The deleterious consequences of inbreeding, especially in the form of inbreeding depression, are well known. However, little is known about how inbreeding affects genome-wide gene expression. Here, we show that inbreeding changes transcription levels for a number of genes. Gene expression profiles of Drosophila melanogaster lines inbred to F ≈ 0.67 at different rates changed relative to those of noninbred lines, but the rate of inbreeding did not significantly affect gene expression patterns. Genes being differentially expressed with inbreeding are disproportionately involved in metabolism and stress responses, suggesting that inbreeding acts like an environmental stress factor.
INBREEDING depression is caused by an increase in the homozygosity of recessive deleterious alleles and loss of overdominance at some loci due to increased homozygosity (Charlesworth and Charlesworth 1987). Empirical studies show the ubiquity of inbreeding depression for traits related to fitness (Charlesworth and Charlesworth 1987; Crnokrak and Roff 1999; Hedrick and Kalinowski 2000; Keller and Waller 2002; Kristensen and Sørensen 2005). However, theory predicts and empirical studies have shown that the level of inbreeding depression is environmental and line specific (Bijlsma et al. 1999; Fowler and Whitlock 1999; Reed et al. 2002; Kristensen et al. 2003; Vermeulen and Bijlsma 2004) and dependent upon the level of genetic load (Ehiobu et al. 1989; Lacy and Ballou 1998; Wang et al. 1999; Day et al. 2003; Reed et al. 2003; Pedersen et al. 2005). Differences between lines in the effects of inbreeding may partly be determined by the rate of inbreeding. If inbreeding is sudden and extreme, the effective population size is strongly reduced, random fixation occurs more often, and selection will have minor impact (Hedrick 1994; Fu et al. 1998; Wang et al. 1999). Consequently, because there are more generations and greater opportunity for selection to act before a given inbreeding level is reached, slower inbreeding is predicted to cause less inbreeding depression than an equivalent level of fast inbreeding (Robertson 1952; Ehiobu et al. 1989; Wang et al. 1999; Day et al. 2003; Reed et al. 2003; Pedersen et al. 2005).
Molecular studies of aging, inbreeding, and environmental stress have been shown to induce similar cellular responses (Kristensen et al. 2002; Pletcher et al. 2002; Sørensen et al. 2003; Girardot et al. 2004; Landis et al. 2004; Pedersen et al. 2005). For instance, the molecular chaperone Hsp70 is found to be upregulated in response to numerous environmental stresses (Feder and Hofmann 1999; Sørensen et al. 2003) and is also upregulated in some inbred lines (Kristensen et al. 2002; Pedersen et al. 2005). However, no studies have investigated how inbreeding affects expression levels of the whole genome. Here, gene-expression profiles of lines inbred to the same level at different rates and noninbred Drosophila melanogaster lines are investigated. On the basis of the obtained results we conclude that inbreeding leads to differential expression of a wide variety of genes with disproportionate representation involved in metabolism and stress resistance.
MATERIALS AND METHODS
Inbreeding procedure and maintenance of the lines:
A genetically diverse mass population of D. melanogaster was founded in August 2002 by mixing 600–700 flies from each of four sets of preexisting populations collected in Denmark, Australia, and The Netherlands. The stocks were all maintained at high population sizes (N > 1000) prior to crossing.
Inbred (fast and slower rate) and “noninbred” control lines were founded from the mass population in December 2002 eight generations after the mass population was founded. Lines with expected equivalent levels of inbreeding (F ≈ 0.67) were obtained by two different rates of inbreeding, either through five generations of full-sib mating (fast rate) or by maintaining a population size of two pairs during nine generations (slower rate). Five independent inbred lines were generated for each of the two breeding regimes. Each inbred line was founded by randomly selecting, respectively, one male and one female (full-sib) and two males and two females (slower inbreeding) from the mass population. Assuming the inbreeding level of the base population to be zero, the expected inbreeding levels were calculated as a measure of coancestry for the full-sib mating (Ft = (1 + Ft-1 − Ft-2)/4) (Falconer and Mackay 1996) and as a measure of genetic drift in the lines inbred by a slower rate of inbreeding (Ft = Ft-1 + (1 − 2Ft-1 + Ft-2)/2Ne) (Crow and Kimura 1970). For the inbreeding procedures, offspring from each line from each consecutive generation were collected as virgins. Four and two pairs were set up per line within the fast inbreeding and the slower inbreeding regimes, respectively, to reduce the extinction of lines throughout the inbreeding procedure. For each breeding regime offspring from one vial were randomly chosen to establish the next generation of inbreeding. However, some lines went extinct through the process of inbreeding; thus to make sure that enough lines reached the expected level of inbreeding, excess lines were set up. Twenty and 15 independent lines were started to make, respectively, the fast and slower inbred lines to compensate for loss of lines as a function of the intensity of inbreeding. Respectively 20 and 10% of the fast and the slower inbred lines went extinct through the inbreeding process. After reaching the desired level of inbreeding, all lines were flushed to minimum sizes of 500 breeding individuals (within two generations) and transferred to bottles. Five “noninbred” control lines, each founded by ∼500 breeding individuals, were established at the time when the inbreeding procedures were initiated. The control lines and the flushed inbred lines were each kept in 10 bottles and within each line flies from all the bottles were mixed in every generation prior to setting up the next generation. The major features of the design used to establish the experimental lines are summarized in Table 1. Throughout and following the inbreeding procedure all flies were maintained in one climate room (25° ± 0.2°, 50% relative humidity, 12/12-hr light/dark cycle).
Sampling of flies and replication:
The inbred flies had a lower productivity, and the density within bottles was therefore lower. To get around this problem, the number of flies was controlled in all generations so that 20, 25, and 30 parental pairs were set up for egg laying before being discarded 24 hr later, within the control, slower inbred, and fast inbred lines, respectively. The numbers of flies emerging from the bottles were not significantly different across the three treatments (control, 356 ± 6, fast inbreeding, 387 ± 11; slower inbreeding, 357 ± 28). Flies were never exposed to strong crowding.
Twenty virgin male flies were collected from each line by sampling four males (<8 hr old) from each of five randomly chosen bottles. Sampling was done by four people in the afternoon within 3 hr. Flies from the different treatments were sampled in rotating order, so that the time of collection and the person collecting were randomized between the treatments. In each bottle the first flies emerging were used. After sampling, flies were immediately frozen in liquid nitrogen (flies <11 hr old). This procedure was followed for all 15 lines (5 controls, 5 fast inbred, and 5 slower inbred lines). A pool of RNA from 20 flies from each replicate line was hybridized to Affymetrix chips.
For RNA purification, 20 virgin male flies from each line were collected and frozen in liquid nitrogen and stored at −80°. Flies were homogenized with a FP-120 Fast Prep bead beater according to manufacturer protocols (Bio-101, Carlsbad, CA) in 1.5 ml Trizol reagent (Invitrogen, San Diego) and 150 μl chloroform. Labeling, hybridization, and staining were performed essentially as described by Dyrskjot et al. (2003). Briefly, double-strand cDNA was prepared from 5 μg of total RNA using the SuperScript Choice system (Life Technologies) according to the manufacturer's instructions except using an oligo(dT) primer containing a T7 RNA polymerase promoter site. Biotin-labeled cRNA was prepared using the BioArray High Yield RNA transcript labeling kit (Enzo). Following the IVT reaction, the unincorporated nucleotides were removed using RNeasy columns (QIAGEN, Valencia, CA).
Array hybridization and scanning:
Fifteen micrograms of cRNA was fragmented at 94° for 35 min in a final volume of 40 μl in a buffer containing 40 mm Tris-acetate pH 8.1, 100 mm KOAc, and 30 mm MgOAc. Next, 260 μl of 6× SSPE-T hybridization buffer (1 m NaCl, 10 mm Tris pH 7.6, 0.005% Triton) was added and the cRNA was denatured by heating to 95° for 5 min. The hybridization mixture was loaded onto the Affymetrix probe array cartridge (Drosophila Genome Array Version 1) and incubated for 16 hr at 45° at constant rotation (60 rpm). The washing and staining procedure was performed in the Affymetrix Fluidics Station. The probe array was exposed to 10 washes in 6× SSPE-T at 25° followed by 4 washes in 0.5× SSPE-T at 50°. The biotinylated cRNA was stained with a streptavidin-phycoerythrin conjugate, final concentration 2 μg/μl (Molecular Probes, Eugene, OR) in 6× SSPE-T for 30 min at 25° followed by 10 washes in 6× SSPE-T at 25°. An antibody amplification step followed, using normal goat IgG as blocking reagent, final concentration 0.1 mg/ml (Sigma, St. Louis), and biotinylated anti-streptavidin antibody (goat), final concentration 3 μg/ml (Vector Laboratories, Burlingame, CA). This was followed by a staining step with a streptavidin-phycoerythrin conjugate, final concentration 2 μg/μl (Molecular Probes, Eugene, OR) in 6× SSPE-T for 30 min at 25° and 10 washes in 6× SSPE-T at 25°. The probe arrays were scanned at 560 nm using a confocal microscope (Hewlett Packard GeneArray Scanner G2500A).
The data were analyzed using programs developed in R, a programming language and developer environment for statistical computing and graphics (http://www.r-project.org/). Preprocessing of expression values was performed using the robust multi-array analysis (GCRMA) algorithm (Irizarry et al. 2003; Wu et al. 2004). In this algorithm, raw intensity values are background corrected on the basis of a model using sequence information followed by a quantile normalization and a robust multichip fit with median polish (Wu et al. 2004). This algorithm combines the strengths of stochastic-model-based algorithms and physical models and has been shown to be superior in accuracy and precision to other normalization methods such as microarray analysis suite, RMA, and PerfectMatch (Wu and Irizarry 2004). To exclude genes that could not be confidently detected in the data analysis probe, sets with less than three present calls within at least one of the three treatments were excluded (a transcript must be represented on at least three chips within the control, fast inbreeding, or slower inbreeding treatments). The filtered gene set contained 8884 transcripts.
Differential expression was assessed using significance analysis of microarrays (SAM) proposed by Tusher et al. (2001). An overall test of significance for a gene was performed using moderated F-statistics in a multiclass analysis. For each gene the three contrasts (control-slow, control-fast, and slow-fast) were tested for differentially expressed gene transcripts. The moderated F-statistic tests whether any of the contrasts are nonzero for that gene, i.e., whether that gene is differentially expressed on any contrast. Specific treatment contrasts (control vs. slow, control vs. fast, and slow vs. fast) were also tested on the basis of a modified t-statistic using the two-class unpaired analysis. Multiple testing was accounted for by controlling the false discovery rate at 20% for both the multiclass and two-class analyses. The SAM analysis was performed as implemented in the R package called siggenes (Schwender 2004).
Groups of genes being differentially expressed were annotated on the basis of the biological process ontology directed by the Gene Ontology (GO) database (Gene Ontology Consortium 2001). The expression analysis systematic explorer (EASE) application on the DAVID homepage (http://david.niaid.nih.gov/david/ease.htm) (Hosack et al. 2003) was used to test for overrepresentation of genes in given annotation categories. EASE scores were calculated for the likelihood of overrepresentation in the annotation categories.
The probability that the overlap of genes being differentially expressed with both types of inbreeding is different from the number expected by chance was calculated by using Monte Carlo simulations. In each simulation the gene list within each treatment was permuted and the overlap of induced genes was determined. A Kolmogorov-Smirnov test (Conover 1971) was used to determine if the distribution of the within-gene variances in gene expression levels differs significantly among the three breeding treatments.
The Affymetrix array contained 13,966 probe sets representing ∼13,000 unique genes and 8884 genes were left after the filtering process.
The within-gene variance in gene expression levels within both inbred treatments was higher than that within the control treatment (control vs. fast, D = 0.354, P < 2.2e-16; control vs. slow, D = 0.296, P < 2.2e-16; see Figure 1). The within-gene variance in gene expression in the fast inbred treatment was higher than that in the slower inbred treatment (fast vs. slow, D = 0.081, P < 2.2e-16; see Figure 1).
The moderated F-statistics in the multiclass analysis revealed 21 genes as being differentially expressed (Table 2, Figure 2). Of these 12 were significantly downregulated within both the fast and the slower inbred treatments and the remaining 9 genes were significantly upregulated within both inbred treatments (Table 2, Figure 2). Several genes were differentially expressed for the contrasts control vs. fast and control vs. slower inbred lines whereas no transcripts were differentially expressed between fast and slower inbred lines (Tables 3–5⇓⇓, Figure 2). The results are summarized in a Venn diagram (Figure 3). Sixty-seven genes are differentially expressed with both control vs. fast and control vs. slower inbreeding. The probability of the observed overlap arising by chance is small (P < 0.00001) under the assumption that all genes are liable to change with treatment. To assess the potential contribution of a smaller gene pool, one can diminish the number of genes that are capable of changing with the two types of inbreeding while simulating the overlap. If the set of genes is reduced to 5000 it is still highly unlikely to observe 67 genes overlapping (P < 0.00001). The 67 genes differentially expressed in both control vs. fast and control vs. slow inbreeding were all either up- or downregulated in both comparisons. Fifty genes were upregulated and 17 were downregulated. Given that all genes are either up- or downregulated within both inbred treatments inbreeding, per se, not the intensity of inbreeding, appears to determine the up or down change in transcript level. Genes that are significantly differentially expressed with fast or slower inbreeding compared to the controls are in the great majority of cases either up- or downregulated with both types of inbreeding (Figure 4). The log twofold change of the 67 genes being differentially expressed was not affected by the type of inbreeding; 34 genes had a fold change that was higher in the control vs. fast inbreeding comparison whereas 33 genes had a fold change that was higher in the control vs. slower inbreeding comparison.
EASE scores for annotation categories with more genes than expected by chance among the genes being significantly differentially expressed with fast and slower inbreeding and with both are given in Table 5. Classes of genes involved in metabolism, immune, and stress responses are overrepresented.
The large number of genes differentially expressed in this study means that a detailed description of the changes on a gene-by-gene basis would be too extensive to list (for a complete list of differentially expressed genes see supplementary material at http://www.genetics.org/supplemental/). However, we here present gene groups defined by function.
Genes involved in stress resistance and metabolism are disproportionately affected by inbreeding (Tables 2–5⇑⇑⇑). Some groups of genes associated with these biological processes are upregulated, whereas a few are downregulated with inbreeding. Most genes being differentially expressed with either fast or slower inbreeding in this study responded in the same direction with both types of inbreeding (Figure 4). Furthermore, the overlap between genes being differentially expressed with both types of inbreeding is much larger than an expected overlap arising by chance, and all 67 genes responded in the same direction. Given the high level of replication (10 inbred lines and 5 control lines) these results indicate that there is a general effect of inbreeding on gene expression patterns.
The data presented here have not been validated by quantitative (Q)RT-PCR or Northern blot. However, Park et al. (2004) showed that the ratios of gene expression obtained from Affymetrix platforms and QRT-PCR analyses are highly correlated (r = 0.93). Park et al. (2004) and Yuen et al. (2002) also showed that results obtained from Affymetrix platforms underestimated the real expression change as detected by QRT-PCR. This latter result shows that our results can be interpreted as being of conservative nature. Moreover, the consistency in direction and magnitude of altered transcript levels (inbreeding vs. control) for the numerous inbreeding replicates suggests that the observed differential transcript levels are not an artifact. A thorough comparison of other published data with ours would require a full treatment of the raw data for common normalization, similar statistical analysis, etc. Nevertheless, we have performed a simple analysis using our own data set and the list of differentially expressed genes in the study by Landis et al. (2004) investigating differentially expressed genes in response to aging and oxidative stress. The present microarray study revealed that a total of 466 genes were differentially expressed with inbreeding whereas Landis et al. (2004) observed that 913 genes were differentially expressed with aging and that 593 genes were differentially expressed with oxidative stress. Monte Carlo simulations were used to test whether the observed overlap between treatments (inbreeding, aging, or oxidative stress) was higher than expected by chance. The overlap was in all comparisons significantly higher than expected by chance even when the number of genes liable to change with treatment was reduced to 5000 (P < 0.001 in all cases). Thirty-four genes are differentially expressed under all of the conditions investigated in both studies (inbreeding, aging, and oxidative stress), including stress response genes [e.g., Hsp83 (Affymetrix probe identifier, AFFYID: 143198_at), Diptericin (AFFYID: 147473_at), and Defensin (AFFYID: 143607_at)] and metabolism genes [e.g., adenosine (AFFYID: 143062_at) and NAD-dependent methylenetetrahydrofolate dehydrogenase (AFFYID: 151767_at)]. Seventy-two and 78 genes overlap in the comparison of aged and inbred individuals or between oxidative-stressed and inbred individuals, respectively. Heat-shock protein- and immune-response genes appear to be differentially expressed in these comparisons [e.g., Hsp60 (AFFYID: 152031_at), Hsc70 (AFFYID: 143191_at), Thor (AFFYID: 153432_at), and Drosocin (AFFYID: 143609_at)]. This indicates that effects of different kinds of stresses (such as inbreeding, aging, and oxidative stress) may bear similarities and those genes being differentially expressed under such conditions may act to maintain homeostasis in organisms exposed to diverse stresses. Clearly these genes are candidate genes that should be investigated in more detail by investigating their protein products and by performing knockout studies to obtain knowledge about their function.
Transcription of a number of genes coding for antibacterial peptides is upregulated with inbreeding (Tables 3–5⇑⇑). These include Defensin (AFFYID: 143607_at), Drosocin (AFFYID: 143609_at), Diptericins (AFFYID, 143443_at; AFFYID, 147473_at) and Thor (AFFYID: 153432_at), which all have well-described antibacterial functions (Bulet et al. 1999; Beutler 2003; Ganz 2003). In this study inbred and noninbred lines were kept under the same laboratory conditions and infection pressure is not expected to differ between control and inbred lines. However, there is a possibility that the inbred lines are more susceptible to infection and that upregulation of antibacterial peptide gene transcripts is a defense mechanism induced due to bacterial infection. Alternatively, the protein products of this group of genes have more general stress resistance functions and are part of a general stress response. In accordance with this idea, upregulated transcription of antibacterial genes has also been observed in response to aging in D. melanogaster and Caenorhabditis elegans (Pletcher et al. 2002; Murphy et al. 2003; Landis et al. 2004) and in response to environmental stress in D. melanogaster (Kayo et al. 2001; Pletcher et al. 2002; Landis et al. 2004).
A number of studies have investigated the association between heterozygosity and disease resistance and susceptibility (Hedrick et al. 2001; Giese and Hedrick 2003; Reid et al. 2003), but no consistent pattern emerges from those studies. Giese and Hedrick (2003) recently showed that noninbred populations of the Gila topminnow (Poeciliopsis occcidentalis) have significantly higher disease susceptibility compared to inbred populations. This result may be understood partly in light of our results showing an upregulation of genes coding for antibacterial peptides in inbred populations, whereby immunity toward pathogens may be increased in inbred populations.
Several gene transcripts coding for molecular chaperones such as Hsp60 (AFFYID: 152031_at), Hsp83 (AFFYID: 143198_at), and the heat-shock protein cognate 1 (AFFYID: 143191_at) are differentially expressed within one or both of the inbred treatments in this study. Heat-shock proteins are well known for their importance as part of the cellular stress response apparatus (Feder and Hofmann 1999; Sørensen et al. 2003), and Hsp70 protein level in inbred D. melanogaster and D. buzzatii lines has previously been shown to be higher than that in noninbred lines (Kristensen et al. 2002; Pedersen et al. 2005).
The reason why different rates of inbreeding were investigated in this study was that slower inbreeding is expected to be less deleterious than faster inbreeding for the same level of inbreeding (Robertson 1952; Ehiobu et al. 1989; Day et al. 2003; Reed et al. 2003; Pedersen et al. 2005). This is because with slower inbreeding there are more generations and greater opportunity for selection to act before a given inbreeding level is reached. The fast inbred lines investigated here have been tested for fertility and heat resistance in another study (Pedersen et al. 2005). Pedersen et al. (2005) showed that fertility, but not heat resistance, was significantly affected by the rate of inbreeding with the slower inbred lines having higher fertility. We hypothesized that purging of deleterious alleles within the slower inbred lines would cause changes in gene expression patterns between the two inbred treatments. More gene transcripts were differentially expressed with slower inbreeding, but there were no genes being significantly differentially expressed between inbreeding treatments (Figure 3). One reason for the higher number of differentially expressed genes in the slower compared to the fast inbred treatment probably is that the variance in gene expression is higher within the fast inbred treatment. This means that for a gene to be significantly differentially expressed a higher-fold change between the control and fast inbred treatments is needed compared to the situation for the control and slower inbred comparison (Figure 1). Another reason for the apparent lack of difference in expression may be that the slower inbreeding treatment investigated here is actually still extreme compared to most situations in nature or in domestic livestock. It would be informative to perform the same experiment on lines being inbred fast and slower than the “slow inbreeding” regime investigated here.
Between-line variance in phenotype is expected to increase with inbreeding (Falconer and Mackay 1996; Lynch and Walsh 1998; Kristensen et al. 2005). However, this is the first experiment showing this on the level of gene expression (Figure 1). The reason why the variance in gene expression is higher with fast compared to slower inbreeding may be that selection is more efficient the slower the rate of inbreeding, while drift will become more important the faster the rate of inbreeding. Given that selection regimes were similar across treatments and lines, this would cause slower inbred lines to have more similar gene expression patterns than fast inbred lines, and therefore the observation made in this study confirms theoretical predictions.
Genetic drift is expected to cause fixation of different genes within the different inbred lines. Given the number of replicate lines and the fact that the set of genes found to be differentially expressed is not a random sample of the gene pool but primarily related to metabolism and stress resistance, we find it unlikely that drift alone can explain our results. This emphasizes that there is a general effect of inbreeding that is an indirect result of the change in genotype frequencies. Perhaps the cumulative fixation of deleterious alleles results in net physiological duress to which the organism responds in a standard manner. Our results also show that transcriptional responses to inbreeding overlap with microarray studies of aging and oxidative stress (Pletcher et al. 2002; Landis et al. 2004). We argue that those genes found to be differentially expressed with inbreeding may be candidate genes for stress resistance in general, and we expect the results to have important implications for other disciplines such as medicine, animal breeding, and conservation biology.
We are grateful to Jesper Dahlgaard, Just Justesen, Morten Mulig Nielsen, and Jesper Givskov Sørensen for fruitful discussions on various parts of this experiment; to Doth Andersen and Bente Devantié for excellent technical assistance; to Stuart Barker, Richard Frankham, Lawrence Harshman, Bob Krebs, and two anonymous reviewers for critical comments on earlier versions of the article; and to the Danish National Research Council for financial support via a centre grant.
Communicating editor: L. Harshman
- Received December 10, 2004.
- Accepted May 30, 2005.
- Copyright © 2005 by the Genetics Society of America