The Cullars Rotation is a century old crop rotation experiment alternating between five crops during the course of a year (cotton, crimson clover, corn, wheat and soybean) (Mitchell et al., 2011). The rotation is composed of 14 plots for the study of the effects of controlled deficiencies of phosphorus (P), potassium (K), nitrogen (N) and sulfur (S), as well as other soil amendments such as lime (for pH control) and addition of a micronutrient mix containing boron (B), zinc (Zn), manganese (Mn), copper (Cu) and iron (Fe) on crop yield (Mitchell et al., 2005) (Figure 1).

Soil samples were obtained from the Cullars Rotation in March 2019. Three topsoil samples from each one of the 14 plots in Cullars Rotation were collected for soil DNA extraction, alongside ~400 g of bulk soil for analysis. The samples for DNA extraction were immediately transported to the lab and kept refrigerated at 4°C until processing. Processing of the 42 soil samples (3 replicates/14 plots) were completed on the same day as sampling. The bulk soil samples were submitted for analysis at Auburn University’s Soil Testing Laboratory the same day. Soil nutrients were obtained using Mehlich I Extraction and analyzed by ICP. Measurements were obtained for calcium (Ca), K, magnesium (Mg), P, aluminum (Al), B, Cu, Fe, Mn, sodium (Na), Zn, Nitrate-Nitrogen (NO₃-N), all in parts per million (ppm). Percent carbon ©, N, and pH were also measured.

For each sample, 5 g of soil was suspended in 25 mL of 0.2% tetrasodium pyrophosphate buffer (pH 8.0), which has been shown to dissociate microbial cells from soil particulates (e.g., clay, sand and silt) while retaining the microbial community structure associated with bulk soil (Portillo et al., 2013; Yuan et al., 2018). Suspensions were homogenized for 10 min, with gentle inversions, followed by 10 cycles of intermittent vortexing for 5 s. The soil slurry was subjected to a low-speed centrifugation (5,000 x g for 3 min) to precipitate larger soil particulates and eukaryotic cells. The supernatant was collected and serially filtered through a 5.0 μm cellulose membrane, then 1.2 μm and 0.45 μm SUPOR® membranes (Pall Corporation, Port Washington, New York) in tandem. The <0.45 μm fraction was then pelleted by centrifugation (20,000 x g for 20 min) and DNA was isolated from the cell pellet using the Power Soil DNA isolation kit (MoBIO Labs) to obtain the <0.45 μm fraction metagenomic DNA (Figure 1). There were three topsoil replicates from each of the 14 plots in the Cullars Rotation, resulting in 42 different samples from which a < 0.45 μm fraction metagenome was generated.

Each <0.45 μm fraction metagenome was uniquely indexed and paired-end sequenced using an Illumina HiSeq platform in one dedicated run (Figure 1). Raw reads from all metagenomes were processed using Trimmomatic v0.39 (Bolger et al., 2014). The remaining reads from all 42 samples were combined and assembled with MEGAHIT v1.2.9 (Li et al., 2015). The reads from each sample were mapped back against the assembled contigs using BWA v0.7.17 (Li and Durbin, 2009). Parallel binning was conducted with MetaBAT2 v2.15 (Kang et al., 2019), MaxBin 2.0 v2.2.7 (Wu et al., 2016) and CONCOCT v1.1.0 (Alneberg et al., 2014). Resulting genomic bins were curated using DasTool v1.1.4 (Sieber et al., 2018). The completeness and contamination of each bin was assessed using CheckM v1.1.3 (Parks et al., 2015) using the standard bacterial marker set and a CPR-specific phylogenetic marker set (Brown et al., 2015). Bins with completeness below 65%, contamination above 10% or strain heterogeneity above 10% were dropped from subsequent analyses. Bins with completeness above 85%, contamination below 5% and strain heterogeneity below 5% were considered high-quality MAGs, while bins with completeness above 75%, contamination below 5% and strain heterogeneity below 5% were considered good-quality MAGs and the remainder were considered draft MAGs. To aid in genome assembly and binning, an additional set of reads from a similar sequencing run of a pilot enrichment for Cullars Rotation’s Plot C was added to the read pool for a second round of assembly and binning following the same steps and criteria above. The total bins obtained were dereplicated using Derep v0.12.8 (Olm et al., 2017) and used in the subsequent analysis.

MAGs were taxonomically classified using two strategies: 16S rRNA gene sequence homology and phylogenetic marker gene voting. First, 16S rRNA sequences were extracted using CheckM v1.1.3 (Parks et al., 2015) and a custom package for CPR 16S rRNA extraction and intron removal (Brown et al., 2015). The MAG-associated 16S rRNA genes were classified against the SILVA database (Quast et al., 2012) (accessed February, 2020) using SINA (Pruesse et al., 2012) through the web client. A homology threshold of 75% was used for phylum classification and 80% for class classification. For the second approach, 13 phylogenetically informative marker genes (Supplementary Table S1) corresponding to ribosomal proteins with similar phylogenetic signal to 16S rRNA (Gregor et al., 2016), were extracted from each bin using CheckM v1.1.3. The markers were annotated against the nr/nt NCBI database. Then each marker “casts a vote” according to their taxonomic annotation. The taxonomy with the higher number of “votes” was assigned to the MAG. Finally, the consensus classification from the two strategies was applied to the MAG. Both strategies used NCBI taxonomy, which is itself based on the taxonomy and nomenclature proposed by Brown et al. (2015).

MAGs classified as CPR were selected for a robust additional classification by constructing a concatenated maximum likelihood tree. A curated collection of reference CPR genomes was constructed by downloading 1,200 CPR genomes from NCBI (accessed October 2019) and filtering out low quality MAGs with the program CheckM v1.1.3 by the following requirements: >75% completeness, < 5% contamination and strain heterogeneity. The set of 13 phylogenetic markers was extracted from each reference genome. For each marker, the Cullars MAG sequences and the references sequences were aligned using MAFFT v7.475 (Katoh and Standley, 2013), trimmed with TrimAl v1.4.1 (Capella-Gutierrez et al., 2009). Then, the marker alignments were concatenated using Phylemon 2.0 (Sanchez et al., 2011) and manually inspected with Geneious 10.3.¹ The concatenated alignment of 13 markers was used to generate a maximum likelihood tree with 1,000 bootstraps using RAxML v8.2.9 (Stamatakis, 2014). Marker genes from the Dependentiae MAGs were used to root the tree. The final tree was visualized and annotated using iTOL (Letunic and Bork, 2019), using Brown et al. (2015) proposed nomenclature. In an effort to also comply with the efforts toward an standardized bacterial classification and taxonomy (Parks et al., 2020), our MAGs were classified according to GTDB (Parks et al., 2022) using GTDB-Tk v2.2 (Chaumeil et al., 2022), and utilized the (Brown et al., 2015) proposed nomenclature for phylogenetic affiliations.

To assess metabolic and nutrient utilization potential for each of our CPR MAGs, each of the sufficient quality MAGs and members of the same taxa from the curated collection of reference genomes, were annotated using DRAM v1.3.5 (Shaffer et al., 2020). This software utilizes a series of databases (UniRef90, PFAM, dbCAN, Refseq viral, VOGDB, and MEROPS peptidase) to annotate genes and reconstruct the organism’s metabolic profile. DRAM annotations were analyzed with special attention to carbohydrate utilization via carbohydrate-active enzymes (CAZymes) (Lombard et al., 2014). Profiles between MAGs and their corresponding references were compared using a 2-way analysis of variance (ANOVA) based on a linear regression model using the lme4 (Bates et al., 2015) and emmeans (Lenth et al., 2019) packages. Output was analyzed in RStudio v. 4.0.3.

The <0.45 μm fraction bacterial community from Cullars soil was reconstructed from the assembled MAGs. To estimate MAG counts, reads from each sample were individually mapped back against each MAG using Bowtie2 (Langmead and Salzberg, 2012) with high sensitivity. Well mapped reads were counted using SAMTools v1.3 (Li et al., 2009). The matrix with read counts per sample was imported into Rstudio. Library sizes were normalized using the cumulative scale sum strategy (CSS) with metagenomeSeq (Paulson et al., 2013). Counts were then normalized by genome length by dividing the total read length by the MAG genome size (bp). Coverages below 1% were discarded as inconclusive and assigned to zero (Levy-Booth et al., 2021). Fractions were rounded up to generate a list of normalized natural counts. These counts were used to reconstruct MAG community and calculate average MAG relative abundance per plot.

Rarefaction curves for each sample were calculated using the package Ranacapa (Kandlikar et al., 2018). For both MAG– and 16S rRNA gene-based community assessments, non-metric multidimensional scaling (NMDS) and correlation analysis (CA) were performed with Bray-Curtis calculated distances using the package PhyloSeq (McMurdie and Holmes, 2013). Differences in community composition across plots were accessed by performing PERMANOVA with 999 permutations and subsequent multivariate homogeneity of groups dispersions (Betadisper) using the package Vegan (Oksanen et al., 2013). The same package was used to perform vector fitting and surface fitting of Cullars Rotation soil parameters with significant correlation to community changes (p < 0.05).

Total DNA was extracted from bulk soil samples which were the same used for CPR cell enrichment. DNA from plots 03, 08, 10 and 11 were used for 16S rRNA gene PCR amplification. The V4 region of the 16S rRNA gene was amplified and paired-end sequenced. Paired-end reads were processed using FLASH (Magoc and Salzberg, 2011) to produce raw tags, which were quality filtered following the QIIME2 protocol (Caporaso et al., 2010). The clean high-quality tags were compared against the SILVA database and chimeras were removed using UCHIME (Edgar et al., 2011). Sequences with >97% similarity were clustered into the same operational taxonomic unit (OTU) using Uparse. Taxonomic annotation of OTU representative sequences was carried out against the SILVA SSU database using Mothur (Schloss et al., 2009).

OTU absolute abundance was normalized by rarefaction, and the subsequent counts were used for alpha and beta diversity analyses. The community comparisons between plots and ordinations were constructed using the same programs outlines above for the MAG communities. Linear discriminant analysis Effect Size (LEfSE) analysis was conducted using the LEfSe software (Segata et al., 2011).

Normalized counts from the 10 most abundant phyla in the16S and MAG abundance data were combined to access the possible relationships between the relative abundance of CPR MAGs and 16S ribotypes from soil bacteria. The resulting count table was used to build the network using the package SPIEC-EASI (Kurtz et al., 2015). The graphs were drawn using the package iGraph (Csardi and Nepusz, 2006) and exported into Cytoscape (Smoot et al., 2011) for visualization.

A previous study of a Cullars Rotation soil metagenome revealed that CPR bacteria were present in this soil in extremely low abundance, but were accessible using a direct cloning and sequencing approach that avoided PCR and 16S rRNA gene-associated biases (Santana-Pereira et al., 2020). Therefore, for each of three samples from the 14 plots at Cullars Rotation, ultra-small cells <0.45 μm were enriched from soil samples to directly sequence and obtain CPR genome bins. From the 124 non-replicated bins, a total of 42 genomic bins were curated ranging from 400 Kb to 3.5 Mb (average 1.17 Mb). We identified 27 MAGs with acceptable quality, which was defined as completeness below 65%, contamination over 10% or strain heterogeneity over 10% (Supplementary Figure S1). These MAGs had an average completeness of 86.8%, with 5 bins later identified as CPR with 100% completeness, and an average genome size of 1.13 Mb in accordance with the reduced genome size characteristic of CPR bacteria. A noticeable exception was MAG DT 16 with a genome size of 3.55 Mb and completeness of 78.6%; however, this relatively larger genome was affiliated with the phylum Gemmatimonadetes which is not known to have a small size but has been previously observed to be capable of transmission through a small sized filter (Portillo et al., 2013).

CPR classification has experienced several revisions (Rinke et al., 2013; Brown et al., 2015; Parks et al., 2017, 2018), highlighting the challenges associated with conducting phylogenetic analyses of uncultured microorganisms (Konstantinidis et al., 2017). Initial studies on CPR taxonomy, based on up ta few thousands of genomes, categorized it as a monophyletic group encompassing several individual candidate phyla, such as Saccharibacteria, and the superphyla Parcubacteria and Microgenomates, each containing numerous candidate phyla (Brown et al., 2015; Anantharaman et al., 2016; Hug et al., 2016). This classification is most congruent to NCBI taxonomy and thus popularized. With the increasing availability of MAGs and genomes, much more comprehensive studies using from 8 to more than 90 thousand genomes made widespread rank normalization and standardization possible, settling on Patescibacteria/CPR as a single phylum (Parks et al., 2017, 2018). This framework was used to develop the Genome Taxonomy DataBase (GTDB) (Chaumeil et al., 2019; Parks et al., 2020), which proposed new nomenclatures and taxonomic ranks for CPR. Thus the previous candidate phyla Saccharibacteria, Parcubacteria and ABY1, correspond to the classes Saccharimonadia, Paceibacteria and ABY1, respectively (Supplementary Table S2). Due to their widespread use and congruence with popular NCBI nomenclature, the names CPR, Saccharibacteria and Parcubacteria will be favored throughout the paper.

A collection of 13 phylogenetically informative marker genes (Supplementary Table S1) and 16S rRNA genes were extracted from Cullars MAGs and were independently compared to the NCBI GenBank nr/nt and SILVA databases, respectively, to form a consensus homology classification of the MAGs. All curated MAGs were classified, with one representative MAG recovered from each of the phyla Actinobacteria, Chlamydiae, Dependentiae and Gemmatimonadetes, and the remaining belonging to the CPR group. Within CPR, there were representatives of the Parcubacteria (16), Saccharibacteria (10) and ABY1 (1) groups (Figure 2). Additionally, CPR MAGs were classified using the GTDB classifier, showing high congruence with the previous classification, since Parcubacteria, Saccharibacteria and ABY1 MAGs, were classified as Paceibacteria, Saccharimonadia and ABY1, respectively. Moreover, GTDB provided finer classification up to family level. Interestingly, a number of MAGs associated with underrepresented taxa in GTDB were recovered from Cullars Soil, such as representatives of families VEQN01 and CAYRIW01, which have only 3 and 1 reference genomes, respectively.

Given the high degree of diversity among members of the CPR group, the candidate CPR MAGs were further classified by constructing a concatenated maximum likelihood tree using 13 phylogenetic markers and 920 reference genomes, allowing the selection of closest relatives for metabolic potential comparisons (Figure 2). The tree confirmed the previous taxonomical assignments. Cullars CPR MAGs affiliated with the Parcubacteria/Paceibacteria class, were all from the order UBA9983_A and distributed across 10 families that correlate to the broad groups Nomurabacteria, Taylorbacteria, Adlerbacteria and Kaiserbacteria in NCBI taxonomy. Additionally, two MAGs (DT 10 and DT 38) clustered with an unclassified CPR MAG to form a deeper branch within Saccharibacteria (Figure 2). These MAGs were classified by GTDB as belonging to the order UBA4664, which carries only 9 references, corresponding to 1.25% of all Saccharibacteria genomes in the database. Indeed, this low representation, whether due to sampling or environmental rarity, might explain the observed deeper branching, although it might be mitigated with the insertion of more reference genomes. The ABY1 group MAG was further classified into the Uhrbacteria group.

The MAGs recovered from Cullars Rotation soil represent three distinct CPR classes including the Saccharibacteria (10), Parcubacteria (16), and ABY1 (1). These MAG genomes (Figure 3) were compared to a curated collection of 124 genome references from NCBI (Saccharibacteria, 43; Parcubacteria, 49; ABY1, 32) to assess their metabolic potential utilizing the metabolic profiling tool DRAM (Supplementary Figure S2).

ABY1, Saccharibacteria, and Parcubacteria MAGs carried partially complete glycolysis, pentose phosphate and reductive pentose phosphate pathways, as consistently observed among these classes of CPR. Additionally, they were found to have nearly complete F-type ATPase. All MAGs were lacking a complete tricarboxylic acid cycle (TCA), subunits of nicotinamide adenine dinucleotide (NADH) dehydrogenase and electron transport chain (ETC), which indicates an anaerobic and fermentative metabolism, although the presence of complete or nearly complete F-type ATPases was widespread. Interestingly, Parcubacteria, namely from the family UBA2103, and Saccharibacteria from Cullars were more likely to carry mostly complete operons coding for cytochrome o ubiquitol oxidase. Cytochrome o functions by catalyzing the oxidation of ubiquinol (CoQH₂) into ubiquinone (CoQ), coupled with the reduction of O₂ into water, supplying energy to pump H+ into the outer membrane contributing to the formation of a proton motive force (Aussel et al., 2014). Interestingly, it was more common among the soil CPR than the references.

The CAZyme database and its enzyme classes were utilized to understand Cullars MAGs carbohydrate metabolic potential, namely glycoside hydrolases (GH, carbohydrate degradation), glycotransferases (GT, carbohydrate biosynthesis), carbohydrate esterases (CE, hydrolysis of carbohydrate esters), carbohydrate-binding modules (CBM, adhesion of carbohydrates), and auxiliary activities (AA, redox coenzymes for lignocellulose degradation) (Supplementary Figure S4). Most Cullars MAGs contained genes for enzymes associated with starch, cellulose/hemicellulose, and chitin degradation with the top hits among the Saccharibacteria and Parcubacteria, but not present in the Cullars ABY1 MAG (Figure 4). Interestingly, only Cullars Saccharibacteria carried the potential for lignin degradation. This continues to support soil CPR as containing the necessary enzymes for degradation of complex plant carbohydrates and soil necromass including bacterial and fungal cell walls.

ABY1 and Saccharibacteria carbon substrate metabolic profiles were not found to be statistically different from that of the reference genomes available on NCBI (p > 0.05). However, Cullars Parcubacteria MAGs did demonstrate a strongly statistically difference from their references (p < 0.001). Parcubacteria was the most represented CPR in our MAGs and particular samples (i.e., DT.28) contained several hits for particular enzyme classes for peptidoglycan and chitin degradation, suggesting these soil microorganisms may be particularly adept for scavenging dead or decaying soil microflora.

Cullars CPR MAGs did not demonstrate any ability to participate in nitrogen or methane cycling. Select Saccharibacteria were predicted to have some capacity for dissimilatory and assimilatory sulfate reduction. Cullars MAGs were predicted to degrade organic forms of nitrogen, particularly, the amino acid methionine in the cysteine biosynthesis pathway. Likewise, Cullars CPR had many predicted peptidases, indicating the ability to degrade peptides from the soil environment.

Sampling the Cullars Rotation soil fertility plots offers a unique opportunity to assess the influence of nutrient availability and soil parameters on the soil CPR community, due to its century-old controlled nutrient deficiency experiments (Mitchell et al., 2011) (Figure 1). To reconstruct the MAG community in each plot, reads from each sample were mapped back against the curated MAGs and normalized to account for variations in library size and genome length (Figure 1). The normalized, filtered counts were used to reconstruct the communities in each of the plots in Cullars Rotation (Figure 5A). Additionally, soil nutrient content for each plot was also measured creating a profile that could be correlated with changes in community structure (Figure 5B).

MAG communities are composed predominantly of Saccharibacteria, followed by Gemmatimonadetes and Dependentiae, even though the latter two have only one representative MAG each. Despite having the higher and most diverse number of representatives, Parcubacteria accounted for a significantly smaller proportion of the community, indicating that taxa dominance was observed in the sampled metagenomes. Interestingly, there was an observed spike in Actinobacteria abundance in Plot 08, which received no pH correction treatment and had a significantly lower pH than other treatment groups.

MAG abundance within the different CPR taxa varied among the different experimental treatment groups (Supplementary Figure S5). In particular, Plot 08 (no pH correction) showed the greatest shift in community structure, with 78% of the relative abundance corresponding to one Actinobacteria and two Saccharibacteria bins, which were dominant exclusively in this plot, while also having the highest number of recovered Parcubacteria genomes. Dependentiae and Gemmatimonadetes genomes were also recovered across all plots, with the later peaking in abundance under increased phosphorous amendment (Plot 05) and dropping to marginal abundance when there was no soil pH correction (Plot 08) indicating a susceptibility to very acidic soils. Otherwise, the abundance of CPR genomes was not significantly different among the different Cullars Rotation plots.

NMDS ordination was used to visualize communities, dissimilarities in relation to fertilization treatment. To further explore the role of soil parameters in community compositions, the measurements of soil macro and micro-nutrients were modeled to the community distance matrix (Figure 5B). Community profiles did vary significantly across plots (PERMANOVA, p = 0.001 and Betadisper, p = 0.8537) with the soil plots with an acidic pH having the most dramatic effect in clustering (p = 0.015). The lack of pH correction coupled with the addition of fertilization treatment resulted in the most acidic pH across all plots for Plot 08 (pH = 4.08 ± 0.08), surpassing even Plot C which received no amendment nor fertilization (pH = 4.58 ± 0.04); however, pH alone could not account for the observed differences (PERMANOVA, p = 0.199), with similar trends not observed for Plot C, which also has a very acidic soil. The environmental fit analysis reveals that the factors with the strongest correlation to the changes in the community were aluminum (Al) and iron (Fe) concentrations (p = 0.001 and p = 0.0001, respectively).

Next, the relationship between the MAG community and overall bacterial community was explored. 16S rRNA gene amplicon data were used to reconstruct the overall bacterial community, recovering representatives of 40 phyla. The overall community was dominated by Proteobacteria, Actinobacteria, Acidobacteria, Bacteroidota, and Firmicutes (Figure 5C). Plot 08 had lower alpha diversity (Supplementary Figure S6) and an overrepresentation of the phylum Chloroflexi. Communities between plots were significantly different (p = 0.001), with Plot 08 also clustering away from other plots in PCA ordinations (Figure 5D). Similarly to the MAG community shifts, an acidic pH was correlated with significant 16S community differences (p = 0.004); however, the only soil parameter with a significant correlation to the community was Al concentration (p = 0.001).

Linear discriminant analysis effect size (LEfSe) was used to identify taxa highly correlated to community shifts in each plot that can act as biomarkers. 18 OTUs were marked as Plot 08 biomarkers at varying taxonomical levels, with the strongest signal belonging to the phylum Chloroflexi, class Ktedonobacteria (Supplementary Figure S7).

The potential relationship between Cullars’ MAGs and members of the soil microbiome was explored by constructing networks between Cullars’ MAGS and soil phyla (Figure 6). The MAGs were correlated with the occurrence of diverse phyla, albeit the weight of the edges was low. No particularly strong relationship was identified between the MAGs or any particular taxa in the overall soil microbiome, even when focusing on OTUs and MAGs distinctively associated with Plot 08, such as MAGs DT35 and DT24, or Plot 08-specific biomarkers such as Chloroflexi or Acidobacteria taxa.