Introduction

Body size clines conforming to Bergmann's rule are widespread, with larger individuals found in higher latitudes. Insect examples include the honeybee (Alpatov 1929), ants (Cushman et al., 1993), ant lions (Arnett and Gotelli 1999), the yellow dung fly (Blanckenhorn and Demont 2004), the housefly (Bryant 1977) and various drosophilids (Coyne and Beecham 1987; Imasheva et al., 1994; James et al., 1995; Karan et al., 1998; Huey et al., 2000; Loeschcke et al., 2000). In Drosophila melanogaster, the occurrence of Bergmann's cline in multiple continents suggests that strong natural selection is acting on body size and/or other correlated traits (Hoffmann and Weeks 2007). Like many other ecological traits, empirical evidence is needed to compare the genes underlying similar adaptive shifts across different species (Hoffmann and Willi 2008).

Although body size in Drosophila is sensitive to environmental factors, the clinal variation in size has a strong genetic basis because patterns persist even when flies are reared in controlled conditions. Linkage analyses by Gockel et al. (2002) and Calboli et al. (2003) on the Australian and the South American D. melanogaster populations have independently shown that wing area, often taken as a proxy for body size, is controlled by quantitative trait loci (QTL) on chromosome arms 3R and 2R. The congruency in QTL positions between the two studies points to a similar genetic mechanism in different continental populations.

However, the QTL on 3R also coincides with the cosmopolitan chromosomal inversion 3RP (In(3R)Payne), which shows a strong clinal pattern and accounts for 30–60% of the size variance (Weeks et al., 2002; Rako et al., 2006). The presence of this inversion polymorphism in mapping populations impedes fine-scale mapping. Linkage disequilibrium is high in the In(3R)Payne region, but meiotic recombination is not entirely prohibitive to mapping. Kennington and co-workers undertook an association study using random individuals from a mid cline population and identified three smaller genomic areas within In(3R)Payne that showed significant association with wing size (Kennington et al., 2006; Kennington et al., 2007). More recently, McKechnie et al. (2010) reported—based on transgenic overexpression, genotype–phenotype association analysis and strong clinal pattern of allelic frequencies—that an insertion/deletion polymorphism in the promoter of Dca (Drosophila cold acclimation) is associated with wing size variation. Hence, Dca represents a gene with relatively large effect on the wing size variation in nature, alongside the contribution from loci within In(3R)Payne (McKechnie et al., 2010).

QTL mapping of wing size in D. melanogaster has suggested overlapping genetic mechanisms in Australia and South America (Gockel et al., 2002; Calboli et al., 2003). However, whether the same genetic basis is also conserved in another species is yet to be established. Similar to D. melanogaster, D. simulans shows a strong latitudinal wing size cline in Australia (Arthur et al., 2008). But unlike D. melanogaster, the absence of major cosmopolitan inversions in D. simulans populations means that genetic mapping at a finer scale can be carried out more easily.

The current study follows an approach similar to Gockel et al. (2002) to investigate the genetic basis of the adaptive wing size variation in D. simulans using inbred lines originating from cline-end populations (Queensland and Tasmania). Mendelian genetic analysis indicated the basic genetic architecture—that large wing size is an autosomal dominant trait, with minor contribution from the sex chromosome, and unequal additive effect from chromosomes 2 and 3. QTL mapping and comparative analysis indicated a genetic basis that is distinct from that reported in D. melanogaster.

Materials and methods

Construction of inbred lines for QTL analysis

Mapping lines used in this study were derived from isofemale lines established from Sorell, Tasmania (latitude 42° 46′ 11′′ S, longitude 147° 34′ 35′′ E) and Maryborough, Queensland (latitude 25° 32′ 44′′ S, longitude 152° 41′ 05′′ E), collected in April 2005. All lines were reared at 19 °C with 70% relative humidity under continuous light in 42 ml vials on 10 ml of Bloomington Drosophila medium adjusted by doubling the quantity of soy flour and replacing light corn syrup with dextrose (http://flystocks.bio.indiana.edu/Fly_Work/media-recipes/bloomfood.htm). We treated the fly media with an antifungal agent (0.14% w/v Nipagin). The isofemale lines underwent 14 generations of inbreeding. In each generation, multiple brother–sister mating cages were set up; once eggs were collected, wings of the parents were measured and only cages with the desired parental phenotypes (smallest 10% for Maryborough and largest 10% for Sorell) were retained to initiate the following generation. Lines L20 and L71 were chosen for QTL analysis for the following reasons: they showed the most extensive between-line size divergence (∼4.4 s.d. values apart), and had minimum within-line wing size variation compared with other line combinations.

Mapping crosses

Owing to genetic dominance (Figure 1), backcrosses to L20 were used for linkage mapping. To establish a male informative cross, an F1 male from the L20 × L71 cross was paired with an L20 virgin to produce the backcross generation. Four such single-paired crosses were used for wing measurement and genotyping, and results were combined to estimate chromosomal contributions. To construct a female informative cross for mapping within chromosomes, an F1 female from L20 × L71 was allowed to mate with an L20 male to produce the backcross progeny. We measured 17 such single-paired female informative pedigrees and selected two, F243 and F257, for initial QTL analysis because they had similar average size (F243=1.52 mm; F257=1.52 mm), variance (F243=0.69 mm; F257=0.68 mm) and brood size (F243=27 males; F257=25 males). Seven additional female informative families (155 male progeny) were subsequently used to verify the major QTL. These additional families were derived from the same isofemale lines that produced F243 and F257.

Figure 1
figure 1

Genetic architecture of wing size variation in D. simulans males. Chart shows the average wing centroid size of the parental lines, F1's and the backcross progeny. The symbol ‘S’ on the horizontal axis represents a chromosome derived from L20, or the small line; ‘B’ represents a chromosome derived from L71, or the big line; ‘Y’ represents the Y chromosome. The order of the letters represents chromosome number. For example, SY:SS:SB=males carrying an X chromosome from the small line, two copies of chromosome 2 from the small line, one copy of chromosome 3 from the small line and the other copy from the big line. Error bars represent s.e. values of the mean wing centroid size based on 15–45 individuals.

Phenotype measurements

One wing from each individual was dissected from the thorax and placed flat on a transparent double-sided sticky tape, between a cover slip and a microscope slide. Wing size was measured by estimating the centroid size using the landmarks reported in Rako et al. (2006). Wing images were captured using a digital camera (PixeLink, Vitana Corp., Ottawa, Canada) attached to a compound microscope (WILD M3B) at × 40 magnification. The images were processed using TPS-Dig version 1.2 software (developed by F. James Rohlf, http://life.bio.sunysb.edu/morph/index.html) on which landmarks 1–8 were placed on specific vein positions. The centroid size of each wing, defined as the square root of the sum of squared distances of the 8 landmarks in the same order, was computed using the CoordGen6d programme (part of the IMP Suite developed by H. David Sheet, http://www3.canisius.edu/~sheets/morphsoft.html) based on the recorded landmark information from the TPS-Dig images. For single-marker regression analysis, because of the difference in phenotypic means and variances among the nine mapping families, centroid sizes were standardised within each family. The resulting z-scores were used for regression analysis.

DNA extraction

Genomic DNA was extracted from adult whole bodies using the chelex method. Single flies were homogenised using the Mixer Mill for 4 min at 25 hz in 250 μl of chelex solution, which was made up of 5% (w/v in distiled water) chelex 100 resin (Bio-Rad; Cat. No. 142–1253, Gladesville, NSW, Australia), 40 μg of proteinase K (Roche; Cat. No. 03115828001, Castle Hill, NSW, Australia) and two 3-mm glass beads (Ajax FineChem; Cat. No. 1700-500G, Taren Point, NSW, Australia). Homogenates were incubated at 65 °C for 30 min followed by a second incubation at 95 °C for 10 min. Samples were centrifuged for 4 min at 20800 r.c.f. and ∼150 μl of the supernatant was transferred to a fresh tube containing 20 μl of 0.1 × Tris EDTA (pH 8). The stock DNA was then diluted 10-fold and 1 μl of the diluted DNA was used in PCR.

Marker design

To develop markers for mapping, we made extensive use of the whole-genome sequences of D. simulans and D. melanogaster. The choice of markers was based primarily on their location in the genome, and to a lesser extent their gene structure. To ensure a good coverage, we opted for evenly spaced markers along chromosomes 2 and 3. Annotated or predicted features such as coding sequences, exons, introns and untranslated regions were obtained from Flybase. D. simulans genomic trace files (in SFC format) were retrieved from National Center for Biotechnology Information (NCBI) trace archives. Sequences were aligned in Sequencher 4.7 (Gene Codes), and primers were designed from highly conserved sequences to flank a small variable region (<150 bases). We imposed such a stringent size limitation to ensure high sensitivity of subsequent high-resolution melt (HRM) analysis. All primers were designed either by hand or using the Primer3 program (http://frodo.wi.mit.edu/primer3/) with optimal melting temperature set at 65 °C. Marker and primer sequence information are summarised in Table 1.

Table 1 Gene markers developed in this study

Genotyping by HRM analysis

All PCR amplifications and HRM analyses were carried out using the Roche LightCycler 480 system (384-well format). The 10-μl PCR reaction contained 1 μl of diluted DNA (see DNA extraction), 1 μl of primer mix (4 μM each), 1 μl of the 10 × reaction buffer, 0.8 μl of dNTP mix (2 mM), 0.4 μl of MgCl2 (50 mM), 0.25 μl of the LightCycler 480 High Resolution Melting Master (Roche; Cat. No. 04909631001), 0.01 μl of IMMOLASE DNA polymerase (Bioline; Cat. No. BIO-21047, Alexandria, NSW, Australia) and UltraPure DEPC-treated water (Invitrogen; Cat. No. 750023, Mulgrave, Victoria, Australia) to make up the remaining volume. Thermocycling conditions were as follows: 95 °C for 10 min, 50 cycles of 95 °C for 5 s, 60 °C for 10 s and 72 °C for 15 s. One fluorescence acquisition was obtained after each 72 °C step. Products were heated to 95 °C for 1 min, cooled to 40 °C for 20 s and raised to 65 °C. As temperature increases gradually from 65 to 95 °C, fluorescence data were acquired continuously. These fluorescence records were used in the HRM analysis using the Genescan module in the LightCycler 480 software package. The fluorescence signals for a given gene were first normalised by defining a 1 °C range before and after the actual melting temperature, or the temperature at which 50% of the amplicons remain double stranded. The melt curves were further normalised to account for slight variation in temperature control across wells. This was done by raising the horizontal threshold level to 5%. Melt curves were then grouped (genotyped) according to similarity in the melt properties.

Linkage analyses

Markers were first tested on the four parents of the mapping crosses (F243 and F257); informative markers were then used to screen the 52 backcross progeny. On the basis of the crossing scheme, two types of melt curves (or genotypes) were expected among these progeny: one corresponds to homozygous L20 and the other to heterozygous L20-L71. A genetic map was constructed de novo based on these 47 markers using MapMaker (Lander et al., 1987). Genotype matrix of all markers and their genetic distances (Mapmaker output) were analysed in Windows QTL Cartographer version 2.5 (Wang et al., 2010), using the composite interval mapping module. Experimental-wise significant threshold level (logarithm of odds=2.5) was determined by 10 000 permutations. Window size for the composite interval mapping analysis was 10 cM; the default composite interval mapping model (model 6) and the forward regression method were applied. A subset of these markers (L208, L209, L218, L219, L222, L311, L318, L319, L325, L328 and L333) were genotyped on 155 additional male progeny from seven other mapping families for regression analysis, that is, regression of the standardised phenotypic values on progeny genotypes (0=homozygous L20; 1=heterozygous L20/L71).

Results

Genetic architecture of wing size variation in D. simulans lines L71 and L20

Comparison of chromosomal configurations and the phenotype in parental lines, F1's and the backcross progeny revealed the basic genetic architecture of wing centroid size (Figure 1). A substantial size difference was observed between the parental lines (t=20.2, d.f.=81, P<0.01). The F1 male progeny from both reciprocal crosses resembled their L71 male parents, suggesting that large wing size is dominant to small size. However, direct comparison of wing size across generations is inappropriate because of potential batch-to-batch variation in food and other subtle environmental differences. Hence, we confirmed this pattern using a specific male informative backcross (discussed below).

The second observation was that all chromosomes contributed to size but to a different extent (Figure 1). The F1 progeny from the two reciprocal crosses differed significantly (t=2.43, d.f.=70, P<0.05) but only slightly (0.030 mm), indicating that the Chr X from L71 had a positive influence on size but its effect was relatively small compared with the autosomes. The effects of Chr 2 and Chr 3 were larger than Chr X. Owing to the lack of meiotic recombination in male Drosophila, four male informative backcrosses—L20♀ × (L20♀ × L71♂) ♂—were used to estimate the relative contribution of each autosome, and to confirm the dominant nature of the size loci. The chromosome composition of each backcross progeny was determined by genotyping two markers (HRM analysis) per autosome: Tor and Pi3K59F for Chr 2; sav and pli for Chr 3 (Table 1). Compared with the L20 homozygous siblings, backcross progeny carrying one copy of chromosome 2 from L71 conferred a 0.068 mm increase in wing size (t=10.9, d.f.=30, P<0.01); progeny carrying one copy of chromosome 3 from L71 were 0.086 mm larger (t=13.8, d.f.=29, P<0.01); and progeny with one copy of each of chromosomes 2 and 3 from L71 gained 0.145 mm (t=23.4, d.f.=36, P<0.01). Hence, wing size was controlled mainly by loci on Chr 2 and 3, and Chr 3 had a greater contribution than Chr 2.

QTL mapping

In all, 47 gene markers covering all four autosomal arms were used for mapping analysis in the female informative crosses. Composite interval mapping revealed one major QTL peak on Chr 3L between markers L311 and L319, explaining 30.6% of the phenotypic variance, and one minor region on Chr 2R between markers L218 and L222, explaining 9.3% of the phenotypic variance (Figure 2). Collectively, these two loci accounted for 39.9% of the total phenotypic variance. No obvious epistatic interaction was detected between these two loci (Figure 3). In addition, three other suggestive QTLs were also evident on the logarithm of odds profile plot (Figure 2). These putative QTL regions accounted for 5.0, 7.3 and 6.6% of the variance. Hence, our QTL mapping captured over half of the phenotypic variance (58.8%) and the two significant QTLs contributed to wing centroid size independently.

Figure 2
figure 2

QTL analysis of wing size variation in D. simulans. QTL analysis using 47 markers covering chromosomes 2 and 3 revealed one major peak on Chr 3L and another on Chr 2R. The percentage phenotypic variance explained by a QTL is indicated above each peak. The y-axis is the logarithm (to the base 10) of odds (LOD) score. The horizontal bar at LOD=2.5 represents the experimental-wise significance threshold at P<0.05 with 10 000 permutations. Triangles on the x-axis indicate marker positions. The horizontal bar underneath each significant QTL represents a one-LOD support interval, with a black dot indicating the peak marker within each interval. The lower graph indicates the additive components along the chromosomes.

Figure 3
figure 3

Absence of epistatic interaction between QTLs on Chr 3L and Chr 2R. Mean wing centroid size in mm is plotted against the four genotypic combinations. ‘L318=1’ on the x-axis denotes the presence of the L71 allele at the L318 locus; likewise for marker L219. Error bars represent s.e. values of the mean wing centroid size.

Single-marker analysis

We tested the robustness of the QTL by genotyping additional 155 male progeny from seven female informative families for 11 QTL-linked markers (L208, L209, L218, L219, L222, L311, L318, L319, L325, L328 and L333). The sample size, therefore, increased from 52 to 207 for most markers (Table 2). Regression analysis showed that marker L318 on Chr 3L had the largest influence on the phenotype (R2=0.318) (Table 2). This was in agreement with the initial QTL scan using 52 individuals. The results were also consistent with the chromosomal substitution data (Figure 1), indicating that Chr 3 had a greater overall contribution to the wing size than Chr 2. However, the rank order of markers on Chr 2L, 3R and 2R was different from earlier results.

Table 2 Summary of single-marker analysis

Comparative analysis between D. simulans and D. melanogaster

Although the chromosome arms (Muller's elements) are highly syntenic between D. simulans and D. melanogaster, genes on each chromosome arm have undergone various degrees of local rearrangements (Clark et al., 2007). Comparison of QTL locations between the two species still requires the use of conserved genetic markers to accurately align the maps. Each gene marker developed for D. simulans in this study has an identifiable orthologue in D. melanogaster (Table 1). Figure 4 summarises the relative positions of the wing size QTLs in D. simulans and D. melanogaster. The comparison suggests that the major QTL on Chr 3L at L318 in D. simulans is at a different genomic location compared with those reported by Calboli et al. (2003) and to Dca in McKechnie et al. (2010). The minor QTL on Chr 2R at L219 does overlap with the D. melanogaster QTL region reported in Calboli et al. (2003), but not the 2R QTL in Gockel et al. (2002). This comparative analysis suggests that D. simulans shares at least one minor QTL location with D. melanogaster, but the major QTLs influencing wing size are different.

Figure 4
figure 4

Comparative mapping of wing size QTLs in D. simulans and D. melanogaster. The four homologous chromosome arms (Muller elements B, C, D and E) are aligned using orthologous markers. A major chromosomal inversion between the two species is evident on Chr 3R, but the linear order of the gene markers used is conserved on other chromosome arms. Solid bars indicate locations of wing size QTL regions. A circle represents the approximate position of each QTL peak in either species. The D. melanogaster QTL information is based on Calboli et al. (2003). The star symbol indicates the location of Dca, the candidate gene for body size variation in McKechnie et al. (2010).

Discussion

We undertook QTL mapping to understand the genetic basis of wing size variation between lines derived from cline-end populations of D. simulans. Five putative genomic regions were identified that influence wing size, which collectively account for 58.8% of the total phenotype variance. However, only two regions were statistically significant (Figure 2).

From a developmental point of view, Drosophila wing morphology is a composite trait influenced by tissue growth, proliferation, differentiation and patterning mechanisms (Neto-Silva et al., 2009). On the other hand, wing size can also be seen as a net outcome of various conflicting selective pressures acting on different parts of the life cycle. Natural selection may operate at the adult stage (Hoffmann et al., 2007), even though the phenotype itself (wing size) is developed during preadult growth stages. Hence, it is believed that a suite of genes can affect wing size.

The number of wing size QTL loci identified in D. simulans in this study is low given that wing morphology (for example, wing shape) is considered to be a typical polygenetic trait (Zimmerman et al., 2000). The low number of QTL found in D. simulans could stem from the small size of the segregating population and/or the low magnitude of the QTL effects. Owing to the inherent limitations and biases associated with QTL mapping (as discussed in Tanksley (1993)), our mapping populations would have the tendency to detect only QTL with large phenotypic effect. Given that N=52 in the initial QTL scan, genes with weak phenotypic effect are likely fall below the logarithm of odds threshold line and not detected statistically. It is also likely that the effect size of the QTL identified might have been inflated because of the Beavis's effect (Beavis 1994; Beavis 1998; Xu 2003). As a result, the R2 (phenotypic variance explained) values might not precisely reflect the true phenotypic impact of the QTL. Although the major QTL on Chr 3L is stable (Figure 2 and Table 2), larger mapping populations (N >500) will be required to more accurately estimate the phenotypic contribution and rank order of the minor QTL.

Nevertheless, our genetic architecture is comparable to those obtained by Calboli et al. (2003), who identified two major QTL regions in D. melanogaster. The major peak identified by Calboli et al. (2003) was later linked to regions within In(3R)Payne by association mapping within a population and through clinal analyses (Kennington et al., 2006; Kennington et al., 2007). Therefore, the mapping procedure used by Calboli et al. (2003) and followed here was successful in identifying genomic areas linked to geographical size variation and isolating regions with large effects on size.

Since their arrival in Australia, both D. melanogaster and D. simulans have successfully adapted to temperate regions and evolved a Bergmann's cline along the east coast (James et al., 1995; Arthur et al., 2008). The advent of whole-genome sequences and QTL mapping in both species offers an opportunity to compare their genetic bases. In the Gockel et al. (2002) study, owing to the presence of In(3R)Payne in the parental lines and insufficient generations (F3), the reduction of meiotic recombination produced a broad QTL plateau, which spans ∼85% (cytological range=65E–99F) of Chr 3. Subsequent mapping using F10 progeny by Calboli et al. (2003) confined the Chr 3 QTL interval to 3R, excluded much of 3L. Because of the higher resolution of the QTL analysis in Calboli et al. we compare our results to this study rather than the Gockel et al. (2002) study. Although potentially sharing a QTL on Chr 2R with D. melanogaster, the most significant QTL in D. simulans is located between markers L311 and L319 on Chr 3L (Figure 4). This genomic region was not implicated in previous studies in D. melanogaster by Calboli et al. (2003) or McKechnie et al. (2010). There is a potential partial overlap in genes affecting size between D. simulans and D. melanogaster, but the major locus on Chr 3 differs. In other words, phenotypic clines in these species probably have a different genetic basis, although further association mapping is required.

Despite being sibling species, D. simulans and D. melanogaster differ in many aspects (David et al., 2004). At the DNA level, interspecies divergence is also obvious based on the genome sequences (Begun et al., 2007; Clark et al., 2007). It is, therefore, perhaps not surprising that the two species can produce a similar geographical wing size divergence via alternative genetic mechanisms.

The Chr 3L QTL interval encompasses ∼3.6 M bases of DNA, containing 482 annotated genes. The Chr 2R QTL on the other hand, spans 4.9 M bases, harbouring at least 721 predicted genes. These estimates are equivalent to 4.2 and 6.3% of the total gene count (11 466) of the D. simulans genome, based on the Begun et al. (2007) annotation. We observed low representation of growth pathway genes within the two QTL regions. Activities of signalling pathways such as the Insulin-Pi3K, Tor, Hippo, in conjunction with hormonal events, are known to alter growth during larval development and hence influence final adult size, as reviewed elsewhere (Edgar 2006; Mirth and Riddiford 2007; Nijhout 2003). Key members of these cellular processes have been implicated as candidate genes for adaptive size control (De Jong and Bochdanovits 2003). The authors summarised the genomic locations of the Insulin–Pi3K pathway as well as metabolic-pathway-related genes in D. melanogaster. The majority of these genes do not fall within the two QTL intervals identified in D. simulans, practically ruling out their direct involvement in size variation in our mapping crosses.

In conclusion, genetic mapping indicates that wing size difference between D. simulans lines is controlled by at least two independent QTL, accounting for a substantial proportion (∼40%) of the size variance. Comparative analysis reveals that the QTL on Chr 3 are different in D. simulans and D. melanogaster, although the other minor QTL (on Chr 2R and 3R) might be shared between the two species. Our results suggest that parallel evolution of Bergmann's clines in different organisms could potentially be driven by distinct genetic mechanisms, but further association mapping across and within populations is required to confirm this pattern (c.f. Kennington et al., 2007). Our QTL map also provides a new starting point for fine-scale mapping and positional cloning of the size loci in D. simulans.