The effects of purifying selection on patterns of genetic differentiation between Drosophila melanogaster populations

Using the data provided by the Drosophila Population Genomics Project, we investigate factors that affect the genetic differentiation between Rwandan and French populations of D. melanogaster. By examining within-population polymorphisms, we show that sites in long introns (especially those >2000 bp) have significantly lower π (nucleotide diversity) and more low-frequency variants (as measured by Tajima's D, minor allele frequencies, and prevalence of variants that are private to one of the two populations) than short introns, suggesting a positive relationship between intron length and selective constraint. A similar analysis of protein-coding polymorphisms shows that 0-fold (degenerate) sites in more conserved genes are under stronger purifying selection than those in less conserved genes. There is limited evidence that selection on codon bias has an effect on differentiation (as measured by FST) at 4-fold (degenerate) sites, and 4-fold sites and sites in 8–30 bp of short introns ⩽65 bp have comparable FST values. Consistent with the expected effect of purifying selection, sites in long introns and 0-fold sites in conserved genes are less differentiated than those in short introns and less conserved genes, respectively. Genes in non-crossover regions (for example, the fourth chromosome) have very high FST values at both 0-fold and 4-fold degenerate sites, which is probably because of the large reduction in within-population diversity caused by tight linkage between many selected sites. Our analyses also reveal subtle statistical properties of FST, which arise when information from multiple single nucleotide polymorphisms is combined and can lead to the masking of important signals of selection.


INTRODUCTION
Natural populations are often divided into subpopulations. Studying the extent to which different subpopulations are genetically differentiated has been of paramount importance in evolutionary genetics, as it provides a way to examine how different evolutionary forces such as genetic drift, natural selection and migration drive changes in the genome (reviewed in Chapter 7 of Charlesworth and Charlesworth, 2010).
Specifically, insights into fundamental processes such as historical demographic changes, (local) adaptation and speciation can be obtained by comparing patterns of genetic differentiation across different genomic regions (Wu, 2001;Weir and Hill, 2002;Charlesworth et al., 2003;Hey and Machado, 2003;Beaumont, 2005;Holsinger and Weir, 2009). For instance, by scanning for loci that show unusually high levels of differentiation relative to the rest of the genome, we can detect loci that are under diversifying selection, whereby different alleles are favoured in different subpopulations (Beaumont and Nichols, 1996;Beaumont and Balding, 2004;Foll and Gaggiotti, 2008;Excoffier et al., 2009). As another example, in a study comparing African and non-African humans, it was found that the X chromosome was substantially more diverged than the autosomes, over and above the null expectation based on the fact that there are four copies of each autosome for every three copies of the X chromosome, which in turn suggests that dispersal in humans may be sex-biased or that the X chromosome may have experienced repeated selection after the divergence of African and non-African populations (Keinan et al., 2009).
Genetic differentiation between subpopulations is often measured by Wright's F ST (Wright, 1951), which is abbreviated as F in this study. F can be defined as the proportion of genetic variation explained by differences in allele frequencies between subpopulations (Charlesworth, 1998;Holsinger and Weir, 2009;Bhatia et al., 2013). F ranges between 0 and 1, which indicate no differentiation and fixed differences between subpopulations, respectively. Various genetic data, for example, single nucleotide polymorphisms (SNPs) and microsatellites, can be used to estimate F, but using statistical procedures that take into account biological properties of the data under consideration (for example, high versus low mutation rate) is vital for acquiring accurate estimates (Weir and Cockerham, 1984;Excoffier et al., 1992;Slatkin, 1995;Nagylaki, 1998;Holsinger and Weir, 2009).
Drosophila melanogaster, a classic model organism for population genetics, offers an invaluable system for studying population differentiation. Despite having a worldwide distribution in the present day, it is believed that the species originated in sub-Saharan Africa (David and Capy, 1988;Stephan and Li, 2007). The colonisation of Europe has been suggested to have taken place about 15 000 years ago (David and Capy, 1988;Stephan and Li, 2007;Duchen et al., 2013). The Americas and Australia were colonised much more recently, possibly in the past few hundred years (David and Capy, 1988;Stephan and Li, 2007;Duchen et al., 2013). By studying patterns of genetic differentiation, investigators have obtained evidence that American populations of D. melanogaster may be formed by admixture between African and European flies (Caracristi and Schlotterer, 2003). Multiple attempts have also been made to identify loci with unusually high F, which may have contributed to local adaptation to different habitats (Turner et al., 2008;Yukilevich et al., 2010;Kolaczkowski et al., 2011;Fabian et al., 2012;Langley et al., 2012;Pool et al., 2012;Campo et al., 2013).
These previous studies of D. melanogaster have mainly focused on determining the evolutionary relationship between subpopulations, quantifying the overall level of differentiation, and detecting genomic regions of interest using outlier scans. However, the role of purifying selection in shaping large-scale patterns of differentiation has not been well characterised, although it has been widely accepted that the majority of new mutations that affect fitness will have detrimental effects (Pal et al., 2006;Eyre-Walker and Keightley, 2007). Supporting this view, it has been estimated that only between 1 and 2% of new nonsynonymous mutations in D. melanogaster are (weakly) positively selected, and about 6% are nearly neutral (that is, |N e s|p1, where N e is the effective population size and s the selection coefficient), and the remaining are deleterious (|N e s|41) Eyre-Walker and Keightley, 2009;Schneider et al., 2011). Thus, we are interested in testing the following predictions based on population genetic theory of subdivided populations (reviewed in Chapter 7 of Charlesworth and Charlesworth, 2010): (i) purifying selection reduces differentiation between populations at functionally important regions; (ii) the level of reduction is positively correlated with the level of selective constraint. Answering these questions will help us better understand the sources of variation in genetic differentiation across the genome, which is important for example, in interpreting results obtained from genome scans (Beaumont and Nichols, 1996;Beaumont and Balding, 2004;Foll and Gaggiotti, 2008;Excoffier et al., 2009).
We will address the questions raised above by making use of the high-quality whole-genome resequencing data published by the Drosophila Population Genomics Project for one French population and one Rwandan population (Langley et al., 2012;Pool et al., 2012). In addition to protein-coding regions, we investigate introns, as previous studies have shown strong evidence that these genomic regions are under substantial selective constraints, probably as a result of the presence of cis-regulatory elements and noncoding RNA genes (Bergman and Kreitman, 2001;Parsch, 2003;Andolfatto, 2005;Haddrill et al., 2005a;Halligan and Keightley, 2006;Casillas et al., 2007;Roy et al., 2010).
Our study proceeds as follows. First, we present an overview of patterns of genetic variation both within and between populations using data from genomic regions where crossing over occurs (crossover (C) regions). We are interested in understanding whether patterns of differentiation at 4-fold degenerate (hereafter 4-fold) sites are affected by selection on codon usage, and whether 4-fold sites and putatively neutral sites in 8-30 bp regions of introns p65 bp (Halligan and Keightley, 2006;Parsch et al., 2010) are comparable with respect to levels of differentiation. These are intended to identify putatively neutral sites which can be used as a reference in the study of the effects of purifying selection on genetic differentiation. We then examine the relationship between K A (nonsynonymous divergence) and diversity/differentiation patterns at 0-fold degenerate (hereafter 0-fold) sites in protein-coding regions, as well as the relationship between intron length and diversity/differentiation patterns in intronic regions. Finally, we compare non-crossover (NC) regions (for example, the fourth chromosome) and C regions regarding differentiation patterns, study the relative contribution of selection and genetic linkage, and examine the correlation between local recombination rates and F at putatively neutral sites.
Sites with a quality score below 31 (equivalent to a PHRED score of 48, and approximately equivalent to one error per 100 kb; see Pool et al., 2012) were also masked.
From the FASTQ files, we extracted protein-coding regions in C regions, which we abbreviate as CDS-C, using gene annotations from FlyBase release version 5.44 (www.flybase.org) and made FASTA files containing all samples (24 alleles). For each D. melanogaster gene with multiple transcripts, we chose one transcript randomly.
Introns belonging to our chosen transcript were extracted, and were further processed by masking any coding regions that form part of the other transcripts. Only introns occurring in C regions were retained (polymorphism data for introns in NC regions were of significantly lower quality, and were therefore excluded).
Protein-coding regions in NC regions of the D. melanogaster genome, abbreviated as CDS-NC, were obtained from Campos et al. (2012). These data included five unlinked regions: N2 (genes located in heterochromatic regions near the centromere of the second chromosome), N3 (genes located in heterochromatic regions near the centromere of the third chromosome), N4 (the fourth chromosome), NXc (X-linked genes near the centromere) and NXt (X-linked genes near the telomere).
For all CDS-C and CDS-NC, we selected D. yakuba as an out-group to avoid any major influence of ancestral polymorphisms on the estimation of sequence divergence, which can potentially create spurious correlations between sequence divergence and recombination (for example, Cutter and Choi, 2010). One-to-one orthologous D. yakuba sequences were obtained from FlyBase (available at ftp://ftp.flybase.net/releases/FB2012_02/precomputed_files/genes/ gene_orthologs_fb_2012_02.tsv.gz). We then performed amino-acid sequence alignments using MAFFT (Katoh et al., 2002). These amino-acid sequence alignments were translated back to nucleotides using custom scripts in PERL to produce in-frame coding sequence alignments that included the 24 D. melanogaster alleles and the D. yakuba reference.
For introns, we used D. simulans as an out-group because we considered D. yakuba too distant for producing reliable alignments, because of the increased prevalence of indels in non-coding regions. We obtained orthologous intronic D. simulans sequences from Hu et al. (2013), which was based on an updated D. simulans genome assembly and careful alignment procedures to preserve gene structures (http://genomics.princeton.edu/AndolfattoLab/ w501_genome_files/alnMSY.tar.gz).
Recombination rate for the midpoint of all alignments was obtained using the Drosophila melanogaster Recombination Rate Calculator v2.3 (Fiston-Lavier et al., 2010) and the high-resolution genetic map published recently by Comeron et al. (2012).

Sequence analysis
For CDS-C and CDS-NC, we calculated K A and K S (the numbers of nonsynonymous and synonymous substitutions per nonsynonymous and synonymous site, respectively) using the kaks() function from the seqinr package (Charif and Lobry, 2007) in R (http://www.r-project.org/), which implements the method of Li (1993). For introns, we calculated divergence (K) to the D. simulans reference using the dist.dna() function in the ape package of R (Paradis et al., 2004), with the 'K80' method (Kimura, 1980). For conducting analyses using polymorphism data in the two D. melanogaster samples, we split CDS-C and CDS-NC into 0-fold degenerate sites and 4-fold degenerate sites by analysing the alignments codon by codon. A codon column was retained if the following requirements were met: (i) data from all individuals were available; (ii) it had at most one SNP. These were to avoid uncertainty of the order of mutations in codons with multiple SNPs. We retained 7235 autosomal and 1150 X-linked CDS alignments for which we had both more than 10 bp of 0-fold sites and more than 10 bp of 4-fold sites.
We split introns into short (p65 bp long) and long (465 bp long) classes, and further trimmed short introns to retain positions 8-30 from the 5 0 end (hereafter the SI sites), to retain sites under the least amount of selective constraint (Halligan and Keightley, 2006;Parsch et al., 2010). This left us with 7483 autosomal and 752 X-linked short introns, and 8851 autosomal and 1869 X-linked long introns. To keep the sample size the same as the CDS data, only intronic sites where data from all individuals were available were retained. This requirement appears to be conservative with respect to detecting SNPs (note also that regions within 3 bp of an indel were also masked by Drosophila Population Genomics Project; see Pool et al., 2012 and http://www.dpgp.org/ dpgp2/DPGP2.html). For instance, p (nucleotide diversity) estimated using data from SI sites that fulfilled the above criterion was 0.0145 in the RG sample, whereas the estimate increased to 0.0164 when we retained, within the same genomic regions, all sites that had data from at least two individuals. However, as we will show in the Results, there is no detectable difference between the SI and 4-fold sites in terms of skewness of allele frequency spectrum (as measured by Tajima's D), average minor allele frequency (MAF), and F ST . Our conservative data filtering procedure is unlikely to bias our analysis of population differentiation. Purifying selection and population differentiation BC Jackson et al Because of the complete linkage between the CDS-NC genes located within a NC region, each non-recombining region of the genome effectively represents a single locus. Therefore all genes within a single NC region were concatenated and analysed as a single gene. In total, we have three autosomal NC regions, which are N2, N3 and N4, and two X-linked NC regions, which are NXc and NXt. These loci were kept intact in the permutation tests used to compare values of summary statistics calculated using data from NC and C regions.
Nucleotide diversity (p), Tajima's D (Tajima, 1989) and relative Tajima's D, (Schaeffer, 2002) were calculated using nuc.div() and a modified version of the tajima.test() function, both from the pegas package in R (Paradis, 2010). Because conclusions drawn from Tajima's D and relative Tajima's D are identical, only the former is presented. Permutation tests were carried out to assess whether these statistics were different between different types of sites. For instance, to investigate whether values of Tajima's D at 4-fold sites and SI sites were comparable, 10 000 pseudosamples were generated by randomly shuffling both the 4-fold and SI sites in the data (SNPs from a single locus were shuffled as a unit), such that in each pseudosample there were similar numbers of SNPs in the two 'site classes' as in the real data.
To assess the effects of selection on codon bias on population differentiation at 4-fold sites, we calculated frequency of optimal codons (Fop) with CodonW (Peden, 1999), using the built-in table of optimal codons for D. melanogaster.

Measuring population differentiation
Levels of differentiation between the two populations of D. melanogaster were measured by Wright's F ST , which is abbreviated as F. We used the definition of Weir and Cockerham (1984), which can be expressed as where p B is the expected divergence between a pair of alleles sampled from two different populations, and p S is the expected within-population diversity (see also Charlesworth, 1998;Keinan et al., 2007). Previous investigations (Maruki et al., 2012;Jakobsson et al., 2013) of the dependence of F on MAF were based on a different definition of F put forward by Nei (1973). We therefore derive the maximum value of Weir and Cockerham's definition of F as a function of MAF. Consider a population divided into two subpopulations. We examine a biallelic locus. The frequency of one of the two alleles in the k-th subpopulation is referred to as p k (k ¼ 1 or 2). It can be shown that (Charlesworth, 1998). Substituting these into Equation (1), F can be rewritten as where d ¼ |p 1 Àp 2 | and s ¼ p 1 þ p 2 . Without loss of generality, we assume that 1psp2. Three properties are of use: (i) MAF ¼ 1-s/2; (ii) 0p2s-s 2 p1 for 1psp2; (iii) 0pdp(2-s) 2 . Rearranging Equation (2), we deduce that Because MAF ¼ 1-s/2, the above inequality is equivalent to max(F)p2 MAF. In Supplementary Figure S1, we display the differences between the upper bounds of F derived here and that obtained in previous studies using Nei's F (Maruki et al., 2012;Jakobsson et al., 2013). It can be seen that F can only assume a very restrictive range of values when MAF is small.
The estimator of F proposed by Hudson et al. (1992) (see also Keinan et al., 2007;Bhatia et al., 2013) was employed: wherep B andp S are estimates of p B and p S obtained from data. Equation (4) can be calculated using information from a single SNP. To combine information from multiple SNPs, the following two methods were used (Weir and Cockerham, 1984;Bhatia et al., 2013): and where S is the number of SNPs, andF ðiÞ ;p ðiÞ B andp ðiÞ S are values of the terms defined in Equation (4) obtained using data from the i-th SNP.
It should be noted that F U gives equal weight to all SNPs, whereas F W gives more weight to SNPs with higher expected levels of polymorphism. In other words, F U is expected to be more sensitive to the presence of SNPs with low MAFs, but F W is dominated by SNPs that are on average more polymorphic. To see this more explicitly, assume that we are combining information from two SNPs (that is, S ¼ 2). We add a subscript j to the symbols defined above to signify the locus under consideration, so that we have p jk , s j and d j ¼ |p j1 Àp j2 |. We further assume that 24s 1 Xs 2 X1. Note that s j are regarded as parameters (for example, a SNP under strong selective constraints is expected to have a larger s (that is, a smaller MAF) than a neutral SNP). Some straightforward algebra leads to the following results: Às j (these results hold when S42; proof not shown). To see the differential sensitivities to SNPs with small MAFs, we define Thus, the behaviour of F W is more akin to that of the more polymorphic SNP (that is, max[F W ] is closer to max[F(s 2 )] than to max[F(s 1 )]). As we will see later, this property of F W can lead to the masking of important signatures of evolution when SNPs with different properties are combined. Table 1 presents summaries of polymorphism patterns for autosomal (A) and X-linked (X) loci situated in genomic regions where crossingover occurs (C regions). For ease of presentation, we will refer to nucleotide diversity, p, calculated using 0-fold sites, 4-fold sites and SI sites (positions 8-30 from the 5 0 end of short introns p65 bp) as p 0 , p 4 and p SI , respectively; a similar notational convention will be used for other statistics. For both A and X, and in both the Rwandan (RG)  (5) and (6) and French (FR) samples, p 0 , Tajima's D 0 (Tajima 1989), and MAF 0 are significantly smaller than the corresponding estimates obtained from 4-fold and SI sites (P permutation o0.001 in all cases), consistent with the well-known fact that most nonsynonymous mutations are deleterious (Pal et al., 2006;Eyre-Walker and Keightley, 2007), and are therefore kept at low frequencies in the population by purifying selection (Kimura, 1983). Previous studies have suggested that SI sites may be neutrally evolving (Halligan and Keightley, 2006;Parsch et al., 2010). In our data set, p SI seems to be somewhat smaller than p 4 , which may be due to the stringent data filtering procedure we employed (see Materials and Methods), or the higher GC content at 4-fold sites compared to intronic sites, which in turn is expected to result in an increased mutation rate in 4-fold sites (Singh et al., 2005;Keightley et al., 2009). There is, however, no statistically discernible difference with respect to either MAF or Tajima's D between 4-fold and SI sites (Table 1; P permutation 40.1 for both A and X).

Genome-wide polymorphism patterns in crossover regions
The FR sample has a lower level of diversity than RG for all three types of sites (Table 1), reflecting a loss of genetic variation induced by population bottlenecks which are believed to have occurred as the species migrated out of Africa (Haddrill et al., 2005b;Li and Stephan, 2006;Thornton and Andolfatto, 2006;Hutter et al., 2007;Duchen et al., 2013). The difference in p 0 between the two populations is somewhat smaller than those observed for p 4 and p SI (for example, on A, p 0 (FR)/p 0 (RG) ¼ 0.83 versus p 4 (FR)/p 4 (RG) ¼ 0.77). This is probably because more 0-fold sites are under strong selective constraint, so that variants at these sites behave almost deterministically, and are therefore less sensitive to demographic changes (for example, Zeng, 2013).
To inspect overall patterns of genetic differentiation between the RG and FR populations, we calculated F ST (abbreviated here as F; see Equation (1) in Materials and Methods), as defined by Weir and Cockerham (1984), using the estimator of Hudson et al. (1992). Two approaches were employed to combine information over multiple SNPs: un-weighted mean F (Equation (5)) and weighted mean F (Equation (6)), which will be referred to as F U and F W , respectively. Because most nonsynonymous mutations are likely to be deleterious, it is expected that levels of population differentiation at these selectively constrained sites should be lower than those at less constrained sites (for example, 4-fold sites) (Barreiro et al., 2008;Maruki et al., 2012). Surprisingly, values of F W 0 , estimated using either the autosomal or X-linked data, are not statistically different from those of either F W 4 or F W SI (Table 1; P permutation 40.1 in all cases). There is also no detectable difference between F W 4 and F W SI (P permutation 40.1 for both A and X). In contrast, F U 0 was found to be significantly smaller than both F U 4 and F U SI (P permutation o0.001 for both A and X), whereas the differences between F U 4 and F U SI remain non-significant (P permutation 40.1 for both A and X). The patterns obtained from F U are therefore more compatible with the a priori expectation that 0-fold sites are on average more constrained than 4-fold and SI sites. We will investigate causes for the lack of difference between F W 0 and either F W 4 or F W SI in a later section. Several differences between A and X are of note (Table 1). Firstly, consistent with previous reports (Caracristi and Schlotterer, 2003;Hutter et al., 2007;Charlesworth, 2012b;Pool et al., 2012;Campos et al., 2013), the X:A ratio in diversity at putatively neutral sites (that is, 4-fold and SI sites) is about 1 in the RG population (p 4 (X)/ p 4 (A) ¼ 1.08 and p SI (X)/p SI (A) ¼ 1.10), higher than the null expectation of 3/4. Secondly, the reduction in diversity in FR is more pronounced for X than A for all three types of sites (for example, p 4 (FR)/p 4 (RG) ¼ 0.41 and 0.77 for X and A, respectively), as reported in previous investigations (Caracristi and Schlotterer, 2003;Hutter et al., 2007). Finally, the extent of population differentiation at both 4-fold and SI sites, as measured by either F U or F W , is significantly higher on the X than on A (P permutation o0.001 for all comparisons). This is probably largely driven by the greater reduction in diversity on the X in non-African populations, as values of D xy , the mean number of nucleotide substitutions between sequences taken from different subpopulations (Nei and Miller, 1990), are comparable between A and X in this study: D xy,4 ¼ 1.65 and 1.64%, and D xy,SI ¼ 1.51 and 1.58%. A systematic examination of possible causes of the apparent differences between A and X is beyond the scope of this study; the interested reader can refer to previous studies of this topic (Charlesworth, 2001;Pool and Nielsen, 2007;Singh et al., 2007;Pool and Nielsen, 2008;Yukilevich et al., 2010;Charlesworth, 2012b;Campos et al., 2013). In what follows, results obtained from A and X will be presented separately.
Limited evidence for selection on codon usage bias affecting patterns of population differentiation at 4-fold degenerate sites To investigate whether selection on codon usage bias affects differentiation patterns at 4-fold sites, we first examined the relationship between F U 4 and Fop, as the latter is well known to be correlated with the intensity of selection on codon usage bias (reviewed in Hershberg and Petrov, 2008;Zeng and Charlesworth, 2009). Considering the large variance of the F estimators and the dearth of SNPs in individual genes, we grouped the genes into equal-sized bins with similar numbers of SNPs at 4-fold sites. As shown in Supplementary Figure S2A, Fop and F U 4 are not correlated on A (Kendall's t ¼ À0.01, P40.1). On the X, some evidence for a weak negative correlation was obtained (Supplementary Figure S2B), but it is not statistically significant (Kendall's t ¼ À0.6, P ¼ 0.13). When F W 4 was considered, no correlation was found on either A or X (Supplementary Figures  S2E and F). To investigate this further, for the genes within each bin on the X, we tested whether F U 4 differed from F U SI statistically. Amongst the six bins, no evidence of a significant difference was found for the first four bins, whereas the differences were marginally significant for the last two bins with highest Fop (P permutation ¼ 0.04 and 0.05, respectively). Similarly, we did not detect any correlation between K S and either F U 4 or F W 4 (Supplementary Figure S2). Overall, there is limited evidence that selection on codon usage bias is strong enough to substantially alter patterns of genetic differentiation at 4-fold sites. Considering that 4-fold and SI sites in C regions are comparable with respect to both MAF and F, in what follows, we will use population differentiation patterns obtained from the two types of site as neutral standards, and will refer to them as putatively neutral sites.
Evolutionarily conserved genes are under stronger purifying selection and have reduced F at 0-fold degenerate sites Genes in C regions were divided into equal-sized bins (with similar numbers of SNPs) based on their K A values between D. melanogaster and D. yakuba. We inspected polymorphism patterns in the RG sample as a function of K A ; a qualitatively identical set of results were obtained using the FR sample (Supplementary Figure S3). On both A and X, K A was found to be significantly positively correlated with both p 0 (Figures 1a and b; A: Kendall's t ¼ 0.989 and Po0.001; X: Kendall's t ¼ 1 and P ¼ 0.009) and Tajima's D 0 (Figures 1c and d; A: Kendall's t ¼ 0.884, Po0.001; X: Kendall's t ¼ 0.867 and P ¼ 0.024). No statistically significant relationship was found when comparing K A with Tajima's D 4 (Figures 1c and d; Kendall's t ¼ À0.2 and À0.333, P40.1, for X and A), although there is a negative correlation between K A and p 4 on A (Figure 1a; Kendall's t ¼ À0.6, Po0.001) (see also  Andolfatto, 2007;Haddrill et al., 2011). In particular, on both A and X, p 0 and Tajima's D 0 approach p 4 and Tajima's D 4 , respectively, as K A increases. In contrast, values of p 4 and Tajima's D 4 , regardless of the K A bin from which they were obtained, remain similar to the values of p SI and Tajima's D SI . These results suggest that 0-fold sites are under stronger constraints than 4-fold and SI sites, and that 0-fold sites in genes with smaller K A are, on average, under stronger purifying selection. We obtained the same results when we used the D. simulans genome as an out-group (Supplementary Figure S4). Figures 2a and b show that evolutionarily conserved genes have significantly smaller F U 0 (A: Kendall's t ¼ 0.663, Po0.001; X: Kendall's t ¼ 0.867, P ¼ 0.02). Again, we obtained the same result when using D. simulans as the out-group (Supplementary Figure S5). The pattern remains statistically significant for autosomes when F W 0 was considered (Supplementary Figure S6). The reduction in F 0 for genes with smaller K A is associated with a strong reduction in MAF 0 (Figures 2c and d) and an increase in the proportion of 0-fold SNPs that are private to one of the two populations (Figures 2e and f), both of which are hallmarks of selection against deleterious mutations (cf., recent findings in humans; Nelson et al., 2012;Fu et al., 2013), and are expected to drive both F U and F W downwards, as shown in Materials and Methods (see also Maruki et al., 2012;Bhatia et al., 2013;Jakobsson et al., 2013). For the 4-fold sites on both A and X, no correlation with K A was observed for F U , F W , MAF and the proportion of private SNPs (Figure 2; P40.1 in all cases based on Kendall's t).
The data presented in Figures 1 and 2 suggest that the lack of difference between F W 0 and either F W 4 or F W SI reported in the previous section is probably because of the fact that F W gives more weight to SNPs with higher expected levels of polymorphism (for example, nearly neutral variants), as we have shown in Materials and Methods. In other words, when all 0-fold sites in C regions were analysed together (Table 1), the effects of purifying selection on a substantial fraction of 0-fold sites were probably masked by those 0-fold sites that are nearly neutrally evolving. Consequently, the overall distribution of F W 0 appears non-distinguishable from those of F W 4 and F W SI . In contrast, F U gives equal weight to all SNPs. Considering that the value of F when calculated using a single SNP is constrained by MAF (see Equation (3) in Materials and Methods), F U is expected to be more sensitive to the action of purifying selection than F W , consistent with the observation reported above. In the Discussion, we will further explore the implications of these statistical properties of F, which arise when information from multiple SNPs is combined.

Longer introns are under stronger selective constraints and are less differentiated
In agreement with earlier findings (Haddrill et al., 2005a;Halligan and Keightley, 2006), longer introns tend to have lower divergence (K) between D. melanogaster and D. simulans (A: Kendall's t ¼ À0.635, Po0.001; X: Kendall's t ¼ À0.486, Po0.001; Figures 3a and b), probably as a result of the presence of functional elements that are subject to purifying selection (Bergman and Kreitman, 2001;Parsch, 2003;Andolfatto, 2005;Haddrill et al., 2005a;Halligan and Keightley, 2006;Casillas et al., 2007;Roy et al., 2010). Here, we report further support for this hypothesis by examining within-population polymorphism patterns as a function of intron length. Consistent with the action of purifying selection, longer introns have lower p (Figures 3c  and d) and more negative Tajima's D (Figures 3e and f) compared with 4-fold and SI sites (similar results were observed in the FR sample; see Supplementary Figure S7). Interestingly, the patterns of divergence and polymorphism level off for introns longer than 2000 bp. Using the RG sample, the values of p and Tajima's D obtained from introns longer than 2000 bp are 0.0072 and À0.5476 for A, and 0.0076 and À0.9013 for X, respectively; all these values are substantially lower than the corresponding values observed at 4-fold and SI sites, but are higher than those obtained from 0-fold sites (see Table 1). Furthermore, the K A values for CDS in C regions between D. melanogaster and D. simulans are 0.015 and 0.018 for A and X, respectively, which are significantly smaller than the values of K for long introns 42000 bp on A and X, which are 0.061 and 0.074, respectively (Mann-Whitney U test, Po0.001). These results imply that long introns, especially those 42000 bp, are more constrained than the 4-fold and SI sites, but probably contain fewer strongly selected sites than 0-fold sites.
Estimates of F W , when calculated using sites from introns more than 65 bp in length, were 0.171 and 0.283 for A and for X, respectively. None of these was found to be statistically different from the corresponding values estimating using 4-fold and SI sites reported in Table 1 (P permutation 40.1 in all cases). F U for introns 465 bp were 0.157 and 0.174 for A and X, respectively, both of which were significantly smaller than both F U SI and F U 4 (P permutation o0.001 in all cases). There is a clear negative relationship between F U and intron length (Figures 4a and b; for A and X, Kendall's t ¼ À0.356 and À0.364; P ¼ 0.010 and Po0.001, respectively), which mirrors that between MAF (or the prevalence of private SNPs) and intron length (Supplementary Figure S8), and is consistent with the expected effect of purifying selection on genetic differentiation between populations. The relationship between differentiation and intron length is weaker when F W was analysed (Supplementary Figure S8; for A and X, Kendall's t ¼ À0.271 and À0.146, and P ¼ 0.05 and 0.16, respectively). These differences between F W and F U can be explained by the fact that fewer sites in introns 465 bp are expected to be strongly selected compared with 0-fold sites. As discussed in the previous section, F W , which tends to reflect differentiation patterns at neutral sites in the data, is less likely to recover signatures of purifying selection compared to F U .
Regions with reduced recombination tend to have higher F It is known that genomic regions that lack crossing over (NC regions) have very different patterns of divergence and polymorphism than those seen in C regions (Haddrill et al., 2007;Betancourt et al., 2009;Arguello et al., 2010;Campos et al., 2012;Campos et al., 2014). In Table 2, we present summary statistics of the NC data pertinent to the current study (see Materials and Methods for a list of the NC regions considered). It can be seen that, for both 0-fold and 4-fold sites, values of F in NC regions are generally higher than those obtained using the same type of site in C regions, regardless of the way in which information from multiple SNPs was combined. Specifically, the average K A to D. yakuba is about 0.05 for the NC loci (Campos et al., 2012). F U 0 calculated using autosomal and X-linked NC data are 0.1817 and 0.3012, respectively (Table 2), higher than the values of 0.1569 and 0.1685 for autosomal and X-linked genes in C regions spanning the same K A values (Figures 2a and b; P permutation ¼ 0.05 for A and P permutation o0.001 for X).
It should be noted that the elevation in F in NC regions is probably caused by an extreme reduction in within-population diversity induced by tight linkage between a large number of selected sites (Table 2; Kaiser and Charlesworth, 2009;O'Fallon et al., 2010;Seger et al., 2010;Zeng and Charlesworth, 2010). This is because F is a relative measure of differentiation (see Equation (1)), and therefore all else being equal, F is expected to be elevated by forces that reduce within-population diversity (that is, p S in Equation (1)), irrespective of whether diversifying selection or reduced gene flow has affected the genomic region under study (Charlesworth, 1998;Noor and Bennett, 2009).
To further examine the effects of selection at linked sites, we inspect the correlation between F at putatively neutral sites and local recombination rates in C regions. Figure 5 presents results based on autosomal loci, where it can be seen that F U 4 is reduced with more frequent recombination (Kendall's t ¼ À0.474, P ¼ 0.004; the data point obtained from the NC regions was not included in the calculation). However, there is no statistically significant relationship between recombination rate and F U SI (Figure 5b; Kendall's t ¼ À0.179 and P ¼ 0.28). Weak negative correlations were also found on the X chromosome for 4-fold and SI sites (Supplementary Figure S9). The patterns remained unchanged when F W was used (Supplementary Figure S10).

DISCUSSION
By using the high-quality data provided by the Drosophila Population Genomics Project, we have found that evolutionary conserved regions (that is, genes with lower K A and longer introns) show clear evidence of more intense on-going purifying selection than less conserved genomic regions, which can be detected by analysing patterns of genetic variation both within and between subpopulations. The negative correlation between p and intron length reported in Figure 3 extends the study by Parsch et al. (2010) who examined a much smaller data set and did not find evidence of such a correlation. Because we did not find support for a correlation between local recombination rate and intron length (Kendall's t ¼ À0.004 and 0.011 for A and X, respectively, and P40.1 in both cases; Supplementary Figure S11) (cf., Carvalho and Clark, 1999;Comeron and Kreitman, 2000), the relationship is unlikely to be driven by the well-known positive correlation between diversity and recombination. It is unclear why the effect of intron length levels off for introns longer than 2000 bp. Analysis of theoretical models (for example, Ometto et al. 2005) and improved annotation of noncoding functional elements (for example, Roy et al. 2010) are both needed to solve this problem. Finally, there is evidence that the severe reduction in within-population diversity in NC regions of the genome induced by tight linkage between selected sites has led to elevated F ST values, but there is limited support for this effect in C regions.
Purifying selection as a major determinant of population differentiation Our analysis reveals (i) a positive correlation between K A and F 0 ( Figure 2) and (ii) a negative correlation between intron length and F calculated using intronic sites (Figure 4). After examining other aspects of polymorphism and differentiation patterns (Figures 1  and 3), we suggest that the observations can be most readily explained by differential intensity of purifying selection acting on different parts of the genome. Similar observations have also been reported in humans (Barreiro et al., 2008;Maruki et al., 2012), suggesting the universal importance of purifying selection as a factor that shapes genetic differentiation between populations.
It should be noted that the above conclusion is not inconsistent with the existence of outlier loci with unusually high F, which may have been caused by diversifying selection (Turner et al., 2008;Yukilevich et al., 2010;Kolaczkowski et al., 2011;Fabian et al., 2012;Langley et al., 2012;Pool et al., 2012;Campo et al., 2013). Our analysis intends to detect forces with large-scale effects (there are typically hundreds of genes in each of the bins in our analysis), and is therefore unlikely to respond to processes that have more localised effects in the genome. In fact, it has been suggested that the number of loci contributing to differences between populations may be relatively small (   example, after taking into account the confounding effects of complex demography and correcting for multiple testing, only four loci had strong statistical support for being driven to high levels of differentiation by diversifying selection between North American and African populations of D. melanogaster (Yukilevich et al., 2010). Furthermore, in line with the low level of linkage disequilibrium in the D. melanogaster genome (for example, Pool et al., 2012), previous genome scan studies have shown that most candidate variants that show evidence of involving in local adaptation only affect differentiation patterns in its immediate neighbourhood, typically on the order of the size of a gene (Kolaczkowski et al., 2011;Fabian et al., 2012). Finally, we have focussed on protein-coding regions and introns, whereas a substantial number of previously found candidate loci fall within intergenic regions. A noticeable exception is chromosome 3R, in which the cosmopolitan inversion In(3R)P is situated. Multiple studies concerning differentiation between various D. melanogaster populations have found that chromosome 3R has a disproportionally large number of candidate loci, especially within the In(3R)P region, and that these candidate variants tend to affect differentiation patterns in a larger genomic neighbourhood (Kolaczkowski et al., 2011;Fabian et al., 2012). To further test the robustness of our results, we repeated the analysis leading to Figure 2a by removing all genes on chromosome 3R, and found that the pattern remains unchanged (Supplementary Figure S12). In summary, it is unlikely that highly differentiated regions driven by adaptive changes have made a substantial contribution to our observations. The relationship between F and recombination As pointed out previously (Charlesworth, 1998;Noor and Bennett, 2009), forces that reduce within-population diversity can lead to elevated F ST values in the absence of diversifying selection and restricted gene flow. Hence, in light of the lack of evidence of adaptive evolution in NC regions of the D. melanogaster genome (Betancourt et al., 2009;Arguello et al., 2010;Campos et al., 2014), the high F values obtained from NC regions is probably a result of the diversity-reducing effect of linkage between selected sites, which is often referred to Hill-Robertson interference or HRI (Hill and Robertson, 1966;Comeron et al., 2008;Sella et al., 2009;Charlesworth, 2012a;Cutter and Payseur, 2013). Within the C regions, although negative correlations between F at putatively neutral sites and local recombination rate, as predicted by the HRI theory, were observed ( Figure 5, Supplementary Figures S9 and S10), these patterns are weak and often non-significant. Langley et al. (2012) also reported weak negative correlations between a different measure of genetic differentiation and fine-scale recombination rates estimated from linkage disequilibrium patterns, but the relationship was inconsistent between chromosome arms and was sometimes weakly positive when broad-scale recombination rates were used.
The weak association between F and recombination in C regions is somewhat surprising given that both p 4 and p SI are clearly positively correlated with local recombination rates in both the RG and FR populations (Supplementary Figure S13). A possible explanation is that, because hitchhiking effects induced by both positive and negative selection can lead to an excess of low-frequency variants at linked neutral sites (Charlesworth et al., 1993;Braverman et al., 1995;Zeng and Charlesworth, 2011), the negative correlation between F and recombination may be weakened, if rare variants are more common in low-recombination regions, as these variants tend to lower F (see Equation (3)). Tajima's D is somewhat more negative in autosomal C regions with reduced recombination (Supplementary Figure S14), but it is hard to determine to what extent this has contributed to the observations in Figure 5 and Supplementary Figure S9, especially when noting that NC regions have more negative Tajima's D and yet higher F ST values. Further research that takes into account HRI, demography and statistical properties of estimators of F (see below) is needed to clarify the matter.
The importance of sampling strategy regarding using F ST to study population differentiation As is the case for other definitions of F ST (Maruki et al., 2012;Jakobsson et al., 2013), Weir and Cockerham's F ST can only take a very restricted range of values when MAF is small (max(F ST )p2 MAF; Supplementary Figure S1). When information is combined across SNPs, the weighted mean F ST (F W ) is likely to be dominated by SNPs that are more polymorphic (that is, those having a higher expected MAF). This can lead to the masking of signals of purifying selection, as we have shown above. Thus, F W may be a better choice when the intention is to ascertain the overall level of genetic differentiation. In this case, as long as the data contain a substantial number of putatively neutrally evolving variants, a reasonably accurate estimate can be obtained, even in the presence of sites under strong selective constraints. In contrast, the unweighted mean F ST (F U ) gives equal weight to all SNPs, and is more responsive to the presence of rare variants (for example, those under purifying selection). These considerations, as well as the recommendations proposed by Bhatia et al. (2013), suggest that care should be exercised when deciding which sampling strategy is most appropriate for the question in hand.

DATA ARCHIVING
All genetic data analysed in this study are publically available.