Global change experiments often observe shifts in bacterial community composition based on 16S rRNA gene sequences. However, this genetic region can mask a large amount of genetic and phenotypic variation among bacterial strains sharing even identical 16S regions. As such, it remains largely unknown whether variation at the sub-16S level, sometimes termed microdiversity, responds to environmental perturbations and whether such changes are relevant to ecosystem processes. Here, we investigated microdiversity within Curtobacterium, the dominant bacterium found in the leaf litter layer of soil, to simulated drought and nitrogen addition in a field experiment. We first developed and validated Curtobacterium-specific primers of the groEL gene to assess microdiversity within this lineage. We then tracked the response of this microdiversity to simulated global change in two adjacent plant communities, grassland and coastal sage scrub (CSS). Curtobacterium microdiversity responded to drought but not nitrogen addition, indicating variation within the genus of drought tolerance but not nitrogen response. Further, the response of microdiversity to drought depended on the ecosystem, suggesting that litter substrate selects for a distinct composition of microdiversity that is constrained in its response, perhaps related to tradeoffs in resource acquisition traits. Supporting this interpretation, a metagenomic analysis revealed that the composition of Curtobacterium-encoded carbohydrate-active enzymes (CAZymes) varied distinctly across the two ecosystems. Identifying the degree to which relevant traits are phylogenetically conserved may help to predict when the aggregated response of a 16S-defined taxon masks differential responses of finer-scale bacterial diversity to global change.
IMPORTANCE Microbial communities play an integral role in global biogeochemical cycling, but our understanding of how global change will affect microbial community structure and functioning remains limited. Microbiome analyses typically aggregate large amounts of genetic diversity which may obscure finer variation in traits. This study found that fine-scale diversity (or microdiversity) within the bacterial genus Curtobacterium was affected by simulated global changes. However, the degree to which this was true depended on the type of global change, as the composition of Curtobacterium microdiversity was affected by drought, but not by nitrogen addition. Further, these changes were associated with variation in carbon degradation traits. Future work might improve predictions of microbial community responses to global change by considering microdiversity.


Microbial responses to global change experiments are often measured by observing shifts in the composition of operational taxonomic units (OTUs) based on 16S rRNA gene sequences (14). However, a high degree of genetic variation exists within such OTUs, variation here referred to as microdiversity (58). It remains unclear whether this variation among microdiverse lineages contributes to differential responses to environmental factors and whether such trait variation manifests at the functional level. For instance, microbiomes with disparate OTU composition can encode similar functions, suggesting few trait differences among 16S-defined taxa (911). On the other hand, taxa with identical 16S rRNA gene sequences can vary in ecologically important traits (1216).
The importance of microdiversity for global change responses will likely depend on the traits involved (17). If an environmental perturbation selects for traits that are phylogenetically conserved at broader taxonomic levels, then the OTU designations would capture the consistent responses of the underlying microdiversity. Conversely, if a change selects on traits that vary at finer taxonomic levels, then 16S-level resolution will mask differential responses within the OTU. Therefore, identifying which bacterial traits are responsible for responses to environmental change and their degree of phylogenetic conservation will determine the level of genetic resolution needed to accurately assess compositional responses (1820).
The Loma Ridge Global Change Experiment (LRGCE) in southern California simulates anthropogenically driven environmental changes expected in this region by experimentally manipulating rainfall and nitrogen in two adjacent plant communities: grassland and coastal sage scrub (CSS). Drought is projected to increase in intensity and frequency (21), while soil nitrogen deposition has steadily increased for decades due to human activities (22). Previous work at this site has shown that drought and nitrogen addition alter both bacterial and fungal community composition in the leaf litter (23, 24). However, these studies focused on the compositional responses at the community level, as assessed by 16S amplicon sequencing. To investigate whether bacterial microdiversity also responded to global change treatments, we revisited this site and focused on the response of one of the most abundant taxa.
Curtobacterium is a genus of Gram-positive, putative aerobic Actinobacteria highly abundant in leaf litter, the topmost layer of soil (25). Though originally known as a plant pathogen (26), Curtobacterium appears to be a major litter decomposer in southern Californian ecosystems and responds to both global change treatments at the LRGCE (6). Within the genus, the 16S region is highly conserved across this genus such that only four Curtobacterium OTUs are delineated at 100% sequence similarity of the hypervariable V4 region of the 16S region. These four OTUs comprise a high degree of genomic diversity, sharing <80% average nucleotide identity and representing at least nine distinct subclades based on genotypic and phenotypic characteristics (27). Using 16S amplicon sequencing, we found that the abundance of Curtobacterium (relative to all other bacterial taxa) increased in response to drought while exhibiting variable responses to nitrogen addition in the grassland and CSS ecosystems (28). Only two 100% OTUs (or ESVs, exact nucleotide variants) were detected, and one of these made up 98.9% of all the Curtobacterium sequences. However, using metagenomic sequencing, we also found that at least six subclades within the genus coexist at the field site (27). Thus, a 16S-level analysis cannot resolve within-Curtobacterium variation in the responses to the global change treatments.
To investigate the response of Curtobacterium microdiversity to simulated global change, we first needed to develop a new genetic marker to delineate and track the finer-scale genetic diversity within the genus. We created Curtobacterium-specific primers of the groEL gene, which after in silico analyses is expected to capture a higher degree of genetic variation than the 16S region. Thus, our first goal was to evaluate the groEL gene as a marker to capture Curtobacterium microdiversity and to distinguish between subclades defined by multiple marker genes. To assess the validity of this method, we compared the relative abundances of subclades derived from our groEL approach to those previously observed using shotgun metagenomes (6).
After initial validation of the groEL marker, we then tested two hypotheses. Based on the high degree of genomic variation within Curtobacterium, we first hypothesized that the previously observed response of Curtobacterium to simulated global changes masks disparate responses of its microdiversity. For instance, although Curtobacterium increased in relative abundance in response to drought overall, this aggregate response may obscure that some subclades responded differentially, with some being more drought tolerant than others. Moreover, we suspected that the response of Curtobacterium microdiversity to specific global change treatments will differ depending on the depth of conservation of the traits involved. In particular, drought responses of soil bacteria generally appear to be more deeply conserved than the 16S OTU level (20), whereas nitrogen assimilation is often highly variable; that is, it is shallowly conserved at the sub-OTU level (29, 30). In this case, Curtobacterium microdiversity would be altered by the nitrogen treatment, but not the drought treatment.
Second, we hypothesized that the response of Curtobacterium microdiversity to global change depends on other traits within the genus, especially those related to resource utilization strategies (31). Therefore, we expected that the microdiversity response to global change would depend on the ecosystem (grassland or CSS) because of differences in the leaf litter substrate. While both ecosystems experience similar abiotic conditions due to their close proximity, their leaf litter differs in polysaccharide content and selects for distinct bacterial communities (28). Genes encoding the breakdown of leaf litter, the carbohydrate-active enzymes (CAZymes), vary in their distribution and diversity among Curtobacterium subclades and are correlated with differential degradation rates of cellulose and xylan, two abundant carbohydrates in leaf litter (6). We thus further explored the role of CAZymes in this hypothesis by comparing the CAZyme diversity of fully sequenced Curtobacterium genomes to the response of Curtobacterium CAZymes from shotgun metagenomic libraries at the LRGCE.


Comparison of Curtobacterium microdiversity by two methods.

The groEL region captured a high degree of Curtobacterium microdiversity (12,228 exact sequence variants, including 6,184 singletons) while also providing a useful taxonomic marker of subclade diversity. For comparison, 16S sequencing of the same samples detected only two Curtobacterium exact sequence variants (ESVs). We also attempted to assess the microdiversity response to the LRGCE treatments using the 16S data. While the first ESV was found in all samples, sometimes as high as 30% relative abundance, the second ESV was found only in four samples, making any further 16S microdiversity analysis unfruitful.
To evaluate the success of groEL as a marker for microdiversity within Curtobacterium, we used a two-pronged approach. First, we generated a phylogeny based on the targeted groEL region. The groEL phylogeny was highly congruent with a multilocus phylogeny based on 21 single-copy core genes (Fig. 1a; Mantel test, z = 0.749, P = 0.001). The genus can be divided into at least nine subclades, all of which could be distinguished by the targeted portion of the groEL gene. The overall topology of the two trees was also highly similar, but there were some differences in branching patterns (Fig. 1a). In particular, the clustering of the strains designated previously as clade II (light yellow in Fig. 1a) was uncertain, as was the placement of the two other strains (Wood-2 and Pine-20).
FIG 1 (a) Phylogenetic trees of the Curtobacterium genus created using 21 single-copy core genes (left [27]) and using only the groEL gene (right). The subclades are denoted by color, and at the top is the outgroup, Frigoribacterium, a closely related genus. Clades with bootstrap support of >70 are noted with a black circle at the node. (b) The relationship between the relative abundance of each subclade from the metagenomic analysis versus the relative abundance from the groEL amplicon analysis across the climate gradient. Each point represents a different sample that is colored by subclade. The overall correlation (dashed line) across all samples is plotted as reported in the top right of the figure. (c) Density plots of the relative abundance of Curtobacterium detected from the climate gradient samples. The groEL amplicon abundances are outlined by a solid line, and the metagenomic abundances, by a dotted line.
Our second test of the groEL marker was to compare the relative abundances of the subclades detected by classification of core genes in metagenomes versus groEL amplicons from the same samples. Approximately 52% of the groEL amplicon reads were classified as Curtobacterium. These sequences revealed a rough correspondence of the relative abundance of Curtobacterium subclades with that assessed by the metagenomic sequences in that the same subclades dominated the samples in both methods. For some subclades, such as II and VA, there was a close match in the relative abundance predicted by the two methods (Fig. 1c). Further, the three relatively low-abundance subclades (IIA, IIB, and IIIA) showed similarly low abundances by both methods; however, the relative abundances of closely related subclades IVA and IVB were reversed by the two methods (Fig. 1c). Lastly, while there was a weak overall correlation (r = 0.209, P < 0.0001), the relative abundances estimated by the two methods within subclades were not correlated (Pearson’s R range from −0.25 to 0.33), except for a significantly positive correlation within subclade IIA (P = 0.011; Fig. 1b). Thus, while the groEL approach provides a high-throughput method to track changes in relative Curtobacterium microdiversity across samples, it is likely not a good quantitative estimator of subclade abundance.

Response to the global change treatments and ecosystem type.

The composition of Curtobacterium microdiversity varied significantly between the drought and ambient treatments, across the two ecosystems, and over time. In particular, the composition of groEL ESVs within the genus responded to the drought treatment, explaining 2.8% of the variation (permutational multivariate analysis of variance [PERMANOVA], P = 0.001; Table 1; Fig. 2). As we know more about subclade distribution, we further analyzed at this level. Summarizing microdiversity by subclade, subclade IVA had a higher relative abundance under ambient conditions, whereas subclades IVB, IVC, and VA were more abundant in the drought treatment (Fig. 3a). In contrast to drought, added nitrogen did not alter groEL ESV composition.
FIG 2 Nondimensional metric scaling (NMDS) plots of Curtobacterium microdiversity (ESVs of groEL amplicons) from the LRGCE samples using Bray-Curtis distance. Vectors on the bottom left represent the direction and strength of correlation with the Curtobacterium subclades. The centroids of each treatment combination are marked by a black circle, and the centroids of all samples in the two ecosystems are marked by an X, showing the greater effect of drought in the grassland.
FIG 3 (a and b) Percent difference (mean ± standard error [SE]) in relative abundance of Curtobacterium subclades in (a) the drought versus ambient rainfall treatments and (b) the CSS versus grassland ecosystems (bottom panel). Subclades are colored as in Fig. 1. The percent change value for drought is calculated as the relative abundance under drought conditions minus the relative abundance under ambient conditions, divided by the relative abundance in drought. A positive value in the drought panel means higher relative abundance under drought conditions.
TABLE 1 PERMANOVA results evaluating the contribution of date of sample collection, ecosystem (grassland or CSS), drought treatment, and nitrogen treatment to variation of Curtobacterium composition (as assayed by groEL ESVs)a,b
SourcedfSSMSPseudo-FP(Perm)Estimated variation (%)
Ecosystem × drought10.
Ecosystem × nitrogen10. 
Ecosystem × date62. 
Drought × nitrogen10. 
Drought × date62. 
Nitrogen × date61. 
Residuals10437.10.3  82.3
df, degrees of freedom; SS, sums of squares; MS, mean squares.
Variation estimates are reported for statistically significant variables (indicated in bold).
Ecosystem (grassland or CSS) contributed the largest amount of variation in Curtobacterium microdiversity, accounting for 10.4% of the total variation (P = 0.001), and the temporal sampling accounted for an additional 1.8% (P = 0.004). These differences across ecosystems were largely driven by particular subclade preferences; for instance, subclade IVA was 150% more abundant in the grassland, while clade II and associated subclades (IIA and IIB) preferred the CSS (Fig. 3b). Ecosystem and drought also interacted to explain 2.7% of the variation in Curtobacterium ESV composition (PERMANOVA, P = 0.002). This interaction can be visualized in an ordination plot by the pronounced effect from the drought treatment on the grass samples but had little discernible effect on the CSS samples (Fig. 2, ecosystem centroids).

Genomic and metagenomic CAZyme variation.

While the global change treatments influenced subclade distributions, the majority of variation was explained by the two ecosystems, likely reflecting differences in litter chemistry. Therefore, we targeted functional traits directly associated with carbohydrate degradation by investigating the CAZyme genomic content across subclades. The number of CAZymes across all Curtobacterium genomes (n = 143) ranged from 61 to 84 genes, with a median of 70 CAZymes per genome. Across subclades, there was significant variation in the total number of CAZymes (see Fig. S1 in the supplemental material; Kruskal-Wallis P < 0.001), as well as significant compositional variation (analysis of similarity [ANOSIM], R = 0.685, P = 0.001). Major CAZymes driving the differences between some of the subclades include GH158, found in subclades IVA and IVB and thought to hydrolyze fungal laminarin, and GH116, found in subclade VA and thought to hydrolyze glucose and xylan (Fig. 4).
FIG 4 Heatmap of normalized mean CAZyme abundances of Curtobacterium clades. Data were normalized to a mean of 0 and a standard deviation of 1. CAZymes and clades were grouped by hierarchical clustering using Ward’s method and Euclidean distances.
Lastly, we investigated CAZyme variation in Curtobacterium within the LRGCE treatments in metagenomic libraries. There was no significant difference in the overall CAZyme count between ecosystem or drought treatment. However, as with groEL ESV composition, the composition of Curtobacterium CAZymes varied significantly between the drought and ambient treatments, across the two ecosystems, and over time. The largest effect was from ecosystem, also consistent with groEL composition, explaining 18.6% of the total variation (Table 2; PERMANOVA P = 0.001). CAZyme composition was highly distinct between ecosystems (Fig. 5; P = 0.001), and this strong division was driven by small but highly consistent differences in the relative abundances of many CAZymes (Fig. S2). There was also variation in dispersion, which can contribute to significant compositional effects in a PERMANOVA test (32), with significantly less dispersion in the grass ambient samples than the CSS ambient, grass drought, or CSS reduced samples (PERMDISP; P = 0.001; Fig. 5b).
FIG 5 (a) Nondimensional metric scaling (NMDS) plot showing the differences in genomic CAZyme content among Curtobacterium strains, colored by subclade. Overlaid are vectors whose direction and magnitude reflect the correlation of the CAZymes to the NMDS axes. (b) NMDS plot of Curtobacterium CAZyme content from metagenomic samples in the experimental plots at the LRGCE. Vectors represent the correlation of CAZymes to the NMDS axes.
TABLE 2 PERMANOVA results evaluating the contribution of ecosystem, drought treatment, and date of sample collection on Curtobacterium CAZyme composition from metagenomic samplesa,b
SourcedfSSMSPseudo-FP(Perm)Estimated variation (%)
Ecosystem × drought1438.4438.310.90.00114.6
Ecosystem × date3447.9149.33.60.0017.8
Drought × date3298.399.42.50.0014.2
Ecosystem × drought × date3261.687.22.20.0016.7
Residuals893,592.440.4  36.3
df, degrees of freedom; SS, sums of squares; MS, mean squares.
Variation estimates are reported for statistically significant variables (indicated in bold).


Closely related bacteria co-occur in the same habitat (30, 33), and such microdiversity suggests that coexistence may be promoted by fine-scale variation in traits (12, 14, 34). Overall, our results indicate that aggregating genetic diversity at the 16S level at least partially masks disparate responses to global change at finer scales of genetic variation. Previous work found that the relative abundance of Curtobacterium increased under drought conditions (6, 28). We find further variability in the response to drought within the genus, supporting our first hypothesis that microdiversity responds disparately to the global change simulations. However, in contrast to our expectation that the response to added nitrogen would be stronger than that to drought, Curtobacterium microdiversity did not respond to nitrogen addition. Nitrogen addition often affects bacterial community composition in other soil systems at the 16S level (3537), including in these same samples from the LRGCE system (28). This result suggests that the response to nitrogen addition may be more conserved than that of drought within Curtobacterium, in contrast to a meta-analysis that found that the phylogenetic depth of response across all bacteria was similar for various simulated global changes (20).
A main reason that we selected the groEL gene as a taxonomic marker is because its phylogeny and (presumably) evolutionary history are well matched to that of the 16S region and other markers. This suggests that it evolves relatively slowly compared to most parts of the genome and is not subject to high levels of horizontal gene transfer (HGT). Indeed, we found that variation in the groEL sequence reflected the phylogenetic distinctiveness of Curtobacterium subclades, indicating that shifts in groEL diversity will generally reflect relative shifts in Curtobacterium microdiversity. A further benefit of using a protein-coding gene over the 16S region to study microdiversity is that much of the genetic diversity will encompass synonymous mutations (in many first and third codon positions). In contrast, the 16S region is famously slow in its evolution, presumably because much of it is under stabilizing selection (38). The high number of synonymous mutations thus allows protein-coding regions such as groEL to resolve finer-scale differences in diversity than even variable regions of the 16S region can. Of course, selection for particular groEL genes in our treatments is possible; the groEL gene codes for a stress response chaperonin, and mutations therein may confer tolerance to different stress conditions (39, 40). However, the vast majority of genetic variability observed in groEL diversity is likely neutral with respect to our treatments and, because HGT of the gene is rare, linked to other regions in the genome that are under selection. Thus, we conclude that targeting groEL to assay shifts in fine-scale genetic composition is a useful approach for microdiversity exploration (16, 4143).
Our second hypothesis, that the response of microdiversity would depend on other traits, was supported by the large ecosystem effect, which mirrors the disparate chemical signatures of the litter substrates between the grassland and CSS ecosystems (31). This is further supported by the interactive effect between treatment and ecosystem, which is larger even than the main drought effect. The microdiversity response may be influenced by selection on other shallowly conserved traits such as CAZymes that might be easily shared via horizontal gene transfer among closely related strains within a genus. For example, subclade IVA encodes a high proportion of chitinases and laminarinases compared to other clades, and the target substrate for these enzymes is likely fungal cell walls, as their major component is β-1,3/1,6-glucan (44). Thus, specialization on fungal necromass versus other substrates may allow for resource partitioning to promote the coexistence of microdiversity. We suspect that these functional traits contribute to environmental distributions; in this case, subclade IVA was more abundant in the grassland, where the fungal to bacterial ratio is significantly greater than in the CSS (31). Supporting this idea, we further found that the composition of Curtobacterium-specific CAZymes responded strongly to drought as well as the ecosystem.
That said, there was not an obvious relationship between shifts in the abundance of particular subclades by treatment or ecosystem and the relative abundance of Curtobacterium CAZymes in metagenomic samples, even though the subclades differ distinctly in their genomic CAZyme content. This lack of correspondence between our results could suggest that we have not adequately isolated the genomic diversity of Curtobacterium and/or that groEL amplification is not a quantitative marker of the subclades. Indeed, groEL amplicon abundances correlated only weakly with Curtobacterium subclade abundances from the metagenomes. However, assigning metagenomic sequences at a fine phylogenetic scale is also challenging, so these may not reflect true subclade abundances.
While we did not observe a clear correlation between taxonomy and CAZyme composition, we did detect a strong response across ecosystems. Notably, the response (assessed by groEL) of Curtobacterium microdiversity to drought in the grassland was much larger than that in the CSS. Bacteria on the more recalcitrant CSS litter may be required to invest most of their cellular energy in degrading the leaves, perhaps revealing a trade-off in microbial life history strategies (4548). However, a similar interactive pattern was not reflected in the composition of Curtobacterium CAZymes, as these genes responded distinctly in both ecosystems, albeit in a different way. This discrepancy may indicate that the sensitivity of Curtobacterium microdiversity to drought depends not only on CAZymes but on additional traits not considered here.
In conclusion, the responses assessed at the 16S level reflect the summation of divergent responses among finer-scale taxa. However, it is important to distinguish between the type of information that the 16S and groEL analyses provide. Because the groEL primers target only Curtobacterium, the analysis does not provide information about how Curtobacterium changes relative to all other bacterial taxa. Thus, the 16S analysis still provides unique and useful information about Curtobacterium, even when only one 16S OTU is effectively present. Separately, the groEL analysis demonstrates that it would be incorrect to extrapolate from the 16S data that all Curtobacterium respond relatively positively to drought. In contrast, the 16S responses to nitrogen addition seem to apply throughout the genus at our field site. The varied responses to simulated drought within a genus of soil bacteria demonstrate the potential importance of finely conserved traits—traits that are otherwise aggregated by the broader taxonomic shifts observed in most microbial global change studies. Of course, it is not feasible to assess every bacterial genus’ response to global change at a finer and finer scale of diversity. However, future investigations from closely studied models such as Curtobacterium could identify traits involved in responses to environmental change, and the genetic resolution is needed to assess them. Such studies could thereby improve predictions about global change responses of the wider microbial community and consequences for soil ecosystems.


Field sites.

The Loma Ridge Global Change Experiment (LRGCE) was established in Irvine, CA, USA (33°44′N, 117°42′W, 365 m elevation) in 2007. Precipitation and nitrogen treatments are applied in two adjacent ecosystems, deciduous shrubland (coastal sage scrub, CSS) and annual grassland (49). The grassland plots are dominated by Avena, Bromus, and Lolium, and the CSS plots are dominated by Artemesia and Salvia (50). The site contains 48 experimental plots, 6.1 by 12.2 m in the grassland, and 18.3 by 12.2 m in the CSS. For this study, we sampled from four replicate plots from four treatments (drought, added nitrogen, drought and added nitrogen, and control) in both the grassland and CSS habitats. Since 2007, the drought treatment (∼50% reduction of annual ambient rainfall) was achieved by intercepting rainfall with a large polyethylene sheet during approximately half of the rainstorms. Nitrogen addition was achieved by adding soluble CaNO3 (60 kg N−1 ha year−1 added) (49).

groEL sequencing and analysis.

To identify a marker gene that enabled the characterization of the microdiversity within Curtobacterium, we initially identified 12 orthologous proteins that could distinguish between Curtobacterium subclades. Based on the feasibility to amplify the flanking regions of these candidate genes, we selected the groEL gene that encodes part of the groEL-groES chaperonin complex (51). Notably, the groEL gene has been used as an alternative to small subunit rRNA genes in other taxa (52). The final primer design targeted a 20-bp region at each end—(Curto-groEL-72F) CGCCGTGAAGGTGACGCTCG and (Curto-groEL-450R) GCCGAGATCGACGCSGTSGC. To test the specificity of the primers, we confirmed amplification of the correct size band from six Curtobacterium isolates spread across the phylogeny. We also tested that we did not amplify the region from leaf litter isolates of the closely related genus Frigoribacterium or other distantly related taxa such Pseudomonas (two genera that are commonly found in the LRGCE leaf litter communities).
Curtobacterium microdiversity was characterized on 219 leaf litter samples collected from 3 areas of each treatment plot at 7 time points (once per season) over a year and a half—22 August 2016, 12 December 2016, 30 March 2017, 28 June 2017, 13 September 2017, 13 December 2017, and 29 March 2018. In the lab, the leaf litter was ground using a coffee grinder, and DNA was extracted from 0.05 g using the ZymoBIOMICS DNA miniprep kit (Zymo Research; Irvine, CA). We added 3 μL of this DNA to a 25-μL PCR cocktail containing 12.5 μL AccuStart II PCR ToughMix 2× reagent (Quantabio; Beverly, MA) and 1 μL each of the groEL primers, and the remainder was made up with nuclease-free water. PCR conditions were as follows: 94°C for 5 min, followed by 27 cycles of 94°C for 45 s, 68°C for 45 s, and 72°C for 45 s, with a final extension of 72°C for 10 min. The amplicons were purified using a 1:1 ratio of AMPure beads and eluted into 20 μL nuclease-free water. Then, 1 μL each of unique i5 and i7 Nextera (Illumina, Inc., San Diego, CA) indices was added to each sample, followed by a second PCR under the following conditions: 94°C for 3 min, followed by 8 cycles of 94°C for 30 s, 55°C for 30 s, and 72°C for 30 s, with a final extension of 72°C for 10 min. Amplified DNA was quantified using a Qubit instrument (BioTek, Winooski, VT) and assessed for quality using a Nanodrop instrument (Thermo Fisher, Waltham, MA). Samples were prepared according to the Nextera XT library preparation kit (Illumina, Inc.). The samples were then pooled equimolarly and cleaned up using a 1:1 ratio of AMPure beads, eluted into 100 μL, and sequenced on an Illumina MiSeq instrument with 300-bp paired-end reads, though the reverse reads were not used due to low quality.
We used a QIIME 2 v2018.11 (53) pipeline to denoise and trim unprocessed FASTQ sequence files with DADA2 (54). To classify groEL amplicons at the subclade level, we first clustered sequences to exact sequence variants (ESVs). To validate that ESVs belonged to Curtobacterium and were not erroneously assigned, we screened the amplicon sequences against a custom database (since we used a non-16S marker gene with custom primers) containing 70 groEL sequences extracted from publicly available genomes of Curtobacterium (27), as well as 10 other major taxa commonly found in leaf litter communities to serve as outgroups. Using a trained classifier from the feature-classifier QIIME plugin (55), we used phylogenetic inference of the groEL gene to place ESVs (56) into a maximum likelihood phylogeny generated from the 70 reference sequences built with RAxML v8.2.12 (57) under the GTRGAMMA model with 100 replicates.
To validate the groEL amplicon approach, we compared our results to previously published shotgun metagenomic data. While the groEL sequences provide a much finer resolution of genetic diversity, we first amplified the groEL region from 96 samples from a regional climate gradient in southern California that included the LRGCE. We compared the groEL diversity with the shotgun metagenomic results at the clade level (five clades originally reported in [27]). In an earlier study, Curtobacterium abundances from the metagenomic libraries were calculated with a multilocus sequence assignment using 21 single-copy genes specific to Curtobacterium (27). Separately, we placed the amplified groEL sequences on a phylogenetic tree of Curtobacterium groEL sequences using the SEPP plugin in QIIME2 (56). Approximately 52% of the ∼10,000 ESVs were assigned to Curtobacterium. Of these, 88% were assigned to the species level, with the remaining 12% only placing to genus level, indicating that most of the genetic diversity was captured by the Curtobacterium genomes in the phylogenetic tree. The samples had a mean of 51% Curtobacterium hits (±25.3% standard deviation), and the non-Curtobacterium reads were removed from further analysis. For both the metagenomes and our amplicon method, the relative abundance of each subclade was calculated by dividing the number of reads assigned to a subclade by the total number of assigned reads in that sample.

Curtobacterium CAZyme composition.

To assess CAZyme diversity within Curtobacterium genomes, we downloaded all available (n = 143) Curtobacterium genomes from NCBI and analyzed their CAZyme profile using the program dbCAN2 (58) with the HMMdb v9 database released on 4 August 2020 and the latest CAZyDB, released on 30 July 2020 (59). dbCAN2 uses a 3-fold approach to identify CAZymes, and we only retained CAZymes that were detected by all three approaches. CAZyme abundances were normalized across genomes by calculating a Z-score with a mean of 0 and a standard deviation of 1, after which we performed an nondimensional metric scaling (NMDS) analysis and plotted the correlation vectors.
From the detected Curtobacterium CAZymes in the genomes, we built a custom database to identify the abundance of Curtobacterium CAZymes in metagenomic sequences from the LRGCE treatments. The sequences were checked for quality using FastQC v0.11.9 (60), after which they were quality filtered and trimmed using Trimmomatic v0.39 (61). We then used DIAMOND BLASTX (62) on the forward reads, taking into account all six putative open reading frames, for comparison to the reference database. To assess the validity of hits, we used CAZymes from the closely related sister genus Frigoribacterium to benchmark our thresholds and set positive hits at an E value of 1e−20 and an identity of >98%. Finally, we performed a reciprocal BLAST search of the putative hits against the full NCBI nonredundant protein database to confirm the reads identified as Curtobacterium.

Statistical analyses.

To assess the contribution of the global change treatments on Curtobacterium composition, we performed a permutational multivariate analysis of variance using PRIMER6 with PERMANOVA+ (Primer‐E Ltd., Ivybridge, UK) (63). The groEL PERMANOVA model included four factors, ecosystem, date of sample collection, precipitation treatment, and nitrogen treatment, as fixed effects. Compositional similarity between samples was calculated using the Bray-Curtis dissimilarity metric after rarefying to a sample depth of 5,000 sequences with 100 repetitions. We ran a type III partial sum of squares PERMANOVA for 999 permutations of residuals. Variation explained was calculated by summing the estimated components of variation for the statistically significant terms and the residuals and dividing each by this total (64). We followed the same process for the metagenomic PERMANOVA model, although there was no nitrogen addition treatment. The Curtobacterium genomic CAZyme content data were normalized by centering to a mean of 0 and a standard deviation of 1 using the scale function in R. ANOSIM and SIMPER analyses to assess differences in subclade CAZyme composition used Bray-Curtis dissimilarities and were performed with the vegan package in R (65).

Data availability.

The raw amplicon and metagenomic reads can be accessed under the NCBI BioProject accession number PRJNA781975. The BioSample accession numbers are SAMN23310455 for the groEL amplicons and SAMN23390019 for the metagenomic reads.


UC Irvine and the LRGCE are located on the ancestral homelands of the Indigenous Kizh and Acjachemen nations.
We thank Alejandra Rodriguez Verdugo, Katrine Whiteson, Kendra Walters, Cynthia Rodriguez, Kristin Barbour, Alberto Barron Sandoval, Joanna Wang, Joia Kai Capocchi, Pauline Uyen Phuong Nguyen, Khanh Thuy Huynh, and Clara Barnosky for their input on analyses and previous drafts and for laboratory help.
This work was supported by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research grants DE-SC0016410 and DE-SC0020382.

Supplemental Material

File (aem.02429-21-s0001.pdf)
ASM does not own the copyrights to Supplemental Material that may be linked to, or accessed through, an article. The authors have granted ASM a non-exclusive, world-wide license to publish the Supplemental Material files. Please contact the corresponding author directly for reuse.


Liu J, Meng Z, Liu X, Zhang X-H. 2019. Microbial assembly, interaction, functioning, activity and diversification: a review derived from community compositional data. Mar Life Sci Technol 1:112–128.
Hutchins DA, Fu F. 2017. Microorganisms and ocean global change. Nat Microbiol 2:17058.
Valencia E, Gross N, Quero JL, Carmona CP, Ochoa V, Gozalo B, Delgado-Baquerizo M, Dumack K, Hamonts K, Singh BK, Bonkowski M, Maestre FT. 2018. Cascading effects from plants to soil microorganisms explain how plant species richness and simulated climate change affect soil multifunctionality. Glob Chang Biol 24:5642–5654.
DeAngelis KM, Pold G, Topçuoğlu BD, van Diepen LTA, Varney RM, Blanchard JL, Melillo J, Frey SD. 2015. Long-term forest soil warming alters microbial communities in temperate forest soils. Front Microbiol 6:104.
Johnson ZI, Zinser ER, Coe A, McNulty NP, Woodward EMS, Chisholm SW. 2006. Niche partitioning among Prochlorococcus ecotypes along ocean-scale environmental gradients. Science 311:1737–1740.
Chase AB, Karaoz U, Brodie EL, Gomez-Lunar Z, Martiny AC, Martiny JBH. 2017. Microdiversity of an abundant terrestrial bacterium encompasses extensive variation in ecologically relevant traits. mBio 8:e01809-17.
Larkin AA, Blinebry SK, Howes C, Lin Y, Loftus SE, Schmaus CA, Zinser ER, Johnson ZI. 2016. Niche partitioning and biogeography of high light adapted Prochlorococcus across taxonomic ranks in the North Pacific. ISME J 10:1555–1567.
García-García N, Tamames J, Linz AM, Pedrós-Alió C, Puente-Sánchez F. 2019. Microdiversity ensures the maintenance of functional microbial communities under changing environmental conditions. ISME J 13:2969–2983.
Yin B, Crowley D, Sparovek G, De Melo WJ, Borneman J. 2000. Bacterial functional redundancy along a soil reclamation gradient. Appl Environ Microbiol 66:4361–4365.
Wohl DL, Arora S, Gladstone JR. 2004. Functional redundancy supports biodiversity and ecosystem function in a closed and constant environment. Ecology 85:1534–1540.
Louca S, Polz MF, Mazel F, Albright MBN, Huber JA, O’Connor MI, Ackermann M, Hahn AS, Srivastava DS, Crowe SA, Doebeli M, Parfrey LW. 2018. Function and functional redundancy in microbial systems. Nat Ecol Evol 2:936–943.
Jaspers E, Overmann J. 2004. Ecological significance of microdiversity: identical 16S rRNA gene sequences can be found in bacteria with highly divergent genomes and ecophysiologies. Appl Environ Microbiol 70:4831–4839.
Brazelton WJ, Ludwig KA, Sogin ML, Andreishcheva EN, Kelley DS, Shen C-C, Edwards RL, Baross JA. 2010. Archaea and bacteria with surprising microdiversity show shifts in dominance over 1,000-year time scales in hydrothermal chimneys. Proc Natl Acad Sci USA 107:1612–1617.
Schloter M, Lebuhn M, Heulin T, Hartmann A. 2000. Ecology and evolution of bacterial microdiversity. FEMS Microbiol Rev 24:647–660.
Acinas SG, Klepac-Ceraj V, Hunt DE, Pharino C, Ceraj I, Distel DL, Polz MF. 2004. Fine-scale phylogenetic architecture of a complex bacterial community. Nature 430:551–554.
Braker G, Zhou J, Wu L, Devol AH, Tiedje JM. 2000. Nitrite reductase genes (nirK and nirS) as functional markers to investigate diversity of denitrifying bacteria in pacific northwest marine sediment communities. Appl Environ Microbiol 66:2096–2104.
Martiny JBH, Jones SE, Lennon JT, Martiny AC. 2015. Microbiomes in light of traits: a phylogenetic perspective. Science 350:aac9323.
Martiny AC, Treseder K, Pusch G. 2013. Phylogenetic conservatism of functional traits in microorganisms. ISME J 7:830–838.
Burke C, Steinberg P, Rusch D, Kjelleberg S, Thomas T. 2011. Bacterial community assembly based on functional genes rather than species. Proc Natl Acad Sci USA 108:14288–14293.
Isobe K, Bouskill NJ, Brodie EL, Sudderth EA, Martiny JBH. 2020. Phylogenetic conservation of soil bacterial responses to simulated global changes. Philos Trans R Soc Lond B Biol Sci 375:20190242.
IPCC. 2021. Climate change 2021: the physical science basis contribution of Working Group I to the sixth assessment report of the Intergovernmental Panel on Climate Change. https://www.ipcc.ch/report/ar6/wg1/.
Galloway JN, Townsend AR, Erisman JW, Bekunda M, Cai Z, Freney JR, Martinelli LA, Seitzinger SP, Sutton MA. 2008. Transformation of the nitrogen cycle: recent trends, questions, and potential solutions. Science 320:889–892.
Matulich KL, Weihe C, Allison SD, Amend AS, Berlemont R, Goulden ML, Kimball S, Martiny AC, Martiny JBH. 2015. Temporal variation overshadows the response of leaf litter microbial communities to simulated global change. ISME J 9:2477–2489.
Amend AS, Martiny AC, Allison SD, Berlemont R, Goulden ML, Lu Y, Treseder KK, Weihe C, Martiny JBH. 2016. Microbial response to simulated global change is phylogenetically conserved and linked with functional potential. ISME J 10:109–118.
Chase AB, Arevalo P, Polz MF, Berlemont R, Martiny JBH. 2016. Evidence for ecological flexibility in the cosmopolitan genus Curtobacterium. Front Microbiol 7:1874.
Wood BA, Easdown WJ. 1990. A new bacterial disease of mung bean and cowpea for Australia. Austral Plant Pathol 19:16–21.
Chase AB, Gomez-Lunar Z, Lopez AE, Li J, Allison SD, Martiny AC, Martiny JBH. 2018. Emergence of soil bacterial ecotypes along a climate gradient. Environ Microbiol 20:4112–4126.
Finks SS, Weihe C, Kimball S, Allison SD, Martiny AC, Treseder KK, Martiny JBH. 2021. Microbial community response to a decade of simulated global changes depends on the plant community. Elementa 9:00124.
Stewart CE, Roosendaal D, Denef K, Pruessner E, Comas LH, Sarath G, Jin VL, Schmer MR, Soundararajan M. 2017. Seasonal switchgrass ecotype contributions to soil organic carbon, deep soil microbial community composition and rhizodeposit uptake during an extreme drought. Soil Biol Biochem 112:191–203.
Moore LR, Rocap G, Chisholm SW. 1998. Physiology and molecular phylogeny of coexisting Prochlorococcus ecotypes. Nature 393:464–467.
Malik AA, Swenson T, Weihe C, Morrison EW, Martiny JBH, Brodie EL, Northen TR, Allison SD. 2020. Drought and plant litter chemistry alter microbial gene expression and metabolite production. ISME J 14:2236–2247.
Anderson MJ. 2006. Distance-based tests for homogeneity of multivariate dispersions. Biometrics 62:245–253.
Thompson JR, Pacocha S, Pharino C, Klepac-Ceraj V, Hunt DE, Benoit J, Sarma-Rupavtarm R, Distel DL, Polz MF. 2005. Genotypic diversity within a natural coastal bacterioplankton population. Science 307:1311–1313.
Kashtan N, Roggensack SE, Rodrigue S, Thompson JW, Biller SJ, Coe A, Ding H, Marttinen P, Malmstrom RR, Stocker R, Follows MJ, Stepanauskas R, Chisholm SW. 2014. Single-cell genomics reveals hundreds of coexisting subpopulations in wild Prochlorococcus. Science 344:416–420.
Swathi AT, Rakesh M, Premsai SB, Louis ST, William KT, Subhash CM. 2013. Chronic N-amended soils exhibit an altered bacterial community structure in Harvard Forest, MA, USA. FEMS Microbiol Ecol 83:478–493.
Fierer N, Lauber CL, Ramirez KS, Zaneveld J, Bradford MA, Knight R. 2012. Comparative metagenomic, phylogenetic and physiological analyses of soil microbial communities across nitrogen gradients. ISME J 6:1007–1017.
Lekberg Y, Arnillas CA, Borer ET, Bullington LS, Fierer N, Kennedy PG, Leff JW, Luis AD, Seabloom EW, Henning JA. 2021. Nitrogen and phosphorus fertilization consistently favor pathogenic over mutualistic fungi in grassland soils. Nat Commun 12:3484.
Wang H, Hickey DA. 2002. Evidence for strong selective constraint acting on the nucleotide composition of 16S ribosomal RNA genes. Nucleic Acids Res 30:2501–2507.
Sen R, Tripathy S, Padhi SK, Mohanty S, Maiti NK. 2014. Sequence polymorphism of GroEL gene in natural population of Bacillus and Brevibacillus spp. that showed variation in thermal tolerance capacity and mRNA expression. Curr Microbiol 69:507–516.
Kim H-W, Wi AR, Jeon BW, Lee JH, Shin SC, Park H, Jeon S-J. 2015. Cold adaptation of a psychrophilic chaperonin from Psychrobacter sp. and its application for heterologous protein expression. Biotechnol Lett 37:1887–1893.
Kent AG, Baer SE, Mouginot C, Huang JS, Larkin AA, Lomas MW, Martiny AC. 2019. Parallel phylogeography of Prochlorococcus and Synechococcus. ISME J 13:430–441.
Hu L, Lu W, Wang L, Pan M, Zhang H, Zhao J, Chen W. 2017. Assessment of Bifidobacterium species using groEL gene on the basis of Illumina MiSeq high-throughput sequencing. Genes 8:336.
Caro-Quintero A, Ochman H. 2015. Assessing the unseen bacterial diversity in microbial communities. Genome Biol Evol 7:3416–3425.
Fesel PH, Zuccaro A. 2016. β-glucan: crucial component of the fungal cell wall and elusive MAMP in plants. Fungal Genet Biol 90:53–60.
Malik AA, Martiny JBH, Brodie EL, Martiny AC, Treseder KK, Allison SD. 2020. Defining trait-based microbial strategies with consequences for soil carbon cycling under climate change. ISME J 14:1–9.
Schimel J, Balser TC, Wallenstein M. 2007. Microbial stress-response physiology and its implications for ecosystem function. Ecology 88:1386–1394.
Carlson RP, Taffs RL. 2010. Molecular-level tradeoffs and metabolic adaptation to simultaneous stressors. Curr Opin Biotechnol 21:670–676.
Strickland MS, Keiser AD, Bradford MA. 2015. Climate history shapes contemporary leaf litter decomposition. Biogeochemistry 122:165–174.
Potts DL, Suding KN, Winston GC, Rocha AV, Goulden ML. 2012. Ecological effects of experimental drought and prescribed fire in a southern California coastal grassland. J Arid Environ 81:59–66.
Allison SD, Lu Y, Weihe C, Goulden ML, Martiny AC, Treseder KK, Martiny JBH. 2013. Microbial abundance and composition influence litter decomposition response to environmental change. Ecology 94:714–725.
Zeilstra-Ryalls J, Fayet O, Georgopoulos C. 1991. The universally conserved GroE (Hsp60) chaperonins. Annu Rev Microbiol 45:301–325.
Vancuren SJ, Hill JE. 2019. Update on cpnDB: a reference database of chaperonin sequences. Database 2019:baz033.
Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet CC, Al-Ghalith GA, Alexander H, Alm EJ, Arumugam M, Asnicar F, Bai Y, Bisanz JE, Bittinger K, Brejnrod A, Brislawn CJ, Brown CT, Callahan BJ, Caraballo-Rodríguez AM, Chase J, Cope EK, Da Silva R, Diener C, Dorrestein PC, Douglas GM, Durall DM, Duvallet C, Edwardson CF, Ernst M, Estaki M, Fouquier J, Gauglitz JM, Gibbons SM, Gibson DL, Gonzalez A, Gorlick K, Guo J, Hillmann B, Holmes S, Holste H, Huttenhower C, Huttley GA, Janssen S, Jarmusch AK, Jiang L, Kaehler BD, Bin Kang K, Keefe CR, Keim P, Kelley ST, Knights D, Koester I, Kosciolek T, et al. 2019. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat Biotechnol 37:852–857.
Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. 2016. DADA2: high-resolution sample inference from Illumina amplicon data. Nat Methods 13:581–583.
Bokulich NA, Kaehler BD, Rideout JR, Dillon M, Bolyen E, Knight R, Huttley GA, Gregory Caporaso J. 2018. Optimizing taxonomic classification of marker-gene amplicon sequences with QIIME 2’s q2-feature-classifier plugin. Microbiome 6:90.
Janssen S, McDonald D, Gonzalez A, Navas-Molina JA, Jiang L, Xu ZZ, Winker K, Kado DM, Orwoll E, Manary M, Mirarab S, Knight R. 2018. Phylogenetic placement of exact amplicon sequences improves associations with clinical information. mSystems 3:e00021-18.
Stamatakis A. 2014. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30:1312–1313.
Zhang H, Yohe T, Huang L, Entwistle S, Wu P, Yang Z, Busk PK, Xu Y, Yin Y. 2018. dbCAN2: a meta server for automated carbohydrate-active enzyme annotation. Nucleic Acids Res 46:W95–W101.
Lombard V, Golaconda Ramulu H, Drula E, Coutinho PM, Henrissat B. 2014. The carbohydrate-active enzymes database (CAZy) in 2013. Nucleic Acids Res 42:D490–D495.
Andrews S. 2010. FastQC: a quality control tool for high throughput sequence data. https://www.bioinformatics.babraham.ac.uk/projects/fastqc/.
Bolger AM, Lohse M, Usadel B. 2014. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30:2114–2120.
Buchfink B, Xie C, Huson DH. 2015. Fast and sensitive protein alignment using DIAMOND. Nat Methods 12:59–60.
Anderson MJ. 2001. A new method for non-parametric multivariate analysis of variance. Austral Ecol 26:32–46.
Anderson MJ. 2017. Permutational multivariate analysis of variance (PERMANOVA), p 1–15. In Wiley StatsRef: statistics reference online. John Wiley & Sons, Ltd., Hoboken, NJ.
Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, Minchin PR, O’Hara RB, Simpson GL, Solymos P, Stevens MHH, Wagner H, Szoecs E. 2019. vegan: community ecology package. R.

Information & Contributors


Published In

cover image Applied and Environmental Microbiology
Applied and Environmental Microbiology
Volume 88Number 622 March 2022
eLocator: e02429-21
Editor: Jennifer B. Glass, Georgia Institute of Technology
PubMed: 35108096


Received: 9 December 2021
Accepted: 23 January 2022
Accepted manuscript posted online: 2 February 2022
Published online: 22 March 2022


  1. bacteria
  2. global change
  3. microdiversity



Department of Ecology and Evolutionary Biology, University of California, Irvine, California, USA
Center for Marine Biotechnology and Biomedicine, Scripps Institution of Oceanography, University of California San Diego, La Jolla, California, USA
S. S. Finks
Department of Ecology and Evolutionary Biology, University of California, Irvine, California, USA
A. A. Malik
School of Biological Sciences, University of Aberdeen, Aberdeen, United Kingdom
C. Weihe
Department of Ecology and Evolutionary Biology, University of California, Irvine, California, USA
S. D. Allison
Department of Ecology and Evolutionary Biology, University of California, Irvine, California, USA
Department of Earth System Science, University of California, Irvine, California, USA
Department of Ecology and Evolutionary Biology, University of California, Irvine, California, USA
Department of Earth System Science, University of California, Irvine, California, USA
Department of Ecology and Evolutionary Biology, University of California, Irvine, California, USA


Jennifer B. Glass
Georgia Institute of Technology


The authors declare no conflict of interest.

Metrics & Citations


Note: There is a 3- to 4-day delay in article usage, so article usage will not appear immediately after publication.

Citation counts come from the Crossref Cited by service.


If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. For an editable text file, please select Medlars format which will download as a .txt file. Simply select your manager software from the list below and click Download.

View Options

Figures and Media






Share the article link

Share with email

Email a colleague

Share on social media

American Society for Microbiology ("ASM") is committed to maintaining your confidence and trust with respect to the information we collect from you on websites owned and operated by ASM ("ASM Web Sites") and other sources. This Privacy Policy sets forth the information we collect about you, how we use this information and the choices you have about how we use such information.
FIND OUT MORE about the privacy policy