Genetic Exchange Across a Species Boundary in the Archaeal Genus Ferroplasma

Study location:

Genomic sequence data were obtained from DNA extracted from a pink biofilm floating on pH 0.7 acid mine drainage and sampled from the five-way location within the Richmond Mine, near Redding, California (supplemental Figure S1 at http://www.genetics.org/supplemental/). As described previously (Tyson et al. 2004), the biofilm community was dominated by bacteria (Leptospirillum groups II and III) and four archaeal organisms from the Thermoplasmatales lineage.

Community genomic data:

Approximately 130 Mbp of shotgun sequence data (3-kb library) from DNA extracted from the biofilm were assembled using the phredPhrap program (Ewing and Green 1998). Assemblies were comparable to those generated using the first 76 Mb of sequence data assembled with the JAZZ program, as described by Tyson et al. (2004). Sequences in the expanded data set were grouped by organism type, as described for the initial 76 Mb of sequence data (Tyson et al. 2004). As noted previously (Tyson et al. 2004), assembly of all sequencing reads from the community generates genome fragments (scaffolds) from both Ferroplasma type I [similar to the Ferroplasma acidarmanus isolate (Edwards et al. 2000)] genome and Ferroplasma type II. These assembled Ferroplasma genome fragments were manually curated to resolve errors arising from repetitive elements. Some gaps were closed by allowing for insertion/deletion of single or small groups of genes in some, but not all individuals (strain variants).

Visualization of variability, cluster analysis, and identification of recombination:

Composite Ferroplasma types I and II genome fragments were assembled from reads from coexisting individuals with similar or identical sequence types and do not capture population-level heterogeneity. To understand the degree to which sequences from these organisms form clusters, it is necessary to determine the extent to which their sequencing reads form a continuum. Genomewide comparison of sequencing reads to each other and to both composite sequences provided a means to investigate the extent of the continuum.

All reads from the community genomic data set were trimmed to remove low-quality bases (below a Phred score of 20) and then aligned to all community scaffolds using BLASTN with a cutoff e-value of 10⁻⁷⁵. Reads better assigned to composite sequences of organisms other than Ferroplasma types I and II were excluded from subsequent analyses. For analysis of clustering, all Ferroplasma types I and II reads were aligned to the F. acidarmanus isolate genome (Edwards et al. 2000) using a very low cutoff e-value (10⁻²⁵). Some reads and parts of reads were not brought into this analysis, often because they are associated with genes not present in F. acidarmanus.

In Figure 1A, the BLASTN output has been converted into a graphical representation showing Ferroplasma type I and II sequence reads (white bars) aligned to a reference genome fragment. Colored ticks within each read represent single-nucleotide polymorphisms (SNPs) relative to the reference sequence. Cyan, pink, magenta, and yellow represent substitutions of the bases A, C, T, and G, respectively. Thinner white lines connect a read to its mate pair (the sequence from the opposite end of the same clone).

Figure 1B illustrates a simpler rendition of data, used to examine larger genomic regions, than the one shown in Figure 1A (also see Figures 2–4). Background colors (shades of brown and violet) are used to highlight the separation of sequences into clusters based on SNP frequency (Figure 1B). For analysis of variation within these clusters, we compared reads with few SNPs (brown background in Figures 1–4 and in supplemental Figure S2 at http://www.genetics.org/supplemental/) to the very similar F. acidarmanus genome. The second cluster (violet background in Figures 1–4 and in supplemental Figure S2) has a composite sequence similar to that of the Ferroplasma type II composite genome. Consequently, for analysis of variation, these reads were aligned to the Ferroplasma type II composite genome fragments. Within clusters corresponding to the Ferroplasma type I and II populations, distinct variant types could be reconstructed on the basis of patterns of SNPs. Variants within populations are illuminated using slightly different background colors (Figure 1B).

Some reads show discrepancies in linkage patterns. For example, the read outlined in red in Figure 1B can be split into two regions where, on the basis of SNP patterns, each region belongs to a different variant type. These mosaic patterns of SNPs, which can be viewed as a skewed distribution of polymorphic sites, are evidence of past recombination (Sawyer 1989). A small black box indicates cases where the recombination point occurs within a read (the transition between sequence types within the read is evident from the pattern of SNPs) (Smith 1999). When reads of a mate pair have different strain types, a recombination point is inferred to be present in the unsequenced part of the clone.

Sequence types related by a recombination event can be identified and reconstructed from overlapping reads (Figure 1C). For the purpose of calculating the nucleotide divergence between sequences related by a recombination event, those sequences that represent the dominant linkage pattern types were defined as “parental” (Figure 1, B and C).

Despite the general separation of sequences into the two distinct Ferroplasma types, regions exist where linkage patterns inferred from sequence similarity indicate recombination between Ferroplasma types I and II. To quantify this form of genetic exchange genomewide, we identified recombination points by investigating the BLASTN (Altschul et al. 1990) alignments of the clones that could not be clearly assigned to one or the other population.

Estimation of recombination frequency as a function of sequence divergence:

We quantitatively evaluated the dependence of recombination frequency on sequence divergence. For each clone showing intra- or interpopulation genetic exchange, we calculated, for the span of the recombinant clone, the nucleotide sequence divergence between the reconstructed parental sequences (as shown in Figure 1C). The full length of the recombinant clone was considered because, since only one recombination point is visible, there is no clear indication as to which end of the clone contains the introduced sequence. In some cases, the comparison region was reduced because one or both of the parental strain fragments did not span the full length of the recombinant clone.

Recombinant clones were counted in increments of 1% sequence divergence for intrapopulation data and 3% for interpopulation data. However, the frequency with which recombination events will be detected at any value of sequence divergence will be determined in part by how common that level of sequence divergence is across the genomes. Specifically, few recombination events are expected at high divergence due to a small fraction of the genomes showing these high divergence levels. Consequently, the numbers of recombination events as a function of sequence divergence (recombination frequencies for defined sequence divergence intervals) were normalized to account for the nonuniform distribution of sequence divergence. The distribution of sequence divergence was determined from sequence comparisons involving all reads.

Relatively few recombination events can be detected below 2% divergence because sequencing errors and recent mutations obscure already subtle differences between strain types. Due to the rarity of highly divergent variants, significant noise levels are anticipated for these cases. For these reasons, we restricted our analyses to divergences between 3 and 11%. Results were fit to an exponential model using maximum-likelihood estimation in the program R (Nelder and Mead 1965; R Development Core Team 2006). Errors in raw recombinant counts were assumed to be normal and were estimated by the maximum-likelihood optimization process.

Quantifying rates of recombination relative to mutation:

The ratio of recombination frequency to mutation rate (ρ/θ) gives an indication of how much more or less likely it is that any nucleotide will change due to recombination than due to mutation (Feil et al. 1999). Because our data differ from those published previously, we derived a new method of calculating this probability. The number of recombinant clones identified (as described above) was used as a proxy for the number of recombination events in the data set. When foreign sequence replaces existing DNA in an organism, only a fraction of positions involved will change (determined by sequence similarity, 2% for intrapopulation recombination events and 18% for interpopulation events). The observed recombination event frequency, r₀, is multiplied by the estimated average recombination block size, h, and the average divergence between sequence types (as a proportion of nucleotides that are different), d, to estimate how often a single base was changed by recombination. The product r₀dh is divided by the estimated mutation rate, μ₀ (the average nucleotide divergence within sequence types), to give the relative frequency of recombination to mutation:

(1)

We used Watterson's segregation site-based calculation of the mutation rate to estimate θ (Watterson 1975). This method assumes no recombination, so the calculation was performed on each strain and averaged over the entire data set.

If interpopulation recombination events are rare and widely separated around the genome, the recombination block size can be determined so long as the event is evident in the majority of cells sampled and the sequence transition is captured at both ends of the block. The average block size for interpopulation recombination was estimated from these observations.

Intrapopulation, determination of the recombinant block size may be more problematic because the blocks are likely to overlap. We used several models that made use of the interpopulation observations and the frequency with which two recombination end points were detected in a single clone to estimate this value.

Given that a block of size x is detected, the probability that it lies entirely within a clone of size C is

(2)

We know the average clone size, C, to be 3.3 kb. Multiplying the above probability by the probability of having a recombination block of size x and summing that product over all possible block sizes gives the expected value of r. Working backward we can estimate the average recombination block size for a given (assumed) distribution of recombination block sizes. In model 1, we performed this calculation assuming a uniform distribution of sizes ranging from 0 to twice the average block size.

In model 2, we calculated the average block size by assuming a normal distribution of block sizes. Eight direct observations of interpopulation recombination block sizes were resampled with replacement 1 million times. Intrapopulation average block size was then estimated for all variances found in these bootstrapped data.

In model 3, we calculate the average length of sequence per observed recombination event. We compare the estimates of all three models with the block size observed for interpopulation events to estimate the likely range of average block sizes for intrapopulation recombination.

Fluorescence in situ hybridization:

Specific oligonucleotide probes were designed to target the 23S rRNA in Ferroplasma type I (fer1_23-850, CCTCACTAGACATCTCC) and Ferroplasma type II (fer2_23-242, CGATCGCCTTTAATCACGG) as previously described (Hugenholtz et al. 2002) and took into consideration accessibility of the target region (Fuchs et al. 2001). Probes were commercially synthesized and 5′ labeled with the fluorochrome fluorescein isothiocyanate (FITC), CY3, or CY5 (Operon, Alameda, CA). Paraformaldehyde-fixed samples were analyzed with both Ferroplasma probes and used in combination with general, bacterial (EUB338mix), and archaeal (ARC915) probes at optimal hybridization and wash conditions, to qualitatively determine the abundance of Ferroplasma type I and Ferroplasma type II. The optimal stringency for the probes was determined empirically using fixed AMD samples and 5% formamide increments from 10 to 50%. F. acidarmanus was used as a positive control for fer1_23-850 and a negative control for fer2_23-242. A Leica DM RX microscope was used for visualization of FISH preparations.

Comparative genomics of composite genome sequences:

On the basis of the number of clones analyzed and an estimate of the number of cells sampled, we obtained one clone for every 1000–10,000 Ferroplasma types I and II individuals present in the sample. A total of 24,184 clones were assigned to Ferroplasma type II (sampling depth of ∼22×) and 6820 clones were assigned to Ferroplasma type I (sampling depth of ∼6×). An additional 1027 clones could not be uniquely assigned to either population. Assembled contiguous DNA fragments (contigs) assigned to Ferroplasma type II were further assembled into a composite genome with only a small number of internal gaps. Assembly made use of mate pair information (i.e., through placement of paired sequences from the ends of clones on separate contigs; supplemental Table S1 at http://www.genetics.org/supplemental/) and accommodated the effects of heterogeneity in gene content within populations. In the refined composite Ferroplasma type II genome, the previously reported JAZZ assembly fragments have been ordered and linked (Tyson et al. 2004). Some small fragments representing strain variants are not included in the composite sequence and were not included in most analyses. The reconstructed scaffolds assigned to Ferroplasma type I contain a 16S rRNA gene identical to the gene in the F. acidarmanus isolate genome (Edwards et al. 2000). These scaffolds cover 93% of the isolate genome and share 98% nucleotide-level sequence identity. For all subsequent analyses, the F. acidarmanus isolate genome was used in place of the composite sequences as the reference for alignments.

The 16S rRNA gene sequence of Ferroplasma type II shares 99.2% identity with Ferroplasma type I (Tyson et al. 2004). The Ferroplasma type II genome encodes 1979 open reading frames (ORFs), compared to 1971 ORFs in the F. acidarmanus genome (Allen et al. 2007). The genomes share 1343 orthologs with an average ∼83% sequence identity. Gene order of the two composite genomes is largely conserved, with an average block size of ∼8 genes (see supplemental Table S1). There are 649 genes in the Ferroplasma type II genome absent in the F. acidarmanus genome compared to 610 in the F. acidarmanus (and a similar number in Ferroplasma type I) genome not present in the Ferroplasma type II genome (Edwards et al. 2000).

Clustering of sequences and evidence of recombination:

Sorting of reads based on sequence variation confirmed the division of Ferroplasma sequence into two genomic clusters that are generally represented by the original Ferroplasma type II JAZZ scaffolds and F. acidarmanus genome. Figure 2A (BLASTN e-value <10⁻²⁵) clearly illustrates this separation across a 3-kb window. Comprehensive analysis showed that these clusters of reads occur genomewide (Figure 2B), with an average nucleotide identity of 82% between Ferroplasma types I and II clusters and 98% identity within each Ferroplasma cluster.

Despite the existence of clear genomic clusters, the data indicate ongoing genetic exchange between Ferroplasma types I and II (e.g., Figure 3 and supplemental Figure S2 at http://www.genetics.org/supplemental/). Of the 1027 clones that could not be uniquely assigned to one or the other type, 709 derive from three regions where the sequence type is essentially the same in all sampled individuals in both populations. In two of these regions, the single sequence type suggests recombination was followed by a sweep that selected for only those genome types with the adaptive gene variant. The other region, containing an integrase and transposases, may be the result of phage insertion into both genomes. Detailed inspection confirmed that a further 171 clones show strong evidence of recombination between the two Ferroplasma species (e.g., Figure 3 and supplemental Figure S2) at 70 locations distributed across the genomes. Twenty-one of the identified locations were associated with transposable elements. In these cases, multiple reads align to the same part of a transposable element in one genome, but their mate pairs align to unassociated genes spread randomly throughout the other genome. Clones aligning to transposable elements were removed from consideration, leaving 139 clones showing evidence of 49 recombination events. Inclusion of the three identical regions brings the total observed recombination events to 52.

Within the two Ferroplasma populations we observed finer-scale clustering to form strain variant groups. Reads within groups share >99.5% nucleotide-level sequence identity whereas strain variants typically share ∼98% sequence identity. There is evidence of extensive recombination between strain variant sequence types. Every few thousand bases across the genomes, a single clone is found to transition from one variant type to another, indicating a population characterized by many mosaic genome types (Figure 4 and Tyson et al. 2004).

By analysis of all sequencing reads assigned to Ferroplasma type I, we identified 1293 (19%) clones showing recombination between strain variants. For Ferroplasma type II we identified 1278 (24%) recombinant clones in reads assigned to seven large genome fragments derived from around the genome (a total of 564,543 bp, sampling approximately the same number of reads as for the lower coverage Ferroplasma type I data set). Because recombination events become more difficult to identify as parental sequence types become more similar, the numbers of recombination events identified in the data sets for each population are minimums.

Decline in recombination with sequence divergence:

Intrapopulation recombination frequencies for both Ferroplasma types showed a clear log-linear relationship with sequence divergence (Figure 5). The coefficient of this relationship for Ferroplasma type I is 29 and for Ferroplasma type II is 44. Between Ferroplasma types, the relationship between the rate of recombination and sequence divergence is less clear due, at least in part, to the small number of recombination events found between populations. However, the data are consistent with log-linear dependence (Figure 5C).

Recombination rates relative to mutation:

The rate of recombination relative to mutation is a useful measure of the relative strengths of two of the major forces shaping the Ferroplasma populations. To calculate this rate, an estimate of the average length of transferred sequence, h, is needed. For interpopulation recombination events, we calculated h to be 9 kb on the basis of direct observation of eight recombination blocks (e.g., Figure 3 and supplemental Figure S2 at http://www.genetics.org/supplemental/). For intrapopulation recombination a direct measurement of h was not possible. We observed that ∼1 in 16 clones contained more than one recombination end point and used this value to estimate h, assuming a uniform (model 1) and a normal (model 2) distribution of sizes. We also estimated the recombination block size using the number of recombination points found relative to the total number of sequenced bases (model 3). In conclusion, we esti mate h to be in the range of 5000–10,000 kb (supplemental Figure S2).

Using the number of recombination events (recombination points/2) per base (5.4 × 10⁻⁵ for Ferroplasma type I and 5.1 × 10⁻⁵ for Ferroplasma type II), an estimate of the mutation rate (0.006 for both populations), and the above range of average recombination lengths, we infer that the probability that a base is changed by intrapopulation recombination rather than by mutation is between 2:1 and 4:1. On the basis of the 52 observed interpopulation recombination events, we determined that the relative frequency of base change within Ferroplasma type I or II due to interpopulation recombination compared to mutation is ∼1:25.

Locations of interpopulation recombination events:

We investigated the distribution of interspecies recombination points in both Ferroplasma genomes. A total of 77% of the events clustered within 500 genes of the origin of replication, where gene order is most conserved (Figure 6). This result suggests that breakdown of synteny further from the origin contributes to the reduction of homologous recombination (Kowalczykowski 2000). Interpopulation recombination events were often associated with regions of local high sequence similarity, such as shown in Figure 3B (vertical shading; also see supplemental Figure S3 at http://www.genetics.org/supplemental/), or tRNAs (Figure 6), supporting the conclusion that short regions of sequence similarity are a critical factor in initiating homologous recombination (Majewski and Cohan 1999). Five clones show evidence for interpopulation recombination within the 23S rRNA gene, a finding that highlights the potential for recombination to skew phylogenetic inferences made from such genes.

Community genomics and population genetics:

Our study demonstrates that unique insights into population genetics of natural microbial systems can be provided by deeply sampled metagenomic data sets. Specifically, analysis of sequence clustering and the frequency of recombination within and between clusters helps to resolve whether Ferroplasma types I and II represent different species. Additionally, the dependence of recombination on sequence identity is found to be the same as in prior experimental studies.

Are Ferroplasma types I and II different microbial species?

Genomewide analysis of sequence heterogeneity reveals two generally discrete clusters, approximately represented by the assembled composite genome sequences of Ferroplasma types I and II. These clusters have only an 87% average nucleotide identity (ANI) genomewide, significantly smaller than the 94% ANI found to correspond with the 70% DNA–DNA hybridization threshold used as part of the standard bacterial species definition (Konstantinidis and Tiedje 2005). The significant difference in metabolic potential (>600 genes are unique to each Ferroplasma type) may indicate adaptation to different niches (Coleman et al. 2006). Direct observation of phenotypic differences is not possible because Ferroplasma type II has not yet been obtained in pure culture. However, fluorescence in situ hybridization using lineage-specific probes for both Ferroplasma types suggests that the relative abundance of each can vary considerably and independently in samples from different environments (data not shown). On the basis of these findings—sequence clustering, the different gene content of these clusters, and evidence suggesting different environmental selection for the two organism types—we infer that these clusters effectively correspond to two species populations.

A general rule for recombination rates:

Our results demonstrate that the rate of genetic exchange in coexisting archaea in a natural microbial community has a log-linear dependence on sequence similarity, genomewide. This finding is consistent with experimental data for a single locus studied in isolated bacteria (Roberts and Cohan 1993; Zawadzki et al. 1995). Thus, as inferred for some isolated bacteria (Roberts and Cohan 1993), the recombination rate in archaea may also be primarily determined by the stability of the heteroduplex complex. Consequently, the log-linear relationship is far more general than previously established and may thus control on the distribution of diversity within and between population clusters.

When divergence within populations becomes sufficiently large, it can be inferred from the general log-linear dependence of recombination rates on sequence divergence that recombination rates will no longer be fast enough to spread newly acquired functions to all members prior to a selection event. In this case, it is tempting to suggest that the ancestral Ferroplasma population could have sympatrically diverged into separate species.

Fraser et al. (2007) recently used simulations to suggest that, under a model free of selection, such a scenario is possible only for organisms where the relationship between recombination and divergence is very steep (coefficient of ∼300). Although the coefficients found in this study (29 and 44) are higher than most of those previously reported (∼8–25) (Roberts and Cohan 1993; Vulic et al. 1997; Majewski et al. 2000), they are far short of this requirement. Thus, in the absence of a more sophisticated model incorporating the effects of selection, the divergence of these species may have involved some degree of physical separation.

It is perhaps not a coincidence that rates for interpopulation recombination involving the most similar sequences are lower than rates for intrapopulation events with the same sequence divergence (Figure 5D, 7–10% nucleotide divergence). These lower rates do suggest the presence of external barriers to genetic exchange. The mechanism of genetic exchange in extremely acidophilic natural archaeal populations is unknown. Uptake of free DNA from the environment is likely inhibited by its acid-promoted degradation. Genetic exchange may involve formation of cytoplasmic bridges, as described for Haloferax volcanii (Tchelet and Mevarech 1989), a conjugative mechanism documented in Sulfolobus species (Schleper et al. 1995; Grogan 1996), or vesicle-mediated genetic transfer (lipofection) (Felgner et al. 1987). For all of these mechanisms, increased physical separation due to partitioning of Ferroplasma types I and II into different microniches within the same biofilm would decrease the incidence of recombination.

The effect of short regions of similarity:

Most interspecies recombination events contained or were located near short regions of high sequence similarity that may affect the analysis of the relationship between recombination rate and sequence divergence. These short sequences are significantly shorter than the regions over which sequence similarity was measured (the length of the clone, usually ∼3 kb, see above). As a result, the measured sequence similarity will be much lower than in the short similar regions that may have contributed to the initiation of recombination. The lower slope (∼7) for the relationship between inter- vs. intrapopulation recombination and sequence divergence (Figure 5) may partly reflect this effect.

High recombination rates:

Our estimate that recombination is two to four times more likely than mutation to alter a single nucleotide is consistent with some recent observations finding higher rates of recombination than once thought, particularly among archaea (Papke et al. 2004; Whitaker et al. 2005). However, high recombination rates do not appear to be universal in archaea (Hallam et al. 2006). Thus frequent recombination within Ferroplasma populations is likely the result of many ecological and genetic factors. For example, like most archaea, both Ferroplasma types I and II lack the well-known mutS and mutL DNA mismatch repair systems and the alternative systems they employ have unknown influence on the relative rates of mutation and recombination.

Implications for species concepts:

The foundation for the biological species concept in eukaryotes is sexual reproduction. In bacteria and archaea, genetic exchange is not necessarily linked to the process of reproduction and can occur between distantly related organism types (Cohan 2002). Accordingly, the applicability of a biological species concept to bacteria and archaea is in question (Lawrence 2002; Nesbø et al. 2006). The reduced rate of genetic exchange seen between recently diverged Ferroplasma types I and II relative to the high rates within each population provides support for the concept that the breakdown of homologous recombination in these archaea serves as a species boundary. We infer a significant difference in the importance of recombination vs. mutation in sequence change within “species” populations (in the range of 2:1 to 4:1) compared to between them (∼1:25). This difference indicates that evolutionary trajectories are much more tightly linked within species groups than between them. Thus, individual populations are diverging.

Recombination between distinct organisms generates mosaic genomes with inconsistent phylogenetic affiliations at different loci around the genome, and it has been suggested that this can preclude meaningful species designations (Nesbø et al. 2006). If large fractions of genomes have inconsistent affiliations, we concur that species designations are not appropriate. This can occur if recombination is very frequent (as is observed within our Ferroplasma populations but not between them) or if recombination is rare but the block sizes are large. Consequently, it may be important to evaluate not just whether interorganism recombination occurs, but also how pervasive it is. We show that this can be achieved using comprehensive genomewide analyses made possible by deeply sampled community genomic data sets. In the case of these Ferroplasma populations (in which ∼2.6% of the genome recently participated in interpopulation recombination), a biological species definition based on the frequency of genetic exchange may be reasonable. This may not be an isolated case because, in many systems, recombination frequencies decrease in a systematic manner as sequence divergence increases. Increasing sequence divergence and sexual isolation will eventually ensure metabolic and ecological differences and distinct species identities.