Genome dynamics in a natural archaeal population

Isolate Versus Population Composite Genome Analysis.

The fer1 isolate genome sequence was assembled from data obtained from small and large insert libraries. General features of the 1.94-Mb genome, classification of the 1,963 ORFs by functional category, and a metabolic reconstruction are reported in supporting information (SI) Table 3 and SI Figs. 6 and 7.

Composite scaffolds and contigs from the fer1(env) population were assembled independently of the fer1 isolate sequence. Assembled fragments were ordered and aligned by using the fer1 genome as a reference sequence to reconstruct ≈92% (1,792 genes) of the fer1 isolate genome (Fig. 1a), with an average read depth of 4.5×. The colinearity plot (Fig. 1b) illustrates a very high degree of conserved synteny between fer1 and fer1(env).

Fig. 1 highlights the differences between the fer1 and the composite fer1(env) genomes. There are 45 genes in fer1(env) not detected in the fer1 isolate and 152 genes in fer1 not present in fer1(env). These gene products are predominantly hypothetical and conserved hypothetical proteins, as well as transposases, integrases, transport-related proteins, products involved in restriction/modification and DNA repair, and glycosyltransferases (see SI Table 4). No missing genes were found in functional categories involved in core energy metabolism, lipid biosynthesis, cofactor biosynthesis, or transcription. Of the unique fer1 isolate genes, 79% occur in distinct genomic blocks (the largest comprising 24 genes; Fig. 2 and SI Fig. 8). Such regions often posses anomalous G+C contents (Fig. 1a) and are often associated with phage-type integrase genes. Thus, we infer that the blocks of unique genes are the result of prophage insertion. The localization of integrated elements at tRNA genes and the presence of partitioned integrase gene fragments flanking the insertion is consistent with a mechanism of site-specific integration typified by the SSV-type archaeal integrated elements (14).

Movement of insertion sequence (IS) elements also occurs fast enough to give rise to strain-level differences (Fig. 1a). There are 94 transposase insertions in fer1 and 51 in fer1(env) representing 10 families of ISs (SI Table 5). The much higher incidence of ISs in fer1 may indicate that these mobile elements proliferated while the organism was in culture. Only 18 of the 43 transposases in common to fer1 and fer1(env) occur in the same location (the 25 rearranged elements are shown in Fig. 1b). A few notably proliferous transposases have been copied into multiple locations, particularly in fer1 (e.g., seven identical copies of one novel transposase). As sources of strain-level diversity, together with the effects IS elements can elicit on gene expression and genome topography (15–20), heterogeneous IS dynamics may alter the adaptive fitness of population members population over short time scales.

Genomic heterogeneity between fer1 and the composite fer1(env) sequence was quantified by determining the number of differences that occur across comparable genomic regions (Table 1). The average nucleotide-level divergence between fer1 and fer1(env) is ≈3%. However, ≈50% of proteins encoded by the composite fer1(env) sequence are identical to the orthologous isolate sequences. The average amino acid-level similarity among orthologs is 99.16%. Overall, blocks of strain-specific genes localized to integrated elements represent the largest input of novel genetic potential within the population. The gain or loss of single or a few genes (mostly hypothetical) is the second most significant source of strain heterogeneity. In contrast, gene duplications (paralogous sequences) are very infrequent events and gene order is almost completely preserved.

Environmental Genotypic Heterogeneity.

Comparison among composite genome scaffolds revealed that coexisting members of the fer1(env) population exhibit differences in gene content. When a subset of individuals contain genes not present in all population members, the same genomic region will be assembled into two or more discrete fragments (10). The genes most commonly responsible for this phenomenon are transposases. In addition to documenting their presence/absence in composite genome sequences, the insertion of transposases into a subset of population members is evidenced by the localized “stretching” of mate-pairs (paired end reads) across genomic regions.

Integrated elements also contribute to genetic heterogeneity between fer1(env) members. Fig. 2 shows that fer1(env) variant 1 is missing three genes comprising a complete type I restriction-modification (R/M) system. In fer1, this region is adjacent to a tRNA^Leu and includes an insertion of 24 putative prophage-associated genes flanked by a partitioned integrase gene. Interestingly, a distinct block of possibly phage-derived genes is also present in the same region in fer1(env) variant 2, suggesting that the same genomic region has been targeted by disparate phage. Consequently, the form of heterogeneity that distinguishes individuals from one another within the fer1(env) population is similar to that which distinguishes fer1 from fer1(env). Whereas the importance of phage as sources of novel genes that differentiate species and strains has been documented (5, 14, 21, 22), our results indicate that the insertion and loss of blocks of genes of probable phage origin occurs fast enough to differentiate the gene content of members of a coexisting population and may give rise to phenotypic differences. These results extend previous findings documenting extensive genomic variation within natural coexisting bacterial populations (23), including the significance of phenotypic differentiation mediated by phage (24), and serve to highlight the dynamic nature of microbial population structure.

Analyses presented above rely upon comparison between the fer1 isolate and composite fer1(env) sequences. However, greater resolution into the fer1(env) population structure can be achieved because each environmentally derived sequencing read likely originated from a different individual within the natural population (9). We evaluated sequence variability in the population on a gene-by-gene basis by aligning all reads constituting the fer1(env) population with the nucleotide sequence of fer1 (average coverage at each locus is ≈4.5×). Overall, 56% of the genes/gene fragments encoded on fer1(env) reads (out of a total of ≈14,000 reads) have the fer1 sequence type (defined as having ≥99.9% nucleotide-level similarity to the fer1 ortholog). At least one read has the fer1 sequence type at ≈80% of loci (see gene content dendogram, SI Fig. 9).

We categorized all gene types in the fer1(env) population as either homogeneous (one sequence type only, 30.4%) or heterogeneous (multiple variants at one locus, 69.6%) (Figs. 3 and 4). Of the homogeneous group, 78.3% of genes are isolate type (427 genes in total), indicating that a significant number of genes possessing the isolate sequence type have persisted over the >4 years since the isolate was obtained. The genome-wide distribution of sequence variation in the fer1(env) population is summarized in Table 2.

To quantify and evaluate the distribution of gene variants within the fer1(env) population, we evaluated several >100-kb regions in detail, verifying that mate-pairs carrying environmental variant sequences anchor them into the same genomic context (see 109-gene region in Fig. 5). Of the fer1(env) heterogeneous loci, 55% possess two variants and 45% possess three or more variants (these products share on average ≈99% amino acid level identity with fer1 isolate orthologs).

Given the ≈4.5× sequencing depth across the fer1(env) population, we cannot determine whether regions are homogeneous because variants were absent or simply present in low abundance. However, our results do show the dominance of sequence types present. In the future, PCR-based analyses together with population genomic sampling will provide the sequence depth necessary to more fully resolve the abundance distribution of gene types in the population.

Mate-pair information (Fig. 5b; see also SI Fig. 10) and direct examination of nucleotide polymorphism patterns within the reads confirm that the population does not consist of a few dominant clonal genome types that have been sampled to different depths of coverage. Several kilobase long regions characterized by only isolate-type genes transition into regions with several gene variants at each locus (Figs. 4 and 5). Multiple linkage patterns between isolate-type and non-isolate type genes can be identified in single regions and some reads capture the transitions directly (Fig. 5 and SI Fig. 11). Thus, it is evident that the population exhibits a mosaic genome structure inferred to have arisen via recombination between distinct yet closely related strain variant sequence types. Frequently, but not exclusively, the site of recombination occurs at the same nucleotide position in multiple reads indicating that the resulting recombinant sequence has been amplified in the natural community, either by natural selection or by genetic drift.

Extensive recombination has been documented in Ferroplasma type II (9). Subsequently, MLST-based analyses provided evidence for genetic recombination in the halophilic archaeon Halorubrum sp (25) and the thermoacidophilic archaeon Sulfolobus islandicus (7). Results documenting recombination in F. acidarmanus suggest that mosaic genome structure may be a common feature of natural archaeal populations.

The restriction subunit of the type I R/M system described above exemplifies the form of sequence heterogeneity that can exist in the fer1(env) population. There are three distinguishable sequence types for the restriction subunit, exhibiting at least five recombination events across the 2,985-bp gene. The fer1(env) population also contains a type I specificity subunit gene with several allelic variants that show evidence of recombination among them. Recombination among individuals with distinct sequence types may provide the capacity to generate a large number of novel, hybrid restriction enzymes. Overall, R/M system dynamics may factor significantly in natural population structure and evolution (26–28) as heterogeneity in R/M activity (diversity in molecular recognition or the complete presence/absence of such systems) may lead to “recombination isolation.” Ultimately, the absence of recombination as a cohesive mechanism hindering genetic divergence among population members may result in an altered evolutionary trajectory (e.g., speciation or extinction).

Population Structure.

An important question is the role of selection in shaping genetic variation in natural populations. It has been postulated that genetic homogenization can be achieved when an adaptive genotype providing a selective advantage sweeps through a population (29, 30). Cohesion can also be maintained by high levels of genetic exchange between population members (31). As suggested here for the fer1(env) population and recently for local populations of Ferroplasma type II (9) and S. islandicus (7), recombination prevents periodic selection from purging total diversity. Because gene variants are not linked to a single genome type, selection removes only the subset of mosaic genome types with less adaptive gene variants. The removal of sequence variation by incomplete selective sweeps may explain the homogeneous genomic regions within the fer1(env) population.

It has been suggested that genetic diversity can be maintained within a clade of very close relatives by one of two mechanisms: (i) the clade contains only a single population, and selection for adaptive genotypes is too weak to purge diversity (23, 32), or (ii) the clade contains multiple ecologically distinct populations (ecotypes) (33). In this latter model, localized chromosomal regions of sequence identity are explained by rare “globally adaptive” mutations, which are beneficial in the context of any ecotype to which they are introduced (29). When these mutations are transferred by recombination between ecotypes, selection may result in homogenization of the recombined region across ecotypes, while leaving the rest of the genome heterogeneous.

We cannot rule out the possibility that some sequence diversity in the fer1(env) population is due to the partitioning of specific mosaic genome types into different niches within the biofilm. In support of the multiple-ecotype model, it is clear that recombinant genotypes have been amplified throughout the clade and that heterogeneous loci fall into discrete clusters. However, as our data indicates that nearly a third of the genes have been involved in recent global selection-like events, it is more likely that the fer1(env) organisms share a single niche, are subject to very similar selective pressures, and thus represent a single population. Furthermore, high rates of recombination effectively blurs the distinction between ecotypes (29). Ultimately, definitive discrimination between either model will require additional genomic monitoring of these populations over space and time.

Evidence for Selection.

The model for the mosaic population structure requires selection to explain the presence of genes with only one sequence type. To seek evidence for selection we modeled synonymous vs. nonsynonymous nucleotide polymorphism patterns in the fer1(env) population, both between population variants and against the respective fer1 ortholog. Pairwise dN/dS ratios (ratios of nonsynonymous to synonymous nucleotide substitutions) among partial and complete orthologous sequences were determined (34). Approximately seventy-nine percent of the 3,745 pairwise comparisons among variable sequence segments exhibited a dN/dS ratio <0.15, indicating that the majority of the genome is under strong purifying selection. No functional bias in distribution of dN/dS values could be detected (data not shown).

Only six cases were identified with dN/dS values >1.0, suggesting that they are under positive selection. Four of these genes encode small hypothetical proteins, one is a conserved hypothetical protein, and one is a predicted siroheme synthase [siroheme, an iron-containing porphinoid, serves as a cofactor for sulfite and nitrite reductases (35)]. In many cases, elevated dN/dS values were associated with regions of genes, rather than entire genes. More sensitive methods for identifying positive selection (e.g., on a codon-by-codon basis) will require reconstruction of complete environmental variant gene sequences and the inclusion of sequences from more distantly related species. Whereas it is expected that many genes would display low dN/dS ratios (e.g., housekeeping genes), these results suggest that stringent functional constraints [including codon usage bias, the selection for preferred synonymous codons (36)] are characteristics of adaptation to this environment. These results are consistent with analyses revealing strong functional constraints in thermophilic species as revealed by global dN/dS analyses among closely related organisms (37).

To seek evidence for the basis of selection, we examined the physical and functional distribution of homogeneous and heterogeneous genes within the fer1(env) data set. Clusters of genes with only one sequence type may be moved in blocks, possibly by recombination events (Figs. 4 and 5). No significant bias in functional gene category representation was detected among homogeneous and heterogeneous regions (SI Fig. 9). Whereas we might predict that genes with a single sequence type may confer metabolic properties that have been important in recent adaptation, greater sequence depth will be required for definitive analysis.

Population Structure in the Context of Environmental and Geologic History.

To understand population dynamics it would be helpful to able to correlate environmental perturbations with their genomic and evolutionary consequences. Geologic and historic records can provide information about the sequence of geochemical changes at the site. Using magnetism recorded in hematite gossans (products of pyrite oxidation), Alpers et al.^†† determined that ore bodies at Iron Mountain were first exposed by erosion ≈700,000–1,000,000 years ago. When mining began ≈150 years ago, the physical fracturing of the system resulted in increased porosity and permeability within the mountain thus accelerating pyrite oxidation. As solution chemistry and temperature are strongly controlled by the dissolution of pyrite, this process resulted in a decrease in solution pH to values <1.0 and increased temperatures to ≈40°C in the contemporary underground mine. Seasonal rainfall-driven fluctuations in temperature and ionic strength occur within the mine and have been shown to be factors in species selection (38). These factors may also be important in strain genotype selection.

Our comparative analyses have revealed genomic heterogeneity that may be understood in the context of these environmental changes. First, we determined that the fer1(env) population consists of mosaic genomes in which a small number of sequence types predominate. One possible scenario is that subsurface populations that were previously evolving in isolation from one another were brought together and began to recombine when different regions of the mountain were physically interconnected by mining. Recombination-driven homogenization of distinct but closely related populations at their geographic boundaries may be a common process in many ecosystems. Perhaps a limited number of mosaic genome type rose in abundance relative to a diverse mosaic pool in a relatively recent clonal expansion (39), possibly induced by environmental changes caused by the onset of mining. Under this scenario, genes with only a single sequence type may have derived from any of the ancestral strains.

As recombination rates are inferred to be very fast relative to nucleotide polymorphism fixation rates (otherwise we would see many closely related variant types), rapid recombinatorial mixing of variants into the mosaic genome pool may maintain population cohesion. However, mechanisms such as phage and IS dynamics contribute to strain-level variation over short time scales via the rapid introduction and distribution of novel genes. Ultimately, accumulated variation may diminish the frequency of genetic recombination, promoting divergence and possibly speciation (31).

A second finding is that almost all genes are under strong stabilizing selection. Environmental sequence types that are nearly identical (≥99.9% similar) to the isolate type are found at most loci despite >4 years difference in sample collection time. The exact preservation of a large amount of DNA sequence over 4 years may provide one constraint on rates of nucleotide polymorphism fixation. Such information is needed to quantify rates of genome evolution. At present, it is only possible to speculate about the relationship between strain divergence and site history. One possibility is that all of the observed nucleotide-level heterogeneity arose in situ over 700,000–1 million years. If this possibility could be confirmed, the relative rates for the accumulation of the different forms of genome change (Table 1) could be converted to absolute rates. Support for this correlation may be obtained by future genomic monitoring of these populations over space and time and through deductions at other sites for which dates of onset of colonization can be determined. This approach may provide a route to the long-term objective of mapping genome change onto environmental change so as to discover how these processes are interconnected.

F. acidarmanus fer1 Genome Sequencing and Assembly.

F. acidarmanus fer1 was isolated from enrichment cultures inoculated with mine water and sediment samples collected in July 1997 from the five-way site within the Richmond Mine (11). Details of the fer1 genome assembly and annotation are presented in SI Materials and Methods.

fer1(env) Genome Assembly and Comparative Analyses.

We compared the isolate fer1 genome with genomic information obtained by shotgun sequencing small insert libraries constructed from DNA extracted directly from an environmental biofilm. The biofilm was sampled in 2002 from the five-way site of the Richmond Mine at the same location from which the fer1 isolate was collected. Tyson et al. (9) used 76 Mb of environmental sequence data to reconstruct large genome scaffolds from five species, one of which was closely related to F. acidarmanus. The partial composite genome reported for this fer1(env) population was previously referred to simply as Ferroplasma type I (9). In this work, an additional 27 Mb of community small insert library sequence not reported in ref. 9 was coassembled with the initial 76 Mb of data. The total community genome data set used in the present study consisted of 141,312 small-insert library reads (average trimmed length 737 bp) for a total of ≈103 Mb of environmental sequence data. Assembly was performed by using PHRAP (minmatch 50; penalty −15; minscore 40) (40) and ATLAS (k-mer size 32 bp, MaxMismatch 10%, band size 7; minmatch 12, minscore 20, penalty −2, forcelevel 10) (41). Newly generated genomic contigs were binned by alignment with existing contigs and scaffolds described in Tyson et al. (9).

Scaffolds and contigs were initially assigned to fer1(env) based solely on assembly and subsequent analyses including % G+C, read depth, and profiling of nucleotide frequencies using TETRA (42). Subsequently, community genome contigs and scaffolds assigned to fer1(env) were verified by comparison with the F. acidarmanus fer1 genome by using BLAST (43) and MUMmer 3.0 software package's NUCmer (44). Approximately 14,000 environmentally derived sequencing reads constituted the genomically sampled fer1(env) population. In total, 172 scaffolds belonging to fer1(env) were used in the comparative analyses totaling >1.79 Mb of assembled sequence. Genome circular diagrams were generated from .embl format genome files exported from ARTEMIS (45) by using “circular diagram.pl” (Pathogen Group, The Wellcome Trust Sanger Institute). Synteny dot plots were generated by alignment of the fer1 genome sequence with all fer1(env) composite contigs and scaffolds by using NUCmer and MUMMERPLOT (44). Orthologous sequences between fer1 and the composite fer1(env) sequences were identified by using the reciprocal smallest distance method described by Wall et al. (46). To determine genome-wide sequence variation in the fer1(env) population, all fer1(env) environmental sequencing reads were aligned to fer1 genes by using BLASTN, and sequence identity was analyzed on a gene-by-gene basis. Multiple alignment outputted blast results were also generated, and mate-pair consistency for all reads localized to a given fer1 gene were verified on a gene-by-gene basis to verify that environmental reads aligned to fer1 genes in the correct genomic context.

Analysis of Selection in the fer1(env) Population.

To screen for evidence of selection within the fer1(env) population and against the orthologous fer1 sequences, the set of predicted fer1 orfs was queried against all reads from the community genome sequence by using BLASTN with a strict E value of e−75. To verify correct assignment of environmental reads to the fer1(env) population, reads comprising independently assembled fer1(env) contigs were extracted and checked for consistency against all blast matches. All alignments were preformed in CLUSTALW. dN/dS ratios for pairwise comparisons between fer1 orfs and orthologous fer1(env) reads were determined by using the method of Nei and Gojobori (34) in SNAP (47). Correct assignment of environmental sequences to the fer1(env) genome and the accuracy of alignments were confirmed for all pairwise comparisons resulting in dN/dS values >0.9 by reconstructing variant alleles and estimating dN/dS values in DnaSP (48).