Phylogenetic conservation of soil bacterial responses to simulated global changes.

Soil bacterial communities are altered by anthropogenic drivers such as climate change-related warming and fertilization. However, we lack a predictive understanding of how bacterial communities respond to such global changes. Here, we tested whether phylogenetic information might be more predictive of the response of bacterial taxa to some forms of global change than others. We analysed the composition of soil bacterial communities from perturbation experiments that simulated warming, drought, elevated CO2 concentration and phosphorus (P) addition. Bacterial responses were phylogenetically conserved to all perturbations. The phylogenetic depth of these responses varied minimally among the types of perturbations and was similar when merging data across locations, implying that the context of particular locations did not affect the phylogenetic pattern of response. We further identified taxonomic groups that responded consistently to each type of perturbation. These patterns revealed that, at the level of family and above, most groups responded consistently to only one or two types of perturbations, suggesting that traits with different patterns of phylogenetic conservation underlie the responses to different perturbations. We conclude that a phylogenetic approach may be useful in predicting how soil bacterial communities respond to a variety of global changes. This article is part of the theme issue 'Conceptual challenges in microbial community ecology'.

Introduction

Soil bacterial communities play critical roles in ecosystem functioning such as n class="Chemical">carbon transformation and stabilization, nutrient and biogeochemical cycling and plant host defence. The composition of these communities is sensitive to a variety of global changes, and such shifts can alter their functioning [1,2]. A predictive understanding of how these communities respond to their environment is therefore of great interest. Among a variety of obstacles, the enormous diversity of soil bacteria creates a challenge for making predictions.

Phylogenetic information might simplify this diversity by offering a structure to its underlying biological variation or traits [3]. In particular, if large phylogenetic clades of bacterial taxa (such as amplicon sequence variants (ASVs) or 97% operational taxonomic units (OTUs)) respond in a similar manner, then one could reduce the number of bacterial groups considered. In such a case, the response would be said to be phylogenetically conserved, defined here as positive or negative responses that are non-randomly distributed across the bacterial phylogenic tree [4]. The more deeply conserved the response (the greater genetic depth at which descendant bacterial taxa show a similar response), the fewer groups that one would need to track as they could be lumped into broader taxa (figure 1a).

Figure 1.

Conceptual framework for the study's three hypotheses. Bacterial taxa respond either positively (blue) or negatively (red) in their relative abundance to a perturbation. (a) The responses to a perturbation might be phylogenetically conserved (left) or random (right). (b) Hypothesis 1: The depth of conservation of the responses (the average phylogenetic depth of clades responding positively or negatively) varies by type of perturbation. Here, the response to perturbation A is more deeply conserved than the response to perturbations B or C. (c) Hypothesis 2: The degree to which a bacterial response is context dependent might depend on the type of perturbation. The degree of context dependency is assessed by comparing the difference between the average depth of responses at individual locations (solid circles) and the depth of responses when the datasets are merged across locations (filled green diamonds). In this hypothetical example, the context dependence of perturbation A is smaller than that of perturbations B and C. (d) Hypothesis 3: The consistency of a clade's response to different perturbations. Clades might respond consistently (responding generally positively or negatively) to multiple perturbations, consistently to just one perturbation, or inconsistently to all.

More generally, such a phylogenetic signal demonstrates a tendency for more closely related taxa to be more similar in their traits than less closely related taxa [5]. Indeed, bacterial traits are consistently phylogenetically patterned [3,6,7], despite that horizontal gene transfer and rapid evolution may break up any pattern [8,9]. Specifically, several recent studies have found that bacterial responses to environmental changes are also phylogenetically conserved [10-13] such that closely related taxa respond more similarly to a perturbation than those that are distantly related. We recently conducted an analysis of soil bacterial response to experimental N additionpan> at 13 locationpan>s across five conpan>tinenpan>ts [14]. The bacterial responpan>se to added N at individual locations was phylogenetically conserved across the tree of life; closely related bacterial taxa on the tree responded more similarly to N addition than distantly related taxa. Further, we found that the phylogenetic pattern of responses was context independent [14]. Phylogenetic clades generally responded in the same direction (positively or negatively) in different locations, even as the baseline environment or surrounding microbial community varied. As a result, we could identify the taxonomy of larger phylogenetic clades that responded to N addition in the same direction in a variety of experimental locations. Yet increased N availability is just one way that soil ecosystems may change in the future [15-17]. We therefore wanted to expand on this prior work to consider whether a phylogenetic approach might be more predictive of bacterial responses for some global changes than others.

To do this, here we re-analysed publicly available data on soil bacterial community composition from field experiments that simulated a variety of global changes, including warming, drought, elevated atmospheric CO2 conpan>cenpan>trationpan>, P additionpan> and increased soil pH. This wider collectionpan> of types of experimenpan>ts allowed us to test three new hypotheses. First, we hypothesized that the responpan>ses to each type of perturbationpan> would vary in the depth at which they are conpan>served (figure 1b). While a responpan>se is not a trait itself, we reasonpan>ed that the degree of conpan>servationpan> is based onpan> the bacterial traits underlying the responpan>se, which do vary [3]. For instance, previous studies suggest that N uptake rate is more deeply conserved than the preference for soil temperature [3,18,19]. However, a systematic comparison of the phylogenetic depth of responses across many perturbation types has not been done.

Second, we hypothesized that the degree to which a bacterial response is context dependent depends on the type of perturbation (figure 1c). We tested this hypothesis by comparing the average depth of responses at individual locations to the depth of responses when the datasets for each treatment type were merged (as in [14]). If the bacterial responses were highly context dependent, then the averaged phylogenetic depth in the merged dataset would be shallower than the individual locations. We reasoned that the importance of the surrounding context might depend on the degree to which a perturbation directly impacted bacterial physiology and growth. For instance, increased drought might directly select for bacteria with drought-adaptation traits such as the ability to accumulate osmolytes, produce exopolysaccharides or form thick cell walls [20-23]. Similarly, P additionpan> might directly select for bacteria with traits for enpan>hanced P acquisitionpan> [24,25]. By conpan>trast, the primary effects of elevated CO2 and warming on soil bacteria might be mediated through the response of the surrounding plant community through changes in litter or root exudate chemistry [26,27]. Thus, the bacterial traits selected for by the perturbation may depend on the baseline nutrient environment at each location. In addition, we expected that variation in the experimental manipulations across studies might increase the context dependency of response. For instance, the warming experiments increased temperature between 1.5 and 5°C for anywhere between 1 year and 20 years; by contrast, the P addition experiments all received the same amount and type of P for 2–4 years.

Finally, we hypothesized that, while a clade's response to different perturbations would not generally be correlated [28,29], some clades might be consistently sensitive to environmental changes and therefore several types of perturbations (figure 1d). This pattern might occur if a clade is adapted to responding positively to environment change––for instance, possessing traits that allow a bacterium to turn on reproduction quickly, persist under stressful conditions or use a wide range of substrates [20,30]. Alternatively, a clade might respond consistently negatively if members are highly specialized to their preferred environment, such that any change in their conditions results in a decline in abundance.

Material and Methods

Study inclusion criteria

We searched for published studies that assessed soil bacterial community composition within global change field experiments. The experiments manipulated the soil environment by warming, rain exclusion (drought), elevating atmospheric CO2 conpan>cenpan>trationpan>, P additionpan> and liming (increasing pH). Each study met the following criteria: (i) published before 1 June 2019, (ii) included at least three replicates for manipulationpan> and conpan>trol (nonpan>-additionpan>) plots, (iii) used high-throughput ampliconpan> sequenpan>cing conpan>taining the V4, V3–V4 or V4–V5 regionpan> of 16S rRNA genes in bacteria, and (iv) sampled from surface soil (within top 15 cm). We identified 27 published studies but excluded 14 because the raw sequence datasets and/or accompanying metadata were not publicly deposited or otherwise obtainable ((table 1; electronic supplementary material, table S1). Ultimately, we included data from six locations for warming, seven locations for drought, six locations for elevating atmospheric CO2 concentration and six locations for P addition (table 1; electronic supplementary material, figure S1). Although we had only one location for liming, we included the data because soil acidification is also a global problem [42] and bacterial pH preference is thought to be a particularly deeply conserved trait [3]. The specific samples used at each location are summarized in electronic supplementary material, table S1.

Table 1.

Characteristics of study locations and number of replicate plots and OTUs used in this study (control, treatment).

perturbation	location	habitat (country)	treatment	duration	no. replicate plotsa	no. OTUs
warming	Che [31]	grassland (China)	plus 1.57°C	6 years	4, 4	2046
	DeAngelis_1 [32]	temperate forest (USA)	plus 5°C	5 years	4, 4	1138
	DeAngelis_2 [32]	temperate forest (USA)	plus 5°C	8 years	4, 4	1062
	DeAngelis_3 [32]	temperate forest (USA)	plus 5°C	20 years	4, 4	1046
	Waghmode [33]	cropland (China)	plus 1.5°C	7 years	6, 6	2298
	Zhang [29]	grassland (China)	plus 2°C	1 year	6, 6	1372
drought	Bastida_1 [34]	forest (Spain)	rainfall exclusion	6 years	6, 6	1631
	Bastida_2 [34]	forest (Spain)	rainfall exclusion	6 years	6, 6	1698
	Bouskill_1 [35]	forest (Puerto Rico)	rainfall exclusion	4 years	3, 3	934
	Bouskill_2	desert (USA)	rainfall exclusion	3–10 months	10, 20	1730
	Fernandes_1 [36]	desert soil crust (USA)	rainfall exclusion	3 years	10, 10	1903
	Fernandes_2 [36]	desert soil crust (USA)	rainfall exclusion	3 years	10, 10	2453
	Zhang [29]	grassland (China)	rainfall exclusion	1 year	6, 6	1612
elevated CO2	Deng [37]	grassland (USA)	368 ppm, 560 ppma	10 years	12, 12	1058
	Raut_1 [38]	grassland (USA)	380–250 ppm, 500–380 ppma	9 years	11, 16	1082
	Raut_2 [38]	grassland (USA)	380–250 ppm, 500–380 ppma	9 years	18, 14	1051
	Raut_3 [38]	grassland (USA)	380–250 ppm, 500–380 ppma	9 years	12, 12	956
	Xia [39]	grassland (New Zealand)	ambient, 475 ppma	12 years	3, 3	708
	Yang [40]	grassland (USA)	ambient, ambient + 275 ppma	4 years	4, 4	1447
P addition	Leff_1 [24]	grassland (Switzerland)	100 kg P ha−1 yr−1 as Ca(H2PO4)2	3 years	6, 3	2069
	Leff_2 [24]	grassland (Australia)	100 kg P ha−1 yr−1 as Ca(H2PO4)2	3 years	6, 3	1307
	Leff_3 [24]	grassland (South Africa)	100 kg P ha−1 yr−1 as Ca(H2PO4)2	3 years	6, 3	1232
	Leff_4 [24]	grassland (South Africa)	100 kg P ha−1 yr−1 as Ca(H2PO4)2	2 years	3, 3	1530
	Leff_5 [24]	grassland (South Africa)	100 kg P ha−1 yr−1 as Ca(H2PO4)2	2 years	6, 3	1545
	Leff_6 [24]	grassland (Australia)	100 kg P ha−1 yr−1 as Ca(H2PO4)2	4 years	6, 3	1275
liming	Guo [41]	paddy (China)	7500 kg ha−1 of CaCO3	>1 year	4, 4	945

acontrol plots, treatment plots.

Sequence processing

Sequence data (FASTQ-formatted raw sequence or FASTA-formatted denoised sequence files) and associated metadata were either shared by the original authors or downloaded from public databases. All of the datasets were sequenced on the Illumina MiSeq or Roche 454 platforms. In general, results from these platforms are comparable [43], and the sequencing error rates of both platforms are low relative to sequence differences between our OTU classifications (see below).

To allow comparisons across studies, we reanalysed the sequence data in a consistent manner. We trimmed all sequences to the V4 region of 16S rRNA genpan>es that corresponpan>ds to the regionpan> amplified with 515F (GTGYCAGCMGCCGCGGTAA)/806R (GGACTACNVGGGTWTCTAAT) primers, after removing the primers and sequences outside of the target region with the Cutadapt toolkit [44]. The UPARSE pipeline [45] was used to merge the paired-end sequences of FASTQ-formatted raw sequence files, conduct quality filtering and cluster the sequences into OTUs. A minimum overlap of 20 bp was set for merging the sequences. A maximum per sequence expected error frequency value of 1.0 was set for quality filtering the sequences. Singleton sequences were removed. Paired-end sequences from all locations within a perturbation experiment (e.g. for all warming studies or for all drought studies) were merged and clustered into OTUs at ≥ 97% sequence similarity, and chimeric sequences were removed at the step of OTU clustering. Taxonomy of the representative sequence of each OTU was assigned within QIIME using the RDP classifier [46] at 80% confidence threshold trained on the latest version of SILVA database (v. 132, https://www.arb-silva.de/download/archive/qiime/). OTUs assigned as chloroplasts or mitochondria, unassigned at kingdom level, and Archaea were removed. The archaeal OTUs were relatively rare and archaeal 16S rRNA genes are known to be preferentially amplified among the primer pairs used [47].

Using a 97% OTU definition (rather than ASVs as in [48]) allowed us to compare the same taxon across many locations, which was key to our analysis. However, we also tested the sensitivity of our results by reanalysing the warming experiments using ASVs with the DADA2 pipeline (v. 1.12) [48]. The values of τD, the mean phylogenetic depth of clades across a tree of life sharing either positive or negative responses to a particular perturbation (figure 1b and see below), were 0.5–0.6 times lower when using the ASV pipeline; this result is to be expected because the abundance of each ASV is generally lower than each OTU and the response ratios of low abundance taxa would be subject to a high degree of noise. However, the τD values based on 97% OTUs and ASVs for each experimental location were strongly positively correlated for both positive (R2 = 0.88) and negative (R2 = 0.83) responses. This tight correlation suggests that our conclusions are robust to the OTU definition used [2].

Overall community composition responses at each location

We first tested whether the perturbations altered overall community composition at each location. Because we obtained a variable number of sequence reads per sample within each location, the sequence data were rarefied to the lowest number of reads per location to account for the variation of samples within each location (table 2). The Bray–Curtis dissimilarity metric was used to calculate compositional differences between each sample from the rarefied OTU tables with the vegan package [49] in the R environment (v. 3.6.0, http://www.R-project.org). We tested for differences in bacterial community composition between manipulation and control plots with a permutational multivariate analysis of variance (PERMANOVA) test with 999 permutationpan>s using the vegan package. Note that our PERMANOVA results could differ from the original study; such discrepancies could arise because we sometimes used only a subset of samples based on the treatment of interest, as well as differences due to sequence processing.

Table 2.

Permutational multivariate analysis of variance (PERMANOVA) and consenTRAIT results comparing control and treatment plots by each study location. Bold values indicate a significant response (p < 0.05); bold italicized values indicate a highly significant response (p < 0.005). The consenTRAIT statistic (τD) is given for both positive and negative responding consensus clades, defined as clades in which >90% of the descendant OTUs show the same direction of response.

perturbation	location	PERMANOVA		consenTRAIT
		no. of rarefied sequences	R2	τD of positive response	τD of negative response
warming	Che [31]	14 577	0.131	0.020	0.021
	DeAngelis_1 [32]	9370	0.286	0.020	0.017
	DeAngelis_2 [32]	6226	0.206	0.018	0.017
	DeAngelis_3 [32]	8285	0.165	0.017	0.019
	Waghmode [33]	13 408	0.218	0.022	0.021
	Zhang [29]	2275	0.100	0.022	0.020
	merging locations	—	—	0.019	0.018
drought	Bastida_1 [34]	9702	0.195	0.019	0.020
	Bastida_2 [34]	9647	0.130	0.019	0.020
	Bouskill_1 [35]	1225	0.168	0.019	0.018
	Bouskill_2	5933	0.045	0.017	0.020
	Fernandes_1 [36]	35 695	0.335	0.018	0.019
	Fernandes_2 [36]	38 303	0.111	0.018	0.021
	Zhang [29]	2228	0.118	0.023	0.022
	merging locations	—	—	0.017	0.019
elevated CO2	Deng [37]	211	0.043	0.018	0.021
	Raut_1 [38]	2729	0.058	0.019	0.022
	Raut_2 [38]	2168	0.067	0.021	0.022
	Raut_3 [38]	2065	0.098	0.022	0.019
	Xia [39]	2054	0.202	0.020	0.018
	Yang [40]	5790	0.179	0.018	0.019
	merging locations	—	—	0.019	0.019
P addition	Leff_1 [24]	27 667	0.123	0.019	0.019
	Leff_2 [24]	16 881	0.191	0.019	0.019
	Leff_3 [24]	17 984	0.150	0.020	0.018
	Leff_4 [24]	18 855	0.146	0.018	0.017
	Leff_5 [24]	15 168	0.191	0.017	0.018
	Leff_6 [24]	18 855	0.069	0.017	0.020
	merging locations	—	—	0.016	0.019
liming	Guo [41]	3230	0.177	0.024	0.020

Permutational multivariate analysis of variance (PERMAn class="Chemical">NOVA) and consenTRAIT results comparing control and treatment plots by each study location. Bold values indicate a significant response (p < 0.05); bold italicized values indicate a highly significant response (p < 0.005). The consenTRAIT statistic (τD) is given for both positive and negative responding consensus clades, defined as clades in which >90% of the descendant OTUs show the same direction of response.

Individual operational taxonomic unit responses at each location

To quantify the response of bacterial OTUs to a perturbation at each location, we used the full sequence dataset (not rarefied) and the DESeq2 package [50] in the platform of the phyloseq package [51] in the R enpan>vironpan>menpan>t. We used DESeq2 in a limited way that differs from its typical use for RNA-Seq data. Specifically, we used it to (i) normalize the sequence counts by sample within a location by replacing the original counts with variance stabilized counts and (ii) calculate the log2-fold ratio of averaged relative abundance in manipulation plots relative to control plots for each taxon. Before calculating the response ratios, we first removed rare OTUs present in less than half of all the plots within a location (electronic supplementary material, table S1), as the response ratios of low occupancy, low abundance taxa would be subject to a high degree of noise. For seven locations where this filtering process left very few OTUs (fewer than 900), we relaxed this criterion to use OTUs present in two or more plots (electronic supplementary material, table S1). Note that we did not use DESeq2 to test for statistical significance (as the program is often used) but exported the normalized log2-fold ratios for further analyses.

Phylogenetic conservation of responses at each location

To assess whether the response to each perturbation was phylogenetically conserved at an individual location, representative sequences of each OTU (the most abundant sequence within each OTU from the experiments) were aligned using the DECIPHER package [52]. A neighbour-joining (n class="Chemical">NJ) phylogenpan>etic tree was inferred with bootstrap analysis (100 replicates) using the phangornpan> package [53]. We thenpan> applied a conpan>senpan>TRAIT analysis [4] (using the castor package [54]) to test whether an OTU's responpan>se to perturbationpan> was related to the bacterial phylogenpan>y.

The consenTRAIT algorithm identifies phylogenetic clades in which the direction of the response is conserved (consensus clades) and calculates the average depth of those clades from a phylogenetic tree. If the response is significantly conserved, then the average phylogenetic depth of those clades is greater than a distribution of randomly distributed responses. The consenTRAIT approach only considers binary traits (here, whether the response is positive or negative). A positive or negative response was assigned for each OTU on the n class="Chemical">NJ phylogenpan>etic tree based onpan> the log2-fold ratio exported from DESeq2. The tree was traversed from the root to the tips, recording the deepest nodes where more than 90% of the descenpan>ding tips (OTUs) shared the same directionpan>al responpan>se (a ‘conpan>senpan>sus’ clade). The genpan>etic depth (the average distance of the node to its descenpan>ding tips) and size (total number of the descenpan>ding tips) of each conpan>senpan>sus clade was calculated. The genpan>etic depth of clades with a single descenpan>ding tip (OTU) was calculated as half the branch lenpan>gth to the nearest neighbour as previously recommenpan>ded [4]. Finally, the mean genpan>etic depth, τD, of the conpan>senpan>sus clades sharing either positive or negative responpan>ses was calculated (electronpan>ic supplemenpan>tary material, figure S2). To assess the statistical signpan>ificance of phylogenpan>etic conpan>servationpan> of responpan>ses, simulated τD values were calculated by randomizing the responpan>ses amonpan>g the tips 1000 times. The probability of phylogenpan>etic conpan>servationpan> (nonpan>-randomness) of the distributionpan> of positive and negative responpan>ses was calculated as the fractionpan> of simulated τD values that were greater than or equal to the observed τD.

We used NJ trees for the consenTRAIT analysis, because the genetic scale of these trees roughly represents sequence dissimilarity. However, to consider whether our results were robust to the phylogenetic reconstruction method, we tested for a correlation between NJ and maximum-likelihood (ML) trees for each perturbation. A ML tree with 100 bootstrap replications was constructed with RAxML v8.0, using the GTR + Gamma distribution model [55] at the CIPRES science gateway (v. 3.3, http://www.phylo.org/index.php/).

Phylogenetic conservation of responses across locations

To assess whether the responses were context dependent, we identified 1364 (warming), 1557 (drought), 1284 (elevated CO2) and 1079 (P additionpan>) OTUs that were presenpan>t in three or more experimenpan>tal locationpan>s. (Note these OTUs were from the pool of non-rare OTUs in each location; electronic supplementary material, table S1.) For each of these widespread OTUs, we averaged the response values across locations. This procedure treats the results of each location equally, regardless of differences in methods and sequencing effort. Unlike a typical diversity metric (e.g. a metric of richness or phylogenetic diversity), the response ratio parameter that we estimate for each taxon in a location should not be biased by sequencing effort, although it will presumably get more accurate with more sequencing. For each perturbation experiment, we then created a NJ tree of these widespread OTUs and performed the consenTRAIT analysis as above with this merged dataset.

Because the responses of widespread OTUs were significantly phylogenetically conserved, we next identified the taxonomy of clades whose response to perturbation was significantly more positive or negative than expected by chance using the RDP classifications with the SILVA database (version 132). n class="Chemical">Notably, this versionpan> of SILVA uses the Genpan>ome Taxonpan>omy Database (GTDB), which classifies taxonpan>omy based onpan> monpan>ophyletic lineages and normalizes taxonpan>omic ranks based onpan> phylogenpan>etic depth [56]. We thenpan> calculated the number of OTUs that had a positive or negative responpan>se at each phylum, class, order, family or genpan>us level. We performed a two-tailed exact test [57] against the equal distributionpan> of positive and negative responpan>ses within each taxonpan>omic group. To compare these results with our previous study about bacterial responpan>ses to n class="Chemical">N addition where the older SILVA database (v. 128) was used [14], we re-identified the taxonomy of those clades with the new SILVA version.

Results

Overall bacterial community composition

After reanalysing the sequence data in a consistent way across all datasets, we first tested whether the perturbations altered overall community composition. Bacterial community composition had a mixed response to the perturbation types. The perturbations significantly altered bacterial composition (PERMAn class="Chemical">NOVA; p < 0.05) in fewer than 50% of the locations (11 of 26 locations), although if locations of marginal significance are included (p < 0.1), this number increases to 60% (16 of 26 locations, table 2). Where composition significantly shifted, the perturbation explained between 5.8 and 28.6% of compositional variation.

Phylogenetic conservation of bacterial responses at individual locations

In contrast with overall composition, bacterial responses were significantly conserved within the majority of locations across all perturbation types. For instance, the responses to soil warming were conserved at all six locations; the mean genetic depth (τD) of the consensus clades (the clades in which more than 90% of the descendant OTUs show the same response direction) for both or either of the responses (positive and negative) was greater than expected given a randomized distribution of responses (permutation test; p < 0.05, table 2). The responses to the other perturbations were phylogenetically patterned in 6 of 7, 5 of 6 and 3 of 6 locations, for drought, elevated CO2 and P additionpan>, respectively. The responpan>se to the single soil liming experimenpan>t was also phylogenpan>etically patternpan>ed (table 2). As we found previously [14], the conpan>senpan>TRAIT results were robust to the phylogenpan>etic reconpan>structionpan> method (neighbour-joining or maximum-likelihood), as expected from the high correlationpan> betweenpan> trees for the differenpan>t locationpan>s and perturbationpan> types (electronpan>ic supplemenpan>tary material, table S2). Notably, the phylogenetic analysis was more sensitive than the community analysis; the locations that displayed overall shifts in bacterial composition were a subset of those that showed significant phylogenetic conservation. This discrepancy likely arises because the phylogenetic analysis gives equal weight to all taxa, whereas the community analysis primarily considers the most abundant taxa.

While the responses to the perturbations were phylogenetically conserved, there was less evidence that the degree of conservation varied by perturbation type, in contrast with our first hypothesis (figure 1b). The mean depth of the consensus clades responding positively or negatively (τD) ranged from 0.017 to 0.024 across all perturbations and all locations (table 2), equivalent to an average sequence dissimilarity in the 16S rRNA genpan>e ampliconpan> of 3.4–4.8% amonpan>g OTUs. These values differed amonpan>g the perturbationpan> types (onpan>e-way ANOVA; p = 0.05, figure 2). This difference was driven by the deeper level of conservation of the elevated CO2 response than of the N addition response (post hoc Tukey test: p = 0.03), whereas the depth of conservation of responses to N addition, P addition, drought and warming experiments substantially overlapped. The soil liming experiment showed the deepest response of any study location, but as the only experiment of its type, we excluded it from the statistical test.

Figure 2.

Mean genetic depth (τD) of consensus clades as calculated with the consenTRAIT algorithm. Open circles indicate τD at each location. Filled circles indicate the average τD at individual locations. Filled (green) diamonds indicate τD of widespread OTUs present in at least three locations. For comparison, the τD values from N addition experiments from [14] are also plotted. The number of locations where the responses were significantly phylogenetically patterned (positive and/or negative) out of the total locations is shown on the left margin (table 2). Different letters on the right margin indicate significant differences (p < 0.05) in the mean τD between perturbations based on one-way ANOVA with Tukey's honestly significant difference. (Online version in colour.)

Mean genetic depth (τD) of consensus clades as calculated with the consenTRAIT algorithm. Open circles indicate τD at each location. Filled circles indicate the average τD at individual locations. Filled (green) diamonds indicate τD of widespread OTUs present in at least three locations. For comparison, the τD values from N additionpan> experimenpan>ts from [14] are also plotted. The number of locationpan>s where the responpan>ses were signpan>ificantly phylogenpan>etically patternpan>ed (positive and/or negative) out of the total locationpan>s is shown onpan> the left margin (table 2). Differenpan>t letters onpan> the right margin indicate signpan>ificant differenpan>ces (p < 0.05) in the mean τD betweenpan> perturbationpan>s based onpan> onpan>e-way ANOVA with Tukey's honestly significant difference. (Online version in colour.)

Context dependence of the bacterial responses

To investigate the context dependence of the responses, we next quantified the average responses of widespread OTUs across locations. Widespread OTUs (present in at least three locations) accounted for more than 50% of the sequences at all locations except one (Bouskill_2) (electronic supplementary material, figure S3). The average responses of these widespread OTUs were also phylogenetically conserved for all perturbations (excluding liming; permutation test; p < 0.05, table 2 and figure 3), indicating that the context of the particular location and experiment did not overshadow the phylogenetic signal of the responses observed at individual locations. However, contrary to our second hypothesis (figure 1c), the degree of context dependence did not vary by perturbation. For all perturbations, the depth (average τD of positive and negative responses) of the clades conserved across locations was exactly 0.001 smaller (only 0.2% difference in 16S rDn class="Chemical">NA sequenpan>ce) than the mean of the individual locationpan>s (diamonpan>d symbols in figure 2).

Figure 3.

Phylogenetic distribution of the averaged responses to (a) soil warming, (b) drought, (c) CO2 elevation and (d) P addition of widespread OTUs (i.e. present at three locations or more). Coloured nodes and lineages show the consensus clades in which more than 90% of the descendant OTUs show the same response direction (blue, positive response; red, negative response). Note that all OTUs in the tree respond, but only consensus clades are coloured for clarity. The outer ring shows the phylum-level taxonomy of OTUs determined using the RDP classifier trained on the SILVA database.

Phylogenetic distribution of the averaged responses to (a) soil warming, (b) drought, (c) n class="Chemical">CO2 elevationpan> and (d) P additionpan> of widespread OTUs (i.e. presenpan>t at three locationpan>s or more). Coloured nodes and lineages show the conpan>senpan>sus clades in which more than 90% of the descenpan>dant OTUs show the same responpan>se directionpan> (blue, positive responpan>se; red, negative responpan>se). n class="Chemical">Note that all OTUs in the tree respond, but only consensus clades are coloured for clarity. The outer ring shows the phylum-level taxonomy of OTUs determined using the RDP classifier trained on the SILVA database.

Correlations of responses among perturbation types

Finally, while the average genetic depth of consensus clades (τD) was relatively shallow for all perturbations, some deep clades (defined by the GTDB taxonomy) responded consistently to some perturbations. For instance, all major phyla present responded consistently to at least one type of perturbation; the direction of the OTU responses within the phylum was significantly more negative or positive rather than random (figure 4; electronic supplementary material, table S3). In particular, more than 80% of OTUs within the phyla Cyanobacteria, Rokubacteria and WPS-2 responded consistently to a perturbation (figure 4; electronic supplementary material, table S3). As one would expect, these patterns were generally even stronger at the class and order levels; for example, the classes Acidobacteriia, Blastocatellia and Holophagae appear to be driving the overall negative response to drought within the phylum Acidobacteria.

Figure 4.

Taxonomic groupings whose response was positive or negative. Responses are coloured when they display significantly more positive or negative responses than expected by chance (two-tailed exact test; p < 0.05). The per cent of OTUs responding positively (blue) or negatively (red) is plotted in the heatmap to the left. The number of OTUs are plotted in the heatmap to the right.

Contrary to our third hypothesis (figure 1d), however, there was little evidence that a clade that responds positively or negatively to one type of perturbation will also respond to other perturbations. In fact, there were relatively few instances where the same taxonomic group responded consistently to more than one perturbation type (figure 4). This result did not seem to be driven by differences in the taxonomic diversity present among the experiments; OTU richness within the taxonomic groups was largely similar across perturbation types (figure 4). For example, many OTUs in the Planctomycetes and Proteobacteria were present in more than one perturbation, but most of the taxonomic groups within these phyla showed consistent responses to only one type of perturbation.

Among the taxonomic groups that responded significantly to two or more perturbations, the responses were often in both directions, showing a positive response to some perturbations and a negative response to others. For instance, the orders Propionibacteriales and Solirubrobacterales within the Actinobacteria responded negatively to warming and positively to N additionpan>. Similarly, the family Sphingobacteriaceae (phylum Bacteroidetes) responpan>ded in opposing directionpan>s to drought (negative) and N addition (positive). The only group that seemed to defy this trend was the phylum Acidobacteria. As a whole, the phylum responded negatively to drought, elevated CO2 and N addition.

Discussion

Our re-analysis of field experiment data provides several indications that phylogenetic information can be used to predict the response of soil bacteria to global changes. Overall, bacterial responses were phylogenetically conserved within the majority of locations across all perturbation types. More specifically, the depth of conservation of these responses did not vary much by perturbation type, inconsistent with our first hypothesis. Further, the degree of context dependence did not depend on the type of perturbation, countering our second hypothesis, but indicating that the results from a variety of perturbation experiments in one location can help predict bacterial community responses in other locations.

Contrary to our third hypothesis, however, there was little evidence of clades that were generally sensitive to perturbations. Indeed, at a broad taxonomic level, most groups responded consistently to only one or two types of perturbations, suggesting that responses to different perturbations are generally due to different traits rather than common traits that broadly increase or decrease a clade's sensitivity to environmental change. A notable example is the response to N and P additionpan>. Although onpan>e might hypothesize that some shared traits might explain the responpan>se to nutrienpan>t additionpan> genpan>erally, onpan>ly onpan>e family within the order Oligoflexales responpan>ded conpan>sistenpan>tly (negatively) to both treatmenpan>ts. This result is conpan>sistenpan>t with previous studies that found that differenpan>t taxa (OTUs) from the same community responpan>ded to differenpan>t perturbationpan>s (e.g. heat shock versus cold shock [28] and warming versus drought [29]). Our study suggests that this trend might hold true for broader taxonomic groups across a range of locations and perturbations.

However, the lack of a correlation between a taxonomic group's sensitivity to one perturbation versus its response to others might ultimately make predictions easier. Further, the result that few taxonomic groups showed consistent responses to more than one perturbation might mean that the interactive effects of multiple global changes on microbial composition might be minimal. On the other hand, the evidence of legacy effects on microbial composition, where the community response to one perturbation influences the response to future perturbations [11,28,35,58], perhaps suggests the opposite. Additional research on the phylogenetic patterns of responses to co-occurring perturbations is therefore needed, as future environmental changes will not occur in isolation [59].

The phylogenetic patterns described above provide several insights into the traits underlying the responses of soil bacteria in the face of a perturbation [60-62]). First, the depth of phylogenetic conservation of all responses fell within a narrow range. The average depth across perturbations (τD = 0.018–0.020) corresponds to a 3.6–4.0% divergence in the 16S rRNA genpan>e ampliconpan>, or approximately the level of a bacterial genpan>us [63]. The onpan>e exceptionpan> was a study that manipulated pH by liming (τD = 0.022), and pH preferenpan>ce appears to be quite deeply conpan>served relative to other traits [3]. All of these values are relatively shallow compared with traits such as methanogenpan>esis (τD = 0.071), but much more conpan>served than organic phosphorus acquisition or carbon substrate utilization, for example. We expect that the response to a perturbation would be governed by a suite of traits that vary in conservation depth. For instance, the response to drought might involve a deeply conserved trait such as spore formation and cell wall type [18] as well as shallower traits such as salt tolerance or biosynthesis of organic osmotic solutes [64].

The patterns in phylogenetic responses also suggest that there is not a common suite of traits involved in responding to most perturbations, as indicated by the lack of correlation in responses among the types of experiments. Instead, the underlying traits responsible for a conserved response seem to be specific to each perturbation type, even though the degree of conservation among the bacterial responses was similar. It is also important to note that the responses to some perturbations might not be due to direct effects of changes in abiotic conditions, but to indirect effects of changes in plant communities and soil and litter chemistry. For instance, warming might increase nutrient availability through the enhanced supply of plant root exudates or decrease nutrient availability through the drying of soil and suppression of litter decomposition [65,66]. Such indirect effects might also explain why it might be easier to observe phylogenetic patterns of bacterial responses in these experimental studies versus along biogeographic gradients [67].

Finally, the phylogenetic responses provide clues about the traits of rare and understudied taxa. For instance, the uncultivated candidate phylum WPS-2, first observed in polluted soil via clone library analysis [68], responded negatively to warming. WPS-2 has a global distribution but most samples containing WPS-2 with high abundance were collected typically from cold environments [69], such as Antarctica soil [70] and Greenland ice [71]. The uncultivated candidate phylum Rokubacteria, first observed in an alpine meadow via clonpan>e library analysis [72], responpan>ded positively to elevated CO2. Metagenomic assemblies suggest that the order Rokubacteriales has the potential for a versatile, mixotrophic metabolism [73]. Elevated CO2 might increase the release of root exudates to soil, enhance nutrient availability and shift the ecological strategy of the soil bacterial community to a higher contribution of fast-growing r-selected taxa [27].

Of course, there are notable limitations to these results, and these caveats point to future areas of research. Perhaps the most importantly, we could only obtain sequence data from a small number of experiments for each type of perturbation. These experimental locations represent a limited number of ecosystems and geographic regions. Indeed, some experiments were located quite close to one another, and all of the elevated n class="Chemical">CO2 and P additionpan> experimenpan>ts were conpan>ducted in grasslands. Thus, the degree to which the bacterial responpan>ses are conpan>text depenpan>denpan>t might increase with the inclusionpan> of additionpan>al experimenpan>ts from a broader range of soil biomes. By conpan>trast, a wider represenpan>tationpan> of soil communities would increase the degree of phylogenpan>etic represenpan>tationpan> across bacterial phylogenpan>y and improve our ability to predict the responpan>ses of particular taxa or clades to particular perturbationpan>s.

Conclusion

The enormous diversity of soil bacteria would seem to overwhelm attempts to make detailed predictions about their responses to global change. Nevertheless, this study provides additionpan>al evidenpan>ce that efforts to compile databases of microbial traits [74], classify responpan>ses to within experimenpan>ts [2] and develop new statistical methods [75] make this a tractable goal. A critical next step, however, is to conpan>nect shifts in community compositionpan> to the functionpan>al processes soil bacteria carry out. Many potenpan>tial functionpan>s (as predicted through genpan>omic traits) are phylogenpan>etically conpan>served [4,74], and new techniques estimating taxonpan>-specific process rates find that these measures are also phylogenpan>etically patternpan>ed [10,19]. Finally, recenpan>t evidenpan>ce suggests that phylogenpan>y has a stronpan>ger effect than the enpan>vironpan>menpan>t onpan> bacterial processes such as growth rate and carbon assimilation [76]. This result suggests that, like their responses in abundance, bacterial process rates might not be overly context dependent. In this case, phylogenetic patterns of bacterial responses and their functional traits could be combined to predict how global changes will alter ecosystem functioning, as has been proposed for both microbial and plant communities [62,77].