Occurrence, Diversity, and Genomes of “Candidatus Patescibacteria” along the Early Diagenesis of Marine Sediments

We investigated CPR in marine sediment cores from three locations: Mohns Ridge (MR) (which is part of the Arctic Mid-Ocean Ridge) core MR-GS14-GC08, North Pond (NP) core NP-U1383E, and organic-rich sediments off the Costa Rica Margin (CR) at site U1379 (core CR- U1379B) (Fig. 1A), where concomitant metagenome and 16S rRNA gene amplicon sequencing data were published (45–47) and are available for detailed analyses of CPR. The seafloor depths of these sites are in the range of 127 to 4,435 m. Sediments at these sites are characterized by a nearly constant rate of sediment deposition but without underlying energy inputs like hydrothermal vents or episodic sediment changes like mass wasting events in trenches. Based on the existing geochemical profiles (see Fig. S1 in the supplemental material) and definitions described previously (45–47), the investigated sediments represent a variety of redox zones ranging from the oxic zone down to the sulfate reduction zone (Fig. 1), although they do not represent a coherent time series of sediment samples formed at a single location. The four sediment horizons from the MR core span the upper oxic zone (10 cm), the oxic-anoxic transition zone (OATZ; 100 cm), the nitrate-ammonium transition zone (NATZ; 160 cm), and the Mn reduction zone (250 cm) (Fig. 1B). Additionally, a total of four oligotrophic sediment horizons from the NP core were included, two of which are located in the upper oxic zone (100 cm and 1,000 cm), and the remaining two represent the transition zones between the oxic and anoxic sediments (i.e., the oxic-anoxic transition zone [at ~2,200 cm] and the anoxic-oxic transition zone [AOTZ, at 2,950 cm below seafloor {bsf}]) (Fig. 1E). Finally, for the CR core (47), four sediment horizons in the interval of 2 to 9 m below the seafloor were included (Fig. 1H). Sulfate should be the dominant electron acceptor in sediments of this interval, due to the combination of the following three observations. (i) Sulfate in the porewater decreased with depth (Fig. 1H), although its depletion depth (~40 m below seafloor [mbsf]) was not covered by the investigated sediment layers of this study. (ii) The sulfur isotope of sulfate (δ³⁴SO₄) increased with depth (48), indicating biological sulfate reduction. (iii) Dissolved Mn in the porewater, the product of Mn oxide reduction, showed a decreasing rather than increasing trend with depth in the investigated sediment interval (Fig. S1C), indicating that Mn oxides were no longer an important terminal electron acceptor. Therefore, this extensive data set allowed us to capture CPR diversity and genomes from a wide range of redox conditions formed in the sediment’s early diagenesis.

We first examined the presence and diversity of CPR in the metagenome sequencing data. By searching and classifying the putative 16S rRNA gene reads using phyloFlash (49), we found that CPR bacteria account for 2.0 to 4.2% of the total microbial communities in MR core GC08 (Fig. 1B), 0.4 to 28% in NP (Fig. 1E), and 0.3 to 1.6% in CR (Fig. 1H). We noted that the variations in the relative abundance of CPR in different sediment layers, especially between those layers in the NP core, were not artifacts caused by different metagenome sequencing depths between samples, because the total 16S rRNA gene reads detected in the metagenome sequencing data sets were comparable (Table S1). In contrast, based on the concomitant 16S rRNA gene amplicon sequencing data, there were only 26 (1.1% of the total recovered operational taxonomic units [OTUs]), 27 (0.4% of the total), and 3 (0.3% of the total) OTUs affiliated with CPR at NP, MR, and CP, respectively. These CPR OTUs accounted for <2% of the total microbial communities of the three sediment sites (Fig. 1B, E, and H), lower than those estimated by metagenome sequencing. The discrepancy between these two methods may be related to the fact that the amplicon sequencing underestimates the relative abundance of CPR because most of the divergent and intron-containing 16S rRNA genes of CPR can evade the detection of universal primers (10). Nevertheless, these results suggested that members of CPR are usually part of the rare community (<5% of the total prokaryotic communities, irrespective of the detection method) yet are widespread (frequently encountered) in marine sediments of various redox zones, which was also supported by reports from marine sediments of other locations such as hadal trench sediments in the Pacific Ocean (41, 42).

Metagenome sequencing also revealed a higher diversity of CPR than did the 16S rRNA gene amplicon sequencing. Based on the taxonomical classification of the 16S rRNA gene reads in the metagenome sequencing data, CPR of the classes “Candidatus Paceibacteria,” “Ca. Gracilibacteria,” “Ca. Microgenomatia,” “Ca. Saccharimonadia,” “Ca. ABY1,” and “Ca. WWE3” were present in the three sediment sites. Among the total 12 sediment layers where metagenome sequencing data are available, “Ca. Paceibacteria” was generally the most abundant CPR class, although the occasional dominances of “Ca. Gracilibacteria” at 22 mbsf of NP (Fig. 1G) and “Ca. Microgenomatia” and “Ca. ABY1” in the three deep layers of CR were also detected (Fig. 1I). In contrast, only members of the classes “Ca. Paceibacteria” and “Ca. Gracilibacteria” (see Fig. S2 for the phylogenetic placements of individual OTUs) were detected by 16S rRNA gene amplicon sequencing (Fig. 1C and F): “Ca. Paceibacteria” was the dominant CPR class (~80% of the total CPR community) in most of the sediment layers of MR core GC08 (Fig. 1C), and members of “Ca. Gracilibacteria” dominated in most of the examined sediment layers of NP (Fig. 1F). The discrepancy in obtained CPR community structure from the two methods again demonstrated that the amplicon sequencing method can probably capture only a subset of the total CPR population in complex communities.

To further characterize CPR in these sediments, we recovered CPR genomes from the existing metagenome sequencing data of these sites. We obtained four CPR MAGs (MR_Bin143, MR_Bin147, MR_Bin1662, and MR_Bin1762) from Mohns Ridge, two (NP_Bin194 and NP_Bin050) from North Pond, and five from sediments of the Costa Rica site (CR_Bin034, CR_Bin039, CR_Bin047, CR_Bin053, and CR_Bin021) (Table 1). These 11 CPR MAGs are considerably small in genome size (0.39 to 1.10 Mbp) and contain 12 to 153 scaffolds (Table 1). These genomes have 462 to 1,264 coding sequences, with the coding density varying in the range of 86.2 to 92.0% (Table 1). Based on the universal 43 single-copy genes of CPR bacteria (10), all except two (MR_Bin1762 and CR_Bin021) were estimated to be of >97.7% completeness with <2.3% redundancy (Table 1). CPR bacteria are known to lack some of the single-copy genes found in non-CPR bacteria (10), so based on the single-copy genes automatically determined for bacteria by CheckM (50), completion estimates were lower (61.0 to 80.3% complete) (Table 1). All of these CPR MAGs except CR_Bin039 have a reconstructed 16S rRNA gene and thus can be regarded as high-quality MAGs.

The novelty of these CPR MAGs was evident from phylogenetic analyses based on (i) the 14 concatenated ribosomal proteins (Fig. 2A) and (ii) the 16S rRNA gene (Fig. 2B). To ensure consistency with the current literature, we adopted the updated genomic distance-based classification scheme of GTDB (51) to name CPR lineages. While six CPR classes were present in the three sediment sites, as revealed by the bulk metagenome sequencing data (Fig. 1), all 11 recovered CPR MAGs were classified as members of novel genera, families, or even orders within the class “Ca. Paceibacteria.” The phylogenetic novelty of these genomes was also supported by the GTDB classification, in which all the CPR MAGs showed <80% average nucleotide identities (ANI) with their most similar genomes included in the GTDB (07-RS207) database, and relative evolutionary divergence (RED) values of 0.74 to 0.83 were calculated. In particular, NP_Bin050 represented a new order (named o__JAHCSD01 in GTDB 07-RS207), with the closest phylogenetic relationship to members of the order Portnoybacterales (Fig. 2A). MR_Bin147 was a member of a new genus in the family GWB1-50-10 in the order UBA6257, an order mainly constituted by MAGs that previously were classified as Jorgensenbacteria recovered from the groundwater environment (20, 52) (Fig. 2A). The other CPR MAGs, belonging to the order “Ca. Paceibacterales,” showed phylogenetic novelties at suborder levels. MR_Bin1662 and MR_Bin1762 formed a new family (we provisionally named it “Candidatus Bathypaceibacteraceae”) (Fig. 2A). Three MAGs from the CR core (CR_Bin039, CR_Bin053, and CR_Bin034) also formed a new family (provisionally named “Candidatus Sedimentipaceibacteraceae”) (Fig. 2A). NP_Bin194 represented a new genus within the RBG-13-36-15 family. In addition, MR_Bin143 represented a new genus in the UBA10102 family, which is composed mainly of genomes previously classified as Wildermuthbacteria. Finally, CR_Bin047 represented a new genus within the family GWA2-38-27 (Fig. 2A).

The phylogenetic novelty of these CPR MAGs was supported by the phylogenetic analysis of the 16S rRNA gene, which was broadly congruent with the concatenated ribosomal protein tree (Fig. 2B). It also showed that close relatives of these CPR MAGs were previously detected by clone libraries in anoxic marine sediments of various locations, including the Gulf of Mexico (53), Shimokita subseafloor sediments (54), Angola Basin sediments (55), and South Pacific Gyre ferromanganese nodules (56), although quantitative information about their abundance at these sites is lacking. The phylogenetic relatedness between CPR MAGs in marine sediments of different oceanographic regions may reflect the habitat or ecological niche preferences or the availability of the hosts of the CPR groups.

After integration of the 16S rRNA gene amplicon sequences into the 16S rRNA gene phylogenetic tree, it was clear that three of the four CPR MAGs recovered from MR were also captured by the 16S rRNA gene amplicon sequencing: MR_Bin1662 showed a 100% match with OTU_152, MR_Bin1762 matched OTU_556, and MR_Bin147 corresponded to OTU_525. It was also clear that the recovered CPR genomes represent only a subset of the CPR diversity revealed by the 16S rRNA gene amplicon sequencing. Despite that most CPR MAGs from NP and CR contained a 16S rRNA gene sequence, none of them was a good match with the amplicon-recovered OTUs from these sediments. The misrepresentation of these novel CPR MAGs in the 16S rRNA gene amplicon sequencing data may be due to the divergence of their 16S rRNA genes because they have one or more mismatches with the forward PCR primer used than those recovered from the MR sediments (see Fig. S3 in the supplemental material). Another possibility is that these CPR MAGs were rare taxa such that stochasticity and local heterogeneity might have led to their evasion in the PCR amplification. “Ca. ABY1” and “Ca. Microgenomatia,” two important CPR classes that prevail in the three investigated sediment cores, especially in the CR core (Fig. 1I), were absent in our genome inventory but had been recovered from Mariana Trench hadal sediments (40). This indicated that the CPR MAGs detected by bulk metagenome sequencing are indeed part of the microbiome of the vast marine sedimentary environment. The mismatch of 16S rRNA gene sequences between the MAGs and the amplicon sequencing also suggested that some yet unknown CPR bacteria remained to be revealed in global marine sediments (for an example, see reference 57) and highlights the necessity of genome recovery for the discovery of novel CPR bacteria.

The presence of the 11 CPR MAGs was confined to anoxic sediment layers at their source locations. Based on the genome coverage calculation, CPR bacteria represented by these MAGs prefer anoxic sediment layers. The four MAGs from MR were present mainly in anoxic sediments in the NATZ or the Mn reduction zone below the oxygen depletion depth (Fig. 3A), the two CPR MAGs from NP mainly existed in the AOTZ (29.50 mbsf), and the five MAGs from CR were detected in the uppermost sequenced sediment horizon (2.0 mbsf) and to a lesser extent in the deepest horizon (~9.0 mbsf) (Fig. 3A). Read mapping across the three sediment sites suggested that the individual CPR MAGs were generally specific to the sediment core of origin (Fig. 3A), except that the five CPR MAGs recovered from CR were also detectable in the AOTZ of core NP-U1383E (Fig. 3A). Similar site-specific diversity of CPR genomes was also recently reported in the groundwater and freshwater lake environment by large-scale metagenomic surveys (23, 58).

As with where they were detected, the growth of the 11 CPR MAGs, inferred from the index of replication (iRep; values of >1 indicate proliferation) (59), occurred only in their primary niches (i.e., the sediment layer where the highest genome coverage was detected for a certain genome). The highest iRep of the four individual MAGs from MR occurred in either the NATZ or the Mn reduction zone, the two from NP in the AOTZ, and the five from CR in the uppermost anoxic zone (Fig. 3B). The calculated iRep (in the range of 1.2 to 2.2) suggested that these CPR MAGs were actively replicating in their primary niches at the time of sampling.

CPR are present in a wide range of redox niches but may not be directly selected by prevailing electron acceptors. Through amplicon sequencing and genome binning, we detected CPR in a variety of redox niches in marine sediments, with the dominant electron acceptors shifting from oxygen to nitrate, Mn oxides, and sulfate. Although none of the CPR MAGs was recovered from oxic sediment metagenomes, CPR 16S rRNA genes were detected in the investigated oxic sediments, especially in cores MR-GS14-GC08 and NP-1383E, which contain extensive oxic zones (Fig. 1). We estimated the absolute abundances of CPR in these two cores as the product of the total cell numbers (i.e., the sum of archaeal and bacterial 16S rRNA gene abundances) (46) and the relative abundances of CPR shown in Fig. 1. In cores of both MR and NP, the estimated absolute abundance of CPR can reach over 10⁵ cells g⁻¹ in the oxic zones (Fig. 4). Among the CPR genomes detected in MR sediments, MR_Bin147 was present in the oxic zone of GC08 (Fig. 3), despite that this genome lacks genes encoding the cytochrome o ubiquinol oxidase (complex IV) involved in oxygen reduction (see results below). Most known CPR organisms were previously detected in oxygen-limited or anoxic environments (33, 60), whereas some CPR organisms were also occasionally reported in oxic groundwater (24, 60–62), freshwater lakes (58), and soils (63). Our survey suggests that the overall CPR population was present in marine sediments throughout the early diagenesis processes even in bulk oxygenated sediments, but individual genomes showed preferences for different anoxic layers in different cores. Whether they live in the presence of oxygen or inhabit microenvironments in sediment particles without oxygen remains unknown. Considering that (i) the recovered CPR genomes seem not to have metabolic pathways sensitive to the external redox condition changes (the switching of the dominant terminal electron acceptors) and (ii) they have their narrow niches (indicated by the restricted distributions of their presence and active replication [i.e., iRep values]), the distribution of CPR bacteria in marine sediments may not be directly selected by the prevailing electron acceptors. The lack of direct environmental dependences of CPR has also been observed in a recent large-scale genome survey in the groundwater environment (24).

Limited energy metabolisms and biosynthesis capacities in marine sediment CPR. The CPR organisms recovered from marine sediments are likely fermentative organotrophs with simplified energy metabolism pathways. All recovered MAGs lack hydrogenases, especially the cytoplasmic bidirectional group 3 [NiFe] hydrogenase that has been proposed to be involved in pumping protons to build the proton motive force and help ATP generation in CPR genomes (33). All marine sediment CPR MAGs, except MR_Bin147, lack ATPases (Fig. 5), indicating that they are not capable of conserving energy by proton motive force. Instead, they may synthesize ATP through partial glycolysis via fermentation and substrate-level phosphorylation (33), similar to the recently characterized CPR bacterium Vampirococcus lugosii (64). This was confirmed by the presence of the glycolysis pathway-related genes in all sediment CPR MAGs recovered in this study, which facilitates the degradation of glucose to produce pyruvate or further to acetate, although some genes that regulate few intermediate steps were missing likely due to the incomplete nature of the CPR genomes (Fig. 5). The important intermediate compound of glycolysis, fructose 6-phosphate, can be provided by the pentose phosphate pathway, which is encoded in most of the recovered CPR MAGs (Fig. 5). We also note the lack of a complete tricarboxylic acid cycle (TCA), despite some CPR MAGs having a small subset of the enzymes in this cycle, likely for biosynthetic purposes. Like their close relatives from other environments (33), the marine sediment CPR MAGs lack a respiratory electron transport chain, as evidenced by the absence of NADH dehydrogenase (complex I) and complexes II to IV of the oxidative phosphorylation pathway (Fig. 5), suggesting that they are nonrespiring.

Although various glycoside hydrolases (GHs) have been detected in some CPR MAGs in groundwater (22, 65) and their expression in some CPR has been detected in subseafloor sediments (66), the GHs detected in our CPR MAGs are very limited. Six MAGs have a GH1 (for hydrolyzing carbohydrate moiety), while GH3, GH57 (starch), GH63, and GH130 (mannose) were also detected in fewer than three CPR MAGs (Fig. 5), suggesting that marine CPR MAGs have very limited saccharolytic capacities. However, the majority of the detected carbohydrate-active enzymes (CAZymes) belong to families 2 and 4, which associated with glycolipid synthesis, similar to CPR genomes detected in groundwater (21). Marine sediment CPR MAGs also lack the genes encoding nitrate reductase that were previously reported in genomes recovered from hadal sediments (43). Although the copper-containing nitrite reductase (NirK) has been previously noticed in some CPR genomes (22, 33), none of the marine sediment CPR genomes contain this gene. Therefore, whether CPR MAGs have an impact on the nitrogen cycle in marine sediments remains unclear.

Most of the sediment CPR genomes lack complete biosynthesis pathways for amino acids, lipids, cofactors, and nucleotides. Among the 20 known amino acids, only the genes for the biosynthesis of lysine and arginine are present (Fig. 5), indicating that they may need to obtain the rest of the amino acids from the external environment or their presumed host cells. Similarly, for nucleotide synthesis, they have only 4 to 14 genes for the metabolisms of purine and pyrimidine (Fig. 5). The CPR genomes have no recognized genes responsible for the synthesis of lipids or cofactors. Also lacking in these sediment CPR genomes are genes for flagellar biosynthesis (only 1 of the 46 required genes is encoded) (Fig. 4), indicating that they are probably nonmotile in the sediment environment. These CPR genomes also lack ABC type transporters and have very limited genes (<5 genes) for the Sec secretion system (Fig. 5). However, like other CPR genomes (33), all sediment CPR MAGs have extensive genes for the biosynthesis of peptidoglycan, suggesting that they may have intact cell walls.

Given the multiple auxotrophies detected in the recovered genomes and small genome sizes, we anticipate an episymbiotic lifestyle for the sediment CPR MAGs similar to that of their nonmarine relatives (34, 64, 67, 68). The lack of ABC transporters in the CPR MAGs may force them to rely on the host cells to obtain the necessary substrates for their metabolism. In MR_Bin147, the genome carries genes for hemolysin synthesis and hemolysin transporter, which could export hemolysin to the surface of the host cells and contribute to the host cell wall and membrane disruption and cell content release (Fig. 6). Similar to CPR genomes in other environments (27, 34, 69, 70), marine sediment CPR genomes harbor genes for type IV pilus synthesis (Fig. 6), which may provide access for the membrane-bound translocation complex to environmental double-stranded DNA (71). In addition, it also has a competence-related integral membrane protein, ComEC, which plays a role in the uptake of host DNA (72) that can be degraded to various restriction endonucleases to provide the nucleotides necessary for growth (Fig. 6).

We attempted to explore the potential prokaryotic hosts of CPR in sediments of the MR core, where CPR organisms were the most abundant among the three sites, by performing a cooccurrence network analysis based on the existing 16S rRNA gene amplicon sequencing data. Our results indicated that MR_Bin147 (represented by OTU_525) (Fig. S4) has a tight association with OTU_59 only, a member of the order Gemmatimonadales in the phylum Gemmatimonadota. Therefore, this analysis suggested that Gemmatimonadales bacteria could be hosts of CPR in MR sediments. It is worth noting that the (relative) abundances of OTU_59 (1.20% to 2.50% of the total communities) were almost 4-fold higher than those of of OTU_525 (01.4% to 0.74%), suggesting that only a fraction of Gemmatimonadota cells may harbor episymbionts in the MR sediments. This is consistent with the observation that “Candidatus Nanosynbacter lyticus TM7x” (Saccharibacterium) cells establish a long-term parasitic association with host cells (Actinomyces odontolyticus XH001) by infecting only a subset of the population (73). Future cultivation efforts are necessary to confirm such prediction and to also explore the hosts of the other CPR bacteria identified in these marine sediments.

Although the 11 CPR MAGs are phylogenetically novel and were confined within their ideal redox niches, they are not geographically exclusive to the three sediment cores investigated in this study. Recently, Zhou et al. (40) reported >500 microbial MAGs from the Challenger Deep sediments of the Mariana Trench, which included 24 “Ca. Patescibacteria” MAGs. Based on the GTDB-tk classification, “Ca. Paceibacteria” was the most abundant class (n = 14), while members of the classes “Ca. ABY1” (n = 4), “Ca. Andersenbacteria” (n = 1), and “Ca. Microgenomatia” (n = 5) were also detected. By adding the “Ca. Paceibacteria” MAGs into our phylogenetic trees (Fig. 2A and B), we found that eight of these “Ca. Paceibacteria” MAGs were affiliated with the novel branches together with our MAGs. In particular, six fell into the family UBA10102 (containing MR_Bin143), one into RBG-13-36-15 (containing NP_Bin194), and the final into the order UBA6257 (containing MR_Bin174) (Fig. 2A). This comparison indicated that the CPR MAGs of these branches may have wider distribution ranges.

We report the abundance, diversity, and metabolic capacities of CPR bacteria in three sediment sites with contrasting environmental contexts, in which the examined sediments span most of the redox zones developed during the early diagenesis of organic matter. In this vast environment, CPR bacteria are generally rare (accounting for <5% in all but one sample) in the total microbial communities and are mainly in anoxic sediments, although ~10⁵ cells g⁻¹ of CPR cells are also present in the oxic zone. 16S rRNA gene amplicon sequencing tends to severely underestimate the diversity and abundance of CPR bacteria in the examined marine sediments, probably due to primer biases. The recovered CPR MAGs are novel at the family or genus level in the “Ca. Paceibacteria” class. All recovered CPR MAGs should be fermentative organotrophs relying on the truncated glycolysis for energy generation but lack complete biosynthesis pathways for amino acids, lipids, cofactors, and nucleotides and therefore are either symbionts or closely dependent on other community members for these key building blocks. The novel CPR genomes recovered in this study lay a foundation for further phylogenomic analysis of CPR in marine sediments. Despite some recent progress from untargeted metagenome sequencing surveys (for examples, see references 39 and 40), the total number of CPR MAGs from marine sediments is still much smaller than that from other natural habitats such as the groundwater (10, 23, 24), freshwater lakes (58), and soil (63). More targeted approaches (with some preselection before sequencing efforts) should be employed in future investigations to fully capture the overall diversity of CPR and other ultrasmall-celled microbes in global marine sediments. Cultivating diverse prokaryotes from native sediments and using cultured hosts to cue the presumed episymbionts in the laboratory (for examples, see references 27 and 44) could be among the promising and scalable ways to coculture CPR and to shed more light on these microbes in the vast habitat of marine sediments. Our study indicates that in global marine sediments, there are still some microbes that are yet uncaptured by single-gene-based surveys.

We analyzed MAGs affiliated with Candidate Phyla Radiation (CPR) in marine sediments recovered from three different locations, the Arctic Mid-Ocean Ridge (MR-GS14-GC08) and North Pond (NP-U1383E) and Pacific coastal sediments (core CR-U1379B) off the Coast Rica Margin (47). The details about sample collection, DNA extraction, 16S rRNA gene amplicon sequencing and analysis, metagenome sequencing, assembly, binning, and refinement were described in previous reports (46, 47). We highlighted the contrasting redox conditions between the three sediment sites based on previously reported geochemical profiles (45–47). For those measurements, sediment porewater samples at discrete depths were extracted using either Rhizon samplers (for Arctic and North Pond cores [74]) or a titanium squeezer (for the Costa Rica core [75]) without air contacts to minimize alterations caused by sampling, handling, and measurement processes. Briefly, oxygen concentrations in sediment cores from the Arctic and North Pond cores were measured directly when the cores were still in the core liner by using a needle-type fiber optic oxygen microsensor (PreSens, Regensburg, Germany). Other porewater constituents, including nitrate, nitrite, ammonium, and dissolved inorganic carbon, were measured by a QuAAtro continuous flow analyzer (SEAL Analytical Ltd., Southampton, UK), using colorimetrical methods specific for each solute (46). Sulfate concentrations in the porewater of CR-U1379B were measured using a Dionex ICS-3000 ion chromatograph.

We leveraged the previously generated shotgun metagenome sequencing data and 16S rRNA gene amplicon sequencing data of these three sites to examine the diversity and distribution of CPR bacteria. For the metagenome sequencing data, 16S rRNA gene reads in the quality-controlled sequences were identified and taxonomically classified using phyloFlash v3.2 beta1 (49), which uses BBMap (76) to identify and then classify putative 16S rRNA gene reads using the SILVA 138 release (77) as the reference database. The fractions of CPR classes in the total prokaryotic communities in the sediment samples were calculated and are shown as pie charts.

All the amplicon sequencing data sets were generated using “universal primers” (Uni519F/806R) (45) that target the variable region 4 (V4) of archaeal and bacterial 16S rRNA genes. Among the nine variable regions, the V4 region contains the second lowest frequency of insertions after the V8 region (10) and therefore was selected as the targeting region in our primer-dependent approach to capture CPR in complex communities, although primer-free methods (for an example, see reference 78) in theory should be better alternatives to characterize novel microbial populations. OTUs (220 bp, clustered at the cutoff of 97% nucleotide similarity) of putative CPR MAGs classified by CREST (79) against the SILVA 138.1 release were extracted from the OTU table. OTUs with <10 total reads across all sediment samples were excluded. Sediment samples with <5 CPR reads (i.e., 0.025% of the total community in GC08 and 0.05% of the total community in NP-U1383E) were excluded and are indicated as “not detected.” The relative abundance of each CPR OTU in each sample was calculated using normalization scaling. OTUs were aggregated at the CPR class level and visualized using bar charts made using the R package ggplot2 (80).

CPR MAGs with identifiable 16S rRNA gene fragments were matched with the short amplicon OTUs by blast alignment using the MAG 16S rRNA gene as the query and the OTU sequences from the corresponding sediment cores as the reference database. For example, if the 16S rRNA gene of MAG_A showed a >99% match with OTU_001 in core NP-U1383E, they were considered matches and the distribution of MAG_A in NP-U1383E was considered to be that of OTU_001. This method is helpful to describe the distribution of MAGs in all the sampled sediment horizons rather than just those of the metagenome sequenced.

The eight metagenome sequencing data sets from MR (n = 4) and NP (n = 4) were coassembled, while the four from CR were also coassembled using MEGAHIT v1.2.9 (81), using k-mers from 21 to 141 (–k-min 21 –k-max 141 –k-step 10 –presets meta-large). Contigs longer than 1,000 bp were binned using MaxBin2 v2.2.5 (82) and MetaBAT v2.15.3 (83), and the highest-quality ones were selected using DAS_Tool (84) with the default parameters. The resulting MAGs were quality assessed using CheckM v1.1.3 (50) and taxonomically classified using GTDB-Tk v2.0.0 (parameters: classify_wf defaults) (85). As a proxy of relative abundance, the coverage of each genome in each sample was determined by read mapping using BBMap v38.76 (76), with a 99% identity threshold (minid = 99, idfilter = 99). To refine and improve the quality of each MAG, reads in the highest-coverage samples were recruited onto the scaffolds of the MAG and reassembled using SPAdes v.3.12.0 (86), with the k-mers of 21, 33, 55, and 77. The resulting assembly was refined using gbtools (87), mainly by plotting coverages of contigs in the samples in which the two highest coverages were determined and outlier contigs were removed. We used CheckM v1.1.3 (50) to assess the quality of the CPR MAGs based on two sets of genes: (i) the single-copy genes that are automatically determined by CheckM for bacteria and (ii) the 43 single-copy genes proposed by Brown et al. (10) that were implemented in CheckM as a specific custom-defined workflow (for detailed steps, see https://github.com/Ecogenomics/CheckM/wiki/Workflows#using-cpr-marker-set).

Genes in the CPR genomes were predicted using Prodigal v2.6.3 (88). Genome annotation was conducted using Prokka v.1.13 (89), eggNOG v2.1.4 (90), and BlastKoala (91) using the KEGG database. The functional assignments of genes of interest were also confirmed using BLASTp (E value < 1e−10, sequence identity > 50%, and shortest alignment rate > 60%) against the NCBI RefSeq database. The metabolic pathways were reconstructed using KEGG Mapper (92). Carbohydrate-active enzymes (CAZYmes) were searched using dbCAN hmm profiles (downloaded April 2019 from http://bcb.unl.edu/dbCAN2/download/Databases/dbCAN-old@UGA/), using an E value of 1e−18 and a coverage of >0.35, per the instructions of the developers.

Phylogenomic analysis was performed separately for the 11 MAGs of CPR recovered. In addition, the 14 MAGs of the “Ca. Paceibacteria” class from the Mariana Trench sediments reported previously (40) were also included. This phylogenomic analysis was based on marker genes consisting of 14 concatenated ribosomal protein genes (rpL2, -3, -4, -5, -6, -14, -16, -18, and -22 and rpS3, -8, -10, -17, and -19). These genes in the genomes were identified in Anvi’o v6.2 (93) by hidden Markov model (HMM) profiles. Sequences of four “Candidatus Doudnabacteria” genomes were used as the outgroup. Sequences of each marker gene were aligned individually using MUSCLE (94), and alignment gaps were removed using trimAl (95), with the mode “automated.” Individually trimmed alignments were concatenated. The maximum-likelihood phylogenetic tree was reconstructed using IQ-TREE v1.6.10 (96) with LG+F+R6 as the best-fit evolutionary model selected by ModelFinder (97) and 1,000 ultrafast bootstraps by UFBoot2 (98).

Likewise, phylogenetic analyses of the 16S rRNA gene were also performed for the CPR MAGs. 16S rRNA gene sequences in the MAGs were annotated using barrnap (99) implemented in the Prokka package. In some cases where more than one copy of 16S rRNA gene fragments was found in a genome, these fragments were assembled by aligning them with other high-quality intron-free CPR 16S rRNA gene sequences in Unipro UGENE (100) and removing the 100% overlap region. These sequences of MAGs were used to identify their close relatives by a BLASTn search in the NCBI database with a similarity threshold of 90%. All sequences were aligned using MAFFT-LINSi (101). Given the known large insertions/introns in CPR 16S rRNA gene sequences (10), the alignment was manually inspected using Unipro UGENE (100) and gaps were manually removed. A maximum-likelihood phylogenetic tree was constructed using IQ-TREE v1.6.10 (96) with SYM+R6 as the best-fit evolution model selected by ModelFinder (97), and 100 ultrafast bootstraps were performed by UFBoot2 (98). Relevant amplicon OTU sequences (220 bp) were added to the phylogenetic tree without changing the overall topology to identify the clades of the OTU sequences.

Index of replication (iRep) values of the recovered CPR MAGs in all the metagenomic sequenced depths were calculated in accordance with the procedure outlined in reference 59. Before the calculations, contigs of <5 kb in each genome bin were removed. The quality (completeness and contamination levels) of the bins assessed by CheckM (50) was not affected by this omission and fulfilled the prerequisite of genome quality of such calculation.

Network analysis was performed for the amplicon sequencing data sets of the MR core GC08. Given that the CPR MAGs were recovered mainly from the anoxic sediments of GC08 (Fig. 3A), the sediment horizons collected from the oxic zone were removed to minimize the effect of niche filtration on network analysis outcomes. We also retained only those OTUs representing >100 reads totally across all samples in this analysis to reduce the effect of spurious correlations caused by low-abundance taxon distribution. The relative abundance data were converted to ratio data using the centered log ratio (Clr) method. The network was inferred using SPIEC-EASI (102) using the “glasso” method and characterized using the igraph package (103). The nodes in the reconstructed network represent OTUs at 97% identity, and their sizes are proportional to their degrees (number of connections), while the edges (i.e., connections) correspond to strong and significant (positive or negative) correlations between nodes. Only positive correlations were plotted, whereas four weak negative correlations existed in the data set but were not shown.

All sequencing data generated in this study have been deposited in the NCBI Sequence Read Archive under project numbers PRJNA529480 (MR), PRJNA489438 (NP), and PRJNA599172 (CR). Raw metagenomic sequencing data are available in NCBI under BioSample numbers SAMN11268098 (GC08_10cm), SAMN11268104 (GC08_100cm), SAMN11268106 (GC08_160cm), and SAMN11268109 (GC08_250cm). The MAGs described in this study are available from NCBI under the accession numbers provided in Table 1.