Genome Evolution in Bacteria Isolated from Million-Year-Old Subseafloor Sediment

We sequenced the genomes of 21 Thalassospira isolates (10 from 3-mbsf and 11 from 6-mbsf sediments) using a hybrid assembly approach consisting of long-read Nanopore sequencing corrected and polished via short-read Illumina technologies at >100× coverage. The mean genome completeness of these hybrid assemblies is estimated to be 99.7% ± 0.3% (mean ± SD), with most being 100% complete and representing the complete chromosome (Table S1). The mean length of the new Thalassospira genomes is 4.71 ± 0.08 Mbp (mean ± SD), with 4,567 ± 107 (mean ± SD) protein-encoding genes, and they are assembled to an average of 12 ± 2 (mean ± SD) contigs (Fig. S1 and Table S1). The genome size and number of protein-encoding genes are similar to those observed within the existing Thalassospira type strains isolated from the surface world (Fig. S1).

Phylogenetic analysis of the full-length 16S rRNA gene sequences from the cultured isolates recovered two clades of subseafloor Thalassospira, affiliated with T. xiamenensis and T. lohafexi, respectively (Fig. 2C). However, phylogenomic reconstruction of the core genomes of Thalassospira type strains and the newly isolated subseafloor strains (1,809 orthologous genes) resolves into three distinct clades of subseafloor Thalassospira isolates (Fig. 3). One clade shares 96 to 97% genome-wide average nucleotide identity (ANI) and 99% 16S rRNA gene sequence identity with T. xiamenensis and T. permensis (Fig. 3). We named isolates in this clade T. xiamenensis strain “Neogene,” after the Neogene eon (2.8 to 23 mya), which covers both estimated ages of the sediment (3 mya and 6 mya) from which the strains in this clade were isolated. The subseafloor genomes in this clade correspond to the 16S rRNA gene ASV1 detected in the in situ frozen sediment core samples (Fig. 2B and C). A second clade contains three isolates from 6-mbsf sediment, which shares 97% genome-wide ANI with T. xiamenensis and T. permensis. Since all isolates in this clade were recovered from ca. 6-mya sediment deposited during the Miocene eon (5.33 to 23 mya), we report them as T. xiamenensis strain “Miocene” (Fig. 3). T. xiamenensis strain Miocene shared 99% 16S rRNA gene sequence identity with these most closely related type strains. The subseafloor genomes in the Miocene clade also correspond to ASV1 detected in the in situ frozen sediment core samples (Fig. 2B and C).

A third clade contains subseafloor Thalassospira cultures isolated only from 3-mbsf sediment and shares 95 to 96% genome-wide ANI with T. lucentensis and T. lohafexi (Fig. 3) and 99% 16S rRNA gene sequence identity with T. lohafexi. The subseafloor genomes in this clade correspond to the 16S rRNA gene ASV2 detected in the in situ frozen sediment core samples (Fig. 2B and C). Based on the genetic distinctness of these isolates, we consider them to be a new candidate species, according to recently provided criteria based on genome-wide ANI (23). Because all isolates of this third clade were recovered from 3-mya sediment, we propose the candidate name “Candidatus Thalassospira pliocenensis,” named after the Pliocene age (2.58 to 5.33 mya) of the deep-sea clay from which they were isolated.

The three subseafloor clades shared 38 to 41% of their genome content (Table S1). However, within each subseafloor clade, the proportion of gene content shared was much higher (>90%) (Fig. 4). This pangenome analysis shows substantial differences in the variable gene content across the three subseafloor clades, providing further evidence that each clade is a genetically distinct population (Fig. 4). Of the genes that are unique to the subseafloor genomes (e.g., core only to the subseafloor genomes), these represented 3 to 6% of the total gene content in the subseafloor genomes (Fig. 4), suggesting that they were acquired before burial. Furthermore, the fraction of the pangenome made up of singleton genes (seen in exactly one genome and no more) is 14% (±4%) of the flexible genome in the type strains (Fig. 4), which is significantly higher than the <0.1% in the subseafloor strains (P = 10⁻⁷ by a two-sided t test). The significantly higher numbers of singleton genes in the type strains are presumably recent acquisitions in the type strain genomes and consistent with a lack of recent gene gain in subseafloor genomes. The higher number of singleton genes in the flexible genomes of the type strains (Fig. 4) is consistent with the longer branch lengths of the type strains in the core genome phylogeny than those of the subseafloor strains (Fig. 3).

The ratio of nucleotide substitutions originating from mutations to those originating from homologous recombination (r/m ratio) can be used to measure the relative effect of homologous recombination on the genetic diversification of populations (24). Due to the physical isolation of individual bacterial cells, reduced cell concentrations, and the reduced availability of extracellular DNA for recombination in subseafloor sediments (25), we hypothesized that r/m ratios would be lower in the subseafloor Thalassospira populations than in the type strains. To test this, we used an established method, ClonalFrameML (26), to calculate the relative ratio of recombination to mutation (R/θ), the mean length of recombined DNA (δ), and the mean divergence of imported DNA (ν) for branch tips (existent genomes) and internal nodes (ancestral states) in the Thalassospira core genome phylogeny, which allows the calculation of r/m [r/m = (R/θ) × δ × ν]. This analysis showed that in the Thalassospira core genome, the r/m value is approximately 10 times lower in the subseafloor core genome (r/m ratio = 0.078) than in the type strains (r/m ratio = 0.71) (Table 1), indicating that homologous recombination plays a much lesser role in the diversification of the subseafloor strains. Because the r/m value estimated for the subseafloor isolates includes three different subseafloor clades, with an ANI differing up to 10% (Fig. 3), the r/m estimate for the subseafloor taxa should be reliable given that there is enough genetic variation to detect recombination. For example, the mean ν value between all of our subseafloor strains (ν = 0.026) (Table 1) is comparable to the values considered sufficient for detecting recombination (26).

Laboratory growth rates of the subseafloor Thalassospira isolates ranged from 0.064 to 0.31 h⁻¹, which were significantly lower than that of the type strain of T. xiamenensis (Fig. 5). These lower growth rates than that of the “surface-world” type strain are consistent with the r/m values of the subseafloor Thalassospira isolates (r/m = 0.078) being anomalously low compared to those of free-living bacteria isolated from the surface world, which have r/m values that range from 0.1 to 64 (24). Concomitant with the 10-fold-lower r/m values than those of the type strains (Table 1), the subseafloor Thalassospira core genomes exhibit far fewer inferred recombination events with imported DNA than the Thalassospira type strains and the ancestral states of the last common ancestors shared with the type strains (Fig. 6).

The relatively low levels of recombination in the genomes of the subseafloor isolates, compared to the type strains, is consistent with a shift away from a recombination-influenced population toward a more clonal population structure since the subseafloor Thalassospira strains diverged from the last common ancestor shared with the type strains. This finding is supported by the pangenome analysis, which shows a significantly (P = 10⁻⁷ by a two-sided t test) lower number of singleton genes in the subseafloor strains than in the type strains (Fig. 4), indicating a lack of recent gene gain in subseafloor genomes compared to the type strains. Clonal growth over millions of years appears to have resulted in a phenotypic change given the relatively lower growth rates of the subseafloor strains (Fig. 5).

We looked for evidence of evolution in the subseafloor strains by investigating the number of pairwise nucleotide-level differences in the subseafloor Thalassospira genomes. We identified tens to thousands of differences within each subseafloor Thalassospira clade (Fig. S2). The differences in subseafloor Thalassospira strains were present in a clade-specific manner (Fig. S3) and consisted of SNPs and insertion-deletion (indel) events. Most (>90%) within-clade differences were represented by indels across all three clades, and all intraclade differences in T. xiamenensis strain Miocene were represented by indels (Fig. S3). While the precise chromosomal locations of SNPs and indels are clade specific, the SNPs and indels often occurred in genes coding for related functions across the different clades. These included indels within genes with predicted annotations involved in flagellar motility (flhB, flhL, fliQ, and fliM), transcription (tetR and fis family transcriptional regulators), cell wall biogenesis (peptidase S41 and peptidoglycan dd-metalloendopeptidase M23), and the transport and metabolism of amino acids and carbohydrates (Fig. S3). This suggests that although the three Thalassospira lineages are evolving independently, the shared environmental setting renders similar traits superfluous across populations and susceptible to genetic changes or loss.

We considered the possibility that the observed nucleotide substitutions and indels occurred within the subseafloor populations during the culture enrichment and cultivation processes. The generation times of Thalassospira (OTU_6) in the incubation measured by quantitative PCR (qPCR) (36 ± 1.5 days [mean ± SD]) indicate an estimated maximum of 15.7 doublings over the enrichment (Fig. S1). Approximately 63 additional generations may be attributable to subsequent cultivation prior to genome sequencing (Fig. S1). At common bacterial mutation or chromosomal deletion rates of 10⁻⁹ to 10⁻¹⁰ mutations bp⁻¹ generation⁻¹ (27, 28), more than 4,000 generations would be required to obtain the observed nucleotide diversity present in the subsurface Thalassospira genomes (Fig. S4). Thus, newly accumulated nucleotide changes in culture are likely insufficient to explain the observed interpopulation nucleotide diversity between the subseafloor genomes (Fig. S4). While it is possible that the observed clonal population structure resulted from the 18-month enrichment, the genomes that we sequenced were not identical, even within clades, which might be expected if clonality was the result of enrichment. For example, the number of pairwise nucleotide differences (including indels and SNPs) in the core genomes of each clade ranged from 0 bp (3_G and 3_H) to 3,602 bp (Fig. S2). At a mutation or deletion rate of 10⁻⁹ per generation, between 10³ and 10⁴ generations would need to accrue to attain this observed genetic diversity (Fig. S4).

Compared to the type strains, the subseafloor Thalassospira genomes exhibit high numbers of pseudogenes (nonfunctional parts of the genome that resemble functional genes) and nonsynonymous substitutions (substitutions that alter the amino acid sequence of a protein). We identified 47.9 ± 8.57 pseudogenes (mean ± SD) in the genomes of subseafloor Thalassospira isolates, which is significantly higher than the number of pseudogenes identified in the type strains (22.1 ± 5.52 pseudogenes [mean ± SD]) (Table 1, Fig. 7, and Fig. S1) (P = 1.5E−10 by a two-sided t test). We investigated whether the genomes from the subseafloor Thalassospira isolates exhibited elevated ratios of nonsynonymous to synonymous substitutions (dN/dS ratios) that might further indicate evolution in clonal populations and relaxed purifying selection across the genome. Similarly, we observed a modest but significant elevation of genome-wide dN/dS ratios in the genomes of the subseafloor Thalassospira strains (0.035 ± 0.006 [mean ± SD]) relative to the type strains (0.022 ± 0.012 [mean ± SD]) (P = 0.0002 by a two-sided t test) (Table 1, Fig. 7, Fig. S1, and Table S2). While the increased dN/dS ratio is consistent with relaxed purifying selection, such small changes are possibly not biologically or evolutionarily meaningful. Additionally, the calculation of dN/dS ratios may be confounded by the high similarity of the subseafloor genomes to each other and to the type strains (see Materials and Methods).

Similar to the nucleotide differences (SNPs and indels) (Fig. S3), the composition of pseudogenes occurred in a clade-specific manner (P = 0.001 by analysis of similarity [ANOSIM]) (Fig. 8). Compared to the type strains, the predicted annotations of the pseudogenes in the subseafloor Thalassospira genomes are skewed toward those involved in transcription, energy conservation, amino acid and carbohydrate metabolism, and flagellar motility (Fig. 8). Moreover, subseafloor genomes have significantly higher numbers of pseudogenes involved in motility (flagellar biosynthesis genes fliN, fliK, and flhO) than the type strains (P = 0.003 by a two-sided t test) (Fig. 9). These flagellar biosynthesis genes are different from the flagellar biosynthesis genes containing intrapopulation indels (fliQ, fliK, flhB, flgH, and flgL) (Fig. S3), but the pseudogenized fliN and fliK genes are syntenic with the indel-containing flagellar genes (Fig. S3), indicating that both indels and pseudogenization have targeted different genes within the flagellar operon.

Burial in the sediment millions of years ago resulted in the physical isolation of subsurface Thalassospira cells and was the impetus for a transition from freely recombining populations with a large N_e to recombination-limited clonal populations with a small N_e that rarely encounter genetically diverse recombination partners that introduce genetic diversity into the population. The combined effects of physical isolation, low cell concentrations, and little or no motility due to low permeability and energy limitation further reduced recombination events. These extreme physical conditions resulted in clonal subseafloor Thalassospira populations that are characterized by the inability to overwrite deleterious mutations through recombination. Our observations of slightly elevated dN/dS ratios and an increased number of pseudogenes (Fig. 7) are evidence the subsurface microbes are experiencing variability in selection, changes that are more likely to be fixed in the absence of homologous recombination. However, because of the small magnitude of the changes, it is challenging to determine whether these observations result from the fixation of deleterious mutations via genetic drift, as is commonly observed in endosymbiotic bacteria (29), or variation in selective pressures resulting from the environmental changes that accompanied burial. For example, some functions appear to be more prone to gene decay than others, such as those involved in flagellar motility. We identified both indels (fliQ, fliK, flhB, flgH, and flgL) (Fig. S3) and pseudogenes (fliN, fliK, and flhQ) (Fig. 8 and 9) in flagellar operons, suggesting that flagellar motility was a slightly beneficial trait in Thalassospira populations before burial but has become a nonessential trait subject to relaxed selection in abyssal clays characterized by very low permeability and an extremely small pore diameter despite their high porosity (18). This observation is further supported by the apparent convergent gene decay in the flagellar biosynthesis operon across the Miocene and Neogene strains (Fig. S5).

Our findings demonstrate that the subseafloor Thalassospira genomes analyzed here have evolved within clonal populations that have experienced a reduction in homologous recombination. Clonal populations can become subject to genetic drift whereby deleterious mutations become fixed and Muller’s ratchet (30) eventually leads to the extinction of endosymbiotic bacterial lineages (29). Our genomes show similar signs of evolution, but the dN/dS ratios observed in the subseafloor Thalassospira genomes are much lower than those observed in endosymbiotic bacteria (31), and subsequent genome reduction was absent. We speculate that the subseafloor Thalassospira isolates have not experienced sufficient generations to observe large changes in dN/dS ratios or genome reduction given that perpetual and severe energy limitation results in increased generation times with depth below the seafloor (2–7).

Subseafloor cell concentrations are low and decrease substantially with increasing depth (1). Because genetic drift has a stronger effect on populations with a small N_e (29), physically isolated microbes experiencing clonal growth and reduced homologous recombination in the deep biosphere may be particularly prone to genetic drift-mediated evolution. Because the subseafloor biosphere contains a large fraction of all bacterial cells on Earth (1), our findings suggest that drift-like evolutionary processes in the absence of homologous recombination may be much more widely distributed in nature than previously thought. Future assessments of homologous recombination and drift in single-cell genomes from uncultured lineages of bacteria and archaea that comprise most subsurface energy-limited communities (32) could be used to assess how widespread this evolutionary mechanism is within the subsurface biosphere.

All samples were taken during expedition KN223 on the R/V Knorr in the North Atlantic, from 26 October to 3 December 2014, at site 11 (22°47.0′N, 56°31.0′W; water depth of ∼5,600 m), via a long-core piston-coring device (∼28 m). Additional details of sampling have been described previously (33). Dissolved oxygen concentrations in the core sections were measured with optical O₂ sensors from the equilibrated core sections and measured with needle-shaped optical O₂ sensors (optodes) (PreSens, Regensburg, Germany) as described previously (33). The dissolved-O₂ data from expedition KN223 are archived and available in the Integrated Earth Data Applications (IEDA) database (https://www.nodc.noaa.gov/archive/arc0127/0178259/2.2/data/0-data/KN223_603663_r2rnav/).

Deep-sea abyssal clay is characterized by very low permeability and an extremely small pore diameter despite its high porosity (18). Deep-sea clay particles have a grain size of <0.2 μm, and the pore space between clay particles is smaller than a bacterial cell, limiting the movement of bacteria through pore space in the clay. Bioturbation can vertically redistribute cells within marine sediment, but bioturbation is restricted to the upper 0.5 m of sediment (34) and thus cannot vertically transport sediment surface material and microbes to depths of 3 and 6 mbsf. Considering the mean sedimentation rate of 1 m per million years, it can be concluded that the bacterial cultures obtained from sediment collected at 3 and 6 mbsf have been physically isolated from the surface world for millions of years.

DNA extractions, qPCR, and 16S rRNA gene sequencing were described previously by Vuillemin et al. (17). In brief, subcores were sampled aseptically with sterile syringes in a UV-sterilized DNA/RNA clean HEPA-filtered laminar flow hood. DNA was extracted from 10 g of sediment transferred into 50 ml of lysing matrix E tubes (MP Biomedicals) containing silica glass beads and homogenized for 40 s at 6 m/s using a FastPrep-24 5G homogenizer (MP Biomedicals) in the presence of 15 ml of preheated (65°C) sterile-filtered extraction buffer (76 vol% 1 M NaPO₄ [pH 8], 15 vol% 200-proof ethanol, 8 vol% MoBio lysis buffer solution C1, and 1 vol% SDS). The samples were incubated at 99°C for 2 min, frozen overnight at −20°C, thawed, and frozen again at −20°C overnight, followed by an additional incubation at 99°C for 2 min and a second homogenization step using the settings described above. After the second homogenization, the samples were centrifuged for 15 min, and the supernatants were concentrated to a volume of 100 μl using 50-kDa Amicon centrifugal filters (Millipore). Coextracted PCR-inhibiting humic acids and other compounds were removed from the concentrated extract using the PowerClean Pro DNA cleanup kit (MoBio). Extraction blanks were included alongside the samples to assess laboratory contamination during the extraction process.

DNA was quantified fluorometrically using a Qubit system with a double-stranded DNA (dsDNA) high-sensitivity kit (Life Technologies). DNA templates were used for qPCR amplifications with updated 16S rRNA gene primer pair 515F (5′-GTG YCA GCM GCC GCG GTA A-3′) and 806R (5′-GGA CTA CNV GGG TWT CTA AT-3′) as previously described (35). All qPCRs were set up in 20-μl volumes with 4 μl of the DNA template, 20 μl of SsoAdvanced SYBR green supermix (Bio-Rad), 4.8 μl of nuclease-free H₂O (Roche), 0.4 μl of primers (10 μM; biomers.net), and 0.4 μl of MgCl₂ and carried out on a CFX-Connect qPCR machine (Bio-Rad) for gene quantification. For 16S rRNA genes, we ran 40 PCR cycles of two steps corresponding to denaturation at 95°C for 15 s and annealing and extension at 55°C for 30 s. All qPCRs were set up in 20-μl volumes with 4 μl of the DNA template. Gel-purified amplicons were quantified in triplicate using QuantiT dsDNA reagent (Life Technologies) and used as a standard. An EpMotion 5070 automated liquid handler (Eppendorf) was used to set up all qPCRs and prepare the standard curve dilution series spanning from 10⁷ to 10¹ gene copies. Reaction efficiency values in all qPCR assays were between 90% and 110% with R² values of >0.95% for the standards. Barcoded PCR for Illumina sequencing was performed using the custom primer dual-indexed approach that targets the V4 hypervariable region of the 16S rRNA gene using dual-indexed 16S rRNA gene primer pair 515F (5′-GTGYCAGCMGCCGCGGTAA-3′)/806R (GGACTACNVGGGTWTCTAAT) (35). Barcoded V4 hypervariable regions of amplified 16S rRNA genes were sequenced on an Illumina MiniSeq platform according to a previously established protocol (35). Bioinformatic processing of these previously published sequence data was described previously by Vuillemin et al. (17) in detail.

Prior to setting up the incubations, the subcores were sampled with sterile syringes using the sample aseptic technique used for DNA extraction. For each sample depth, 7 g of abyssal clay was placed into sterile 20-ml glass flasks and incubated with 4 ml of sterile artificial seawater composed of either H₂¹⁸O (97% atomic enrichment) or unlabeled artificial seawater. Vials were crimp sealed, with an oxygenated headspace of approximately 10 ml, and incubated at 8°C. The artificial seawater was different from the porewater at depth because there was no added nitrate, but there was also no added ammonia, which should be similar to in situ conditions where ammonia is generally below the limit of detection (17). Oxygen was measured continuously throughout the incubations using noninvasive fiber-optic measurements as described previously (36). Small fluctuations in the oxygen measurements in the killed-control and experimental incubations were likely due to temperature fluctuations of the incubator itself (±1°C) since the noninvasive fiber-optic oxygen sensor spots are temperature sensitive (36). Oxygen consumption was detectable over 18 months in slurries consisting of sediment and sterile artificial seawater (see Fig. S1 in the supplemental material), suggesting the presence of actively respiring microbes. The rate of oxygen consumption was calculated by subtracting the starting concentration of O₂ from the final concentration of O₂ and dividing this by the number of days for the incubation.

We used quantitative stable isotope probing (qSIP) to measure the atoms percent ¹⁸O enrichment of actively growing microbial taxa as described previously (17, 37). In brief, after 7- and 18-month incubations, DNA was extracted and subjected to cesium chloride (CsCl) density gradient centrifugation. The same 16S rRNA gene 515F/806R primers (described above) were used for qPCR (described above) to determine density shifts in the peak DNA of buoyant density (BD) for each incubation. 16S rRNA gene amplicons from each fraction resulting from density gradient fractionation were Illumina sequenced as described previously (17). To identify contaminants that may have entered during the fractionation process, we also included in the sequencing run extraction blanks from the SIP fractionation. OTUs containing sequences from extraction blanks were removed. Excess atoms percent ¹⁸O enrichment were calculated for each OTU (including OTU_6, corresponding to the subseafloor Thalassospira strains) according to the equations for quantifying per-OTU atomic enrichment.

The number of doublings for the Thalassospira OTU (OTU_6) detected at the 18-month time point was calculated using the qPCR-normalized relative abundance of the 16S rRNA genes at 0 and 18 months. The number of generations was calculated as X = 2ⁿX₀, Where X and X₀ are the final and initial cell densities, respectively, and n is the number of generations.

After 18 months of incubation in sterile ¹⁸O-labeled artificial seawater, 25 μl of the slurry was plated onto solid medium (10 mg/ml yeast extract and 8 mg/ml agar in artificial seawater [30 mM MgCl₂·6H₂O, 16 mM MgSO₄·7H₂O, 2 mM NaCO₃, 10 mM KCl, 9 mM CaCl₂, 450 mM NaCl]), and after 2 days of incubation in the dark at room temperature, abundant colonies were observed growing on the surface of the petri dishes (Fig. S1). No colonies were observed to grow on control petri dishes that received 25 μl of the ¹⁸O-labeled artificial seawater slurry incubated for 18 months using starting material from the autoclaved sediment (killed controls). This indicated that the colony-forming bacteria were from the sediment and not contaminants introduced during the experimental setup of the incubations. We attempted to culture chemoheterotrophic microbes directly from the collected sediment samples under the same conditions, but no CFU were observed on the petri dishes, even after several months of incubation. Thus, long-term incubation of the sediment at 8°C simply in the presence of added water apparently stimulated the activity of many subseafloor bacteria to a point at which they were able to grow on the surface of a petri dish.

A total of 21 colonies were picked from petri dishes containing colonies from the 3-mbsf and 6-mbsf slurries (10 colonies picked from the 3-mbsf slurry and 11 colonies picked from the 6-mbsf slurry). These colonies were streaked individually onto separate new petri dishes, and a single colony was picked from each of the separate petri dishes (representing the original CFU from the enrichment) and grown in sterile liquid medium (10 mg/ml yeast extract and 8 mg/ml agar in artificial seawater [30 mM MgCl₂·6H₂O, 16 mM MgSO₄·7H₂O, 2 mM NaCO₃, 10 mM KCl, 9 mM CaCl₂, 450 mM NaCl]). Single colonies were then grown in liquid medium. A portion of each of these colonies was used for DNA extraction and genome sequencing, and the remaining volume was frozen as glycerol stocks. Growth rates were determined in experiments with 20-ml crimp-sealed glass flasks containing 0.1 ml of the glycerol stock inoculated into 10 ml liquid medium (10 mg/ml yeast extract in artificial seawater [30 mM MgCl₂·6H₂O, 16 mM MgSO₄·7H₂O, 2 mM NaCO₃, 10 mM KCl, 9 mM CaCl₂, 450 mM NaCl]) with a 10-ml headspace and gentle shaking. The optical density (OD) at 600 nm was measured once every 30 min with a spectrophotometer, and growth rates were calculated from the exponential phase using the formula μ = [ln(X) − ln(X₀)]/(t − t₀), where X and X₀ are the final and initial optical densities (at 600 nm), respectively, and t and t₀ are the times at which those optical densities were made (final and initial, respectively). Under these conditions, the growth rates of the subseafloor Thalassospira strains were similar across all three clades and ranged from 0.064 to 0.31 h⁻¹ (Fig. S5). The growth experiments for all subseafloor and type strains were performed alongside one another in the same bacteriological growth medium formulation, at the same temperature, in identically constructed and prepared flasks, and at the same shaking speed.

Because bacteria can evolve on laboratory experimental timescales (10, 11), we considered the possibility that all diversification and evolution happened during the 1.5-year enrichment. Using the qPCR-based estimate of the doubling time of the subseafloor Thalassospira OTU (OTU_6) in the incubation, which was 36 (±1.5) days, the number of doublings with this rate over this time period would be 15.67. After the enrichment, we spread cells onto petri plates. Assuming that each colony attained 3.3E9 cells per colony (38) from a single cell, we calculate another 31.62 generations per colony (using X = 2ⁿX₀), for a total of 47.29 (15.67 + 31.62 = 47.29) generations per isolate; repeating this isolation with a single cell and colony formation a second time prior to genome sequencing, we arrive at 47.29 + 31.62 = 78.91 generations. According to “Drake’s rule” (28), bacteria experience on average 1 mutation per 300 genomes replicated; thus, the amount of nucleotide diversity (tens to thousands of nucleotide differences) (Fig. S2) that could be accumulated during the incubation is insufficient to explain the observed divergence among the three subseafloor populations. We thus conclude that the interpopulation nucleotide diversity resulted from mutations that were acquired after they were buried. We acknowledge that some cells likely did not survive during the 18-month incubation, which may cause an underestimation of the number of generations. Moreover, upon growth in a colony, additional numbers of generations are passed. Therefore, while our calculated number of generations is an underestimate, we are within an order of magnitude of the real value.

DNA was extracted from the isolates grown in liquid culture until the end of exponential phase as described above. After reaching stationary phase, cultures were pelleted via centrifugation, and the supernatant was decanted. The cell pellets were resuspended in a preheated (65°C) sterile-filtered extraction buffer (76 volume % 1 M NaPO₄ [pH 8], 15 volume % 200-proof ethanol, 8 volume % MoBio lysis buffer solution C1, and 1 volume % SDS), added to lysing matrix E tubes (MP Biomedicals) containing silica glass beads, and homogenized for 40 s at 6 m/s using a FastPrep-24 5G homogenizer (MP Biomedicals). The samples were centrifuged for 15 min, and the dissolved high-molecular-weight DNA in the supernatant was concentrated to a volume of 100 μl using 50-kDa Amicon centrifugal filters (Millipore). The concentrated extract was cleaned of proteins and other nongenomic DNA organic matter using the PowerClean Pro DNA cleanup kit (MoBio). Extraction blanks were added alongside the samples to assess laboratory contamination during the extraction process. Genomic libraries were prepared using the Nextera XT DNA library prep kit (Illumina). Quality control and quantification of the libraries were performed on an Agilent 2100 Bioanalyzer system using high-sensitivity DNA reagents and DNA chips (Agilent Genomics). Genomic libraries were sequenced to a depth of ca. 100× coverage using a high-output paired-end 2-by-150 sequencing reagent kit (Illumina).

In addition to Illumina sequencing, the high-molecular-weight genomic DNA was sequenced using the Nanopore MinION system. Sequencing libraries for the MinION system were prepared using the ligation sequencing kit (Oxford Nanopore Technologies) according to the manufacturer’s instructions. Barcoded libraries were sequenced on the MinION system using a Flongle R9 flow cell and base called and demultiplexed using the MinIT system with ont-minit-release v19.12.5 and ont-guppy-for-minit v3.2.10 for base calling (Oxford Nanopore Technologies).

A hybrid assembly was performed using both the short (Illumina)- and long (Nanopore)-read sequencing data using Unicycler (v.0.4.0), which uses de novo-assembled Illumina data from SPADES to polish the de novo-assembled contigs obtained from Nanopore data using RACON (39). We used the default settings available in Unicycler to coassemble the Illumina data (forward and reverse fastq files) and Nanopore data (single fastq file), which consisted of the following commands: ./unicycler -1 isolate-6-1_illumina_R1_001.fastq -2 isolate-6-1_illumina_R2_001.fastq -l isolate-6-1_Nanopore.fastq –spades_path ./SPAdes-3.13.0-Darwin/bin/spades.py –racon_path ./racon/build/bin/racon –bowtie2_build_path ./bowtie2-build –bowtie2_path ./bowtie2 -o out. The combined assemblies of Illumina and Nanopore data resulted in a relatively low number of contigs (9 to 12 per genome) and a predicted genome completeness of 100% for nearly all genomes (Table S1). Genome completeness was determined using CheckM (40), using the family Rhodospirillaceae as the reference group (which contains all species within the genus Thalassospira). Genomes were annotated using RASTk (41).

The core genome was defined as the set of orthologous genes that were shared in all subseafloor and “type strain” Thalassospira genomes. NCBI genome accession numbers for the type strain genomes of Thalassospira (those isolated from surface-world habitats) are as follows: NZ_JAATJD010000001.1 for T. tepidiphila, NZ_JPWA00000000.1 for T. xianhensis, NZ_ATWN01000001.1 for T. lucentensis, NZ_FTON01000032.1 for T. xiamenensis, NZ_AMRN01000001.1 for T. profundimaris, NZ_FNTU01000002.1 for T. permensis, NZ_JFKB00000000.1 for T. alkalitolerans, NZ_CP031555.1 for T. indica, NZ_CP024199.1 for T. marina, NZ_PGTS00000000.1 for T. povalilytica, NZ_NXGX00000000.1 for T. lohafexi, NZ_JFKA00000000.1 for T. mesophila, and NZ_JRJE00000000.1 for T. australica. Orthologous genes between the subseafloor and type strains were defined as those sharing >30% amino acid similarity to the collective suite of genes within the type strain Thalassospira xiamenensis M-5. T. xiamenensis M-5 was chosen as the reference genome for this purpose because it is the only publicly available genome of a cultivated Thalassospira isolate that is completely closed and represents a single chromosome and a 190-kb plasmid. A total of 1,809 orthologous genes (defined as being present as only a single copy in each genome) were identified that are carried by all Thalassospira strains that had >30% sequence similarity to genes within the T. xiamenensis M-5 genome. We define these genes as orthologs, as opposed to paralogs, because they are present as only a single copy in each genome. Paralogs were identified in Thalassospira genomes by performing similarity searches for genes (BLASTp queries) that had similarity to a T. xiamenensis gene (BLASTp subjects) and then removing subjects (potential core genes) that had hits from multiple queries from the same query genome. Only T. xiamenensis genes (BLASTp subjects) that had similarity from a single BLASTp query across all Thalassospira strains were kept as orthologs. This should select only single-copy genes (genes that have not been duplicated within the genomes). Because the genes have homology between the different genomes and are present as a single copy across all genomes compared, we interpret them to be orthologs. Each of these 1,809 orthologous genes was individually aligned between all Thalassospira strains using MUSCLE (42), and the 1,809 individual alignments were then concatenated into a single core genome alignment for the subsequent phylogenomic analysis using Geneious Prime (version 2019.2.1). After the concatenation of all core genes, the total size of the core genome alignment was 1,817,073 nucleotide characters, with 34 taxa (21 subseafloor strain and 13 type strain taxa). A maximum likelihood phylogeny was created using PhyML (43) with a generalized time-reversible (GTR) model of evolution and 100 bootstrap replicates, which was implemented in SeaView (44). The resulting phylogenetic tree and the concatenated core genome alignment were used as inputs for subsequent ClonalFrameML and dN/dS ratio (ratio of nonsynonymous to synonymous mutations) analyses. Full-length 16S rRNA gene sequences were amplified from the isolated subseafloor bacteria from DNA extracted from cell pellets (using the protocol described above) using PCR (30 cycles of 95°C for 30 s, 60°C for 30 s, and 72°C for 2 min) with universal PCR primers 27F (5′-AGAGTTTGATCMTGGCTCAG-3′) (45) and 1391R (5′-GACGGGCGGTGTGTRCA-3′) (46) and Sanger sequenced, as a quality control measure. Quality control and assembly of the Sanger-sequenced 16S genes (forward and reverse sequences, each ca. 850 bp in length) were performed in CLC Genomics Workbench using the default base-calling quality and de novo assembly settings for Sanger sequence data. Phylogenetic analyses were performed on the Sanger-sequenced 16S sequences with maximum likelihood methods as described above.

The contributions of mutations and recombination to the genomic diversity in the concatenated core genome alignment, the number of recombination events (imports) per genome, and the positions of recombination hot spots were investigated using ClonalFrameML (26). This approach allows the assessment of homologous recombination within the core genome only (e.g., genes shared between all compared genomes) and does not assess nonhomologous gene acquisition (e.g., horizontal gene transfer) or gene loss events. Nucleotides unaffected by homologous recombination are referred to as unimported, and nucleotides subject to recombination are referred to as imported (26). ClonalFrameML provides the relative ratio of recombination to mutation (R/θ), the mean length of recombined DNA (δ), and the mean divergence of imported DNA (ν). These results were used to calculate the relative contribution of recombination versus mutation to the overall genomic diversity (r/m), using the formula r/m = (R/θ) × δ × ν. ClonalFrameML was performed in three separate runs, containing a core genome alignment that contained (i) all genomes, (ii) only the subseafloor genomes, and (iii) only the type strains. The resulting r/m values from these three groups (presented in Table 1) were then used to interpret the relative importance of mutations compared to recombination in the separate groups (e.g., type strains versus subseafloor strains). The r/m analysis provided by ClonalFrameML is capable of identifying recombination between closely related strains where the mean ν value can be as low as 0.031 (26). The mean ν value between our closely related subseafloor strains is slightly lower but comparable (ν = 0.026) (Table 1). The mean divergence of imported DNA that we describe is like those described in the original ClonalFrameML methods paper, suggesting that the r/m values that we calculated can detect recombination between the closely related subseafloor strains.

In addition to calculating sites and rates of recombination in the core genome, ClonalFrameML also estimates the ancestral sequences at internal nodes of the clonal genealogy and any missing base calls in the observed sequences. The reconstruction of ancestral sequence states is performed using maximum likelihood, and the ClonalFrame model can be thought of as a hidden Markov model (HMM) when the ancestral and descendant genomes for each branch of the clonal genealogy have been observed or reconstructed (26). The hidden state of the HMM records whether each nucleotide was subject to recombination or not on the branch connecting the two genomes. We acknowledge that drawing inference under the resulting ancestral recombination graph is a notoriously complex statistical problem (26). Instead, here, we use ClonalFrameML only to assess within-group recombination (e.g., between species within the genus Thalassospira), and thus, our analysis cannot assess the influence of external recombination (from species outside the genus Thalassospira).

The dN/dS ratio in the core genome alignment (global ω ratio) was estimated using HyPhy v2.2.4 (47) and applying the adaptive branch site random-effects likelihood (aBSREL) approach (48) to all branches in all subfamilies. The dN/dS calculation is based on an in-frame codon alignment of the core genome, together with the corresponding core genome phylogeny. As described previously by Smith et al. (48), aBSREL tests for each branch in the phylogeny (e.g., each genome in the core genome alignment) and whether a proportion of sites have evolved under positive selection and will infer the optimal number of ω classes (dN/dS) for each branch. aBSREL takes into consideration that different branches may feature evolutionary patterns with differing complexities and hence may be better modeled by more or fewer ω classes. To this end, aBSREL uses AICc (small-sample AIC [Akaike information criterion]) to infer the optimal number of ω for each branch in the core genome phylogeny. Tests for statistical significance of branches under positive selection (e.g., genomes with a significant dN/dS ratio) in aBSREL are based on likelihood ratio test (LRT) distributions, whereby a null model is defined in which no positive selection is allowed on a branch, and the LRT is used to determine whether the null model can be rejected (48). The aBSREL algorithm that calculates the global ω (dN/dS) ratio does not provide accurate results for positive selection when highly similar sequences are included in the in-frame codon alignment (48). Therefore, we calculated the dN/dS ratios from the in-frame codon core genome alignment with aBSREL 10 separate times, with each run containing one selected genome from each of the three subseafloor clades: (i) T. xiamenensis strain Neogene, (ii) T. xiamenensis strain Miocene, and (iii) “Ca. Thalassospira pliocenensis.” The results for each of these separate aBSREL runs are summarized in Table S2. By running aBSREL in this manner (multiple times using only one representative of highly similar subseafloor clades), statistically significant dN/dS values were obtained for each subseafloor genome (Table S2). The 10 different aBSREL run output files produced from the HyPhy program (.json files) were then visualized in HyPhy Vision (http://vision.hyphy.org/aBSREL), a Web graphical user interface that allows visualization of branch site selection, P values, and global ω (dN/dS) classes for each branch (genome) from the in-frame codon core genome alignment.

Low dS values may produce unreliable dN/dS ratios. aBSREL does not calculate separate dN or dS values. We investigated how the magnitude of dS in the subseafloor strains compared to that of the type strains using the FitMG94 workflow in HyPhy (https://github.com/veg/hyphy-analyses/tree/master/FitMG94) to fit the Muse-Gaut model of DNA sequence evolution (49) to the core genome alignment, which reports the number of synonymous and nonsynonymous substitutions per nucleotide site for each branch on the tree. The FitMG94 workflow also uses a corrected empirical estimator (CF3x4) that provides improved estimates of several parameters in the evolutionary model (50). The results of FitMG94 showed that there was no significant difference in the magnitude of dS in the type strains (0.02 to 0.36) compared to the subseafloor strains (0.001 to 0.07) (P = 0.59 by a two-sided t test).

All subseafloor and type strain Thalassospira genomes were analyzed in Anvi’o v6.2 using the pangenome workflow (51). Briefly, each genome was converted into an Anvi’o contig database. Genes were functionally annotated using eggNOG v5.0 (52) with eggNOG-mapper (53) and imported back to each genome’s Anvi’o contig database. Genome storages were generated using anvi-gen-genomes-storages, and anvi-pan-genome was deployed with the parameters –min-bit 0.5, –mcl-inflation 10, and the flag –use-ncbi-blast. The Anvi’o pangenome database and a summary of gene clusters are stored in Figshare (https://figshare.com/s/06ba1287a00ab01a1ee).

We estimated the number of pseudogenes within the genomes using two programs, Psi-Phi (54) and DFAST (55). Psi-Phi uses a conservative criterion considering a pseudogene only when it lost >20% of its original length and enhances pseudogene recognition among closely related strains both in annotated regions by identifying incorrectly annotated open reading frames (ORFs) and in intergenic regions by detecting new pseudogenes (54). Psi-Phi classifies pseudogenes as either identified pseudogenes or possible but potentially not pseudogenes. To be conservative, we considered only genes identified as pseudogenes from Psi-Phi and did not consider those flagged as potential pseudogenes. As a second check of the pseudogene content, we searched genomes for pseudogenes using DFAST (55). The estimated number of pseudogenes per genome was then taken as an average of the numbers detected using both methods (Psi-Phi and DFAST). On average, Psi-Phi identified a higher number of pseudogenes per genome (57 ± 10) than DFAST (32 ± 4), but the variation between methods for the same genome was consistent (average variation = 27; standard deviation of averages = 7). This minimal variation between individual genomes indicates that biases inherent to the pseudogene prediction methods affected the different genomes equally and thus allow for a pseudogene comparison between the genomes.

Data are publicly available through NCBI BioProject accession number PRJNA473406. The raw genomic sequence reads are available in the SRA under BioSample accession numbers SAMN17168194, SAMN17168195, and SAMN17168196. The hybrid assemblies (contigs) for each of the subseafloor isolate genomes are available through the LMU Open Data database (https://doi.org/10.5282/ubm/data.233). The 16S data from the frozen core are available in the SRA under BioSample accession numbers SAMN10929403 to SAMN10929517. Figures and output files from the pangenomic analysis in Anvi’o are available at Figshare (https://doi.org/10.6084/m9.figshare.13372619). Additional data related to this paper may be requested from the authors.