Assessment of the degree to which gene expression is additive and heritable has important implications for understanding the maintenance of variation, adaptation, phenotypic divergence, and the mapping of genotype onto phenotype. We used whole-genome transcript profiling using Agilent long-oligonucleotide microarrays representing 12,017 genes to demonstrate that gene transcription is pervasively nonadditive in Drosophila melanogaster. Comparison of adults of two isogenic lines and their reciprocal F1 hybrids revealed 5820 genes as significantly different between at least two of the four genotypes in either males or females or across both sexes. Strikingly, while 25% of all genes differ between the two parents, 33% differ between both F1's and the parents, averaged across sexes. However, only 5% of genes show overdominance, suggesting that heterosis for expression is rare.
THE null hypothesis in quantitative genetics is the infinitesimal model, that there are an infinite number of genes, that the individual effects of genes are vanishingly small, and that gene effects are predominately additive (Fisher 1918). While recent work has challenged the assumption of large numbers of genes per trait and the accompanying assumption that allelic effects are small (Mackay 2001), surprisingly little evidence has been gathered to either support or contradict the assumption of additivity, particularly at the molecular level. The assumption of additivity follows as a parsimonious inference from the observation that phenotypic traits tend to blend in the offspring of divergent parents and is also built into standard mathematical models of complex traits. Yet no genome-wide assessment of the degree of additivity of gene effects has been reported. Gene expression profiling provides one mechanism for performing such a study, recognizing that the path from genotype to phenotype passes through transcript abundance. To address this issue in Drosophila, we combined moderate repetition with an experimental platform that provides high measurement repeatability, so as to quantify subtle differences in gene expression.
Both theoretical and empirical studies in a variety of species suggest that with a half dozen replicates of several genotypes (or other treatments), differences in transcript abundance as small as 1.2-fold can be detected reliably using analysis of variance or Bayesian statistical methods (Kerr et al. 2000; Wolfinger et al. 2001; Efron and Tibshirani 2002). Three studies comparing isogenic strains of Drosophila melanogaster have demonstrated that of the order of 10–20% of all transcripts significantly differ in abundance between adults or developmental stages (Jin et al. 2001; Meiklejohn et al. 2003; Rifkin et al. 2003), and a comparable fraction of the transcriptome varies between specific tissue samples of individual fish, mice, and humans (Pavlidis and Noble 2001; Oleksiak et al. 2002; Cheung et al. 2003; Hsieh et al. 2003). What has not been assessed directly is the additivity or lack thereof of these differences, namely the degree and direction of dominance and epigenetic effects on gene expression.
We approached this question by comparing gene expression profiles of 1-week-old reproductively mature male and female flies of two strains, OregonR and 2b, as well as their reciprocal F1 offspring. These two well-studied strains (Nuzhdin et al. 1993) differ in reproductive and other life history traits and have been used to map QTL for a wide range of traits from bristle number and wing shape to metabolism and fitness in a set of 98 derivative recombinant inbred lines (Gurganus et al. 1998; Zimmerman et al. 2000; Wayne et al. 2001; Montooth et al. 2003). The strains show considerable sexual dimorphism for a range of traits, some of which may show more extreme divergence than that observed between two randomly chosen inbred lines of D. melanogaster. Also, gene expression profiling of midpupae of the two parent strains has been used to narrow down a QTL interval affecting ovariole number to a set of candidate genes (Wayne and Mcintyre 2002).
The expression platform that we present is a custom-designed long-oligonucleotide microarray printed in situ by Agilent Technologies (Hughes et al. 2001). The array provides nearly comprehensive coverage of the version 3.0 release of the Berkeley Drosophila Genome Project annotation of the Celera genome sequence of D. melanogaster (Celniker et al. 2002). A complete loop design (Churchill 2002; Simon et al. 2002) involving 28 microarrays with seven replicates of each of eight genotype and sex samples was analyzed with gene-specific mixed models using the SAS Microarray Solution. Surprisingly, twice as many genes were observed to differ in gene expression between F1 flies and their parents as between the two parent strains. We describe this observation in the context of modes of dominance and relative to the contribution of each chromosome and class of gene ontology, as well as the differences between the sexes.
MATERIALS AND METHODS
Each microarray was synthesized by Agilent Technologies using phosphoramidite chemistry to synthesize 60 mers in situ on glass slides (Hughes et al. 2001). Probes were designed using release 3.0 of the BDGP annotation of the Celera D. melanogaster genome sequence, and sequences will be available to the public through Agilent as part of the annotation of the commercial release. The 21,929 Drosophila probes represent 12,559 genes with 13,631 unique probes; 3112 of the probes are represented by between two and four replicates that are randomly distributed across the layout, and 3189 of the genes are represented by duplicate probes, 1072 of which are alternate probes. Most probes were designed against a common exon, while the alternate probes typically correspond to an alternate splice variant. Average expression levels derived from seven measurements of each sex and genotype were extremely highly correlated between duplicated probes (R2 > 0.98; data not shown) whereas alternate probes gave uncorrelated measurements, as expected. Oligonucleotide sequences and gene identities are supplied as supplementary information, as is an SAS file containing all fluorescence intensity measurements, at http://statgen.ncsu.edu/ggibson/SupplInfo/SupplInfo3.htm.
The eight genotype and sex combinations (OregonR, 2b parents and their reciprocal F1's, O2b and 2bO; 7- to 10-day-old mated males and females) were each represented by seven hybridizations on 28 microarrays, according to a loop design (Churchill 2002; Figure 1). Each sample was generated by extracting total RNA from several hundred flies combined from multiple (5–10) rearing vials; three independent amplified RNA-labeling reactions were performed for each fluorescent dye (Cy3 or Cy5, using Agilent kit G2554A according to manufacturer's instructions). To minimize the source of technical error due to labeling, these were pooled and then split into three or four separate hybridization mixes for each of the 16 sex, genotype, and dye samples. The loop design ensured that each sex and genotype was contrasted directly against each other sex and genotype and that the dyes were balanced for sex, genotype, and sex-by-genotype contrasts. Parental and F1 flies were both reared at a density of ∼100–150 larvae per bottle on standard cornmeal medium (1 kg cornmeal, 200 g yeast, 100 g agar, 700 ml molasses, and 14 liters water, plus tegosept and propionic acid) supplemented with live yeast paste. Adults were maintained in vials on a 12-hr light-dark cycle, separated by sex between days 4 and 7 after eclosion, and collected over a 3-hr midafternoon window to minimize circadian influences, by snap freezing in 70% ethanol at day 7–10. Total RNA was extracted using Trizol reagent. All hybridizations were performed at the Agilent facility in Palo Alto, California, according to standard protocols.
Slides were scanned using an Agilent scanner (G2565BA), and all analyses were performed on background-subtracted, loess-transformed log base 2 fluorescence intensity measurements. Two-step mixed model analysis of variance (Wolfinger et al. 2001) was performed with SAS Microarray Solution software. Each of the 56 hybridization samples was subjected to linear normalization in PROC MIXED with the model where i (= 1–28) specifies the microarray number, j specifies the dye(channel), and Yij are the log2-transformed expression intensities. The residuals from this model were then analyzed in 21,929 separate mixed models by probe using a model of the form where k specifies the genotype (ORE, O2b, O2b, or 2bO), l is the sex (male or female), and Array is the only random effect. The significance of each of the 12 specific contrasts of the four genotypes within each sex was then determined using the ESTIMATE option in PROC MIXED (SAS code is available on request). Standardized least-squares means for each sex and genotype were subjected to two-way hierarchical clustering using Ward's method in JMP to generate Figure 2. A total of 15,165 probes showed a difference between the sexes less than the Bonferroni significance cutoff of 2.3 × 10−6 (0.05/21,929 probes) and these correspond to 8775 genes after accounting for duplicate probes and alternate transcripts.
Probes were excluded from further analysis if none of the 12 pairwise genotype-within-sex contrasts showed a significant difference that was <1.9 × 10−7 [0.05/(21,929 probes × 2 sexes × 6 pairwise comparisons)]. This gave a list of 9557 potentially significant probes, corresponding to 5820 genes. The probes were sorted by gene, and the specific probe for each gene that gave the highest mean significance across comparisons was chosen for subsequent analysis of inheritance patterns. The justification for this is that repeated measures were extremely consistent. This analysis excludes comparison of alternate splice variants. For describing patterns of inheritance, we used mean least squares of the transformed data to estimate transcript abundance. As demonstrated in Ranz et al. (2004), slightly different formulation of the mixed models and choice of significance threshold or fitting of a Bayesian approach will alter the categorization of a small percentage of genes. Our analysis thus focuses on the pattern of variation across the genome, rather than on results for specific genes, and on fold changes in transcript abundance, rather than on P-values, given that all of the genes in the analysis show a highly significant genotype effect across the four parental and F1 genotypes.
Genes were classified into the inheritance patterns in Figure 4 and Table 1 by computing the difference between least-squares mean expression level estimates for each genotype within sex. Differences >1.25-fold (a difference in log2 intensity of 1.32 or more) typically corresponded to a significance of 10−5 and were used to assign genes to categories according to the following intuitive rules using an Excel macro: additive genes are those where both F1 genotypes are 1.25-fold greater than one parent and less than the other; overdominant (underdominant) genes are 1.25-fold greater (less) than both parents; dominant high (or low) genes are not significantly different from one parent but greater (less) than the other; genes expressed in the F1 more like the same sex in the parental line that was female (or male) are not significantly different from the female (male) parent of the same sex, but greater than both the reciprocal F1 and the other parent of the same sex; and genes that are high (low) in one F1 are greater (less) than each of the other three genotypes (and in some cases, one of the other genotypes might also have been significantly different from the others). Some of these categorizations are somewhat arbitrary, but our analysis tends toward conservative estimation of the number of genes showing nonadditive inheritance patterns. Approximately half of the significant genes did not fit one of the 16 patterns above. These generally showed differences only between the two most extreme genotypes: many of these may also be additive or show partial dominance that does not satisfy the conservative definitions above. The remaining genes (7519 in females, 8461 in males) did not show significant genotype differences in either sex. Functional categorization of all genes by manually curated FlyBase Gene Ontology (GO) category is indicated in the supplementary table online at http://www.genetics.org/supplemental/.
Sex-specific generational differences in the transcriptome:
To estimate how many genes show sex and genotype differences, we assessed for each gene separately the significance of the difference between genotype means, relative to the variance in measurement within each genotype. Technical error was minimized by pooling multiple labeling reactions in each hybridization mix. Of the 12,017 genes on the array, 8775 are sexually dimorphic in adults after correction for multiple comparisons, and 5820 showed a significant difference between at least one of the four genotypes in one of the sexes. Two-way hierarchical clustering by gene and by genotype/sex combination reveals which genes share similar patterns of inheritance and provides an overall view of the extent of departure from additivity (Figure 2). The eight rows corresponding to the genes that show significant differences between at least one genotype in one sex cluster at two levels. The deepest branch separates males and females, consistent with previous observations that at least half of the whole adult transcriptome in flies is differentially expressed in the reproductive tissue (Jin et al. 2001; Arbeitman et al. 2002; Parisi et al. 2003). In both sexes, a second branch separates the reciprocal F1's from the two parental strains. This is a broad indication that transcription is more divergent between the F1 and parents than between the parents and, hence, that departures from additivity are apparent for a substantial fraction of the transcriptome.
Genotype differences are visualized as volcano plots in Figure 3 that contrast significance against magnitude of difference in transcript abundance. Points above the dashed horizontal line in each plot exceed the conservative Bonferroni significance cutoff of 1.9 × 10−7 (see materials and methods for details). Approximately 25% of the probes differ between the two parents in either sex, 20% between the reciprocal F1 females and 9% between the reciprocal F1 males, but an average of 37% between females of either F1 and either parent or 28% between F1 and parental males. For comparison, two-thirds of the probes differ between sexes. The percentage of genes differing between inbred lines and sexes is larger than that reported in the above-mentioned studies, but is consistent with their findings given our higher level of replication and the reduced technical error of the Agilent arrays: classification of differences in significance is very much a function of experimental design. Dye effects were small relative to biological effects, with only one gene showing a more than twofold sensitivity to whether Cy3 or Cy5 was used in the labeling reaction, although 7% of the probes were significant for the dye effect after Bonferroni correction.
Clusters of genes with similar expression profiles tend to identify patterns of dominance. Each of 16 common modes of inheritance is plotted on a standardized scale of fold change in relative transcript abundance in Figure 4, showing the mean value of each genotype as a boldface line and the range of values as fine lines. Genes showing additive gene expression have transcript abundance estimates intermediate between those of the two parents, whereas dominance is indicated where transcript abundance of the two F1's is similar to one or the other parent. Overdominance is the situation in which transcript abundance is at a higher level in both reciprocal F1's than in both parents, whether or not the parents differ in expression, while underdominance is the reverse situation, with transcript abundance lower in both reciprocal F1's than in the two parents. Female parent-like effects are defined as the case in which the F1 abundance resembles that of the same sex in the line that donated the female parent, and male parent-like effects are defined as the case in which abundance resembles that of the same sex in the line that donated the male parent. (Note that these are not quite the same as maternal or paternal effects.) Numerous examples where just one of the F1 crosses had an extreme high or low expression level were also observed. The three numbers above each plot indicate the number of genes showing the inheritance pattern in females, in males, and in both sexes.
Several features of the analysis of genotype differences stand out. First, 25% more genes are differentially expressed in females than in males, and the total fraction of significant Drosophila genes is around one-third of the genome. Second, <2% of the differentially expressed genes are additive in the sense that both F1's are at least 1.25-fold different from both parents. Many of these actually show partial dominance, as indicated in histograms of the ratio of dominance to additivity coefficients for each sex in Figure 5. Third, there is asymmetry in the direction of full dominance. For example, the 2b parent has higher transcript abundance in more than twice as many genes in males as in females (see numbers above each panel in Figure 4). Fourth, 5% of the genome shows over- or underdominance. Almost half of these effects are observed in both sexes, implying that they may be acting in the soma rather than in germline. Fifth, a remarkably large number of genes are significantly higher in F1 females derived from a female 2b parent than in the other genotypes. Sixth, only half of the differentially expressed genes show these standard patterns of expression: another 2247 genes in females and 2012 in males have less obvious patterns of differential expression between just two of the genotypes. The implications of these findings are discussed below.
Sex biases in nonadditivity of transcription:
Next we asked whether sex chromosomes and autosomes contribute equally to the patterns of inheritance. Disregarding the genetically depauperate Y chromosome, D. melanogaster has six chromosome arms, five of which (X, 2L, 2R, 3L, and 3R) each carry between 15 and 25% of the 13,500 or so genes, while the short chromosome 4 has only ∼100 genes. Table 1 suggests that there is a weak bias against genes on the X chromosome showing genotype-specific differences: only 85–90% of the expected number of X chromosome genes are significantly different for any of the common patterns of inheritance given the number of X-linked genes included on the array, whereas each mode of variation for loci on the autosome arms is represented approximately in proportion to its gene density on the microarray. This pattern can be explained by noting that the X chromosome is hemizygous every other generation. Selection is thus more efficient at removing deleterious and fixing beneficial recessive alleles, which should reduce the amount of standing genetic variation due to nonadditive alleles. However, the paucity of variation for expression on the X appears to stem from underrepresentation of additive, not dominant, variation: very few of the “additive” genes are X-linked (7 of 130 in females, 0 of 137 in males).
The tendency toward over- or underdominance appears to be strongly sex biased. Table 2 shows that among genes that are at least 2-fold different between the sexes, and 1.25-fold different between parents and F1, female-biased genes (expression higher in females) are more likely to be expressed at a higher level in the parents than in the F1, while the reverse is true for male-biased (expression higher in males) genes. This result holds in both sexes, but is much stronger in males: 80% of all female-biased genes are downregulated in F1 males, and 63% of male-biased genes are upregulated. The corresponding percentages in F1 females are 57 and 56%. A similar result is observed in the interspecific hybrid females derived from a cross of D. melanogaster and D. simulans (Ranz et al. 2004) where it has been attributed to a combination of underdevelopment of the ovary and loss of control of male-biased expression. It would seem, though, that there are intrinsic tendencies toward misexpression even in intraspecific hybrids. An alternative explanation might be that expression is at the appropriate level in heterozygous flies, but deviates in isogenic stocks due to deleterious recessive alleles. This in turn implies that misregulation of male-biased genes more frequently causes downregulation, while misregulation of female-biased genes causes upregulation.
Functional clustering of differentially expressed genes:
Careful examination of the list of genes within each pattern of inheritance highlights numerous instances of possible coregulation of functionally related genes. Some of these are located close together on a chromosome as observed by Spellman and Rubin (2002): for example, a cluster of five odorant-binding proteins in cytological intervals 55 and 56 that show paternal inheritance of expression level from 2b. In other cases, similar patterns of regulation are observed in genomically distant but functionally related genes, e.g., five unlinked ubiquitin conjugating enzymes on different chromosome arms show underdominant expression. A similar phenomenon may be observed for genes involved in a common biochemical or cellular process; e.g., of 19 differentially expressed genes involved in EGFR signaling, 11 show underdominance and 6 overdominance. It will be interesting to determine whether such relationships hold over a series of crosses involving diverse lines.
We also asked whether any categories of functional annotation had obvious trends with respect to sex dimorphism or particular patterns of inheritance. Table 1 lists the percentage of genes in each of 14 categories based on the GO annotation of the fly genome. Most categories are equally represented in the two sexes, but one deviation stands out. Transcripts encoding proteins involved in translation are overrepresented in the genes that are differentially expressed across genotypes in females, as expected since the nurse cells must prepare the oocyte for rapid protein synthesis upon fertilization. Many of these genes are also strongly female biased. A particularly interesting class of these genes is a collection of 13 aminoacyl-tRNA ligase mRNAs that show significantly elevated levels specifically in F1 females derived from the cross of 2b females by OregonR males.
Classical quantititive genetic theory is formulated without regard to the molecular events that relate genotype to phenotype. Although it has been known for 50 years that this relationship is mediated through transcription, translation, and protein interactions, only recently have the tools been developed that allow the additivity and dominance of gene activity to be addressed in a systematic manner at the transcriptional or proteomic level. Fundamentally, the assumption and very general observation that gene effects are additive intuitively implies that transcript abundance itself is also additive. Dominance at the transcript level but not the trait level would imply modification of the transcriptional effect by other genes, which while plausible is not consistent with basic quantitative genetic models. Similarly, for dominant traits, differences between reciprocal F1's or the two sexes at the transcript level are not generally expected. While we have not attempted in this study to relate transcriptional to phenotypic variation, our results are nevertheless challenging because they indicate extensive nonadditivity at the transcription level, implying that the assumption of a straightforward mapping of genotype onto phenotype is not justified.
Another way of considering this relationship is to consider that genetic polymorphism induces molecular polymorphism, some summation of which results in phenotypic polymorphism at the trait level. The simplest conception of additive effects is that some polymorphisms affect transcript abundance, some affect protein abundance, some affect catalytic rates, and so forth, so that in moving from genotype to trait we would expect a gradual increase in the level of variance at each successive biochemical level. The finding that there is a higher coefficient of variation associated with gene effects at the transcriptional or protein level than at the phenotypic level would not be consistent with this simple conception. It would likely imply that additivity arises as a result of the averaging of a complex web of nonadditive biochemical effects. We suggest that concerted attention to variability at all levels of biological organization is now warranted, as it has applied and theoretical implications.
As one example, consider that when two inbred strains of any species are crossed together, it is not uncommon for several traits to show improved performance or heterosis. The two most often cited explanations for this phenomenon are unmasking of the effects of deleterious recessive alleles and summation of the dominant effects of multiple loci brought together in the progeny. Neither of these explanations requires over- or underdominance of transcript abundance; nor do they predict it. On the other hand, heterosis and its negative counterpart, outbreeding depression, are so important in plant and animal breeding and so central to evolutionary theory that it is worthwhile to characterize the effect of out-crossing on the expression of genes. We find very little evidence of heterosis for gene expression (only 5% of the transcriptome), although there are multiple other departures from additivity and many instances of altered gene expression in one of the reciprocal F1's or in which the two F1's resemble alternate parents. The percentages of genes in each class may be surprisingly low to some readers and high to others, highlighting the fact that we have little empirical or theoretical evidence on which to base expectations. Comparison of heterosis for expression of specific genes with heterosis for traits regulated by the genes is an obvious research priority.
Whether a transcript shows some degree of dominance for abundance is likely to be influenced by the relative contributions of cis- and trans-acting factors to expression of the gene. If transcription is controlled predominantly by cis-regulatory regions of a particular gene, then transcript abundance might often be expected to be additive in the absence of transvection. By contrast, transcription factors that vary in activity level between the parents are more likely to interact to produce a range of degrees of dominance. Dominance for low-level transcription might be caused by increased repressor activity or, alternatively, by haploinsufficiency for a transcriptional activator. Extreme expression levels in either F1 or similarity to expression in the same sex of either parental line is more difficult to explain, although X linkage and genomic imprinting may contribute. However, genomic imprinting has recently been documented in Drosophila (Maggert and Golic 2002), but only for the Y chromosome. Maternal genetic effects are more commonly considered in relation to embryogenesis, but differences in gene expression of adults might also arise as a result of initial discontinuities in the embryo persisting and even amplifying throughout development.
Little theoretical work is available to guide interpretation of departures from additivity of expression. Arguments from biophysical principles relating to the cooperativity of transcription factors (Gibson 1996) as well as the structure of regulatory networks (Omholt et al. 2000) imply that epistasis (sensu molecular) and dominance are likely to be natural properties of transcriptional regulation. Recent modeling of development with systems of linear equations (Von Dassow et al. 2000; Meir et al. 2002) demonstrates that regulatory networks can be robust to considerable variation in parameters representing transcript and protein abundance, but has not addressed the causes and consequences of dominance. Empirically, the combination of QTL mapping with expression profiling in recombinant inbred lines of mice (Schadt et al. 2003) has demonstrated that cis-acting variation contributes 25% or more of the variation for expression of 15% of genes in the liver. This study also found a number of hotspots for the location of eQTL (QTL affecting gene expression), and these are likely to correspond to the location of trans-acting regulatory genes. A similar result was obtained in crosses between two yeast strains (Yvert et al. 2003). Our results imply that inclusion of F1 individuals in studies mapping gene expression differences could greatly assist in the interpretation of sources of variation affecting divergence among lines.
From an evolutionary genetic perspective, documentation of extensive nonadditivity of gene expression has at least two important implications. First, by challenging the assumption that genotype maps directly onto phenotype through transcript abundance, it calls for more attention to the prevalence of epistasis as a pervasive aspect of genetic architecture. The possibility that genotypes produce additive variation for visible traits irrespective of nonadditive variation for transcription and, conversely, that cryptic variation can be hidden at the transcriptional level should be considered explicitly by theoreticians and experimental biologists working with different organisms. Second, the widespread differences between F1 and parents suggest caution in relation to the use of highly inbred lines to quantify levels of intraspecific variation, at least for gene expression. The degree of differentiation among inbred lines reported here is similar to that documented by others (Jin et al. 2001; Rifkin et al. 2003), but more extensive sampling in the form of a diallel cross will be required to assess the generality and covariance of nonadditive gene expression in Drosophila.
We thank Russ Wolfinger for providing a beta version of the SAS Microarray Solution software. This work was funded by National Institutes of Health award R24-GM65513.
Communicating editor: D. Weigel
- Received January 19, 2004.
- Accepted May 7, 2004.
- Genetics Society of America