Rokubacteria in Northern Peatlands: Habitat Preferences and Diversity Patterns

Rokubacteria is a phylogenetic clade of as-yet-uncultivated prokaryotes, which are detected in diverse terrestrial habitats and are commonly addressed as members of the rare biosphere. This clade was originally described as a candidate phylum; however, based on the results of comparative genome analysis, was later defined as the order-level lineage, Rokubacteriales, within the phylum Methylomirabilota. The physiology and lifestyles of these bacteria are poorly understood. A dataset of 16S rRNA gene reads retrieved from four boreal raised bogs and six eutrophic fens was examined for the presence of the Rokubacteriales; the latter were detected exclusively in fens. Their relative abundance varied between 0.2 and 4% of all bacteria and was positively correlated with pH, total nitrogen content, and availability of Ca and Mg. To test an earlier published hypothesis regarding the presence of methanotrophic capabilities in Rokubacteria, peat samples were incubated with 10% methane for four weeks. No response to methane availability was detected for the Rokubacteriales, while clear a increase in relative abundance was observed for the conventional Methylococcales methanotrophs. The search for methane monooxygenase encoding genes in 60 currently available Rokubacteriales metagenomes yielded negative results, although copper-containing monooxygenases were encoded by some members of this order. This study suggests that peat-inhabiting Rokubacteriales are neutrophilic non-methanotrophic bacteria that colonize nitrogen-rich wetlands.


Introduction
The Rokubacteria is one of the deeply branching lineages of prokaryotes, which are often addressed as members of the rare biosphere [1][2][3]. The existence of this bacterial group was first recognized by Lipson and Schmidt [4], who retrieved the corresponding 16S rRNA gene sequences (GenBank accession numbers AY192281 and AY192282) from an alpine meadow of the Colorado Rocky Mountains. These sequences clustered together with two other environmental 16S rRNA gene sequences (AF234132 and AF428647), which were earlier obtained from an Australian arid soil and a lake sediment in China, and formed a deeply rooting, phylum-level bacterial lineage, which was originally referred to as the Candidate phylum SPAM (for Spring Alpine Meadow) [4]. The name Candidatus Rokubacteria was proposed in 2016 for the phylogenetic lineage defined by the metagenome bin CSP1-6, which was assembled from the sediment sample of the Rifle DOE Scientific Focus Area in Colorado [2]. This lineage was equally divergent (82% 16S rRNA gene sequence similarity) from the Deltaproteobacteria, Firmicutes, Nitrospirae and the Candidate phylum NC10. The fact that Candidatus Rokubacteria represents the same phylogenetic lineage as SPAM was noticed by Becraft et al. [5]. Notably, in treeing analyses that imply normalization of evolutionary distances, Candidatus Rokubacteria formed a common group with the candidate phylum NC10 [6,7]. Therefore, based on the results of the comparative genome analysis, this bacterial group was recently defined as the order-level lineage, Rokubacteriales, within the phylum Methylomirabilota, the class Methylomirabilia [7]. The second order within this class is the Methylomirabilales (also referred as NC10), which accommodates bacteria that couple anaerobic oxidation of methane with the reduction of nitrite to dinitrogen [8].
Despite the wide distribution of Rokubacteria, all attempts to isolate these microorganisms to date remain unsuccessful. All insights into the characteristic features of these enigmatic organisms were made by the cultivation-independent molecular approaches, including single-cell genomics and metagenomics [3,5]. These analyses showed that rokubacterial genomes are large (6)(7)(8) and have high GC content. At the same time, the cell sizes of these bacteria assessed using the fluorescence-activated cell sorting are small, 0.3-0.4 µm in diameter [5]. The presence of large genomes in small cells may imply extensive DNA packaging or dormancy. Notably, Rokubacteria are characterized by an intriguing high genomic heterogeneity among individuals, with no environments identified to date with near-clonal populations. This feature may present certain challenges to future studies [5].
Several recent studies have illuminated the possible metabolic and ecological properties of Rokubacteria [3,16,18]. Anantharaman and co-authors discovered that these bacteria contain some of the earliest evolved dissimilatory sulfite reductases (DsrAB) [18]. Rokubacterial DsrAB-encoding genes represent a novel deep-branching lineage on the tree among the sulfate/sulfite-reducing microorganisms. Besides dsrAB genes, Rokubacteria also possess apr, sat, and qmo genes that are required for sulfate reduction or sulfite oxidation [18]. In addition, a number of novel putative isopropanol dehydrogenases affiliated with the Rokubacteria were identified while studying wetland sediments [16]. Therefore, it cannot be ruled out that Rokubacteria may utilize ethanol and/or isopropanol as electron donors, which suggest their potential contribution to alcohol cycling in wetlands [16]. The most intriguing hypothesis regarding the metabolic capabilities of the Rokubacteria, however, was published by Kroeger and coauthors, who assembled and examined a rokubacterial metagenome-assembled genome (MAG), named Amazon-R-15-13, from Brazilian Amazon rainforest soil [3]. One of the unique features of Amazon R-15-13 was the presence of a gene operon with a structural similarity to that of particulate methane monooxygenase (pMMO), a key enzyme of methanotrophic metabolism. The corresponding enzyme from Amazon R-15-13, however, clustered together with monooxygenases from two non-methanotrophic actinobacteria, Nocardioides luteus and Streptomyces theroautotrophicus, and revealed only a very distant relationship to true pMMO from characterized methanotrophic bacteria.
This study was initiated in order to examine the diversity patterns of the Rokubacteriales in boreal wetlands and to verify the hypothesis regarding the presence of methanotrophic capabilities in these bacteria in incubation experiments. All currently available metagenomes of the Rokubacteriales were also examined for the presence of pMMO-encoding genes as well as the genes encoding evolutionary related copper-containing monooxygenases.

Analysis of Microbial Diversity in Boreal Wetlands
The comparison of microbial diversity patterns in two types of boreal wetlands, acidic razed bogs and neutral eutrophic fens, was performed using the 16S rRNA gene sequence dataset retrieved by Dedysh et al. (2021) and deposited in Sequence Read Archive (SRA) under the accession numbers SRR11280489-SRR11280524 (Bioproject PRJNA610704). The set of examined wetlands included four raised bogs (Shichengskoe, Piyavochnoe, Barskoe and Alekseevskoe) and six eutrophic fens (Shichengskoe, Piyavochnoe, Rodionskoe, Ileksa, Povreka and Charozerskoe) located in the Vologda region of European North Russia, within the zone of middle taiga. Detailed descriptions of the sampling sites and the sampling procedure are reported elsewhere [17,19]; characteristics of the examined peatlands are also provided in Table 1.

Peat Sampling for Incubation Studies
To examine the response of peat-inhabiting representatives of the Rokubacteriales to methane availability, one additional batch of peat samples was collected in June 2021 from the eutrophic fen site of the mire Shichengskoe (59 • 56 56 N, 41 • 16 59 E). The samples were transported to the laboratory, homogenized and used for the determination of methane oxidation activity, incubation experiments and molecular analyses.

Methane Incubation Experiments
Weighted portions of peat (10 g, water content 91%) sampled from the fen Shichengskoe were placed in 160 mL glass flasks, which were then sealed hermetically. Methane (CH 4 ) was injected in the flasks up to the concentration of~1000 ppm and the flasks were incubated at room temperature for 24 h. Samples of the gas phase (0.5 mL) were taken from the flasks periodically and analyzed for CH 4 concentration on a Kristall 5000 chromatograph (Khromatek, Yoshkar-Ola, Russia). The experiments were made in triplicate. The rates of methane oxidation by the samples were calculated in mg CH 4 g −1 of wet peat h −1 . After these activity measurements, methane was added to the flasks up to the concentration of 10% (v/v) in the headspace. The flasks were then incubated at room temperature for 4 weeks. At the end of each week, the flasks were flushed with ambient air using a sterile filter (0.22 µm) to remove remaining CH 4 and accumulated CO 2 , and methane was re-injected up to the concentration of 10% (v/v), in order to keep high methane availability during the whole incubation period. Portions of peat from each incubation flask were taken before and after incubation with CH 4 and used for the molecular analysis.

Bioinformatic Analyses of Microbial Community Structure
Both the datasets of 16S rRNA gene sequences from Dedysh et al., 2020 and from the current study were analyzed with QIIME 2 v.2019.10 (https://qiime2.org) [22]. The DADA2 plugin was used for sequence quality control, denoising and chimera filtering [23]. Operational taxonomic units (OTUs) were clustered applying the VSEARCH plugin [24] with open-reference function using the Silva v. 132 database [25,26] with 97% identity. Taxonomy assignment was performed using BLAST against the Silva v. 132 database with 80% identity.

Correlations between Microbial Groups and Chemical Properties of Peat
Graph Pad Prism v. 6.0.1 (GraphPad Software, San Diego, CA, USA, www.graphpad. com) package was used to calculate Pearson correlation coefficients between abundances of the taxonomic groups and chemical properties of peat samples. Tests were considered significant if they had a p-value < 0.05.

Analyses of the Publicly Available Rokubacterial Metagenomes
All available rokubacterial MAGs were retrieved from the GenBank using taxid and combined into the single database. The entire pool of OTUs identified in this study was compared against this Rokubacteria database using BLAST (v. 2.9.0) [27] (accessed on 15 November 2021) with a sequence identity threshold of 97%. The constructed database was annotated using Prokka [28]. Amino acid sequence of the putative monooxygenase from the assembly Amazon-R-15-13 was taken as a query for tblastn searches against the Rokubacteria database with a sequence identity threshold of 30%. All valid blast hits were included into the phylogenetic analysis. A genome-based tree for members of the Rokubacteria was reconstructed using the Genome Taxonomy Database and GTDBtoolkit (https://github.com/Ecogenomics/GtdbTk (accessed on 16 November 2021)) [29], release 04-RS89. The maximum likelihood phylogenetic tree was constructed using MegaX software [30]. Functional annotation of CDSs in MAGs was performed on the basis of the results of BLASTP searches [27] against the NCBI (National Center for Biotechnology Information) nonredundant database and the KEGG (Kyoto Encyclopedia of Genes and Genomes) database [31].

Microbial Community Composition in Boreal Wetlands
The pool of 16S rRNA gene reads examined in this study included 1,024,783 sequences (average length,~440 bp), which were retrieved from the peat samples and retained after quality filtration, denoising and removing chimeras. The microbial community composition in raised bogs was clearly different from that in fens ( Figure 1). Thus, archaeal populations in raised bogs were represented by members of the Euryarchaeota and Taumarchaeota, while fens were colonized by the Euryarchaeota, Crenarchaeota and Nanoarchaeaeota. The patterns of bacterial diversity in raised bogs and fens were also different. Raised bogs were dominated by the Acidobacteria (from 30.0 ± 0.4% to 42.4 ± 2.9% of total reads retrieved from the Alekseevskoe and Barskoe fens, respectively, mean ± SE). The next numerically abundant bacterial groups in bogs were the Planctomycetes (11.0 ± 2.7% and 19.2 ± 0.4% in the Piyavochnoe and Barskoe), Proteobacteria (10.7 ± 2.0% and 15.7 ± 2.8% in Alekseevskoe and Piyavochnoe) and Verrucomicrobia (11.3 ± 1.2% and 14.8 ± 0.5% in Barskoe and Piyavochnoe). By contrast, the bacterial communities in fens were dominated by the Proteobacteria (from 21.1 ± 3.2% to 30.6 ± 2.1% of total reads retrieved from the Piyavochnoe and Rodionskoe fens, respectively) and Chloroflexi (from 15.7 ± 1.8% to 30.2 ± 3.7% of total reads retrieved from the Shichengskoe and Ileksa fens, respectively). Other major groups were the Acidobacteria (7.2 ± 0.1% to 22.3 ± 1.2% retrieved from the Ileksa and Shichengskoe) and Patescibacteria (6.0 ± 0.5% to 18.8 ± 2.3% in Povreka and Piyavochnoe). Notably, members of the Methylomirabilota were detected in all studied fens but were absent from raised bogs (Figure 1). These bacteria were especially abundant in Charozerskoe and Shichengskoe fens (4.2 ± 0.3% and 3.3 ± 0.3% of total reads) but represented only a minor group in fens Povreka and Ileksa (0.3 ± 0.1% of total reads in both wetlands). The relative abundances of Rokubacteria determined in our study are in agreement with previous reports. Bacteria of this group were low in abundance at almost every site where they were identified [11,32], and often were represented by a single 16S rRNA gene sequence. The two so far reported exceptions were the microbial community in a grass root zone in the Angelo Coast Range Reserve, California, where Rokubacteria constituted~10% of all prokaryotes [33] and the microbial community in virgin soils of southern Vietnam tropical forests, where the relative abundance of these bacteria reached 6% in the lower soil horizons [14].

Diversity and Most Abundant OTUs of Methylomirabilota
Since Rokubacteria were only detected in eutrophic fens, our further analyses were focused on the wetlands Shichengskoe, Piyavochnoe, Rodionskoe, Ileksa, Povreka and Charozerskoe. The pool of Methylomirabilota-affiliated 16S rRNA gene fragments obtained in our study included 7837 reads. The highest numbers of reads were received from Shichengskoe and Charozerskoe fens (3458 and 2647, respectively) and the lowest from Ileksa and Povreka (247 and 244). Taxonomic analysis revealed that peat samples from Ileksa and Povreka fens were mostly represented by uncultured bacteria from the family Methylomirabilaceae of the order Methylomirabilales (98% and 100% of all Methylomirabilotaaffiliated reads). On the contrary, samples from Shichengskoe and Charozerskoe fens were dominated by uncultivated members of the order Rokubacteriales (96% and 93%, respectively). In peat samples from the fens Piyavochnoe and Rodionskoe, the proportions of Rokubacteriales-/Methylomirabilales-like reads were as follows: 59/41% and 55/45%, respectively (upper panel in Figure 1).
In total, 137 OTUs of Methylomirabilota determined at 97% sequence similarity were identified in our study. Of these, 1 OTU was present in all fens, 3 OTUs were detected in four samples, 3 OTUs were detected in three samples, 8 OTUs were present in two samples and the remaining 122 OTUs were specific for each fen (Figure 2). The list of the most abundant OTUs comprising ≥ 0.1% of all reads retrieved from the corresponding fen is given in Table 2 and is shown in the phylogenetic tree in Figure 3. The most abundant OTUs were represented by uncultivated members of the Rokubacteriales (OTUs No 1, 3). However, OTU 8, which was present in all studied fens, belonged to the family Methylomirabilaceae. Overall, among the fourteen most represented OTUs, three OTUs belonged to the family Methylomirabilaceae and eleven OTUs affiliated with the order Rokubacteriales (Table 2, Figure 3).  Table 2. The OTUs affiliated with Methylomirabilaceae and Rokubacteriales are indicated by red and blue, respectively.  All rokubacterial MAGs available in the GenBank were combined into the single database, which also contained MAGs of higher quality from GTDB. OTUs retrieved from the examined fens were compared against this database using BLAST with a species-level sequence identity threshold of 97%. Nine blast hits against GTDB metagenomes affiliated with MAGs assembled from grassland soil and aquifer well ecosystems (Figure 4). The remaining blast hits were represented by the GenBank MAGs from meadow soil, wetland sediment, freshwater lake and hot spring sediment (not shown).

Correlations between Chemical Properties of Peat and the Abundance of Methylomirabilota
To reveal the influence of several chemical characteristics of peat on the relative abundances of Methylomirabilota groups, Pearson correlation analysis with further significance test was performed. Relative abundance of Rokubacteriales correlated positively with pH, total nitrogen (TN), Ca and Mg availability, and negatively correlated with total organic carbon (TOC), sulfate, Fe and P availability ( Figure 5). Notably, a strong positive correlation with pH, TN, Mg and Ca availability, as well as a negative correlation with TOC, was observed for the representatives of Methylomirabilaceae ( Figure 5).

Microbial Community Response to Incubation with Methane
To examine the response of peat-inhabiting Rokubacteriales to methane availability, peat samples from the fen Shichengskoe were incubated with 10% (v/v) CH 4 for 4 weeks. Incubation with methane resulted in a pronounced increase of methane-oxidizing activity of peat samples from 0.20 ± 0.02 to 0.40 ± 0.09 mg CH 4 g −1 h −1 . The bacterial community composition in these samples was analyzed before and after incubation with methane. The original methanotroph community was composed of representatives of the order Methylococcales, particularly uncultivated Crenothrix bacteria, and members of the family Beijerinckiaceae ( Figure 6). The strongest positive response to methane availability was detected for members of the order Methylococcales, which were represented mostly by Methylovulum spp., Methylomonas spp., Methyloglobulus spp. and Methylomagnum spp. (Figure 6). In contrast, the relative abundances of the Beijerinckiaceae and Crenothrix declined after incubation with methane, which was also true for members of the Methylomirabilaceae. No statistically significant response of the Rokubacteriales to CH 4 availability was observed in this incubation experiment.

Search for the Presence of Genomic Determinants of Methanotrophy in Rokubacterial MAGs
To verify the hypothesis regarding the occurrence of methanotrophic capabilities in Rokubacteria, we inspected all metagenomes of these microorganisms available in the Gen-Bank for the presence of genes encoding pMMO, a key enzyme of aerobic methanotrophy. The latter is encoded by the pmoCAB gene cluster; the pmoA is used as a functional gene marker for detecting aerobic methanotrophs in the environment [34,35]. The pMMOs belong to the superfamily of Cu-containing membrane-bound monooxygenases (CuMMOs), which also includes ammonia monooxygenases (AMO) and a number of short-chain alkane and alkene monooxygenases [36,37]. The genes coding for CuMMOs are usually organized in the same way as pMMO (xmoCAB gene cluster) [37]. It should be taken into account, however, that the search for the presence of xmoCAB operons in rokubacterial MAGs has some limitations due to the lack of complete rokubacterial genomes.
The amino acid sequence (WP_008359136.1) deposited in GenBank as a putative MMO from the MAG Amazon-R-15-13 [3] was used as a query against the database of rokubacterial metagenomes. Only one blast hit displayed a reliable identity threshold (89%). This amino acid sequence belonged to the rokubacterial MAG AR15 retrieved from a meadow soil (GCA_003220675) (Figure 7). Both XmoA sequences were organized in clusters with XmoB and XmoC. The evolutionary relationships of rokubacterial XmoA sequences among other members of the CuMMO superfamily were determined by means of phylogenetic analysis (Figure 7). Notably, XmoA sequences from rokubacterial MAGs did not cluster with PmoA sequences of true methanotrophs from the Alphaproteobacteria, Gammaproteobacteria, Verrucomicrobia and candidate phylum NC10 [8,38,39]. Instead, two rokubacterial XmoA sequences clustered together with XmoA sequences from several members of the Actinobacteria, which are known for the ability to utilize short alkanes but not methane. Three of these actinobacteria, namely Nocardioides sp. strain CF8 and Mycolicibacteria rhodesiae strains NBB3 and NBB4, possess putative CuMMOs (pBMOs), which function as butane monooxygenases [40,41]. The fourth actinobacterium in this clade, Rhodococcus sp. ZPP, is capable of growth on propane, which is dependent on the presence of CuMMO genes [42]. Thus, clustering of both rokubacterial XmoA sequences with those from Actinobacteria suggests that Rokubacteria may potentially act as short chain alkane oxidizers. Although CuMMOs display clear specialization for one particular substrate, they can co-oxidize several substrates, including methane and alkanes containing up to five carbons [43]. We, therefore, looked also for the occurrence of C1-metabolism genes in two selected rokubacterial MAGs, AR15 and Amazon-R-15-13. Rokubacteria bacterium AR15 possesses PQQ-dependent methanol dehydrogenase (MDH), which displays 57% identity to that in Methylococcus capsulatus Bath. MDH is most likely absent from Amazon-R-15-13, since the best blast hit revealed only 27% identity to MDH sequence in AR15. The genes coding for formate dehydrogenase are present only in AR15. Formaldehyde assimilation may potentially function via the serine cycle. Completeness of the latter is difficult to confirm due to the incompleteness of both MAGs. However, phosphoenolpyruvate carboxylase, a functionally important enzyme, is absent from both AR15 and Amazon-R-15-13.
On the other hand, utilization of short chain hydrocarbons is supported by the presence of genes encoding alcohol and aldehyde dehydrogenases as well as genes coding for the pathway of β-oxidation of fatty acids. The fatty acids produced by aldehyde oxidation are further metabolized by β-oxidation, generating Acyl-CoA, which enters the tricarboxylic acid cycle [44]. Therefore, AR15 and Amazon-R-15-13 most likely represent bacteria that utilize short chain hydrocarbons.

Conclusions
In summary, our study identified neutral eutrophic fens as one of the common habitats for bacteria of the order Rokubacteriales. Relative abundance of these microorganisms correlated positively with pH, total nitrogen, Ca and Mg availability. No response to methane availability was detected for the Rokubacteriales, while a clear increase in relative abundance was observed for the conventional Methylococcales methanotrophs. The search for methane monooxygenase encoding genes in the publicly available Rokubacteriales metagenomes yielded negative results. Some members of this bacterial order may possess Cu-containing membrane-bound monooxygenases, but the exact substrate specialization of these enzymes in Rokubacteriales remains unclear and requires experimental verification. Apparently, peat-inhabiting Rokubacteriales are non-methanotrophic bacteria that colonize nitrogen-rich wetlands. Funding: This study was supported by the Russian Science Foundation (project 21-14-00034) and the Ministry of Science and Higher Education of the Russian Federation.

Data Availability Statement:
The raw data generated from sequencing of the 16S rRNA gene fragments before and after methane incubation studies have been deposited in the NCBI Sequence Read Archive (SRA) under the accession numbers SAMN23339836 and SAMN23339837, respectively (BioProject PRJNA782125).

Conflicts of Interest:
The authors declare no conflict of interest.