Gene flow between Drosophila yakuba and Drosophila santomea in subunit V of cytochrome c oxidase: A potential case of cytonuclear cointrogression.

Introgression is the effective exchange of genetic information between species through natural hybridization. Previous genetic analyses of the Drosophila yakuba-D. santomea hybrid zone showed that the mitochondrial genome of D. yakuba had introgressed into D. santomea and completely replaced its native form. Since mitochondrial proteins work intimately with nuclear-encoded proteins in the oxidative phosphorylation (OXPHOS) pathway, we hypothesized that some nuclear genes in OXPHOS cointrogressed along with the mitochondrial genome. We analyzed nucleotide variation in the 12 nuclear genes that form cytochrome c oxidase (COX) in 33 Drosophila lines. COX is an OXPHOS enzyme composed of both nuclear- and mitochondrial-encoded proteins and shows evidence of cytonuclear coadaptation in some species. Using maximum-likelihood methods, we detected significant gene flow from D. yakuba to D. santomea for the entire COX complex. Interestingly, the signal of introgression is concentrated in the three nuclear genes composing subunit V, which shows population migration rates significantly greater than the background level of introgression in these species. The detection of introgression in three proteins that work together, interact directly with the mitochondrial-encoded core, and are critical for early COX assembly suggests this could be a case of cytonuclear cointrogression.

Speciation is the process by which one interbreeding population is split into two descendant populations that, in time, acquire reproductive isolation (RI) (Mayr 1963). When RI is not complete, sympatric species can still effectively exchange genetic information through natural hybridization, also called introgression or introgressive hybridization (Anderson and Hubricht 1938; Riley 1939). Initially, it was believed that hybridization between different animal species was a rare accident. We now know, however, that hybridization is much more common in nature than previously suspected and introgression has been studied in many animal systems, including primates (Ferris et al. 1983; Hey 2003; Becquet and Przeworski 2007; Hey 2010a), mice (Geraldes et al. 2008; Teeter et al. 2008; Staubach et al. 2012; Janousek et al. 2015; Liu et al. 2015), birds (Ellegren et al. 2012; Alcaide et al. 2014; Poelstra et al. 2014), butterflies (Bull et al. 2006; Dasmahapatra et al. 2007; Mallet et al. 2007; Putnam et al. 2007; Giraldo et al. 2008; Pardo‐Diaz et al. 2012; The Heliconius Genome 2012; Martin et al. 2013; Nadeau et al. 2014; Kronforst and Papa 2015), mosquitoes (Besansky et al. 1994; Marsden et al. 2011; Weetman et al. 2012; Clarkson et al. 2014; Neafsey et al. 2015; Norris et al. 2015), and Drosophila (Machado et al. 2002; Hey and Nielsen 2004; Llopart et al. 2005b; Bachtrog et al. 2006; Kulathinal et al. 2009; Nunes et al. 2010; Garrigan et al. 2012; Herrig et al. 2014).

The study of introgression is important for at least three reasons. First, understanding introgression and the mechanisms behind it is paramount in determining the evolutionary independence of taxa. Second, introgression can be a source of species’ adaptations (Arnold 1997; Arnold 2007). Through adaptive introgression, species can potentially benefit from genetic variation that arose in another species (Castric et al. 2008; Kim et al. 2008; Rieseberg 2011; Song et al. 2011; Pardo‐Diaz et al. 2012; Staubach et al. 2012). Additionally, it has been suggested that the acquisition of new traits may be much more common via introgression than by de novo mutation, particularly with respect to traits that depend on multiple substitutions (Rieseberg 2009). Finally, gene products that function well in their own species’ genetic background can generate incompatibilities when interacting with products from another species and result in hybrid sterility and inviability (i.e., DM incompatibilities) (Bateson 1909; Dobzhansky 1937; Dobzhansky and Queal 1938; Muller 1942; Coyne and Orr 2004). The study of introgression can provide insights into the genomic location of DM incompatibilities, local adaptation, or alternatively, uncover candidates for cointrogression.

In regions of sympatry where species meet and hybridize, the potential for introgression exists. These hybrid zones are “natural laboratories” where new genic combinations are generated and tested by selection (Hewitt 1988). Unfortunately, occurrences of these hybrid zones in genetically amenable model systems, like Drosophila, are not very common. A hybrid zone of two Drosophila species (D. yakuba and D. santomea) was discovered by Daniel Lachaise on the volcanic island of São Tomé and later characterized (Lachaise et al. 2000; Llopart et al. 2005b, a). This hybrid zone is located on the slopes of Pico de São Tomé where D. yakuba is found at low elevations, below 1450 meters, while D. santomea stays in shady regions covered by the rainforest, between 1150 and 2024 m (Lachaise et al. 2000; Llopart et al. 2005a). At middle elevations, the two species meet and hybridize creating a hybrid zone. Hybrids have been reported at a frequency of approximately 1–3%, as estimated by the presence of pigmentation patterns intermediate between the two parent species (Llopart et al. 2002) and confirmed by SNP genotyping (Llopart et al. 2005a). This frequency is ∼100 times higher than the reported frequency of hybridization between D. pseudoobscura and D. persimilis, another common pair used to study speciation (Dobzhansky 1973). D. yakuba and D. santomea began diverging genetically approximately ∼400,000 years ago (Cariou et al. 2001; Llopart et al. 2002; Llopart et al. 2005b) but are still able to produce fertile hybrid females (Coyne et al. 2004), which makes them capable of backcrossing in nature (F1 hybrid males are completely sterile), and are therefore an ideal system for studying introgression.

Previous work based on the analysis of 47 complete mitochondrial genomes indicates that mtDNA introgression is recurrent in this species system (Llopart et al. 2014). Studies have shown that the D. santomea mitochondrial genome has been completely replaced by that of D. yakuba (Llopart et al. 2005b; Bachtrog et al. 2006; Llopart et al. 2014). This replacement is of particular interest because it raises the possibility that nuclear‐encoded proteins interacting with mitochondrially encoded proteins may have cointrogressed. The proteins encoded in the mitochondrial genome are heavily involved in oxidative phosphorylation (OXPHOS), a metabolic pathway that is responsible for the production of most of the energy in an animal's cells and is remarkably important to cell survival (Smeitink et al. 2001). Five protein complexes, which are all localized in the lipid bilayer of the inner mitochondrial membrane, play key roles during OXPHOS. While Complex II is composed solely of proteins encoded by nuclear DNA (nDNA), the four remaining complexes are composed of interacting proteins that are encoded in both nDNA and mtDNA (Tripoli et al. 2005). Complex IV of OXPHOS, often referred to as cytochrome c oxidase or COX, is of special interest because its core is composed entirely of mitochondrial‐encoded subunits (I, II, and III) and is surrounded by multiple nuclear‐encoded proteins with which it must interact (Saraste 1999) (Fig. 1). These functional interactions in COX have been proposed to result in cytonuclear coadaptation in some species (Blier et al. 2001; Rand et al. 2004; Osada and Akashi 2012).

Figure 1

Cytochrome c oxidase (COX) subunits of oxidative phosphorylation (bovine). Mitochondrially encoded subunits (I–III) form the core of COX. Nuclear subunits Va and Vb showed evidence of cointrogression in D. yakuba and D. santomea. Adapted with permission from http://www.genome.jp/kegg/pathway.html, map00190 (Kanehisa and Goto 2000).

Although the mitochondrial genome represents only a small fraction of the genetic material of an organism, the multiple interactions between mitochondrial and nuclear partners provide unique opportunities to study DM incompatibilities and their effects in hybrids (Rand et al. 2004; Burton and Barreto 2012). For example, the F2 hybrid males of the parasitoid wasps Nasonia giraulti and N. vitripennis experience increased mortality and this hybrid breakdown is also associated with disruption of mitochondrial function and a reduction in ATP production (Breeuwer and Werren 1995; Ellison et al. 2008; Niehuis et al. 2008; Gibson et al. 2013). Interpopulation hybrids of the copepod Tigriopus californicus show variable degrees of reduced fitness that map to mitochondrial‐nuclear incompatibilities (Edmands and Burton 1999; Burton et al. 2006; Ellison and Burton 2006, 2008; Burton and Barreto 2012). A recent study of the genomic trajectory of hybrid swarms in this same species shows a decrease in cytonuclear mismatch over time, which is consistent with coordinated changes of mitochondrial and nuclear interactors (Pritchard and Edmands 2013). In Drosophila, the search for DM incompatibilities between species involving the mitochondrial genome has produced mixed results, possibly reflecting the examination of different species/strains (Sackton et al. 2003; Montooth et al. 2010; Hoekstra et al. 2013; Meiklejohn et al. 2013).

Here, we evaluate levels of introgression between D. yakuba and D. santomea in the COX nuclear genes using multilocus approaches based on maximum‐likelihood population genetics methods. Our results revealed significant gene flow from D. yakuba to D. santomea for the entire COX complex and more specifically for COX subunit V. Most important, the population migration rate estimated for this subunit is significantly greater than the background migration rate obtained from 25 genes unrelated to OXPHOS and randomly distributed across the genomes of the two species (Llopart et al. 2005b). This supports the notion that introgression in subunit V cannot be fully explained by baseline levels of gene flow between D. yakuba and D. santomea. The signature of gene flow in subunit V is not limited to a single gene. We detected the same distinct pattern typical of introgression in the three nuclear genes encoding subunit V (CoVa, CoVb, and CG11043): no fixed difference between species and abundant shared variation, a pattern not seen in any of the other COX loci analyzed. Because these three genes are part of the same main subunit, which interacts with the mitochondrial‐encoded core early on in COX assembly, we suggest that this pattern is the result of coordinated events of nuclear and mitochondrial cointrogression.

Materials and Methods

FLY LINES STUDIED

We established a total of 33 isofemale lines of Drosophila in the laboratory, including 16 D. yakuba, 16 D. santomea, and one D. teisseri. Allopatric samples of D. yakuba included two lines from Tanzania and one from Zimbabwe. The D. yakuba sample also included lines collected from regions of sympatry or parapatry with D. santomea, eight from the São Nicolau waterfall (São Tomé) and five from or near the hybrid zone of Pico de São Tomé. All D. santomea samples were collected from regions of sympatry with D. yakuba and included four from Rio Quija in southwest São Tomé, five from Bom Sucesso, the field station on Pico de São Tomé, and seven from or near the hybrid zone. The D. teissieri line was established from a female collected in Brazzaville, Republic of Congo. All flies were maintained in the laboratory at 24°C on standard cornmeal‐yeast agar medium and kept in 12 hour light/dark cycles.

LOCI ANALYZED

We analyzed patterns of polymorphism and divergence in the 12 genes encoding COX nuclear subunits (14 regions total because we included the 5’ and 3’ noncoding flanking regions of the CoVa gene; Table 1). We chose COX due to the intimate nature of protein–protein interactions between the nuclear‐encoded subunits and the three mitochondrial‐encoded subunits composing the complex. We also focused on this complex because it has the highest ratio of mitochondrial to nuclear subunits out of all the OXPHOS complexes, and it is the rate‐limiting step of OXPHOS. While ten subunits are found in COX in mammals, so far only eight have been identified in Drosophila, which lacks VIIb and VIII nuclear‐encoded subunits (Tripoli et al. 2005) (http://www.mitocomp.uniba.it/). We designed primers for each of these regions using the D. yakuba reference genome sequence obtained from the 12 Drosophila Genome Project (Clark et al. 2007) and Primer3Plus (Untergasser et al. 2007). Thirteen of the regions under analysis are autosomal and one is X‐linked (Table 1).

Table 1

Gene names, chromosomal locations, and subunit associations of nuclear genes of OXPHOS COX

Gene name in	Gene name in	Location in	Cytogenetic	Location in
D. melanogaster	D. yakuba	D. melanogaster	map1	D. yakuba	Subunit2
CG10396 3	GE15061	2R	41E5	Un	IV
	GE14691	2R	—	Un	—
CoIV	GE13288	2L	38A8	2R	IV
CoVa	GE24204	3R	86F9	3R	Va
CoVb	GE13966	2L	26E3	2L	Vb
CG11043	GE13977	2L	26E3	2L	Vb
CG30093	GE11738	2R	52D3‐D4	2R	VIa
levy	GE14301	2R	59E3	2R	VIa
CoVIb	GE17355	X	18E5	X	VIb
Cype	GE25385	2L	25D6	2L	VIc
CG9603	GE24854	3R	84F13	3R	VIIa
CG18193	GE25872	3R	84F13	3R	VIIa
CoVIIc	GE19330	2R	46D8‐D9	3L	VIIc

1cytogenetic map in D. melanogaster.

2Subunit association of OXPHOS COX (http://www.mitocomp.uniba.it/).

3Gene CG10396 is duplicated in the D. yakuba genome project sequence (Clark et al. 2007).

According to the D. yakuba genome project sequence, CG10396 is actually duplicated (denoted GE15061 and GE14691 in D. yakuba) and we were able to detect the duplication in every single D. yakuba line we studied. Whether these two copies remain functional is currently unknown. In contrast, we found no evidence of the CG10396 duplication in any of the D. santomea lines investigated. Due to the lack of divergence information with D. santomea, the duplicated copy of CG10396 (i.e., GE14691) was excluded from all our analyses. In addition, one line of D. santomea shows a premature stop codon in CG30093, resulting in a protein two amino acids short.

DNA EXTRACTION, SEQUENCING, AND DESCRIPTION OF INTER‐ AND INTRASPECIFIC VARIATION

We isolated DNA from single male flies from each line using the Qiagen DNeasy Blood and Tissue kit (Qiagen, Valencia, CA) and amplified the regions of interest using PCR. PCR products were then purified using the Wizard Magnesil PCR clean‐up system (Promega, Madison, WI). Sequencing reactions were performed using the Big Dye 3.1 chemistry (Applied Biosystems, Foster City, CA), cleaned‐up using a 6.25% Sephadex column, and run in an ABI 3730 DNA Analyzer (Applied Biosystems, Foster City, CA). To maximize accuracy, both DNA strands were sequenced and assembled into single contigs using Sequencher 5.0 (Gene Codes, Ann Arbor, MI). A pool of haplotypes was then reconstructed for each locus containing heterozygous sites using the PHASE 2.0 program excluding indels (Stephens et al. 2001). For each line a single haplotype was selected for further analysis from the pool based on highest confidence as calculated by PHASE 2.0 (Stephens et al. 2001).

To calculate the number of synonymous (d) and nonsynonymous (d) changes per site and d/d, we used the codeml program within the PAML v4.5 package (Yang 2007). We calculated the average nucleotide frequencies from third codon positions (F3×4), and we assumed a single d/d for all the lineages (Model 0) and among sites (NSites = 0). Gene trees were generated using MEGA 5 (Tamura et al. 2011), with bootstrap values (Felsenstein 1996) based on 10,000 replicates. Analyses of polymorphism were performed using the program DnaSP 5.1 (Librado and Rozas 2009).

DETECTION OF INTROGRESSION

To detect introgression between D. yakuba and D. santomea, we used the multilocus method based on the Isolation‐with‐Migration (IM) model (Nielsen and Wakeley 2001; Hey and Nielsen 2004) implemented in the program IMa2 (Hey and Nielsen 2004, 2007). This program generates simulations of gene genealogies using a Markov‐Chain Monte‐Carlo (MCMC) approach to obtain posterior distributions of population demographic parameters (Nielsen and Wakeley 2001; Hey and Nielsen 2004). The IM model assumes that one panmitic ancestral population divides into two descendant populations that may (or may not) have exchanged genes after their split. To determine whether a population model with gene flow fits better to the data than an alternative model with 0 migration, IMa2 performs likelihood ratio tests based on the maximal posterior probability for a migration parameter and the posterior probability when that parameter is set at zero (Nielsen and Wakeley 2001).

The IM model assumes that all detected variation is neutral (i.e., not affected by directional or balancing selection), no intralocus recombination and free recombination between loci (Hey and Nielsen 2004). To test the assumption of neutrality we applied Tajima's D (Tajima 1989), McDonald‐Kreitman (MK) (McDonald and Kreitman 1991), and Hudson‐Kreitman‐Aguade (HKA) (Hudson et al. 1987) tests. Tajima's D statistic rejected neutrality in only 2 out of 28 tests performed (14 regions analyzed in two species) (Table S1). None remained significant after multiple‐test correction. There was an overall trend toward negative values of Tajima's D, which indicated a deficit of polymorphisms at intermediate frequencies (Tajima 1989). These results potentially indicated a recent population expansion (Tajima 1989) and/or purifying selection (Fu 1997; Fay and Wu 1999; Gordo et al. 2002), in agreement with previous findings in the same pair of species (Llopart et al. 2005b; Bachtrog et al. 2006). MK tests detected no significant departure from neutral expectations (McDonald and Kreitman 1991) (Table S2). The mulitlocus HKA test (https://bio.cst.temple.edu/∼hey/software/software.htm#HKA) showed consistency in the polymorphism to divergence ratio across loci (P = 0.16; Table S3). From these findings combined, we conclude that overall there is no evidence of positive or balancing selection in our dataset. To address the issue of no intralocus recombination, IMa2 was run additionally using the largest sequence block of each region showing no evidence of intralocus recombination assessed using the four gamete type test in DNASP 5.1 (Librado and Rozas 2009).

IMa2 was run for all 14 COX regions, as well as groups of 3–5 loci corresponding to each main nuclear subunit of COX (IV, V, VI, and VII). We also used IMa2 to calculate migration rates for a subset of 25 OXPHOS‐unrelated genes randomly distributed throughout the D. yakuba/D. santomea genomes (Llopart et al. 2005b). These genes provide a background signal of introgression in these species. For all datasets the burn‐in period was at least 1 × 10+6 steps under the HKY model (Hasegawa et al. 1985) with 150 Markov chains (a = .99 b = .75). Prior values were set following the recommendations in the IMa2 documentation (Table S4). Mixing was assessed by monitoring the Effective Sample Sizes (ESS) and parameter trends over a run. Large ESS and lack of long‐term parameter trends (zero autocorrelations) indicate good mixing. To ensure convergence, two independent sets of simulations were obtained, each run for at least 1 × 10+5 genealogies. Posterior probabilities of the different parameters were calculated using both runs combined. To convert model parameters to more familiar demographic quantities (e.g., the population migration rate, 2N, where N is the effective population size and m is the migration rate per locus and per generation) we followed Won and Hey (2005). Model parameters and converted units are shown in Table S5.

To obtain posterior probability distributions of migration rates additional to those provided by IMa2, we applied the program MIMAR, which can incorporate intralocus recombination (Becquet and Przeworski 2007). The program GenCo was used to estimate intralocus recombination for each gene (Gay et al. 2007) (Table S6). For each dataset we performed at least three MIMAR runs with different seeds, each carried out independently utilizing a burn‐in period of 1 × 10+5 steps and at least 1 × 10+6 steps after burn‐in. The ranges of the prior distributions for MIMAR were U [0.002, 0.05] for the population mutation rates (in per base pair units), U [2.5 × 10+6, 4.5 × 10+6] for splitting times (in number of generations), and logarithm from U [4.54 × 10‐5, 2.72] for population migration rates. Contrary to IMa2, MIMAR requires information from a third species used as outgroup. We used D. erecta as outgroup instead of D. teissieri for two reasons. First, it has been proposed that D. teissieri has exchanged genes with D. yakuba–D. santomea, at least for the mitochondrial genome (Llopart et al. 2005b; Bachtrog et al. 2006). Gene flow between D. teissieri and D. yakuba can severely affect the inference of the ancestral/derived allele upon which MIMAR is based. Second, the whole genome of D. erecta is available as part of the 12 Drosophila genome project and the species constitutes a valid outgroup for D. yakuba and D. santomea (Clark et al. 2007). For the 14 COX loci and the 25 random genes, MIMAR runs produced acceptance rates that were extremely low (7.78 × 10‐6–3.27 × 10‐5) and thus considered unacceptable based on MIMAR documentation. MIMAR runs with individual main subunits showed good acceptance rates for subunits IV (0.19 – 0.21), VI (0.12 – 0.13), and VII (0.07 – 0.14), and acceptance rates that were marginally acceptable for subunit V (0.014 – 0.017).

Results and Discussion

INTER‐ AND INTRASPECIFIC VARIATION

To obtain an overall view of nucleotide variation associated with the 12 nuclear genes of COX in D. yakuba and D. santomea, we investigated between‐species divergence and within‐species polymorphism levels. Table 2 shows the number of nonsynonymous (d N) and synonymous (d S) changes per site between D. yakuba and D. santomea for each gene. The average d N of all 12 genes was consistent with previous estimates for this same pair of species (Llopart et al. 2005b) suggesting that the nuclear genes of COX are not particularly conserved when considered as a class. Levels of nucleotide polymorphism were also obtained. The average nucleotide heterozygosity at synonymous sites is 1.65% for D. santomea and 2.32% for D. yakuba (Table 3), indicating only a modest reduction in polymorphism in D. santomea relative to D. yakuba consistent with previous findings (Llopart et al. 2005b; Bachtrog et al. 2006). Because D. santomea is an island endemic and D. yakuba can be found throughout sub‐Saharan Africa, the two species are predicted to have very different census population sizes based on their current geographical distribution, and population genetics theory would predict a more substantial difference in intraspecific variation. This same trend of similar levels of polymorphism, however, has been seen previously in other closely related species with apparently large differences in census population size, and has been attributed, at least partially, to ancestral polymorphism and very recent species split (Kliman et al. 2000; Ramos‐Onsins et al. 2004; Llopart et al. 2005b; Leffler et al. 2012).

Table 2

Average divergence between D. yakuba and D. santomea at nonsynonymous (d) and synonymous (d) sites of OXPHOS COX nuclear genes

Locus1	dN	dS	dN/dS
CG10396	0.003	0.02	0.15
CoIV	0	0.03	0
CoVa	0	0.06	0
CoVb	0	0.03	0
CG11043	0.003	0.06	0.05
CG30093	0.01	0.20	0.05
levy	0	0.18	0
CoVIb	0	0.010	0
Cype	0	0.03	0
CG9603	0.005	0.09	0.06
CG18193	0.004	0.06	0.07
CoVIIc	0	0.05	0
All	0.002	0.07	0.03

1Gene symbols from D. melanogaster.

Table 3

Polymorphism data for D. yakuba and D. santomea in the nuclear genes of OXPHOS COX

Locus1	Species	n 2	π 3	θ w 4	θ syn 5	θ sil 6	θ N 7	Syn8	Rep9
CG10396	D. santomea	16	0.08	0.19	0	0.40	0.074	0	1
	D. yakuba	14	0	0	0	0	0	0	0
CoIV	D. santomea	15	0.80	0.99	1.17	1.67	0	5	0
	D. yakuba	16	0.41	0.49	2.07	0.84	0	9	0
CoVa	D. santomea	14	0.54	0.82	2.55	1.60	0	9	0
	D. yakuba	15	0.53	0.87	2.50	1.43	0.18	9	2
CoVa3’	D. santomea	15	1.16	1.15	‐	1.20	‐	‐	‐
	D. yakuba	15	0.59	1.12	‐	1.12	‐	‐	‐
CoVa5’	D. santomea	15	1.03	1.34	‐	1.34	‐	‐	‐
	D. yakuba	15	1.25	1.84	‐	1.84	‐	‐	‐
CoVb	D. santomea	15	1.71	2.11	2.72	2.64	0.11	7	1
	D. yakuba	14	1.78	1.90	0.80	2.5	0	2	0
CG11043 10	D. santomea	14	1.43	1.55	4.14	1.90	0.55	16	7
	D. yakuba	14	1.22	1.53	4.14	1.81	0.27	17	5
CG30093	D. santomea	13	1.27	1.47	2.81	1.88	0.91	6	6
	D. yakuba	15	2.60	2.68	9.07	4.12	0.57	20	4
levy	D. santomea	16	0.42	0.69	1.90	1.06	0	5	0
	D. yakuba	16	0.92	1.08	4.19	1.70	0	11	0
CoVIb	D. santomea	16	0.11	0.20	1.21	0.26	0	2	0
	D. yakuba	16	0.42	0.64	0.61	0.87	0	1	0
Cype	D. santomea	16	0.75	0.86	0.99	1.11	0.18	2	1
	D. yakuba	14	1.23	1.34	2.06	1.78	0	4	0
CG9603	D. santomea	16	0.76	0.98	1.85	1.47	0.15	4	1
	D. yakuba	16	0. 2	0. 34	0	0.47	0	0	0
CG18193	D. santomea	14	1.35	1.31	0.43	1.77	0.67	1	5
	D. yakuba	16	0.13	0.23	0.41	0.29	0.26	1	1
CoVIIc	D. santomea	16	0.33	0.53	0	0.67	0	0	0
	D. yakuba	14	0.28	0.33	1.93	0.42	0	3	0

1Gene symbols from D. melanogaster.

2Number of sequences analyzed.

3Nucleotide diversity (%).

4Watterson's (1975) estimate of heterozygosity (%).

5Heterozygosity (%) at synonymous sites.

6Heterozygosity (%) at silent (synonymous and noncoding) sites.

7Heterozygosity (%) at nonsynonymous sites.

8Number of synonymous polymorphisms.

9Number of amino acid replacement polymorphisms.

10A single codon with multiple changes at a site was excluded.

PATTERNS OF VARIATION IN GENE TREE TOPOLOGY

Introgression violates the simple bifurcating model of species divergence upon which phylogenetic analyses are based, and as such, is expected to leave a phylogenetic signal. We obtained gene trees for each locus individually using Neighbor‐Joining and maximum‐likelihood methods (Figs. S1 and S2). As expected for closely related species, we detect a variety of gene tree topologies, ranging from two reciprocally monophyletic clades with good bootstrap support (e.g., gene levy) to single groups with sequences from both species clustering together (e.g., CoVa) (Fig. 2). The observed heterogeneity in gene trees across loci is consistent with either introgression between species, incomplete lineage sorting, or high sequence conservation across species due to strong functional constraints. Additional analyses are thus required to determine whether there is evidence of introgression between D. yakuba and D. santomea in the COX loci.

Figure 2

Gene trees reconstructed with the Neighbor‐Joining algorithm and bootstrap values based on 10,000 replicates for levy and CoVa. Bootstrap values higher than 70% are shown. D. santomea, D. yakuba, and D. teisseri are shown with filled circles, shaded squares, and open triangles, respectively.

EVIDENCE OF GENE FLOW FROM D. YAKUBA TO D. SANTOMEA IN COX NUCLEAR GENES

To investigate whether there has been genetic exchange between D. yakuba and D. santomea in COX, we first determined the number of shared polymorphisms and fixed differences in each nuclear gene under study (Table 4). There was substantial heterogeneity in the distribution of shared polymorphisms and fixed differences across loci, which suggests the possibility of introgression at some but not all of the genes investigated. Shared variation can occur from (1) introgression (i.e., gene flow), (2) parallel mutations, or (3) retention of ancestral polymorphism due to incomplete lineage sorting in closely related species (Clark 1997). To rule out parallel mutation, we calculated the expected number of shared polymorphisms produced by independent convergent mutations in each of the two species (Kliman et al. 2000). In all cases, the number of observed shared polymorphisms was substantially greater than that expected from parallel mutation (Table 4). The other evolutionary source of shared variation in very closely related species is ancestral polymorphism, which can persist due to incomplete lineage sorting. However, the genealogies under introgression and incomplete lineage sorting are expected to have different branch lengths and topologies. Population genetic models, like the IM model (Nielsen and Wakeley 2001; Hey and Nielsen 2004) implemented in the software IMa2, can differentiate between introgression and incomplete lineage sorting using likelihood methods (Hey and Nielsen 2004, 2007).

Table 4

Shared and fixed variation in nuclear genes of OXPHOS COX

Locus1	L(bp)2	Sx1 3	Sx2 4	F5	SS 6	Main subunit
CG10396	942	4	0	7	0(0)	IV
CoIV	1138	31	16	6	0(0.44)	IV
CoVa	924	12	12	0	8(0.16)	V
CoVa3’	858	18	18	0	8(0.38)	V
CoVa5’	850	30	42	10	1(1.48)	V
CoVb	1487	64	61	0	22(2.63)	V
CG11043	1805	53	53	0	25(1.56)	V
CG30093	883	17	48	6	7(0.92)	VI
levy	816	11	22	5	5(0.30)	VI
CoVIb	789	4	14	2	1(0.07)	VI
Cype	932	18	28	2	4(0.54)	VI
CG9603	785	23	6	9	1(0.18)	VII
CG18193	850	26	4	8	1(0.12)	VII
CoVIIc	892	12	7	9	0(0.09)	VII

1Locus symbols from D. melanogaster.

2Length in base pairs.

3Polymorphisms exclusive to D. santomea.

4Polymorphisms exclusive to D. yakuba.

5Differences fixed between species.

6Number of shared polymorphisms (Expected number of shared polymorphisms from parallel mutation).

To determine whether the patterns of intra‐ and interspecific variation in the nuclear genes of COX could be fully explained in absence of gene flow between D. yakuba and D. santomea, we applied the likelihood ratio (LLR) tests (Nielsen and Wakeley 2001; Hey 2010b) implemented in IMa2 (Hey and Nielsen 2004, 2007). Simulations yielded contained posterior probabilities for all demographic parameters estimated, suggesting that our dataset contained sufficient information (Fig. 3). The LLR tests indicated that a population model in which there was gene flow from D. yakuba to D. santomea fits significantly better to our data (LLR = 7.07, P < 0.01 and LLR = 7.22, P < 0.01 for seeds 1 and 2, respectively) than a model with 0 migration. In contrast, there was no evidence of gene flow in the opposite direction (LLR = 1.03, P > 0.05 and LLR = 1.11, P > 0.05 for seeds 1 and 2, respectively). Qualitatively similar results were obtained but with significance in both directions when we used the largest block of each region with no evidence of intralocus recombination (LLRseed1 = 11.87, P < 0.001 and LLRseed2 = 12.71, P < 0.001 from D. yakuba to D. santomea, and LLR seed1 = 7.28, P < 0.01 and LLR seed2 = 7.25, P < 0.01 from D. santomea to D. yakuba). We conclude that the nuclear genes encoding COX show the molecular signature of introgression from D. yakuba into D. santomea.

Figure 3

The marginal posterior probability distributions for speciation parameters scaled by the neutral mutation rate and obtained by fitting the IM model to a 14‐region dataset using IMa2. (A) Population mutation rates, θ = 4N, where N is the effective population size and u is the neutral mutation rate per locus and per generation. (B) Rates of migration for each gene copy per mutation event, m = m/u, where m is the migration rate per locus and per generation. (C) Time since population split in terms of mutations, t = t u y, where t is time since population splitting in years and u y is the neutral mutation rate per locus and per year.

GENE FLOW IN MAIN SUBUNIT V

To assess whether introgression occurs in clusters or is scattered throughout the different genes of COX, we performed independent IMa2 runs for each main nuclear‐encoded subunit (IV, V, VI, and VII; Table 4). The number of regions associated with each main subunit ranges from two (subunit IV) to five (subunit V, 3 genes and two noncoding regions). A recent study evaluating the accuracy of IMa2 has shown that, even when analyzing only two loci, the false‐positive rate is smaller than 0.01 for species with divergence time greater than 2N generations, where N is the effective population size (Cruickshank and Hahn 2014). Based on the analysis of 25 loci unrelated to OXPHOS and randomly distributed across the D. yakuba/D. santomea genomes (Llopart et al. 2005b), the two species split from their common ancestor ∼4N generations ago (Table S5). This corresponds to ∼0.43 Mya, in agreement with estimates from other studies (Cariou et al. 2001; Llopart et al. 2002; Llopart et al. 2005b; Bachtrog et al. 2006). Subunit V, consisting of Va and Vb, showed a significant rejection of an isolation model (i.e., gene flow = 0) in favor of an alternative model with gene flow from D. yakuba to D. santomea (LLRseed1 = 5.11, 0.01< P < 0.05 and LLRseed2 = 4.89, 0.01< P < 0.05) but not in the opposite direction (LLRseed1 = 0, P > 0.05 and LLRseed2 = 0.0, P > 0.05). In contrast, nucleotide variation in main subunits IV, VI, and VII is compatible with a model with no gene flow between species (all LLR = 0, P > 0.05). Equivalent results were obtained when the largest block with no evidence of intralocus recombination was used (LLR = 3.33 – 3.42, P < 0.05 for migration from D. yakuba to D. santomea in subunit V and LLR = 0 – 0.23, P > 0.05 in all the other subunits). Altogether these results suggest that the strongest signature of introgression in COX resides in subunit V, the innermost polypeptide in the mitochondrial membrane that interacts directly with the three mitochondrially encoded subunits (I, II, and III) (Fig. 1).

All three genes composing main subunit V (i.e., CoVa, CoVb, and CG11043) show the same pattern of variation: no fixed differences between species and many shared polymorphisms (Table 4). One could reason that shared polymorphisms in a group of functionally connected genes may indicate relaxed selection, which would result in abundant ancestral polymorphism. However, several observations argue against this possibility in subunit V. First, long‐term relaxed functional constraints would produce fixed differences between species, but there are none. Second, d estimates for the three genes of subunit V analyzed in the phylogeny of D. yakuba, D. teisseri, D. erecta, D. simulans, D. sechellia, and D. melanogaster are not greater than those for other COX genes, indicating no evidence of relaxed functional constraints (Mann–Whitney U test z = 0.83, P = 0.48; Table S7). Note that in this analysis we did not use D. santomea, as introgression with D. yakuba would tend to reduce d. Altogether these observations suggest that the fraction of shared polymorphisms present in the three genes of COX subunit V is not due to relaxed selection.

The shared variation in all genes of subunit V supports the notion that introgression is not concentrated on a single gene but spread throughout its three components. The probability of detecting introgression in a nuclear gene in these species is ∼0.1 (5/49 genes) based on this study and two additional large‐scale polymorphism surveys (Llopart et al. 2005b; Bachtrog et al. 2006). If we assume that this probability is similar across the genome, our observation of all three components in one main subunit showing abundant shared variation and no difference fixed between species is unlikely to be the result of unrelated introgression events (P = 0.001, binomial test). Altogether our results suggest the possibility that the three genes of subunit V may have cointrogressed with the D. yakuba mitochondrial genome.

SUBUNIT V IS A CANDIDATE FOR CYTONUCLEAR COINTROGRESSION

Mitochondrial introgression appears to be recurrent in the D. yakuba and D. santomea system in which an early invasion of the D. yakuba mitochondrial genome fully replaced the D. santomea native haplotypes (Llopart et al. 2014). In addition, a previous survey of nuclear loci randomly distributed throughout the genomes reported limited but detectable introgression between the two species (Llopart et al. 2005b). We hypothesize that genes involved in cytonuclear cointrogression will show migration rates significantly greater than the background signal from random genes in D. yakuba and D. santomea. To test this hypothesis, we inspected posterior probability distributions of migration rates for each of the main COX subunits. Only the curve corresponding to migration rates from D. yakuba to D. santomea for subunit V has a nonzero peak (Fig. 4). Similar posterior probability distributions for population migration rates from D. yakuba to D. santomea were obtained using MIMAR, which incorporates intralocus recombination (Fig. S3).

Figure 4

The marginal posterior probability distributions for migration rates (m = m/u) from D. yakuba to D. santomea (black line) and in the opposite direction (gray line) obtained by fitting the IM model to datasets of 2–5 regions corresponding to each main subunit using IMa2 (see legend of Fig. 3 for a detailed explanation of the units of m).

Subunit V shows an estimate of the population migration rate from D. yakuba to D. santomea that is incompatible with 0 migration (2Nm = 0.40, 95%HPD: 0.079 – 3.98) and is more than one order of magnitude greater than the migration rate estimated from random genes (2Nm = 0.014, 95%HPD: 0.002–0.11). We further applied Mann–Whitney tests to random subsamples from the two posterior probability distributions of 2Nm using the number of genes within each category. This study indicates significantly greater migration rates for subunit V than for the random genes, with 99.88% and 98.7% of the subsamples showing significant Mann–Whitney tests at critical values of P = 0.05 and P = 0.01, respectively. In addition, we also applied a nonparametric test. We find that the ratio of shared to exclusive polymorphisms, which is an index of introgression, is significantly greater in the loci associated with subunit V than in all the other 34 genes combined (25 random genes and nine genes unrelated to COX subunit V) (Mann–Whitney test for small sample sizes, exact P = 0.007; Table S8). Equivalent results were obtained when the coding sequence of CoVa was combined with the 5’ and 3’ flanking sequences into a single region (Mann–Whitney test for small sample sizes, exact P = 0.011; Table S8), and when only autosomal genes were analyzed (Mann–Whitney test for small sample sizes, exact P = 0.020). All these results combined indicate that (i) the significant introgression in subunit V is not due to a single region, and (ii) the signal of introgression in subunit V stands out from the background signal of introgression in D. yakuba and D. santomea nuclear genes.

To gain insight into the potential molecular causes of the introgression of main subunit V, we investigated its role in OXPHOS. The basic functions of subunits Va and Vb (Cox6p and Cox4p respectively in yeast) have been previously characterized in mouse and yeast, and more recently in Drosophila (Khalimonchuk and Rodel 2005; Fontanesi et al. 2006; Galati et al. 2009; Ren et al. 2010). Additionally, the specific protein–protein interactions of these subunits have also been characterized. These subunits are unique for several reasons. First, they are the innermost subunits of COX, and are extrinsic to the lipid bilayer (Khalimonchuk and Rodel 2005). These subunits are also unique because of their role in the formation of the COX holoenzyme. Subunit Vb (genes CoVb and CG11043) is the primary nuclear subunit involved in assembly of COX (Khalimonchuk and Rodel 2005; Fontanesi et al. 2006; Galati et al. 2009) and interacts with all three mitochondrially encoded subunits. More specifically, subunit Vb has been shown to interact intimately with mitochondrial protein Cox1 to initiate COX assembly (Tsukihara et al. 1996; Khalimonchuk and Rodel 2005). This is the first interaction that takes place in COX assembly, and once the Cox1/CoVb subcomplex is formed, the remaining subcomplexes of the holoenzyme come together to form COX (Galati et al. 2009). Subunit Vb also physically interacts tightly with mitochondrialy encoded Cox3 (Tsukihara et al. 1996), is specifically regulated by the mitochondrial content of cardiolipin, and in turn enhances binding affinity of Cox2 for cytochrome c, which initiates COX activity (Galati et al. 2009). Cardiolipin is a mitochondrion‐specific phospholipid required for COX activity (Khalimonchuk and Rodel 2005). Previous studies have also shown that a loss of subunit Vb prevents assembly of COX (Galati et al. 2009). Likewise, subunit Va (gene CoVa) has been shown to be a major regulator of COX by blocking inhibition of COX activity, and also plays a very important role in COX assembly and regulation. After the primary Cox1/CoVb interaction takes place, the CoIV/CoVa heterodimer is necessary for association with Cox2. This heterodimer then interacts with Cox1 (Fornuskova et al. 2010).

We observe similar patterns of shared variation and no fixed difference between species in the three proteins that form COX subunit V (CoVa and CoVb components). This observation, together with the complete replacement of the D. santomea mitochondrial genome by that of D. yakuba (Llopart et al. 2005b; Llopart et al. 2014) and the unique roles of CoVa and CoVb in COX assembly, opens the possibility that the three genes of subunit V may have cointrogressed together with the mitochondrial genome. This is also supported by population migration rates for the entire subunit V that are significantly greater than the background signal of introgression in this system. However, none of the genes of subunit V shows the molecular fingerprints of a recent “hard” selective sweep (Maynard Smith and Haigh 1974; Kaplan et al. 1989; Sella et al. 2009), which has been invoked in some cases of adaptive introgression (Song et al. 2011; Pardo‐Diaz et al. 2012; Staubach et al. 2012; Brand et al. 2013; Hedrick 2013; Clarkson et al. 2014; Huerta‐Sanchez et al. 2014; Llopart et al. 2014). The examination of the distribution of polymorphisms across lines suggests that introgression in the genes of subunit V is not incomplete or partial, for all variation is shared and there is no evidence of diverged haplotypes (Fig. S4). Note that, because our sample of D. yakuba includes sequences from allopatric populations of Africa mainland, our finding suggests that the variation present today in the three genes of subunit V consists of actually D. yakuba haplotypes, as it is the case for the mtDNA. We put forward the possibility that multiple D. yakuba nuclear haplotypes of subunit V cointrogressed into D. santomea capturing a fraction of the D. yakuba standing variation, which would easily explain the presence of shared polymorphisms and no fixed differences (Fig. 5C). At least qualitatively, this scenario would produce patterns of variation within the recipient species reminiscent of “soft” sweeps (Hermisson and Pennings 2005; Pritchard et al. 2010; Orr and Unckless 2014), as seen in Mytilus (Fraisse et al. 2014). Although quantitative studies of how much gene flow between species is necessary for this scenario to be realistic, the high frequency of hybridization of D. yakuba and D. santomea in nature today (∼3%) (Llopart et al. 2005a) and the fertility of F1 hybrid females (Coyne et al. 2004) suggests that abundant introgression may be allowed in some areas of the genome.

Figure 5

Schematic representation of different cases of introgression between species. (A) Adaptive introgression with partial replacement of variation in the recipient species. (B) Adaptive introgression with complete replacement of variation in the recipient species. (C) Introgression of multiple haplotypes from the donor species capturing a significant fraction of the standing variation. Note that in both (A) and (C) shared variation will be detected in the recipient species.

Conclusions

Our multilocus analyses indicate that a population model with gene flow from D. yakuba to D. santomea fits significantly better to the patterns of nucleotide variation in COX nuclear loci than an alternative model of isolation. More specifically, there is strong evidence that the genes composing major subunit V (genes CoVa, CoVb, and CG11043) are the main source of gene flow, with migration rates from D. yakuba to D. santomea significantly greater than the background signal of introgression detected in these two species. Because all three proteins show the same pattern of shared variation and lack of fixed differences between species, are part of the same main subunit, and interact with the mitochondrial‐encoded core early on in COX assembly, we interpret this finding as a potential case of cointrogression. Based on what we know about the molecular function of these subunits, we suggest that cointrogression of these genes contributed to maintain proper activity of cytochrome c oxidase and ultimately OXPHOS.

Table S1. Tajima's D test for selection in OXPHOS COX loci.

Table S2. McDonald–Kreitman test for selection in OXPHOS COX genes.

Table S3. HKA test for selection in OXPHOS COX genes.

Table S4. Values of uniform priors used in IMa2 runs.

Table S5. Parameters and converted demographic quantities obtained with IMa2.

Table S6. Estimates of recombination by GenCo.

Table S7. Number of nonsynonymous (d) and synonymous (d) changes per site obtained with PAML in COX genes analyzed in the phylogeny of D. yakuba, D. teisseri, D.erecta, D. simulans, D. sechellia, and D. melanogaster.

Table S8. Exclusive, fixed and shared variation in D. yakuba and D. santomea for the COX loci and random genes.

Figure S1. Gene trees reconstructed with the Neighbor‐Joining algorithm with bootstrap values based on 10,000 replicates.

Figure S2. Gene trees reconstructed with the Maximum‐Likelihood algorithm with bootstrap values based on 10,000 replicates.

Figure S3. Posterior probability distributions of the population migration rates from D. yakuba to D. santomea in the different COX subunits obtained with MIMAR.

Figure S4. Nucleotide variation (singletons not shown) in the three genes (CoVa, CoVb, and CG11043) of COX subunit V.

Click here for additional data file.