We tested for differences in community characterization based on 16S RNA vs. DNA biomarker sequence patterns in the rhizosphere of B. stricta, a perennial herb native to Wyoming. Seeds for this study were originally collected from the Snowy Range Mountains (41.32971759902109 N, -106.50515422710646 W), and grown for one generation in the greenhouse to increase seed numbers and minimize maternal effects. Prior to planting, seeds were surface sterilized, by rinsing for 1 min in 70% ethanol 0.1% Triton 30% RO water mixture, then rinsed in RO water, then rinsing for 12 min in a 10% bleach 0.1% Triton 90% water mixture. Seeds were then rinsed a final three times with RO water, before being placed on sterile filter paper for ease of planting (adapted from Lundberg et al., 2012).

Surface-sterilized seeds were planted into pots with a mixture of field and potting soils. For the field soil, we collected soil from unvegetated sites adjacent to a field location with a native B. stricta population, referred to hereafter as the Crow Creek field site (CRW). Field soil was sieved to 4 mm to remove large debris, then autoclaved three times for 30 min with the soil being mixed between autoclaving steps. This soil was then mixed with autoclaved potting soil [Redi-Earth Potting Mix (Sungro Horticulture, Agawam, MA, United States)] in a 9:1 ratio; we included a small percentage of potting mix because its greater water-holding capacity relative to the field soil improves the overall rate and synchrony of seed germination. This soil mixture was next inoculated with a 4% v/v of non-autoclaved field soil inocula. This approach of autoclaving all soil and then applying a field soil inoculation was used to ensure that detected microbes were those native to B. stricta and not derived from the potting mix.

Before field transplanting, seeds germinated in 2″ mesh net pots (2″ Inch TEKU Net Slit Pots for Hydroponic Aeroponic Use) and were allowed to grow for 4 weeks under greenhouse conditions (UW Laramie Research and Extension Center, Laramie, WY) with ambient day/night light and temperature cycles. In addition to pots planted with seeds, which would be used to characterize the rhizosphere microbiome, we also prepared soil-filled pots without plants, which would be used to estimate the microbiome of bulk soil or unvegetated microsites, referred to hereafter as the “control soil.” Pots were placed in a randomized checkerboard array into tray blocks, such that no pot was directly adjacent to another. Plastic covers were placed over all trays to retain humidity and promote germination. All materials for planting, such as bench tops, trays, pots, and covers were bleached and rinsed before use. Pots were individually watered, initially via subirrigation and after 2 weeks of growth via overhead misting. Covers were removed 2 weeks after germination. Three weeks after germination, all pots were acclimated to the outside environment via 2-h field exposures. Plant germination time, rosette size, and true leaf number were measured weekly to estimate plant performance.

In June of 2019, 4 weeks after germination, plants were transferred to the Crow Creek field site, which has a naturally occurring population of B. stricta. The Crow Creek (CRW) field site was located in the Medicine Bow-Routt National Forest in south-eastern Wyoming (41.227318 N, -105.383343 W), with an elevation of ~2,560 m. Six 26 cm by 140 cm plots were cleared of plants and debris. Mesh pots were randomized into each plot and planted roughly 10 cm apart in two rows of 24 pots. Mesh pots were removed from filters and cups and placed directly in field site plots, to minimize root damage while transplanting and to facilitate collection of plant rhizospheres. Rosette size and true leaf number were measured weekly to estimate plant performance. All pots were watered every other day using RO water, and checked for insect damage.

Four weeks after being transplanted to the field site, soil and plant samples were harvested. Rhizosphere soil, bulk soil, and control soil samples were collected over the course of 3 days, from July 15 2019 to July 17 2019 between 1:30 p.m. and 2:30 p.m. At each collection time point, pots with and without plants were randomly harvested and bulk soil was collected directly from within each plot. Samples from the three treatments types were handled as follows: The large samples of soil collected in the field from (1) the control pots without plants and (2) the bulk soil were stored in whirlpacks (Whirl-Pak^® Bags); smaller subsamples of these soils were taken and stored in 2 mL test tube for later nucleic acid extraction. (3) For collection of rhizosphere soil from the pots with plants, mesh pots were removed from the soil, with care taken to ensure that roots which had grown out of the mesh pots were damaged minimally. Plants were then removed from the pots and shaken to remove excess loose soil. Soil that adhered closely at approximately 1 mm to the root mass was considered to be rhizosphere soil. Roots and adhered soil were separated from plant leaves and stem using flame sterilized scissors, and placed in falcon tubes containing approximately 200 mL of PBS buffer (200 ul silwet, 900 mL RO water, 100 mL 10 × PBS) (adapted from Bulgarelli et al., 2013). Samples were stored on ice in the field. Immediately upon returning to the lab, bulk and control soil samples were transferred to a-80°F freezer. For subsequent sample processing, falcon tubes containing rhizosphere soil and plant root tissue were defrosted, then lightly vortexed to remove adhered soil from root tissue. Root tissue was removed using sterile forceps. The soil slurry was then vacuum filtrated through a 0.2 nm filter. The resulting soil and filter were transferred to a 0.5 mL centrifuge tube, and flash frozen with liquid nitrogen, before being stored long term at −80°C.

During the extraction process, we took several steps to minimize biases that might artificially inflate the differences between RNA- and DNA-derived community composition measurements. Methodological biases were minimized through the use of technical replicates, and the simultaneous extraction of nucleic acids from a single soil sample, so each paired DNA and RNA community profile derives from an identical collection (Moeseneder et al., 2005; Gentile et al., 2006; Morgan et al., 2010; McCarthy et al., 2015). Some methodological differences were unavoidable, such as RNA but not DNA having a reverse transcription step (Zhen et al., 2015). Control, bulk, and rhizosphere soil RNA and DNA samples were simultaneously extracted using the methods described in the RNeasy PowerSoil Total RNA Kit and RNeasy PowerSoil DNA Elution Kit (Qiagen, 2017). After extraction, RNA was reverse transcribed to cDNA using the QuantiTect Reverse Transcription Kit (Qiagen, 2009). Negative control blank samples were included for extractions and reverse transcription. Samples were stored at −20°C until further processing.

Positive control (a ZymoBiomics mock community) and negative control (blank-H₂O) samples were included in library preparation. 16S rRNA-gene DNA and 16S rRNA-transcript cDNA amplicons were amplified using the 515–806 (Walters et al., 2016) primer pair to amplify the V4 region of the 16S rRNA locus. Two technical replicates were made for each sample; non-blank samples were normalized to a standard concentration of 10 ng/ul. Kapa HiFi Hot Start polymerase, Kapa HiFi Hot Start buffer and reagents, and HPLC grade water were used during PCR. PCR conditions for the first round were: 95° for 3 min; followed by 15 cycles of 98° for 30 s, 62° for 30 s, and 72° for 30 s; with a final 72° elongation step for 5 min and a 4° hold. PCR products were cleaned using AxyPrep MagBead magnetic beads (Axygen; Union City, CA, United States). PCR conditions for the second round were: 95° for 3 min; followed by 19 cycles of 98° for 30 s, 55° for 30 s, and 72° for 30 s; with a final 72° elongation step for 5 min and a 4° hold. Products from the second round of PCR were also cleaned using AxyPrep MagBead magnetic beads. Library amplification success was confirmed using a Bioanalyzer fragment analyzer (Agilent; Santa Clara, CA, United States). PCR amplicon libraries were sequenced by Psomagen (Rockville, Maryland, United States) on an Illumina NovaSeq 6000 using 2 × 250 paired-end sequencing.

A custom perl script (created by C. Alex Buerkle) was used to demultiplex sequence data, and unique reads were dereplicated using vsearch v.2.9.0 (Edgar, 2010). Identified reads were clustered using the “cluster_unoise” (Quast et al., 2015) algorithm and a 99% similarity threshold, and sequences that occurred 12 or more times were considered a potential OTU. Chimeric sequences were removed using “uchime3_denovo” algorithm (Edgar, 2010), and the resulting OTU were used to make an OTU table using the “usearch_global” algorithm. Taxonomy was assigned to each OTU using vsearch v 2.15 and the Silvia v123 (Quast et al., 2015) reference database, with a minimum bootstrap confidence of 80%. Computing was performed using the Teton Computing Environment at the Advanced Research Computing Center, University of Wyoming, Laramie.¹

Reads in negative control samples that occurred in low frequency (i.e., less than 100 reads per sample) in other samples were considered contamination and removed from the data set. Reads that did not assign to the kingdom bacteria were removed. Samples with less than 30,000 total reads were removed from analysis, and samples were filtered so each RNA derived sample had a corresponding DNA sample that was extracted from the same collected soil. OTU that occurred in less than 10% of all samples with less than 2 reads were removed. To minimize uninformative noise, OTU with less than 100 total reads across all samples were grouped together into a single taxonomic group. A total of 29,667 OTU consisting of 908,171 reads (1.2% of total reads; an average of 30 reads per OTU) were merged into a single OTU and labeled “Low Abundance OTU Group.”

The final data set consisted of 90 total samples; 36 rhizosphere soil, of which 18 were DNA derived and 18 RNA derived, 26 control soil samples, of which 13 were DNA derived and 13 RNA derived, and 28 bulk soil samples, 14 of which were DNA derived and 14 RNA derived. These samples contained 17,896 unique OTU, and 33,732,811 total reads (64% RNA reads: 36% DNA reads). Notably, selecting an even number of samples at random from each treatment (rhizosphere, control, or bulk) did not alter the conclusions of the analysis, and we therefore present results based on all 90 samples.

Samples were rarified to 32,997 reads per sample, and alpha diversity was estimated using the phyloseq package “estimate_richness.” Statistical differences between community nucleic acid derivation and soil type were determined using ANOVA (analysis of variance). To calculate beta diversity, we used the package Phyloseq v1.30.0 (McMurdie and Holmes, 2013). Reads assigned to an OTU in a sample were divided by the total number of reads per sample to calculate within site proportional abundances. MDS plots were used to visualize Bray-Curtis pairwise dissimilarities of community data. Significant differences between DNA and RNA derived communities, as well as the three soil types within and between communities were determined using PERMANOVA (pairwise adonis testing with Bonferroni correction). The corncob package (Martin et al., 2020) was used to estimate differential abundances between the DNA and RNA-derived communities, and between the DNA-derived and RNA-derived rhizosphere, bulk, and control soil samples, using the absolute abundance of reads, rather than the proportional abundances. Differentially abundant taxa with p-values < 0.01 were considered statistically significant.

T-tests of 16S-ratios and differential abundance analysis of the RNA-derived communities are two methods which can be used to determine elevated levels of RNA-derived reads, and therefore protein synthesis potential (PSP), between groups. In this study, we compared both analyses for overarching patterns of PSP in the microbial communities (see Supplementary Tables S1, S2). 16S-ratios were generated by normalizing the number of RNA reads (representing the number of ribosomes in a cell) by the number of DNA reads (representing the number of gene copies that code for ribosomes in a cell), using method 3 outlined in Bowsher et al. (2019). Though 16S-ratios are a common method of estimating overarching patterns of potential microbial activity in a community (Steven et al., 2017), we focused primarily on analyzing the RNA generated data, due to biases associated with 16S ratios such as (1) the within cell variability 16S rRNA-gene copy number in different taxa, causing over estimations of the 16S-ratio, which cannot be corrected for (Schaechter et al., 1958; Cooper and Helmstetter, 1968; Klappenbach, 2001; Lee et al., 2009; Franklin et al., 2013; Louca et al., 2018), and (2) using a single arbitrary cutoff point of 16S-ratio activity, which is problematic in diverse microbial communities and may include dormant microbes in analysis (Jones and Lennon, 2010; Blagodatskaya and Kuzyakov, 2013; Blazewicz et al., 2013; Steven et al., 2017). Mean 16S-ratios were calculated for each genus (Bowsher et al., 2019), in order to be compared to the genus-level differential abundance analysis output from corncob. We again used corncob to estimate RNA-based differential abundance. Though there are some concerns associated with directly using RNA as an indicator of PSP (outlined in depth in Blazewicz et al., 2013) corncob corrects for errors inherent in microbial data analyses. Corncob correlates taxa to covariates of interest and infers possible taxa presence in samples with small sequencing depth, with calculated variance around that inference, thereby reducing sampling bias (Willis, 2019; Martin et al., 2020). It is likely that this analysis better reflects which genera have a higher number of RNA reads compared to DNA, and therefore which genera have higher PSP in the different soil types. This study therefore focuses on reporting the RNA-derived community results, with the 16S-ratio analysis as Supplementary material.

Hypothetical functional gene profiles were created using the package Tax4Fun (Aßhauer et al., 2015), the output of which shows what percentage of the sample being analyzed is associated with a known gene and function. The RNA-derived genera that were significantly different in abundance among the soil environments in the corncob analysis were included in the Tax4Fun analysis, as were genera with a mean 16S-ratio value over one deemed significantly different between soil types using t-tests. For the Tax4Fun analysis, taxonomy of the OTU were reassigned using the Silva123 library (Quast et al., 2015) to ensure data compatibility. Gene profiles were compared to the KEGG Orthology database (Kanehisa et al., 2016) to determine the metabolic pathways with which the identified genes were associated. Mean Tax4Fun abundances were calculated for each genus in the different sample types, and the differences between genera identified as enriched vs. depleted were calculated to determine whether a hypothetical gene was more abundant in a certain soil treatment compared to the other. For example, in taxa identified as significantly different between rhizosphere vs. control soils, relative gene abundances were calculated for genera depleted in the rhizosphere and genera enriched in the rhizosphere soil. Genes with a positive value were more abundant in the “enriched (corncob −log2(FoldChange) > 0, or 16S-ratio p < 0.01 & t > 0)” rhizosphere soil category, and genes with a negative value were more abundant in the “depleted (corncob −log2(FoldChange) < 0, or 16S-ratio p < 0.01 & t < 0)” rhizosphere soil category (Figure 1).

After quality filtering and removal of chimeras, an average of 374,809 reads per sample were retained with 17,896 unique taxa in 90 samples (45 DNA, 45 RNA). Of those taxa 17,851 were present in the RNA community (16 taxa unique to RNA), and 17,880 in the DNA community (45 taxa unique to DNA).

Using rarefied communities, alpha diversity did not differ significantly between the communities characterized by DNA and RNA, in richness, Shannon, or Simpsons diversity on average between the bulk soil or control treatments. However, between the DNA and RNA, rhizosphere soils differed significantly from one another in Shannon (p < 0.01**) and Simpsons diversity (p < 0.1*) (Figure 2), indicating that DNA vs. RNA characterization capture different diversity among member vs. potentially active microbes.

Visualization of the relative abundance of some representative phyla between the DNA and RNA derived communities showed differering patterns of PSP vs. microbial membership. Bacteriodetes was not significantly differ between communities, indicating an equal amount of PSP for the number of cells present in the community (Figures 3A, 4). Two phyla (p: Gemmatimonadetes, p: Saccharibacteria) were prevalent in the DNA-derived community compared to the RNA-derived community (Figures 3A,C,D), indicating little to no PSP for how present the two phyla were in the community. Proteobacteria showed significantly greater number of RNA- to DNA-derived reads, indicating a high level of PSP for number of cells present in the community (Figures 3A, 4). Visualization of the absolute abundance of phyla between the DNA and RNA derived communities show 21 phyla as being significantly different in abundance between the RNA and DNA communities, with 13 of those significantly more abundant in the RNA-derived community and the remaining 6 more abundant in the DNA-derived community (Figure 4). Visualization of the 16S-ratio PSP compared to the RNA-derived community (which indicates PSP) showed similar patterns of phyla abundance, though 16S-ratios showed a much higher PSP of p: Cyanobacteria (Figure 3B), despite making up less than 1% of both the total RNA- and DNA-derived community (Figures 3B,C).

Based on Bray-Curtis analysis of the normalized data, beta diversity was significantly different between the overall DNA- and RNA-characterized communities (p < 0.001***; R2 = 0.104) (Figure 5A). Differential abundance analysis via corncob identified 11 phyla (137 genera) that differed in abundance based on whether they were characterized via DNA or RNA (p < 0.01). Of these, 71 genera were more abundant in the RNA- compared to the DNA-derived community estimation (Figure 1A).

The three soil environments compared in this study were rhizosphere soil, control-bulk soil with a similar physical matrix as the soils that plants were grown in, and bulk soil taken directly from the study site.

Some expected artifacts of RNA- and DNA-derived community analysis were found, though these made up <1% of the total reads used for analysis and the three soil types do not significantly differ in OTU unique to each soil type. In the DNA-derived community reconstruction 62 taxa and 3,286 reads (<1% of DNA reads) were unique to the rhizosphere soils, 12 taxa and 64 reads (<1% total DNA reads) were unique to control soils, and 256 taxa and 12,016 reads (<1% of total DNA reads) were unique to the bulk soils. In the RNA-derived community 38 taxa and 227 reads (0.1% of RNA reads) were unique to the rhizosphere soils, 76 taxa and 629 reads (<1% total RNA reads) were unique to control soils, 245 and 15,765 taxa and reads (<1%% of total RNA reads) were unique to the bulk soils.

Within the DNA-derived community, significant differences were seen between the bulk and rhizosphere soils (p < 0.1*) and control and rhizosphere soils (p < 0.01**) in richness, between the bulk and rhizosphere soils (p < 0.1*) and control and rhizosphere soils (p < 0.1*) in Shannon Diversity Indices, and between the bulk and rhizosphere soils (p < 0.0001***) and bulk and control soils (p < 0.01**) in Simpsons Diversity Indices (Figure 2). Based on Bray-Curtis analysis of normalized data from the DNA-derived community, the three soil environments were significantly different from each other (p < 0.001***, R2 = 0.398) (Figure 5D). Differential abundance analysis via corncob identified 144 genera that differed significantly between rhizosphere and bulk soils, with 128 genera being more abundant in rhizosphere than bulk soils (p < 0.001***) (Figure 1B). Between rhizosphere and control soils, 83 total genera were differentially abundant, with 73 being more abundant in the rhizosphere than control soils (Figure 1C). Between the control and bulk soils, 155 total genera were differentially abundant, with 135 genera more abundant in the control than bulk soils (Figure 1D).

Within the RNA-derived community composition, no significant differences were seen between the three soil communities in terms of richness, Shannon Diversity, or Simpson Diversity (p = 0.185) (Figure 2). Based on Bray-Curtis analysis of normalized data from the RNA group, the three soil environments were significantly different from each other (p < 0.001***, R2 = 0.4169) (Figure 5C). Differential abundance analysis via corncob identified 168 genera that differed significantly between rhizosphere and bulk soils, with 120 genera being more abundant in rhizosphere than bulk soils (p < 0.01) (Figure 1B). Between rhizosphere and control soils, 113 total genera were differentially abundant, with 79 more abundant in the rhizosphere than control soils (Figure 1C). Between the control and bulk soils, 179 total genera were differentially abundant, with 129 genera being more abundant in the control than bulk soils (Figure 1D). T-tests of significant differences between 16S-ratios of genera between soil types conflicted with corncob analysis of RNA-derived communities (Supplementary Table S1).

The genera of the RNA-derived community characterization that showed differential abundances between the rhizosphere and control soil environments were assigned hypothetical functional profiles using the software Tax4Fun (Aßhauer et al., 2015), as were the genera deemed by convention significantly different between soil environments using 16S-ratio analysis T-tests (Figure 6; Supplementary Table S2). Tax4fun predicts functional capabilities of microbial communities using read abundance of 16S datasets, and returns what percentage of a community is associated with a particular gene, which relates to a metabolic capability.

In rhizosphere vs. control soils, profiles of hypothetical metabolism between the corncob and 16S-ratio t-test genera do follow some similar patterns. For example, the corncob analysis indicated a higher proportion of the community of genera enriched in the rhizosphere is dedicated to amino acid and carbohydrate metabolism compared to other forms of metabolism (Figure 6A), and the 16S-ratio t-test agrees, though it reverses which of these categories are more prevalent (Figure 6B). Genera derived from corncob analysis (Figure 6A) also indicate that lipid metabolism is prevalent in the rhizosphere compared to control soils, which we do not see in the 16S-ratio derived soils (Figure 6B). Genes associated with the metabolism of various compounds such as galactose were highly abundant in the rhizosphere soils compared to control soils, in both the corncob and 16S-ratio t-test analysis. Genes associated with fatty acid degradation were only visibly prevalent in the corncob genera analysis (Figure 6).