Hybridization breaks species barriers in long-term coevolution of a cyanobacterial population

Previous analysis of 16S amplicons identified two main clusters similar to OS-A and OS-B’, but also several others at frequencies below 1% (Rosen et al., 2018). But whether clusters within 16S amplicons correspond to clusters at the whole-genome level in such a rapidly recombining population is not clear (Rosen et al., 2015; Miller and Carvey, 2019). Therefore, we first compared each single-amplified genome (SAG) to Synechococcus OS-A and OS-B’. Almost all SAGs (330 out of 331) comprised two clusters with an average nucleotide divergence of ∼2% from one of the reference genomes and ∼15% from the other. Note that the divergence between clusters is well above the 5% threshold that is commonly used to distinguish different species (Nayfach et al., 2016; Jain et al., 2018; Garud et al., 2019). We labeled the cluster closer to OS-A by α (128 SAGs) and the other by β (202 SAGs), to distinguish them from the reference genomes. Surprisingly, one SAG (approximately 70% complete) was ∼15% diverged from both OS-A and OS-B’ and formed a separate cluster which we labeled γ (Figure 1A). Clustering methods based on different metrics gave identical results (Appendix 2—figures 2 and 3).

Figure 1

Download asset Open asset

To investigate the extent that these clusters characterize the distribution of genomes throughout the two springs, we used single-cell genomes as references and recruited reads from a large collection of metagenome samples to determine the abundance of α, β, and, particularly, γ across different environments (see Appendix 3). We found that 95% of reads were part of a main cloud of sequences with ≤5% nucleotide divergence from one of the three main clusters. We used the fraction of reads within each main cloud to infer the abundances of each cluster in the metagenomic samples. The abundance of γ was highly correlated with α (Appendix 3—figures 3 and 4) and increased rapidly with temperature in the Mushroom Spring samples (Figure 1B). One interpretation of this pattern is that α, β, and γ represent distinct species adapted to different temperatures. Alternatively, genotypes within the same cluster could be adapted to different temperatures or other environmental conditions. The relative abundances of the clusters would then reflect the complex eco-evolutionary processes that maintain the diversity in the population.

Previous analysis using a much more limited dataset proposed that β formed a quasisexual population, in which rapid recombination led to random assortment of alleles within the cluster (Rosen et al., 2015). Motivated by these results, we investigated the impact of recombination on the diversity within the two main clusters.

We quantified recombination rates in α and β using a standard measure of linkage disequilibrium, ${\sigma }_d^2$. Given a pair of sites with alleles $\mathtt{a}\mathtt{/}\mathtt{A}$ and $\mathtt{b}\mathtt{/}\mathtt{B}$, ${\sigma }_d^2$ is defined as

where $⟨\cdot ⟩$ is the average across pairs of sites. In asexual populations with conventional demographic drift, ${\sigma }_d^2=5/11$ independently of the location of the sites. Recombination events affecting one but not both sites can unlink mutations, leading to a decrease in ${\sigma }_d^2$. Because sites further apart have more opportunities to become unlinked, ${\sigma }_d^2$ decreases with the distance between sites.

In β, the linkage disequilibrium ${\sigma }_d^2$ decreased by around a factor of 10 over distances of ∼1 kbp (Figure 2A). This agreed quantitatively with earlier results based on deep amplicon sequencing from Rosen et al., 2015, which showed a decrease by a factor of ∼1 over distances of up to 300 bp, using a different but related measure of linkage, conventionally denoted as $⟨r^2⟩$. Over genome-wide length scales, linkage was small (${\sigma }_d^2=0.06$) and close to a fully unlinked control (${\sigma }_d^2=0.04$). We found similar results for $\alpha$, but with a quantitatively slower decrease in linkage by a factor of ∼8 over ∼1 kbp. Nevertheless, genome-wide linkage values were still close to the fully unlinked control (${\sigma }_d^2=0.14$ for the data vs ${\sigma }_d^2=0.12$ for the control). As a consequence of rapid recombination, we found that diversity in marker genes, such as 16S rRNA, was poor predictors of whole-genome diversity (Figure 2B). Further details and analyses can be found in Appendices 4 and 8. Appendix 8 also provides estimates of the recombination parameters for each species and compares these to estimates obtained from previously proposed methods (Didelot and Wilson, 2015; Lin and Kussell, 2019). Altogether, these results are consistent with each cluster representing a quasisexual population and imply that phylogenies inferred from whole genomes do not reflect the evolutionary history of the population.

Figure 2

Download asset Open asset

While linkage decay in α and β was similar, we found diversity patterns at the gene level were qualitatively different (Figure 3A). The typical diversity at synonymous sites in β was ${\pi }_S\approx 0.04$, with ∼2× variation across genes (Figure 3A), consistent with previous results (Rosen et al., 2018). In contrast, ${\pi }_S$ in α varied wildly by over two orders of magnitude (Figure 3A), ranging from ${\pi }_S\approx 0.2$ in some genes to others that did not have a single mutation. Note that because here we only used synonymous sites, these variations cannot be explained by differences in the strength of purifying selection or the accumulation of non-synonymous mutations during rapid adaptation (Neher, 2022). Importantly, the low-diversity segments of the α genome spanned multiple samples and springs and are thus unlikely to be the result of a very recent clonal expansion as is often seen in human commensal and pathogenic bacteria (Didelot and Wilson, 2015; Garud et al., 2019; Sakoparnig et al., 2021; Weimann et al., 2024).

Figure 3

Download asset Open asset

Beyond the low-diversity genes in $\alpha$, we found genomes in both species containing segments very similar or identical to homologous regions from the other species. We performed extensive checks to verify that these segments were not assembly artifacts or cross-contamination between reaction wells during whole-genome amplification (Appendix 10). Such segments could represent accessory genes transferred between species, as has been widely reported in both Synechococcus and other species (Bhaya et al., 2007; Smillie et al., 2011). Alternatively, they could be the result of homologous recombination with other species within the core genome. We refer to the later process as hybridization throughout. The rest of our analysis aims to quantify the extent of hybridization and its impact on the evolution of the population. We began by identifying the core genes shared by α and β.

We defined core genes by constructing groups of orthologous genes (orthogroups) and identifying those that were likely present in all genomes. Specifically, we clustered gene sequences in several steps using standard methods based on similarity (Appendix 1). We used phylogenetic distances to identify and separate distant paralogs within the same orthogroup. To assign clusters corresponding to each species, we analyzed the abundance distributions of orthogroups in α and β separately. Because of the incomplete nature of our genomes (average coverage ≈0.3, with a range of ≈0.05–0.95; see Appendix 1), the assignment has to be done statistically. Distributions for both species were bimodal, with a peak near zero corresponding to flexible genes and another at high abundances corresponding to core genes. We used a binomial fit to determine single-copy genes that were likely present in all genomes and found 1825 core α and 1737 core β orthogroups. Of these, 1448 orthogroups were shared by both α and β and make up the Synechococcus core genes. We focus the rest of our analysis on these core genes.

We first performed a quantitative comparison of the synonymous diversity between α and β across the core genes. We hypothesized that the low-diversity segments of the α genome were the result of mutation accumulation since their common ancestor. If mutation rates across genes are similar, we expect the fraction of polymorphic sites $f_p$ for each gene to follow a Poisson distribution. Empirically, we found that for $f_p<0.05$, the data closely matched the Poisson expectation with a single fit parameter (Figure 3B). This group, which we call the ancestral backbone of$\alpha$, contained around half of the core genes (675 out of 1448). Note that this result is distinct from the distribution of divergences between species, which was shown previously to also follow a Poisson distribution (Rosen et al., 2018). The remaining half (773 out of 1448) had up to ten times higher fraction of segregating sites and were inconsistent with the backbone polymorphism rate. As our subsequent analysis shows, the high diversity in these genes is the result of hybridization (see Figures 4—7 and subsequent discussion). We therefore refer to this group as the hybrid α genes. In β, the distribution of $f_p$ had a single peak, but a simple Poisson process was a poor fit to the data (Figure 3C). This suggests a substantial fraction of the diversity is not simply due to mutations but is introduced through recombination with other strains or species. Consistent with this interpretation, we found the distribution of synonymous diversities in the hybrid α overlapped that of β, and both were two orders of magnitude higher than in the α backbone (Figure 3D). These results suggest that hybridization had a large impact on the genetic diversity of both α and β. We therefore developed a method to identify hybrid segments directly and test this hypothesis.

Figure 4

Download asset Open asset

Figure 5

Download asset Open asset

Figure 6

Download asset Open asset

Figure 7

Download asset Open asset

We used hierarchical clustering to identify and characterize hybridization patterns between species. If hybridization is common, as in our case, inferences based on phylogenetic trees would not provide reliable estimates of pairwise divergences. We therefore used pairwise nucleotide distances directly. The distribution of pairwise distances was bimodal, with a clear gap around $d=0.075$, consistent with our comparisons to the reference genomes (Figure 1). Based on this observation, we defined species clusters for each gene using a simple cutoff at $d_c=0.075$ (Appendix 1—figure 2) together with additional constraints to ensure that clusters were well separated and consistent with asexual divergence from a common ancestor (Appendix 1).

Detailed analysis of species clusters revealed evidence of hybridization in roughly half of the genome. We found two distinct patterns of hybridization. Around 75% of core genes (1078 out of 1448) had two to three clusters of sizes roughly proportional to abundances of α, β, and γ in the samples (Appendix 5—figure 1). The majority of these genes (844 out of 1078) had clusters, each of which contained sequences from a single species (Figure 4A), as expected if they were asexually diverged. We call these non-hybrid loci. A significant minority (234 out of 1078) contained a mixture of sequences from different species (Figure 4B), representing recent transfers of whole genes between species. We call these simple hybrid loci. The remaining 25% of core genes (370 out of 1448) formed either a single cluster or had one very large cluster containing most sequences (Figure 4C). As we show below, these large clusters result from complex patterns of hybridization within single genes. We thus call these mosaic hybrid loci. In around a quarter of these loci, we found evidence of complete or near-complete gene sweeps, but in most, the divergence between α and was typical. Almost half (297 out of 773) of high-diversity α genes were mosaic hybrids. To better understand the dynamics of hybridization in the population, we analyzed both hybridization patterns in detail, starting with the simple hybrid loci.

We used a simple set of heuristic criteria to assign the transfer direction of simple hybrids. If species evolved by asexual divergence from a common ancestor, we expect each orthogroup to contain two to three clusters of sizes roughly proportional to the species abundances. For example, based on the average coverage and the proportion of β in our samples (Figure 1A), we expect the largest cluster to contain ∼80 sequences from β genomes. In some cases (45 out of 1063), we also find a minority of sequences from $\alpha$ genomes in this β cluster, which we infer to be the result of a hybridization event from β into α (Figure 4B). An alternative explanation would be a transfer from a third species that underwent a full selective sweep through β and only a partial sweep through α. This scenario would imply that hybridization is even more common than we infer, so our method is conservative. We infer transfers from α into β in a similar way (Figure 4B). When the third cluster is present and contains γ, we assign all other sequences in the cluster to hybridization with γ. Because we only have a single γ sequence, these clusters are difficult to classify. We cannot distinguish between transfer from and to γ, and it is possible that some simple hybrid loci involving γ may actually be mosaic hybrids. When the third cluster is present but does not contain γ or in cases where we find a fourth cluster, we assign those sequences to hybridization with an unknown species $X$. Using this procedure, we systematically determine the frequencies of different hybridization events in α and β.

We found simple hybrids were present in 5−10% of core genes in both species (Figure 5A, B). Specifically, we found 185 simple hybrid loci in α and 68 in β. Note that we conservatively excluded contigs that were fully hybrid from this analysis, so our numbers likely underestimate hybridization rates in the population (Appendix 5).

In both species, simple hybrids were scattered along the genome and contained alleles from different donor species (shown in different colors in Figure 5A, B), implying multiple hybridization events occurred independently, as expected when recombined segments are much shorter than the genome. Around 80% of α hybrid alleles were present in multiple cells, while in β the fraction was lower, at around 50%. The frequencies of simple hybrids were highly variable, ranging from singletons to frequencies close to 100%, consistent with a continual hybridization scenario. Surprisingly, around half of hybrids in β and a fifth of those in α were from unknown donors. However, a systematic search over the 34 metagenomic samples did not find evidence of other Synechococcus species in either spring (Appendix 3—figures 5–7). These sequences could be the result of hybridization with species that were formerly present in the caldera and are now either extinct or at very low abundances. Alternatively, the hybrid sequences could have resulted from hybridization with, as yet, undiscovered species from other springs. But regardless of the exact scenario, overall, these results suggest the current population is descended from multiple cyanobacterial species that hybridized with each other.

Is the genomic mixing between species the result of neutral processes or does positive selection for beneficial hybrids play an important role? Distinguishing between these alternatives is difficult in general, so we instead focused on cases where beneficial hybrid genes swept through the entire host species, that is, full gene sweeps. Full gene sweeps would greatly reduce the synonymous divergence between α and β, making them much easier to detect. Visually, we observed the synonymous divergence between α and β contained several divergence troughs, in which the average diversity was well below typical values. We identified the genes contained in these troughs by imposing a simple cutoff $d_{trough}=0.2$ on the synonymous divergence between α and β, which we chose close to the inflection point of the synonymous divergence distribution across genes (data not shown).

A systematic search revealed 91 very low divergence loci, spread across more than half a dozen divergence troughs (Figure 5C). Comparing alignments of typical genes and divergence through genes revealed a dramatic difference in fixed substitutions between species, which can only be explained by a recent selective sweep. The largest divergence trough contained 27 genes, representing around 25% of divergence trough genes, most of which were part of a nitrogen-fixation pathway. In addition, the nitrogen-fixation pathway is one of the largest genomic segments in which the gene order between α and β is conserved (Bhaya et al., 2007), suggesting all genes were transferred simultaneously. The most likely explanation is that at least one of the ancestors acquired the ability to fix nitrogen after it colonized the Yellowstone caldera and that this gave it a substantial selective advantage. The genetic diversity within the pathway genes gave further clues about the hybridization process. The diversity within the trough was similar in both species (${\pi }_S=0.052$ in α and ${\pi }_S=0.039$ in β) and was close to the average diversity within β (${\pi }_S=0.042$), suggesting the β ancestor already had the ability to fix nitrogen when it colonized the caldera. However, the trough diversity was more than an order of magnitude larger than the diversity within the α backbone (Figure 6A, B). This large discrepancy strongly suggests α acquired the pathway from β through multiple independent transfers after they colonized the Yellowstone caldera.

We performed a similar analysis of diversity within simple hybrid genes to determine the direction and timing of transfers. We focused our analysis on α, where we had more data. Presence of hybrid genes at high frequencies suggests the hybridization process may be driven by selection. Two distinct patterns of diversity could result. If a single allele was transferred into α and then underwent a selective sweep, we would see very little diversity among the hybrid alleles. This is known as a hard sweep. Alternatively, if multiple hybrid alleles were transferred into $\alpha$ while undergoing a sweep, the diversity would be comparable to that in the host species. This is the pattern we observed in the divergence troughs and is known as a soft sweep.

We found evidence that both hard and soft sweeps contribute to hybridization in the population (Figure 6C). Around 15% of α simple hybrids (31 out of 153 genes) had no diversity at all, regardless of the donor species, consistent with hard sweeps driving their spread. In the remaining 85% (122 out of 153 genes), diversity within α was 1–2 orders of magnitude higher than the α backbone diversity. As a result, the diversity could not be due to mutation accumulation since the common α ancestor, but must have originated from multiple independent transfers—that is, a soft sweep. We further verified the presence of soft sweeps within α by comparing the diversity of individual genes. In most cases, the diversity within α at these loci was similar to β, with a handful of cases of α having significantly less, consistent with a subset of β alleles sweeping through the population simultaneously (Figure 6C).

Selective sweeps of individual genes could also lead to mixing of other genomic segments from the donor species into the host population through genetic hitchhiking. If the diversity within the two hybridizing species is low, this process would result in blocks of tightly linked SNPs consisting of primarily two haplotypes, one from the host and another from the donor (Figure 7A). Anecdotally, we frequently observed such patterns in mosaic hybrid genes (Figures 4C and 7A). We developed a systematic method to search for SNP blocks based on clustering the SNP linkage matrix. Briefly, we search for consecutive SNPs that have perfect linkage (see Appendix 6). We used α hybrid genes to calibrate our method and chose very conservative parameters to avoid false positives. Consistent with this, only around 17% of SNPs from the α hybrid genes belonged to blocks.

Despite the very conservative thresholds used in our algorithm, we found over one thousand SNP blocks in each species, scattered across hundreds of genes, including in highly conserved ribosomal genes and genes in which distinct species clusters could be assigned (Figure 7BC; Appendix 6—figures 1 and 2). Most blocks were short, with typical lengths around 70 bp in α and around 30 bp in β (Appendix 6—figure 1). The divergences between SNP block haplotypes were similar to that between species, consistent with hybridization (Appendix 6—figure 3). To confirm their hybrid origin, we mapped the haplotypes from each block to the two other species in our sample. We found 66% (406 out of 605) of SNP blocks within α and 15% of those within β (128 out of 776) contained a haplotype identical to consensus of other species (Figure 7B). Note that we only considered blocks in which all three species were present for this comparison. In addition, we found a comparable number of blocks that were statistically similar but did not match the other species α, β, and $\gamma$ (Figure 7B, C). The most likely explanation is that these blocks resulted from hybridization with other species, similar to the whole-gene hybrids found previously (Figure 5).

How much linkage is maintained between blocks? Given the relatively high divergences between species, one might expect epistatic interactions between blocks to constrain the genotypic structure of the population. Conversely, if recombination is dominant, each genome would be a mosaic of segments from the two ancestral species. To answer this question, we calculated the linkage coefficients between haplotypes for each pair of blocks and compared their distribution to a fully unlinked control (see Appendix 6 for details). Within the β species, the two distributions were almost identical, implying that the genomes are consistent with a random mixture of block haplotypes, as was previously proposed (Rosen et al., 2015). Within $\alpha$ the linkage coefficients were statistically larger than those predicted from the unlinked control (Figure 7D), but the overall value was still low, with $⟨r^2⟩\approx 0.12$ for data compared to $⟨r^2⟩\approx 0.05$ for unlinked control. These results are consistent with the slower rate of linkage decrease within $\alpha$ compared to β discussed previously. We also found similar results for the linkage between hybrid gene alleles (Appendix 4). Thus, rather than finding clonal subpopulations, every one of the ∼300 cells in our sample has a unique combination of SNP block haplotypes.

In this section, we outline the method we used to assign genomes to each species cluster. We also verify the robustness of the species assignment by using three different methods and show that they give identical results.

Comparisons of genomes from bacterial populations often reveal clusters at different length and divergence scales (Jain et al., 2018). But the interpretation of such clusters is controversial (Cohan, 2002; Fraser et al., 2009). There are three main possible interpretations for such patterns. One interpretation is that the clusters represent the bacterial analog of eukaryotic species, which are sexually and ecologically isolated from each other. In practice, a whole-genome divergence cutoff of 5% is frequently used to define bacterial species (Jain et al., 2018). An alternative interpretation is that clusters are the result of clonal expansions of particular strains from a much more diverse and possibly more genetically homogeneous population. The third interpretation is that the clusters are the result of uneven sampling from a larger and more diverse population. The unique combination of having a large collection of genomes, with minimal compositional bias, from a geographically isolated population allows us to test these alternatives. We will mainly aim to distinguish between the first two interpretations and the third one by checking whether there are robust genomic clusters in our sample. More discussion on the evidence concerning the first two hypotheses can be found in Appendices 3 and 5.

We used three different methods to check whether the genomes form robust clusters: the phylogenetic tree of the 16S rRNA sequences, average genome-wide divergence from the two reference genomes OS-A and OS-B’, and patterns of divergences between triplets formed by the SAG and the two reference sequences at different loci. The phylogenetic tree of the 16S is the simplest and still the most commonly used method to characterize the diversity of natural microbial communities (Yarza et al., 2014). Because of its practical importance, we included it for comparison with the other two methods. However, in a highly recombining population like the Yellowstone Synechococcus, it is not clear to what extent the 16S divergence accurately reflects the whole-genome diversity. To address this question, we used the divergences across the entire genome as a comparison. The main limitation of using summary statistics such as the average divergence across the entire genome is that it obscures heterogeneities along the genome that could have important ecological and evolutionary implications. There are different ways to address this concern. Here, we used an extension of a method first used in Rosen et al., 2018, which provides a coarse description of the variation in the diversity patterns across different loci. We present the results from each of these analyses next.

We first constructed the phylogenetic tree of 16S sequences. Predicted 16S rRNA sequences (117 from 331 filtered SAGs) were extracted and aligned using MAFFT (v7.453) (Katoh and Standley, 2013). Two sequences (SAGs MA02M11_3 and MuA02C9) were incomplete and were removed from the analysis. The remaining (115) sequences were hierarchically clustered using the average linkage (UPGMA) method implemented in Scipy. Results are shown in Appendix 2—figure 3, with the colored background showing the species assignment using whole-genome divergences (see below) for comparison. We found perfect agreement between the whole-genome assignment and the largest three clusters in the 16S phylogenetic tree.

The average divergence within the α cluster is approximately half that of β ($1\cdot {10}^{-3}$ vs $2\cdot {10}^{-3}$). However, a significant fraction of the average divergence within β was due to two outliers (MuA1L16 and MusA1J8). After removing the outliers, the average divergence within β decreased to $1.7\cdot {10}^{-3}$. Closer investigation of the alignments revealed that both outlier β sequences contained a short segment of three SNPs that were identical to all α sequences and different from all other β sequences, at positions 1067–1072. In addition, MuA1L16 contained a larger segment of almost 200 bp (positions 390–570) containing 10 SNPs and an insertion that were different from all other β and were identical to the consensus $\alpha$. Apart from these segments, the outlier sequences were typical of the β species. This pattern is similar to the SNP blocks discussed in the main text and in Appendix 6 and therefore is likely the result of hybridization between α and β.

We compared the 16S rRNA clusters to whole-genome divergences from the two reference genomes, OS-A and OS-B’. We aligned all of the predicted genes from each SAG to all of the coding DNA sequences (CDS) from each reference genome using BLAST (v2.13.0+) (Altschul et al., 1990). For each SAG, we calculated the divergence to each reference genome ($d=1-\text{pident}/100$) and calculated the average over the best hits across all genes. The results for all SAGs with more than 100 hits to both references are shown in Figure 1A. Three clear clusters were observed, which we labeled as α, β, and γ. The species of each SAG was assigned as the label corresponding to the cluster it belonged to in Figure 5A.

The previous two metrics do not reflect possible heterogeneity in the patterns of diversity along the genome. To characterize this variation on a gene-by-gene basis, we extended a method that was previously used to analyze amplicon data from the Mushroom Spring in Rosen et al., 2018. We refer to this method as gene triplet analysis. This analysis exploits the separation of scales between and within species clusters. At most loci and in most cells, we observed previously in order to coarse grain the high-dimensional matrix of pairwise divergences into a small number of discrete patterns. Specifically, we take triplets of sequences from the SAGs, OS-A, OS-B’, at loci where all three are homologous and calculate the divergences between all three pairs of sequences. These divergences typically form a triangle as shown in Appendix 2—figure 1. Note that it is possible to have triplets of divergences that violate the triangle inequality, which cannot be represented as triangles in a two-dimensional space. In our data, such cases were rare, but this may not be generally true, particularly for sequences that are not clearly separated and which are the result of extensive recombination. Note, however, that the classification into distinct patterns given below can still be done even in these cases. To reduce the dimensionality of the data further, we define a main cloud diameter $d_M=0.1$ for the diversity within species. We then classify the shape of the triangle into the 5 distinct patterns shown in Appendix 2—figure 1, depending on whether each side of the triangle is longer or shorter than $d_M$. For example, denote the divergence triplet by $(d_{XA},d_{XB},d_{AB})$, representing the divergences between the SAG sequence and OS-A, the SAG sequence and OS-B’, and OS-A and OS-B’, respectively. Then the first pattern in Appendix 2—figure 1A (labeled ‘XA–B’) contains loci where $d_{XA}<d_M$, $d_{XB}\ge d_M$, and $d_{AB}<d_M$, the second pattern (labeled ‘XAB’) contains loci where $d_{XA}<d_M$, $d_{XB}<d_M$, and $d_{AB}<d_M$, and so on. The cutoff of $d_M=0.1$ is the same as the one used in Rosen et al., 2018 and provides a good separation between species at a gene by gene basis. The results of our analysis, however, do not depend on the specific value of $d_M$ as discussed at the end.

For each SAG, we assigned all genes that were orthologous to OS-A and OS-B’ to one of the five patterns described above. The three rows in Appendix 2—figure 1 show typical examples of an α (Appendix 2—figure 1A), a β (Appendix 2—figure 1B), and the γ (Appendix 2—figure 1C) SAG. The majority of genes across each row in Appendix 2—figure 1 are found in the XA–B, XB–A, and X–A–B patterns, respectively. This is consistent with the species clusters we saw in the genome-wide average divergences from before. However, we also find a significant minority of genes forming other patterns, most notably XAB, representing genes where all three sequences are within the main cloud ($d<d_M$). Some of these genes were contained in the divergence troughs discussed in the main text. In the case of the γ SAG, the fact that X–A–B is the dominant pattern rules out the possibility that γ is a hybrid resulting from the mixture of $\alpha$ and β. Such hybrids have been observed after recent clonal expansions in acid mine drainage microbial communities (Tyson et al., 2004; Lo et al., 2007), but do not appear to be present at substantial frequency in the Yellowstone community. Instead, it appears that γ is a distinct species that diverged from the other two around the same time α and β diverged from each other.

Using the pattern described above allows us to compare different SAGs to each other and allows us to test the robustness of the species assignment based on the average genome-wide divergence. Specifically, we assigned to each SAG a gene triple fingerprint $f=(f_{XA-B},f_{XAB},f_{X-A-B},f_{X-AB},f_{XB-A})$, representing the fraction of genes with each pattern (note that $\sum _i{f}_i=1$). We then compared the fingerprints of each SAG with the species assignment based on the average divergence from OS-A and OS-B’ from Figure 1A. We found that the fingerprints formed three clear patterns that exactly matched the assignments based on the average divergence (Appendix 2—figure 2). We tested this assignment for different values of $d_M$ ranging from 1% to 20% and found that the fingerprint patterns were robust across the entire range (data not shown). These results show that the classification into distinct species is robust and consistent across the three methods described.

Appendix 2—figure 1

Download asset Open asset

Appendix 2—figure 2

Download asset Open asset

Appendix 2—figure 3

Download asset Open asset

To estimate abundances of the three main species α, β, and γ in the Octopus and Mushroom Springs, we recruited metagenomic reads from 34 samples taken from both springs to SAGs from each species. We chose the one γ SAG (OA3L13, ≈70% coverage) along with representative SAGs from α (MA02D15, ≈77% coverage) and β (OcA2H14, ≈96% coverage) as reference genomes. Hybridization on length scales longer than the typical read size (150 bp in our case) can lead to errors in our estimates of species abundances. To minimize these errors, we only used core genes where three distinct species clusters could be identified and there was no evidence of hybridization in the single-cell data (see Appendices 1 and 5). This resulted in a set of ≈400 genes for each reference. We aligned all the metagenomic reads to these target gene sequences using BowTie2 (v2.3.4, default parameters) (Langmead and Salzberg, 2012) and extracted reads aligning to each of the three reference sequences using Samtools (v1.7.0) (Li et al., 2009) and bedtools bamtofastq command (v2.26.0) (Quinlan and Hall, 2010).

We assigned the recruited reads to each species cluster if their divergence from their respective reference was less than the ‘main cloud’ radius $d_M$. We determined $d_M$ self-consistently from the data as follows. For each reference, counted the number of reads that had divergences less than $d_M$ from that reference sequence and more than $d_M$ from the other two reference sequences, for $d_M$ from 1% to 15%. If the population consists of species clusters with typical diversity $\pi$, the number of reads within the cloud will increase until $d_M\approx \pi$. As $d_M$ approaches the typical divergence between species $d\approx 15\mathrm{\%}$, more reads will start to map to multiple reference sequences, so the number of reads uniquely assigned to each species cluster will decrease. Therefore, the number of reads uniquely assigned to each species cluster should have a maximum as a function of $d_M$. We indeed found this to be the case for all three species (data not shown). Empirically, we found that $d_M=5\mathrm{\%}$ was close to the maximum for all three species, so we used this value throughout our analysis.

We found that the average depth of recruited reads across genes varied within the same sample. However, the variations across samples were highly correlated, which confirmed that most reads aligning to our chosen genes came from genomes that were part of the same species cluster rather than from genes transferred across species (Appendix 3—figure 3). We found that the genes that recruited the largest numbers of reads had read depths that were much more similar within the same sample compared to other genes. Interestingly, we found much more variation in depth between genes within the γ cluster compared to α and β (Appendix 3—figure 3). One likely explanation for this is that γ contains much more diversity compared to either α or β.

Appendix 3—figure 3

Download asset Open asset

We estimated the abundance of each species cluster by averaging the depth of the 100 genes with the highest average depth across all samples for each cluster. We found that γ was more than an order of magnitude more abundant than our detection threshold given by the metagenomic read depth in 32 out of the 34 metagenome samples (Appendix 3—figure 4). Moreover, its abundance in the two samples for which we had single-cell and metagenome data (MSe4 and OSM3) was close to the detection threshold of the single-cell samples, which was around 2%. This explains why we only obtained a single γ cell from the single-cell samples after our quality controls. Finally, we found that the abundance of γ was highly correlated with α across all samples and was only present at abundances higher than $O(1)\mathrm{\%}$ in the Mushroom Spring sample taken from 65°C. This suggests that γ is, together with α, primarily adapted to growth at higher temperatures.

Appendix 3—figure 4

Download asset Open asset

Are other Synechococcus species present in Mushroom or Octopus Spring apart from α, β, and γ? To answer this question, we again used the 34 samples from which we have high-depth metagenomes. We used BLAST with highly sensitive parameter choices to minimize the number of true Synechococcus reads not captured. For computational efficiency, we focused on the 100 genes with highest coverage depth in the metagenomes across all three references.

To determine how much of the population is in the α–β–γ main cloud, we used the same three SAGs as in the previous subsection as references and recruited reads using BLAST (-task blastn -evalue 1e-3 -word_size 6 -num_threads 32). Using these parameters, we could recruit reads belonging to the dominant Roseiflexus species in our samples, which is in a different phylum from Synechococcus. We were therefore confident that any reads from other potential Synechococcus species would be included in our analysis. To limit the analysis to just Synechococcus species, we removed all reads that were more than 20% diverged from all of our references from subsequent analysis.

Because of the large amount of data in this analysis, we focus on just the Mushroom Spring samples taken from different temperatures (MS temperature series). We recruited a total of 4,522,833 reads with an average depth per gene of 4980. We calculated the minimum divergence from the three reference sequences for each read and found 94% within 5% divergence of one of them (Appendix 3—figure 5A). For comparison, we determined the typical divergences between species by dividing the reference sequences into 150 bp segments (the same size as our reads) and calculating the divergences between all three pairs of species (Appendix 3—figure 5B). The middle 95-percentile of the combined distributions extended down to a divergence of around 5% (eight SNPs), consistent with our cutoff for the intraspecies divergences at these loci. These results are consistent with all of the Synechococcus reads being from one of the three species and set an upper bound on the relative abundance of any other species at below 4%.

To check how sensitive our detection method is to the presence of other species, we used the Mushroom Spring sample from 65 in which γ had a relatively high abundance. To mimic the presence of another species, we compared the previous distributions to the distribution of the minimum divergence from α and β only (dashed red line in Appendix 3—figure 5A). We found a 19% decrease in cumulative number of reads within 5% from the references in this case, which was very close to the 21% abundance of $\gamma$ we estimated in the previous subsection based on the average depth across loci. This shows that our method is sensitive to the presence of other species and illustrates the benefits of combining single-cell and metagenome analysis for identifying unknown species.

Appendix 3—figure 5

Download asset Open asset

To further verify that no other Synechococcus species were present in our samples, we closely examined the diversity of 16S sequences. We used only single reads in our analysis instead of trying to assemble full sequences to prevent misinference caused by hybridization. Because the diversity along the 16S sequence is highly heterogeneous, variation between read recruitment locations can dominate the variation between species. We controlled for this by limiting the analysis to reads that aligned within a given region of the 16S. For 150 bp reads, we found that using segment sizes of 200 bp provided the best compromise between increasing depth and minimizing variance due to alignment location. To find the most informative region, we calculated the nucleotide diversity at each site ${\mathcal{H}}_i=2⟨f_i(1-f)⟩$, where $f_i$ is the minor allele frequency at site $i$, and averaged it over all 200 bp segments along the 16S (Appendix 3—figure 6A). We found the diversity varied by more than a factor of four, consistent with the presence of conserved and variable regions. The diversity was nearly identical when using all three species as reference genomes, suggesting our recruitment captured all sequences within the main cloud. The highest average diversity was close to 400–600 bp segments, which we chose for further analysis.

We calculated the distribution of divergences from the γ allele for all reads recruited from this segment for MS temperature series. We found temperatures of 60°C and below had two distinct peaks close to Appendix 3—figure 6B 4% and 8% divergences. For comparison, we calculated the range of divergences from α and β for all 150 bp subsegments within the 400–600 bp segment and found these matched with the two peaks we observed. At 65°C, we also found only two peaks, one close to zero and another close to 4%, consistent with our previous estimates showing that this sample mainly consisted of α and γ cells.

Appendix 3—figure 6

Download asset Open asset

Do peaks correspond to single-species clusters or are they amalgamations of several distinct species with similar divergences from each other? To answer this question, we constructed a multiple-sequence alignment of all of the 16S reads recruited across the 34 metagenome samples and performed a hierarchical clustering on the pairwise divergence matrix using average linkage. Pairwise divergences were calculated by padding the reads with gaps to ensure all sequences were 200 bp in length and calculating the nucleotide divergence ignoring any gaps. For the MS temperature series, we found 97% of reads were in one of the three clusters associated with α, β, or γ (Appendix 3—figure 7A–D). Consistent with our previous results, samples below 65°C had two main clusters close to OS-A and OS-B’ (Appendix 3—figure 7A–C), while at 65°C they were close to γ and OS-A (Appendix 3—figure 7D). Aggregating across all 34 metagenome samples, we found that α and β were the two dominant clusters, with γ being the next largest (Appendix 3—figure 7E, F). We also observed a large difference in the overall abundances of α and β across the two springs, but the very large heterogeneity across samples we found previously (Appendix 3—figure 4) suggests this is likely due to the particular samples used in this study and may not necessarily be an accurate reflection of the overall abundances in the two springs.

Appendix 3—figure 7

Download asset Open asset

In this section, we summarize our analysis of linkage disequilibrium within α and β. The structure of this section is as follows. We begin by reviewing definitions of two commonly used ways to quantify linkage disequilibrium. We then compare the rate of decay of linkage disequilibrium within and between the two species. We next compare the decrease in linkage in the data to theoretical predictions from the neutral model and show that while the decay is similar to the neutral model prediction, the overall shape of the linkage decay curves cannot be fully explained by the classical neutral model. We then show that within α, there is substantial heterogeneity in the rate of linkage decay in low-diversity regions compared to high-diversity ones. Surprisingly, we find that linkage decreases more rapidly in low-divergence regions, suggesting the absence of a clonal frame for α even on short time scales. The final subsection presents a technical validation of our method for processing alignments showing that our results are robust to filtering for outlier gene alleles.

We used two distinct but related methods for quantifying linkage disequilibrium, ${\sigma }_d^2$ and $r^2$. Both metrics are related to the correlation coefficient between alleles at two loci but differ in their normalization. For a single pair of loci A and B with alleles A/a and B/b, we can define the correlation coefficient between their respective frequencies as

Using this expression, we can define a statistical measure of linkage as a function of distance along the genome by averaging over pairs of loci with a given separation $x=|x_A-x_B|$, where $x_A$ and $x_B$ are the positions of the two loci along the genome. The two most common ways to perform the average are

and

For a neutral panmitic and well-mixed population, both metrics decrease as ${\sigma }_d^2(x)\sim r^2(x)\sim 1/x$. For ${\sigma }_d^2$, there also exists an analytic expression (Ohta and Kimura, 1969):

where $N$ is the population size, and $R$ the recombination rate per base pair per generation. When recombination is rare and occurs via crossover $R(x)\approx R_{0}x$, and the expression is consistent with the scaling ${\sigma }_d^2(x)\sim 1/x$ mentioned previously. While numerical and theoretical studies have shown that $r^2$ is more sensitive to low-frequency alleles compared to ${\sigma }_d^2$ (Hudson, 1985; Song and Song, 2007), our analysis did not reveal significant differences between two metrics. We therefore focus here on ${\sigma }_d^2$ which is more commonly used and easier to compare to theoretical predictions.

We calculated ${\sigma }_d^2(x)$ and $r^2(x)$ for α and β separately and together. Both procedures are similar, so we focus on the analysis for all α and β sequences. Unless otherwise stated, all analyses were performed on alignments of individual genes and the results were averaged across genes. Using this method, we can easily obtain linkage information over distances of up to around 1 kbp. Over longer distances, the results can be dominated by specific recombination events in a small number of unusually long genes and are therefore less useful for quantifying genome-wide statistical patterns. We thus perform most of our analysis up to distances of around 1 kbp, except when calculating typical genome-wide linkage as discussed below.

Alignments for each orthogroup were first trimmed to remove segments with excess numbers of gaps. We then determined a consensus sequence by taking the most common codon sequence at each location along the reading frame, excluding gaps. Each nucleotide was labeled 0 if it was a gap, 1 if it was the same as the consensus, and 2 if it was different from the consensus. While this does not distinguish between biallelic and multi-allelic sites, the latter are much less common than the former and should not affect our conclusions. Finally, we calculated the two linkage coefficients using Equations 3 and 4. To estimate linkage across distances longer than a single gene, we performed the average across SNPs spanning a pair of genes and then averaged the values across a random sample of 1000 pairs of genes.

The linkage calculations within each species were performed in an analogous way, with the exception of a two-step filtering process of the alignments at the beginning. We describe this procedure in detail at the end of this section.

Appendix 4—figure 1

Download asset Open asset

Our analysis showed that within both α and β, linkage quickly decreased by approximately an order of magnitude over a distance of ∼1 kbp (orange circles and blue squares in Appendix 4—figure 1). However, when considering the population as a whole, we find linkage disequilibrium is maintained across the whole length of the genome (black diamonds in Appendix 4—figure 1). Within β, the linkage at distances of 1 kbp was identical to the genomewide linkage between random pairs of genes. Within α, the decrease was slower, but followed a similar pattern. The value of the genome-wide linkage was slightly higher than at 1 kbp, likely because of the reduced sample size at long genomic distances caused by poor genome coverage. Consistent with this, we found that the genomewide-linkage values in both species were close to completely unlinked controls obtained by randomly permuting each column of the alignments (${\sigma }_d^2=0.14$ vs ${\sigma }_d^2=0.12$ in α and ${\sigma }_d^2=0.06$ vs ${\sigma }_d^2=0.04$ in β). Overall, these results suggest that the dynamics leading to the reassortment of alleles are similar in the two species and only differ in the length scale over which alleles become unlinked. We discuss this in more detail in the next section.

The fact that both species have similar levels of genomewide linkage despite the sample size of α being less than half of that of β shows that our estimates for the genomewide linkage are not significantly affected by the sample size effect. We further validated this by taking random samples from the β main cloud equal to the size of the α main cloud at each locus and recalculating the linkage curves. We found the results for the full β main cloud and the subsampled main cloud were almost identical (Appendix 4—figure 2). This shows that the sample size does not have a large effect on our results and cannot explain the differences in linkage decay between the two species.

Appendix 4—figure 2

Download asset Open asset

The rapid decrease in linkage sets important limits on how much information about the diversity in the population is contained in samples of marker genes. To illustrate this, we performed a linear regression between the whole-genome divergences and the 16S rRNA divergence. Consistent with our previous analysis, we found slightly stronger linkage in α compared to β, but in both cases, the correlations were weak and the 16S divergences explained less than 10% of the genome-wide divergences (Appendix 4—figure 1B). These results show that frequent recombination in the population makes marker genes, such as the 16S rRNA, poor predictors of the genome-wide diversity.

To get a more quantitative understanding of the effects of recombination on the genetic diversity in our population, we compared our previous results to theoretical predictions from population genetics.

At short distances, the prediction from Equation 5 is that ${\sigma }_d^2=5/11$. However, in our data, the value of ${\sigma }_d^2$ for nearby sites was very close to one. While the Yellowstone Synechococcus have a strong mutational bias toward GC nucleotides (Rosen et al., 2018), Equation 5 is independent of the relative rates of mutations. As a result, we conclude that the neutral model cannot explain the linkage curves we observe.

Appendix 4—figure 3

Download asset Open asset

While the overall level of ${\sigma }_d^2$ in our data is inconsistent with the neutral prediction, we can still ask whether the rate of decay with distance has a similar functional form. To do this, we followed Garud et al., 2019 and fit the linkage curves to the following functional form:

where $\rho =2NR$ is the population recombination rate. We performed the fit by choosing $A$ such that the above expression matched the value of ${\sigma }_d^2$ in the data for neighboring sites ($x=1$). We then fitted the curves by eye by adjusting the value of $\rho$.

We found that this functional form could fit the data very well for both species up to distances of $x\approx 300$ bp (Figure 2 in main text). For longer distances, the linkage decrease was slightly slower than predicted by Equation 6. One possible interpretation of this slowdown is that longer distances have a significant overlap with the distribution of lengths of recombined segments $L_r$, beyond which we expect linkage to reach the genomewide level (Garud et al., 2019). This suggests that recombined segments as short as 300 bp might play an important role in the recombination process within the population.

We further validated that the linkage curves have the same functional form by rescaling distances between sites by the recombination rate $\rho$ inferred previously. We found a very good collapse of the linkage curves for the two species when plotted against the rescaled distances (Appendix 4—figure 3B). This strongly implies that the dynamics leading to the decay of linkage are the same in both species, with the differences between them being length scale over which sites are effectively linked.

Our analysis so far was limited to either distances up to the length of typical genes (∼1 kbp) or across the entire genome. Extending our analysis over distances in between these two extremes is difficult due to several technical challenges. The first is the fact that typical coverages for our SAGs are around 30%. Together with a typical contig size of ∼10 kbp, this means that our sample size quickly decreases with distance. The second challenge is that at longer distances, gene order might vary across genomes, making it difficult to define a genomic distance at the population level.

We calculated the linkage disequilibrium between multi-SNP features, such as SNP blocks or hybrid alleles, using Equation 2 as in the case of biallelic SNPs. In the case of SNP blocks shown in Figure 7D in the main text, we calculate the frequencies for the four haplotypes within SAGs that are covered at both loci. The unlinked controls (represented by solid lines in Figure 7D) were obtained by a random permutation of the alleles at each locus independently, while keeping the SAGs that are covered at each locus fixed. This procedure preserves any correlations in genome coverage between nearby loci and provides the most accurate estimate of $r_{AB}^2$ for a completely unlinked genome.

The linkage between hybrid gene alleles was calculated similarly as for the SNP blocks, except for the genotype assignment. At each locus, we labeled the most abundant allele using upper case letters ($A,B,\dots$) and grouped all of the remaining ones into a single allele labeled by a lower case letter ($a,b,\dots$). The linkage matrix $r_{AB}^2$ and the unlinked control were calculated the same as previously described. The results for α were $⟨r_{AB}^2⟩=0.04\pm 0.15$ compared to $⟨r_{AB}^2⟩=0.02\pm 0.11$ for the unlinked control. For β, the results were $⟨r_{AB}^2⟩=0.01\pm 0.09$ compared to $⟨r_{AB}^2⟩=0.009\pm 0.08$ for the unlinked control. Note that the average linkage disequilibrium was close to the unlinked controls in both species, with α having slightly higher linkage. These agree with those from the analysis of SNP blocks shown in the main text.

Appendix 4—figure 4

Download asset Open asset

For each species, we first selected all of the sequences that belonged to SAGs assigned to that species (see Appendix 2). We then calculated the consensus sequence for this group of sequences, as described above, and removed all sequences with divergences greater than a main cloud cutoff $c$. This is a simple way to remove sequences that are the result of gene-level hybridization (see Appendix 5 and Rosen et al., 2015 for a similar procedure). We chose the main cloud cutoff for each species by calculating the linkage curves for different values and choosing the largest value for which the curves did not noticeably change (Appendix 4—figure 4). Based on this criterion, we use $c=0.05$ for α and $c=0.1$ for β throughout.

This section describes the method we used to quantify hybridization at the gene level. Our method uses a set of simple heuristics to assign labels to the species clusters defined in Appendix A based on the composition of species to which the sequences belonged. We then identified hybridization events by searching for discrepancies between the gene and whole-genome labels as explained below.

We first assigned labels to the species-level clusters constructed as described in Appendix 1. Based on the gene triplet results presented in Appendix 2, we expect most sequences to be contained in two main clusters, close to OS-A and OS-B’, respectively, with at least one sequence forming a separate cluster whenever the orthogroup was covered in the γ SAG. Based on this expectation, we considered all of the clusters from each orthogroup in decreasing order of their sizes and labeled them as described in the main text.

We next determined which orthogroups are part of the core genome. We defined a gene as a core gene if it is likely present in all of the genomes in our sample. The average coverage for each SAG in our sample was approximately 30%, so it is not possible to directly determine whether some genes are always present. Instead, we use a statistical measure of the coverage to assign orthogroups to the core genome. Specifically, we used the rank-coverage plot across all orthogroups and SAGs (equivalent to the complementary cumulative distribution) to determine an appropriate cutoff for the core genome. The rank-coverage plot showed three distinct regions. First, there were around a dozen orthogroups that had unusually high abundances, with the number of sequences in some cases being greater than the number of SAGs in our sample. These are likely multiple-copy genes and IS elements and were excluded from further analysis. Then there was a large plateau comprising over 1000 orthogroups with coverages in the range 100−150, followed by a sharp decrease (on a logarithmic scale) in the coverage. We found that setting a coverage threshold for the presence at 20% of the number of SAGs from each species approximately matched the shoulder between the plateau and the sharp decrease, and so we used this as a criterion for assigning core genes.

The final step in our algorithm is to catalog hybridization events, which we do by simply searching for gene sequences whose labels were different from the whole-genome species assignment. Analysis of these results showed a large number of transfers, mainly from β into α, which came from contigs that did not contain any genes labeled as the host species. One explanation for the presence of such contigs is that they represent larger genomic segments flanked by repetitive regions (such as IS elements) that are transferred and integrated into the host genome. Long repetitive genomic segments can lead to contig breaks during assembly, which would explain the complete absence of genes belonging to the host cluster. Alternatively, such contigs could be the result of contamination with external DNA. However, the contamination hypothesis does not explain why there seem to be more transfers from β into than vice versa. We also did not find any evidence of differences in the number of such transfers between the sample containing only $\alpha$ SAGs (OS2005, see Appendix 10) and samples containing mixtures of and β. While this evidence suggests these transfers are very likely the result of hybridization, we cannot rule out the possibility of contamination. We therefore excluded all hybridization events on contigs that did not contain any genes labeled as the host genome from the analysis. The results of our analysis are summarized in Figure 5A, B in the main text.

Appendix 5—figure 1

Download asset Open asset

In this section, we describe the method we used to discover SNP blocks and the statistical analyses we performed.

In order to identify linkage blocks originating from hybridization, the raw alignments from the pangenome construction were processed as follows. We divided core OGs into distinct species clusters, as explained in Appendices 1 and 5, and removed any sequences from SAGs assigned to a different species than the cluster. This was done to avoid false positive blocks due to the presence of recently transferred whole-gene alleles from different species. For OGs with mixed-species mosaic clusters, all sequences in the cluster were initially included. Alignments were then trimmed and codons that contained a high frequency of gaps were removed. Finally, we defined a main cloud of sequences for α and β comprising sequences within a divergence diameter $d_M=0.1$ from the consensus for each species and removed all sequences outside it.

We used a linkage-based approach to search for evidence of hybridization within genes. Visual inspection of the alignments in the highly diverse α genes and in most β genes revealed segments with high SNP density grouped into two haplotypes in strong linkage disequilibrium. We hypothesized that these SNP blocks were the result of intragenic hybridization, and we developed an algorithm to systematically identify them.

The algorithm can be summarized by the following steps:

1. Read orthogroup alignment

2. Pre-process alignment

   1. Trim alignment edges and remove columns with high fraction of gap codons

   2. For each species, infer consensus sequence and remove sequences outside a predetermined diameter around the consensus

   3. Replace all SNPs with absolute frequencies $n<n_c=4$ with the consensus for their respective columns

3. For each species $s$, subsample sequences from cells assigned to $s$

4. Calculate the linkage matrix $r_{ij}^2$ for all pairs of SNPs

5. Identify groups of consecutive SNPs with $r_{ij}^2=1$, up to a rounding error of $ϵ={10}^{-5}$, containing at least $l_{min}=5$ SNPs

The inclusion of step 2c accounts for the possibility of more recent low-frequency SNPs emerging after the initial hybridization took place. Note that typical alignments depths are ∼30 for α and ∼70 for β and are considerably larger than $n_c$.

We applied the above algorithm to the alignments of all of the core genes and obtained a total of 1000 blocks within $\alpha$ and 1483 blocks within β. Most blocks were short, but those in $\alpha$ were on average more than twice as long as those in β ($⟨L⟩=68$ in α vs $⟨L⟩=29$ in β; Appendix 6—figure 1A). Similarly, α blocks contained more SNPs on average compared to β, although here the difference was less pronounced ($⟨l⟩=11.3$ in α vs $⟨l⟩=7.9$ in β; Appendix 6—figure 1B).

Appendix 6—figure 1

Download asset Open asset

To identify hybrid haplotypes, we first determined the consensus sequences for the two block haplotypes, together with the other two species within the block segment. We then labeled any haplotypes whose sequence was identical to one of the other two species as a hybrid segment with the corresponding species as the donor. These results are shown in Appendix 6—table 1.

To control for other processes that could potentially produce similar patterns of linkage disequilibrium, we compared the distributions of divergences between haplotypes and block lengths for cases where a hybrid haplotype could be identified directly to cases we could not. No obvious difference in the distributions could be seen between the two cases, consistent with our interpretation that the blocks without matches to other species in our sample are also the result of hybridization.

Appendix 6—figure 2

Download asset Open asset

We further verified that the divergences between block haplotypes were consistent with what we would expect from hybridization. We calculated the fraction of nonsynonymous and synonymous nucleotide substitutions, $p_N$ and $p_S$, respectively, between the consensus sequences of each haplotype (see Appendix 7). We then removed blocks with either $p_N$ or $p_S$ were above a threshold of 0.6. This was done to avoid large errors in estimating the divergence, which we calculated using the Jukes–Cantor correction $d=-(3/4)\ln (1-4p/3)$ (Jukes and Cantor, 1969). The overall number of such blocks was small (46 out of 1442 in α and 143 out of 2008 in β) so their removal does not affect our conclusions about the typical blocks. We found that the distributions of both the nonsynonymous and synonymous divergences were similar to those previously reported from comparisons of OS-A and OS-B’ (Rosen et al., 2018), and consistent with our hybridization hypothesis (Appendix 6—figure 3).

Appendix 6—figure 3

Download asset Open asset

In this section, we present our analysis of the genetic diversity within individual genes in α and β. In the first subsection, we define the different metrics we used and present the variations in diversity along the genome of both species. In the next four subsections, we compare the diversity of different genes and use the α backbone to roughly infer the time scales of the evolutionary process. In the last subsection, we present a scenario for the evolution of the population.

We quantified the nucleotide diversity across all genes in both α and β using the following procedure. Alignments were trimmed for excess gaps and sequences were assigned to each species as described in Appendix 1. We then define the nucleotide diversity in the usual way as

where $I$ is an index for sequences in the alignment, $S$ is the sample size, equal to the number of elements in $I$, and ${\pi }_{ij}$ is the fraction of nucleotide substitutions between i and j, excluding gaps. To calculate synonymous (${\pi }_S$) and non-synonymous diversities (${\pi }_N$) we replace ${\pi }_{ij}$ with the fraction of synonymous (${\pi }_{ij}^S$) and non-synonymous substitutions (${\pi }_{ij}^N$), respectively. We calculated ${\pi }_{ij}^S$ and ${\pi }_{ij}^N$ using two different methods: the fraction of fourfold degenerate (4D) and onefold degenerate (1D) sites, respectively, at which i and j differ and using the Nei–Gojobori method (Nei and Gojobori, 1986). We found that both methods gave similar results and used the one that was most convenient for each analysis.

To calculate the divergence between species, we first determined the proportion of sites that were substituted as follows:

where $I_{\alpha }$ and $I_{\beta }$ are the indices of sequences from α and β in the alignment, and $S_{\alpha }$ and $S_{\beta }$ are number of sequences from each species, respectively. The fraction of pairwise substitutions ${\pi }_{ij}$ is the same as above. We then use the Jukes–Cantor correction (Jukes and Cantor, 1969) to correct for the effect of multiple substitutions at a single site

The synonymous ($d_S$) and non-synonymous divergences ($d_N$) were calculated in a similar way by replacing ${\pi }_{ij}$ with ${\pi }_{ij}^S$ and ${\pi }_{ij}^N$ in Equation 8, respectively. Finally, to calculate the diversity within each species separately, we take the sum in Equation 7 over $I_{\alpha }$ or $I_{\beta }$ only.

There are two plausible explanations for the heterogeneity within α found in Figure 3 (main text). The low-diversity genes could be due to local sweeps similar to the ones that formed the divergence troughs shown in Figure 5C in the main text. Under this scenario, the high-diversity genes would represent the typical diversity within α that has accumulated since its common ancestor. Since the high-diversity genes are as diverse as typical β genes, this scenario would imply α and β each had common ancestors at similar times in the past. Alternatively, the low-diversity genes could represent the typical diversity within α, with the high-diversity genes being the result of hybridization with other species.

To distinguish between these two possibilities, we developed a method to separate low diversity from high-diversity genes within α. Specifically, we calculated the fraction $f_l^S=P_l^S/L_l^S$ of 4D sites at locus l that were polymorphic within the α alignment. Here, $P_l^S$ is the number of synonymous polymorphisms and $L_l^S$ the number synonymous sites at that gene. By using 4D sites, we naturally control for differences in gene conservation. If mutation rates across loci do not vary too much, we expect the distribution of $f_l^S$ across genes to be given by a Poisson distribution. The mean of this Poisson distribution can be fitted from the fraction of sites that are not polymorphic

from which we get the following estimate for the mean:

When we compared this prediction to the data for α, we found that the Poisson distribution gave an excellent fit for around half of the core genes (675 out of 1448) with $f_l^S\lesssim 0.05$. In the remaining core genes (773 out of 1448), $f_l^S$ was much larger than predicted from the Poisson (Figure 3B in main text). Note that the agreement in the inset of Figure 3B is based on a single fitting parameter. Based on these results, we chose a threshold of $f_l^S\le 0.05$ to define the $\alpha$ ancestral backbone, with all other genes labeled as hybrid.

To validate our partition, calculate the synonymous diversity ${\pi }_S$ within each group. We found that our criterion produced a clear separation in the synonymous diversity, with low-diversity genes having a mean ${\pi }_S\approx 1.7\times {10}^{-3}$, compared to ${\pi }_S\approx 3.4\times {10}^{-2}$ for the high-diversity genes. As a result, we estimate that more than 95% of the genetic diversity within $\alpha$ is the result of hybridization. This result shows that the molecular clock estimate of the time since the most common ancestor based on $T_{MRCA}={\pi }_S/(2\mu \tau )$, where μ is the mutation rate per division and $\tau$ the generation time, should be interpreted with caution. In the case of $\alpha$ the estimate using the average ${\pi }_S$ across the whole genome would be around an order of magnitude less than the true value, even without considering how mutation rates might change over such long time scales.

We performed a similar analysis for β but did not find any large peak near $f_l^S=0$ (Figure 3C in main text). Instead, the distribution of $f_l^S$ was close to unimodal but with significantly fatter tails compared to the Poisson. Given the excellent fit to the Poisson for α in the region where SNPs are purely due to mutations, this suggests that other mechanisms play an important role in generating the diversity within β. The most parsimonious explanation is that the hybridization process we observed in α has covered most of the original β genome by this point. While it is not clear what the distribution of $f_l^S$ should be in this case, it is plausible that it need not be Poisson (Sakoparnig et al., 2021). Further evidence in favor of this hypothesis is the presence of SNP blocks across a large fraction of β core genes, as discussed in the main text and Appendix 6.

The large number of SNP blocks found throughout the genomes of both α and β shows that hybridization is an important source of genetic diversity for the Synechococcus population. However, because of the requirement for perfect linkage within SNP blocks, the total number of SNPs that are the result of hybridization is likely much higher than the one we infer from our algorithm. To obtain quantitative estimates of the total contribution from hybridization to the diversity, we use α as a control. As shown previously, only ≈5% of the total diversity in the high-diversity loci from α can be attributed to mutations. These loci contained 603,588 sites and 64,462 SNPs, of which 10,391 SNPs were found inside SNP blocks. Thus, only around 17% of the SNPs that are likely the result of hybridization, or hybrid SNPs, are included in the SNP blocks.

Within the β core genome, we have 1,591,476 sites and 146,906 SNPs, of which 14,819 were within SNP blocks. Assuming a similar ratio of hybrid SNPs to SNPs in SNP blocks as in $\alpha$, we estimate that around 90,000 or approximately 60% of SNPs in β are the result of hybridization. This is consistent with β having evolved over longer time compared to α. However, our knowledge of how multi-locus SNP patterns change over the course of evolution is limited, so the estimates for β should be interpreted with caution. In particular, if larger segments continue to be broken up due to recombination, it is likely that the fraction of hybrid SNPs within blocks would decrease over time for a constant minimum block size $l_{min}$ (see Appendix 6). Thus, the overall contribution of hybridization to the β genetic diversity may be even higher.

Combining genetic diversity analysis with evidence from the geological record of the Yellowstone caldera allows us to make inferences about the evolutionary history of the Synechococcus population. If we assume that most of the fixed variants between species are the result of mutation accumulation, then we can use the typical synonymous divergence between species to infer the time since the most recent common ancestor for the whole population. Typical bacterial mutation rates for natural populations are around $\mu \sim {10}^{-9.5}$ per bp per generation (Chen et al., 2022). Using this rate as a baseline and assuming generation times of $\tau \sim 1-10$ days gives an estimate of $T_{MRCA}\sim {10}^{6.5}-{10}^{7.5}$ years (Rosen et al., 2018).

For our argument, the lower bound is the more important one since it determines whether we can exclude the single-colonizer hypothesis. A $T_{MRCA}$ below 10^6.5 years would imply a shorter division time or higher mutation rate than our estimates. Our lower bound of 1 day for the division time is around the same as the fastest division rate observed in the laboratory. If the population is in steady state, as it appears to be, the division time sets the maximum generation time, which can be attained if the death rate is equally high. If the death rate is lower than this, other factors would limit cell division and the generation times would be longer. Taken together, we can be fairly confident that a generation time of 1 day is a hard lower bound.

The mutation rate is less certain since we do not have any reliable measurements and experimentally measured mutation rates might not be a good predictor of substitution rates over evolutionary times (Rocha, 2018). We can instead use other bacteria to estimate a range of values. In wild-type E. coli, mutation rates of $\approx 2\cdot {10}^{-10}$ per bp per generation have been measured in laboratory conditions (Lee et al., 2012), but mutator strains with rates that are up to 100 times higher are also known (Sniegowski et al., 1997). In natural populations, rates as high as $4\cdot {10}^{-9}$ have been inferred using longitudinal sequencing (Moran et al., 2009). If we use 10⁻⁸⁵ mutations per bp per generation as an upper bound, we would infer a lower bound of $T_{MRCA}\sim {10}^{5.5}$ years, which very narrowly overlaps with the time since the formation of the Yellowstone caldera. Thus, while it is not inconceivable that the three species could have diverged within the caldera, such a scenario would require all the relevant parameters to be at the limit of currently known values. In addition, such a scenario would require additional hypotheses to explain why the strains initially diverged into distinct species, with seemingly little genetic exchange between them (Rosen et al., 2018), and then only much later started to hybridize. Taken together, current evidence suggests the most likely scenario is that α, β, and γ diverged from a common ancestor before the formation of the caldera around 640,000 years ago.

Appendix 7—figure 1

Download asset Open asset

By using the diversity within $\alpha$ as a proxy for time, we can obtain estimates on the time scales of spatial mixing that influences the diversity within the population. To do this, we calculated the pairwise distances in the low-diversity regions of the $\alpha$ genome between pairs of cells in each spring separately and pairs of cells from different springs. We found that all three distributions were similar (Appendix 7—figure 1). The mean of the distribution for cells Mushroom Spring was slightly lower than those from the Octopus spring (${\pi }_S=8.6\cdot {10}^{-4}$ vs ${\pi }_S=1.1\cdot {10}^{-3}$) and this difference could not be explained by differences in sample sizes between the two springs. However, both values were very close to the mean of cells from different springs (${\pi }_S=1.0\cdot {10}^{-3}$). This suggests that the time scale of spatial processes relevant for the long-term evolution of the population, ${\tau }_m$ is much shorter time scales than $T_{MRCA}$ for α. Taking $T_{MRCA}\sim 3\cdot {10}^3$ years as a reference point sets an upper bound of ${\tau }_m\ll 3\cdot {10}^3$ years.

We propose the following scenario that is consistent with our previous discussion. First, we argued that a small number (3–6, but possibly more) of Synechococcus populations that had been evolving separately for $\sim {10}^6-{10}^7$ years independently colonized the Yellowstone caldera after its formation 640,000 years ago. After colonization, they started to hybridize, leading to a gradual mixing of their genomes (Figure 8). This can be most clearly seen in the α species, where the diversity in core genes that do not show signs of hybridization is more than an order of magnitude less than those that do. Moreover, genes that have undergone hybridization account for the majority of diversity within α. Within β, we do not observe such heterogeneity, but we find SNP blocks are even more prevalent than in α (Appendix 6—table 1). The simplest explanation for this pattern is that β has been undergoing hybridization for longer, and by now, all of its core genome has been mixed with other species. Consistent with this explanation, we find that only a small fraction of SNP blocks within β can be traced back to α (Appendix 6—table 1). Most of the β SNP blocks cannot be traced to either α or γ and are likely the result of a combination of hybridization with species that are not currently within the Lower Geyser Basin, have abundances lower than our detection threshold, or have gone extinct.

There are several important parameters that determine the impact of recombination and hybridization on the evolution of the population. First are the rates of recombination and mutation in the population. We denote the population recombination rate by $\rho =2TR$ and the population mutation rate by $\theta =2T\mu$, where T is the average time since the most common ancestor, $R$ is the recombination rate per site, and μ is the mutation rate per site. Note that in the panmitic Wright–Fisher model, T in units of generation time is equal to the population size $N$, but this is not generally true (Neher and Hallatschek, 2013; Weissman and Hallatschek, 2014; Birzu et al., 2018; Birzu et al., 2021). A second important parameter is typical length of recombined segments, which we denote by $L_r$. Finally, there is the typical divergence of recombined segments. In this Appendix, we use the nucleotide distance $p$ or the synonymous distance $p_S$, as defined in Appendix 7. We estimated these parameters for the hybridization process, which we argued previously is the dominant source of diversity in the population.

Before describing the analysis, it is important to note that while some of the simple statistics discussed here may be captured by simple neutral models, the underlying dynamics may be far more complex. Indeed, the analysis in the main text shows extensive evidence of selective sweeps affecting the genetic diversity across a broad range of genomic length scales, contrary to simple neutral theory. Moreover, while our analysis shows that—over time scales implied by our sequencing depth—the Synechococcus population is effectively well-mixed within the caldera, the actual population is obviously not. Synechococcus cells live in a dense microbial mat spread across thousands of springs separated by distances up to ∼10² km (Ward et al., 1998). As a result, we expect the long-term evolution of the population to be described by effective parameters, which may have little relation to recombination rates and population sizes directly measurable in the field or in the laboratory (Neher and Hallatschek, 2013; Weissman and Hallatschek, 2014; Birzu et al., 2018; Birzu et al., 2021). Understanding the relationship between the measurable parameters describing recombination, selection, and spatial dynamics of individual cells and the effective parameters at the population scale is an important and interesting topic for future work.

Many different methods have been proposed to estimate these or other closely related parameters (Croucher et al., 2015; Didelot and Wilson, 2015; Lin and Kussell, 2019). All of the methods we are aware of are based on models that make various assumptions that may not be valid for any given population, such as ours. In our case, an accurate estimate of these parameters requires a quantitative model of the hybridization process and is beyond the scope of the paper. Here, we instead used the results of our previous analysis to get rough estimates of these parameters and compared them to two other commonly used methods.

We considered each species separately, starting with α. For most evolutionary models, θ can be estimated from the nucleotide diversity π (Neher and Hallatschek, 2013). As in the main text, we used only synonymous sites to control for differences in the strength of purifying selection or the accumulation of non-synonymous mutations during rapid adaptation. Based on our previous analysis, we estimated the population synonymous mutation rate ${\theta }_S\approx 2\cdot {10}^{-3}$ (Appendix 7). Because of the clear distinction between the ancestral backbone and hybrid genes in α, we expect this estimate to be reasonably accurate. A rough estimate of $\rho$ can be obtained from the linkage disequilibrium curves (Figure 2 in main text). Our previous analysis showed that $\rho \approx 0.03$ and $\rho \approx 0.12$ gave good fits for α and β, respectively. Taking the ratio, we found a ratio between recombination and mutations of $\rho /{\theta }_S\approx 15$—much higher than previous estimates.

Estimating the mutation rate in β is more difficult because there are no genes unaffected by hybridization. We instead used the fraction of SNPs found in SNP blocks to estimate the diversity purely due to mutations. In Appendix 7, we estimated that around 40% of SNPs are the result of mutations. Multiplying with the synonymous diversity gave ${\theta }_S\approx 0.02$. Using this estimate, we found $\rho /{\theta }_S\approx 6$, which is lower than in $\alpha$ but still significantly higher than values that have been reported previously.

The other recombination parameters can also be inferred from our analysis. As we showed previously, the typical divergence of hybrid segments is given by the synonymous distance between species $p_S\approx 0.4$. We estimated the average size of recombined segments from the linkage curves. The recombination rate between two sites on the genome increases linearly with the separation x for distances shorter $L_r$. In the standard neutral model, this leads to the scaling ${\sigma }_d^2(x)\sim 1/x$ at intermediate distances (see Appendix 4). For $x>L_r$, the recombination rate saturates to the genomewide recombination rate $\sim \rho L_r$, which results in deviations of the linkage decrease from the predicted ${\sigma }_d^2(x)\sim 1/x$ scaling. In α, we observe a deviation around $x\sim 500$ bp. Note that we verified that this plateau cannot be explained by the smaller finite sample size in $\alpha$. A similar argument gives a shortened segment length $L_r\sim 300$ bp in beta. While we cannot directly test that the plateau is not due to sample size, the small effect of subsampling we observed earlier suggests it is an unlikely explanation (Appendix 4—figure 2). Note that in both species, our estimated genomewide recombination rate $\rho L_r$ is much greater than one, consistent with the genomewide linkage being close to the unlinked control, as discussed in the main text.

We compared our analysis to standard methods. ClonalFrameML and other related software are perhaps the most widely used methods to infer bacterial recombination (Didelot and Wilson, 2015). Such methods assume that recombination rates are rare enough that they can be described by a Poisson process on an asexual tree, similar to mutations. We refer to the assumption as the clonal frame hypothesis. ClonalFrameML uses the inferred phylogenetic tree from pairwise distances to identify discrete recombination events and assign them to branches on the tree. The recombination rate, genetic divergence, and length of recombined segments are then inferred from these events.

Note that the clonal frame hypothesis is inconsistent with the rapid decrease in linkage we observed (Figure 2 in main text). Nevertheless, it is useful to see how the method performs on our data.

We used a subset of high-coverage SAGs to create whole-genome alignments for α and β, which we then analyzed with ClonalFrameML. We first selected SAGs more than 75% complete, which resulted in a dataset of 2 α and 12 β SAGs (Appendix 1—figure 1). For each species, we chose the genes that mapped to OS-A or OS-B’ and created concatenated alignments of the entire genome in the order of the references. Missing genes or sequences were replaced by gaps. We constructed the phylogenetic tree using FastTree and ran ClonalFrameML on the alignments and trees with default parameters.

Appendix 8—figure 1

Download asset Open asset

For both species, ClonalFrameML inferred a ratio of recombination to mutation rates $\rho /\theta \approx 1$ ($R/\theta$ in Didelot and Wilson, 2015) and a fraction of nucleotide differences for imported segments of $p=0.07-0.09$ ($\nu$ in Didelot and Wilson, 2015). Note that both θ and $p$ reported by ClonalFrameML refer to nucleotide rather than synonymous differences. The values of $p$ are around half of the distance between species $p\approx 0.15$ and are inconsistent with our metagenome analysis which showed no genomes present in OS or MS at these divergences. The main difference between the species was in the average length of recombined segments, which was almost 10 times higher in α ($L_r=365$) compared to β ($L_r=45$).

A detailed comparison of the results of ClonalFrameML with our analysis is beyond the scope of this paper. Nevertheless, we found several inconsistencies in the assignment of recombined segments by ClonalFrameML. Because of the lack of data in α, we focused this analysis on β. First, we found that the distribution of recombined segment lengths was inconsistent with the exponential distribution assumed by ClonalFrameML and had a mean that was almost double the inferred $\delta =45$ (Appendix 8—figure 1). Second, manual validation of a small sample of randomly chosen recombined segments revealed identical patterns present in other cells that were not assigned to the recombination event (data not shown). We suspect that these misassignments are due to the inconsistency of the observed segments with the asexual tree.

The failure of ClonalFrameML to fit the data should not be surprising. As we showed, the low levels of linkage at whole-genome scales in α and β are inconsistent with the existence of a clonal frame assumed by the model. In addition, manual validation of some of the recombined segments detected revealed spurious segments caused by gaps in coverage due to the incomplete nature of the SAGs. These results caution against the naive application of methods such as ClonalFrameML to SAG datasets in general, but also to other populations for which the clonal frame hypothesis may not hold.

Appendix 8—figure 2

Download asset Open asset

An alternative method for fitting the clonal frame model was proposed recently by Lin and Kussell, 2017; Lin and Kussell, 2019. Instead of identifying recombination events directly, they use the correlations between SNPs to infer the recombination parameters statistically. Specifically, they model a local (sample) population that evolves according to a neutral model with recombination. They assume that the population recombines with a larger and more diverse ‘pool’ that evolves according to the same model but with different parameters. Under this model and assuming the recombination to mutation rate ratios are the same in the sample and pool populations, they show that the conditional probability that two sequences differ at a site given that they differ at a second site l bp away, $P(l)$, depends only on population mutation rate (θ in our notation), the population recombination rate ($\rho$) and the average length of recombined segments ($L_r$). By fitting the analytical expression for $P(l)$ from the model to the data, all of the model parameters can be inferred.

We generated alignments of the high-coverage α and β SAGs using Mauve (Darling et al., 2004) and followed Lin and Kussell, 2019 to fit $P(l)$. Note that unlike in ClonalFrameML, we used only synonymous sites for the mcorr analysis, so the divergences cannot be directly compared between the two methods. We found the model fit the data for $\alpha$ reasonably well up to two distances around 300 bp, with the exception of the first two points. However, the model was not able to capture the slower decay in correlations at longer genomic distances (see also Appendix 4). In β, the fit at short distances was considerably poorer (Appendix 8—figure 2). In particular, correlations between nearby SNPs were ∼2× higher in the data compared to the model. The increased correlations at short distances could be due to genetic hitchhiking during gene sweeps or could be a sign of subtle forms of population structure within β. Distinguishing between these scenarios is an interesting topic for future work.

Compared to ClonalFrameML, the recombined fragment lengths inferred by mcorr were larger (especially for β) and were reasonably close to our estimate. In addition, the mcorr estimate of the diversity due to mutations (θ in our notation) agreed remarkably well with our estimates in both species. Similar to ClonalFrameML, the synonymous divergence of hybrid segments was between the divergence between species and the diversity within species and was inconsistent with our metagenome analysis. Finally, the ratio of recombination and mutation rates was similar to the one inferred from ClonalFrameML and more than an order of magnitude smaller than our estimates. These results are summarized in Appendix 8—table 1. We leave a more careful analysis of the differences as a topic for future work.

In this section, we describe how we quantified associations between genomic features and populations from different samples. Our primary focus was on differences in hybridization patterns, so we analyzed the variations in hybrid gene and SNP block haplotype frequencies across different samples and locations. Given the small number of samples in our study, we limited the number of statistical tests we performed in order to avoid false positives. For each of the two main species, α and β, we tested for the effect of sample location and sample composition separately, as described below.

We considered two distinct hypotheses that could lead to associations between block haplotypes and environmental conditions from which the sample was taken. First, we considered if the community composition was associated with specific haplotypes by comparing samples in which a single species was present to samples which contained a mixture of $\alpha$ and β. For α, we chose the 48 cells from pure-α sample OS2005 as the test group and the 45 α cells from MS2004, MS2006, and OS2009 as the control. For β, we chose the 44 cells from the pure-β sample MS2005 as the test group and the 46 β cells from MS2004 and MS2006 as the control. We denote these two tests as ‘α-composition’ and ‘β-composition’ tests, respectively. Second, we considered the hypothesis that sample location was associated with specific haplotypes. For this test, we divided all cells from each species into two groups labeled MS and OS, based on the location from which the sample was taken. These tests are denoted by ‘α-location’ and ‘β-location’, respectively.

We used the following permutation test to find block haplotypes that were significantly associated with each environmental variable. We first removed all blocks in which the minority haplotype was below an arbitrary cutoff frequency, set at 5 sequences. This was done in order to avoid tests which were statistically underpowered. We then calculated the two-sided p-value for the association for each block individually using Fisher’s exact test.

Because the sampling depth for each block is variable due to variations in coverage, standard methods that correct for the false-discovery rate are not well-suited in our case. We use the following procedure to determine significant associations. We randomized the composition of the test and control cell groups for each test, while keeping the total number within each group fixed. We then recalculated the p-values using the same procedure as before and recorded the lowest p-value across all blocks $p_{rand}^{min}$. We used the ratio $p_{rand}^{min}/p$ between the minimum p-value in the randomized control and the p-value in the test group as a measure of significance of the association. We arbitrarily set a threshold value of $p_{rand}^{min}/p_t=10$, above which associations were deemed highly significant.

The results of these tests are summarized in Appendix 9—table 1.

Overall, we found very few blocks that were associated with either the sample composition or location. This is consistent with the results from the main text (Figure 7D), which showed that the block haplotypes are close to unlinked even on distances comparable to the whole genome. The largest number of associations was with the sample composition in the $\alpha$ species, consistent with the previously mentioned results, which showed a small but statistically significant increase in linkage within the α species compared to the β. Nevertheless, the overall fraction of blocks whose haplotypes were strongly associated with the sample composition was ≈1.2% and still small. We therefore concluded that neither the composition of the population within samples nor the location has a large effect on the hybridization patterns revealed by the SNP blocks. A complete table with the results of these tests can be found in the supplementary data tables.

We also examined variations in the rate of transfers of whole-gene alleles across species. The lower overall number of such events compared to the SNP blocks makes statistical tests for associations much more likely to yield false positives in this case. As a result, we limited the analysis to a qualitative comparison of the number of different types of transfers across samples shown in Appendix 9—figure 1. We did not observe any systematic variations with either sample location or sample composition for either species.

Appendix 9—figure 1

Download asset Open asset

Beyond the processing steps outlined in Appendix 1, we also explicitly verified that the hybridization patterns we observe cannot be explained by sequencing and genome assembly artifacts. Specifically, we tested if misassembly of reads from distinct alleles of the same gene could explain the presence of hybrid sequences. We considered two main sources of such reads: cross-contamination from neighboring reaction wells during whole-genome amplification (WGA) and the presence of paralogs horizontally transferred from other species.

Distinguishing between hybridization and artifacts is hard because both can lead to similar signals. Here, the structure of data provides us natural controls. First, α has genes with low diversity and little sign of hybridization and high-diversity genes with extensive hybridization. We can use former to establish a baseline for genome coverage and check for deviations in the later. If the hybridized genes are caused by misassemblies of actual between α genes and contamination from β or between different alleles of the same gene within α cells, we should see an increase in the coverage at these loci. Second, the variability in species composition across samples provides controls for cross-contamination between reaction wells. Specifically, OS2005 contained 48 cells, all classified as α, while MS2004, MS2006, and OS2009 together contained 45 α cells amplified on plates with mixtures of α and β cells in similar proportions (Appendix 10—figure 1). We will refer to the former as the control SAGs sample and the latter as the test SAGs sample. We used these two groups as a measure of the contribution of cross-contamination from neighboring β cells to the observed hybridization in α.

Appendix 10—figure 1

Download asset Open asset

We first divided the core α OGs into test and control groups as follows. We first divided OGs into low- and high-diversity groups based on the probability that synonymous sites are segregating in the sample using a cutoff at 0.05. From the low-diversity OGs, we chose the subset that were grouped into α and β, or α, β, and $\gamma$ clusters during the pagenome construction (approximately 300 each out of a total of ∼800). These OGs were nearly clonal in the α sample and showed little to no evidence of hybridization in our entire data. These OGs will form the control OG group. From the high-diversity OGs, we selected those where all sequences were grouped together in a single mixed-species cluster (∼300 out of a total of ∼600). These include the cases with the most extensive hybridization and thus formed the test OG group.

We validated that core gene coverage can be described by a random process in which each gene has a fixed probability $p$ of being present in the assembly. We fitted the baseline probability $p$ using the low-diversity genes from the OS2005 α cells. Many of these genes were clonal within the entire α sample and showed no evidence of hybridization with the other species in our sample. They thus provide a good baseline for the expected coverage in the absence of contamination. While the coverage among different SAGs was highly variable (Appendix 10—figure 2A), we found that the distribution of coverages of different genes was well fitted by a Gaussian distribution with the mean and variance given by assuming each gene is an independent Bernoulli random variable with success probability $p\approx 0.29$. We tested the independence assumption by calculating the covariance between the presence of different genes (Appendix 10—figure 2C). Under our very simple model, any correlations between the coverage at different loci would be purely due to statistical fluctuations caused by finite sample sizes. The distribution of diagonal ($C_{aa}$) and off-diagonal ($C_{ab}$) components of the covariance matrix can then be calculated exactly given $p$. We compared the values obtained from the data with these predictions and found a very good agreement (Appendix 10—figure 2D).

Appendix 10—figure 2

Download asset Open asset

Having validated our gene presence model, we can now use it to perform simple but quite powerful tests for different artifacts. Consider the extreme case where all of the hybrid genes we observe are the result of misassemblies of reads from two distinct alleles, one from $\alpha$ and one from either β or γ, from the same reaction well. Note that for this argument, it is unimportant whether the alternative allele is due to contamination or genuine horizontal transfer that resulted in two distinct, unmixed copies of the gene. In either case, the probability of observing the gene in our assembly would be $P_{obs}=p(2-p)\approx 0.5$, which is significantly higher than for a true single-copy gene. We tested for such an increase by comparing the distribution of coverages among α cells in the OS2005 sample of high-diversity genes that formed a single mixed cluster across all sequences in our data (see Appendix 1) and low-diversity genes with distinct α and β or α, β, and γ clusters. We found the two distributions were very similar to each other (Appendix 10—figure 3A), with a small but statistically significant increase in the mean coverage at low-diversity genes with distinct species clusters ($p_{control}\approx 0.3$ vs $p_{test}\approx 0.28$). The difference is in the opposite direction from the one expected from the misassembly hypothesis, but it corresponds to less than an extra copy of a gene on average and is unlikely to have any consequence for our conclusions.

Note that the misassembly hypothesis above would also predict a specific pattern for the sequences within α, with a large fraction $2(1-p)/(2-p)\approx 0.83$ having either the pure α or the pure β alleles, and $p/(2-p)\approx 0.17$ containing a mixture. We do not find evidence of such patterns in our data. More importantly, the 0.17 fraction is an upper bound on the maximum frequency of mixtures that can be explained through misassembly, as any lower amounts of contamination would lead to fewer chimeric sequences. Note that both of the above tests take advantage of the fact that the typical genome coverage is low.

Using a similar method to the one described above, we can also check for contributions from cross-contamination between reaction wells on the formation of chimeric sequences. Here, we compared the coverage distribution of mixed-species orthogroups in the mixed $\alpha -\beta$ samples described previously with the control group from the all α sample OS2005. If cross-contamination with β was a significant contributor to the presence of mixed-species orthogroups, one would expect a large difference between the two distributions. We did observe a large systematic difference between the two groups in the overall coverage $p$ across all genes (data not shown), which was likely due to batch effects in the amplification and sequencing. We therefore compared the mean-centered distributions between the two groups and found no difference between them and the fit to the independent presence-absence model from before (Appendix 10—figure 3B).

Appendix 10—figure 3

Download asset Open asset

Beyond the indirect presence of hybridization given by the presence of mixed-species orthogroups, we also tested for any effect of cross-contamination on the presence of linkage blocks. We applied the same algorithm described in Appendix 6 to the subsample of 48 α cells from OS2005 and compared the results with those from random samples of α cells from samples with mixtures of α and β cells. We found that the number of blocks detected in OS2005 was similar to the bootstrapped distribution from mixed $\alpha -\beta$ samples (Appendix 10—figure 4A). We also found the distribution of block lengths and divergences between block haplotypes to be similar as well (Appendix 10—figure 4BC).

Appendix 10—figure 4

Download asset Open asset

Based on the evidence presented above, we concluded that cross-contamination between reaction wells and genome assembly artifacts does not contribute significantly to the presence of mixed-species orthogroups or linkage blocks. While we cannot definitively exclude that some sequences we observe are the result of amplification or assembly errors, the evidence above shows that they are rare enough to not have a significant effect on our main conclusions.

The funders had no role in study design, data collection, and interpretation, or the decision to submit the work for publication.

We thank Benjamin H Good, Alana Papula, and Freddy Bunbury for helpful comments and discussions. GB and DSF gratefully acknowledge support from the Simons Foundation via a Postdoctoral Fellowship Award 730295 to GB and Sabbatical Fellowship to DSF. This work was also supported by NSF Grants PHY-1607606 and PHY-2210386 to DSF. HSM was supported by the National Institutes of Health grant R01-AI-100947 to Mihai Pop. DB acknowledges support from the NSF/BIO-BBSRC collaborative research grant 1921429, NSF MIM grant 2125965, Joint Genome Institute Community Sequencing Project 503441 and 509352, and the Carnegie Institution for Science. Samples processed for single-cell genome amplification were collected using NP Park Permits: YELL-5494 to David Ward (multi-year), YELL-5660 to DB (2007–2008), and YELL 5694 to DB (2007–2009). The work (proposal: https://doi.org/10.46936/10.25585/60001132) conducted by the U.S. Department of Energy Joint Genome Institute (https://ror.org/04xm1d337), a DOE Office of Science User Facility, is supported by the Office of Science of the U.S. Department of Energy operated under Contract No. DE-AC02-05CH11231.

1. Preprint posted: June 7, 2023
2. Sent for peer review: August 11, 2023
3. Reviewed Preprint version 1: November 7, 2023
4. Reviewed Preprint version 2: October 21, 2025
5. Version of Record published: December 31, 2025