An attempt to select non-genetic variation in resistance to starvation and reduced chill coma recovery time in Drosophila melanogaster

Phenotypic variance observed in quantitative traits is classically additively split into genetic (G), environmental (E) and G by E (G×E) interaction components. Among these, only the additive contribution to genetic variance is transmitted to the next generation. Non-genetic inheritance is suggested to explain part of the phenotypic variance that is observed in nature (Salinas et al., 2013). Epigenetic marks play an important role in this non-genetic inheritance as DNA methylation patterns or chromatin conformation can be transmitted across generations (Herman and Sultan, 2016). In addition, it is becoming clear that maternally or paternally transmitted small RNAs can play an important role in the maintenance of gene expression patterns (Watanabe et al., 2011; Conine et al., 2013; Holoch and Moazed, 2015). Moreover, microbes may also participate in non-genetic inheritance, often exhibiting vertical transmission and affecting host gene expression (Vastenhouw et al., 2006).

While it is clear that non-genetic inheritance contributes to phenotypes, despite some attempts, it has remained difficult to quantify this contribution. For example, the part of phenotypic variation that can be explained by changes in DNA methylation was estimated using epiRILs (epigenetic recombinant inbred lines) in Arabidopsis thaliana (Johannes et al., 2009; Zhang et al., 2013; Cortijo et al., 2014). The authors showed that a significant percentage of variance can be explained by the epialleles. Using the same epiRILs, the non-genetic heritability of several traits (such as leaf area or flowering time) was estimated to be low but significant (Zhang et al., 2013; Kooke et al., 2015). Epialleles were thus demonstrated to play a part in the evolution of the organisms (Zhang et al., 2013; Kooke et al., 2015). No equivalent experiments have been done in animals, including in Drosophila.

With this work, we intended to test the hypothesis that non-genetic inheritance can play a part in phenotype and that it can be selected. If this is the case, it will provide us with new mechanisms to understand how species adapt to different environments. To test this hypothesis, we performed artificial selection experiments on Drosophila melanogaster inbred lines, which harbour low levels of genetic variability. Drosophila melanogaster is one of the best studied model organisms in quantitative genetics, and a great number of selection experiments have been performed for a large number of traits (Harshman and Hoffmann, 2000). Here, we selected for a decrease in chill coma recovery time (CCRT) and an increase in starvation resistance (as survival time, ST). Indeed, CCRT and ST are often recorded as displaying high levels of heritability (Ayrinhac et al., 2004; Hoffmann et al., 2005), and the underlying mechanisms are starting to be elucidated (Slocumb et al., 2015; Hardy et al., 2018). If non-genetic inheritance is not transmitted to the next generation, the selection procedure should fail. Because the selection procedure can be considered a stress condition, we estimated developmental instability using bristle number and fluctuating asymmetry measures, which can indicate a break of the canalization process. We also measured the expression of candidate genes known to be involved in the stress response, such as Hsp, and genes implicated in the epigenetic pathways involved in transposable element (TE) silencing.

As previously shown in other organisms, we found that response to selection was strongly genotype dependent (Groot et al., 2017; Herman and Sultan, 2016). However, and contrary to what has been shown in plants (Cortijo et al., 2014), the extent of the response to selection was weak in our D. melanogaster lines. It was previously shown that morphological and fitness alterations can occur through the selection process and remain after selection relaxation in Drosophila (Rutherford and Lindquist, 1998; Sollars et al., 2003). This was not an obvious result in our experiment. However, we did detect expression changes for some genes following the selection process, such as thor, Hsp27 and ago3.

Samples of natural populations of D. melanogaster were collected from a single population at Gotheron, France (44°56′0″N, 04°53′30″E) in June 2014 using fruit bait. Thirty isofemale lines were established directly from gravid females from the field. Brother–sister matings were subsequently performed for 30 generations, resulting in 30 inbred lines, which are thought to harbour very low intra-line genetic variability (Falconer and Mackay, 1996). We randomly chose three of these lines to continue with the experiments, denominated 6.6, 10.1 and 15.4. An outbred line was built from one virgin pair sampled from each of the 30 isofemale lines. The progeny arising from these flies constituted generation 0. This outcome provided a control line displaying initial genetic variability, which we expected to respond to selection. Indeed, this initial genetic pool was made of a variety of alleles, either favourable or unfavourable to starvation resistance or chill coma recovery time. The selection process applied on this outbred line allowed the removal of unfavourable alleles, and therefore made the genetic pool evolve in the intended direction. As inbred lines display non-genetic variability, whereas the outbred line displays both genetic and non-genetic variability, response to selection was expected to be stronger in the outbred line than in the other lines. Flies were maintained in the laboratory at 24°C in a standardized culture medium for Drosophila.

We only considered 2–5 day old flies, as previous studies found that CCRT depended on age (David et al., 1998). Flies were first sexed on ice within a 5°C chamber, and 50 females were then transferred into empty plaques (one female per well) and placed in chambers containing melting ice. After 16 h, individuals were promptly removed from cold to room temperature (24°C), and CCRT was measured individually by recording the time until the fly could stand on its legs (Gibert et al., 2001).

Flies were put into starvation vials (1.5% agar medium) with no nutritional value, only allowing the flies to obtain ingestible moisture (Bubliy and Loeschcke, 2005; Harshman et al., 1999; MacMillan et al., 2009). The 2–5 day old females were sampled using an insect aspirator, without anaesthesia. Five tubes were established, each containing 10 females and kept at 25°C with 70% relative humidity. The number of dead flies was recorded three times per day.

Artificial selection was applied during 10 generations, without relaxation, on samples of 50 females.

For CCRT decrease, a selection pressure of 20% was applied, i.e. the first 10 recovering females were used as the breeders for the next generation. Ten other pairs were randomly selected from a pool of flies to make up the control lines. Each selected line had its own control, maintained in the same conditions (temperature, density and culture medium), except that it was subjected to neither chill nor starvation, with the aim of minimizing the effects of microenvironmental variations. Individuals that took longer than 120 min to recover from chill coma were excluded from the analysis.

For starvation resistance, approximately 50% of surviving females (L50) from the five replicates were used to make the subsequent generation. As selection was performed only on females, at each generation, males were randomly taken from the pool of males in the same numbers as females. After each treatment, flies were placed into vials containing fresh food and survival was measured 24 h later. All control lines were kept with the same number of flies as the selected lines at each generation. Thus, we ensured that density was not a variable to be considered in the analyses.

All measures were recorded 10 generations after the end of the selection process; that is, at generation 20. Thirty females per line were analysed, for a total of 480 individuals. All visible morphological alterations (deformed body parts) were recorded (Table S1). We also recorded the number of scutellar bristles (SCT), which is considered a strongly canalized trait (Rendel, 1959; Sgrò et al., 2010). Indeed, changes in SCT number would be reflective of alterations of canalization.

We counted the number of sternopleural bristles from both the left and right sides of individual flies, and estimated FA as the absolute value of their difference (Van Valen, 1962; Palmer and Strobeck, 1986).

Twenty 2–4 day old mated females were placed into four vials (five flies per tube) to lay eggs for 48 h, and then placed into new vials for 48 h. Vials were maintained at 24°C in a standardized culture medium. Hatching adults were counted daily for 8 days. We used the total number of adults in the progeny as a proxy for the fitness of the line.

We chose to work on ovaries because it is the female tissue that is most closely related to trans-generational transmission. Pools of 70 pairs of ovaries from 2–5 day old mated females were dissected in 1% PBS and split into three biological replicates. Total RNA was extracted using the RNeasy kit (Qiagen) plus QIAzol Lysis Reagent and subsequently treated with DNase (DNA-free kit; Ambion); 1 μg of total RNA was then converted into cDNA using Superscript III reverse transcriptase (Invitrogen) primed with oligo(dT)₂₀.

cDNAs were 50-fold diluted and then quantified using SYBR green quantitative PCR (qPCR) in a LightCycler apparatus (Roche Diagnostics). We chose genes that have been shown to be implicated in the canalization process and stress response [Hsp26, Hsp27, Hsp68, Hsp83, technical knockout (tko)] and in epigenetic regulation [argonaute 3 (ago3), piwi, zucchini (zuc), vasa (vas), thor, Methyltransferase 2 (dnmt2), Suppressor of variegation 3-9 (suvar3-9), modifier of mdg4 (modmdg4), oskar (osk), Helicase at 25E (Hel25E)], and the 412 TE. All these genes are expressed in the ovaries. We tested three genes for use as normalization genes: Ribosomal protein L32 (rp49), 18S rRNA (18S) and tubuline (tub). As rp49 had the lowest coefficient of variation across samples, we normalized the whole dataset by rp49 levels. Three biological replicates were obtained for each condition, and reactions were performed in triplicate. Standard curves were calculated using serial dilutions of the cDNAs. An efficiency value between 1.8 and 2.0 was maintained. For each gene, standard curves were used to convert Ct values into absolute concentrations, which were subsequently divided by rp49 absolute concentrations in each sample. Primer sequences are given in Table S2.

To measure intra-line nucleotide diversity, we focused on one intron of the Adh gene. This gene is particularly well studied in Drosophila for diversity measures (McDonald and Kreitman, 1991). Working on the intronic sequence gave us access to neutral diversity. We performed individual DNA extractions on five females per line using the NucleoSpin 96 Tissue Kit (Macherey Nagel, Düren, Germany) following the manufacturer's instructions. We PCR amplified a 453 bp region in the first intron of the Adh gene (FBgn0000055) (see Table S2 for primer sequences). PCR products were directly sequenced using the Sanger procedure.

Estimates of the realized heritability for both traits were computed for each line by regression of the cumulative selection response (as a deviation from the control) on the cumulative selection differential based on the data from generation 0 to 10 (Falconer and Mackay, 1996). Because the resulting selected lines were derived from inbred lines with different genetic backgrounds, we computed heritability estimates for each line separately.

All statistical analyses were performed using R. Data were analysed using mixed models with line as a random effect. Models were implemented in R using the lmer function of the lme4 package.

At each generation during the 10 generations of selection, we recorded the survival time of the first three dead flies in each vial. We could not wait until longevity was recorded to choose individuals that gave birth to the next generation, so we decided to make measurements on the same number of flies for each line to limit bias. We assessed the significance of the selection effect by comparing the null model [survival_time∼(1|line)] with the complete model [survival_time∼generation+(1|line)] (likelihood ratio test, d.f.=1). Next, we analysed intra-line behaviours using linear models (survival_time∼generation). The strength of the selection effect was estimated by the corresponding slope.

We recorded CCRT from generation 1 to 10 for the selected lines. Contrary to our expectations, CCRT increased instead of decreasing. Therefore, we started to record CCRT in control lines from generation 6, and, for unknown reasons, we noticed that CCRT also increased in control lines. To test the existence of a response to selection, we only used data obtained from generations 6 to 10, assessing the significance of the treatment by generation interaction by comparing the complete model [recovery_time∼treatment×generation+(generation|line)] with the model without an interaction term [recovery_time∼treatment+generation+(generation|line)] (likelihood ratio test, d.f.=1). A negative interaction term for the selection treatment is expected if recovery time increases less in the selection condition compared with the control. We also analysed intra-line behaviours using linear models (recovery_time∼treatment×generation). The effect of the treatment by generation interaction corresponds to the difference between the respective slopes for the selection and control treatments.

We assessed the significance of the selection effect by comparing the null model [expression∼(1|line)] with the complete model [expression∼treatment+(1|line)] (likelihood ratio test, d.f.=1). The corresponding P-values are provided in Fig. 3. We determined per-line global patterns of variation between selection and control conditions using paired Wilcoxon tests on gene results.

We assessed the significance of a global effect of the selection procedure by comparing the null model <FA or fitness>∼(1|line) with the complete model <FA or fitness>∼treatment+(1|line) (likelihood ratio test, d.f.=1). At the intra-line level, we compared the selection and control treatments using t-tests for FA and Wilcoxon tests for fitness estimates.

We used three inbred lines of D. melanogaster and one outbred line to test for selection response in two traits: reduced CCRT and starvation resistance (2759 and 3881 flies analysed, respectively). We used the outbred line as a positive control for response to selection as it is known that the traits of interest display a genetic basis (Falconer and Mackay, 1996). As inbred lines display non-genetic variability, whereas the outbred line displays both genetic and non-genetic variability, response to selection was expected to be stronger in the outbred line than in the other lines.

Globally, we found weak, although significant, responses to our selection protocols. In the starvation assay, survival time significantly increased over time (χ²=15.5, P=8×10⁻⁵; generation effect=0.44). In the CCRT assay, the generation by treatment interaction was significant (χ²=30.5, P=3×10⁻⁸), with a negative effect (−2.5) due to selection, indicative of a relative shortening of recovery time along generations. However, in both selection experiments, we noticed a high inter-line variability, which we describe in more detail below.

As expected, the strongest response to selection for starvation resistance was observed for the outbred line (slope=1.09, P=2×10⁻⁴) (Fig. 1A). Line 10.1 also responded to selection (slope=0.63, P=0.003). However, no significant response to selection could be detected in lines 6.6 (slope=−0.07, P=0.75) and 15.4 (slope=0.10, P=0.48).

We observed that CCRT increased significantly less across generations in the selected line than in the control for the outbred line and 10.1 (slope difference=−2.89, P=1×10⁻⁶ and slope difference=−5.78, P<10⁻⁷, respectively). The effect was not significant for line 15.4 (slope difference=0.42, P=0.62). Missing data in line 6.6 prevented us from comparing slopes regarding the same generations (Fig. 1B). We also noticed a high variability of the response across and within generations. Such oscillations are frequently observed in experimental selection protocols (Falconer and Mackay, 1996).

Strikingly, we noticed an absolute decrease of CCRT only in the selected outbred line. For both traits (ST and CCRT), we noticed that the outbred line displayed a strong response to selection, as expected. Inbred line 10.1 also showed some response to both selection processes. In contrast, line 15.4 did not significantly respond to selection for either of the traits. Line 6.6 was insensitive to starvation resistance selection. This illustrates the strong line effect in response to selection in inbred lines.

Broad heritability values of cumulative realized heritability (ΣR) for CCRT and ST were estimated per line (Table 1).

As expected, the outbred line presented the highest heritability estimates (0.20 for CCRT and 0.16 for ST). The heritability estimates for CCRT and ST were low and consistent with those estimated for most physiological or behavioural traits in outbred populations (Mousseau and Roff, 1987; Roff, 1997).

We analysed 30 females per line in both selection experiments and recorded the number of flies with visible morphological alterations (including deformed scutellar bristles). We did not detect any increase in the number of aberrant phenotypes in selected lines compared with control lines (Fisher's exact tests; Table S3). We also recorded the number of scutellar bristles and considered as an aberrant phenotype any number different from four. We did not detect any increase in the number of aberrant scutellar bristles following the selection process, except in the case of line 6.6 for CCRT selection (Fisher's exact tests; Table S3).

We could not detect any effect of selection on FA using the complete dataset (CCRT assay: χ²=0.99, P=0.319; ST assay: χ²=1.12, P=0.291). However, at the intra-line level, we detected a significant increase in FA of sternopleural bristle numbers for line 6.6 for both selection experiments (t-tests; CCRT: P=0.036; ST: P=0.031) (Fig. 2).

Fitness was not altered after the selection experiments (CCRT assay: χ²=1.18, P=0.278; ST assay: χ²=0.02, P=0.890), except for a significant decrease in line 6.6 selected for reduced CCRT (Wilcoxon test, P=0.03; Fig. 2B) and in line 15.4 selected for ST resistance (Wilcoxon test, P=0.03) (Fig. 2A).

We quantified gene expression for a set of genes that could be involved in response to selection, genome stability or both. Our set included genes associated with stress, epigenetics and the Piwi-interacting RNA (piRNA) pathway. These genes were chosen because a previous experiment showed that their expression levels were modified after chill coma and starvation stress in D. melanogaster and Drosophila simulans species (B.F.M., J.S.-O., M.F. and C.V., unpublished data). Colinet et al. (2010) also showed an increase in the expression Hsp27 and other Hsp genes after a cold shock. We found that the selection process had an effect on the expression levels of some genes: thor in the CCRT assay (P=0.025), and thor (P=0.007), Hsp27 (P=0.017), modmdg4 (P=0.029) and ago3 (P=0.016) in the starvation assay (Fig. 3). Additionally, line 6.6 showed a significant global decrease in gene expression levels for cold treatment (Wilcoxon signed-rank test, V=135, P<10⁻⁴), while line 10.1 and the outbred line showed an overall upregulation (Wilcoxon signed-rank tests, V=16, P=0.005, and V=12, P=0.002, respectively).

We sequenced a portion of the first intron of the Adh gene to assess genetic diversity within lines. All sequences retrieved from line 6.6 (16 sequences) and line 15.4 (15 sequences) were 100% identical. Among the sequences retrieved from line 10.1 (17 sequences), we found two variants that differed by five single-nucleotide polymorphisms (SNPs) and a short indel. Among the sequences retrieved from the outbred line, only eight were analysable because of a high proportion of heterozygotes, and these matched the same two variants as recorded in line 10.1. Such genetic diversity is expected in the outbred line.

In this study, we performed artificial selection experiments on inbred lines of D. melanogaster to affect resistance to chill exposure and nutrient restriction. In the absence of genetic variation, we expected selection to be inefficient, unless non-genetic variability can also be inherited. As far as we know, this is the first study of selection on traits performed on inbred lines of insects of a natural origin. Apart from this study, and using D. melanogaster transgenic lines, Ciabrelli et al. (2017) recently successfully performed divergent selection for eye colour determined by alternative epialleles, as defined by differential levels of H3K27me3 on a mini-white construct.

We observed responses to selection only in lines for which we cannot exclude genetic variation. This prevents us from concluding that we managed to select non-genetic variation. As expected, the response to selection was more intense in the case of the outbred line than in any one of the inbred lines. Heritability estimates were low and consistent with those estimated for most physiological or behavioural traits in natural populations (Mousseau and Roff, 1987; Scheiner and Callahan, 1999; Scheiner et al., 2000; Gerken et al., 2016). Such low estimates were expected in inbred lines.

We established that the selection processes had an effect on gene expression for some genes involved in the stress response, epigenetics and TE control. It has already been shown that direct and indirect responses to selection can affect several sets of genes with different pleiotropic effects (Mackay, 2014). In the present study, we detected upregulation of modmdg4 and ago3 in the CCRT assay, downregulation of thor in both assays, and downregulation of Hsp27 in the starvation assay.

We noticed variability across lines in the intensity of the response to selection, which we propose is partly related to the extent and nature of the epigenetic variability specific to each line. Such line variability was also observed regarding gene expression measures: the selection process induced changes in expression profiles; however, these changes were line dependent, as expected for independent inbred lines (England et al., 2003).

In addition, we noticed a large variability of responses across and within generations (Fig. 2). Such oscillations are frequently observed in experimental selection protocols and are particularly observed in populations with low genetic load (Falconer and Mackay, 1996), which is probably not the case for the inbred lines. Morgante et al. (2015) showed that a large phenotypic plasticity existed within lines from a natural population of D. melanogaster, which they called micro-environmental plasticity. Our data suggest that such within-line variability exists, which could be due to variability in the non-genetic component of phenotypic plasticity, but we cannot exclude residual genetic variability. Indeed, despite a large number of sibling mating crosses, inbred lines probably carry residual polymorphism. Ciabrelli et al. (2017) performed deep-sequencing of the genomes of their inbred lines and reported hundreds of thousands of polymorphisms in each line. However, based on sequence analyses, they claim that these differences do not account for the observed phenotypic differences. Here, the data that we have on a limited genetic loci suggest that genetic variability may be somehow higher for line 10.1 than for lines 6.6 and 15.4. Thus, we cannot exclude the possibility that the larger response to selection observed in 10.1 compared with the other inbred lines is related to a higher level of residual genetic variability.

Several studies have described stress-induced variation through natural and artificial selection (Badyaev, 2005), but very few have suggested that artificial selection could be considered a stress (see Belyaev, 1979; Belyaev et al., 1985; Trut, 1999, for selection on tame behaviour in the silver fox). In ongoing work in our group, we observed that selection experiments in Drosophila populations presented occurrences of aberrant phenotypes, significant FA index and alterations of canalized traits such as the number of scutellar bristles (B.F.M., unpublished data). Therefore, artificial selection may lead to fitness decreases and impairment of buffering mechanisms. We suspect alteration in the buffering mechanism following the selection protocol only in strains 6.6 and 15.4. In line 6.6, both selection procedures led to an increase in the number of aberrant scutellar bristles, an increase in FA levels and a decrease in fitness. The fitness of line 15.4 decreased following starvation selection. No morphological alterations were detected for both 10.1 and the outbred line. One possible explanation is that the potential of response to selection, regardless of its mechanism, did not reach its plateau for both traits. Therefore, the environmental perturbations applied did not overcome the buffering mechanisms, while the cryptic genetic variation revealed by strain 6.6 indicates a possible rupture in some canalization processes (Dworkin, 2005).

Indeed, we noticed that lines displaying the weakest response to selection were also those that had impaired buffering mechanisms and fitness decreases. This was mostly line 6.6 but also line 15.4 to a more limited extent. In contrast, line 10.1 responded to selection and showed no buffering impairment nor fitness decrease. In addition, line 6.6 displayed an overall decrease in the expression levels of most studied genes following the selection process, while line 10.1 mostly upregulated the studied genes. Together, all these elements may indicate that line 6.6 is unable to trigger stress response pathways, which results in an inability to respond to selection and deleterious buffering impairments. These data should be considered as possible clues, as we did not measure the expression of genes – namely, stress response genes – after a stress stimulus.

We hypothesized that selection could act upon non-genetic inheritance (e.g. upon an epigenetic methylation pattern or a chromatin structure), which introduces a conceptual similarity between non-genetic and genetic inheritance. Non-genetic variability could arise randomly (e.g. as an epimutation) and subsequently be exposed to selection so that it follows similar dynamics to those of ordinary genetic variants. Many studies address questions about the prevalence of non-genetic effects in natural conditions and their inheritance mechanisms, such as epigenetics (Cubas et al., 1999; Vaughn et al., 2007; Bossdorf et al., 2008; Bossdorf and Zhang, 2011; Zhang et al., 2013; Cortijo et al., 2014; Tricker, 2015; Heard and Martienssen, 2014; Ciabrelli et al., 2017). Contrary to what is known from plants, our results demonstrate that selection of non-genetic variation is not straightforward in Drosophila. One inbred line showed a response to starvation selection, but this line appears to have residual standing genetic variation. We observed a line-dependent weak response to selection, accompanied by some changes in gene expression and buffering mechanisms. For example, the expression level of ago3, a gene involved in TE silencing through the piRNA pathway, was affected by the selection process. We speculate that these differences were maintained after the release of selection, which would indicate a transmission across generations. The transmission of epigenetic marks could explain the final phenotype after selection.

We note that, to accurately disentangle genetic from non-genetic compounds on the phenotypes of lines subjected to selection, we must be able to strictly control genetic variance. As seen with our data, this approach is nearly impossible. In future experiments, a knowledge of the genomic sequences before and after selection will be fundamental. With this work, we hope to provide some clues on the difficulties of clearly demonstrating the ability to select for non-genetic variability. The current advances in high throughput sequencing should help us to delve further into these issues.

We thank E. Desouhant, I. Amat and P. Gibert for helpful discussion.