The evolution of larger size in high-altitude Drosophila melanogaster has a variable genetic architecture.

Important uncertainties persist regarding the genetic architecture of adaptive trait evolution in natural populations, including the number of genetic variants involved, whether they are drawn from standing genetic variation, and whether directional selection drives them to complete fixation. Here, we take advantage of a unique natural population of Drosophila melanogaster from the Ethiopian highlands, which has evolved larger body size than any other known population of this species. We apply a bulk segregant quantitative trait locus mapping approach to 4 unique crosses between highland Ethiopian and lowland Zambian populations for both thorax length and wing length. Results indicated a persistently variable genetic basis for these evolved traits (with largely distinct sets of quantitative trait loci for each cross), and at least a moderately polygenic architecture with relatively strong effects present. We complemented these mapping experiments with population genetic analyses of quantitative trait locus regions and gene ontology enrichment analysis, generating strong hypotheses for specific genes and functional processes that may have contributed to these adaptive trait changes. Finally, we find that the genetic architectures indicated by our quantitative trait locus mapping results for size traits mirror those from similar experiments on other recently evolved traits in this species. Collectively, these studies suggest a recurring pattern of polygenic adaptation in this species, in which causative variants do not approach fixation and moderately strong effect loci are present.

Introduction

Well into the genomic era, considerable debate persists over the types of genetic architectures that underlie adaptive evolution. For example, it is unclear how polygenic adaptive phenotypic changes tend to be—genes of major effect on adaptive traits are often reported (e.g. van’t Hof et al. 2011; Miller et al. 2014), and in the case of local adaptation, these are more likely to overcome the homogenizing force of migration (Yeaman and Whitlock 2011). However, it is possible that most adaptive events may instead involve large numbers of small-effect changes (Pritchard and Di Rienzo 2010; Rockman 2012). It is also unclear how often adaptive variants are selected as newly occurring mutations (e.g. Linnen et al. 2009), vs selection on standing genetic variation after an environmental change (e.g. Colosimo et al. 2005). In the latter case, the detection of “soft sweeps” is a distinct and more challenging exercise than for classic “hard sweeps” (Pennings and Hermisson 2006). It is also unclear how often adaptive variants actually reach fixation, vs remaining polymorphic due to factors such as traits reaching a new optimum or threshold value, changes in selective pressures, balanced equilibria such as heterozygote advantage, or ongoing migration (Stephan 2016; Höllinger et al. 2019; Thornton 2019; Barghi and Schlötterer 2020; Barghi et al. 2020; John and Stephan 2020).

Population genomic scans for natural selection provide some insight into the genetic basis of adaptive evolution, identifying large numbers of loci with signals of recent positive selection, and estimating the frequency at which different functional categories of sites are targeted. However, the biological basis of natural selection at these loci is usually not clear from genetic variation alone, and the properties of adaptive mutations may depend on the biological process (e.g. morphological vs physiological changes; Carroll 2008; Liao et al. 2010). Therefore, an essential complement to population genomic scans is detailed experimental case studies of the genetic basis of specific adaptive phenotypic changes, to gain a clearer and more nuanced understanding of how natural selection operates at the genetic level.

The molecular and evolutionary genetics model Drosophila melanogaster provides an efficient system for illuminating the genetic basis of evolutionary change, in part because of its ease of laboratory study, its well-developed molecular genetic toolkit, and its compact and well-annotated genome. D. melanogaster expanded from a warm ancestral range in southern-central Africa to occupy diverse worldwide environments (Sprengelmeyer et al. 2020). Latitude and especially altitude gradients allow the comparison of geographically proximate, closely related populations from contrasting environments. Phenotypic differences between genetically similar populations provide ideal raw material for studies of evolution at the genetic level, because the power of population genetic scans for local selection is maximized, and once the relevant genes are identified, the number of plausible causative mutations that differ between populations may be limited.

Size is a fundamental organismal quality. In D. melanogaster and other Drosophilids, larger body size is correlated with cooler latitudes (David et al. 1977; Gilchrist and Partridge 1999) and may provide a fitness advantage in cool environments (McCabe and Partridge 1997; Reeve et al. 2000). Instead of a direct effect of size on thermal tolerance (Drosophila are small enough to be virtually isothermic with their environment), higher larval density in the tropics may select for earlier pupation, leading to smaller adults, while in cooler regions viability selection may favor larger, more robust adults (Partridge and French 1996).

In other Drosophilid species, larger flies are also found at higher altitudes (Stalker and Carson 1948; Norry et al. 2001), but this phenomenon was not extensively studied in D. melanogaster until recently (Louis et al. 1982; Collinge et al. 2006). In the past decade, a unique highland Ethiopian population of D. melanogaster was found to be the largest known naturally occurring members of this species, with particularly enlarged wings (Pitchers et al. 2013; Klepsatel et al. 2013; Klepsatel et al. 2014; Fabian et al. 2015; Lack et al. 2016a, 2016b). The increase in wing size is associated with lower wing loading and could benefit flies in highland environments that are persistently cool (limiting the speed of wing movement) and feature thinner air (providing less resistance against fly wings). A plastic decrease in D. melanogaster wing loading that occurs at low developmental temperatures is associated with improved flight performance in cold (Frazier et al. 2008). Other studies have found that wing loading does not necessarily predict flight performance under different pressure or temperature conditions (Dillon and Frazier 2006; Hoffmann et al. 2007).

Comparing wing length between a highland Ethiopian population and a low-altitude ancestral range population from Zabmia, phenotypic differentiation (Q = 0.985) greatly exceeded genetic differentiation (genome-wide F = 0.151), implying that directional selection acted on wing length or a pleiotropically correlated trait (Lack et al. 2016a). The species is only estimated to have occupied the Ethiopian highlands about 2,700 years ago (Sprengelmeyer et al. 2020), or roughly 40,000 fly generations ago (based on 15 generations per year; Turelli and Hoffmann 1995; Pool 2015). In light of an effective population size on the order of 1 million for this lineage (Sprengelmeyer et al. 2020), the evolution of larger size has occurred on a recent population genetic time scale (∼0.01 autosomal coalescent units).

There has been some progress on understanding the tradeoffs and mechanisms involved in this population’s size evolution. Compared to a low-altitude Zambian population from the ancestral range, Ethiopian flies lay fewer but larger eggs, which develop into larger adults without prolonging the larval growth phase (Lack et al. 2016b). Ethiopian size changes were found to involve increases in cell size (likely a function of increased somatic ploidy; Smith and Orr-Weaver 1991; Edgar and Orr-Weaver 2001) as well as cell proliferation (Lack et al. 2016b). The evolution of larger wings in Ethiopian D. melanogaster was accompanied by a decanalization of wing development, implying that ancestral buffering mechanisms had been disrupted in the course of adaptive trait evolution (Lack et al. 2016a).

The genetic basis of Ethiopian size evolution has not been investigated. Outside Africa, initial progress has been made to identify genes underlying latitude-size clines in D. melanogaster outside Africa. In Australia, Dca and srp are potential contributors to wing and body size differences, respectively (Lee et al. 2011; Chen et al. 2012). Functional experiments (Carreira et al. 2009), association testing (Jumbo-Lucioni et al. 2010; Lafuente et al. 2018; Watanabe and Riddle 2021), and experimental evolution (Turner et al. 2011) have suggested that many genes could influence within-population body size variation. However, quantitative trait locus (QTL) mapping of size differences between high- and low-latitude populations has suggested a few major loci, with uneven chromosomal contributions not predicted by a highly polygenic model (Calboli et al. 2003). Hence, the polygenicity of body size variation may depend on whether diversity is examined within populations where stabilizing selection may predominate, or between populations where adaptive phenotypic evolution is suspected.

In this study, we aim to understand the genetic architecture of adaptive trait evolution, using the Ethiopian population’s thorax and wing size changes as model traits. Here, thorax length represents a proxy for overall body size, whereas wing length represents a trait that has particularly evolved in this population (Lack et al. 2016a). We focus on the polygenicity of trait evolution and genetic predictability within a population. We perform bulk segregant analysis to ascertain QTLs that are involved in thorax and wing size trait evolution. We also use population genetic statistics and gene ontology (GO) enrichment to find the evidence of local adaptation and to identify candidate genes for future functional investigation.

Materials and methods

Bulk segregant analysis

Here, we sought to map the larger thorax and wing size of a high-altitude Ethiopian population (from Fiche, Ethiopia, “EF,” 9.81°N, 38.63°E, alt. 3,070 m) compared to a low-altitude ancestral range population (Siavonga, Zambia, “ZI,” 16.54°S, 28.72°E, alt. 530 m). Four independent crosses were conducted, each involving a unique pair of EF and ZI strains: EF8N × ZI403N, EF15N × ZI366N, EF73N × ZI418N, and EF86N × ZI274N. These lines had been inbred for 8 generations from wild-caught isofemale lines, though inbred lines still contain significant levels of residual heterozygosity (Lack et al. 2016c). These lines had been found to be homozygous for common inversions—In(1)A, In(1)Be, In(2L)t, In(2R)NS, In(3L)OK, In(3L)P, In(3R)K, In(3R)Mo, and In(3R)P—by a combination of genomic data analysis (Lack et al. 2016c) and PCR (Lack et al. 2016a). All flies used were raised at 20°C on medium prepared in batches of 4.5 l water, 500 ml cornmeal, 500 ml molasses, 200 ml yeast, 54 g agar, 20 ml propionic acid, and 45 ml tegosept 10% (in 95% ethanol), at room temperature.

To determine which regions of the genome harbor the causative variants responsible for the evolution of larger thorax and wing size, bulk segregant analysis was performed to detect QTL. Four different population cages were started—1 for each of the Ethiopia-Zambia crosses mentioned. Each population cage is 28 cm × 14 cm × 15 cm and has 14 vials containing the above medium. In each population cage, reciprocal crosses were established between 8 inbred parental individuals of each strain (Zambia and Ethiopia). From each reciprocal cross, 125 F1 offspring of each sex were used to establish the second generation. For the duration of the experiment, nonoverlapping generations were maintained at ∼1,200 individuals (Fig. 1). Adult flies were allowed to lay eggs on the food for 1 week before being removed. The food vials were replaced when adult flies in the cage were 7–10 days old. At the 16th generation, 600 3–5-day-old female flies from each population cage were measured as described below. For each trait, thorax size and wing size, the flies were placed into pools constituting the 10% smallest (N = 60) and 10% largest (N = 60) individuals, with the remaining individuals discarded. The QTL detection power of this experimental design is expected to be higher for QTLs explaining more than 15% of the parental strain trait difference (Pool 2016—Fig. 7A).

Fig. 1.

The bulk QTL mapping experimental design is illustrated. As further described in the Materials and Methods, F1 offspring of reciprocal crosses were allowed to interbreed in a relatively large population without selection until the F15 generation, at which point 600 females were sorted to obtain the top and bottom 10% for a size trait for sequencing. This design allows a large number of unique recombination events to take place, which should improve mapping performance.

Body size

To measure thorax and wing size, we followed the protocol described in Lack et al. (2016a). Thorax size measurements were in 3–5-day-old adult females. From each mapping cross, females were photographed with a digital camera attached to a stereo dissecting microscope (AmScope SM-4BX), and thorax length was measured from the base of the anterior humeral bristle to the posterior tip of the scutellum (see Lack et al. 2016b, Fig. A2B). For wing size, we also examined 3–5-day-old adult females from each of the mapping crosses. For 5 females per cross, a wing was removed and photographed at 509 magnification using a digital camera attached to a compound microscope (Olympus BH-2). The length and depth of each wing were then measured using ImageJ version 1.48 (http://imagej.nih.gov/ij/; last accessed January 5, 2022), and we measured a straight line drawn from the intersection of the anterior crossvein and L4 longitudinal vein, to where the L3 longitudinal vein intersects the wing margin (see Lack et al. 2016b, Fig. A2A). For depth, we measured a straight line from the intersection of the L5 longitudinal vein and the posterior wing margin, passing through the intersection of the posterior crossvein and L4, and terminating at the anterior wing margin. For wing area, we imaged individual wings using the “wing grabber” apparatus described by Houle et al. (2003), wing area was determined by outlining each wing using ImageJ version 1.48 (http://imagej.nih.gov/ij/; last accessed January 5, 2022), and the reported area for each cross is the mean of the 5 wings.

Genome preparation

We sequenced the genomes of pooled samples (N = 30 individuals) for the parental lines and 2 such pools for each of the large- and small-size groups (0–5% and 5–10% extremes for each direction, summing to N = 60 total for each extreme). Genomic DNA was obtained using a chloroform extraction and ethanol precipitation protocol. The DNA was fragmented with a Bioruptor sonicator (Diagenode), and paired-end libraries with ∼300-bp inserts prepared using NEBNext DNA Library Prep Reagent Set for Illumina (New England Biolabs no. E6000L). Each library’s concentration and quality was analyzed with an Agilent 2100 Bioanalyzer (Agilent Technologies, Inc.). The prepared libraries were sequenced at UW-Madison Biotechnology Center on the Illumina HiSeq 2000 platform. Having concluded that the full 10% extremes would best be analyzed together (Pool 2016), we merged reads from the 0–5% and 5–10% pools (similar numbers of reads were obtained from these pools in each case) before proceeding with the analysis.

Genome alignment

All the raw data that passed the Illumina filters were processed using a Perl-scripted pipeline. Reads from each sequenced genome were mapped to the D. melanogaster reference genome (release 5.57) obtained from Flybase (www.flybase.org), with the default parameters in BWA ver. 0.6.2-r126 (Li and Durbin 2009). Using Stampy ver. 1.0.21(Lunter and Goodson 2011), the BAM files were then remapped. With samtools ver. 0.1.18 (Li et al. 2009), reads were filtered for a mapping quality of 20 and for proper pairs. The BAM files were further processed by removing unmapped reads and sorted by coordinate, and PCR duplicates were marked using Picard ver. 1.109 (http://picard.sourceforge.net; last accessed January 5, 2022). To improve the alignment around indels, we used GATK ver. 3.2 (McKenna et al. 2010). The average depth of coverage per genome was calculated for the parental lines and the low- and high-tolerant lines (Supplementary Table 1).

QTL mapping

Synchronized mpileup files for the aligned genomes were created with the PoPoolation2 ver. 1.201 software package (Kofler et al. 2011). The 2 large (and 2 small) pools from a given cross were then combined with a custom perl script. Ancestry difference (ad) was then calculated with each biallelic SNP (Bastide et al. 2016). Ancestry difference estimates the difference between the proportion of the large-fly pool’s sequencing reads carrying an allele from the large (Ethiopia) parental line and that same proportion from the small-fly pool. It was estimated as: where p is the frequency of the major allele in the large parent, p is the small parental allele, f is the frequency of the large parent allele in the large pool of F16 offspring, and f is that same allele’s frequency in small F16 offspring. The 5 chromosomal arms (X, 2L, 2R, 3L, and 3R) were divided into windows based on SNP density (Lack et al. 2015) which created 2,728, 3,131, 2,357, 2,956, and 2,935 windows, respectively, each roughly 8.4 kb in size on average. Only sites that had a parental strain frequency difference of ≥0.25 were used in the analysis, to avoid noisy ancestry proportion estimates from SNPs with modest frequency differences between parental strains. A simulation-based inference for BSA mapping (SIBSAM) was performed (Pool 2016) to identify significant QTLs and calculate their confidence intervals and effect sizes. SIBSAM is able to evaluate both primary QTL peaks and flanking secondary QTL peaks, evaluating whether ragged peaks contain significant evidence for more than 1 QTL. Forward simulations incorporate recombination in multiple individuals for multiple generations, selection on phenotype in the final generation with additive gene effects, plus environmental variance, and then the sampling of sequence reads to obtain ad.

Genetic differentiation and GO enrichment analysis

QTLs identified in the previous step will contain many genes that may or may not be involved in the evolution of these traits. To help identify the causative genes within the significant QTLs for thorax and wing size, window F and maximum SNP F per window (“SNP F”), and the haplotype statistic χ (Lange and Pool 2016) were analyzed. Genomes from Zambia (n = 197) and Ethiopia (n = 68) were used from the Drosophila Genome Nexus (Lack et al. 2015). The 2 F statistics, quantifying allele frequency differences between populations, were used because window F should be sensitive to locally adaptive selection sweeps with larger linkage effects such as hard sweeps, while SNP F should be sensitive to narrower sweep signals such as those associated with selection on standing variation that result in soft sweeps (Pennings and Hermisson 2006). The χ compares length of identical haplotype blocks among individuals in one population vs another.

For our local adaptation genome-wide scans, each of the 5 chromosomal arms (X, 2L, 2R, 3L, and 3R) were divided into windows based on SNP density (Lack et al. 2015), which created 2,728, 3,131, 2,357, 2,956, and 2,935 windows, respectively, each roughly 8.4 kb in size on average. To narrow down potential candidate genes, a chromosomal arm quantile outlier approach was used to identify genes with an extreme population genetic signal. For each statistic and each chromosome arm separately, we defined a window’s quantile as the proportion of all windows on that chromosome arm with an equal or greater value for that statistic. We classified windows that were in the top 2.5% quantile for any of the 3 statistics as outlier windows. We then grouped neighboring outlier windows together into outlier regions, since they may reflect the same instance of local adaptation. To form an outlier region, a maximum of 2 nonoutlier windows were allowed between 2 outlier windows. Genes associated with outlier windows (overlapping them or the nearest gene in either direction) were retained for subsequent analysis.

We performed a GO enrichment analysis to identify potential functional categories that may contribute to the contrasting phenotypes found between the Zambia and Ethiopia populations. The outlier genes that were identified in the significant QTL regions were used for window-based GO enrichment analysis (Pool et al. 2012). A GO enrichment analysis was conducted for both thorax and wing size. A P-value was calculated based on the probability of observing a given number of outlier genes from a GO category. P-values were obtained from permutation in which outlier regions were randomly reassigned 10,000 times.

Results

We used bulk segregant analysis to perform QTL mapping for both thorax and wing length using 4 different unique between-population crosses. Each mapping population used individual inbred strains from an ancestral range Zambia population with smaller thorax and wing length, and from the high-altitude Ethiopia population that has evolved larger thorax and wing length. In our bulk segregant analysis, offspring of reciprocal crosses were allowed to interbreed for 16 nonoverlapping generations without selection at a large population size (N ≈ 1,200). After the 16th generation, 600 adult females were measured for both thorax and wing length and the top and bottom 10% of individuals were grouped for pooled genomic sequencing (Fig. 1; Materials and Methods). SIBSAM (Pool 2016) was then used to identify primary and secondary QTL peaks, along with their estimated effect sizes and genomic confidence intervals.

For thorax length, 4 Ethiopia × Zambia mapping crosses revealed a total of 12 significant peaks (Fig. 2 and Supplementary Table 2). The EF8N cross had 1 significant peak with an estimated effect size of ∼17%. EF15N had 2 significant peaks, each having an estimated effect size of ∼15%. EF73N had the most significant peaks with a total of 5, and these had estimated effect sizes that ranged between 12% and 20%. EF86N had 4 significant peaks with estimated effect sizes between 13% and 16%.

Fig. 2.

Significant QTL peaks for 4 Ethiopia/Zambia thorax length crosses. A point for each ∼8-kb window corresponds to the average difference in ancestry from the larger parental strain between the large and small F16 pools (y-axis). Significant primary or secondary QTL peaks are denoted with an arrow. The significance threshold for primary peaks is approximately 0.17.

For wing length, these same 4 crosses revealed a total of 33 significant peaks (Fig. 3 and Supplementary Table 3). EF8N had a total of 12 significant peaks, with estimated effect sizes that ranged between 7% and 24%. EF15N had 3 significant peaks, with estimated effect sizes that range between 16% and 24%. EF73N had a total of 10 significant peaks, with estimated effect sizes that ranged between 6% and 27%. EF86N had 8 significant peaks, with estimated effect sizes that ranged between 11% and 25%.

Fig. 3.

Significant QTL peaks for 4 Ethiopia/Zambia wing length crosses. A point for each ∼8-kb window corresponds to the average difference in ancestry from the larger parental strain between the large and small F16 pools (y-axis). Significant primary or secondary QTL peaks are denoted with an arrow. The significance threshold for primary peaks is approximately 0.17.

In general, very different QTL landscapes were observed between independent Ethiopia/Zambia crosses (Fig. 4). In some cases, QTLs do overlap between crosses, which may reflect either chance (different QTLs located close together) or else genuine sharing of causative variants underlying thorax and/or wing size. We identified QTL overlap when the QTL peak of 1 cross overlaps with the genomic confidence interval of another cross. For thorax length there are no regions between the 4 Ethiopia crosses where a QTL peak overlapped with another peak’s genomic confidence interval (Fig. 4). However, for wing length between the 4 crosses, there were 12 regions where QTL peaks overlapped with genomic confidence intervals involving 12 of the 33 QTLs (Fig. 4). Within these overlapping peaks, there are no overlap between all 4 crosses.

Fig. 4.

The locations of significant QTLs on the 5 euchromatic chromosome arms of D. melanogaster. The colors indicate for 4 Ethiopia strains used in mapping crosses for a) thorax length and b) wing length. The width of each box indicates the 90% CI of each QTL. Intervals that are less than 10 kb in width are marked with triangles. Dotted gray lines indicate Mb increments.

Some differences in the significant QTLs between crosses could represent chance detection of a shared QTL in some crosses but not others. However, with this experimental design, we expect to have >90% power to detect a QTL with 20% effect size (Pool 2016). Hence, at least for several of the strongest of the QTLs detected here, their absence in other crosses is likely to reflect real differences in genetic architecture.

Potential targets of local adaptation within QTL regions

Regions of the genome where the Zambia and Ethiopia populations greatly differ in their genetic variation may harbor genes involved in these adaptive traits. We used 3 population genetic statistics, window F, maximum SNP F within a window, and the haplotype statistic χ to identify possible candidate genes for body and wing size evolution within the significant QTLs. Using 3 different statistics is advantageous due to the differing power each statistic has in detecting local adaptation, depending on whether selective sweeps are complete or incomplete, or hard vs soft (Lange and Pool 2016). A quantile approach was used to identify population genetic outlier regions within QTLs that had one of the 3 statistics with a quantile of below 0.025. These local adaptation candidate regions are typically much narrower than QTLs, and hence, there can be multiple outlier regions per QTL, but most outlier regions are associated with just a few genes (Supplementary Tables 4 and 5). There are many genes within these outlier regions with no known role in either thorax or wing size. However, there are also genes known to be involved in size regulation.

For thorax length, genes corresponding to QTLs and population genetic outliers that are known to be involved in growth included bbc (Liu et al. 2014), ct (Thumm and Kadowaki 2001), msn (Kadrmas et al. 2004; Carreira et al. 2009), RasGAP1 (Dworkin and Gibson 2006), scyl (Reiling and Hafen 2004), spi (Nagaraj et al. 1999), srp (Bánréti et al. 2012; Chen et al. 2012), and tara (Bejarano et al. 2008). Of these, bbc and RasGAP1 provide examples of loci with promisingly narrow F peaks at the SNP level (Fig. 5), which may merit targeted investigation by future studies. msn and RasGAP1 are also within wing QTLs and are therefore relevant to the analysis described below as well.

Fig. 5.

At 2 candidate genes identified for thorax length evolution (bbc and RasGAP1), small numbers of SNPs show the highest F values between Ethiopia and Zambia. Depicted above is the gene transcript, while the x-axis indicates kb position along the relevant chromosome arm (release 5).

Within the outlier regions for wing length, these genes included caps (Milán et al. 2001; Carreira et al. 2009; Jumbo-Lucioni et al. 2010; Watanabe and Riddle 2021), dally (Takeo et al. 2005; Carreira et al. 2009), Dlish (Wang et al. 2019), Dronc (Verghese et al. 2012), fj (Villano and Katz 1995), Gbp1 (Koyama and Mirth 2016), Pka-C3 (Dworkin and Gibson 2006; Jumbo-Lucioni et al. 2010), and salr (Wang et al. 2017). Dlish and Pka-C3 are examples of genes with individual SNPs having high F values (Fig. 6). Functional testing will be needed to establish if genetic variants found within these genes are indeed responsible for the associated phenotypes.

Fig. 6.

Peaks of SNP F center on 2 candidate genes identified for wing size evolution (dlish and Pka-C3), showing elevated genetic differentiation between Ethiopia and Zambia at these genes. Depicted above is the gene transcript, while the x-axis indicates kb position along the relevant chromosome arm (release 5).

GO enrichment

We conducted individual GO enrichment analysis for thorax and wing length. We used only the genes found in the outlier windows located within significant QTL regions from the 4 crosses. Full results are presented in Supplementary Tables 6 and 7. Categories with raw P-values below 0.001 included some that are either known or potentially involved in body and wing size. For thorax size, the top categories included: negative regulation of Ras protein signal transduction (Prober and Edgar 2002), regulation of protein polymerization (Fernández et al. 2011), and brahma complex (Krupp et al. 2005). For wing size, the top categories included: neurogenesis (Rutledge et al. 1992), ubiquitin protein ligase binding (Cornell et al. 1999), cellular amino acid catabolic process (Zinke et al. 1999), cellular response to anoxia (Heinrich et al. 2011), transmembrane transport (Bartscherer et al. 2006), and negative regulation of proteolysis (Lee et al. 2001). Some of these functional processes might underlie Ethiopia size adaptation, while others may be driven by unrelated trait evolution in this high-altitude population.

Discussion

We employed quantitative and population genetic strategies to investigate the genetic architecture of adaptive size evolution in our highland Ethiopia population. Our bulk segregant analysis revealed that between the 4 crosses, thorax size has 12 associated QTLs with moderate-to-large effect (∼13–20%). However, between the 4 crosses wing size has 33 QTLs with small-to-large effects QTLs (∼6–27%). A greater ability to detect wing length QTLs than thorax length QTLs may reflect the greater magnitude of the population difference in this trait (Lack et al. 2016a).

One striking result was the lack of QTL overlap between crosses for either thorax or wing size. Between the 4 thorax crosses, there is no overlap. This is especially notable given that we have almost have very high power to detect QTLs with effect size of 20% (Pool 2016) and yet the QTL on chromosome arm 2L with ∼20% effect size is not present in any other cross. For wing size, there was overlap in only 12 of the 33 QTL regions and no overlap between all 4 crosses. The QTLs with the 3 largest effect sizes of over 25% are present in only one cross. The low QTL overlap between crosses could reflect persistent genetic variation at causative loci in the Ethiopian and/or Zambian populations. Given that the Ethiopian population appears to have experienced directional selection for larger size and still maintains similar genetic variance for size traits as Zambia (Lack et al. 2016b), we suggest that some favored size variants have not reached fixation in the Ethiopian population. There are multiple reasons why favored alleles might not fix, including the Ethiopian population reaching its new optimum or threshold trait value (especially if ample standing variation means that not all large alleles needed to fix), heterozygote advantage, or ongoing adaptation. Indeed, simulation and theory have shown that depending on the genetic architecture of an adaptive trait, nonfixed causative variants may be the norm (Stephan 2016; Höllinger et al. 2019; Thornton 2019; Barghi and Schlötterer 2020; Barghi et al. 2020; John and Stephan 2020).

Our conclusions of persistent variability underlying an evolved trait mirror similar results for pigmentation (Bastide et al. 2016) and for ethanol resistance (Sprengelmeyer and Pool 2021) in this same population and others, all from mapping experiments with similar design and scale. With 5 traits now examined (ethanol resistance, abdominal background color, abdominal stripe width, thorax length, and wing length), consistent patterns are starting to emerge. First, these traits each average at least a few detectable QTLs per cross, with means ranging from 3 to 8.25 (Table 1). Second, there is notably little QTL peak overlap between parallel mapping crosses involving different strains from the same populations, with QTL overlap proportions ranging from 0% to 35% (Table 1). For each of these traits, there are moderately large effect QTLs, associated with very high detection power (Pool 2016), which are not present in other crosses. Hence, at the population level, it is fair to say that each of these traits is at least moderately polygenic and involves nonfixed differences between populations. Third, moderately strong QTLs are consistently present in any given cross, with the average QTL effect size ranging from 13% to 19% (Table 1), although undetectable smaller effects may be present as well.

Table 1.

The results of bulk QTL mapping experiments for 5 different traits.

Trait	Avg. QTLs per cross	Pairwise QTL overlap (%)	Avg. effect size (%)
Ethanol resistance	8	35	14
Thorax length	3	0	15
Wing length	8.25	36	17
Stripe width	3	8	19
Background color	5.67	18	19

All mapping used the same experimental design described in the Materials and Methods, aside from minor variation in the number of generations of interbreeding (15–20). Both the pigmentation stripe and pigmentation background data are from Bastide et al. (2016), while the ethanol results are from Sprengelmeyer and Pool (2021). Listed are the number of significant QTLs for each mapping population, the proportion of QTLs that overlap between parallel crosses from the same 2 populations, and the average QTL effect size across all mapping crosses.

Results compatible with the above findings have also been obtained from other experimental scenarios and from other species. For example, Barghi et al. (2019) studied the adaptation of replicate D. simulans populations to high temperature and found that the resulting heterogeneity in genomic responses among replicates was well-modeled by a scenario of genetic redundancy in which not all favored variants are needed to achieve an adaptive change. An experimental selection study involving desiccation resistance in D. melanogaster also observed considered heterogeneity among replicates (Griffin et al. 2017). Among human populations, the strongest known examples of local adaptation involve nonfixed genetic differences, including lactase persistence in northern Europeans and some African populations (e.g. Tishkoff et al. 2007; Itan et al. 2009), high-altitude adaptation in Andeans/Ethiopians/Tibetans (e.g. Bigham et al. 2009; Yi et al. 2010; Huerta-Sánchez et al. 2013), and enhanced diving ability in the Bajau (Ilardo et al. 2018). Likewise, shifts in natural stickleback populations have also involved nonfixed genetic changes, with pleiotropy suggested to play a role (Rogers et al. 2012). Similarly, pleiotropic effects of insecticide resistance variants may help explain why they are often in natural populations of D. melanogaster and other species (e.g. Catania et al. 2004; Clarkson et al. 2021). We also note that a population genetic analysis of African populations of D. melanogaster found ample evidence for incomplete sweeps (Vy et al. 2017).

Polygenic adaptation may have diverse outcomes, depending in part on the number of segregating variants at the onset of selection that affect a trait, as well as the magnitudes of their effect on the trait relative to the shift in trait optimum. While each of the traits summarized above might be described as “polygenic,” it is worth considering the type of polygenic adaptation that these mapping studies imply. The persistently variable genetic basis of these evolved traits may suggest a scenario of abundant standing genetic variation prior to selection for each of these traits. In light of the consistent presence of moderately strong QTLs for these traits, such standing variation may have included relatively large effect loci, which would experience relatively stronger directional selection during the trait’s evolution. An abundance of standing variation is consistent with the large population size and high genetic diversity of this species (e.g. Sprengelmeyer et al. 2020). Further studies will be needed to quantify the models of polygenic adaptation that experiments such as ours indicate, and to assess whether such persistent variability is a widespread outcome of trait evolution not only in this species but also across the tree of life.

Data availability

All raw sequence data have been deposited in the NIH Short Read Archive, with accession numbers given in Supplementary Table 1. The scripts used for SIBSAM can be found at: http://github.com/JohnEPool/SIBSAM1; last accessed January 5, 2022.

Supplemental material is available at G3 online.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.