Introgression in the Drosophila subobscura--D. Madeirensis sister species: evidence of gene flow in nuclear genes despite mitochondrial differentiation.

Species hybridization, and thus the potential for gene flow, was once viewed as reproductive mistake. However, recent analysis based on large datasets and newly developed models suggest that gene exchange is not as rare as originally suspected. To investigate the history and speciation of the closely related species Drosophila subobscura, D. madeirensis, and D. guanche, we obtained polymorphism and divergence data for 26 regions throughout the genome, including the Y chromosome and mitochondrial DNA. We found that the D. subobscura X/autosome ratio of silent nucleotide diversity is significantly smaller than the 0.75 expected under neutrality. This pattern, if held genomewide, may reflect a faster accumulation of beneficial mutations on the X chromosome than on autosomes. We also detected evidence of gene flow in autosomal regions, while sex chromosomes remain distinct. This is consistent with the large X effect on hybrid male sterility seen in this system and the presence of two X chromosome inversions fixed between species. Overall, our data conform to chromosomal speciation models in which rearrangements are proposed to serve as gene flow barriers. Contrary to other observations in Drosophila, the mitochondrial genome appears resilient to gene flow in the presence of nuclear exchange.

In early speciation models, a severe reduction of gene flow initially facilitated by geographic isolation (i.e., allopatry) was required to maintain distinctness between coexisting clusters (Dobzhansky 1937; Mayr 1942, 1963; Coyne 1992; Coyne and Orr 1998, 2004). However, the observation that approximately 25% of plant and 10% of animal species hybridize in nature suggests that gene flow is not only possible, but likely to occur between closely related species (Mallet 2005). Studies using large datasets, together with newly developed statistical tools, have allowed us to recognize more complex models of speciation wherein certain genomic regions are extensively divergent, whereas others experience the homogenizing effects of gene-sharing (Wu 2001; Pinho and Hey 2010; Yang 2010). These scenarios portrait species’ genomes as mosaics and present a new line of inquiry questioning the amount of gene exchange allowed between two species while preserving their evolutionary independence (Hey and Pinho 2012).

The study of gene flow between species is relevant for at least two additional reasons beyond shedding light on evolutionary independence. The first is the identification of genes responsible for reproductive isolation between species (i.e., “speciation” genes). Situations where secondary contact and hybridization follows a period of geographic isolation represent “natural laboratories,” where the effects of genes and genic interactions on hybrid fitness can be assessed in their natural environment (Barton and Hewitt 1981b; Hewitt 1988; Harrison 1990). Genes responsible for species-specific traits or reproductive isolation will not be able to cross species boundaries, whereas genes with no deleterious effects on hybrid fitness will move more freely across hybrid zones and species’ genetic pools (Barton and Hewitt 1981a; Hewitt 1988; Harrison 1990; Navarro and Barton 2003a; Hey 2006; Teeter et al. 2008). By studying gene flow in hybrid zones, we gain insight into the genomic location of factors likely responsible for species identity and reproductive isolation (Teeter et al. 2008; Garrigan et al. 2012). The second reason to study gene flow between species is that it can serve as a source of adaptive genetic variation (Arnold 1997, 2006; Arnold and Martin 2009; Rieseberg 2009). Species can access genetic variation that originated in another species through hybridization, with the possibility of benefitting from nonnative adaptive mutations (Castric et al. 2008; Kim et al. 2008; Rieseberg 2011; Song et al. 2011; Pardo-Diaz et al. 2012; Staubach et al. 2012). Traits that depend on multiple substitutions at the same or different genes are particularly more attainable through introgression of preexisting variants than through de novo mutations (Rieseberg 2009). This latter scenario can be particularly relevant for species with small effective population sizes that have limited levels of standing variation available and reduced effectiveness of natural selection.

To date, gene flow between different species of Drosophila has been investigated in the obscura and melanogaster groups of the subgenus Sophophora and in the Drosophila repleta group: between D. pseudoobscura and D. persimilis (Machado et al. 2002), D. simulans and D. mauritiana (Ballard 2000; Kliman et al. 2000; Nunes et al. 2010; Garrigan et al. 2012), D. simulans and D. sechellia (Kliman et al. 2000; Garrigan et al. 2012), D. yakuba and D. santomea (Llopart et al. 2005; Bachtrog et al. 2006), and D. mojavensis and D. arizonae (Counterman and Noor 2006). Although mitochondrial gene exchange accompanies variable degrees of nuclear gene flow in all the Sophophora cases above, additional studies are needed to further examine this correlation. We use a multilocus approach to determine whether the sister species Drosophila madeirensis and D. subobscura show evidence of gene flow in nuclear and mitochondrial DNA (mtDNA). Both species belong to the obscura group and are believed to have diverged allopatrically between 0.6 and 1.0 million years ago (Mya; Ramos-Onsins et al. 1998). Drosophila madeirensis is endemic to Madeira, a small volcanic island located approximately 580 km west of Morocco in the Atlantic Ocean (Monclus 1984). This species is invariantly associated with the Laurisilva forest, a relic of the subtropical forests of the Tertiary era that once covered the entire Mediterranean region. Drosophila subobscura can be found throughout the Palearctic region and has recently also invaded the Americas (Ayala et al. 1989; Rozas et al. 1990; Mestres et al. 2001). Due to secondary invasion, D. subobscura currently inhabits Madeira Island where its range overlaps with that of D. madeirensis. Because F1 hybrid females are fertile in laboratory conditions (Khadem and Krimbas 1991; Papaceit et al. 1991), their overlapping geographic ranges open the possibility of natural hybridization and gene flow.

In a recent study, Khadem et al. (2012) investigated nucleotide variation patterns at six nuclear regions and found no evidence of gene flow between D. madeirensis and D. subobscura. However, all investigated regions were located very close to inversion breakpoints or along a chromosomal inversion that is currently fixed between D. madeirensis and D. subobscura and, therefore, their results are consistent with the role of inversions limiting gene flow (Rieseberg 2001; Navarro and Barton 2003b; Faria and Navarro 2010; Stevison et al. 2011). Because there are large chromosomal areas free of inversions fixed between D. madeirensis and D. subobscura (Larruga et al. 1983; Krimbas and Loukas 1984; Papaceit and Prevosti 1991), it remains unknown whether gene flow occurs between the two species at other genomic locations. To investigate this possibility, we analyzed polymorphism and divergence in 26 regions randomly distributed throughout the genome of D. madeirensis and D. subobscura, including mtDNA and the Y chromosome. Using phylogenetic analyses and population genetic methods that incorporate the stochastic variance of the coalescent process under an isolation-with-migration (IM) model (Hey and Nielsen 2004), we find evidence of gene flow between the two species at some nuclear loci, while the mitochondrial genome has remained distinct.

Materials and Methods

DROSOPHILA LINES

Thirty-three isofemale lines were established: 16 D. madeirensis, 16 D. subobscura, and one line of a more distantly related species endemic to the Canary Islands, Drosophila guanche, used as outgroup. All lines of D. madeirensis were established from individual females collected in Ribeiro Frio, Madeira Island, in an area of Laurisilva forest where D. subobscura is also found. Drosophila subobscura is distributed throughout the Palearctic region of the globe and our sample included lines from three geographic locations. Five D. subobscura lines come from Ribeiro Frio and represent a population in sympatry with D. madeirensis. Eight lines were collected in a pine forest near the small town of Camacha (Madeira Island) and represent a parapatric population. The final three lines, collected in Heidelberg (Germany), were obtained from the Drosophila Species Stock Center (https://stockcenter.ucsd.edu) and represent an allopatric population that shares no gene flow with D. madeirensis.

LOCI INVESTIGATED AND DNA SEQUENCING

We obtained DNA sequences for 26 regions in 15–16 lines of D. madeirensis, 13–16 lines of D. subobscura, and one line of D. guanche (Tables 1 and 2). Drosophila subobscura has five major acrocentric chromosomes, designated with vowels A, U, E, J, and O, and a dot chromosome (F). Chromosome A is largely homologous to the D. melanogaster X chromosome, U to 2L, E to 2R, J to 3L, O to 3R, and F to 4 (Lakovaara and Saura 1982; Segarra and Aguade 1992; Segarra et al. 1996; Papaceit et al. 2006). The regions studied include three genes located on the A chromosome, 16 genes on the autosomes (three on U, four on E, six on J, and three on O), five genes on the mitochondrial genome, and four genes on the Y chromosome. As the complete genome sequences of D. madeirensis and D. subobscura are not yet available, chromosomal locations of the nuclear regions under study were inferred based on the position of the orthologous regions in D. melanogaster and D. pseudoobscura, and on the conservation of Muller's elements across the genus Drosophila (Muller 1940; Sturtevant and Novitski 1941; Loukas et al. 1979; Schaeffer et al. 2008). Precise locations of genes within each chromosome arm are unknown.

Table 1

Divergence Data Summary

		Drosophila madeirensis–D. subobscura			D. madeirensis–D. guanche			D. subobscura–D. guanche
Locus	cds (bp)1	dN2	dS3	d4S4	dN2	dS3	d4S4	dN2	dS3	d4S4
Autosomal
Cch1	618	0.13	5.84	7.02	0.67	9.04	9.90	0.78	10.06	11.11
CG3609	609	1.14	5.29	4.75	1.29	15.09	16.05	0.24	14.46	15.99
CG8560	702	1.07	11.25	10.04	3.07	18.84	16.60	2.59	18.00	16.46
CG8768	498	0.78	5.99	6.49	1.86	14.06	16.42	2.14	12.17	13.83
Ef1β	564	0.08	2.25	2.04	0.00	2.84	2.21	0.08	2.53	1.95
eIF-5A	408	0.00	0.93	0.23	0.00	0.20	0.23	0.00	0.72	0.00
Hsp22	477	0.36	6.46	6.51	0.34	10.09	8.13	0.53	8.90	5.63
Jon65A	648	3.15	1.82	2.14	5.69	7.27	6.87	4.92	7.57	6.71
Nmt	648	0.00	2.59	3.48	0.00	7.25	7.73	0.00	7.24	6.98
Prosβ1	609	0.03	2.73	3.12	0.43	6.97	8.20	0.40	6.47	7.16
RpL32	402	0.00	1.27	1.70	0.00	1.81	3.13	0.00	1.16	2.39
RpLP1	336	0.00	2.29	3.43	0.00	1.63	2.33	0.00	1.86	3.24
RpS12	372	0.00	2.01	1.76	0.34	2.75	4.27	0.00	2.70	3.11
RpS26	342	0.00	0.40	0.40	0.00	0.93	1.66	0.00	1.04	1.88
RpS7	333	0.00	2.67	1.56	0.00	5.94	4.97	0.00	4.90	3.47
Vha26	375	0.00	0.65	0.40	0.00	0.97	0.25	0.00	1.30	0.14
X-linked
CG5004	384	0.02	1.50	1.24	0.02	2.16	4.18	0.00	2.38	4.87
sesB	519	0.00	0.04	0.00	0.47	1.60	2.14	0.47	1.56	2.14
YP3	531	0.54	0.57	0.27	1.74	3.48	4.80	1.33	3.42	4.54
Y-linked
kl2	693	1.28	2.04	1.21	1.55	4.22	7.47	1.03	5.25	6.79
kl3	285	0.00	3.08	4.50	0.53	3.94	4.11	0.53	5.01	4.16
Ory	351	0.97	2.15	2.14	1.29	6.04	4.35	0.97	3.83	2.16
Ppry	381	0.87	3.95	3.92	0.74	4.98	4.94	1.47	6.10	6.04
Average	11085	0.55	3.22	3.22	1.10	6.47	6.45	0.90	6.32	6.45
mtDNA	1446	0.79	10.46	9.66	1.40	18.60	19.62	1.21	15.78	16.56

Length of coding sequences (cds).

Number of nonsynonymous substitutions per 100 sites.

Number of synonymous substitutions per 100 sites.

Number fourfold degenerate substitutions per 100 sites.

Table 2

Polymorphism Data Summary

Locus	Species	n	Pol.	Rep.	Syn.	πN1	πSil2	Length (bp)
Autosomal
Cch1	Drosophila madeirensis	16	19	1	14	0.10	2.76	680
	Drosophila subobscura	16	30	2	20	0.09	4.89
CG3609	D. madeirensis	16	17	3	8	0.08	0.98	746
	D. subobscura	15	22	3	14	0.09	1.83
CG8560	D. madeirensis	16	34	7	28	0.46	4.89	711
	D. subobscura	16	31	4	27	0.19	5.11
CG8768	D. madeirensis	16	19	3	11	0.13	2.05	563
	D. subobscura	16	36	11	16	0.86	5.74
Ef1β	D. madeirensis	16	11	0	8	0.00	1.09	679
	D. subobscura	16	17	3	11	0.16	1.48
elF-5A	D. madeirensis	15	35	0	3	0.00	0.89	1273
	D. subobscura	16	23	0	0	0.00	0.55
Hsp22	D. madeirensis	16	15	5	5	0.17	1.26	570
	D. subobscura	16	8	1	6	0.04	0.48
Jon65A	D. madeirensis	16	8	1	1	0.05	0.91	747
	D. subobscura	16	12	3	7	0.24	0.64
Nmt	D. madeirensis	16	11	0	10	0.00	1.41	670
	D. subobscura	16	16	0	16	0.00	2.04
Prosβ1	D. madeirensis	15	8	2	5	0.06	0.64	694
	D. subobscura	13	24	0	11	0.00	4.34
RpL32	D. madeirensis	16	12	0	5	0.00	1.11	591
	D. subobscura	16	10	0	4	0.00	0.84
RpLP1	D. madeirensis	16	10	0	1	0.00	0.64	630
	D. subobscura	15	7	0	5	0.00	0.59
RpS12	D. madeirensis	16	8	0	2	0.00	0.33	762
	D. subobscura	13	10	0	2	0.00	0.86
RpS26	D. madeirensis	16	12	0	3	0.00	0.79	640
	D. subobscura	16	9	0	2	0.00	0.49
RpS7	D. madeirensis	16	7	0	2	0.00	0.30	673
	D. subobscura	16	7	0	4	0.00	0.23
Vha26	D. madeirensis	16	17	0	2	0.00	1.18	720
	D. subobscura	15	12	0	4	0.00	0.66
X chromosome
CG5004	D. madeirensis	16	10	1	3	0.04	0.64	698
	D. subobscura	16	7	0	3	0.00	0.46
sesB	D. madeirensis	16	5	0	1	0.00	0.38	752
	D. subobscura	16	2	0	0	0.00	0.13
YP3	D. madeirensis	16	13	3	4	0.23	0.71	682
	D. subobscura	16	12	1	4	0.06	1.06
Y chromosome
kl2	D. madeirensis	16	2	1	1	0.02	0.08	695
	D. subobscura	16	0	0	0	0.00	0.00
kl3	D. madeirensis	16	0	0	0	0.00	0.00	362
	D. subobscura	16	0	0	0	0.00	0.09
Ory	D. madeirensis	16	0	0	0	0.00	0.00	430
	D. subobscura	16	0	0	0	0.00	0.00
Ppry	D. madeirensis	16	0	0	0	0.00	0.00	437
	D. subobscura	16	1	1	0	0.04	0.10
mtDNA	D. madeirensis	16	31	5	20	0.08	0.65	2039
	D. subobscura	16	13	2	11	0.02	0.19

Nonsynonymous nucleotide diversity per 100 sites.

Synonymous and noncoding nucleotide diversity per 100 sites.

We designed polymerase chain reaction (PCR) primers using a variety of sequence sources and the software Primer3Plus (Untergasser et al. 2007). In most cases, primers were designed on sequences obtained initially from a D. subobscura cDNA library (J. M. Comeron, pers. comm.). Primers for the mitochondrial regions and the nuclear gene RpL32 (known as rp49 in D. subobscura) were based on published sequences of D. subobscura (Rozas and Aguade 1990; Volz-Lingenhohl et al. 1992), while Y-chromosome primers were designed on D. guanche sequences (Carvalho and Clark 2005). All regions under study were PCR-amplified using approximately 25 ng of total DNA isolated from individual male flies using the DNeasy kit for Blood and Tissue (Qiagen, Valencia, CA). PCR products were purified using the Wizard MagneSil clean-up system (Promega, Madison, WI). We used the Big Dye 3.1 chemistry (Applied Biosystems, Foster City, CA) to perform cycle sequencing reactions, which were purified with a 6.25% sephadex column and run in an ABI 3730 DNA Analyzer (Applied Biosystems). Both DNA strands were sequenced directly (no prior cloning). Sequences were inspected by eye and assembled using Sequencher 4.10.1 (Gene Codes, AnnAbor, MI). Multiple sequence alignments were obtained using the program Clustal-W (Larkin et al. 2007) and the program Phase 2.1 (Stephens and Scheet 2005) was used to reconstruct haplotypes from sequences with heterozygous sites. We screened all lines for the presence of the endosymbiont Wolbachia using diagnostic wsp PCR primers (Mateos et al. 2006) as well as positive and negative controls and found no evidence of infection in any of our lines. All sequences newly reported in this study were deposited in GenBank/EMBL/DDBJ database libraries under accession numbers KC622351–KC623163.

ANALYSIS OF POLYMORPHISM AND DIVERGENCE

Analyses of polymorphism were conducted using the program DnaSP v.5 (Librado and Rozas 2009). We reconstructed gene genealogies using the neighbor-joining algorithm (Saitou and Nei 1987) implemented in MEGA 5 (Tamura et al. 2011). The degree of confidence for each branch point in gene genealogies was assessed by the bootstrap method with 1000 replications (Felsenstein 1985). To estimate the number of synonymous (d), nonsynonymous (d), and fourfold degenerate (d4) changes per site, we applied Kumar's method as implemented in MEGA 5. These estimates were corrected for multiple hits at a site using the Kimura-2-p method (Kimura 1980). Divergence times were estimated following the approach proposed by Tamura et al. (2004). Briefly, this approach takes into account selective constraints associated with synonymous codon usage bias to estimate mutation distances (d) as a measure of neutral divergence according to the formula , where C is the base composition skew (BCS; Tamura and Kumar 2002), and η is the fraction of mutations eliminated per unit BCS. In Drosophila madeirensis, D. subobscura, and D. guanche, BCS is 0.044 and η is 3.33 (Tamura et al. 2004). Finally, we have assumed a nuclear mutation rate of 1.1 × 10−8 mutations per site and per year (Tamura et al. 2004) and five generations per year (Mestres et al. 2001).

DETECTION OF GENE FLOW

To determine whether D. madeirensis and D. subobscura have experienced gene flow after their initial split, we explored posterior probabilities of demographic parameters under an IM model using Markov chain Monte Carlo (MCMC) simulations of gene genealogies, as implemented in the program IMa2 (Hey and Nielsen 2004, 2007). According to the IM model, an ancestral population splits in two descendant populations, after which gene exchange between the two descendant populations may have occurred. The IM models used in the present study have the following demographic parameters: two (or three) population sizes, two migration rates, and the time since population splitting. The model assumes that the variation within the dataset is neutral (i.e., not affected by directional or balancing selection), free recombination between loci, and no intralocus recombination. We evaluated whether there is evidence in our dataset of possible violations of these assumptions. Overall, we found little evidence of positive selection. MK test (McDonald and Kreitman 1991) reveled only a significant excess of amino acid replacements fixed between species relative to synonymous changes in Jon65A (12 fixed replacements, two fixed synonymous changes, four polymorphic replacements, and 17 polymorphic synonymous changes; G = 16.33, P = 5×10−5). There is no evidence of positive selection for any locus when the Hudson-Kreitman-Aguade (Hudson et al. 1987) test is applied and limited evidence using Tajima's D test (Tajima 1989); no departure remained statistically significant after correcting for multiple tests (Table S1). To accommodate the assumption of no intralocus recombination, we constructed a dataset using the largest nonrecombining block of each gene, as determined using the four-gamete test (Hudson and Kaplan 1985). In addition to the original dataset containing all positions, one of only the largest nonrecombining block for each region, and one without positively selected amino acid replacements in Jon65A, we also created a dataset in which the gene with the strongest signal of gene flow (RpS26) was removed, and one that excluded the mtDNA. We conducted exploratory IMa2 runs with different numbers of chains and heating terms until we achieved good MCMC mixing as determined by large numbers of effective sample sizes (ESS) and the lack of long-term trends (zero autocorrelations). To ensure convergence, we completed two independent IMa2 runs, each with its own burn-in period and random number seed, and verified that they produced very similar estimates of the different demographic parameters. All parameter estimates were based on the combined results from two independent runs. We used 100 Markov chains (a = 0.96 and b = 0.90) with at least one million steps of burn-in and 100,000 genealogies per locus. The upper bounds of prior distributions of population mutation rates (θ), migration rates (m), and time since the populations split (t) were set following IMa2 manual: θ from U(0, 18.1), m from U(0, 0.55), and t from U(0, 7.23). Because our dataset includes mtDNA loci, we used the Hasegawa-Kishino-Yano mutation model (Hasegawa et al. 1985). To obtain estimates of the migration rates for each locus individually, we also used the program MIMAR (Becquet and Przeworski 2007). We conducted three independent runs, each with a different random seed and at least 4.5×106 steps and 5×105 burn-in steps. We estimated recombination rates using the program GenCo (Gay et al. 2007) setting gene conversion rates to zero. None of our estimates of recombination rates were significantly different than zero as determined by log-likelihood ratio tests.

Results and Discussion

DATING SPECIATION IN THE D. SUBOBSCURA COMPLEX

To obtain an overall view of speciation in the D. subobscura complex and provide context to the gene flow analysis, we first investigated the origin and divergence of the three closely related species D. madeirensis, D. subobscura, and D. guanche. Initial studies based on allozymes indicated unresolved phylogenetic relationships in the subobscura complex (Cariou et al. 1988). Pairwise comparisons of our 26 regions, however, show that levels of divergence at synonymous (d) and fourfold degenerate sites (d) between D. madeirensis and D. subobscura are consistently lower than in comparisons involving D. guanche, suggesting that the former two species share a more recent common ancestor and are therefore considered sister species (Table 1). To estimate time since population split, we used mutation distances (d; see Materials and Methods) and took into account ancestral polymorphism. The weighted average nucleotide diversity at fourfold degenerate sites is 0.015 in D. subobscura. If we assume that current levels of polymorphism in the Palearctic D. subobscura closely represent intraspecific variation in the ancestral species, the net d between D. subobscura and D. madeirensis for nuclear genes is 0.020 (d = 0.017) and between D. guanche and D. madeirensis/D. subobscura is 0.061 (d = 0.052). These estimates allow us to date the split between D. madeirensis and D. subobscura to ∼0.92 Mya and between D. guanche and D. madeirensis/D. subobscura to ∼2.75 Mya. Unlike D. simulans, D. mauritiana, and D. sechellia, which show a phylogeny with a basal polytomy (Kliman et al. 2000; McDermott and Kliman 2008; Garrigan et al. 2012), our results support previous findings indicating that D. madeirensis and D. subobscura are sister taxa and that D. guanche is a clear outgroup to these species (Gonzalez et al. 1990; Molto et al. 1994; Ramos-Onsins et al. 1998). The number of generations since the split between D. subobscura and D. madeirensis is approximately the same as the number of generations since the split of D. yakuba and D. santomea (Llopart et al. 2002; Llopart et al. 2005). This system is also formed by one mainland and one island endemic and shows evidence of gene flow as well (Llopart et al. 2005).

INTRASPECIFIC NUCLEOTIDE VARIATION

Overall, nucleotide diversity at silent (i.e., synonymous and noncoding, πSil) sites in autosomal and X-linked genes is reduced by ∼30% in the island endemic relative to the Palearctic species, with average estimates of 0.012 (95% CI: 0.0073–0.017; neutral coalescent simulations with no recombination) in D. madeirensis and 0.017 (95% CI: 0.010–0.024; neutral coalescent simulations with no recombination) in D. subobscura (Table 2). Although a reduction was anticipated in D. madeirensis relative to D. subobscura, it is not as large as the extreme difference in current census population sizes based on the geographic distribution of both species. This nonlinear relationship between census population sizes and levels of diversity was discovered in early surveys of polymorphism and is known as Lewontin's “paradox of variation” (Lewontin 1974; Leffler et al. 2012). One possible resolution of this paradox put forward by Leffler et al. (2012) is that species with small census population sizes, such as D. madeirensis, experience less genetic draft (i.e., weaker effects of selection at linked neutral sites) and therefore smaller effects of variation-reducing selection (Gillespie 2001). Another possibility is that in very closely related species with large differences in census population sizes, particularly endemics, the paradox can be partially explained by ancestral polymorphisms that predate the separation of the two species (Kliman et al. 2000; Ramos-Onsins et al. 2004; Llopart et al. 2005; McDermott and Kliman 2008; Nowell et al. 2011). Introgression with a close relative may also be a contributing factor to the maintenance of intraspecific variation (Ramos-Onsins et al. 2004; Llopart et al. 2005; Li et al. 2010). All these factors may play a role explaining the similar levels of observed diversity in these two species.

We also investigated silent nucleotide diversity on the X chromosome and autosomes. After correcting for potential differences between mutation rates using d, we observe an X/A ratio of diversity significantly smaller than 0.75 (expected under neutrality) in D. subobscura (0.0095/0.019; χ2 = 5.79, P = 0.016) but not in D. madeirensis (0.0099/0.013; χ2 = 2.3×10−4, P > 0.05). Similar probabilities are obtained without correction for possible differences in mutation rates (χ2 = 8.91, P = 0.0028 for D. subobscura and P = 0.07 for D. madeirensis). The result in D. subobscura resembles those seen in analyses of North American and cosmopolitan populations of D. melanogaster (Andolfatto 2001; Hutter et al. 2007; Langley et al. 2012), D. simulans (Begun and Whitley 2000; Begun et al. 2007), and D. pseudoobscura (Haddrill et al. 2010), as well as mammalian species (Baines and Harr 2007; Carneiro et al. 2010). The exception among Drosophila populations has been African D. melanogaster with an X/A ratio close to 1, significantly higher than the expected 0.75 under neutrality (Andolfatto 2001; Hutter et al. 2007; Singh et al. 2007; Langley et al. 2012).

In populations where background selection plays a significant role, an X/A ratio larger than 0.75 (as in African D. melanogaster) could be explained by a higher effective recombination rate on the X relative to autosomes (Charlesworth 2012b). In contrast, background selection is expected to play minor roles in populations with overall higher recombination rates, as seen in D. pseudoobscura and D. simulans, thus, predicting X/A ratios closer to 0.75 (Charlesworth 2012b). In D. subosbscura, accurate measures of crossover rates for all chromosomes are lacking, but partial genetic maps indicate higher rates than in D. melanogaster (Loukas et al. 1979; Pegueroles et al. 2010), which would explain limited effects of background selection. Notably, we observe an X/A ratio substantially lower than 0.75, which suggests that additional factors may be playing a role in reducing πSil on X-linked loci. These factors include (i) recent demographic changes in population size that differentially affected X-linked and autosomal variation (Hutter et al. 2007; Pool and Nielsen 2007) or (ii) a higher rate of adaptive evolution on the X than on autosomes (i.e., faster-X effect; Charlesworth et al. 1987; Begun and Aquadro 1993; Begun and Whitley 2000; Andolfatto 2001; Hutter et al. 2007; Singh et al. 2007; Vicoso and Charlesworth 2009; Stephan 2010; Connallon et al. 2012; Mackay et al. 2012). Additional multilocus or genomewide studies are required to distinguish between these possibilities.

Finally, we examined levels of polymorphism in mitochondrial and Y-linked loci (Table 2). Overall, patterns of intraspecific nucleotide variation in regions of no recombination are consistent with the effects of linked selection caused by background selection (Charlesworth et al. 1993a; Charlesworth 2012a) and/or hitchhiking (Maynard Smith and Haigh 1974; Kaplan et al. 1989; Wiehe and Stephan 1993; Stephan et al. 2006; Stephan 2010). As expected, the Y chromosome shows severely reduced πSil in both D. madeirensis and D. subobscura, with values close to those observed in some populations of D. melanogaster (Zurovcova and Eanes 1999; Larracuente and Clark 2013). In contrast, mitochondrial πSil is not as dramatically reduced, likely reflecting that, besides linkage effects, the mutation rate of mtDNA in Drosophila is higher than the nuclear mutation rate (Moriyama and Powell 1997; Haag-Liautard et al. 2008). Interestingly, intraspecific variation in the D. subobscura mitochondrial genome is threefold lower than that in D. madeirensis. This pattern could be explained by (i) a recent selective sweep (Braverman et al. 1995; Simonsen et al. 1995) in D. subobscura, and/or (ii) a reduction of the linkage effects of deleterious mutations in the mtDNA of the island endemic. Our analysis of the frequency spectrum of polymorphic sites shows a significant excess of rare variants in D. subobscura as estimated by Tajima's D statistic (Tajima 1989; D = −1.94, P < 0.05 from coalescent simulations), which yields some support for the recent selective sweep scenario. However, this result must be interpreted with caution as we detect an overall trend toward negative Tajima's D in many nuclear genes as well, which is consistent with purifying selection and/or population expansion (Fu 1997; Fay and Wu 1999; Gordo et al. 2002; Table S1).

VARIATION ACROSS LOCI IN GENE TREE TOPOLOGY

Gene flow between species is expected to affect some, but not all, genes in the genome and thus produce heterogeneity in gene tree topologies (Hey 2001; Machado and Hey 2003). To gauge whether divergence between D. madeirensis and D. subobscura has accumulated at similar rates in different regions of the genome, we obtained genealogies for all regions analyzed (Figs. 1 and S1). Ef1β, RpL32, and RpS26 produced gene trees in which sequences of both species intermingle. Support for gene flow between D. madeirensis and D. subobscura is particularly strong in RpS26, which shows the largest number of shared polymorphisms even though both nonsynonymous and silent sites are fairly conserved (Fig. 1A). This is also consistent with this gene showing the lowest FST (Hudson et al. 1992; Table S1). RpS26 is located on chromosome U, which shows the gene arrangement U1+2 in D. madeirensis (Krimbas and Loukas 1984; Papaceit and Prevosti 1991). This arrangement is present at a high frequency (∼60%) in current populations of D. subobscura from Madeira (Larruga et al. 1983), thus, providing opportunity for gene flow. Additional support for gene flow is observed in CG8768 and Prosβ1, which form two reciprocally monophyletic groups (bootstrap values >90%) that do not conform to the species phylogeny, as the D. madeirensis clade includes several sequences of D. subobscura (Fig. 1B and C). This suggests that gene flow from the former species to the latter has occurred in these regions of chromosome E (homologous to Muller's element C and arm 2R in D. melanogaster) that are homosequential on Madeira (Larruga et al. 1983; Krimbas and Loukas 1984; Papaceit and Prevosti 1991). The remaining 11 autosomal gene trees show either two distinct clades of sequences corresponding to D. madeirensis and D. subobscura or monophyly in only one of the two species (Fig. S1). The large variation in tree topology detected across loci may be attributable to a limited number of phylogenetically informative sites, likely due to functional constraints at protein level, gene flow, and/or incomplete lineage sorting (Kliman et al. 2000; Machado et al. 2002).

Figure 1

Neighbor-joining trees reconstructed from total nucleotide variation. Bootstrap values were obtained after 1000 replicates and values higher than 85% are shown. Drosophila madeirensis, D. subobscura, and D. guanche are indicated with filled circles, open squares, and filled triangles, respectively. (A) RpS26, (B) CG8768, (C) Prosβ1, and (D) mtDNA.

In contrast to the autosomes, we found high levels of genetic differentiation between the two species on sex chromosomes, which is consistent with the presence of two inversions fixed between D. madeirensis and D. subobscura on the X chromosome (Krimbas and Loukas 1984; Papaceit and Prevosti 1989, 1991; Khadem et al. 2011). These results are in agreement with reduced rates of gene flow on sex chromosomes compared to autosomes and with previous findings in multiple species (Tucker et al. 1992; Dod et al. 1993; Payseur et al. 2004; Llopart et al. 2005; Macholan et al. 2007; Geraldes et al. 2008a, 2008b; Teeter et al. 2008; Garrigan et al. 2012). The leading explanation for limited gene flow in sex chromosomes is the disproportionally large effect of the X chromosome on hybrid male sterility (Coyne and Orr 1989; Coyne 1992; Llopart 2012; Meisel et al. 2012), which is also detected in these species (Khadem and Krimbas 1991), and the potential contribution of the Y chromosome (Coyne 1985; Coyne et al. 2004; Geraldes et al. 2008b; Sackton et al. 2011). The smaller effective population size of the sex chromosomes that results in faster sorting of lineages could also be a contributing factor.

The mitochondrial gene tree shows also two distinct reciprocally monophyletic groups (Fig. 1D; bootstrap >90), clearly indicating that no gene flow has occurred recently in the mitochondrial genome between these two species. In Drosophila, this observation is unexpected because when the molecular signature of gene flow is detected, it is often seen in the mitochondrial genome (Ballard and Whitlock 2004). The reasons for the typical preferential introgression of mitochondrial over nuclear genes remain poorly understood, but several hypotheses have been proposed. One possibility is that the lack of physical linkage with nuclear DNA, the small genome size, and the strong functional constraints of mitochondrial protein-encoding genes could facilitate the exchange of the mitochondrial genome between species (Barton and Jones 1983; Ferris et al. 1983; Powell 1983; Ballard 2000; Sota and Vogler 2001; Nabholz et al. 2013). A second possibility is that within the framework of Hill–Robertson effects and Muller's ratchet (Lynch and Gabriel 1990; Charlesworth et al. 1993b; Comeron and Kreitman 2002; Gordo et al. 2002; Hadany and Comeron 2008; Charlesworth 2012a), mitochondrial genomes with fewer deleterious mutations in closely related species could be selected for and rapidly replace the “damaged” native haplogroups (Rieseberg 2009). However, there are also evolutionary reasons for mtDNA to remain distinct in the face of species hybridization. For instance, species-specific cytonuclear interactions may have already evolved at the time of secondary contact, thus, prohibiting mtDNA gene flow (Blier et al. 2001; Rand et al. 2004; Dowling et al. 2008; Burton and Barreto 2012; Osada and Akashi 2012; Meiklejohn et al. 2013). In D. madeirensis and D. subobscura, significant linkage disequilibrium has been detected between certain mtDNA haplotypes and chromosomal arrangements (Castro et al. 1999; Oliver et al. 2002), which could have minimized the likelihood of mtDNA introgression via hybridization in these species.

SIGNIFICANT POSTSPLIT GENE FLOW HAS OCCURRED FROM D. SUBOBSCURA INTO D. MADEIRENSIS

To determine whether D. madeirensis and D. subobscura have experienced gene flow and obtain estimates of demographic parameters, we fit our multilocus data to an IM model using the software IMa2 (Nielsen and Wakeley 2001; Hey and Nielsen 2004). Although we sampled both Heidelberg and Madeira populations of D. subobscura, we found no consistent evidence of genetic differentiation (Hudson 2000; Table S2), thus, we simulated an IM model with only two descendant populations (D. madeirensis and D. subobscura). Multiple independent runs show well-defined distributions of the posterior probabilities of all parameters, indicating that our data contain sufficient information to parameterize the model confidently. To estimate the time since the split between D. madeirensis and D. subobscura, we examined the marginal posterior probability distribution of the parameter t, which shows a sharp well-defined peak at 3.344 (Fig. 2A). We can convert this maximum likelihood (ML) estimate of time to years following Won and Hey (2005) and using the mutation rate (substitution rate) obtained from the comparison between D. guanche and D. subobscura/D. madeirensis. We estimated the split between D. guanche and the two sister species to have occurred ∼2.75 Mya (see “Dating speciation in the D. subobscura complex”), which corresponds to an average net divergence per locus of 23.12 and results in a rate per locus and per year of 4.20 × 10−6. We determined that the ancestor of the sister species colonized Madeira ∼0.80 Mya (95% highest posterior density (HPD) interval: 0.63–1.0 Mya), well in agreement with the time obtained using net nucleotide divergence at fourfold degenerate sites (see above). We also obtained ML estimates of the population mutation rate parameter, θ (Fig. 2B) and found them to be very similar in both species (θmad = 7.76, 95% HPD interval: 6.63–9.22 in D. madeirensis and θsub = 7.32, 95% HPD interval: 6.11–8.95 in D. subobscura). The precise estimates of these parameters may be slightly affected by population expansion (i.e., negative values of Tajima's D; Table S1; Strasburg and Rieseberg 2010). However, the similarity between species suggests that both D. madeirensis and D. subobscura have comparable effective population sizes, which is consistent with our observed estimates of nucleotide diversity (see above).

Figure 2

Smoothed marginal posterior probability distributions for (A) time since population split, (B) effective population sizes, and (C) migration.

Examination of the patterns of intra- and interspecific nucleotide variation across loci indicates heterogeneity in the number of shared and exclusive polymorphisms as well as in the number of fixed differences (Table 3). This heterogeneity strongly suggests the presence of gene flow at some regions of the genome but not others. ML estimates of the migration rates (m) are ms>m = 0.015 (95% HPD interval: 0.00082–0.041) from D. subobscura to D. madeirensis and mm>s = 0.0069 (95% HPD interval: 0–0.36) in the reverse direction (Fig. 2C). To determine whether these estimates were significantly different than zero, a log-likelihood-ratio test based on the comparison to a condition where m is set at zero was applied (Hey and Nielsen 2007; Hey 2010). Although ML estimates of the migration rates are nonzero in both directions, only gene flow from D. subobscura to D. madeirensis is significantly different than zero (LLRrun1 = 4.431 and LLRrun2 = 4.653, P < 0.05). In the reverse direction m = 0 is compatible with our data (LLRrun1 = 0.137 and LLRrun2 = 0.331, P > 0.05). We conclude that there is detectable, albeit limited, flow of nuclear genes between the two sister species after they split from the ancestral population, with a particularly strong signal of gene movement from D. subobscura to D. madeirensis. The significant migration from D. subobscura into D. madeirensis is maintained when (i) mitochondrial sequences are excluded from the dataset (LLR = 3.541, P < 0.05), (ii) only the longest nonrecombining blocks are used (LLR = 8.227, P < 0.05), (iii) codons with positively selected amino acid replacements in Jon65A are excluded (LLR = 4.04, P < 0.05), and (iv) when the ancestral population is not confined to be equal to the largest descendant population (LLR = 3.828, P < 0.05; Fig. S2).

Table 3

Shared and Fixed Variation between Drosophila madeirensis and D. subobscura

Chromosome	Locus	Shared Polymorphisms1	Fixed Differences	Exclusive to D. madeirensis	Exclusive to D. subobscura
U	CG3609	0 (0.524)	14	17	23
	RpLP1	0 (0.078)	11	7	7
	RpS26	7 (0.016)	0	5	2
E	CG8768	4 (0.909)	0	16	32
	elF-5A	0 (0.610)	9	37	21
	Ef1β	1 (0.221)	0	10	15
	Prosβ1	2 (0.199)	1	6	23
J	Jon65A	0 (0.145)	14	9	12
	CG8560	8 (0.911)	6	27	24
	Nmt	1 (0.224)	1	10	15
	Hsp22	0 (0.211)	8	15	8
	RpS12	0 (0.115)	29	8	11
O	Cch1	6 (0.476)	2	12	27
	RpL32	1 (0.154)	0	13	7
	RpS7	0 (0.073)	2	7	7
	Vha26	0 (0.300)	6	18	12
A	CG5004	0 (0.100)	6	10	7
	sesB	0 (0.013)	2	5	2
	YP3	1 (0.176)	2	12	10
Y2		0 (0.003)	28	2	3
mtDNA3		1 (0.193)	39	31	12

The number of expected parallel mutations are indicated in parentheses (Kliman et al. 2000).

Partial sequences of kl2, kl3, ORY, and PprY were combined into a single Y-linked region.

Partial or complete sequences of COXIII, ND3, ND5, Cytb, and ND1 were combined into a single mitochondrial region.

Our results are at odds with a recent study that reported no clear evidence of introgression between these species (Khadem et al. 2012). Two fundamental differences can reconcile these two studies. First, Khadem et al. (2012) investigated six nuclear regions, while in the present study we analyzed 26 nuclear and mitochondrial regions and thus have greater statistical power to detect gene flow. Second, and possibly more important, the regions studied by Khadem et al. (2012) were all associated with the O3 inversion, which is monomorphic in D. madeirensis and not present in D. subobscura. (D. subobscura shows the complex inversion O3+4). Chromosomal inversions fixed between species have been proposed to prevent gene flow and serve as “traps” for speciation genes (Noor et al. 2001; Rieseberg 2001; Sota and Vogler 2001; Navarro and Barton 2003b; Noor et al. 2007; Stevison et al. 2011). In contrast, our study includes loci randomly distributed through the genome, several of which are located on chromosomes that are colinear in D. madeirensis and D. subobscura and thus represent regions where gene flow will not be restricted.

RPS26 IS THE PRIMARY, BUT NOT THE SOLE, SOURCE OF INTROGRESSION

To determine which genes may be “trespassing” species boundaries, we applied the software MIMAR (Becquet and Przeworski 2007). The gene RpS26 shows an unusual pattern, with a very high estimate of the migration rate between D. madeirensis and D. subobscura (Fig. 3). This result is in agreement with those found when assessing levels of shared polymorphisms and gene trees of the individual loci (see above). Standard neutrality tests reveal no evidence of positive selection operating in this gene suggesting that gene flow in RpS26 does not represent a case of recent adaptive introgression, or the acquisition of beneficial alleles from closely related species (Arnold 1997, 2006; Whitney et al. 2006; Arnold and Martin 2009; Rieseberg 2009). To investigate the effect of RpS26 on the significant gene exchange discovered in our multilocus analysis, we excluded this gene from our dataset. Upon removal of RpS26, no significant gene flow was detected in either direction (all LLR < 0.85, P > 0.05). This result suggests that RpS26 appears to be the locus with the strongest signature of gene flow between the two sister species. However, nonzero estimates of the migration rates remain slightly more likely than zero estimates (ms>m = 0.0030, P = 0.049 vs ms>m ∼ 0, P = 0.036), which suggest that other loci may contribute as well to the overall signal obtained from the multilocus analysis.

Figure 3

Maximum likelihood (ML) estimates of migration rates obtained with MIMAR for individual nuclear regions.

Conclusions

Overall, our multilocus data fit better with a speciation model that incorporates migration between D. madeirensis and D. subobscura than with an alternative model of complete isolation. This observation indicates that in spite of abundant polymorphic inversions on all major chromosomes of D. subobscura, gene flow has occurred after the split of the two sister taxa. The strongest signal of gene flow is detected in autosomal regions that show the same arrangement in D. madeirensis as current Madeira populations of D. subobscura. There is no evidence of gene flow on sex chromosomes that is possibly associated with the disproportionally large X effect on hybrid male sterility seen in these species. In addition, X-linked loci reveal a significant deficit of silent nucleotide diversity in D. subobscura. This pattern contrasts with the X/A ratios seen in some species of the melanogaster subgroup and raise the interesting possibility that the X chromosome may experience higher rates of adaptive evolution than the autosomes, likely mediated by higher recombination rates. Although larger datasets are required to further study this possibility, we put forward that species in the “obscura” group, such as D. subobscura, are ideal systems to study the interplay of recombination, dominance, and natural selection. Finally, the mitochondrial genome appears to be resilient to gene flow in D. madeirensis and D. subobscura despite clear evidence of movement of nuclear genes, a trend that does not conform to the norm in Drosophila.