Drosophila melanogaster shows clinal variation along latitudinal transects on multiple continents for several phenotypes, allozyme variants, sequence variants, and chromosome inversions. Previous investigation suggests that many such clines are due to spatially varying selection rather than demographic history, but the genomic extent of such selection is unknown. To map differentiation throughout the genome, we hybridized DNA from temperate and subtropical populations to Affymetrix tiling arrays. The dense genomic sampling of variants and low level of linkage disequilibrium in D. melanogaster enabled identification of many small, differentiated regions. Many regions are differentiated in parallel in the United States and Australia, strongly supporting the idea that they are influenced by spatially varying selection. Genomic differentiation is distributed nonrandomly with respect to gene function, even in regions differentiated on only one continent, providing further evidence for the role of selection. These data provide candidate genes for phenotypes known to vary clinally and implicate interesting new processes in genotype-by-environment interactions, including chorion proteins, proteins regulating meiotic recombination and segregation, gustatory and olfactory receptors, and proteins affecting synaptic function and behavior. This portrait of differentiation provides a genomic perspective on adaptation and the maintenance of variation through spatially varying selection.

THE amount and genomic distribution of polymorphism may be influenced by genetic drift, by mutation-selection balance, and by various forms of positive selection such as spatially varying selection, heterozygote advantage, or negative frequency-dependent selection. Spatially varying selection can generate allele-frequency differences between populations in spite of gene flow and lead to local adaptation, which is of particular interest as an intermediate step between intra- and interspecies variation (Felsenstein 1976; Endler 1977; Barton 1983).

Drosophila melanogaster has been a model system for investigating the forces maintaining polymorphism for many decades. Indeed, differentiation along latitudinal clines in this species is one of the most thoroughly documented cases of spatially varying selection (Singh and Rhomberg 1987; Hale and Singh 1991). D. melanogaster originated in Africa and was introduced to Australia and the Americas in historic times (Lachaise and Silvain 2004). Multiple phenotypes, including development time (James and Partridge 1995), diapause incidence (Schmidt et al. 2005), body size (De Jong and Bochdanovits 2003), and temperature tolerance (Hoffmann et al. 2002) show genetically determined variation along latitudinal clines. Several allozyme (Oakeshott et al. 1982), DNA (Gockel et al. 2001; Sezgin et al. 2004), and inversion polymorphisms (Knibb 1982) also show clinal variation. Most putatively neutral markers show little correlation with latitude within continents, supporting the idea that clinal variation is often caused by natural selection (Hale and Singh 1991; Berry and Kreitman 1993; Long and Singh 1995; Gockel et al. 2002; Kennington et al. 2003). The observations that clinal variation is strongly associated with easily measurable phenotypes affecting fitness (Eanes 1999; Gockel et al. 2002; Calboli et al. 2003; Norry et al. 2004; Kennington et al. 2007) and that clines for a number of phenotypes and genetic variants appear to have been independently established on multiple continents (De Jong and Bochdanovits 2003) also strongly support the proposition that many clinal variants are under spatially varying selection and that the biology of temperate and tropical populations may be quite different.

Nevertheless, these data represent a small and highly biased picture of the phenotypes and genes influenced by such selection and provide a fragmentary view of adaptive polymorphism in this genetic model system. A comprehensive genomic description of clinally differentiated genes and chromosomal regions would provide an estimate of the distribution of differentiation across the genome, which is a prerequisite for estimating the fraction of the genome influenced by local adaptation. This description would also provide new insights into the biological processes under selection and stimulate investigation into the functional biology of adaptive variants. Here we report estimates of differentiation with respect to latitude at ∼3 million markers across the D. melanogaster genome. We used an array genotyping technique, “single feature polymorphism mapping,” which requires no SNP discovery phase (Winzeler et al. 1998; Borevitz et al. 2003; Turner et al. 2005; Werner et al. 2005; Gresham et al. 2006). This approach allowed us to map differentiation between northern and southern populations of D. melanogaster from the east coasts of North America and Australia, which have previously been shown to be part of latitudinal clines (Singh and Rhomberg 1987; Hale and Singh 1991; Kennington et al. 2003; Anderson et al. 2005). The resulting portrait of genomic differentiation generates candidate loci for phenotypes previously known to vary clinally and implicates interesting and unsuspected biological processes in gene-by-environment interactions.


Fly collections:

Australian flies were collected in 2004 and described in Anderson et al. (2005). Southern flies for the hybridization experiment were collected at Sorell (northern Tasmania, 41.14° S), and Miller's Orchard (southern Tasmania, 42.46° S); northern flies were collected in northern Cardwell (Queensland), Cairns, and Cooktown (15.28° S and 16.54° S). U. S. flies were collected by Paul Schmidt in 2000. Northern flies were collected at Rocky Ridge Orchards in Bowdoinham, Maine, and Biscay Farms in Walpole, Maine. Southern flies were collected at the “Fruit and Spice Park” north of Homestead, Florida, and at the source orchards for the “Robert is Here” fruit stands south of Homestead on Route 1. Additional Australian fly lines used for sequencing (latitudes 18.16° S, 29.59° S, 32.10° S, 33.57° S, and 37.39° S) were also described in Anderson et al. (2005).

DNA isolation:

We used isofemale lines from the above locations as follows: 16 lines (northern United States), 16 lines (southern United States), 17 lines (northern Australia), and 15 lines (southern Australia), pooling flash-frozen females from all lines within a region. DNA extraction of these pooled females was conducted differently for U. S. and Australian experiments. For the U. S. populations, nuclei were isolated by centrifugation, and DNA was isolated by phenol chloroform extraction followed by ethanol precipitation. For the Australian populations, DNA was extracted with QIAGEN (Valencia, CA) DNeasy kits, followed by 1:1 phenol chloroform clean up and ethanol precipitation.

Fragmentation and labeling:

DNA samples (all at a volume of 30 μl) were fragmented with a mix of DNase I (Promega, Madison, WI), One-Phor-All buffer (Amersham Biosciences), and acetylated BSA (Invitrogen, San Diego); the amounts added per sample were 4 μl 10× One-Phor-All, 0.14 μl acetylated BSA, and ∼0.085 μl of DNase per microgram of DNA. We created a master mix at these ratios, fragmented DNA, and then adjusted the amount of DNase in the master mix to optimize performance. Fragmentation was performed in a thermocycler at 37° for 16 min, 99° for 15 min, and 12° for 15 min [we found a large effect of different thermocyclers on the variance of fragment size, with a MJ Research PTC-200 performing much better than an Applied Biosystems (Foster City, CA) 2720]. Fragmentation performance was assessed by running 3 μl of fragmented DNA in a high-percentage agarose gel, with goals of a 50-bp mean fragment size, maximum intensity, and minimal variance in fragment size. When a master mix was found to perform well, all samples were fragmented with the same master mix and run in a gel together, and samples that appeared most similar were selected for labeling. Labeling was performed by adding 2 μl of biotin-N6-ddATP (Enzo) and 3 μl of RTdT enzyme (Promega) to each sample. RTdT was diluted from 30 to 15 units/μl enzyme before use by mixing RTdT enzyme, RTdT 5× buffer, and water at a ratio of 5:1:4. Labeling was accomplished in a thermocycler at 37° for 90 min, 99° for 15 min, and 12° for 5 min. Samples were stored at −20° until array hybridization. As with the fragmentation step, extreme care was taken to minimize variation between samples within an experiment: all samples were labeled from a single master mix at the same time.

Array hybridization:

We hybridized four replicate arrays per population (16 total arrays) at the Affymetrix core facility at the University of California at Davis Genome Center. The U. S. and Australian samples were hybridized differently. For the U. S. samples, the hyb cocktail was based on Affymetrix protocols for SNP genotyping: 132 μl 5 m tetramethylammonium chloride (TMAC); 2.2 μl 1 m Tris, pH 7.8; 2.2 μl 1% Tween 20; 2.2 μl 50 mg/ml BSA; 2.2 μl 10 mg/ml herring sperm DNA; 3.67 μl B2 control oligo; 60 μl labeled DNA target; and water to 220 μl total volume. The fluidics protocol used was DNAARRAY_WS4_450. For the Australian arrays, the Affymetrix cocktail for measuring a whole-transcript double-stranded target was used: 125 μl 2× hybridization mix; 17.5 μl DMSO; 4.17 B2 control oligo; 60 μl of labeled DNA target; water to 250 μl. The fluidics protocol used was FS450_0001. The 2× hybridization mix includes the BSA and herring sperm and uses morpholineethanesulfonic acid (MES) (as in the expression arrays) instead of TMAC. We recommend the protocol used for the Australian experiment because the U. S. protocol occasionally generated chips with conspicuous spatial artifacts (which were discarded). For the Australian samples, ∼7.5 μg of DNA were used per microarray, while ∼10 μg were used in the U. S. experiment (as estimated using λDNA gel standards and a Nanodrop spectrophotometer). All chips within each experiment were hybridized to equal amounts of DNA.

Data extraction:

Arrays were scanned at the Affymetrix core facility according to the manufacturer's specifications, which generates a binary file with intensities for each probe (a .cel file). These files were converted to text with a program available at the Affymetrix website. We used NCBI megablast to identify array probes with a single perfect match in version 4.3 of the D. melanogaster reference genome. Of the 6,553,600 probes on the array (which include mismatch and control oligos), 2,971,885 were retained for analysis. We used the D. melanogaster version 5 assembly, which contained significantly more assembled heterochromatic DNA than version 4.3, to identify probes that were unique perfect matches to the heterochromatic regions of chromosomes 2, 3, and X; this yielded an additional 32,256 probes. Thus, the total number of analyzed probes was 3,004,141.


The goals of normalization were to control for heterogeneous average signal intensities across chips, to control for spatially nonrandom patterns of signal intensity within chips (Borevitz et al. 2003), and to remove the correlation between mean and variance of signal intensity for each probe. To spatially normalize our data, we divided each array into 1600 subarrays of 64 × 64 probes. We log transformed raw intensity values and then divided the intensity of each oligo by the median intensity of the local 64 × 64-probe window. This 64 × 64-probe window was defined using all the probes on the array, so that all windows were the same size. However, medians were calculated using only the unique probes described above. We expect most probes to hybridize equally well in all populations. Thus, the expected proportion of probes that are more intense in either population should be ∼50%. After normalizing the arrays using the local 64 × 64 probe medians, 51.4% of probes had higher mean hybridization intensity in the southern United States. However, in Australia, 56.9% of probes were more intense in the southern United States (average difference in probe intensities were 3.1 × 10−4 in America, 3.5 × 10−3 in Australia). Therefore, we further normalized these data using quantile normalization in R (Gautier et al. 2004). After quantile normalization, 50.1 and 49.9% of probes were more intense in the southern United States and Australia, respectively, with the average difference in the means <4 × 10−10 on both continents. The heterochromatic regions from the v5 annotation were added to the analysis later; arrays were renormalized to include these additional 32,356 probes, with only the data for these new probes retained.


We estimated geographic differentiation between northern and southern populations at each probe using a two-tailed t-test on normalized signal intensities; all measurements were done separately for the Australian and U. S. experiments. The resulting P-values were converted into false discovery rates (FDR) by dividing the frequency distribution of P-values into 1000 bins and calculating a q-value for each bin (the q-value is the estimate of the FDR for any given P-value; Storey and Tibshirani 2003). Under the assumption of a uniform distribution of truly null features, the expected frequency of null features in each bin is 2,971,885 euchromatic probes × 0.001. Our q-value for each bin is simply this value divided by the number of observations in each bin (q-values >1 were considered to be 1). This q-value is different from that of Storey and Tibshirani in that it is the expected FDR of a given range of P-values, rather than the FDR of all P-values less than a given value.

With an array probe for every 40 bp of the genome, on average, differentiated regions will often span multiple probes. To detect windows of differentiated probes, we measured the average q-value along each chromosome arm in windows of 20, 50, and 100 probes, moving 10% of the window size at each measurement. We generated null distributions for these overlapping windows by permuting the locations of probes on each chromosome arm 50 times. We estimated FDR for windows by dividing the expected number of windows of a given value (from the permutations) by the observed number for each window size; each chromosome arm was treated separately. Assembled heterochromatic regions are relatively small and experience little crossing over, so linked selection could lead to differentiation of entire heterochromatic regions relative to the genomic average. Significance of the fourth chromosome and heterochromatic regions of chromosomes 2, 3, and X was therefore determined by permuting all probes in the genome to generate a genomic average. The expected degree of overlap between continents was calculated under the naive assumption of independence between base pairs and is simply the percentage of the eukaryotic genome significant in Australia (1.36%) times the percentage significant in the United States (0.43%) times the number of base pairs in the eukaryotic genome (1.17 × 106).

Copy number variants (CNVs) at different frequencies in DNA pools may generate directionality in signal intensities between northern and southern samples. For example, a window corresponding to a duplicated region that is at high frequency in the north and low frequency in the south should have greater signal intensity in the north at all probes that are performing well. We conservatively considered a window to be a candidate copy number variant if 45 probes in a 50-probe window or 80 probes in a 100-probe window all had higher means in the same population; the permutation-based 0.001 FDR threshold for this test for each chromosome arm was equal to or less than these values for all chromosomes.

We used a gene ontology (GO) database (amigo.geneontology.org/cgi-bin/amigo/go.cgi) to determine whether genes that overlap differentiated windows are nonrandomly distributed with respect to function; a gene was considered to overlap a differentiated window if any portion of its transcript was within the window. We used only terms associated with more than five genes. Significance was determined by comparing the observed proportion of genes in each GO term to an expected number, which is the binomial sampling probability of sampling an equal or larger number of genes in each category, given the number of significant genes in the genome. Because regions of the genome are significant in our analysis, and some genes with related functions are clustered in the genome, these probabilities should be interpreted with caution. To partially account for this effect, we sampled regions of the genome that exactly matched the distribution of significant regions 250 times and determined the number of genes for each GO term in these simulated data. We consider a GO term to be significant when P < 0.05 for both the binomial sampling test and the sampling-of-genomic-regions test.


To investigate how inferences made from chip data reflect underlying patterns of DNA sequence variation, we sequenced three types of regions: significantly differentiated regions, random regions, and nonsignificant regions (chosen nonrandomly). The significant regions are those that were significant in the sliding-window or single-probe analyses. The random regions were determined using a random number generator; regions were discarded if they overlapped a significant region or a TE insertion site. We also sequenced a set of regions that were not significant, but were not chosen randomly, to investigate particular aspects of the data. For example, some sequenced loci were near significant regions and others overlapped a single probe with a very low q-value that was not located within a significant window. DNA was isolated from each isofemale line, and primers were designed from the D. melanogaster genome sequence using Primer Select (Lasergene). All PCR products were cloned using TOPO TA kits (Invitrogen), and one clone from each isofemale line was sequenced in both directions. All sequences have been deposited in GenBank under accession nos. EF670724EF672031 and EU476191EU476879. See Table 4 for sequence lengths, sample sizes, and summary statistics at each locus. Sequence traces were assembled and edited in CodonCode (http://www.codoncode.com), and summary statistics were calculated in DNAsp (Rozas et al. 2003). For comparisons of probe q-values to FST, FST was estimated for each polymorphic base that overlapped the 25-bp probe using the method of Hudson, Slatkin, and Maddison (Hudson et al. 1992); probe FST values are the mean FST of all polymorphic bases within a probe. Wilcoxon tests, χ2 tests, and Spearman rank correlations were computed in R; due to the presence of ties, rank correlation P-values are not exact.

For clinal analysis, we obtained sequence data from five additional Australian population samples at Or22b, TfIIB, and dmrt93B. Five to 19 alleles per population were used; the 16.54° S and 15.28° S latitude populations were combined and considered 15.91° S for all loci, and the 32.10° S and 33.57° S latitude populations were combined and considered 32.84° S for Or22b because of low numbers of sequenced alleles from each population.

Quantitative real-time PCR:

Genomic DNA was prepared from 10 northern Australian lines and 9 southern Australian lines. DNA concentrations were estimated and adjusted such that concentrations of all samples varied by no more than twofold. Quantitative real-time PCR (qPCR) primers were designed for the putative large duplication segregating on chromosome arm 3R and tested on genomic DNA; primers for two control genes, RpL32 and Acp53, were taken from Fiumera et al. (2005). For each DNA sample, four replicate qPCR reactions were carried out for each of the three primer pairs. The four replicates were evenly split between two plates (i.e., two replicates for each DNA on each plate). We used 1 μl of DNA plus 11 μl of the SYBR master mix (Applied Biosystems: 6 μl SYBR, 4 μl dH2O, 0.5 μl 25 μm forward primer, and 0.5 μl 25 μm reverse primer). We used both no-DNA controls and dH2O controls. Samples were run using the default settings on an Applied Biosystems 7900HT. Data collection was performed in SDS version 2.1 and output was exported to Excel. The amount of putative copy-number DNA relative to the two control genes was determined using the 2−ΔΔCT method (Livak and Schmittgen 2001). Significance of relative amounts of putative copy-number variant DNA between northern and southern Australian population samples was determined using the Wilcoxon test.


We used Affymetrix D. melanogaster tiling arrays, which have a 25-bp probe for approximately every 40 bp of the genome, to identify regions showing genetic differentiation between northern and southern D. melanogaster populations in the United States and Australia. Because the hybridization intensity of genomic DNA to a microarray depends on sequence similarity, differentiated loci can be mapped through array hybridization (Winzeler et al. 1998; Borevitz et al. 2003; Gresham et al. 2006). If a polymorphism with different frequencies in two populations overlaps a probe on the array, DNA derived from these populations may hybridize to the array with different affinities (Turner et al. 2005). Although this approach has several limitations, including variable sensitivity across the genome and, potentially, nonlinear change in hybridization intensity with increasing sequence mismatches (Zhang et al. 2003), it provides a useful way of measuring differentiation at a very large number of markers. We isolated DNA from pools of ≥15 isofemale lines from Maine (northeastern United States), Florida (southeastern United States), Queensland (northeastern Australia), and Tasmania (southeastern Australia) and hybridized each pool to four replicate arrays. We first present an analysis of differentiation at individual probes from the tiling array, followed by sliding-window analysis. Our goal here is not to formally test a population genetics model, but rather to quantify the observed distribution of differentiation across the genome. This allows us to identify the most unusual regions in the genome from populations near the endpoints of known clines and to generate hypotheses regarding genes and phenotypes harboring adaptive variation.

Differentiated probes:

We used a t-test to measure differentiation between temperate and subtropical populations on each continent separately for the 2,971,885 array probes with single perfect matches to the D. melanogaster reference genome. The resulting distribution of P-values is shown in Figure 1A. On the basis of the excess of probes with low P-values, we can estimate the number of probes that overlap differentiated DNA (Figure 1A and materials and methods). Assuming a uniform distribution of truly null features, 316,123 and 189,447 euchromatic probes overlap detectable differentiation in Australia and in the United States, respectively. The greater skew in the distribution of P-values in Australia suggests that the average differentiated probe is more differentiated between northern and southern Australian than between northern and southern U. S. populations (Figure 1A). For the majority of further analyses, we converted these P-values into q-values, which provide estimates of the false discovery rates for each range of values (see materials and methods and Figure 1B).

Figure 1.—

Probe differentiation. (A) Histogram of t-test P-values. Only P-values <0.50 are shown; the remainder of the distribution is essentially flat. Red, Australian data; black, U. S. data. (B) The relationship of P-value to q-value, with Australian data in red. (C) Mean probe FST vs. probe q-value. (D) Proportion of probes with FST of zero vs. probe q-value; black triangles, U. S. data; red circles, Australian data.

Probes that are differentiated on both continents are excellent candidates for loci under spatially varying selection. There are 13,747 euchromatic probes in Australia with P-values <0.001 (FDR = 21.3%) and 5838 euchromatic probes with P < 0.001 in the United States (FDR = 50.7%); 104 of these probes are shared (P < 0.001 on both continents; Figure 1A). Although the false discovery rate of these probes is high in the U. S. analysis, the combined observation of P < 0.001 on both continents should increase the likelihood that they are truly differentiated. The vast majority of these probes (hereafter, “golden probes”) are not located in regions of reduced crossing over and are much farther apart than the scale of linkage disequilibrium in D. melanogaster [(Andolfatto and Wall 2003) average is 911 kb for probes on the same chromosome arm, with only one pair of probes within 500 bp of each other]. This suggests that golden probes mark many independent observations of small regions of differentiation. If differentiation at these loci was due to stochastic demographic forces or sampling variance, alternative alleles should not show parallel intensity differences with respect to climate on both continents. For example, a probe that is significantly differentiated by chance on both continents and has lower hybridization intensity in Maine should have lower hybridization in Queensland or Tasmania with equal probability. In fact, 80% of the golden probes have hybridization patterns that indicate that mismatched allele frequencies vary in a consistent manner with respect to latitude on the two continents (sign test P < 0.001). We sequenced 10 amplicons containing a golden probe: 8 in which hybridization intensity varied in parallel, and 2 in which we expected allele-frequency variation in opposite directions on the two continents with respect to latitude. The sequence data from these amplicons (see below) supported the array-based inferences for all 10 probes, with no apparent false discoveries.

Each of the four major autosomal arms in D. melanogaster segregates a paracentric inversion found throughout most of the species range. These inversions are more common in tropical populations (Knibb 1982), with the major chromosome arm 3R inversion, In(3R)P, showing particularly strong clinal patterns in Australia (Anderson et al. 2005; Kennington et al. 2006). The 104 golden probes are significantly overrepresented in regions spanned by the four cosmopolitan inversions (48 observed, 30 expected, χ2 = 4.15, P = 0.042). This pattern is driven entirely by In(3R)P (26 observed, 12 expected, χ2 = 5.16, P = 0.023), as the golden probes are not overrepresented in inverted regions when In(3R)P is omitted from the analysis (22 expected, 22 observed). Thus, although some of these 104 probes may be differentiated because of linkage to In(3R)P, the majority are likely located near a site under spatially varying selection. A list of the locations and corresponding annotations for these 104 probes is available as supplemental data.

Correlation between array data and DNA sequence differentiation:

To determine how probe q-values relate to DNA differentiation, we estimated FST between northern and southern populations for 452 25-bp probes (from 17 PCR amplicons) in the United States and for 572 probes (from 23 PCR amplicons) in Australia. Golden probe loci were not included in this analysis because they should have better-than-expected FDR due to the replicated observations of differentiation on two continents.

We found a strong, highly significant correlation between probe q-value and FST: Spearman's rank correlations are −0.306 in Australia (P = 6.8 e-14) and −0.255 in the United States (P = 3.8 e-8). The average FST increases with lower q-value (Figure 1C), and the proportion of probes with an FST of zero decreases with lower q-value (Figure 1D). These data clearly indicate that sequence differentiation between populations is revealed by these array experiments.

The array data indicate that we should observe more differentiation in Australia than in the United States (because of the greater excess of low P-values; Figure 1A). Our sequence data show the expected pattern, as the proportion of differentiated probes and the magnitude of differentiation of the most significant probes are far greater in Australia. Because our U. S. and Australian DNA samples were hybridized to the arrays under slightly different conditions (see materials and methods), comparisons between continents only on the basis of the array experiment should be interpreted cautiously. Nevertheless, the correlation between q-value and FST is similar in the two experiments, supporting the proposition that differentiation is greater between the Australian populations than between the U. S. populations. This could be due, in part, to the broader latitudinal range spanned by our Australian sample (15° S–43° S) than by our U. S. sample (25° N–43° N). Steeper ecological gradients in Australia, greater effective gene flow in the United States, or differences in the genetic composition of the introduced populations on the two continents could also contribute to different patterns of clinal variation on the two continents.

Sliding-window analysis:

Although the probe-by-probe analysis reveals many differentiated loci, limitations of the array data motivate analyses that combine information from multiple, adjacent probes. This should dramatically increase power because the scale of linkage disequilibrium in D. melanogaster is often >25 bp (Haddrill et al. 2005; Ometto et al. 2005). We used a sliding-window analysis to identify genomic regions associated with lower-than-expected q-values, given the distribution of q-values on each chromosome arm (see materials and methods). Because the physical scale and magnitude of genetic differentiation is unknown, we used two significance thresholds to generate comprehensive (FDR < 0.05) and more conservative (FDR < 0.001) lists of significantly differentiated 20-, 50-, and 100-probe windows (which span ∼800, 2000, and 4000 bp on average). Most significant 20- and 50-probe windows were contained within significant 100-probe windows. For example, 89% of significant DNA at FDR < 0.001 in Australia is significant in the 100-probe windows, with the remaining 11% composed of small, highly differentiated regions, which were not detectable with the larger window size. The union of significant overlapping windows of all sizes generated a list of nonoverlapping differentiated regions; regions composed of windows significant at FDR < 0.05 and FDR < 0.001 are hereafter referred to as regions significant at these thresholds. These thresholds merely identify outliers rather than regions having excess differentiation relative to the expectation under a population genetics model.

The window analysis confirms that a much larger proportion of the genome is significant in Australia: U. S. populations have only 31% as many base pairs in significant regions at both FDR thresholds. Although many regions are significant on only one continent (typically Australia), significant regions are shared between continents much more than would be expected by chance (150 kb are significant on both continents at FDR < 0.001, while only 8.7 kb are expected; χ2 = 125,807, P < 2.2e−16). Although this expected value is probably an underestimate of the true null hypothesis (because power varies across the genome and base pairs are not independent), this analysis provides genomic support for previous conclusions that many genes vary clinally on multiple continents and are likely to be experiencing spatially varying selection (Singh and Rhomberg 1987; Hale and Singh 1991).

Table 1 reports the number of base pairs in significant regions and the degree of overlap of significant regions between continents for each chromosome arm. Much of the heterogeneity between chromosome arms is attributable to a few large regions. For example, large centromere-proximal regions of the X chromosome appear to be differentiated in Australia but not in the United States. Results from the window analysis support results from the probe-level analysis in that significant regions were overrepresented on chromosome arm 3R. At the conservative FDR threshold, nearly half (47.3%) of the differentiated regions significant on both continents are on 3R (70,962 bp observed, 35,752 expected; χ2 = 28,271, P < 2.2e−16). An example of the window analysis (50-probe windows) for 3R can be seen in Figure 2. The breakpoints of the known clinal inversion In(3R)P are highly significant on both continents. Overall, however, significant regions are not dramatically more abundant within the inverted region (Figure 2), consistent with previous reports suggesting that levels of linkage disequilibrium in D. melanogaster are not consistently high throughout In(3R)P (Kennington et al. 2006). The only other inversion for which precise breakpoints are known is In(2L)t (Andolfatto and Kreitman 2000). Regions near the breakpoints of this inversion are not significantly differentiated in our analysis, perhaps because the cline for In(2L)t is much shallower than the cline for In(3R)P. We detected numerous, previously known clinal loci in our window analysis, including Est6 in Australia (Coppin et al. 2007), Hex-C in Australia (Eanes 1999), Gdh in the United States (Eanes 1999), hsr-omega in Australia (Anderson et al. 2003), UGPase in the United States (Sezgin et al. 2004), cpo on both continents (P. Schmidt, personal communication), and period in Australia (although this is a portion of this gene that is not known to be clinal; Sawyer et al. 1997) in addition to the In(3R)P breakpoints on both continents.

Figure 2.—

Summary of chromosome arm 3R. Windows of 50 probes significant at FDR < 0.001 are shown in red (Australia) and black (United States). Each circle is a window; significant regions appear as vertical lines because they are small relative to the scale of a whole chromosome (with the exception of the large ace region). The breakpoints of In(3R)P are indicated, and genes that overlap some regions are shown in black type for regions significant on both continents and in red type for regions significant in only Australia. All labeled regions are also significant in the copy number variant analysis, except the dmrt93b region. Locations of the most significant shared probes are shown as purple triangles.

View this table:

Number of base pairs and percentage of each chromosome arm in differentiated regions

The window analysis is complementary to the single-probe analysis. While the latter may detect very short, strongly differentiated regions, the window analysis combines information from multiple probes, thereby providing more power to infer differentiation. In spite of this dual approach, inherent limitations of array performance mean that many differentiated regions will not be detected. An illustrative example is the Adh locus, at which clinal differentiation is limited to a single coding SNP and noncoding indel mutation (Berry and Kreitman 1993). Indeed, we do not detect Adh in our window analysis. Furthermore, we would be unlikely to detect these two mutations in our single-probe analysis because both are 1 bp away from a probe on the tiling array. Other known clinal loci were missed as well, such as Pgm, Gpdh, Treh (Sezgin et al. 2004), and hsp70 (Bettencourt et al. 2002). Although additional arrays could increase power, resequencing efforts will be required for complete descriptions of geographic variation in the species.

The physical extent of differentiated regions linked to sites experiencing spatially varying selection depends on several population genetics parameters, including recombination rate, age of alleles, strength of selection, and amount of gene flow. The vast majority of differentiated regions that we are able to detect with this method are several hundred to several thousand base pairs long (Table 2). For regions containing at least one gene, the mean number of genes per significant region is fewer than two, so these data offer considerable opportunity to identify genes that are targets of spatially varying selection. Additionally, a number of regions that are significant on both continents contain no genes (n = 7 at FDR < 0.001) and are candidates for selection on noncoding sequences (Table 2). However, technical factors such as probe density and probe sensitivity may also affect estimates of the size of differentiated regions, and sequence data will be required for more accurate estimates. A complete list of significantly differentiated regions is available as supplemental data.

View this table:

Summary statistics of differentiated regions

Despite the small average size of most significant regions, several are many kilobases in length. For example, seven regions are >40 kb at the FDR <0.05 threshold. Two of these regions, together totaling >250 kb, are nearly contiguous and are located near the centromere of the X chromosome in Australia; another large region differentiated on both continents is near the centromere on 2R. The large size of these centromere-proximal regions may result from extensive linkage disequilibrium associated with low recombination rates. However, three differentiated regions >40 kb are not located near centromeres. One lies on chromosome 3R in Australia [outside In(3R)P] and contains acetylcholine esterase (ace), part of timeout, and seven unnamed genes. Two additional large regions, each of which contains a cluster of chorion protein-coding genes, are found in Australia. These two regions contain 9 of the 13 chorion protein-coding genes in the D. melanogaster genome (as well as 14 other genes), suggesting strong selection on this gene family. Although the ace region is also significant in the U. S. analysis, the aforementioned chorion clusters are not, again illustrating that some of the most significant regions in Australia are not significant in the United States. Large differentiated regions in areas of normal recombination may often be due to differentiated copy-number variation, which is discussed below.

Sequenced regions:

To investigate how inferences made from array data reflect underlying patterns of DNA sequence variation, we sequenced many loci from a subset of the lines used in the array experiment. Golden probes, significantly differentiated regions from the window analysis, randomly chosen regions, and other regions that were chosen nonrandomly were all investigated. This final category includes regions that contained a probe with a low q-value on only one continent, regions located near (but not within) significant windows, and regions that were sequenced on both continents despite the fact that significant differentiation was observed on only one continent. For three regions, we extended our sampling to additional Australian populations to confirm that differentiation between cline endpoints is consistent with clinal variation.

Sequence data from 10 amplicons containing a golden probe revealed differentiation overlapping the golden probe on both continents in 9/10 cases (Table 3). In the final instance (fas1), the sequenced amplicon contained a single common SNP. It was differentiated in the expected direction, but was located 22 bp away from the golden probe. We suspect that this SNP may have affected hybridization of the array probe through “proximity” effects by altering DNase fragmentation or the behavior of fragmented DNA during hybridization. Most of these loci are located in noncoding regions: differentiation may be easier to detect in these regions due to higher polymorphism, which leads to multiple differentiated SNPs or indels within the probe. An example coding region, Gr93c, is shown in Figure 3. A clear peak of strong differentiation located around the golden probe contains nine SNPs, three of which are nonsynonymous. These nonsynonymous SNPs are among the most differentiated sites in the gene and show parallel temperate vs. tropical frequency differences in the United States and Australia. Such SNPs are attractive candidates for functional investigation. Investigation of differentiated noncoding variants may provide insight into selection on gene regulation. For example, parallel differentiation at intronic sites in cp309, a centrosome component required for sperm motility (Kawaguchi and Zheng 2004; Martinez-Campos et al. 2004); at intronic sites at CG13466, a predicted microtubule-related protein; and at an intergenic region near spn-D, which functions in female meiosis (Abdu et al. 2003), provide opportunities for investigating the functional and population genetics consequences of regulatory variation.

Figure 3.—

Differentiation at Gr93c. Each SNP with FST > 0 is shown as a circle, with the three nonsynonymous SNPs as solid circles; lines are running averages of FST in 50-bp windows (Australia in red, United States in black). Thin vertical lines delineate the golden probe site. The two exons of Gr93c are shown as boxes.

View this table:

Frequencies of sequence variants at golden probe loci

Table 4 shows summary statistics for all 41 amplicons sequenced from Australia and all 32 amplicons sequenced from the United States (many of these regions were also used to analyze the correlation between FST and probe q-value, but additional information is gained by analyzing these data at the scale of amplicons). The randomly chosen loci provide an unbiased sample of population genetic variation in the United States and Australia (although with only seven loci, it is limited in scope). These loci show low differentiation between northern and southern populations on each continent (average FST = 0.054 in Australia, 0.036 in the United States) and little differentiation between continents, with the exception of R10. This X-linked locus shows very little differentiation within continents (Table 4), but remarkably, it has 9 fixed differences between continents (n = 19 for each continent; FST = 0.782). In agreement with previous investigations (Hale and Singh 1991), these seven random loci show no more between-continent differentiation for temperate than for tropical comparisons (average between-continent tropical vs. temperate FST = 0.147; average between-continent tropical vs. tropical or temperate vs. temperate FST = 0.142; Wilcoxon P = 0.66), which supports the hypothesis that the clines were established independently.

View this table:

DNA polymorphism and population differentiation in Australia and the United States

The conventional view of clinal variation is that it often represents the spread of a derived temperate-adapted allele to high frequency in populations farther from the equator (Sezgin et al. 2004). Indeed, we observed several possible cases of this phenomenon, as indicated by much lower heterozygosity in temperate populations (e.g., Or22, TfIIB, CG13466, and α-tubulin 67C). In Australia, polymorphism is lower in temperate populations in significant regions but not in random regions (Wilcoxon P = 0.003, 0.018 for π and θ, respectively, for significant regions). However, in many cases of strong differentiation, there is little difference in heterozygosity between temperate and tropical populations, and other significant regions have lower polymorphism in tropical populations. Further analysis regarding the histories of putatively selected alleles in such regions will be required to illuminate this issue. Finally, in Australia, the frequency spectrum is skewed toward rare alleles in temperate populations for significant but not for random regions (Wilcoxon P = 0.0002 for Tajima's D in significant regions). Differentiated regions in the United States show fewer systematic differences in heterozygosity or the frequency spectrum, but the trend in the U. S. data for the frequency spectrum agrees with the trend observed in Australia, with all golden probe loci having lower Tajima's D in Maine. There are individual loci in the United States that show major differences between northern and southern populations, such as the 5′-UTR of star1, which has a very negative skew in the frequency spectrum and a reduction in haplotype diversity in Florida, while CG6947 shows a negatively skewed frequency spectrum and reduced polymorphism in Maine (this locus was sequenced because it contained one of the most significant probes in the genome).

As was true for the single-probe analysis, the amplicon data reveal a good correspondence between array-based and sequence-based estimates of differentiation. However, several significant regions had low FST. In some cases, this is because the differentiated region is smaller than the size of the amplicon. An example of one such narrow region of differentiation is shown in Figure 4; this region surrounds the first exon of dmrt93B, which is involved in sex determination (Zhu et al. 2000). Many of the golden probe loci have low (or zero) average FST across the sequenced amplicon, but still contain a differentiated SNP at the differentiated probe site: these loci are very attractive candidates for follow-up study, contrary to what one would be led to believe on the basis of the amplicon FST estimates in Table 4. Despite these successes, our sequencing revealed several dramatic false discoveries (such as DmsR-2 in Australia). Although we expected <1 in 1000 of the regions significant at FDR < 0.001 to be false discoveries in the technical sense, some regions may also be false discoveries in the biological sense. For example, partial homology with differentiated DNA elsewhere in the genome may affect probe intensity mapping to some loci. Two apparent false discoveries, CG7738 and CG9153 (Table 4), are unusual compared to the other significant sequenced regions in that they contain no probes in the lowest q-value categories. Rather, they are significantly differentiated because they possess an unusual concentration of intermediate q-values (the most significant probe overlapping the CG7738 amplicon has a q-value of 0.68, but only 5 of 22 probes overlapping this amplicon have a q-value of 1.00). A possible explanation for this result is that our array experiment (which provides estimates of frequency differences based on a large pooled sample) detected many slightly differentiated SNPs at these loci, but this differentiation is too slight to be detected in our small sample of sequenced alleles.

Figure 4.—

Differentiation at dmrt93b. Microarray differentiation is shown for windows of 20 probes. Australia, in red, is more differentiated; the United States is shown in black. Sequence differentiation (FST) is shown as triangles for the United States and circles for Australia (each point corresponds to a 200-bp sliding window, moving 100 bp at a time). Differentiation overlaps the first exon of dmrt93b (box) and part of the first intron.

On the basis of the observed relationship between amplicon FST estimates and q-values, we designed an additional ad hoc analysis, which identifies highly differentiated loci with fewer false discoveries. We made all probe q-values >0.6 equal to 1.0, such that our window analysis would only detect clustering of the most differentiated probes. This modification significantly increased the Spearman's rank correlation between amplicon FST and average q-value in the U. S. data: considering all q-values, ρ = −0.33, P = 0.17; considering q > 0.6 = 1.0, ρ = −0.59, P = 0.008. The Australia correlation also increased: considering all q-values, ρ = −0.38, P = 0.046; considering q > 0.6 = 1.0, ρ = −0.46, P = 0.014. Thus, this may be a good approach for finding the most differentiated regions with the lowest false discovery rate. Alternatively, including all q-values may be the most powerful way to delineate all significantly differentiated regions for genomic analyses such as the gene ontology analysis presented below.

Despite the demonstrated usefulness of these data, which vastly increase the number of potentially clinal loci, we emphasize that the FDR estimates for regions are inaccurate, in contrast to the FDR estimates for individual probes, which appear to be very accurate. We therefore consider our three significance thresholds for regions (FDR < 0.05, FDR < 0.001, and the most conservative analysis just discussed) to provide three relative levels of confidence, which can be used by the community for exploratory analysis of clinal differentiation. An example of how these data may be used to gather evidence bearing on spatially varying selection of genes or classes of genes can be found in the accompanying article analyzing variation in chromatin-remodeling genes (Levine and Begun 2008, this issue). Results from all thresholds from the array analyses are available as supplemental data.

Clinal samples:

Although the array analysis included only populations from the endpoints of known clines, we expect that most of the differentiated regions detected between these populations will prove to be clinal. The presence of known clinal loci in our data supports this assertion. Previous investigations of cosmopolitan D. melanogaster populations, especially in Australia, have found little other differentiation between populations [although one previous study found that Tasmanian populations were differentiated from southern mainland populations, additional studies did not confirm this result (Agis and Schlotterer 2001; Kennington et al. 2003)]. To investigate this issue, we sequenced three differentiated Australian loci (Or22b, TfIIB, and dmrt93B) in Australian populations from additional latitudes (see materials and methods). All 11 common polymorphisms in the Or22b locus and 37/43 common polymorphisms in the TfIIB locus are significantly (P < 0.05) clinal, with r2 values up to 0.996. Of the 20 polymorphisms falling within the most differentiated window of dmrt93B (Figure 4), 16 had significant r2 values (P < 0.05). When the allele frequencies for all polymorphisms in each locus are averaged, all three loci remain significantly clinal (Figure 5). We did not detect a break between the Tasmanian populations used in the array analysis and the other southern Australian populations. These data support the notion that many of the highly differentiated regions identified in the array analysis harbor clinal variation. The analysis presented in the accompanying article (Levine and Begun 2008) also supports the hypothesis that strongly differentiated regions are often clinal.

Figure 5.—

Clinal differentiation. We sequenced Or22b (circles), TfIIB (squares), and dmrt93B (triangles) from additional latitudes to determine if differentiated regions were clinal. The average allele frequency for all polymorphisms in each locus is shown; multiple polymorphisms in complete linkage were considered once when taking the average. The r2 values for these averages vs. latitude are 0.982, 0.907, and 0.868 for Or22b, tfIIB, and dmrt93B, respectively (all P < 0.005).

Differentiation of functional categories:

Given that most significant windows are associated with one or two genes, functional annotations corresponding to genes in significant windows can be used to investigate whether certain GO terms are associated with differentiated regions (see materials and methods). The results are summarized in Table 5, while the full list of significant GO terms is available as supplemental data. Genes are not distributed randomly along chromosomes with respect to function, so that physical clustering of related genes could inflate the observed significance of GO terms. However, all of the GO terms in Table 5 include genes from different locations in the genome, which supports the hypothesis that these functions are under selection. We also assessed significance of GO terms using two complementary strategies, one of which explicitly considered the clustering of related genes (see materials and methods). In this section, we highlight several examples to illustrate the nonrandom nature of differentiation with respect to gene function.

View this table:

Summary of significant GO terms

Regions that are differentiated on both continents contain a dramatic overrepresentation of several gene families, including UDP glucotransferases and cytochrome P450s (cyps). A region on 3R containing five genes with glucuronosyltransferase activity and two unlinked related genes (CG10170 and CG30438) is significant on both continents, leading to a highly significant overrepresentation of genes with glucuronosyltransferase activity, polysaccharide metabolism, and toxin response (Table 5). Cytochrome P450 (cyp) genes are overrepresented in shared significant regions as well: Cyp12a4 and Cyp12a5 are in one of the most significant regions of the genome (Figure 2). Five cyp genes are significant in additional shared regions: Cyp6a17, Cyp6a22, and Cyp6a23 are in a significant cluster on 2R, Cyp6g1 is also on 2R, and Cyp6t1 is on the X (Cyp6a17 and an additional gene on 2R, Cyp12d1, each contain 1 of the 104 golden probes). These cytochromes and transferases lead to the significant overrepresentation of genes involved in steroid metabolism and defense response (which also includes the significant genes hsc70-2, cactin, and turandot C).

The nonrandom distribution of gene functions in regions significant only in Australia provides evidence that these functions are under spatially varying selection on that continent. For example, GO terms associated with the chorion (the outer layer of the eggshell) are overrepresented in Australia. The six chorion proteins that are significant in the 3L and X chorion clusters in Australia represent 75% of the genes in the “structural constituent of chorion” GO term. Additionally, fs(1)N, which is known to function in chorion formation, is differentiated in Australia, contributing to the significant overrepresentation of the “insect chorion formation” GO category. Three other genes, mip130, rbf, and dup, are all involved in chorion gene amplification and are significant in Australia. Finally, three additional chorion genes, Cp7Fa, Cp7Fb, and Cp7Fc, are significant in the X-linked chorion cluster but are not associated with any GO terms. Interestingly, none of these 13 chorion-related genes overlaps a significantly differentiated region in the United States. Because the chorion clusters are polytenized in some tissues in females, it is possible that these regions are significant because of geographic variation in the degree of polytenization in females or different proportions of polytenized tissues in these females rather than from different allele frequencies at these genes.

Of the nine known acetylcholine receptors in the genome, portions of three (in addition to ace itself) are significant on both continents and an additional three are significant only in Australia. These data, and the fact that many other significant genes are involved in synaptic transmission, lead to a significant overrepresentation of a suite of synapse-related GO terms. The known examples of insecticide resistance alleles associated with acetylcholine metabolism in flies and mosquitoes (Fournier et al. 1989; Weill et al. 2002, 2003) and the aforementioned significant Cyp6g1, which is also implicated in insecticide resistance (Chung et al. 2007; Daborn et al. 2007), suggest that selection pressures associated with insecticide or other environmental toxins vary with latitude.

In the United States, the cell cycle/DNA damage checkpoint genes mei-41 and grapes (Larocque et al. 2007) were both significant; the genes are unlinked, and differentiation of both is unlikely to be due to chance (P < 0.04 in GO analysis). The differentiated region of mei-41 included only one nonsynonymous variant; the ancestral allele was present in all other sequenced Drosophila genomes, which supports the idea that the derived allele is functionally important. These genes were not significant in Australia, and sequencing verified that they are differentiated only in the United States. Differentiated genes thought to play a role in double-strand break repair, recombination, or chromosome segregation in Australia include α-tubulin 67C (Matthies et al. 1999), egl (Carpenter 1994), tefu (Oikemus et al. 2006), and glu (Steffensen et al. 2001), with α-tubulin 67C and tefu among the most differentiated loci sequenced (glu sequence data are presented in Levine and Begun 2008). The verified golden probe near spn-D could be related to these functions as well (Abdu et al. 2003). In general, clinal variation could plausibly be associated with temperature variation (directly or indirectly) or with other ecological variables correlated with latitude. Alternatively, the genetics of D. melanogaster populations may vary with latitude as a result of clines for chromosome inversions, which are major modifiers of recombination that could result in selection on genes that regulate crossing over. In tropical populations, autosomes are frequently heterozygous for paracentric chromosome inversions, which should lead to greater frequencies of achiasmate meioses and higher rates of nondisjunction (for autosomes). Indeed, females heterozygous for inversions on all four autosomal chromosome arms suffer a considerable drop in fertility (Stalker 1976). In the absence of data on phenotypic variation associated with these genes, several hypotheses for selection are plausible; at the present time, we can only speculate.

The wide range of overrepresented functions illustrates the pervasive effect of the environment on D. melanogaster biology and clearly indicates that functional analysis of natural variants can contribute to annotation of the Drosophila genome. Differentiation of sets of related genes or genes that interact with one another provide an attractive opportunity to study recent adaptation in a network context.

Copy number variation:

The three large differentiated regions not located near centromeres (the ace region and the two chorion clusters; see above) have an interesting feature in common: most probes in each region have higher mean intensities in one population. For example, the chorion cluster on 3L contains 1106 probes in the 58-kb significant region: 91% of these probes have higher mean intensity in northern Australia than in southern Australia. The simplest explanation for this result is copy number variation—either a duplication with respect to the reference sequence is present at higher frequency in northern Australia or a deletion with respect to the reference sequence is at higher frequency in southern Australia. We scanned the genome for CNVs by recording the sign of the intensity difference for each probe and comparing this distribution to a null expectation from permuted data (see materials and methods). We found many differentiated regions in which most probes (45 probes in a 50-probe window or 80 probes in a 100-probe window; see materials and methods) had a higher mean intensity in one population, consistent with copy number variation. In the United States, 26 (4.7%) regions significant at the FDR < 0.05 threshold met this stringent criterion. In Australia, 22.3% of our significantly differentiated regions were also significant in this copy number variant analysis. However, other biological or technical factors could generate autocorrelation in mean probe intensity along chromosomes. Consider, for example, a region with many differentiated SNPs and strong linkage disequilibrium. If the reference strain matches a haplotype occurring at higher frequency in the south, then all SNP-containing probes would be more intense in the south. Extrachromosomal DNA or differential degrees of chromosome polytenization in two populations could also generate autocorrelated patterns of probe intensity. For probes not containing a SNP, technical factors might create slight differences in mean intensity between populations that are correlated along chromosomes. For example, the nonrandom spatial distribution of GC content across Affymetrix arrays, combined with the existence of spatial artifacts not entirely corrected in our normalization, could result in nondifferentiated probes with spatially correlated intensities (for nondifferentiated probes, all variance is technical variance, so even if this variance is very small it will affect this test). Despite these caveats, it seems very likely that some fraction of our candidate CNVs is real and that the number of geographically differentiated CNVs in the genome is considerable.

To begin investigating these candidate CNVs, we used qPCR (see materials and methods) to confirm that the largest differentiated candidate CNV, which spans nine genes and >40 kb, is a polymorphic duplication at higher frequency in northern Australia than in southern Australia (P = 0.01). This includes ace, which is also contained within a candidate CNV in the United States. However, the candidate U. S. CNV does not include the other genes, suggesting the possibility that Australian and U. S. CNVs originated independently (Labbe et al. 2007). We also validated, using PCR, that the single most significant euchromatic region in Australia is a 2098-bp deletion (Figure 6). This deletion, which includes portions of two olfactory receptors, Or22a and Or22b, is present in 9 of 12 northern chromosomes sequenced but in 0 of 11 southern chromosomes sequenced. It may have resulted from unequal crossing over between the 5′ regions of these two genes, which are 95% identical. In deleted chromosomes, the first exon and intron of Or22a remain, along with exons 2–4 of Or22b. The deleted chromosomes therefore contain a potentially functional chimeric receptor. In addition to the deletion, there are five nonsynonymous SNPs, all of which are polymorphic in the north and monomorphic in the south; several are nearly fixed differences, with frequencies of 0% in the south and 95% in the north. Finally, two SNPs and a 16-bp insertion in the 3′-UTR are present in 85–90% of northern chromosomes but absent from southern chromosomes. Or22a and Or22b are coexpressed in the ab3A neuron, which is sexually dimorphic and responds to a broad range of odors (Dobritsa et al. 2003). Functional analyses indicate that Or22a accounts for the full odor response of this neuron, with no obvious role of Or22b (Dobritsa et al. 2003). Functional investigation of the differentiated alleles detected in this experiment may shed light on the potential selective forces involved in creating such strong differentiation at this locus in Australia.

Figure 6.—

Differentiation at Or22a and Or22b. Australian data are in red, the U S. data are in gray, and the genes are boxes. Circles and dots denote individual probes, and the running averages of 50-probe windows are also shown. We verified the presence of a deletion that includes all but the first exon of Or22a and only the first exon of Or22b; the breakpoints of this deletion are indicated by vertical black lines. Most probes have q-values of 1, except in the deletion in Australia, where most probes have the lowest q-value, 0.213.

Several candidate CNVs overlap previously identified CNVs: mth8 overlaps a polymorphic deletion (A. D. Kern and D. J. Begun, unpublished results) and the differentiated CNV shared between continents at cyp6a17 is coincident with a known polymorphic duplication (C. Robin, personal communication). The breakpoints of In(3R)P are candidate CNVs in our data, and these regions are known to be associated with duplications and deletions (Matzkin et al. 2005). These data support the idea that many CNVs are influenced by spatially varying selection in D. melanogaster. A complete list of candidate CNVs is available as supplemental data.


We also analyzed heterochromatic regions of the version 5.0 genome assembly. These include the assembled portions of the chromosome 2, 3, and X centromeres, which are divided into right and left arms for chromosomes 2 and 3 (e.g., 2Rhet, 2Lhet), and the entire fourth “dot” chromosome. These regions are large, gene poor, located near centromeres, experience little crossing over, and harbor low levels of sequence variation within populations. However, despite the low levels of variation, our data suggest that these regions contain functionally important variants at different frequency in northern and southern populations. 2Rhet, 3Lhet, and 3Rhet are very differentiated on both continents (Table 1). The fourth chromosome shows extreme differentiation in Australia, with 62% of the chromosome significant at 0.001 FDR, while the United States shows no evidence of differentiation on the chromosome. Overall, these data indicate that parts of the centromeres of chromosomes 2 and 3 are very differentiated between temperate and subtropical populations, while conclusions regarding the centromere of the X are less clear, with only one small significant region in Australia (Table 1). Although we used only probes with one perfect match in the reference genome, these probes may still be repetitive in nature, either because they are misassembled in the current reference or because the reference sequence happens to possess only one copy. Indeed, all of the significant regions in the centromeric heterochromatin are also significant as candidate copy number variants (although it is even more likely that we could mistake differentiated SNP haplotypes for CNVs in regions with high linkage disequilibrium). Functionally important heterochromatic variants could be associated with heterochromatic genes. Alternatively, the variation documented here could reflect functional differences among natural D. melanogaster centromeres, perhaps associated with variation in meiotic segregation (Novitski 1955).


The list of differentiated probes and genomic regions presented here adds considerably to the list of candidate genes harboring adaptive polymorphisms in D. melanogaster. Most regions that are differentiated in parallel on both continents are likely influenced by spatially varying selection. The dramatic overrepresentation of some functionally related groups of genes in differentiated regions implicates selection acting on these loci and functions as well. These data provide an excellent opportunity to study how sets of related genes adapt to the environment (including their interactions with each other). They also provide information that may be useful in refining the annotations of known elements and dissecting the functions of unknown elements in the genome. In addition to considerable population genetic and functional analysis motivated by the array data, much descriptive work remains to be done. For example, a substantial portion of the differentiation between populations may be due to copy number variants, but we did not fully characterize this variation. We are also uncertain about whether differences between continents reflect biological differences (different selection pressures or different adaptive responses to the same selection pressures) or differences in population genetics parameters (e.g., the physical scale of linkage disequilibrium) that make regions detectable on only one continent. This is an especially interesting question, as some of our most dramatic differentiation was found only in Australia (e.g., Or22, chorion genes, cytoskeleton genes, and the entire fourth chromosome).


We are especially grateful to Ary Hoffman and Paul Schmidt for sharing fly stocks. We also thank Sergey Nuzhdin, Alisha Holloway, Chuck Langley, and many others for advice and support. Finally, we thank Umbreen Arshad, Liz Milano, Matt Rolston, and Heather Lindfors, and others in the Begun lab for technical and laboratory assistance. This work was supported by seed money from the Center for Population Biology at the University of California, Davis, by National Science Foundation predoctoral fellowships to T.L.T. and M.T.L., and by National Institutes of Health grant GM071926 to D.J.B.


  • Sequence data from this article have been deposited with the EMBL/GenBank Data Libraries under accession nos. EF670724EF672031 and EU476191EU476879.

  • Communicating editor: D. M. Rand

  • Received October 22, 2007.
  • Accepted March 8, 2008.


View Abstract