We sampled 11 males of P. syfanius and P. maackii across a geographic transect (Figure 2A, dashed lines) covering both pure populations (in the sense of being geographically far from the hybrid zone) and hybrid populations. We also include four outgroup species, two of which have chromosome-level genome assemblies (P. bianor, P. xuthus; Lu et al., 2019; Li et al., 2015), while the other two (P. arcturus, P. dialis) are new to this study. All samples were re-sequenced to at least 20× coverage across the genome and mapped to the genome assembly of P. bianor. Among sampled local populations, P. syfanius inhabits the highlands of Southwest China (Figure 2A, red region), whereas P. maackii dominates at lower elevations (Figure 2A, blue region) (see Figure 2—figure supplement 3 and Figure 2—source data 4 for the complete distribution). The two lineages form a spatially contiguous hybrid zone at the edge of the Hengduan Mountains (Figure 2B) with individuals exhibiting intermediate wing patterns (Figure 2A: purple dot, corresponding to population WN in Figure 2B). Consistent with previous results (Condamine et al., 2013), assembled whole mitochondrial genomes are not distinct between the two lineages (Figure 2C, Figure 2—source data 3), suggesting either that divergence was recent, or that gene flow has homogenized the mitochondrial genomes. However, the two species are likely adapted to different environments associated with altitude, as several ecological traits are strongly divergent (Kashiwabara, 1991; Figure 2—figure supplement 1). Similarly, between pure populations (KM and XY in Figure 2B), relative divergence ($F_{ST}$) is also high across the entire nuclear genome (Figure 2D, Figure 2—source data 1). The $F_{ST}$ on autosomes averages between 0.2–0.4, and on the sex chromosome (Z-chromosome) it reaches 0.78. A highly heterogeneous landscape of $F_{ST}$ is accompanied by numerous islands of elevated sequence divergence ($D_{XY}$) and reduced genetic diversity ($\pi$) scattered across the genome (Figure 2—figure supplement 4, Figure 2—source data 2). Since animal mitochondrial DNA generally has higher mutation rates than the nuclear genome (Haag-Liautard et al., 2008), its low divergence between the two species are likely the result of gene flow. Overall, despite ongoing hybridization, genomes of the pure populations of P. syfanius and P. maackii are strongly differentiated.

Figure 2 with 4 supplements see all

Download asset Open asset

Figure 2—source data 1

Estimated $F_{ST}$ between populations KM and XY (50 kb windows with 10 kb increments).

https://cdn.elifesciences.org/articles/78135/elife-78135-fig2-data1-v2.zip

Download elife-78135-fig2-data1-v2.zip

Figure 2—source data 2

$D_{XY}$ and $\pi$ estimated for populations KM and XY (non-overlapping 10 kb windows).

https://cdn.elifesciences.org/articles/78135/elife-78135-fig2-data2-v2.zip

Download elife-78135-fig2-data2-v2.zip

Figure 2—source data 3

The phylogeny and alignment among mitochondrial genomes.

https://cdn.elifesciences.org/articles/78135/elife-78135-fig2-data3-v2.zip

Download elife-78135-fig2-data3-v2.zip

Figure 2—source data 4

Results and data for MaxEnt species distribution models.

https://cdn.elifesciences.org/articles/78135/elife-78135-fig2-data4-v2.zip

Download elife-78135-fig2-data4-v2.zip

Figure 2—source data 5

Sample information.

https://cdn.elifesciences.org/articles/78135/elife-78135-fig2-data5-v2.zip

Download elife-78135-fig2-data5-v2.zip

A natural question is whether genomic differentiation is associated with barrier loci and reproductive isolation. In other words, can $F_{ST}$ variation be attributed to gene flow variation between sister species? We suspect that barrier loci likely exist, because sequence variation between pure populations suggests that elevated $F_{ST}$ is associated with reduced $\pi$ and elevated $D_{XY}$ across autosomes (Figure 3A), as expected for hybridizing species (Irwin et al., 2018). The Z chromosome (sex chromosome) has the highest $D_{XY}$ and the lowest $\pi$ among all chromosomes (Figure 2—figure supplement 4), another characteristic of hybridizing species with barriers to gene flow (Kronforst et al., 2013).

Figure 3 with 1 supplement see all

Download asset Open asset

Figure 3—source data 1

Estimated entropy ($S_w$, $S_b$), $D_{XY}$, $\pi$, $F_{ST}$ on 20 segments per chromosome (including the Z chromosome).

https://cdn.elifesciences.org/articles/78135/elife-78135-fig3-data1-v2.zip

Download elife-78135-fig3-data1-v2.zip

To test for the presence of barrier loci, we augment the analysis with the sequences of four individuals from the population closest to the center of the hybrid zone (Population WN). We investigate whether differences in ancestry variation provide additional evidence for barrier loci in this hybrid zone. The underlying logic is that barrier loci will simultaneously:

1. Reduce linked $\pi$ in pure populations (Ravinet et al., 2017);

2. Elevate linked $D_{XY}$ between pure populations (Ravinet et al., 2017);

3. Elevate pairwise linkage disequilibrium in hybrid zones (Barton, 1983);

4. Enrich linked ancestry from one lineage in hybrid zones (Sedghifar et al., 2016).

Effects 3 and 4 can be bundled together as “reduced ancestry randomness” around barrier loci because both are expected if intermixing of segments of different ancestries within a genomic interval is prevented. For effects 3 and 4, because of small sample sizes, estimating site-specific statistics such as pairwise linkage disequilibrium is untenable. However, our high-quality chromosome-level reference genome enabled accurate estimation of local ancestry. As a remedy for small sample sizes, we employ two entropy metrics borrowed from signal-processing theory to quantify ancestry randomness in local regions along chromosomes (see Materials and methods). By dividing chromosomes into segments, we can extract indirect information about effects 3 and 4 at the expense of reduced genomic resolution. The proposed metrics, $S_b$ and $S_w$, correspond to the randomness of ancestry between and within individual diploid genomes from a local population (Figure 3B). For a cohort of ideal chromosomes with uniform recombination and marker density, if ancestry is independent between homologous chromosomes, both $S_b$ and $S_w$ increase with reduced local ancestry correlation and more balanced parental contribution (Figure 3C).

For a given autosomal segment, we then calculate $\pi$, $D_{XY}$, and $F_{ST}$ between pure populations, as well as entropy metrics $S_b$ and $S_w$ in hybrid individuals for the same segment. To investigate whether effects 1–4 are all present in our system, the joint range among entropy, $\pi$, and $D_{XY}$ is shown in Figure 3D, which suggests that low ancestry randomness (low entropy) is likely associated with reduced $\pi$ within species and elevated $D_{XY}$ between species. To further quantify such association, we estimated Pearson’s correlation coefficients ($\rho$) between entropy and the latter statistics (Figure 3E). These associations are strongly significant ($Z$-scores > 3). Consequently, reduced ancestry randomness in hybrids (effects 3 & 4) coincides with classical patterns of barrier loci between pure populations (effects 1 and 2). This analysis is not sufficient to exclude all alternative hypotheses. For instance, we cannot entirely exclude the possibility that patterns are driven by low-recombination regions ($S_w,S_b\downarrow$) that experience linked selection ($\pi \downarrow$) also have elevated mutation rates ($D_{XY}\uparrow$). Nonetheless, this alternative seems most unlikely, as low-recombination regions should typically be less rather than more mutable (Lercher and Hurst, 2002; Jensen-Seaman et al., 2004; Yang et al., 2015; Arbeithuber et al., 2015; Liu et al., 2017). Overall, adding information from hybrid populations strengthens the evidence for barrier loci acting across autosomes.

The Z chromosome was excluded from this analysis as it likely differs in mutation rate or effective population size (Presgraves, 2018), but its ancestry in hybrid individuals either retains purity or resembles very recent hybridization (long blocks of heterozygous ancestry, Figure 3—figure supplement 1). The Z chromosome has low ancestry randomness, and it also has the highest level of divergence (Figure 2D), both of which suggest strong barriers to gene flow between P. syfanius and P. maackii on this chromosome.

A hint that substitution rates differ between the two species comes from site-pattern asymmetry (but this asymmetry alone is insufficient to establish the existence of divergent rates). We focus on two kinds of biallelic site patterns. In the first kind, choose three taxa (P₁,P₂,O₁), with P₁=syfanius, P₂=maackii, while O₁ is an outgroup. Assuming no other factors, if substitution rates are equal between P₁ and P₂, then site patterns (P₁,P₂,O₁)=(A,B,B) and (B,A,B) occur with equal frequencies, where “A” and “B” represent distinct alleles. This leads to a $D$ statistic describing the asymmetry between three-taxon site patterns:

where $n$ is the count of a particular site pattern across designated genomic regions. A significantly non-zero ${\displaystyle D_3}$ indicates strongly asymmetric distributions of alleles between P₁ and P₂.

In a second kind of site pattern test, we compare four taxa (P₁,P₂,O₁,O₂) and calculate a similar statistic between site patterns (A,B,B,A) vs (B,A,B,A):

Classically, a significantly non-zero ${\displaystyle D_4}$ suggests that gene flow occurs between an outgroup and either P₁ or P₂, thus it is widely used to detect hybridization (the ABBA-BABA test) (Durand et al., 2011; Hibbins and Hahn, 2022). Here, ${\displaystyle D_4}$ is used more generally as an additional metric of site pattern asymmetry.

We compute both ${\displaystyle D_3}$ and ${\displaystyle D_4}$ on synonymous, nonsynonymous, and intronic sites, with P. syfanius and P. maackii samples taken from pure populations (KM and XY). For each type of site, we progressively exclude regions with local $F_{ST}$ below a certain threshold, and report $D$ statistics on the remaining sites in order to show the increasing site-pattern asymmetry in more divergent regions (Figure 4). ${\displaystyle D_3}$ is significantly negative regardless of outgroup or site type for most $F_{ST}$ thresholds, and ${\displaystyle D_4}$ is also significantly negative for most outgroup combinations when computed across the entire genome ($Z$-scores are shown in Figure 4—figure supplement 1), proving that site-patterns are strongly asymmetric between P. syfanius and P. maackii. Importantly, the direction of asymmetry is nearly identical across all outgroup comparisons. This asymmetry cannot be attributed to batch-specific variation as all samples were processed and sequenced in a single run, and variants were always called on all individuals of P. syfanius and P. maackii. Sequencing coverage is normal for most annotated genes used in the analysis (Table 1), suggesting that asymmetry is not due to systematic copy-number variation that could affect variant calls. Nonetheless, two independent processes of evolution could explain observed asymmetric site patterns. In hypothesis I (Figure 5A, left), asymmetry is generated via stronger gene flow between P. syfanius and outgroups, leading to biased allele-sharing. In hypothesis II, site pattern asymmetry is due to unequal substitution rates between P. syfanius and P. maackii, which is further modified by recurrent mutations in all four outgroups (Figure 5A, right). We test each hypothesis below.

Figure 4 with 1 supplement see all

Download asset Open asset

Figure 4—source data 1

, ${\displaystyle D_4}$ statistics with their $Z$-scores.

https://cdn.elifesciences.org/articles/78135/elife-78135-fig4-data1-v2.zip

Download elife-78135-fig4-data1-v2.zip

Figure 5

Download asset Open asset

Figure 5—source data 1

Concatenated gene trees based on 50 kb windows and the underlying support for each binary split.

https://cdn.elifesciences.org/articles/78135/elife-78135-fig5-data1-v2.zip

Download elife-78135-fig5-data1-v2.zip

Figure 5—source data 2

Results of the signed-rank test on branch elongation in P. maackii.

https://cdn.elifesciences.org/articles/78135/elife-78135-fig5-data2-v2.zip

Download elife-78135-fig5-data2-v2.zip

We first test hypothesis I by phylogenetic reconstruction using SNPs in annotated regions. We construct local gene trees for each 50 kb non-overlapping window for all samples, including the four outgroup species. As biased gene flow with outgroups should rupture the monophyletic relationship among all P. syfanius +P. maackii individuals, the fraction of windows producing paraphyletic gene trees can be used to assess the potential impact of gene flow. However, almost all gene trees show the expected monophyletic relationship (Figure 5B). This conclusion is independent of the level of support used to filter out genomic windows with ambiguous topologies (see Materials and methods). Consequently, hardly any window shows a phylogenetic signal of hybridization with outgroups. Nonetheless, branch lengths in reconstructed gene trees suggest possible substitution rate divergence between the two species: Among highly supported monophyletic trees (bootstrap support ≥ 95), P. maackii (in populations YA, BJ, XY) is always significantly more distant than P. syfanius from the most recent common ancestor of P. maackii +P. syfanius (Figure 5C, significance levels reported via a Wilcoxon signed-rank test for each pair of individuals, see Figure 5—source data 2). Second, the direction of allele sharing in the $D$ statistics is unanimously biased towards P. syfanius. If hypothesis I were true, hybridization with outgroups is required to take place mainly in the highland lineage. There is no evidence to support why the highland lineage should receive more gene flow based on current geographic distributions, as outgroups P. xuthus and P. bianor overlap broadly with both P. maackii and P. syfanius, while outgroup P. dialis is sympatric only with P. maackii (Condamine et al., 2013). Although it is possible that historical and modern geographic distributions differ, and archaic gene flow might have occurred preferentially with P. syfanius, it should still leave some phylogenetic signal of introgression. Overall, we find little evidence for biased hybridization required by hypothesis I.

One might worry that by rejecting hypothesis I, we also throw doubt on widely accepted conclusions of the ABBA-BABA test for gene flow in other systems (that a significantly nonzero ${\displaystyle D_4}$ implies hybridization with outgroups) (Durand et al., 2011). However, in the next section, we show why ${\displaystyle D_3}$ and ${\displaystyle D_4}$ are fully consistent with hypothesis II, and so this is just a special case where the ABBA-BABA test produces a false positive for gene flow.

In hypothesis II, divergent substitution rates between P. maackii and P. syfanius interact with recurrent mutations in outgroups to produce asymmetric site patterns. To understand its effect on $D$ statistics, consider a simplified model of recurrent mutation (Figure 5D,E), where a site in outgroup $i$ mutates with probability ${\displaystyle p_i}$, producing the same derived allele with probability $c$. When averaged across the genome, $c$ can be treated as a constant, and ${\displaystyle p_i}$ increases with distance to the outgroup. In the absence of gene flow, for three-taxon patterns, recurrent mutations modify ${\displaystyle D_3}$ by a factor of approximately $(1-2c{p}_1)$ (Figure 5D, see Materials and methods), and observed ${\displaystyle D_3}$ will thus be positively correlated with ${\displaystyle p_1}$. For four-taxon patterns, it can be shown that observed ${\displaystyle D_4}$ is always negative due to larger probabilities of recurrent mutation in more distant outgroups (Figure 5E, see Materials and methods). Assuming no significant contribution of incomplete lineage sorting (Maddison and Knowles, 2006), the value of ${\displaystyle D_4}$ becomes more negative with increasing ${\displaystyle \mathrm{\Delta }{p}_i/\mathrm{\Sigma }{p}_i=(p_2-p_1)/(p_2+p_1)}$.

To test these signatures, we employ estimated node heights of outgroups in the mitochondrial tree (Figure 2C) as proxies for outgroup distance, and hence for the relative probability of recurrent mutation (${\displaystyle {\rho }_i}$). In line with expected signatures, we find that observed ${\displaystyle D_3}$ indeed increases with node height (Figure 5F, left), and observed ${\displaystyle D_4}$ decreases with (Δ Node height/∑ Node height; Figure 5F, right). Thus, the directions and magnitudes of both $D$ statistics are congruent with hypothesis II. As hypothesis II naturally predicts unanimously negative ${\displaystyle D_3}$ and ${\displaystyle D_4}$ as well as their relative magnitudes among different outgroup combinations, it is more parsimonious than hypothesis I. Hence, divergent substitution rates likely exist between P. syfanius and P. maackii.

Having tested for the existence of divergent substitution rates, we now explore how they become mixed by gene flow using a coalescent framework. In particular, we will assess whether the conceptual picture in Figure 1 can be recovered quantitatively in the butterfly system. As gene flow is ongoing between the two lineages, consider two haploid populations of size $N$ exchanging genes at rate $m$ (Figure 6A). This simple isolation-with-migration (IM) model at equilibrium has relative divergence (Notohara, 1990)

Figure 6 with 3 supplements see all

Download asset Open asset

To quantify the signature of mixed rates, consider the asymmetry in observed numbers of substitutions (circles in Figure 6A). Take a pair of sequences from two populations. Let $⟨n_1⟩$ be the expected number of derived alleles exclusive to the sequence in population 1, and let $⟨n_2⟩$ be the same expected number in population 2. Their ratio is defined as observed rate ratio:

Evidently, $r=1$ is the symmetric point where both sequences have the same number of derived alleles. Further, let the actual substitution rate in population 1 be ${\mu }_1$, and the actual substitution rate in population 2 be ${\mu }_2$. The ratio between the two actual rates are defined as the true rate ratio:

At migration-drift equilibrium, observed rate ratio ($r$) and observed divergence ($F_{ST}$) are related by the following formula parameterized singly by ${\displaystyle r_0}$ (Figure 6B, see Materials and methods):

which translates into

These formulae indicate that unequal substitution rates are more mixed in regions with lower genomic divergence. Equation 7 is surprising, because it reveals that the remaining fraction of substitution rate divergence ($(r-1)/(r_0-1)$) is almost the same as $F_{ST}$, which corresponds to the fraction of genetic variance explained by population structure (Wright, 1949). In Figure 6B, the relationship between $r$ and $F_{ST}$ is still largely linear when one species evolves three times as fast as its sister species ($r_0=3$), and so Equation 7 might be robust under biologically realistic rates of substitutions between incipient species. Using extensive simulations (Figure 6—figure supplement 1 to Figure 6—figure supplement 3), we show that the full formula (Equation 6) is also robust in a number of equilibrium population structures.

To test whether such predictions are met, we calculate $r$ on synonymous, nonsynonymous, and intronic sites partitioned by their local $F_{ST}$ values between pure populations (KM & XY), and recover a similar monotonic relationship across most $F_{ST}$ partitions (Figure 6C). For introns, the observed relationship between $r$ and $F_{ST}$ is a near perfect fit to Equation 6 (Figure 6D, squares), with an estimated ${\displaystyle r_0}$ = 1.837. For synonymous (Figure 6D, circles) and nonsynonymous sites (Figure 6D, triangles), Equation 6 also provides an excellent fit in regions with low to intermediate $F_{ST}$. Estimated ${\displaystyle r_0}$ for synonymous sites is 1.813, close to that of introns, while nonsynonymous sites have a considerably lower ${\displaystyle r_0}$ of 1.633. If introns and synonymous substitutions are approximately neutral, we infer that neutral substitution rates are about 80% greater in the lowland species.

A critical step in our analysis is to associate genomic islands with barrier loci. Conventionally, information on increased association between barrier sites in hybrid populations comes from two-locus linkage disequilibrium (Slatkin, 1975). Empirical studies on hybrid populations frequently use such statistics to strengthen the evidence for barriers to gene flow (Knief et al., 2019; Wang et al., 2022). Alternatively, if phased haplotypes can be sequenced, the length of ancestry tracts (Sedghifar et al., 2016) or the density of ancestry junctions (Janzen et al., 2018; Wang et al., 2022) also carry information on barrier loci. All these methods break down with small numbers of unphased samples, as forced phasing can produce a large number of false ancestry junctions. However, conventional notions such as “ancestry tract length” are not definable for unphased local ancestry. The dilemma forces us to consider more robust statistics that carry information on ancestry association even in unphased data.

The development of both entropy metrics follows three intuitions:

• Vectorization: ancestry is a categorical concept, and it should map to a signal containing the contribution of each source population.

• Conservation law: the total ancestry from all sources at a particular locus is conserved.

• Highly (auto)correlated signals have a concentrated spectrum (low spectral entropy).

Since entropy carries information about the degree of ancestry correlation, it will decrease in regions of low recombination and high genetic relatedness. If genetic relatedness is produced by inbreeding, it should affect the entire genome in a similar way, and between-individual entropy ($S_b$) will be similar across different parts of the genome. However, if inbreeding is so severe that $S_b$ is globally zero, it will not be an informative metric. It is worth noting that within-individual entropy $S_w$ shares a conceptual similarity to the wavelet transform of ancestry signals (Pugach et al., 2011; Sanderson et al., 2015). These entropy metrics will be particularly useful for window-based genomic analysis of ancestry correlation with limited sampling, and they are also compatible with larger sample sizes and more than two parental species. (The formulae of entropy with two parental species are presented in Materials and methods, and mathematical details are discussed in Appendix 1.) Nonetheless, this entropy approach assumes local ancestry can be accurately inferred, which will be challenging for studies with low-coverage sequencing, non-chromosomal genome assembly, or lacking reference populations.

Interestingly, a higher substitution rate in the lowland lineage P. maackii is congruent with the evolutionary speed hypothesis (Rensch, 1959), where evolution accelerates in warmer climates. Our finding echoes the results of many previous studies (Gillooly et al., 2005; Wright et al., 2006; Lin et al., 2019; Ivan et al., 2022). Without measuring the per-generation mutation rates in both species, it is unclear what mechanisms cause increased substitution rates, but the lowland lineage typically has an additional autumn brood that is absent in P. syfanius (Takasaki et al., 2007; Figure 2—figure supplement 2). Warmer temperatures in lowland habitats might also increase spontaneous mutation rates in ectothermic insects (Waldvogel and Pfenninger, 2021). Both mechanisms could produce increased substitution rates in the lowland species. It is surprising that estimated ${\displaystyle r_0}$ between synonymous sites and introns agree with each other under such a coarse framework, and estimated ${\displaystyle r_0}$ for nonsynonymous sites is considerably lower. Less asymmetric rates for nonsynonymous substitutions could perhaps be explained by the nearly neutral theory (Ohta, 1993), which argues that many nonsynonymous mutations are mildly deleterious, and selection is more efficient in suppressing them in larger populations and slowing down substitutions. In the field, the lowland species often appears in larger numbers with well-connected habitats, while the highland species faces a highly heterogeneous landscape of the Hengduan Mountains, which could lead to differences in effective population sizes required by our expectation under the nearly neutral theory.

An additional result from our study is that divergent substitution rates might produce spuriously non-zero ${\displaystyle D_4}$ statistics when combined with recurrent mutations, which could increase the false positive rate of the ABBA-BABA test (Durand et al., 2011). This phenomenon has been suspected in humans (Amos, 2020), and is certainly a theoretical possibility (Hibbins and Hahn, 2022), but has not been tested in most empirical studies.

A wide class of introgression tests targeting gene flow between outgroups and a pair of taxa are based on site pattern information. Hahn’s ${\displaystyle D_3}$ computes absolute sequence divergence between groups in a triplet of species (Hahn and Hibbins, 2019), and will be affected by unequal substitution rates in a similar way to our ${\displaystyle D_3}$. Martin’s ${\displaystyle f_d}$ computes the same numerator as our ${\displaystyle D_4}$ (Martin et al., 2015), so it could also produce false positives under similar situations. Related statistics include $D_p$ (Hamlin et al., 2020), $D_f$ (Pfeifer and Kapan, 2019). A general guideline for site-pattern based statistics is that the focal pair of taxa are closely related such that substitution rates do not differ, while outgroups should not be too distant to minimize recurrent mutations (Hibbins and Hahn, 2022). However, whether these assumptions are met in empirical studies is worthy of investigation, and our system provides a counterexample even between sister species with ongoing hybridization.

A separate pitfall might occur if introgression tests based on site-patterns are applied to genomic windows to locate regions introgressed from outgroups. In our case, $F_{ST}$ peaks have the most asymmetric substitution rates between P. maackii and P. syfanius, thus they will most likely be associated with false-positive $D$ statistics. This could lead to the incorrect interpretation that some barrier loci (“speciation genes”) are introgressed from outgroups – a popular hypothesis in adaptive introgression and hybrid speciation, see Edelman and Mallet, 2021.

To this end, we speculate that using appropriate substitution models to infer gene tree topology will perform better in assessing the impact of introgression with outgroups.

In the gray zone of incomplete speciation, interspecific hybrids bridge between gene pools of divergent lineages (Mallet, 2005). We here demonstrate a similar role of hybridization in coupling and mixing differing substitution rates. Divergent rates of substitution carry information about outgroups, while divergence based on allele frequency differences does not. Preserving divergent substitution rates is a stronger effect than maintenance of allele frequency differences, because divergence of allele frequencies is a prerequisite for rate preservation. This dependency can be coarsely quantified across the genome by the relationship between observed rate ratio ${\displaystyle r}$ and relative divergence $F_{ST}$ in an equilibrium system of hybridizing populations (Equation 6). At migration-drift equilibrium, it is not surprising that divergent substitution rates are associated with relative divergence. In Figure 6A, when coalescence occurs rapidly compared to gene flow, most substitutions separating individuals are species-specific. However, when gene flow is faster than coalescence, individuals will carry substitutions that occurred in both species. This could have important implications, because preserving lineage-specific substitution rates as measured by $r$ might not require low absolute rates of gene flow. Instead, reducing effective population sizes via recurrent linked selection might achieve a similar result in populations at equilibrium ($N\downarrow \Rightarrow F_{ST}\uparrow \Rightarrow r\uparrow$).

The theory in its present form has several limitations. First, mutations follow the infinite-site model (Kimura, 1969), so that reverse mutations, double mutations, and substitution types are not accurately reflected in the estimated ratio between species-specific substitution rates. Second, population structure is assumed at equilibrium, whereas real data could carry footprints from non-equilibrium demographic processes (e.g. secondary contact) (Hey, 2010). Third, there could be considerable population structure within each species contributing to elevated $F_{ST}$ but not necessarily to rate-divergence. This effect could be seen in simulated stepping-stone models (Figure 6—figure supplement 3) and will underestimate the level of rate divergence. Fourth, substitution is stochastic across the genome, and accurately estimating observed rate ratio $r$ relies on averaging substitution numbers across a large number of sites. This poses a problem if sites linked to high $F_{ST}$ regions are rare (i.e. few genomic islands). Lastly, as the theory is built upon the neutral coalescent, it is best suited for studying behaviors of neutral sites.

Nonetheless, the monotonic relationship between $r$ and $F_{ST}$ (i.e. larger sequence divergence is associated with more asymmetric substitution rates) might be qualitatively robust regardless of the aforementioned caveats. For instance, even for hybrid zones formed by recent secondary contact, reducing the absolute rate of gene flow by barrier loci in principle also keeps divergent rates from mixing, simply because it prevents substitutions accumulated in the allopatric phase from flowing between species.

In conclusion, our study characterizes several genomic consequences of the rate-mixing process when molecular clocks tick at different speeds between hybridizing lineages. This process provides new information on reproductive isolation but also leads to pitfalls in interpreting results of popular introgression tests. As this phenomenon is neglected in most studies of hybridization and speciation, its full scope awaits further investigation in both theories and empirical systems.

Museum specimens with verifiable locality data of all species were gathered from The University Museum of The University of Tokyo (Harada et al., 2012; Yago et al., 2021), Global Biodiversity Information Facility (The Global Biodiversity Information Facility, 2021b; The Global Biodiversity Information Facility, 2021a), and individual collectors (Figure 2—source data 5). Records of P. maackii from Japan, Korea and NE China were excluded from the analysis, so that most P. maackii individuals correspond to ssp. han, the subspecies that hybridizes with P. syfanius. Spatial principal component analysis was performed on elevation, maximum temperature of warmest month, minimum temperature of coldest month, and annual precipitation, all with 30s resolution from WorldClim-2 (Fick and Hijmans, 2017). The first two PCAs, combined with tree cover (Hansen et al., 2013), were used in MaxEnt-3.4.1 to produce species distribution models that use known localities to predict occurrence probabilities across the entire landscape (Phillips et al., 2017). Outputs were trimmed near known boundaries of each species. See Figure 2—figure supplement 3 for the final result.

Eleven males of P. syfanius and P. maackii, with one male of P. arcturus and one male of P. dialis were collected in the field between July and August in 2018 (Figure 2—source data 5), and were stored in RNAlater at –20 °C prior to DNA extraction. E.Z.N.A Tissue DNA kit was used to extract genomic DNA, and KAPA DNA HyperPlus 1/4 was used for library preparation, with an insert size of 350 bp and 2 PCR cycles. The library is sequenced on a Illumina NovaSeq machine with paired-end reads of 150 bp. Adaptors were trimmed using Cutadapt-1.8.1, and subsequently the reads were mapped to the reference genome of P. bianor with BWA-0.7.15, then deduplicated and sorted via PicardTools-2.9.0. The average coverage among 13 individuals in non-repetitive regions varies between 20× and 30×. Variants were called twice using BCFtools-1.9 – the first including all samples, used in analyses involving outgroups, and the second excluding P. arcturus and P. dialis, used in all other analyses. The following thresholds were used to filter variants: ${\displaystyle 10N<}$ DP ${\displaystyle <50N}$, where $N$ is the sample size; QUAL > 30; MQ > 40; MQ0F < 0.2. As a comparison, we also called variants with GATK4 and followed its best practices, and 93% of post-filtered SNPs called by GATK4 overlapped with those called by BCFtools. We used SNPs called by BCFtools throughout the analysis. Mitochondrial genomes were assembled from trimmed reads with NOVOPlasty-4.3.1 (Dierckxsens et al., 2017), using a published mitochondrial ND5 gene sequence of P. maackii as a bait (NCBI accession number: AB239823.1). We also used the following published mitochondrial genomes (NCBI accession numbers): KR822739.1 (Papilio glaucus), NC_029244.1 (Papilio xuthus), JN019809.1 (Papilio bianor). The neighbor-joining mitochondrial phylogeny was built with Geneious Prime-2021.2.2 (genetic distance model: Tamura-Nei), and we used 10⁴ replicates for bootstrapping. The reference genome of P. xuthus was previously aligned to the genome of P. bianor and we used this alignment directly in all analysis (Lu et al., 2019).

Given a species tree {{P₁,P₂},O}, where P₁ and P₂ are sister species and O is an outgroup, if mutation rates are equal between {P₁,P₂}, and no gene flow with O, then on average the number of derived alleles in P₁ should equal the number of derived alleles in P₂. Let ${\displaystyle \mathcal{S}}$ be a collection of sites, ${\displaystyle f_s}$ be the frequency of a particular site pattern at site ${\displaystyle s\in \mathcal{S}}$. “ABB” be the pattern where only P₂ and O share the same allele, and “BAB” be the pattern where only P₁ and O share the same allele, then the three-species ${\displaystyle D_3}$ statistic is calculated as

where ${\displaystyle \mathcal{S}}$ is always limited to sites without polymorphism in the outgroup O. This statistic is in principle capturing the same source of asymmetry as the statistic proposed by Hahn and Hibbins, 2019, although their version uses divergence to the outgroup instead of frequencies of site-counts. Similarly, the four-species ${\displaystyle D_4}$ statistic, which considers species tree {{{P₁,P₂},O₁},O₂} and site patterns ABBA versus BABA (Durand et al., 2011) is calculated as

where ${\displaystyle \mathcal{S}}$ is always limited to sites without polymorphism in the second outgroup O₂. The significance of both tests was computed using block-jackknife over 1 Mb blocks across the genome. Additionally, we estimated rate ratio as follows. First we restricted to sites where all outgroups are fixed for the same ancestral allele to dampen the influence of recurrent mutation. Then, for each site, sample one allele at random from each focal lineage. Calculate the probability of observing a derived allele in P₁ but not in P₂, and the probability of observing a derived allele in P₂ but not in P₁. The rate ratio is computed as the ratio between the two probabilities. Explicitly, let $I(\cdot )$ be the identity function, and ${\displaystyle f_s}$ be the frequency of the derived allele, then:

Its standard error was estimated using 1 Mb block-jackknifing. We excluded P. xuthus from the outgroups to increase the number of informative sites when using this formula.

In this section, we calculate observed ${\displaystyle D_3}$ and ${\displaystyle D_4}$ assuming that incomplete lineage sorting contributes insignificantly to both statistics. If incomplete lineage sorting is present, it will not create new bias (numerators are on average unchanged), but will likely dampen existing bias (inflating denominators).

As substitutions are independent along each lineage, we can mute recurrent mutations in outgroups and generate them afterwards. For three taxa with gene tree {{P₁,P₂},O₁}, before recurrent mutations, there are ${\displaystyle n_1}$ sites with pattern (B,A,A), and ${\displaystyle n_2}$ sites with pattern (A,B,A). If substitution rate is higher in P₂, we have ${\displaystyle n_2>n_1}$, so the true value of ${\displaystyle D_3}$, written as ${\widehat{D}}_3$, is always negative:

Next, recurrent mutations in O₁ occur at each site with an average probability ${\displaystyle p_1}$, and with an average probability $c$, ancestral alleles from affected sites in O₁ are converted to the same derived allele in {P₁,P₂}. $c$ will be independent of ${\displaystyle n_1}$ and ${\displaystyle n_2}$, as long as substitutions between {P₁,P₂} are only different in rates, but not mutation types. Hence, two possible mutation paths exist:

The expected site counts, after recurrent mutations, become

Using the new expected site counts in ${\displaystyle D_3}$ statistics produce the following value:

Since ${\widehat{D}}_3$ is negative, it grows approximately linearly near small values of ${\displaystyle p_1}$. (The full equation is still monotonic in ${\displaystyle p_1}$.)

Similarly, for four-taxon statistics, before recurrent mutation, there are two types of sites: (A,B,B,B)- ${\displaystyle n_1}$; (B,A,B,B)- ${\displaystyle n_2}$. Suppose the average probability of recurrent mutation is ${\displaystyle p_1}$ in O₁, and ${\displaystyle p_2}$ in O₂, and the conversion probability of each recurrent mutation into derived alleles of {P₁, P₂} is $c$ for both outgroups. Using the same procedure, one can show that

Since recurrent mutations occur more frequently in distant outgroups, ${\displaystyle p_2>p_1}$. Because ${\widehat{D}}_3$ is negative, we have ${\displaystyle D_4<0}$.

Local gene trees were estimated using iqtree-2.0 (Minh et al., 2020) on 50 kb non-overlapping genomic windows with options -m MFP -B 5000. Only SNPs from annotated regions (synonymous sites +nonsynonymous sites +introns) across all individuals were used. For diploid individuals, heterozygous sites were assigned IUPAC ambiguity codes and iqtree assigned equal likelihood for each underlying character, thus information from heterozygous sites is largely retained. This is crucial as we are interested in the branch length of inferred trees. Option -m MFP implements iqtree’s ModelFinder that tests the FreeRate model to accommodate maximum flexibility with rate-variation among sites (Kalyaanamoorthy et al., 2017). We also used UltraFast Bootstrap to calculate the support for different types of splits in each window (the -B 5000 option; Hoang et al., 2018). In each window, we extracted the support for monophyly among P. maackii +P.syfanius directly from the output of UltraFast Bootstrap, and we define the support for paraphyly among P. maackii +P.syfanius as (100 - the support for monophyly). For each level of support, we filtered out genomic windows where both the support for monophyly and the support for paraphyly drop below the given level. The remaining windows were considered informative.

We construct a continuous-time coalescent model as follows. Both populations have $N$ haploid individuals, gene flow rate is $m$, and coalescent rate is $N^{-1}$ in each population. In the equilibrium system, as we track both haploid individuals backward in time, there are six distinct states: ${\displaystyle (1|2),(2|1),(1,2|),(|1,2),(0|),(|0)}$, where 1 and 2 represent two individuals prior to coalescent, 0 is the state of coalescent, and $(\cdot |\cdot )$ shows the location of each lineage. Its transition density $\mathbf{B}(t)$ satisfies ${\partial }_t\mathbf{B}=\mathbf{AB}$, where $\mathbf{A}$ is given as

Let $S_{i|j}(T)$ be the mean sojourn time of an uncoalesced individual inside population ${\displaystyle i}$ during $0\le t\le T$, conditioning on the individual being taken from population $j$ at $t=0$. Assuming the infinite-site mutation model, let ${\mu }_i$ be the substitution rate in population ${\displaystyle i}$, observed rate ratio $r$ is thus

where (due to symmetry)

Let $r_0={\mu }_2/{\mu }_1$, and since $F_{ST}={(1+4Nm)}^{-1}$, we have

Software ELAI (Guan, 2014) with a double-layer HMM model was used to estimate diploid local ancestries across chromosomes. An example command is as follows: elai-lin -g genotype.maackii.txt -p 10 -g genotype.syfanius.txt -p 11 -g genotype.admixed.txt -p 1 -pos position.txt -s 30 -C 2 -c 10 -mg 5000 --exclude-nopos.

Note that -mg specifies the resolution of ancestry blocks, thus increasing its value will increase the stochastic error of incorrectly inferring very short blocks of ancestry. To control for uncertainty, we estimated repeatedly for 50 times. All replicates were used simultaneously in finding the correlation coefficients between entropy and other variables. Results from an example run is in Figure 3—figure supplement 1.

Here we introduce concisely the data transformation framework for calculating the entropy of local ancestry. The mathematical detail of this approach is presented in the appendix.

To characterize the average autocorrelation along an ancestry signal at a given scale $l$, define the following scale-dependent autocorrelation function

where $z(l)$ is understood as a periodic function such that $z(\xi +l)=z(\xi +l-L)$ whenever the position goes outside of $[0,L]$. The Wiener-Khinchin theorem guarantees that $z(l)$’s power spectrum $\zeta (f)$, which is discrete, and the autocorrelation function $A(l)$ form a Fourier-transform pair. Due to the uncertainty principle of Fourier transform, $A(l)$ that vanishes quickly at short distances (small-scale autocorrelation) will produce a wide $\zeta (f)$, and vice versa. So the entropy $S_w$ of $\zeta (f)$, which measures the spread of the total ancestry into each spectral component, also measures the scale of autocorrelation. In practice, $\zeta (f)$ is the square modulus of the Fourier series coefficients of $z(l)$, and we fold the spectrum around $f=0$ before calculating the within-individual entropy $S_w$. The formula used in the manuscript is

where $Z_n$ are the Fourier coefficients from the expansion $z(l)={\sum }_{n=-\mathrm{\infty }}^{+\mathrm{\infty }}{Z}_n{e}^{i2\pi nl/L}$. To speed up the Fourier expansion, we could densely pack equally-spaced markers that sample a continuous ancestry signal into a discrete signal, which then undergoes Fast Fourier Transform (FFT). The spectrum of FFT (discrete and finite) approximates the continuous-time Fourier spectrum (discrete and infinite), and entropy also converges as marker density increases.

As ancestry configuration is far from random around barrier loci, it will also influence the correlation of ancestry between different individuals at the same locus. For a genomic region experiencing strong barrier effects, two individuals could either be very similar in ancestry, or very different. This effect can be quantified by first calculating the cross-correlation $C_{j,j^{\prime }}(l)=z_j(l)\overline{z_{j^{\prime }}(l)}$ at position $l$ between individuals $j$ and $j^{\prime }$, and then averaging across a genome interval: $c_{j,j^{\prime }}=\frac{1}{L}{\int }_0^L{C}_{j,j^{\prime }}(l)dl$. The $J\times J$ dimensional matrix $\mathbf{C}$ with entries $c_{j,j^{\prime }}$ describes the pairwise cross-correlation within the cohort of $J$ individuals. We also have $c_{j,j}\equiv 1$ as each individual is perfectly correlated with itself. The matrix $\mathbf{C}$ is Hermitian, so it has a real spectral decomposition with eigenvalues ${\lambda }_j$ that satisfy ${\sum }_j{\lambda }_j/J=1$. This process is very similar to performing a principal component analysis on the entire cohort of individuals, and ${\lambda }_j/J$ describes the fraction of the total ancestry projected onto principal component $j$. If many loci co-vary in ancestry, the spectrum $\{{\lambda }_j\}$ will be concentrated near the first few components. Similarly, we use entropy to measure the spread of the spectrum, and hence the between-individual spectral entropy is defined as

The ancestry of a hybrid depends on the pure reference populations. The following assumptions are used throughout this section:

• There are a finite number of pure reference populations indexed by $k\in \{1,2,\mathrm{\cdots },K\}$. In most cases, we are only interested in $K=2$, for instance, a hybrid zone between two lineages.

• The chromosome has so many sites so that it can be treated as a contiguous rod with length $L$. $l\in [0,L]$ is the index of positions on a chromosome.

• The species’ ploidy is ${\displaystyle n_p}$ (completely phased data ⇔ $n_p=1$). When not specified, we assume the data is unphased, so that the ancestry on a particular chromosome always refers to the ancestry on a collection of ${\displaystyle n_p}$ homologous chromosomes.

• Ancestry can be inferred. In practice, regions with low density of informative SNPs will make inference difficult.

As "ancestry" is a categorical variable (i.e., there is no intrinsic order among the reference populations), to quantify the correlation of ancestry along a chromosome, we need to map the contribution of each ancestry category to some numeric value that can be used to calculate correlation.

The first choice is to use the probability of each category. For unphased data, let $p_k(l)$ be the fraction of the ${\displaystyle n_p}$-ploid chromosome at position $l$ coming from reference population $k$. The vector-valued function $\mathbf{B}(l)={(p_1(l),p_2(l),\mathrm{\cdots },p_K(l))}^{\top }$ completely describes the ancestry on a single chromosome. By definition, ${\sum }_k{p}_k(l)\equiv 1$. This is the conservation of total ancestry at any position. However, by this representation, the correlation ${\mathbf{B}}^{\top }(l_1)\mathbf{B}(l_2)$, which is the product of ancestries at different positions, will have units of the square of probabilities. If a locus is to be compared with itself, then the correlation will be

Ideally, we want the correlation of ancestry of a locus with itself to be the same regardless of its ancestry configuration. Equation (24) does not meet this criteria as it depends on $\mathbf{B}$. Here we borrow some concepts from signal-processing theory. The conservation of total ancestry means that ancestry is more analogous to the "energy" of a signal, which is in units of [squared-signal], rather than the signal itself. Therefore, we could instead define a second type of ancestry representation by taking the square roots of probabilities.

The hybrid index is defined as the average ancestry in a hybrid from a given reference population over a set of loci. It is usually calculated for each individual or for each chromosome. In our notation:

To characterize the scale of correlation, the distance between loci is important because correlation often drops while distance between loci increases. The information about the spatial scale of correlation is retained when the full spectrum of correlation is considered within a single individual.

A second source of correlation arises when we consider the relationship between individuals. To convey the main idea, let’s compare between a set of ancestry signals from an inversion and a set of ancestry signals from a regular region subject to normal recombination. For the inversion, since recombination between chromosomes of different ancestry is often completely suppressed, the ancestry signal along the inversion will be close to constant. When two haploid individuals are compared at the inverted region, they are either different everywhere in terms of ancestry, or completely the same. Comparing between diploid individuals is similar, although the difference is more fine-grained due to the presence of the heterozygotes. However, in a collinear region subject to a regular rate of recombination, any ancestry signal will switch stochastically between states. In the latter case, the similarity between two ancestry signals will be similar across multiple pairwise comparisons. It will be helpful to think in the principal component (PC) space spanned by all individuals’ ancestry signals. Ancestry signals from an inversion will form several tight clusters in the PC space, while those from a regular region will form a single cloud with a larger dispersion.

Both sources of correlation, and hence both types of randomness, can be measured using entropy in information theory.