Genome Wide Association Studies of early �tness traits in Drosophila melanogaster unveil plasticity and decoupling of different aspects of phenotype

Great efforts have been sustained to explain the relationships between genotype and phenotype for developmental �tness traits through the study of their genetic architecture. However, crucial aspects of functional architecture in�uencing the maintenance of genetic variability, and thus the capacity for evolutionary change, are still unexplored. Here we performed Genome-wide Association Studies for phenotypic variability, plasticity and within-line canalization at two temperatures for Larval Developmental Time (LDT), Pupal Developmental Time (PDT), Larval Viability (LV), Pupal Viability (PV), and Pupal Height (PH) in lines derived from a natural population of Drosophila melanogaster. Results suggest changes in genetic networks linked to resource acquisition and allocation underlying variability for all traits. However, we found low genetic pleiotropy between traits and for different aspects of phenotype (means, plasticity, within-line canalization) within each trait. Their genetic bases were also temperature-specic: we found no variants showing an effect for the same trait at both temperatures. Moreover, a genetic decoupling between larval and pupal traits was con�rmed, as there were no candidate variants signi�cantly associated to phenotypic variability for the same trait across stages. We found evidence of genetic antagonistic pleiotropy for several loci affecting larval and pupal traits. The high degree of modularity at various levels would allow for the independent evolution of distinct aspects of the phenotype in different environments and ontogenetic stages. This may explain why genetic variation for these adaptive traits is not extinguished in natural populations and may entail important implications regarding evolvability.


Introduction
Understanding the genetic, epigenetic and environmental factors affecting natural variation in life history traits is one of the major goals of evolutionary biology (Stearns 1992, Roff 2002, Goldberg et al. 2007, Flatt 2011).Indeed, the study of the genetic bases underlying phenotypic variation for adaptive traits and their responses to environmental conditions (phenotypic plasticity) allows one to elucidate the phenotype-genotype map through the analysis of their genetic architecture (Mackay 2001, Houle et al. 2010, Goddard et al. 2016).Moreover, characterizing the change in developmental systems that underlie adaptive traits and their con gurations in different environments would allow the comprehension of evolutionary dynamics beyond the gene level, wherein phenomena such as pleiotropy, epistatic interactions, phenotypic plasticity and canalization play a major role (Hansen 2006).
Most studies on evolutionary biology have focused on late tness components and when early tness components have been addressed, it has been mostly with the intention to investigate trade-off relationships with adult traits (Roff 2002, Flatt 2011).Certainly, few studies have investigated the genetic basis of phenotypic variation in early tness traits (Lavagnino et al. 2008, Mensch et al. 2010, Petino Zappala et al. 2018, 2019).However, the assessment of developmental adaptive traits has been characterized as fundamental for the understanding of evolutionary dynamics (Hansen 2006, Müller 2007), as the focus of evolutionary developmental biology lies in the generation and maintenance of variability affecting the biological processes that impact reproductive success, through the characterization of developmental pathways and their changes during ontogeny and phylogeny (Raff 2000, Hall 2003).
Drosophila melanogaster represents an excellent model to investigate the genetic basis of phenotypic variation for early tness traits.The species presents a complex life cycle that progresses through immature (embryo, larvae and pupae) to mature (adult) stages and includes complete metamorphosis.The probability of reaching the next life cycle stage (i.e., viability) is affected by genetic and several non-genetic factors such as density (larval crowding), nutritional status, temperature, presence of pathogens or predators (Barker and Podger 1970, Kolss et al. 2009, Trotta et al. 2016, Petino Zappala et al. 2019).The sequence from one ontogenetic stage to the next is temporally organized, physiologically controlled and dependent on environmental factors, such as nutrient availability, temperature, and the changes in daylight, which in turn may affect the timing of certain developmental events (Mirth et al. 2005(Mirth et al. , 2014; Mensch et al. 2008Mensch et al. , 2010, Di Cara and King-Jones 2013, Nijhout et al. 2014, Petino Zappala et al. 2018).Changes in the time elapsed between the embryonic and adult stage (commonly known as developmental time [DT]) is thought to have a great impact on tness on a species that forages on decaying substrates.In this sense, a faster development would correlate with a higher preadult survival (Partridge and Fowler 1992, Nylin and Gotthard 1998, Narasimha et al. 2015).On the other hand, the reduction of adult tness as a consequence of the acceleration of DT, in a pattern called Fast Development Syndrome, illuminates the direct connection of the preadult and adult stages through energetic trade-offs (Chippindale et al. 2004, Yadav andSharma 2014).Preadult traits such as DT are not usually decomposed in their larval and pupal components, despite the fact that larval and pupal stages are characterized by different behaviors, developmental processes and genetic pathway activities (Arbeitman et al. 2002).Previous work has demonstrated that Larval and Pupal DT (LDT and PDT) can be decoupled; this heterochronic relationship shows a genetic basis and may respond to selection (Mensch et al. 2010, Petino Zappala et al. 2018).Certainly, collapsing both components into an "uni ed" egg-to-adult DT would represent a loss of valuable information, as it would hinder the understanding of the diverse processes occurring during different life stages.It has also been shown that the consequences of this ontogenetic modularity may affect evolutionary potential in holometabolous insects (Yang 2001).Therefore, here we will consider the genetic basis for DT and Viability separately for larval and pupal stages.This will also allow us to analyze whether this decoupling can be explained by genetic antagonistic pleiotropy (i.e.differences in sign of the effect of each loci at each stage) and/or by loci presenting effects of different magnitude at each ontogenetic stage (Huang et al. 2020).Finally, Pupation height (PH), de ned as the distance from the surface of the rearing medium to the center of the pupa, is usually an indicator of larval motility, adequate wandering and feeding behavior (Sokolowski and Hansell 1983) which is related to phototaxis (Schnebel and Gross eld 1986) and negative gravitaxis (Fowler and Montell 2013).Despite its dependence upon environmental conditions, especially larval crowding (Mueller and Sweet 1986), it is known to be determined by a polygenic basis and is related to other developmental and behavioral traits affecting reproductive success.As in the case of preadult Viability, PH can be affected as a result of the Fast Development Syndrome (Chippindale et al. 1997(Chippindale et al. , 2004)).Studying the trade-offs and the degree of pleiotropy between these adaptive traits is relevant to determine the maintenance of genetic variability and possible restrictions on their evolution (Davidowitz et al. 2004).
While traditionally trait means have been the major focus of evolutionary biology, other aspects of the phenotype may be relevant to the establishment of ecological strategies to deal with environmental heterogeneity and therefore for evolutionary dynamics (Dewitt 2016(Dewitt , Ørsted et al. 2018).Although phenotypic plasticity, the ability of a genotype to produce different phenotypes according to environmental conditions (Schlichting and Pigliucci 1998), may not always respond to an adaptation, it is expected to be selected for traits that present different phenotypic optima across frequently experienced environments (DeWitt and Scheiner 2004, Ayrinhac et al. 2004, Flatt 2005) when the relationship between phenotype and tness is predictable (Simons 2014).Phenotypic plasticity can be coupled to varied degrees of phenotypic robustness within environments (Simons 2014).In relatively homogeneous or predictably changing conditions, it is expected that adaptive traits will become canalized in a population (i.e.phenotypic variability will be restricted around an optimal phenotypic value for such environmental conditions as a result of stabilizing selection).This would occur at least in part through the xation of variants and the establishment of epistatic interactions promoting phenotypic robustness against minor environmental or genetic perturbations (Rice 1998, Flatt 2005, Ørsted et al. 2018).As some genetic variability will remain "unexpressed", this phenomenon allows for the accumulation of genetic variability within a population which is not subject to natural selection.However, an environmental perturbation may cause changes in epistatic interactions, affecting individual robustness and the degree of canalization of the phenotype at the population level through the liberation of this cryptic genetic variation.In dealing with unpredictable and/or heterogeneous environments, developmental instability could be considered as a bet-hedging strategy (Simons and Johnston 1997), maximizing phenotypic diversity between genetically similar or identical individuals so that part of the population is suited to the changing conditions (Herman et al. 2014), which could provide an opportunity for adaptive evolution (Paaby and Gibson 2016).The degree of phenotypic canalization within genotypes (or lines) at each environment can be estimated by means of the Environmental Variation Coe cient (CVE).Different degrees and combinations of phenotypic plasticity and within-line canalization patterns of a trait can be selected allowing populations to deal with the challenges posed by environments with predictable and unpredictable components (DeWitt and Scheiner 2004, Simons 2014).Both properties of a trait emerge from its genetic architecture; and result in the conditional neutrality of some genetic variants, because their phenotypic effect can be concealed from natural selection in some conditions.Thus, these phenomena are relevant for the maintenance of genetic variability, which increases invasiveness and evolvability (i.e. the ability to foster evolutionary change) (Dewitt

2016, Paaby and Gibson 2016).
Environmental temperature is an important factor affecting organismal development in ectotherms (Angilletta et al. 2004, Hoffmann andWeeks 2007).Indeed, adaptations to different temperature regimes, which may be associated to changes in life history traits, have been reported for D. melanogaster populations (Hoffmann and Weeks 2007, Folguera et al. 2008, Trotta et al. 2016).Therefore, assessing phenotypic plasticity and canalization patterns related to thermal variability for this model organism is key to understand adaptations and responses to environmental change.
In this paper, we performed Genome-Wide Association Studies (GWAS) using 39 of the original "core 40" lines derived from the Drosophila Genetic Reference Panel (DGRP) whose genomes are sequenced (Mackay et al. 2012).With these lines, we had previously demonstrated a "decoupling" for DT at the larval and pupal stages (Petino Zappala et al. 2018) and that variation in preadult viability also presents a stage and thermal speci city (Petino Zappala et al. 2019).In this work, we wanted to expand on our previous studies by corroborating the changes in the genetic bases of these traits across stages and temperatures, while also analyzing different aspects of phenotype (i.e., means and within-line canalization).In this sense, we expected to obtain speci c candidate loci underlying phenotypic variability for each trait at each stage and temperature combination.As we considered here the different traits in the same analysis, including their within-line canalization patterns, we were allowed to explore the degree of genetic modularity and pleiotropy associated to each aspect of the phenotype.We consider that decomposing adaptive phenotypes in their different components while also accounting for changes in environmental factors is key to further understand the dynamics of genetic architecture, the maintenance of genetic variability in natural populations and thus, the evolution of evolvability.

Materials And Methods
The core set of 40 DGRP lines were acquired from the Bloomington Drosophila Stock Center.Brie y, they were obtained by 20 generations of full-sib mating from isofemale lines collected from the Raleigh, NC, USA population, and their genome has been fully sequenced (Mackay et al. 2012, Huang et al. 2014).Stocks were maintained under standard culture conditions (cornmeal-dextrose-agar medium, 25°C, 60-75% relative humidity, 12:12hr light:dark cycle).39 of these lines were used for the nal GWA studies.

Experimental design
For each temperature, 300 pairs of sexually mature ies corresponding to each line were placed for 8 hours in separate oviposition chambers.Eggs were allowed to hatch and batches of 30 rst-instar larvae were transferred to culture vials containing a standard cornmeal, yeast and agar culture medium.Larvae were reared at different controlled temperatures: 17°C and 25°C ± 0.5°C, the other culture conditions (i.e., medium, relative humidity, lightdark cycle) being the same at both temperatures.For each of the 39 lines at both temperatures there were 4 replicate vials each containing 30 larvae.Every 12 hours, we recorded formation of new pupae and their position as well as the number of emerged adults.
We quanti ed ve early tness components (Figure 1): Larval Developmental Time (LDT), Pupal Developmental Time (PDT), Larval Viability (LV), Pupal Viability (PV), and Pupal Height (PH) and their CVEs in all DGRP lines reared at both temperatures as previously stated in Petino Zappala et al. (2018,2019).Viability (LV and PV) was estimated as the proportion of individuals surviving each stage.In the case of LV the viability is relative to the number of rst-instar larvae seeded in culture vials (30) whereas PV was calculated as the proportion of emerged adults relative to the number of formed pupae in each replicate vial.LDT was estimated as the time elapsed since the estimated egg-laying time (at half the oviposition interval) until pupae emergence, averaged by replicate vial.PDT was calculated as the difference between average total DT (since the estimated egg-laying time until adult emergence) and average LDT for each replicate vial.Finally, PH was measured with a caliper as the distance between the center of each pupa and the surface of the culturing medium.Replicates were run on different days to randomize environmental variation.
Coe cients of Environmental Variation (CVEs) per line for all traits were estimated as CVE = σ/μ, where μ stands for the line's mean for the trait and σ represents the standard deviation between replicates of each line (Mackay and Lyman, 2005).

Statistical analyses
To study changes in phenotypic canalization, we performed ANOVAs on the CVEs by line for all traits with the model Y = μ + T + E, where T stands for Temperature ( xed factor) and E represents within-temperature variance.Genetic correlation analyses were performed for all early tness components within each temperature and for the same component between temperatures for both means and CVEs.Also, a correlation analysis was performed between the mean and the CVE for each early tness component within each temperature.Means for each Line-Temperature combination were used in all correlation analyses.As we utilized the same values for different correlation analyses, Bonferroni correction for multiple tests was applied (p B =0.0042).

Genome-wide Association studies
We performed the analysis for each early tness components utilizing line means for phenotypic scores and 1453947 polymorphic markers of the DGRP population publicly available in the SNP (Single Nucleotide Polymorphism) data located in dgrp2.gnets.ncsu.edu;DGRP Freeze 2.0 (Mackay et al. 2012, Huang et al. 2014).
Analyses were done by means of the available web application, using as input the list of mean phenotypes for each line and temperature, instead of using male and female phenotypes for each line.This pipeline adjusts line phenotypic means for the effects of Wolbachia infection and 5 major chromosomal inversions (In(2L)t, In(2R)NS, In(3R)P, In(3R)K, In(3R)Mo).The adjusted phenotypes are tted using a linear mixed model Y = μ + M + G + E, where μ is the population mean, M is the xed effect of the polymorphic marker, and G is a polygenic term with its covariance among lines determined by the genomic relationship matrix (Huang et al. 2014).The web application returns a list of SNPs signi cantly associated to phenotypic variability (adjusted and non-adjusted for the effect of Wolbachia and inversions) at each temperature and for the interaction term (thermal plasticity).
Polymorphisms with FDR<0.1 (Benjamini and Hochberg, 1995) were retained as candidates affecting the corresponding trait.In the case that two or more genes lied within 1kb of a candidate variant, only the one closest to it was retained as a candidate gene.The same criterion was applied to regulatory regions associated to candidate variants.

Enrichment analyses
The Flymine platform (Lyne et al. 2007), which gathers data from different databases, was used to analyze Gene Ontology (GO) terms that were enriched in the list of candidate genes, by testing the null hypothesis that the different terms were represented in the list of candidate genes in the same proportion as in the genome.Holm-Bonferroni correction was applied to set the critical p-value and results were normalized according to gene length in the platform settings.

GWAS for Larval-Pupal phenotypic decoupling
Given the decoupling observed between larval and pupal phenotypic variability, a second complimentary set of GWAS was performed within each temperature for both stages of DT and Viability, wherein the signi cant SNPs for the interaction term indicate loci associated to this phenomenon.Analyses were performed through the dgrp2.gnets.ncsu.edu;DGRP Freeze 2.0 web application, in this case with a list of larval and pupal phenotypes for each line within each temperature.These SNP candidates were not included in the enrichment analyses but were analyzed to assess their effect on the particular phenotype for the larval and pupal stages.

Quantitative genetics of early tness components
We reared 9360 rst-instar larvae resulting in 7883 pupae from which we measured LV, LDT and PH, and 6820 adults allowing to evaluate PV and PDT.As previously reported (Petino Zappala et al. 2018, 2019) all traits showed moderate to high heritabilities (Additional File 1).Both mean and CVE parameters were estimated for all traits analyzed at two different rearing temperatures, 17°C and 25°C (Additional File 2, Additional File 3).The comparison of the genetic (Pearson and Spearman) correlation matrices between temperatures for traits' means (Table 1) revealed that correlation patterns were mostly temperature speci c.This and the fact that no correlation between temperatures for each trait was signi cant suggest that the genetic bases involved in phenotypic variation for these traits depend on the rearing temperature.The analyses of these matrices at 17°C and 25°C indicate that only LDT and LV exhibited signi cant and conserved correlations at both temperatures.Moreover, no correlation was found between ontogenetic stages for the same trait (LV/PV and LDT/PDT) at either temperature.Our results also revealed signi cant differences between temperatures for the CVEs of PDT and PV, indicating a decanalization of phenotypes at 17°C in both cases for the genotypes within this population (Figure 2, Table 2).DT traits were in average more robust (showed lower CVEs) than Viability traits and PH.As was found for traits' means, there were differences between CVE patterns at both temperatures (Table 3), indicating that the genetic bases underlying the traits' canalization patterns were not conserved across temperatures for these lines.On the other hand, correlation analyses between means and CVEs showed signi cant and negative values for both Viability traits and PH (Table 3), although it is worth mentioning that in the case of Viability this result could be an artifact caused by the upper limit of the variable.

Genome-wide Association Studies
We performed GWA Studies on LDT, PDT, PH, LV, and PV means and their CVEs by line (Additional File 4) on 39 lines of the DGRP to identify genomic regions harboring variants that contribute to natural phenotypic variation in ies reared at 17°C, 25°C, and to phenotypic plasticity between these two temperatures.We used a mixed model that accounted for any effects of Wolbachia infection and chromosomal inversions (Huang et al. 2014).There were no signi cant associations of any early tness component with Wolbachia infection at either temperature, but 2 of the 5 polymorphic inversions present in the subset of lines affected variability for LDT after Bonferroni correction (Table 4).
We identi ed a total of 924 variants, excluding duplicates (i.e., variants obtained for more than one trait), that were associated with variation on early tness traits, considering both types of variables (i.e., means and CVEs [Table 5, Additional File 4]).The number of variants affecting means and CVEs differed signi cantly among the traits analyzed (χ 2 4 =651.92p<10 -5 ); the result was still signi cant (χ 2 4 =544.96p<10 -5 ) even if we eliminated variants that are in Linkage Disequilibrium (LD) considering a cutoff of 50bp (Mackay et al. 2012), suggesting that the genetic bases affecting traits' means and CVEs are independent.
From the total genetic variants detected, 38% exhibited signi cant results for traits' means and 59% affected CVEs, whereas only 3% was simultaneously signi cant for traits' means and CVEs.Interestingly, 98.9% of SNPs or InDels affecting traits' means were temperature-speci c, wherein a majority of the variants (221 polymorphisms, representing 58%) affected traits' means in ies reared at 25°C.On the other hand, 492 (86.5%) of candidate polymorphisms for CVEs showed temperature-speci c effects, although in this case the opposite pattern was found, since most polymorphisms (485, 98.5%) presented an effect on CVEs only at 17°C.Interestingly, no single SNP or InDel exhibited signi cant effects at both temperatures for the same trait, either for means, CVEs or a combination of both (means and CVEs).
The GWAS revealed 27 pleiotropic candidate polymorphisms (i.e., polymorphisms that were obtained simultaneously for different traits; matches between candidates for plasticity and the same trait at either temperature were not computed as pleiotropic).If 50bp was considered as the limit for LD (Mackay et al. 2012), the number of pleiotropic candidate polymorphisms was lowered to 19, all of them corresponding to LV and CVE of PH (Additional File 5).Moreover, according to pairwise LD heatmaps of those traits (Figure 3), all pleiotropic loci corresponded to 2 regions in LD, probably associated to polymorphic inversions present in the population (Huang et al. 2014).
For the candidate polymorphisms obtained in GWA studies, we found inverse relationships between the effect sizes and Minor Allele Frequency (MAF) for all traits, such that the less common variants had greater effects, either positive or negative (Additional Files 4 and 6).There were positive effect sizes for both Larval and Pupal Viability, wherein individuals homozygous for the major allele survived in a major proportion than individuals homozygous for the minor allele; at least for Larval Viability, these were comparable between temperatures.On the other hand, a negative effect size pattern was observed for variants affecting LDT and PH at 25°C.Minor frequency alleles of variants associated to variability for CVEs were associated to a decanalization of phenotypes for all traits, with candidate polymorphisms for DT traits showing the lower effects overall.436 candidate polymorphisms were also associated to regulatory regions, of which 406 (93.1%) were binding sites corresponding to 19 transcription factors and several hotspots (i.e.binding sites to multiple transcription factors) (Additional File 5).Moreover, 4 of these transcription factors, invected (inv), knot (kn), senseless (sens), and bric a brac 1 (bab1) were found as well as candidate genes.
A total of 333 genes were located within <1kb of candidate variants affecting early tness traits assessed in this study (Table 5, Additional File 4).The number of genes involved in natural phenotypic variation for each trait varied from 1 for PDT to 110 for the CVE of PV.No single pattern was detected for the number of genes affecting the combination of traits and temperatures.The frequencies of site classes for gene-related candidate polymorphisms did not differ from those corresponding to the whole polymorphic markers database (χ 2 4 = 7.68, p=0.1 (Mackay et al. 2012), with only 3.9% of the variants associated to polypeptide changes in gene products (Additional File 4).
Our data indicated that most of the genes were trait-speci c, considering that only 6% (20 candidate genes) were pleiotropic (Additional File 7).Of these 20 genes, only 2 (Semaphorin 1a [Sema1a] and luna) were pleiotropic affecting mean and CVE for the same trait (LDT).This number amounted to only 0.6% of the 333 candidate genes detected in this study, suggesting a decoupling between the genetic bases of CVEs and means for all traits.
Interestingly, PV means and CVEs, that presented signi cant correlations at both temperatures (Table 3), showed no pleiotropic genes or variants.The genes Ecdysone-inducible gene L1 (ImpL1) and Sema1a exhibited the highest pleiotropic effects, as each one was associated to natural phenotypic variation for 3 traits (LV/CVE LDT/CVE PH and LDT/CVE LDT/CVE PV, respectively) while 18 genes were pleiotropic for 2 traits.No genes presented natural variants simultaneously associated to variation in means for different traits.Moreover, surprisingly, none of the candidate genes was involved in phenotypic natural variation for the same trait at both temperatures.
According to enrichment analyses (Additional File 8), the complete list of candidate genes presented a larger proportion relative to the whole genome in terms related to development and function of the Central Nervous System (such as neurogenesis or locomotion).Moreover, genes related to ontogenetic transitions, growth, metabolism, feeding behavior, tissue repairing and immunity -generally, resource acquisition and allocationseemed to underlie variability for all assessed traits, their plasticities and microenvironmental robustness at both temperatures (Additional File 4).
Finally, to evaluate loci statistically associated to differences across life stages, we performed the complementary GWAS for DT and Viability at each temperature.Results indicate antagonistic pleiotropic effects for most of the signi cant loci for both traits (Figure 4, Additional File 9).Interestingly, almost all signi cant loci at 25º show the same effect signs (with minor frequency alleles associated to longer larval and shorter pupal DTs and with lower larval and higher pupal Viability).Signi cant loci at 17º showed no consistent effects, although all of them present antagonistic pleiotropy across stages.

Discussion
In recent years it has come to attention that studying the genetic architectures of developmental adaptive traits and their variations through ontogeny and phylogeny is essential to understand evolutionary dynamics (Mackay 2001, Hansen 2006).However, despite the recognition of the complexity in the relationship between genotype and phenotype for these traits, in which phenotypic plasticity, robustness and modularity at different levels hold important implications for the maintenance of genetic variability and for evolutionary dynamics, most genetic analyses devote little attention to these issues.Here we focused on Developmental Time, preadult Viability and Pupation Height in D. melanogaster, which contribute to preadult tness.We performed genetic analyses considering their phenotypic variability, plasticity and within-line canalization for two rearing temperatures in different stages along ontogeny, in order to understand the developmental processes behind the observed variation and the complex relationships between these tness components.We detected signi cant phenotypic variability for these traits and a strong genotype-environment interaction in this subset of lines of the DGRP population; a substantial amount of this variation could be explained by genetic factors.We also detected a reduction in phenotypic robustness at 17°C relative to 25°C, which was signi cant for DT traits and PV.It may be noted that although 17°C cannot be considered as a stressful temperature for the species, this nding was expected for a population adapted to a subtropical climate, which could have evolved epistatic interactions favoring a canalization of adaptive traits around phenotypic optima at higher temperatures.
In this sense, a change in the genetic networks at lower temperatures would imply the loss of this buffering.The apparent uncoupling between the genetic architectures of different traits, their plasticity and within-line canalization and also their temperature and stage speci city, which was previously observed in these and other By means of GWA Studies in DGRP lines, whose genomes are sequenced, we found 924 polymorphisms and 333 genes as candidates for underlying phenotypic variability for all traits' means and CVEs at both temperatures.It is worth noting that, even when several candidate polymorphisms were found for most traits, the number of assessed lines limited the detection of candidate variants to those with a MAF of 10% or higher; less common variants would have a greater phenotypic effect but were undetectable with the sample size used.Indeed, we found that alleles with a lower MAF were associated with greater effects, in this case corresponding to longer LDT, lower viabilities and higher CVEs for all traits.The latter was also found by Harbison et al. (2013) in a GWAS for sleep traits; the presence of decanalizing alleles in low frequencies can be explained because some individuals bearing them would present the same phenotype as those carrying the more frequent alleles, allowing them to escape the effects of stabilizing selection.However, compared with viability and PH, decanalizing alleles showed lower effect sizes for DT traits, as expected for characters under stabilizing selection (Gibson andWagner 2000, Gibson 2009).Regardless of this limitation, the focus of our work lies in characterizing how the genetic bases of pre-adult traits vary according to environmental conditions, ontogenetic stage and for different components of phenotype; therefore, we do not claim to have fully characterized the genetic architecture of the traits assessed.
Most of the candidate variants found were found in noncoding regions, a common result in GWAS studies in Drosophila and other organisms (Maurano et al. 2012, Harbison et al. 2013, Massey and Wittkopp 2016) which signals the importance of genetic variability affecting gene expression on phenotypic variation.Moreover, the abundance of variants located near or within the same gene on binding sites for only a few transcription factors suggests that individual variants did not affect the global patterns of expression, but entailed a temporally or spatially restricted effect.This modular characteristic of transcription factor action is already a recognized phenomenon (Wagner et al. 2007) which could also explain the independent effect of different loci within the same gene over different traits that we have found in this work.
The phenotypic analyses and GWA studies performed here point towards a complex and contingent genetic basis for developmental adaptive traits: as was suggested by the low cross-environment genetic correlations and the lack of conserved signi cant correlations between traits, most variants related to phenotypic variability were trait, temperature and developmental stage speci c, and/or were only related to a particular aspect of phenotype (means, CVEs, plasticities for means or CVEs).Moreover, instead of a general mechanism underlying robustness for all traits, we found that genetic bases for CVEs were trait-speci c, a result that is consistent with previous ndings by Morgante et al. (2015).
Even when we found pleiotropic genes, the speci c polymorphisms underlying variability for each trait were usually located in different regions.This phenomenon indicates a great degree of modularity at the genetic level, as distinct variants within the same genes may become selected independently.Also, their effect on phenotype may be probably affected by epistasis and the alteration of genetic networks by external or even internal conditions, which means that the same genetic variants are not necessarily expected to in uence phenotype in different genetic backgrounds, environments or along ontogeny (Mensch et  Interestingly, when analyzing variants that could account for the decoupling of larval and pupal traits, widespread antagonistic pleiotropic loci were found.While almost all of the candidate loci obtained at 25ºC were associated with longer larval and shorter pupal DTs and with lower larval and higher pupal viabilities, the effects of the candidate loci affecting both traits at 17ºC were not consistent in sign and more pronounced.
The study of the functional architecture of developmental traits could also explain how they affect each other.
Statistical correlations between the developmental adaptive traits assessed here have been found in several occasions (Casares and Carracedo 1987, Chippindale et al. 2004, Narasimha et al. 2015), and some were also observed here for the DGRP lines.However, it has come to attention that relationships between adaptive developmental traits are complex, statistical correlations do not hold universally (Chippindale et al. 2003), and according to our work, they may vary within the same population, depending on the environmental factors and ontogenetic stages.According to enrichment analyses, changes in the genetic networks affecting processes related to Central Nervous System development accounted, at least in part, for the phenotypic variation found in these lines.Several genes related to pathways (e.g.Insulin Receptor, Ecdysone Receptor, Octopamine Receptor or Short Peptide F Receptor pathways) controlling resource acquisition and/or allocation (i.e., metabolism, feeding, growth) were also obtained as candidates for all traits.Moreover, some of the candidate genes have been reported as pleiotropic in regard to the pathways they participate in, which raises the possibility that they mediate crosstalks between different cascades and coordinate biological functions.These complex genetic and environmental factors affecting the genetic networks behind developmental processes would explain phenotypic variation in all the traits considered but also the intricate relationships found between them.
Our results entail signi cant implications concerning evolvability for adaptive developmental traits.Modularity, either at the genetic, ontogenetic or anatomical level, limits the loss of genetic variability and favors adaptive radiation (Yang 2001, Hansen 2006, Mäkinen et al. 2018).Given the disparate behaviors and processes that characterize developmental stages and the varying challenges posed by different environments, the possibility of a somewhat independent evolution for traits corresponding to distinct periods of ontogeny in response to environmental heterogeneity would imply a greater phenotypic exibility in the face of changing conditions and limit the loss of genetic variability (Yang 2001, Del Pino et al. 2012, Petino Zappala et al. 2018).In this sense, even when loci presenting antagonistic pleiotropy across ontogenetic stages were found, they were temperature speci c; such trade-offs are expected to maintain natural phenotypic variability in the context of heterogeneous environmental conditions and epistatic effects (Huang et al, 2020).So does the decoupling between the genetic bases of means, plasticities and within-line canalization for adaptive traits when phenotypic optima vary across environments.If these aspects of phenotype are affected by different genetic variants which can evolve independently, different combinations of them may arise and become selected as ecological strategies as populations deal with environmental heterogeneity (Simons and Johnston 1997, DeWitt and Scheiner 2004, Simons 2014, Dewitt 2016).
All these factors favor the maintenance of genetic variability within natural populations and thus increase evolutionary potential (Hansen 2006, Schlichting, 2008).We therefore argue that, in order to understand the evolutionary dynamics underlying biological diversity, forthcoming studies should account for these phenomena emerging from the genetic architecture of developmental adaptive traits.
Our ndings suggest that the search for speci c variants or genes universally associated with complex traits or diseases may be futile.However, focusing on the genetic networks or processes in which candidate genes are involved would probably provide a more representative picture which could be extrapolated to different contexts.
Our results support this hypothesis, and future works on lines derived from other natural populations could con rm it if the same genetic networks and developmental processes also underlie phenotypic variability for these traits through the effects of different genetic variants.Means and standard deviations for within-genotype Coe cients of Environmental Variation of all traits at both temperatures.Signi cant differences for the same trait between temperatures are marked as * p<0.05, *** p<0.001.

Declarations
If plasticity is affected by genotype or genotype-environment interaction exists (Wijngaarden and Brake eld 2001, Callahan 2005), variation in phenotypic plasticity can determine different local strategies, avoiding intraspeci c competition (Del Pino et al. 2012, Forsman and Wennersten 2016).Plastic responses can also increase invasiveness (Bull 1987, Davidson et al. 2011) and limit the loss of genetic variability (see below).
Although previous work allowed to analyze phenotypic variability for these traits in natural populations or their responses to arti cial selection (Casares and Carracedo 1987, Welbergen and Sokolowski 1994, Chippindale et al. 2004, Yadav and Sharma 2014, Narasimha et al. 2015), and some studies have been carried out to study their genetic bases (Mensch et al. 2008, 2010; Riedl et al. 2007, Petino Zappala et al. 2018, 2019), these are, to our knowledge, the rst GWA analyses for these traits performed on lines derived from a natural population, providing a more representative sample of natural genetic variability (Rockman 2012, Gasch et al. 2016).Moreover, our analyses also assessed pre-adult tness traits in different environments and took into account within-line canalization patterns, while simultaneously decomposing preadult traits into their larval and pupal components.

Table 1 .
This work was supported by grants PICT-2012-0640 and PICT-2016-2256 from Agencia Nacional de Promoción Cientí ca y Tecnológica (FONCyT, PICT) and PIP 112-201101-00717 from Consejo Nacional de Ciencia y Técnica (CONICET).Con icts of interest: The authors declare no con icts of interests.Summary of the results of the correlation (Pearson / Spearman) analyses for all traits.

Table 2 .
Summary of the partial ANOVAs between temperatures for the CVEs of all traits.p-values are shown for each ANOVA for the effect of Temperature on CVEs of all traits.§ not signi cant after Bonferroni correction for multiple tests (p B =0.0042, see text for more details)

Table 3 .
Summary of the results of the genetic correlation analyses (Pearson / Spearman) for CVEs of all traits.Spearman correlation estimates are shown for each trait between temperatures (in bold on the diagonal), between each pair of traits within 17°C (above the diagonal) and between each pair of traits within 25°C (below the diagonal).Last two rows contain values for the correlation between mean and CVE for each trait.Signi cant values are indicated as *p<0.05,***p<0.001,§ not signi cant after Bonferroni correction for multiple tests (p B =0.0042, see text for more details).

Table 4 .
Summary of the ANOVAs for the effects of inversions and Wolbachia infection status.shown for the effect of Wolbachia infection and presence of inversions on all traits at both temperatures.Signi cant results are indicated as *p<0.05,**p<0.01,***p<0.001,§ not signi cant after Bonferroni correction for multiple tests (p B =0.0042).

Table 5 .
Summary of the GWAS results.