Stage-Specific Effects of Candidate Heterochronic Genes on Variation in Developmental Time along an Altitudinal Cline of Drosophila melanogaster

Background Previously, we have shown there is clinal variation for egg-to-adult developmental time along geographic gradients in Drosophila melanogaster. Further, we also have identified mutations in genes involved in metabolic and neurogenic pathways that affect development time (heterochronic genes). However, we do not know whether these loci affect variation in developmental time in natural populations. Methodology/Principal Findings Here, we constructed second chromosome substitution lines from natural populations of Drosophila melanogaster from an altitudinal cline, and measured egg-adult development time for each line. We found not only a large amount of genetic variation for developmental time, but also positive associations of the development time with thermal amplitude and altitude. We performed genetic complementation tests using substitution lines with the longest and shortest developmental times and heterochronic mutations. We identified segregating variation for neurogenic and metabolic genes that largely affected the duration of the larval stages but had no impact on the timing of metamorphosis. Conclusions/Significance Altitudinal clinal variation in developmental time for natural chromosome substitution lines provides a unique opportunity to dissect the response of heterochronic genes to environmental gradients. Ontogenetic stage-specific variation in invected, mastermind, cricklet and CG14591 may affect natural variation in development time and thermal evolution.


Introduction
Understanding the ontogenetic trajectories of life-history traits is necessary to explain and predict constraints and variation in the evolution of characters. In this context, any variation in the state of individuals prior to reproduction will be preserved to some degree at the reproductive stage. In this sense, the time elapsed from the embryo to the reproductive phase, commonly known as developmental time (DT), is directly related to an individual's reproductive success [1].
In holometabolous insects like fruit flies, which occupy ephemeral habitats, the impact of DT on fitness is further exaggerated [2]. Although life history models assume that natural selection should favor the maximization of growth rate, genetic correlations with pre-adult viability [2,3], adult body size [4], longevity [5] and starvation resistance [6] constraint the evolution of DT in the model organism D. melanogaster.
Changes in the timing of developmental processes -heterochrony -could account for many evolutionary changes among populations and species. However, the first criterion for heterochonic evolution is the existence of natural variation for developmental time. Several studies have shown latitudinal clinal variation for DT in Drosophila melanogaster populations of Australia [7] and South America [8,9]. Also, Folguera et al. [8] reported that egg-to-adult DT was longer in southern localities and lowland populations, a pattern that leads to the hypothesis of adaptive thermal evolution.
DT is a quantitative trait determined by multiple segregating genes that are sensitive to temperature variation [10]. It has been also shown that heterochronic genes are mainly involved in metabolic and neurogenic pathways, such as Hippo, Notch and insulin signaling pathway [10,11]. However, it is largely unknown whether these pathways affect variation in DT in natural populations. The resolution of this question may uncover the genetic architecture of this ecologically relevant trait and, hopefully, would provide insights into the evolutionary processes governing natural variation. Interestingly, In(2L)t and In(2R)R, two cosmopolitan polymorphic inversions of the second chromosome, affect DT [9,12], and several mutants showing the strongest effects on DT mapped to the second chromosome [10]. Here, we assessed the effects of four genes mapping to the right arm of the second chromosome and that affect metabolism and neurogenesis on natural variation of DT. invected and mastermind are involved in the development of the nervous system [13,14]. cricklet is a gene encoding a carboxylesterase involved in the metabolism of the juvenile hormone [15]. CG14591 has no clear association with any biological process.
Genetic complementation testing is a potent tool to investigate the contribution of individual genes to natural variation in quantitative complex traits [16]. Furthermore, this experimental approach is particularly accurate when employing the chromosome substitution line technique, since natural chromosomes can be extracted from wild organisms and placed into a common genetic background. In this sense, Drosophila brings an impressive toolkit for dissecting multiple interacting loci with individually small effects that affect quantitative developmental traits [16].
We report the results of a survey of variation in DT using second chromosome substitution lines derived from sampling localities that lay along latitudinal and altitudinal gradients from several populations of Drosophila melanogaster of western Argentina. We show that variation in DT among populations is not only correlated with altitude but also with certain features of the thermal regime. Clinal variation in DT over short geographic distances and its association with climatic variables suggests that among-population differentiation is the consequence of different selective pressures along the altitudinal gradient. Also we show that allelic variants at candidate heterochronic genes invected, mastermind, criclket and CG14591 contribute to altitudinal variation in DT with stage-specific effects.

Generating second chromosome substitution lines
Isofemale lines were founded by rearing the progeny of gravid Drosophila melanogaster females collected in six localities along latitudinal and altitudinal gradients in Western Argentina (February 2004 andFebruary 2005). The geographical location, latitude, longitude, altitude and climatic information for each population are given in Table 1. All lines were maintained by full-sib mating for 10 generations on cornmeal-molasses-agar medium under standard conditions of 2561uC, 70% humidity and a 12-h light: 12-h dark cycle. After 10 generations, a single second chromosome was extracted from each isofemale line and substituted into the genetic background of an isogenic Canton-S B strain by standard techniques using balancer chromosomes ( Figure 1). To construct second chromosome substitution lines, one male of each isofemale line were crossed to w; Cy/Canton-S B; Sb/Canton-SB females [generation 1 (G1)]. A single w; Cy/+ 2; Sb/+ 3 male from the progeny of each cross was crossed to w; Cy/Sp; Canton-S B females. Next, w; Cy/+ 2; Sb/Canton-S B males were crossed to w; Cy/Sp; Canton-S B females (G3). Females and males of genotype w; Cy/+ 2; Sb/Canton-S B were intercrossed at G4, and the Cy and the Sb balancers were eliminated in the next generation (G5), obtaining an isogenic second chromosome substitution line with genotype w; +2; Canton S B. By means of this protocol we generated 50 second chromosome substitution lines isogenic for one wild derived chromosome in an otherwise isogenic background common to all lines. An average of 9 lines per population was obtained.

Developmental time assay
For each substitution line, 300 pairs of sexually mature flies were placed in egg collecting chambers for 8 hours. Eggs were allowed to hatch and batches of 30 first-instar larvae were transferred to culture vials containing 10 ml of cornmeal-agar-molasses medium (4 replicate vials per line). Emerged flies from each vial were collected every 12 hours and sorted by sex. We estimated DT as the time elapsed since the transfer of first-instar larvae to the vials until adult emergence. All vials were kept in an incubator at 2561uC, under a 12-h light: 12-h dark cycle and at 70% humidity.
Forward stepwise multiple regression analyses were performed using geographic and climatic variables separately to test for clinal variation in developmental time. The geographic variables considered were altitude and latitude; and the climatic variables were minimal mean temperature, maximal mean temperature, mean temperature and thermal amplitude.

Quantitative complementation tests
Quantitative complementation tests were performed to examine the contribution of invected, cricklet, CG14591 and mastermind to natural variation in DT among the most divergent second chromosome substitution lines. Thus, from our set of isogenic lines we selected those exhibiting the most extreme phenotypes for DT (the fast and slow developing lines). These lines were crossed individually with lines carrying P [GT1] element mutated alleles in one of the four candidate genes identified by Mensch et al. (2008) [10] and, in parallel, with a P [GT1] element-free insertion line with the same genetic background (Canton-S B). The genotypes of the F1 progeny of the crosses were m/+ i and Canton-S B/+ i , respectively, where m is a mutant allele in one of the candidate genes derived from a P [GT1] element insertion line and + i represents a wild derived allele of the candidate gene. We used P [GT1] insertion lines BG00846, BG01339, BG01672 and BG01902 that carry mutant alleles for invected, cricklet, CG14591 and mastermind, respectively [10]. We evaluated variation in developmental time among genotypes by means of a four-way analysis of variance (ANOVA) according to the mixed model: where L, G and S are the fixed cross-classified effects of line (second chromosome substitution line), genotype (Canton-S B or m) and sex, respectively. R stands for the among replicate effect (random), and E is the error (within vial variance). The first criterion that must be met for quantitative failure of complementation is a significant Line by Genotype interaction. However, quantitative failure of complementation can be explained either as due to allelism or epistatic interactions [17,18]. Thus, in order to identify allelic failure to complement, the second criterion to be met is that variation among lines in the Canton-S B background is not greater than variance among lines in the mutant background [17,19]. If both criteria are met, it can be inferred that a candidate gene affects natural variation in developmental time.
In all complementation tests, DT was partitioned into larval developmental time (LDT) and pupal developmental time (PDT), involving the duration of all events that occurred before and after pupation, respectively. All statistical tests were performed using the STATISTICA package (StatSoft, Inc. 1999, 2001).

Patterns of genetic variation among natural substitution lines
A large amount of variation in DT was found among natural substitution lines, indicating that D. melanogaster second chromo-some harbors genetic variation in natural populations for this fitness related trait ( Figure 2). Among the lines showing the fastest phenotypes were Güemes fast (GUE F, 232 hours) and Lavalle fast (LAV F, 233 hours), whereas the lines that took more time to reach the adult stage were San Blas slow (SB S, 339 hours) and Uspallata slow (USP S, 353 hours).
Regression of DT on geographic variables revealed a significant and positive association between developmental time and altitude (r = 0.5, p,0.0001) ( Figure 3) and none with latitude (p.0.05). Also, we detected a positive and significant association between DT and thermal amplitude (r = 0.26, p,0.0001) ( Figure 4). In addition, altitude and thermal amplitude are positively correlated in the array of populations analyzed (r = 0.83, p,0.05) implying that greater daily thermal amplitude is likely to be found in highland populations. This pattern of genetic variation prompted us to explicitly test the hypothesis that altitude may be the main factor accounting for among population differentiation. Thus, we divided the complete set of populations sampled in two groups. The first group included sampling sites situated above 1000 meters above sea level: Uspallata, San Blas and Chilecito, which we to refer as highland populations from here on. The second group includes three lowland localities: Güemes, Lavalle and Neuquén which are located below 800 meters above sea level. Comparisons between these two groups showed that substitution lines from highland populations developed slower than lowland populations in males (F 1,185 = 6.59; p = 0.01) and females (F 1,189 = 7.72; p = 0.006). Indeed, DT substitution lines derived from highland and lowland sampling sites differed, on average, by 12 hours in both sexes.
In order to quantify within population genetic variation, we performed an ANOVA for each sampling locality (Table 2). In all populations, the Line factor was highly significant implying that the second chromosome harbors genetic variation for developmental time in all natural populations. Females developed significantly faster than males in Neuquén, Chilecito and Uspallata populations. The Line by Sex interaction term was not significant in any population.

Candidate genes for altitudinal variation in developmental time
Since we detected great genetic differentiation between highland and lowland populations, we decided to test whether variation in candidate genes previously identified [10] could account, at least partially, for the clinal pattern. To this end we selected substitution lines with the most divergent DT for four independent quantitative complementation tests (QCTs) employing lines bearing heterochronic mutant alleles in a set of four genes, mastermind, invected, cricklet and CG14591, that mapped in the second chromosome. Thus, natural genotypes that differed in DT by more than five days were used in four independent QCTs. Our  results showed a significant line by genotype interaction in each QCT suggesting a failure of complementation in all four candidate genes (Table 3). Likewise, in each analyses, L 6 G 6 S terms were not significant (Table 3) implying that genetic variation affects both sexes in similar way. In addition, in each QCT, variance among wild chromosomes was similar between mutant and Canton-S B backgrounds (Table 4), meaning that failure to complementation can mainly be attributed to allelism, rather than epistasis, as the genetic mechanism generating the phenotypic differences among lines.
Genetic variation in candidate genes was revealed by the different responses of natural chromosomes when combined with mutant alleles or the functional Canton-S B allele. Effects on DT of different wild chromosomes over mutant and Canton SB backgrounds are expressed as ratios for each QCT ( Figure 5). Most ratios showed positive values indicating that wild chromosomes had longer DT on a mutant background than on Canton-S B, the only exception being line SB S in the QCTs for invected, mastermind and CG14591. Nevertheless, in the QCT for CG14591, we also identified a line (GÜ E F) with a DT ratio close to zero, suggesting that this chromosome did not contribute to natural genetic variation for this gene. All in all, our results not only confirm that mastermind, invected, cricklet and CG14591 are quantitative candidate genes for developmental time expression, but also that these genes exhibit natural allelic variants that contribute to natural variation in DT.

Stage-specific variation in developmental time
In all QCTs, we further dissected egg-to-adult developmental time by measuring the duration of the larval (LDT) and pupal (PDT) stages. For all candidate genes, we found that the L 6 G interaction was significant for LDT but not for PDT (Table 5 and  Table 6, respectively), indicating that segregating genetic variation mostly affected the time spent in premetamorphic stages rather than the duration of metamorphosis itself. As it can be observed in Figure 5, LDT accounts for the largest proportion of egg-to-adult  DT genetic variation in these lines (blue and red bars, respectively). Once again, variance among wild chromosomes was similar between mutant and Canton-S B backgrounds (Table 4), supporting the idea that allelism is the best explanation for the genetic variation observed rather than epistasis.

Variation in DT in second chromosomes and candidate genes
Developmental time is a fundamental life history trait closely linked to fitness, a short juvenile period is beneficial since it reduces the risks of predation and/or infections and has the advantage of earlier reproduction [1]. Life history models assume that natural selection should favor the maximization of growth rate. However, genetic variation within and between populations has been frequently identified for developmental time, suggesting that DT in a particular environment is not necessarily the maximum and that rapid growth rates may have an associated fitness cost.
Our survey of DT variation in substitution lines bearing different wild second chromosomes, in an otherwise identical genetic background, revealed a considerable amount of variation. The lines exhibiting the most extreme scores for DT differed by as much as 5 days. Moreover, our survey shows that chromosome 2 harbors an important portion of natural genetic variation affecting DT in all populations analyzed. In fact, the contribution of differences among chromosome 2 substitution lines to egg-to-adult DT variation ranged from 59% to 87% of the total phenotypic variance in the set of populations studied ( Table 2).
Also interesting is that the line by sex interaction did not contribute to variation in DT in any of the localities, suggesting the lack of genetic variation for sexual dimorphism of DT. Interestingly, Mensch et al. (2008) [10] demonstrated a null input of the genotype by sex interaction to total phenotypic variance using a set of lab lines, as did the study of Van't Land et al.(1999) [9] showing clinal variation in South American populations. The absence of genetic variation for sexual dimorphism in developmental time has often been observed in other drosophilids as well [20,21].
Regression analyses of DT on geographic and climatic variables indicated that altitudinal clines and the association with thermal amplitude reported in Folguera et al. (2008) [8] could be partially accounted for by allelic variants in genes located in the second chromosome. The associations between DT and altitude and thermal amplitude are congruent since daily mean thermal amplitude is likely to be more common in highland than in lowland sampling sites [8]. Further analysis comparing highland (after grouping localities situated above 1000 meters above sea level) and lowland (below 1000 meters above sea level) populations showed that DT in highland substitution lines was, on average, 12 hours longer. One of the major challenges to understand the genetic basis of complex adaptive traits (i.e. developmental time) has been the difficulty in mapping individual genes affecting natural genetic variation because differences in genetic background can profoundly affect complex phenotypes [22,23]. Nevertheless, some successful efforts [24][25][26][27] allowed the identification of candidate genes for diverse traits. It is important to note that most of these studies used chromosomes carrying deficiencies [22] that were  induced in heterogeneous collections of genetic backgrounds. In contrast, our study was performed with chromosome substitution lines that share the same genetic background either among them and with P element insertion lines that we used to characterize heterochronic candidate genes [10]. Failure of complementation in QCTs revealed that invected, mastermind, cricklet and CG14591 have natural allelic variants that affect DT. To our knowledge this the first time that natural allelic variation is described for these genes as well as the (partial) characterization of the genetic basis of DT variation in natural populations. Moreover, the similar variances of natural allelic variants when combined with mutant and Canton S B alleles, in all QCTs, suggest an additive effect of allelic variants of mastermind, cricklet, CG14591 and invected on DT. However, non additive effects due to interactions with other second chromosome loci affecting DT cannot be ruled out. Notice, for instance, that the genotype Lavalle/mutant, exhibited faster larval developmental time than the Uspallata/mutant genotype, when the mutant allele is contributed by the P-element insertion line for invected ( Figure 6). On the other hand, the Lavalle chromosome prolonged and Uspallata chromosome speeded up development when combined in genotypes with mutant alleles contributed by the P-element insertion lines for mastermind, cricklet and CG14591 ( Figure 6). These results suggest that complex genetic interactions among several developmental pathways are affecting DT. Noteworthy, two of the genes, invected and mastermind are expressed early in development and also during the larval and pupal stages [14,28]. The former, a paralog of the pair-rule gene engrailed, is a transcription factor that is expressed in neuroectodermic cells [13], while mastermind is a key component of the Notch signaling pathway, which determines cell fate and regulates pattern formation. Likewise, cricklet (CG9858, also known as EstR) encodes a carboxylesterase that has been proposed to be associated with Juvenile hormone (JH) functions [15]. Yet, recent studies have shown poor kinetics towards various forms of JH [29]. Finally, no clear association with any biological process has been reported for CG14591. Our study provides, therefore, the first record concerning functional significance of CG14591 in D. melanogaster.

Holometabolous ontogeny and the nature of variation in DT
As a complex trait, DT displays considerable genetic correlations with other life history and morphological characters.
Particularly, it has been suggested that DT may be involved in a trade-off with body size, since attaining a large size may imply a longer feeding period, and thus a longer DT. Thus, variation in growth patterns, growth rate and the duration of the growth period, affect the age and size at maturity, implying that the study of juvenile growth may be crucial to understand life history evolution [30]. In holometabolous insects, body size and development time are controlled by three factors: i) growth rate, ii) the time required to attain a critical size (related to the cessation of secretion of Juvenile hormone) and iii) the interval to cessation of growth (related to the timing of ecdysteroid secretion) which corresponds to the onset of pupation and metamorphosis [31]. Early in the last larval stage, larvae reach a critical size (at which point they have acquired enough nutrients to complete development), stop feeding and search for a pupation site. However, the cessation of feeding does not occur immediately, and there is a delay between the attainment of critical size and pupation. Thus, total developmental time can be conceived as the sum of the duration of successive discrete phases: embryonic, larval and pupal stages that end with the emergence of the imago. Our results illustrate that allelic variation in the four candidate genes studied only affected the time interval between the egg stage and pupation, which we call Larval Developmental Time (LDT) but not the interval between the pupal stage to adult emergence, Pupal Developmental Time (PDT). This contrasting pattern of variation along the ontogeny provides an initial insight to understand the genetic architecture underlying developmental time. In this sense, our results concur with the idea that variation in LDT may no have direct effect on the duration of the pupal stage [32]. Instead, since resource acquisition and growth are confined to the larval stages, LDT variation might directly affect adult body size [33]. This might be the case of invected alleles that not only affected larval developmental time but also influenced various adult body size traits [34].
The recognition that a single trait, DT, may be subdivided in ''subtraits'' that are not under the control of the same genes imply that the duration of larval and pupal stages responded to different selective pressures; a matter that may have broad implications in evolutionary genetics, since responses to selection may crucially depend on the genetic variance/covariance of individual ''subtraits'' [32]. In this context, thermal selection is usually invoked as a causal factor for the frequently reported geographic clines of body size in Drosophila [7,8,11] however whether thermal selective pressure exerts a direct influence on LDT (reversed Bergman rule) or, actually, the target of selection is adult body size itself (Bergman rule) remains unclear. In this sense, a recent study conducted in our lab has shown that LDT and PDT respond differentially to thermal variation (Folguera et al., unpublished results). Noteworthy, mastermind, cricklet and CG14591 were also shown to be candidates for the plastic responses of DT to temperature variation [10], meaning that the effects of each gene depend on ambient temperature along ontogeny. All in all, thermal selection along the altitudinal cline may influence the large genetic variation in metabolic and neurogenic genes that affect LDT. In summary, our results are consistent with the hypothesis that egg-to-adult development time should not be considered as the sum of additive ''sub-traits''. Indeed, as development is a complex time-dependent process, a system genetics approach would be appropriate in order to test the hypothesis of stage-specific developmental time variation at the genomic level. To this end, our study provides particular candidate genetic pathways to be tested in future studies taking advantage of the abundant genomic information available in D. melanogaster and allied species [35,36].