Mitonuclear interactions and introgression genomics of macaque monkeys (Macaca) highlight the influence of behaviour on genome evolution

(a) Genomic consequences of the macaque social system

Papionin monkeys are compelling subjects for the study of interaction networks that span different genomic compartments because their sex-specific behaviours establish strong predictions. In these monkeys, females tend to be philopatric whereas males tend to leave their natal groups (e.g. [12]). Evolutionary repercussions of extreme female philopatry are evidenced by ancient mitochondrial coalescence times that predate speciation of rhesus macaques [13], pigtail macaques (Macaca nemestrina) [14], olive baboons (Papio anubis) [15] and possibly Sulawesi macaques [14]. By contrast, a deep mitochondrial coalescence time is not observed in hamadryas baboons, Papio hamadryas, but this species has an atypical social system for papionins with more frequent female migration [16].

High intraspecific diversity and interspecific divergence of mitochondrial DNA leads to the prediction that genetic incompatibilities may arise between nuclear-encoded genes that interact with highly diverged mitochondria-encoded genes with whom they have not co-evolved (reviewed in [17]). This could lead to positive selection for compensatory mutations in nuclear-encoded genes that interact with mitochondria-encoded genes. The sharp geographical structure of mitochondrial DNA also has the potential to foster co-adaptation of genes involved with mitonuclear interactions to local environmental conditions (reviewed in [17]). Together these factors could disfavour the exchange of genes involved with mitonuclear interactions, both within and among species [5].

We set out to test these expectations using genomic data from 29 individuals from eight species of macaque, including pigtail macaques (M. nemestrina; six individuals) from Sumatra and Borneo and seven species of Sulawesi macaque (23 individuals). Macaca nemestrina is a member of the silenus group of macaques, which is sister to the Sulawesi macaques [13,18]. These taxa are separated by the Makassar Strait which runs between Borneo and Sulawesi, and corresponds to a sharp biogeographical boundary that is demarcated by Wallace's Line [19]. Sulawesi macaque species are allopatrically distributed, and parapatric species hybridize in nature (e.g. [20]). Mitochondrial genomes of these species diversified more than 2 Ma and are substantially diverged ([14], electronic supplementary material). We first focused on 211 nuclear-encoded genes that interact with mitochondria-encoded genes or mitochondrial DNA with an aim of testing whether population differentiation, polymorphism, metrics of natural selection, and runs of homozygosity in genomic regions that contained these genes differed from regions with other genes, or regions with no genes. Then, with an aim of examining how gene flow is specifically influenced by these 211 genes or more generally by sex-specific demography, we tested whether patterns of gene flow differed in regions containing these genes compared to other genomic regions, and we also evaluated population differentiation and gene flow on the X chromosome and the autosomes.

(a) A subset of N_interact genes has atypically high genetic differentiation between species

Samples in this study have highly diverged mitochondrial genomes (electronic supplementary material, tables S1 and S2, figure 1; electronic supplementary material, [14]). We therefore expected high interspecies differentiation and genomic signatures of positive selection in nuclear-encoded proteins that interact with mitochondria-encoded proteins, rRNA enzymes or mitochondrial DNA (i.e. OXPHOS, ARS2, MRP and REP genes—hereafter N_interact genes) [21]. To test this, we began by quantifying F_ST in 100 kilobase pair (kb) genomic windows in 10 pairwise comparisons between the five species with genomic data from at least three individuals each, as defined in Methods, with the expectation that this statistic (which quantifies population differentiation on a scale from 0 to 1) should be higher in windows that contain N_interact genes (hereafter N_interact windows). We considered the autosomes and the X chromosome separately because F_ST of the X is significantly higher than the autosomes (see below).

• Download figure
• Open in new tab
• Download PowerPoint

Consistent with our expectations for mitonuclear interactions, in all 10 pairwise comparisons, the average F_ST of autosomal N_interact windows (which contain the start site of at least one N_interact gene; see the electronic supplementary material) was significantly higher than autosomal windows that encode only other genes (p < 0.05 for each of 10 permutation tests). These findings are also apparent in density plots which illustrate that F_ST of N_interact windows consistently have more density over higher F_ST values as compared to genomic windows that contain only other genes, or compared to windows containing no genes (electronic supplementary material, figure S1), and also in boxplots of these data (electronic supplementary material, figure S2).

N_interact windows contained a higher number of genes on average (3.0 genes/window) compared to other windows that contained only non-N_interact genes (1.6 genes/window; hereafter non-N_interact windows). One possible explanation for the higher F_ST therefore is that N_interact windows typically have more genes than genomic windows which contain only other genes; this aspect of genome complexity was not considered in the permutation tests. To explore this possibility, for each of the 10 pairwise comparisons a linear regression was performed with an interaction term between two predictor variables, the number of genes in each 100 kb window (density) and whether (1) or not (0) the window was an N_interact window (N_interact), with F_ST as the response variable: F_ST ∼ N_interact * density. In all pairwise comparisons, density was significantly positively correlated with F_ST (p < 0.0001 for each test) but N_interact and the interaction terms were not significant (p > 0.05 for each test). This suggests that the higher F_ST of all N_interact windows is largely attributable to gene density rather than to the presence of N_interact genes (e.g. pink and grey dots in figure 2 overlap extensively), an observation that is not consistent with our expectations.

• Download figure
• Open in new tab
• Download PowerPoint

Although the set of 209 autosomal N_interact windows is not significantly more differentiated than other genomic windows with genes after controlling for gene density, it is still possible that natural selection against mitonuclear incompatibilities could cause a subset of N_interact windows to be extreme outliers for differentiation. To test this possibility, we performed an outlier analysis (electronic supplementary material) to identify windows where F_ST is much higher than expected after considering the positive relationship with density (upper F_ST outliers). For each of the 10 pairwise comparisons, approximately 3% of the non-N_interact windows were upper F_ST outliers. However, in these same 10 comparisons, the proportion of N_interact windows that were upper F_ST outliers was two to three times higher (range: 6–10%, figure 2; electronic supplementary material, table S3), a significant difference for all comparisons (p < 0.01, binomial tests). For each of the 10 pairwise comparisons, the proportions of lower F_ST outliers in non-N_interact windows were all approximately 1%. Unexpectedly, we found that the proportions of lower F_ST outliers were also higher in N_interact windows (range: 4–7%), which was also significant for all comparisons (p < 0.01, binomial tests). Overall, these observations are consistent with a subset of N_interact genes being atypically highly differentiated after controlling for gene density, even though F_ST of all N_interact windows was not significantly higher than other gene containing windows after controlling for gene density. But they also indicate that a subset of N_interact windows is atypically undifferentiated, an unexpected result that could be a consequence of purifying selection in regions with high gene density. Recognizing that these pairwise F_ST comparisons are not independent, it is nonetheless interesting that several N_interact windows were consistently upper F_ST outliers across multiple comparisons (electronic supplementary material, table S3 and Results).

Only two of the 211 N_interact genes are on the X chromosome (NDUFB11, COX7B) but F_ST of their N_interact windows was not significantly different from other X-linked genomic windows that carry only other genes (p > 0.05, permutation tests).

(b) A subset of N_interact genes is subject to atypically strong positive selection

As a complement to the analyses of F_ST, we examined other population-specific statistics that are sensitive to natural selection: pairwise nucleotide diversity (π), Tajima's D [22], and Fay and Wu's H [23]. These statistics are expected to have lower (and sometimes negative for D and H) values in positively selected genomic regions compared to other genomic regions; more details about these statistics, including possible influences of population size change, are provided in the electronic supplementary material. We found that after accounting for relationships with gene density, there was a significant excess of autosomal N_interact windows that were lower outliers for all three of these statistics for all five species considered (p < 0.01, binomial tests; electronic supplementary material, Results and tables S4, S5). In addition, polymorphism was atypically low in multiple species in the flanking regions of several N_interact windows—most prominently the ones containing MRPS21 and NDUFAF3 (electronic supplementary material, figure S3), and several of these lower outliers were also upper F_ST outliers (electronic supplementary material, table S3).

To explore the effects of the size of the genomic window on these results, we repeated the analysis of F_ST, π, Tajima's D, and Fay and Wu's H using a window size of 30 kb. The findings were essentially identical to those from the 100 kb windows (electronic supplementary material, Results).

If natural selection repeatedly favours compensatory mutations in N_interact genes that maintain compatibility with the rapidly evolving mitochondrial genome, we expected that these genes would be embedded in long runs of homozygosity (ROHs) that are generated by selective sweeps. Consistent with this expectation, the average lengths of ROHs containing N_interact genes (and possibly other genes) was 1.9–4.6 times longer than the average length of ROHs that contained only other genes, and 8.2–17.0 times longer than ROHs that did not contain any genes (electronic supplementary material, table S6). These differences are individually significant for all species (p ≤ 0.05, permutation tests; electronic supplementary material, table S6) except M. brunnescens (p = 0.35) and M. maura (p = 0.07). The atypically long ROHs that contain N_interact genes are clearly evident in density plots (electronic supplementary material, figure S4), and the longest ROH in most species contained at least one N_interact gene (electronic supplementary material).

We used linear models to account for the relationship between gene density and ROH length, and the results indicate that the relationship between N_interact genes and ROH length is nuanced for most species, with exceptions discussed in the electronic supplementary material. In general, at low to moderate gene densities, ROHs with N_interact genes tend to be much longer than those without. However, this relationship with ROH length varied among species as gene density increases (figure 3). For ROHs with high gene densities (e.g. greater than 20 genes/100 kb), those with N_interact genes tend to be similar or even shorter than those with only other genes. However, less than 1.2% of the ROHs had a genes density greater than 20 genes/100 kb and greater than 80% of the ROHs had a gene density below 5 genes/100 kb. Overall, this indicates a strong and pervasive positive correlation between N_interact genes and ROH length for the vast majority of observed ROHs after controlling for gene density.

• Download figure
• Open in new tab
• Download PowerPoint

In addition to influencing patterns of polymorphism and differentiation, natural selection on N_interact genes has the potential to influence protein evolution. To test this possibility, we analysed the rate ratio of nonsynonymous to synonymous substitutions per site (dN/dS) in N_interact genes and non-N_interact genes. The mean dN/dS of N_interact genes is significantly higher than the mean dN/dS of non-N_interact genes (p < 0.001, permutation test, electronic supplementary material), which is consistent with N_interact genes being under atypically relaxed purifying selection and/or atypically intense or frequent positive selection as compared to other genes.

(c) N_interact genes in two admixed individuals match expectations for selection on mitonuclear interactions

Using a hidden Markov model [24], we identified two M. tonkeana individuals—both collected near hybrid zones—that exhibited evidence of extensive introgression (figure 4). These individuals allowed us to test the prediction that selection against mitonuclear incompatibility could influence introgression. For individual PF626, which carries M. tonkeana mitochondrial DNA [14], we had a one-sided expectation that there would be fewer introgression blocks that contained N_interact genes from M. maura than expected by chance. For individual PF511, which carries introgressed mitochondria from M. hecki [14], we expected the opposite (i.e. an excess introgression blocks that carried N_interact genes from M. hecki).

• Download figure
• Open in new tab
• Download PowerPoint

In individual PF626, only two homozygous introgression blocks were identified that carried a total of eight genes out of 16 049 annotations; neither carried any N_interact genes. While this was consistent with our expectation, this small sample did not constitute a significant dearth of homozygous introgression blocks carrying N_interact genes (p = 0.902; permutation test). In individual PF511, there were 15 introgression blocks in total, and 12 of these contained a total of 54 genes out of the 16 049 annotated genes (0.3%). We therefore expected less than one of the 211 autosomal and X-linked N_interact genes would be on homozygous introgression blocks if there were no enrichment of N_interact genes. Instead, three N_interact genes were present on homozygous introgression blocks: NDUFB5, MRPL1 and MRPL47. Enrichment of N_interact genes on these introgression blocks in this individual was significant (p = 0.036; permutation test). Phylogenetic analyses provided further confirmation that the three N_interact genes in individual PF511 were indeed introgressed from M. hecki (electronic supplementary material).

(d) Introgression genomics point to sex-differences in demography in Macaca nemestrina

Population structure is higher on the X chromosome than the autosomes in these macaques (electronic supplementary material, figure S5). In each of these genomic compartments, we used two introgression statistics—Patterson's D [25] and f_dM [26]—to assess gene flow between pairwise combinations of one species (P₃) and either of two other species (P₁, P₂) whose phylogenetic relationships are (((P₁, P₂), P₃), O), where O is an outgroup taxon. A positive statistic arises from an excess of derived variants shared between P₂ and P₃ as compared to those shared between P₁ and P₃; negative values indicate the opposite. Non-zero values are potentially consistent with an inference of gene flow (e.g. for positive values, between P₂ and P₃), but can also arise from other scenarios that do not involve gene flow, such as ancestral population structure [27,28]. A rationale for the phylogenetic framework of these tests is provided in the electronic supplementary material, Results and figures S6–S8. Significant differentiation among these species is supported by admixture analysis [29], which favours a model with nine (autosomes) or eight (the X) populations (electronic supplementary material, figure S9).

Our analyses recover no evidence for introgression between M. nemestrina and the Sulawesi macaques after Sulawesi was colonized on the X or the autosomes. When P₁ and P₂ are defined to be any pair of Sulawesi macaques or groups of Sulawesi macaques that are reciprocally monophyletic, and P₃ is defined to be M. nemestrina from Borneo (parameterization 1), the estimated introgression statistics for the X and autosomes are not significantly different from zero (figure 5; electronic supplementary material, figure S10).

• Download figure
• Open in new tab
• Download PowerPoint

We considered an alternative phylogenetic parameterization for introgression statistics that also spans the Makassar Strait. When P₁ and P₂ are defined to be M. nemestrina from Sumatra and Borneo, respectively, and P₃ is defined to be any or all Sulawesi macaque species (parameterization 2), the introgression statistics for the X are consistently significantly positive, and those for the autosomes are either significantly negative, not significant or significantly positive depending on which individual or population is used as P₂ (figure 5; electronic supplementary material, figure S11). The differing results from parameterization 1 informed our interpretation of those from parameterization 2 and provided insights into how behaviour influenced population structure on the X and autosomes of M. nemestrina (see Discussion).

(e) Introgression genomics point to sex-differences in demography on Sulawesi

Across three hybrid zones on Sulawesi, introgression statistics are significantly positive for the autosomes and most analyses of the X chromosome (figure 6). These results thus support gene flow in the hybrid zones between M. nigrescens and M. hecki, M. hecki and M. tonkeana, and M. maura and M. tonkeana, with the highest level of gene flow between M. nigrescens and M. hecki. Across these three hybrid zones, Patterson's D was more extreme for the X chromosome than the autosomes in several analyses. f_dM mostly did not differ substantially between the X and autosomes, with one exception discussed in the electronic supplementary material. Analysis of individual genomes using a hidden Markov model [24] also recovered strong evidence of introgression between parapatric species of Sulawesi macaque across these three hybrid zones and identified substantial geographical heterogeneity in the extent of gene flow (figure 4; electronic supplementary material, figures S12–S15 and Results).

• Download figure
• Open in new tab
• Download PowerPoint

(a) Mitonuclear interactions in macaques

Using genomic data from 29 individuals from eight macaque species, we evaluated the expectation that deep divergences of mitochondrial genomes could lead to positive selection for compensatory mutations in N_interact genes, and negative selection against mitonuclear incompatibilities in admixed individuals. This expectation was supported by an excess of N_interact windows that were extreme upper F_ST outliers, lower π outliers and lower Tajima's D and Fay and Wu's H outliers (electronic supplementary material, tables S3–S5, and file S1), and a higher average rate of protein evolution of N_interact genes compared to non- N_interact genes (electronic supplementary material, Results). In addition, ROHs that contain N_interact genes were typically far longer than those that did not (figure 3; electronic supplementary material, figure S4 and table S6). In one M. tonkeana individual with introgressed mitochondrial DNA (PF511), there was a small but significant enrichment of introgressed N_interact genes from the source species of her mitochondrial genome (M. hecki).

Compensatory mutations in N_interact genes could mitigate the effects of Muller's ratchet in the mitochondria (i.e. in non-recombining genomes, the stepwise decrease in fitness over time associated with loss of the least deleterious genome [30]). Another nonexclusive possibility is that some N_interact genes are adaptively co-diverging with mitochondrial genomes [17]. Some N_interact genes (e.g. MRPS21, NDUFAF3) appear to have been subject to independent selective sweeps in multiple species, which suggests that there is variation among N_interact genes in the susceptibility of mitonuclear interactions to incompatibility mutations, in the degree to which interacting partners in the mitochondria and nucleus adaptively co-evolve, or both. This is also highlighted by our analyses that generally do not find significant signatures of positive selection for all N_interact genes, but consistently detect a substantial proportion of these genes that are outliers for multiple metrics of positive selection and/or population differentiation.

Of the four functional N_interact categories, three (ARS2, MRP and REP) interact with non-translated RNA genes or DNA (i.e. mitochondrial tRNA, ribosomal RNA and the control region), whereas OXPHOS genes interact with protein from translated mitochondrial genes. Thus, conflict resolution and co-adaptation across N_interact categories involve pairwise interactions between mitochondrial DNA, RNA and protein. In many cases, it is unclear exactly which portions of the nuclear and mitochondria-encoded genes have direct contact, e.g. [31]. For example, NDUFAF3 does not co-immunoprecipitate with the mitochondria-encoded ND1 protein even though both of these proteins are in the Q module of OXPHOS complex 1 [31], and NDUFAF3 is crucial for the incorporation and stability of ND1 in this complex [32]. Some aspects of compensatory evolution therefore may be driven by indirect interactions between mitochondria-encoded and autosomal-encoded factors. In other species, mitonuclear interactions have been maintained by compensatory mutations in N_interact genes, the recruitment of novel subunits, and gene duplication followed by subfunctionalization or neofunctionalization [21,33].

(b) Introgression genomics and sex-specific demography

Our analyses recovered no evidence of introgression between macaques on either side of the Makassar Strait in the autosomes or the X (parameterization 1, figure 5; electronic supplementary material, figure S10). This observation underscores the formidable biogeographical barrier represented by Wallace's Line within a widespread genus whose distribution straddles this line. Analysis with an alternative parameterization (parameterization 2) evidences geographically structured variation in M. nemestrina that was strongly influenced by sex-specific demography (figure 5). More specifically, a negative Patterson's D for autosomes when individual PM665 was set to be taxon P₂ suggests that there is an ancestry component that is shared between M. nemestrina from Sumatra and the Sulawesi macaques that is not present or less frequent in individual PM665, even though this individual originated from Borneo, which is geographically far closer to Sulawesi than Sumatra. However, positive introgression statistics are recovered for the autosomes when other Borneo individuals were set to be taxon P₂ (figure 5), which indicate that the ancestry component that is shared between the autosomes of Sulawesi and M. nemestrina from Sumatra is also present in some M. nemestrina on Borneo (more prominently in Borneo individuals PM664, GumGum and Sukai than in the Sumatra individual Ngsang, and at a similar level in Borneo individual PM1206 and Ngsang). This ancestry component was not detected with admixture analysis of individual genomes (electronic supplementary material, figure S9) and thus appears to be a small portion of these genomes.

This pattern in the autosomes contrasts sharply with a significantly more positive Patterson's D with parameterization 2 for the X chromosome (and parallel results for f_dM), which uniformly point to a prominent shared ancestry component between genomes from Borneo and Sulawesi macaques that is not present or rarer in the genome from Sumatra. Together these differing results in the autosomes and X suggest that genetic structure in the X is more strongly influenced than the autosomes by ancient population connections between Borneo and Sulawesi that reach back to more than 2 Ma to when Sulawesi was initially colonized by an ancestor of M. nemestrina from Borneo [14]. When data are concatenated [18] or weightings of alternative phylogenies are considered (electronic supplementary material, figures S6 and S8), there is strong support for monophyly in the autosomes and the X chromosome of M. nemestrina in Sumatra and Borneo with respect to the Sulawesi macaques. However, mitochondrial genomes of M. nemestrina from Borneo is more closely related to that of the Sulawesi macaques than they are to mitochondrial DNA from M. nemestrina from Sumatra [14], which mirrors the results from the X chromosome when analysed using parameterization 2 (figure 5). These differing relationships in the X and autosomes are consistent with male-mediated migration between Borneo and Sumatra M. nemestrina after macaques colonized Sulawesi coupled with female philopatry. In the past, overland migration was frequently possible between Borneo and Sumatra during times of lower sea level [34].

The large-X effect posits that the X chromosome plays an oversized role in species incompatibilities compared to the autosomes [35]. A prediction of the large-X effect is that there should be less gene flow of X-linked than autosomal-linked variation between species. This aspect of the large-X effect is similar to predictions associated with mitonuclear incompatibility in the sense that both may limit introgression of mostly (or entirely) female-inherited genomic compartments (the X and the mitochondria). Surprisingly because female philopatry might be expected to amplify the expectations of a large X effect by decreasing gene flow of this chromosome, our analyses do not identify the X chromosome as having substantially less introgression than the autosomes on Sulawesi. In fact, Patterson's D statistic (though not f_dM) across Sulawesi hybrid zones sometimes suggests higher rather than lower gene flow on the X than the autosomes across three hybrid zones on Sulawesi (figure 6). Perspective into these results on Sulawesi is gained by turning to the results from M. nemestrina, where parameterization 2 consistently recovered more extreme, positive introgression statistics on the X than the autosomes, even though no gene flow occurred between Borneo and Sulawesi after Sulawesi was colonized based on parameterization 1 (figure 5). We therefore suspect that the positive Patterson's D on the X is owing to ancient intraspecific population structure within the X chromosome that has persisted over evolutionary time scales (millions of years). Persistent ancestral geographical structure at the margins of species ranges thus may be conflated with estimates of introgression on the X chromosome, with a counterintuitive consequence of increasing the magnitude of Patterson's D on the X, even when gene flow is modest or absent. Because f_dM uses a dynamic estimator to tabulate patterns of genetic variation from the source (donor) population [26,28], this statistic appears to be less sensitive to the effects of ancient geographical structure on the X than Patterson's D (figure 6). Taken together, these results evidence demographic effects on Patterson's D and suggest this statistic should be evaluated in tandem with other metrics, such as divergence and genetic diversity of introgression blocks, and/or other introgression statistics such as f_dM [28]. Additional discussion of these findings in the context of biodiversity conservation is provided in the electronic supplementary material.

(a) Genetic samples and molecular data

Data analysed in this study include whole-genome data from 29 georeferenced monkeys (18 female and 11 male) from eight macaque species from Southeast Asia (electronic supplementary material, table S1; figure 1) that were sequenced to medium to high depth (range 11.0–40.0X; electronic supplementary material, figure S16). Three of these nuclear genomes were previously reported [18] and complete mitochondrial genomes of all of these samples were also previously reported [14]. For analyses that required an outgroup, we also included genomic data from a male P. anubis baboon (SRR1515161). These data were aligned to a reference genome from a rhesus macaque, genotyped and filtered as detailed in the electronic supplementary material.

(b) Nuclear-encoded proteins that interact with mitochondria-encoded proteins

We quantified summaries of differentiation, polymorphism, natural selection and introgression for genomic windows containing genes, with a focus on four categories of N_interact proteins that interact with mitochondria-encoded proteins [21], with additional detail about genes in each category provided in the electronic supplementary material. Gene acronyms of these N_interact genes are provided in the electronic supplementary material, file S1.

As detailed in the electronic supplementary material, we calculated F_ST, π, Tajima's D statistic [22], and Fay and Wu's H [23] in non-overlapping 30 and 100 kb genomic windows. Permutation tests, linear regressions and outlier analyses were used to evaluate whether these metrics were more extreme in N_interact windows as compared to non-N_interact windows, or windows that contained no genes. We estimated the rate ratio of non-synomymous to synonymous substitutions per site across a phylogeny of these samples and used permutation tests to evaluate whether this rate ratio was higher in N_interact genes than non-N_interact genes (electronic supplementary material). We additionally identified ROHs and used permutation tests and linear regressions to evaluate whether ROHs that contained N_interact genes tended to be longer than those that only contained non-N_interact genes, or those that contained no genes. All permutation tests were performed using scripts written in Perl; all linear models were fitted using the lm() function in R [40].

(c) Analysis of population structure and gene flow

To explore the effects of sex-biased migration within species in the X chromosome and autosomes, F_ST was calculated within each species for which more than one genome was sequenced. For M. nemestrina, F_ST was calculated between the Sumatra specimen and all samples from Borneo, and also between the two northern and three southeastern Borneo specimens. For several Sulawesi species (M. nigra, M. nigrescens, M. hecki, M. tonkeana and M. maura), specimens were divided into two populations based on their geographical distributions for F_ST calculations. Confidence intervals for F_ST were estimated using the weighted block jackknife approach [41].

Patterson's D statistic [25] (not to be confused with Tajima's D) and f_dM [26] were used to explore the effects of sex-biased migration between species in the X chromosome and autosomes. Confidence intervals were estimated using the weighted block jackknife approach [41] with non-overlapping 100 000 bp genomic regions that were weighted by the sum of the numbers of variable positions that were used in each calculation in each window. Additionally, introgression blocks were inferred from individual genomes using a hidden Markov model [24] as detailed in the electronic supplementary material.