Genome Streamlining, Proteorhodopsin, and Organic Nitrogen Metabolism in Freshwater Nitrifiers

To map free-living picoplanktonic (here defined as cells that pass through a 1.6-μm filter) nitrifiers across the Great Lakes, we searched our recent 16S rRNA data sets for known nitrifying taxa (34). We detected putative AOB in the genus Nitrosospira (Betaproteobacteria, family Nitrosomonadaceae) and AOA in the genus Nitrosarchaeum (family Nitrosopumilaceae), along with putative NOB in the genera “Ca. Nitrotoga” (Betaproteobacteria, family Gallionellaceae) and Nitrospira (family Nitrospiraceae). We did not detect 16S rRNA amplicons from Nitrosococcus, Nitrococcus, Nitrospina, or Nitrobacter. The highest relative abundances of picoplanktonic nitrifiers were observed in deep samples from eastern Lake Erie and Lake Ontario (9 to 24% of total amplicons), compared to 2 to 14% in Lakes Michigan, Huron, and Superior. Lakes Erie and Ontario also have higher cell concentrations and higher surface chlorophyll (see Data Set S1 at https://doi.org/10.6084/m9.figshare.15130350.v4). The relative abundance of nitrifiers was negatively correlated with photosynthetically active radiation (PAR; Spearman’s rho = −0.89, P < 2.2e−16) and reached a maximum below the depth of 1% PAR in each lake, up to 20% of amplicon sequences (see Fig. S1a in the supplemental material). The relative abundances of ammonia- and nitrite-oxidizing taxa were strongly correlated (Spearman’s rho = 0.918, P < 2.2e−16) (Fig. S1b). Picoplanktonic nitrifiers were rare (<0.1% relative abundance) in bottom water samples from the southern basin of Lake Huron (HU15M) and the western basin of Lake Erie (ER91M); these two stations are the shallowest in our data set and have relatively high light penetration to the bottom (∼1% PAR). Chlorophyll a concentration was also negatively correlated with the relative abundance of nitrifiers (Spearman’s rho = −0.677, P < 1.7e−7) (Fig. S1c). These findings are consistent with previous work demonstrating photoinhibition of nitrification (37–40), as well as potential competition with phototrophs for ammonium (41).

The taxonomic assemblage of nitrifiers differed across lakes and even among stations within a lake (Fig. 1; see Data Set S1 at https://doi.org/10.6084/m9.figshare.15130350.v4), in association with variable productivity and nitrogen availability. Surface ammonium is typically below 300 nM except in Lake Erie, where it is severalfold higher and spatially variable; nitrate, on the other hand, is very high across the lakes but lowest in Erie due to biological uptake (42, 43). Few measurements of urea exist, but it can exceed ammonium (44) (see Data Set S2 at https://doi.org/10.6084/m9.figshare.15130350.v4). AOB (Nitrosomonadaceae) were observed across all lakes. In contrast, AOA (Nitrosopumilaceae) sequences exceeded 0.5% relative abundance only at the three deepest stations (SU08M, MI41M, ON55M), where the ratio of AOB to AOA ranged from 10:1 to 1:3. We found pronounced shifts in the dominant NOB across stations (Fig. 1), and all stations except those in Lake Ontario showed strong dominance (greater than 10-fold) of either “Ca. Nitrotoga” (family Gallionellaceae) or Nitrospira. Nitrospira was the only nitrite oxidizer detected in Lake Superior and the dominant nitrite oxidizer in parts of Lake Michigan (MI41M, MI18M). In contrast, “Ca. Nitrotoga” was the only nitrite oxidizer observed in Lake Erie and the dominant nitrite oxidizer in Lake Huron and at the shallowest station in Lake Michigan (MI27M). Within each taxon, a single 16S rRNA oligotype dominated the AOA, “Ca. Nitrotoga,” and Nitrospira, while several oligotypes of Nitrosomonadaceae shifted abundance across samples (Fig. S2), consistent with ecotypic diversity as discussed below.

We reconstructed 15 genomes of the AOB Nitrosospira, substantially expanding genome descriptions for this genus (45–47). Based on a phylogenomic tree, free-living Great Lakes Nitrosospira falls into two major clades, both of which are distinct from published species; each of these clades also includes metagenome-assembled genomes (MAGs) recovered from Lakes Biwa and Baikal, suggesting novel globally distributed freshwater lineages (Fig. S3). One clade, which we call NspGL1, has a highly reduced genome (median, 1.42 Mb) and low G+C content (40.7%). The second clade was resolved into three subclades (denoted NspGL2a, -2b, and -3) (Fig. S3) based on phylogeny and average nucleotide identity (ANI), all with small genome sizes of 1.45 to 1.68 Mb and 50% G+C content (Fig. 2; and see Data Set S3 at https://doi.org/10.6084/m9.figshare.15130350.v4). Compared to 86 reference Nitrosomonadaceae genomes, Great Lakes Nitrosospira genomes are not only smaller (median estimated complete genome size for the reference = 3.21 Mb, for GL = 1.45 Mb) (Table 1) but also have shorter intergenic spacers, fewer paralogs, fewer pseudogenes, and fewer sigma factors (Table 1 and Fig. S4; see Data Set S4 at https://doi.org/10.6084/m9.figshare.15130350.v4), consistent with genome streamlining to reduce resource demands (48). Based on short read mapping, these subclades are ecologically distinct: NspGL1 and NspGL2b—with the smallest genomes—are the dominant AOB in the upper oligotrophic lakes, while NspGL2a is abundant only in Lake Ontario and NspGL3 is abundant only in Lake Erie (Fig. 2). Hereafter, we refer to these subclades as ecotypes due to their phylogenetic and ecological divergence.

We next compared gene content between our Great Lakes Nitrosospira and 86 Nitrosomonadaceae reference genomes (see Data Set S5 at https://doi.org/10.6084/m9.figshare.15130350.v4). On average, Great Lakes Nitrosospira genomes encode far fewer two-component signal transduction systems (NspGL = 5 to 8, mean reference = 19), transposases (NspGL = 0 to 7, mean reference = 39), motility genes (NspGL = 0 to 4, mean reference = 52), pilus and secretion genes (NspGL = 2 to 9, mean reference = 27), and defense-related genes (NspGL = 4 to 11, mean reference = 39). They also lack functions related to biofilm formation such as polysaccharide matrix production (e.g., pel genes) and extracellular protein targeting (exosortase and PEP-CTERM motifs). Our 15 new Nitrosospira MAGs have high estimated completion (median, 98.6%), and therefore it is unlikely that these gene absences can be entirely attributed to incomplete assemblies. This overall picture of gene content in Great Lakes AOB contrasts with that of Nitrosospira isolates from soil (45, 47) and even of oligotrophic Nitrosomonas isolates (49) and is consistent with a passive planktonic lifestyle in extremely low-nutrient systems.

We next compared metabolic potential among Great Lakes AOB ecotypes to understand their ecological preferences for upper lakes (NspGL1, NspGL2b), Lake Ontario (NspGL2a), and Lake Erie (NspGL3). Surprisingly, all seven NspGL1 MAGs encode proteorhodopsin, a light-driven proton pump that supports bacterial energy production (50, 51). They also carry the genes necessary to synthesize its chromophore retinal, including those encoding 15,15′-beta-carotene dioxygenase (blh), lycopene cyclase (crtY), phytoene dehydrogenase (crtI), phytoene synthase (crtB), and GGPP synthase (crtE) (52, 53) (Fig. 3a). We also identified proteorhodopsin in a single-cell amplified genome representing NspGL1 from Lake Michigan (Fig. 3a) (99.8% ANI with NspGL1 MAGs), demonstrating that it is not an artifact of metagenome assembly. To our knowledge, this is the first example of a nitrifier with proteorhodopsin. All NspGL1 proteorhodopsins share residues H95, D127, and E138 along with a short beta-turn (G111-P116) between helices B and C, which are characteristic features of proteorhodopsin as distinct from sensory and other rhodopsins (54), and the presence of leucine at position 135 suggests green light tuning (55) (Fig. 3b). All of the genes in this module have highest similarity to homologs from Polynucleobacter but are flanked by Nitrosomonadaceae-like genes, suggesting recent horizontal gene transfer (Fig. 3a). The predicted NspGL1 proteorhodopsins cluster with Polynucleobacter, Methylopumilus, and other freshwater Betaproteobacteria in supercluster III as defined by MicRhoDE (56) (Fig. 3c). We compared the homologous genome regions in two highly similar MAGs from Lakes Biwa and Baikal (Fig. S5); these contigs lack the proteorhodopsin module but appear to flank a variable region where the contig assembly ends. A proteorhodopsin photosystem could support survival of NspGL1 in the presence of sunlight, which has been shown to inhibit ammonia oxidation (37, 57). In the upper lakes where NspGL1 is abundant, light penetration is high well below the thermocline in stratified periods (58), and deep-water taxa are seasonally advected to the surface by water column mixing (34). In addition to proteorhodopsin, NspGL1—but not the other three ecotypes of Great Lakes Nitrosospira—encodes a class I cyclopyrimidine dimer photolyase, which uses light energy to repair UV-induced DNA damage, and carries the catalase-peroxidase katG, suggesting that the NspGL1 ecotype is adapted to relatively shallow depths in the water column (Fig. 2).

Great Lakes Nitrosospira genomes carry a reduced, ecotype-specific complement of nitrogen metabolism genes compared to reference AOB (Fig. 2) (gene absences were verified as described in Materials and Methods). All are presumed to have the core ammonia oxidation enzymes ammonia monooxygenase and hydroxylamine dehydrogenase; these genes were assembled and binned as expected in some MAGs and were manually identified on short unbinned contigs in other cases (see Data Set S6 at https://doi.org/10.6084/m9.figshare.15130350.v4) (see Materials and Methods). Surprisingly, all Great Lakes Nitrosospira MAGs lack the copper protein nitrosocyanin, whose precise function is unknown but so far has been found in all described AOB except one member of the Nitrosomonas oligotropha clade (49). Based on the expanded set of genomes analyzed here, the lack of nitrosocyanin likely extends beyond the Great Lakes MAGs to closely related freshwater and marine strains, along with additional members of the N. oligotropha clade (Fig. S3); its absence may be related to the divergence of these clades. Only NspGL1 and NspGL2b encode NO-forming nitrite reductase (NirK), which confers nitrite tolerance (59); this result is surprising, given that these two clades dominate the upper lakes, where productivity and reduced N are lowest. None of the Great Lakes ecotypes encode NO reductase (NorCBQD), and NspGL1 lacks cytochrome P460 family proteins, both of which are common in AOB and implicated in nitrogen oxide metabolism (18, 49). Nitrogen acquisition is also distinct among Great Lakes AOB: NspGL1 lacks an apparent ammonium transporter but has urease structural and accessory genes (ureABCEFG) and a high-affinity urea transporter (urtABCDE). Further, all Great Lakes ecotypes encode a high-affinity amino acid transporter (livFGHM); these genes are rare (<5%) in reference genomes and could supply reduced nitrogen and/or organic carbon. Finally, NspGL1 and NspGL3 have genes for producing cyanophycin, an intracellular storage compound for nitrogen (47, 60). Together, the distinctive gene complements present in Great Lakes Nitrosospira illustrate the variability and adaptability of AOB gene content, even across a connected freshwater habitat.

As with nitrogen metabolism, carbon metabolism is also distinct between Great Lakes and reference AOB and among Great Lakes ecotypes (Fig. 2). Unlike most reference AOB, Great Lakes Nitrosospira AOB lack two key enzymes of the oxidative pentose phosphate pathway, glucose-6-phosphate dehydrogenase and 6-phosphogluconate dehydrogenase. All ecotypes except Erie-specific NspGL3 also lack genes for glycogen synthesis and degradation, suggesting that they are unable to store and access this carbon reserve. The key enzyme for carbon fixation, RuBisCO, has evolved several kinetically distinct forms whose distribution likely reflects ecological pressures (61). NspGL1 and NspGL3 both contain form IA RuBisCO, while NspGL2a and NspGL2b contain form IC RuBisCO (Fig. 2) (61, 62). NspGL1 genomes also possess an alpha carboxysome-like cso operon, similar to Nitrosomonas eutropha C91 (62), though our draft assembly lacks the expected carbonic anhydrase (csoS3/csoSCA). Carboxysome-associated RuBisCO may allow NspGL1, the ecotype most strongly adapted to energy and nutrient limitation, to more efficiently fix CO₂ by minimizing the wasteful oxygenation reaction and reducing the cellular nitrogen allocation to RuBisCO (61). The ranges of kinetic properties observed in other autotrophs for form IAq (found in NspGL3) and form IC (found in NspGL2a and NspGL2b) overlap, and therefore more work is needed to understand the fitness advantages, if any, that this RuBisCO diversity confers on Great Lakes nitrifiers.

We reconstructed three similar genomes (>99% ANI) of Nitrosarchaeum (NarchGL) (Fig. S6) from three separate samples (two from Lake Superior, one from Lake Ontario), consistent with our low observed 16S rRNA diversity for Thaumarchaeota. These NarchGL genomes are very similar (∼99% ANI) to two genomes reconstructed from Lake Baikal, located thousands of kilometers away (63). Their next closest relatives are also from freshwater environments, and phylogenetic clustering suggests that salinity is an important driver of divergence throughout the Nitrosopumilaceae (Fig. S6). As a group, the Thaumarchaeota tend to have smaller genomes, lower G+C content, higher coding density, and fewer paralogs and pseudogenes than nitrifying bacterial taxa; even within this group, NarchGL genomes fall below the 30th percentile in size and have significantly fewer paralogs than average (Table 1 and Fig. S4) (see Data Set S4 at https://doi.org/10.6084/m9.figshare.15130350.v4). By using our reconstructed genomes as probes for metagenomic read recruitment, NarchGL were detected in Lakes Superior, Michigan, and Ontario; they represented roughly one-third of ammonia oxidizers in the mid-hypolimnion of station SU08M (Fig. 2).

NarchGL genomes share nearly 90% of their predicted proteins with close relatives, including “Ca. Nitrosarchaeum limnia”; at the same time, they show distinctive patterns in gene content that pinpoint the key selective pressures of deep lakes (Fig. 2) (see Data Set S5 at https://doi.org/10.6084/m9.figshare.15130350.v4). All three NarchGL genomes encode urease and a urea transporter, implicating urea as a vital source of nitrogen for energy and/or biosynthesis. Consistent with phosphorus scarcity in much of the Great Lakes (64), NarchGL genomes encode high-affinity transport systems for phosphate and potentially phosphonates, though we did not identify a phosphonate lyase. In addition to both CRISPR/Cas enzymes cas1 and cas4, NarchGL genomes contain several phage proteins, suggesting that viral infection and integration events may be common. DNA photolyases, which have been found in epipelagic clades of marine Thaumarchaeota (7), are present in all genomes representing low-salinity Nitrosarchaeum including NarchGL, consistent with sunlight exposure due to high water clarity (58) and/or annual mixing in the Great Lakes. NarchGL also lack the common tRNA modification 4-thiouridylation (indicated by KEGG Orthology K04487 and Pfam PF02568-PF18297 [65]); we propose that the absence of this modification, which is susceptible to near-UV radiation (65), is also related to sunlight exposure.

Genomes of NarchGL reveal striking reductions in environmental sensing, response, and regulatory functions, relative to most other Nitrosarchaeum and Nitrosopumilaceae. NarchGL genomes encode 9 to 12 domains representing the general archaeal transcription factors TATA binding protein (TBP; Pfam PF00352) and transcription factor B (TFB; Pfam PF00382 and PF08271), compared to 21 in “Ca. Nitrosarchaeum limnia.” NarchGL genomes lack common domains found in two-component systems that transmit environmental signals to control gene expression or protein activity (Pfam domains PF02743, PF00672, PF00512, and PF00072; NarchGL = 0, “Ca. N. limnia” = 53 to 54 copies per genome). Further, they are depleted in ArsR family transcription factors (PF01022; 0 copies in NarchGL versus 2 to 3 in “Ca. N. limnia”), P-II proteins for regulation of nitrogen metabolism (PF00543; 1 copy per NarchGL genome versus 5 in “Ca. N. limnia”), and other potential regulatory domains (CBS PF00571, 5 in NarchGL versus 18 to 19 in “Ca. N. limnia”; USP PF00582, 1 in NarchGL versus 15 in “Ca. N. limnia”). This extremely limited regulatory capacity in NarchGL stands in sharp contrast to closely related “Ca. N. limnia” and instead parallels the oceanic minimalist “Ca. Nitrosopelagicus brevis” (9).

Despite the broad distribution of “Ca. Nitrotoga” in freshwater systems and beyond, only six genomes are available, derived from rivers heavily impacted by urban and agricultural influence, a wastewater treatment plant, and coastal sediment (27, 66, 67). Hence, the metabolic and phylogenetic diversity of this group is virtually unexplored. We reconstructed six new MAGs of “Ca. Nitrotoga,” which form two clusters with >99% ANI within each cluster and ∼97% ANI between clusters (NtogaGL1a and NtogaGL1b) (Fig. S3). These new “Ca. Nitrotoga” MAGs are far smaller than published genomes (median, GL = 1.44 Mb, reference = 2.61 to 2.98 Mb) and have shorter intergenic regions, fewer sigma factors, and fewer paralogs (Table 1 and Fig. S4) (see Data Set S4 at https://doi.org/10.6084/m9.figshare.15130350.v4), consistent with genome streamlining (48). They also have distinctive gene content (see Data Set S5 at https://doi.org/10.6084/m9.figshare.15130350.v4). They lack functions such as motility and chemotaxis, pilus biogenesis, and DNA repair (mutLS). Great Lakes “Ca. Nitrotoga” genomes also encode markedly fewer two-component systems for sensing and responding to environmental cues than four river-derived genomes (NtogaGL = 2 to 6 per genome versus 30 to 35 in four reference genomes; reference strain KNB has 7). Compared to reference genomes, NtogaGL genomes have fewer defense-related genes (restriction-modification, toxin-antitoxin, and CRISPR-Cas systems; mean, NtogaGL = 11 versus 39 for references), and transposases (mean, NtogaGL = 3 versus 19 for references). While incomplete assembly of hypervariable genome regions may explain some of these absences, the overall genome properties are consistent with a relatively stable low-nutrient environment and planktonic lifestyle.

The reduced genomes of NtogaGL1a/b help clarify core features of the genus “Ca. Nitrotoga,” along with accessory functions that may enable local adaptation in specific populations. To date, sequenced “Ca. Nitrotoga” genomes including NtogaGL1a/b encode similar electron transport pathways, including NADH dehydrogenase-complex I, succinate dehydrogenase-complex II, and alternative complex III, along with high-affinity cbb3-type cytochrome oxidases, suggesting adaptation to low-oxygen conditions. They also share the Calvin cycle for carbon fixation, a complete tricarboxylic acid (TCA) cycle, and an evolutionarily distinct nitrite oxidoreductase (NXR) from other NOB (27, 66, 67). All “Ca. Nitrotoga” genomes to date also share transporters for amino acids and peptides, potential sources of C and/or N. “Ca. Nitrotoga” can also potentially access reduced sulfur compounds for energy via sulfite dehydrogenase, suggesting metabolic flexibility beyond nitrite oxidation.

Beyond these similarities, the small genomes of NtogaGL1a/b are distinct from previously described “Ca. Nitrotoga” in many ways. NtogaGL1a/b lack NiFe hydrogenase to use hydrogen as an energy source. They also lack nitrogen metabolism functions, including assimilatory nitrite reductase (nirBD) and nitrite reductase to NO (nirK). Based on gene content, NtogaGL1a/b appear unable to use hexoses like glucose, since they lack the glycolytic enzyme phosphofructokinase and the Entner-Doudoroff pathway, similar to Nitrobacter winogradskyi (68). Consistent with this, they also lack genes for storage and breakdown of glycogen. All but one of the NtogaGL1a/1b genomes encode cyanate lyase (cynS), which is found in other NOB but not in “Ca. Nitrotoga” to date (25, 69, 70). The cynS gene, adjacent to glnK-amtB for ammonium sensing and transport, likely functions in N assimilation, as recently described for Nitrospinae (71). While cyanase has been shown to mediate reciprocal feeding between some NOB and ammonia oxidizers (25), it remains to be seen whether such an interaction occurs in the free-living (<1.6-μm) size fraction and dilute environment sampled here. Notably, cyanase from NtogaGL1a/1b, along with predicted Nitrospirae proteins from Lake Baikal and soil, form a distinct phylogenetic cluster from most nitrifier cyanase proteins observed to date (Fig. S7).

The two ANI-based clusters we detected, NtogaGL1a and NtogaGL1b, appear to be phylogenetically and ecologically distinct ecotypes. Based on short read mapping, NtogaGL1b dominates Lake Erie, while NtogaGL1a dominates all other “Ca. Nitrotoga”-containing samples (Fig. 4). We found several metabolic genes that differentiate the two ecotypes. Lake Erie-specific NtogaGL1b genomes share a region encoding thiosulfate dehydrogenase (tsdA), cytochromes, transport of sulfur-containing compounds, lactate dehydrogenase (ldh), a two-component system, and a Crp-family transcription factor (Fig. S8). This region may be involved in oxidizing thiosulfate as an energy source and sensing and responding to redox changes that accompany seasonal hypoxia in Lake Erie. The corresponding region in NtogaGL1a encodes an integrase and photolyase, consistent with greater DNA photodamage in the more transparent waters of Lakes Michigan, Huron, and Ontario, where NtogaGL1a is abundant.

We reconstructed six closely related genomes of Nitrospira (∼99% ANI) (Fig. S9), representing the predominant NOB throughout Lake Superior and in parts of Lakes Michigan and Ontario (Fig. 4) (see Data Set S1 at https://doi.org/10.6084/m9.figshare.15130350.v4). These genomes, which we refer to as NspiraGL, fall within lineage II (Fig. S9), which is broadly distributed across soil, freshwater, and engineered habitats (20); however, genome analyses to date have focused on strains from wastewater and engineered systems, leaving major blind spots. NspiraGL share core features of Nitrospira metabolism, including a periplasm-facing NXR that is advantageous under substrate-limiting conditions, multiple cytochrome bd-like oxidases, and the reverse TCA cycle for carbon fixation (69). However, as with “Ca. Nitrotoga,” the Nitrospira genomes we reconstructed in the Great Lakes are markedly smaller than published reference genomes (median for NspiraGL = 1.83 Mb, median for reference = 3.72 Mb), with higher coding density and fewer paralogs, sigma factors, and pseudogenes (Fig. S4) (see Data Set S4 at https://doi.org/10.6084/m9.figshare.15130350.v4), consistent with genome streamlining theory (48). Compared to 75 lineage II Nitrospira reference genomes, NspiraGL have reduced capacity for environmental sensing (two-component systems: NspiraGL = 7, mean reference = 26), transport (NspiraGL = 76 to 83, mean reference = 140), defense (NspiraGL = 7 to 8, mean reference = 26), and transposition (NspiraGL = 0 to 2, mean reference = 15) and lack pilus or flagellar motility (see Data Set S5 at https://doi.org/10.6084/m9.figshare.15130350.v4). NspiraGL encode just five sigma factors, compared to 18 in Nitrospira moscoviensis. Further, NspiraGL genomes encode a single NXR, while N. moscoviensis carries five copies that are differentially regulated (26, 72). NspiraGL also lack the glnE gene for glutamine synthetase (GS) adenylyltransferase, suggesting that GS activity is not repressed by this mechanism. Together, these features suggest limited regulatory and ecological flexibility, consistent with a relatively constant, oligotrophic environment.

Compared to other Nitrospira, NspiraGL exhibit limited energetic flexibility but can access diverse nitrogen sources (Fig. 4). We predict that NspiraGL are unable to grow on hydrogen or formate as alternative energy sources (23, 26), as they lack NiFe-hydogenase and formate dehydrogenase. The glycolysis and oxidative TCA cycles appear to be incomplete, lacking phosphofructokinase and citrate synthase, respectively; this suggests a limited capacity for organic carbon utilization. NspiraGL lack nirK, encoding NO-forming nitrite reductase, which is found in a majority of reference genomes. To obtain nitrogen for biosynthesis, NspiraGL encode a high-affinity nitrate/nitrite/cyanate transporter (nrtABC), assimilatory nitrite reductase (nirA), and cyanase (cynS), along with amt family ammonium transporter. Although none of the NspiraGL MAGs include urease (ureCBA), one does contain urease accessory proteins (ureEFGD) and two contain a urea transporter (urtABCD), suggesting incomplete assembly of the urea utilization pathway. As with “Ca. Nitrotoga,” we suggest that cyanase, along with urease where present, functions in nitrogen assimilation rather than cross-feeding, given the dilute environment and free-living planktonic cells.

Beyond energy, carbon, and nitrogen metabolism, we discovered striking differences between NspiraGL and reference Nitrospira related to DNA repair. NspiraGL encode two additional photolyase-related proteins, along with a class I cyclopyrimidine dimer (CPD) photolyase found in most reference Nitrospira taxa (Fig. 5). Photolyases use blue light energy to repair DNA lesions caused by UV radiation (73). The two additional genes in NspiraGL are adjacent and share best hits with Betaproteobacteria, suggesting recent horizontal transfer (Fig. S10). One likely encodes an FeS-BCP photolyase, which repairs (6-4) dipyrimidine lesions (74, 75). The other shares an FAD-binding domain with photolyases, but the C-terminal region has no recognizable domains (Fig. 5). This protein is widespread in aquatic bacteria and has not been functionally characterized, though an actinobacterial homolog was suggested to be involved in light sensing and regulation (76). Beyond photolyases, NspiraGL also encode uracil-DNA glycosylase (UNG), which removes misincorporated uracil from DNA. Uracil results from deamination of cytosine, which can occur spontaneously or be induced by NO (77). In addition to the photolyases and UNG that repair DNA lesions, NspiraGL encode translesion DNA polymerase V (umuCD), which enables replication to proceed past lesions. Together, these genes indicate that members of Nitrospira in the Great Lakes experience significant DNA damage, including UV-induced damage that also requires light for the repair process, in hypolimnion waters with high transparency (58) and/or during seasonal mixing.

Other major differences between NspiraGL and reference Nitrospira genomes are related to reactive oxygen species (ROS). Surprisingly, despite their oxic habitat, NspiraGL lack superoxide dismutase (SOD), monofunctional catalase (katE), and bacterioferritin, which limits the Fenton reaction by sequestering free iron. However, all six NspiraGL MAGs, but few reference genomes (7% of 75), have recently acquired bifunctional catalase-peroxidase katG; interestingly, we also observed katG in Great Lakes “Ca. Nitrotoga” and Nitrosospira (Fig. 2 and 4). The absence of SOD suggests that NspiraGL does not produce damaging levels of endogenous superoxide, perhaps because NspiraGL lack the major respiratory and nonrespiratory flavoproteins that produce ROS in other SOD-containing Nitrospira taxa (78). Unlike superoxide, H₂O₂ can cross membranes and is known to be produced by both photooxidation of dissolved organic matter and dark heterotrophic activity (79). The lakes where NspiraGL dominate have high water clarity (58) and low productivity and are fully oxic, consistent with abiotic photochemistry as the primary source of exogenous ROS; this stress may have selected for katG as a defense. NspiraGL also lack cytochrome c peroxidase, which is found in 70 of 75 reference genomes; this protein is proposed to function in anaerobic respiration of H₂O₂ (80), and therefore its absence in NspiraGL is consistent with a constant oxic environment. Together, these results indicate that members of Nitrospira in the Great Lakes face distinct ROS pressures that have shaped their gene content.

The Laurentian Great Lakes harbor nitrifiers that are phylogenetically related, but markedly different in genome size and functional capacity, from their well-studied relatives inhabiting wastewater systems, soils, and even other freshwater systems. By examining the entire nitrifier assemblage at once, we detected common features across taxa that illuminate the selective pressures faced by microbes in deep lakes. All the lineages we describe show small genome sizes (1.3 to 1.7 Mb), reduced capacity for environmental sensing and response, and adaptation to a passive (i.e., nonmotile) planktonic lifestyle, features which have not been previously associated with AOB, Nitrospira, and “Ca. Nitrotoga.” Within the AOB Nitrosospira, we found ecotypes with a gradient of genome reduction that maps onto their habitats’ trophic gradient: from NspGL1 (1.4 Mb, low GC, upper lakes) to NspGL2b (1.5 Mb, upper lakes) to NspGL2a (1.6 Mb, Lake Ontario) to NspGL3 (1.7 Mb, Lake Erie) (Fig. 2). The thaumarchaeal NarchGL have a markedly reduced regulatory capacity like the open ocean strain Nitrosopelagicus brevis (9). The NOB NspiraGL have genomes 50 to 60% smaller than the genomes of described Nitrospira taxa and dominate the deeper more oligotrophic basins, while “Ca. Nitrotoga” favors shallower, more productive basins. The emergence of Lake Erie-specific ecotypes of both Nitrosospira (NspGL3) and “Ca. Nitrotoga” (NtogaGL1b) demonstrates how distinct this habitat is compared to the other lakes. Importantly, our findings here represent planktonic cells in the smallest size fraction (<1.6 μm); it is likely, especially in Lake Erie, that particle-associated nitrifiers may be abundant and genetically distinct.

Nitrifiers inhabiting the transparent waters of the upper Great Lakes show distinctive adaptations to light, including diverse photolyases, ROS detoxification, and even proteorhodopsin. This discovery is surprising, given that nitrifiers are rare in the surface mixed layer of the Great Lakes (Fig. S1) and that photoinhibition of ammonia oxidation and nitrifier growth is well documented (37, 40, 57). We propose that proteorhodopsin could be used to augment energy metabolism when ammonia oxidation is photoinhibited and/or ammonia oxidation is substrate limited. Water clarity has increased over the past several decades in Lakes Michigan and Huron, now surpassing that of Lake Superior (58). High light penetration along with seasonal mixing likely exposes deep-water cells to damaging levels of light and oxidative stress. Future cultivation and physiological studies should examine photoinhibition and potential phototrophy in Great Lakes nitrifiers.

Our work unveils new clues about the ecological and evolutionary potential of nitrifiers in their natural freshwater habitat. This collective nitrifier diversity undoubtedly influences the cycling of carbon and nitrogen across this ecosystem, and future work will explore the differential contributions to nitrification by the distinct lineages we described here. Understanding what controls the diversity of nitrifiers and other key functional groups, and the consequences of this diversity for biogeochemistry, are essential for forecasting the effects of rapid environmental change across the large lakes of the world (e.g., see reference 81) and predicting impacts on the critical ecosystem services they provide (82).

Water samples were collected from the Laurentian Great Lakes aboard the R/V Lake Guardian, during the biannual Water Quality Surveys conducted by the U.S. EPA Great Lakes National Program Office (83). Station information is provided in Data Set S8 at https://doi.org/10.6084/m9.figshare.15130350.v4. Data presented here were collected in April and August 2012. Samples were collected using a conductivity-temperature-depth (CTD) rosette sampler (Sea-Bird Scientific) at the surface (2 m), deep chlorophyll maximum (if present), the mid-hypolimnion (depths ranging from 19 m in Lake Erie to 200 m in Lake Superior) (see Data Set S1 at https://doi.org/10.6084/m9.figshare.15130350.v4), and near the bottom of the water column (10 m above the lake bottom at most stations, 1 m above bottom at shallow stations). For each sample, 5 to 8 L of water was prefiltered through a GF/A glass fiber filter (Whatman 1820-047; nominal pore size, 1.6 μm) to exclude eukaryotic phytoplankton and particle-associated microbes, and cells were collected on 0.22-μm Sterivex filters (Millipore SVGP01050). Filters were stored at −80°C. For dissolved nutrient analysis, 0.22-μm filtrate was collected in 125-mL acid-clean high-density polyethylene (HDPE) bottles (Nalgene) and stored at −20°C. Samples for single-cell amplified genomes (SAG) were collected in August 2014. For each sample, 1 mL of raw water was incubated with 100 μL of glycerol-TE buffer (20 mL 100× Tris-EDTA [TE], pH 8, plus 100 mL glycerol plus 60 mL water; final concentrations after sample addition are 10 mM Tris, 1 mM EDTA, 5% glycerol) for 10 min in the dark and then flash frozen in liquid nitrogen and stored at −80°C until processing.

CTD profiles, water chemistry, and chlorophyll a data were collected by the U.S. EPA according to standard protocols (84) and retrieved from the Great Lakes Environmental Database (https://cdx.epa.gov/) for 2012 and 2013. In addition, we measured dissolved nitrogen species from August 2013 samples. Ammonium concentrations were measured using the OPA method in a 96-well plate (85). Nitrate and nitrite concentrations were measured using the Griess reaction method in a 96-well plate (86). Urea concentrations were measured in a 24-well plate using a colorimetric reaction (87).

The full 16S rRNA amplicon data set was described by Paver and colleagues (34). Here, we focus on data from the V4‐V5 region (primers 515F‐Y and 926R [88]), collected in 2012 in tandem with metagenome samples from select stations. We classified sequences using the Silva v.132 database (89) and the method of Wang et al. (90) as implemented by mothur (91). Sequences classified to each detected family of nitrifiers (ammonia oxidizer families Nitrosomonadaceae and Nitrosopumilaceae; nitrite oxidizer families Gallionellaceae and Nitrospiraceae) with a mothur-assigned confidence score above 90 were delineated into taxonomic units using minimum entropy decomposition with a minimum substantive abundance of 10 (92).

One station per lake in Lakes Superior, Michigan, Huron, and Ontario and two stations in Lake Erie were selected for metagenome sequencing. Spring 2012 metagenome samples were collected from the surface, and summer 2012 metagenome samples were collected from the mid-hypolimnion (depths listed in Data Set S1 at https://doi.org/10.6084/m9.figshare.15130350.v4). DNA was extracted using a modified phenol-chloroform extraction protocol (34), and libraries were prepared according to the Illumina TruSeq protocol. Samples from spring 2012 were sequenced at the Joint Genome Institute using Illumina HiSeq (2 × 150 bp). Samples from summer 2012 were sequenced at the University of Chicago Functional Genomics Core Facility using Illumina HiSeq 2500 (2 × 250 bp).

To confirm the presence of proteorhodopsin, we analyzed a single-cell amplified genome from Nitrosospira collected from Lake Michigan and sequenced by the Joint Genome Institute. Quality filtered reads were downloaded from Joint Genome Institute (JGI) IMG/ER and normalized using bbnorm.sh with a target of 100 and a mindepth of 2. Normalized reads were assembled using SPAdes 3.1.11 in single-cell mode (93) with flags –sc and –careful. Resulting scaffolds were annotated identically to MAGs as described below.

Raw reads for spring surface samples were quality controlled at the Joint Genome Institute, using bbduk.sh for adapter trimming (ktrim = r, minlen = 40, minlenfraction = 0.6, mink = 11, tbo, tpe, k = 23, hdist = 1, hdist2 = 1, ftm = 5) and quality filtering (maq = 8, maxns = 1, minlen = 40, minlenfraction = 0.6, k = 27, hdist = 1, trimq = 12, qtrim = rl). Raw reads for summer hypolimnion samples were adapter trimmed, quality filtered, and interleaved using bbduk (parameters: ktrim = r, mink = 8, hdist = 2, k = 21, forcetrimleft = 10, forcetrimright = 199, minlen = 150) using BBTools suite version 35.74 (https://sourceforge.net/projects/bbmap/). Separate assemblies of quality-filtered reads were carried out for each metagenome using metaSPAdes 3.1.11 –meta mode using default k sizes of 21, 33, and 55 (94). To enable binning based on sequence coverage, forward and reverse reads were merged using bbmerge in BBtools, using qtrim2 = r trimq = 10,13,16 and adapter = default. Merged short reads were then mapped onto each assembly using bowtie2 2.2.9 in –sensitive mode (95), and this coverage information was used to bin assembled contigs. Binning was performed using MetaBAT2 2.12.1 (96), Binsanity 0.2.6.3 (97), and CONCOCT 1.0.0 (98) using default parameters. The resulting bins were scored, aggregated, and dereplicated using DAS_Tool 1.1.1 (99), followed by manual curation using Anvi’o 4.0 (100). We assessed genome completion and contamination of manually curated bins using CheckM 1.1.0 lineage_wf (101), and all new MAGs presented here are greater than 70% complete with less than 10% contamination (see Data Set S3 at https://doi.org/10.6084/m9.figshare.15130350.v4). Potential nitrifiers were screened by searching for ammonia monooxygenase, hydroxylamine oxidoreductase, and nitrite oxidoreductase genes within reconstructed genomes using blastp 2.5.0 (102). For bins where any of these genes were detected, we identified bacterial single-copy core genes (103) or archaeal single-copy core genes (104) using HMMER (105), as implemented in Anvi’o. Single-copy core genes were queried against proteins predicted from bacterial and archaeal genomes in RefSeq (NCBI) (106), and taxonomic identity of these core genes was ascertained based on a least common ancestor approach using a 0.1% window around the bit score of the best hit using KronaTools 2.7.1 (107). Taxonomic assignment was further validated using GTDB-tk 1.0.0 (108). Grouping of MAGs into clades and subclades based on ANI was carried out using fastANI 1.1.0 (109). Genome characteristics for each genome group were calculated as the median of those values for the group. Estimated complete genome size was calculated for MAGs and for references in the pangenome analysis using CheckM (101) completion and contamination, as follows: estimated = actual × [(1 − contamination)/completion]. To quantify the abundance of each clade/ecotype across samples, we used competitive mapping of merged short reads using bowtie2 in sensitive mode against all nitrifier MAGs, summing up the mapped read counts across all MAGs in a given clade/ecotype and dividing by the total mapped nitrifier reads in a sample; these values are shown in bubble plots (Fig. 2 and 4).

Reference genomes were obtained from GenBank (accession numbers listed in Data Set S4 at https://doi.org/10.6084/m9.figshare.15130350.v4). The full pangenome analyses included all the genomes listed therein, but we only report results from the subset of genomes most closely related to our MAGs. This subset consists of 86 Nitrosomonadaceae, 5 “Ca. Nitrotoga,” 78 Nitrosopumilaceae within Thaumarchaeota, and 75 Nitrospira genomes that fall within lineage II. Reference genomes were treated consistently with GL MAGs, with de novo gene calling by prodigal 2.6.3 (110) via Anvi’o. Unless otherwise noted, default settings were used for all software. Genes were annotated using InterProScan 5.30–69.0 (111), GhostKOALA (112), and eggnog-mapper 1.0.3 against the bactNOG database (113). Gene cluster analysis was carried out using the Anvi’o pangenome pipeline (114), using blastp to determine sequence similarity, ITEP to eliminate weak similarity (115), and MCL to cluster, using a minbit of 0.5, MCL inflation of 2, and minimum gene occurrence of 1 (116). Sigma factors were tallied by identifying gene clusters annotated with the following PFAMs: PF00309, PF03979, PF00140, PF04542, PF04539, PF04545, and PF08281. Pseudogene counts were retrieved where available from NCBI PGAP annotated genomes (117). Paralog counts are reported as the number of gene clusters with more than one gene per genome. Intergenic spacers were calculated using bedtools complementBed function (118). Coding fraction is defined as the summed length of all protein-coding genes divided by the estimated total genome length. Prokka 1.14.5 (119) was used to generate GenBank format files from MAGs and SAGs, and genoPlotR 0.8.9 (120) was used to generate initial gene neighborhood maps.

The NspGL1 proteorhodopsin sequence was inserted into the MicRhoDE rhodopsin tree using pplacer (121) through the MicRhoDE Galaxy pipeline (56). We then constructed a more targeted phylogenetic tree using aligned reference sequences of supercluster III from MicRhoDE, filtered to exclude fragments shorter than 220 amino acids. To this alignment, we added NspGL1 sequences using MAFFT 7.310 (122) along with high-similarity sequences from NCBI nr that were not present in MicRhoDE. The tree was inferred using RaxML 8.2.12 with model PROTGAMMALG (123). The tree was visualized in iTOL (124), and more distant clusters were collapsed for clarity.

A cyanase phylogenetic tree was created using sequences drawn from querying NtogaGL cyanase against NCBI nr using blastp, as well as sequences from references 2, 25, and 125. Sequences were aligned using MAFFT (122), and the tree was inferred using RaxML 8.2.12 with model PROTGAMMALG (123). The tree was visualized in iTOL (124), and branches were colored based on the taxonomy of the parent genome.

Photolyase-related proteins in GL MAGs were identified by searching for the following features: KEGG Orthology K01669, NCBI Clusters of Orthologous Genes COG0415, Pfams PF03441, PF00875, PF04244, Superfamilies SSF48173, and SSF52425. Reference proteins (n = 56) spanning the previously defined families of photolyases and cryptochromes (126) were obtained from UniProt, along with aquatic bacterial sequences described by Maresca and colleagues (76). The reference sequences were aligned using MAFFT (122), and sequences from GL MAGs were added using the MAFFT –addfragments option. The tree was estimated using IQ-TREE 2 1.6.11 (127) and visualized using iTOL (124).

Nitrospirae, Thaumarchaeota, Gallionellaceae, and Nitrosomonadaceae genomes were downloaded from GenBank (NCBI) (128) and included in the phylogenomic trees for their respective family. Phylogenomic analyses were carried out within Anvi’o. Briefly, single-copy core genes were extracted as described above, individually aligned at the protein level using muscle (129), and concatenated for each genome. Concatenated alignments were trimmed using Gblocks 0.91b (130) and analyzed by RAxML 8.2.12 (123) to create a phylogenetic tree using the PROTGAMMALG model and 50 bootstraps. Trees were visualized in iTOL (124).

We used several approaches to validate the presence of proteorhodopsin in assembled Nitrosospira genomes, to rule out the possibility of chimeric assemblies from different species. We note that proteorhodopsin-containing contigs were independently assembled and binned together with core Nitrosospira contigs from seven different samples (i.e., each sample was assembled and binned separately, rather than coassembled). In five of seven cases, proteorhodopsin and retinal biosynthesis genes were assembled together with core Nitrosospira genes on the same contig. To rule out a systematic reproducible error in assembly and/or binning, we compared these seven MAGs to a single-cell amplified Nitrosospira genome (SAG) from Lake Michigan, obtained as part of another project with the JGI. This SAG was processed through JGI’s standard decontamination pipeline and manually investigated to ensure lack of contamination. We found no evidence of contaminating core genes, as all core genes had best hits to either Nitrosospira or more generally Nitrosomonadaceae in nr. SAG contigs were matched to homologous contigs from NspGL1 MAGs to determine if any SAG contigs were unique using FastANI 1.1.0 (109) with –visualize flag. All contigs from this Nitrosospira SAG were found within an NspGL1 MAG. Bandage 0.8.1 (131) was used to manually inspect the assembly graph around the contig that contained the NspGL1 Nitrosospira proteorhodopsin to ensure that the assembled contig did not represent a chimeric contig or inappropriate scaffolding. We verified that a single, unique path exists from the beginning to the end of the NspGL1 contig containing proteorhodopsin (Fig. 3). Further, we verified that consistent coverage across this contig existed by mapping short reads from the original sample using bowtie2 (95) and viewing results using Integrated Genomics Viewer 2.7.0 (132). A closely related assembly of the same genomic region from Lake Biwa did not show evidence of proteorhodopsin; to confirm this difference between the Lake Biwa and Great Lakes MAGs, we mapped reads from Lake Biwa (133) (BioProject PRJDB6644) onto the assembled contig described above using bowtie2 (95). This analysis demonstrated that while a large fraction of the NspGL1 contig in question recruited reads from Lake Biwa at high identity (98 to 99%), starting upstream of proteorhodopsin and retinal biosynthesis, this contig no longer recruited reads from Lake Biwa.

Despite recovery of 15 high-completion MAGs in NspGL1/2a/2b/3, many of these MAGs lacked key nitrification genes in amo and hao operons. This was largely due to the fact that amo and hao operons were often assembled on small contigs below the minimum size cutoff we imposed for binning contigs. Difficulty in assembling these contigs was likely due in part to the several amo and hao operons with extremely high identity to one another in each genome, a phenomenon which has been observed in other Nitrosospira genomes (18). Manual assembly graph inspection with Bandage (131) supported this hypothesis, as did assessment of abundance of short reads associated with amo operons from NspGL and comparison of abundance of short reads associated with core gene rpoB from NspGL, using ROCker (134). Still, an exemplar MAG from at least one representative of each ecotype (NspGL1/2a/2b/3) was found with both amo and hao operons. Further, manual inspection of unbinned contigs confirmed that amo and hao operons existed on contigs in every sample from which a MAG for a particular ecotype was recovered. That is, for every time that an NspGL1 MAG was recovered from a sample, we were able to determine that an amo and hao operon which could be affiliated with NspGL1 existed, even if it was not correctly binned. Affiliation for these unbinned key nitrification genes was carried out by alignment of amoAB and haoAB sequences to amoAB and haoAB sequences correctly binned in NspGL ecotypes. This process was also carried out for two NtogaGL1a MAGs for nxrAB, which were poorly assembled in those two samples. Data Set S6 (at https://doi.org/10.6084/m9.figshare.15130350.v4) summarizes the presence of genes related to nitrification and nitrogen metabolism across all our MAGs.

Metagenome-assembled genomes typically comprise tens or even hundreds of contigs, and this fragmented nature makes it impossible to say with certainty whether a particular gene is truly absent. To substantiate our claims of gene absence based on MAGs, we used several lines of evidence. First, we note that our MAGs have high estimated completion (median of 96.4%, mean of 94.3%), based on the presence of universal core gene markers. Second, for all new lineages described here except NspGL3, we assembled multiple similar MAGs independently from different samples, and we inferred gene absences only if the absence was replicated in multiple assemblies. Together, these two factors provide strong support for cases where a missing gene would be expected to occur in a region of predominantly core genes; however, these factors are less informative for cases where a missing gene might occur in a genomic island, because we have no way of assessing the completion of regions lacking core genes, and islands tend to have systematic poor assemblies across samples. A third line of evidence that we considered is chromosome organization: if a single gene is deleted from an otherwise conserved region of synteny, then this deletion should be apparent in a gene neighborhood diagram (e.g., see Fig. S5, S7, S8, and S10 in the supplemental material). Unfortunately, in many cases, our MAGs are too dissimilar from reference genomes and share little synteny with them, so this approach is not always informative.

We used a fourth approach based on quantitative analysis of short reads to verify gene absences. If a suspected missing gene were actually present in the population, but failed to assemble and/or bin with the rest of the genome, then it should be detectable in the unassembled short reads. The frequency of a gene in the population can be estimated from its abundance in the short reads, compared to the abundance of core marker genes in the short reads. We implemented this approach as follows. We searched unassembled short reads for each gene of interest that we identified as absent from MAGs (e.g., nitrosocyanin) using tblastn. Short reads with significant similarity were then filtered by best-hit taxonomy to the appropriate nitrifier group (i.e., Nitrosomonadaceae, “Ca. Nitrotoga,” Nitrospira, Thaumarchaea). These filtered short reads were enumerated and length normalized [1,000 × (number of short reads/length of the target gene of interest)]. The same procedure was repeated for genes expected to be present in every cell (e.g., amoAB, hao, nxrAB, ribosomal protein genes) for comparison. If a putative missing gene (based on MAGs) has near-zero detection in the short reads, we can be confident that the gene is truly missing (or has undetectable sequence similarity, or was so recently acquired from another lineage that its best hit points to a different taxon). In contrast, if a putative missing gene is detected in the short reads, then the gene may be present in genomes related to our MAGs but was unassembled/unbinned, or the gene may be present in another lineage of nitrifiers that is not represented by our MAGs. Short read-based quantification of select genes is presented and described in Data Set S7 and Supplemental Text (both at https://doi.org/10.6084/m9.figshare.15130350.v4).

All statistical comparisons were carried out in R version 3.5.3 (135), and plots were generated using ggplot2 3.2.0 (136). Code and data files are available at bitbucket.org/greatlakes/gl_nitrifiers.

The metagenome-assembled genomes presented here are available via NCBI BioProject PRJNA636190. 16S rRNA data are available at NCBI BioProject PRJNA591360. Metagenomes sequenced by JGI are available at https://genome.jgi.doe.gov under project ID 1045056, 1045059, 1045062, 1045065, 1045068, and 1045071. The single-cell amplified genome is available at http://img.jgi.doe.gov/ under IMG Genome ID 3300033241. Raw reads are available in NCBI SRA (SRR14240538–SRR14240543) or through JGI with the project IDs listed above.