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 transformationpan> and stabilizationpan>, nutrient and biogeochemical cyclinpan>g and plant host defence. The compositionpan> of these communpan>ities is sensitive to a variety of global changes, and such shifts can alter their funpan>ctionpan>inpan>g [1,2]. A predictive unpan>derstandinpan>g of how these communpan>ities responpan>d to their environpan>ment is therefore of great inpan>terest. Amonpan>g a variety of obstacles, the enormous diversity of soil bacteria creates a challenge for makinpan>g predictionpan>s.

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 addition at 13 locations across five continents [14]. The bacterial response 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 n class="Chemical">CO2 conpan>centrationpan>, P additionpan> and inpan>creased soil pH. This wider collectionpan> of types of experiments allowed us to test three new hypotheses. First, we hypothesized that the responpan>ses to each type of perturbationpan> would vary inpan> 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 unpan>derlyinpan>g the responpan>se, which do vary [3]. For inpan>stance, previous studies suggest that N uptake rate is more deeply conpan>served than the preference for soil temperature [3,18,19]. However, a systematic comparisonpan> of the phylogenetic depth of responpan>ses across many perturbationpan> types has not been donpan>e.

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 addition might directly select for bacteria with traits for enhanpan>ced P acquisition [24,25]. By contrast, the primary effects of elevated pan> class="Chemical">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 concentration, P addition anpan>d liminpan>g (inpan>creasinpan>g pH). Each study met the followinpan>g criteria: (i) published before 1 June 2019, (ii) inpan>cluded at least three replicates for manpan>ipulation anpan>d control (non-addition) plots, (iii) used high-throughput amplicon sequencinpan>g containpan>inpan>g the V4, V3–V4 or V4–V5 region of 16S rRNA genes inpan> bacteria, anpan>d (iv) sampled from surface soil (withinpan> top 15 cm). We identified 27 published studies but excluded 14 because the raw sequence datasets anpan>d/or accompanpan>yinpan>g metadata were not publicly deposited or otherwise obtainpan>able ((table 1; electronic supplementary material, table S1). Ultimately, we inpan>cluded data from six locations for warminpan>g, seven locations for drought, six locations for elevatinpan>g atmospheric pan> class="Chemical">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 genes that corresponds to the region 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 permutations 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

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 n class="Chemical">phyloseq package [51] inpan> the R environpan>ment. We used DESeq2 inpan> a limited way that differs from its typical use for RNA-Seq data. Specifically, we used it to (i) normalize the sequence counpan>ts by sample withinpan> a locationpan> by replacinpan>g the originpan>al counpan>ts with variance stabilized counpan>ts and (ii) calculate the log2-fold ratio of averaged relative abunpan>dance inpan> manipulationpan> plots relative to conpan>trol plots for each taxonpan>. Before calculatinpan>g the responpan>se ratios, we first removed rare OTUs present inpan> less than half of all the plots withinpan> a locationpan> (electronpan>ic supplementary material, table S1), as the responpan>se ratios of low occupancy, low abunpan>dance taxa would be subject to a high degree of noise. For seven locationpan>s where this filterinpan>g process left very few OTUs (fewer than 900), we relaxed this criterionpan> to use OTUs present inpan> two or more plots (electronpan>ic supplementary material, table S1). Note that we did not use DESeq2 to test for statistical signpan>ificance (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 (NJ) phylogenetic tree was inferred with bootstrap analysis (100 replicates) using the phangorn package [53]. We then applied a consenTRAIT analysis [4] (using the castor package [54]) to test whether an OTU's response to perturbation was related to the bacterial phylogeny.

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 NJ phylogenetic tree based on 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 descending tips (OTUs) shared the same directional response (a ‘consensus’ clade). The genetic depth (the average distance of the node to its descending tips) and size (total number of the descending tips) of each consensus clade was calculated. The genetic depth of clades with a single descending tip (OTU) was calculated as half the branch length to the nearest neighbour as previously recommended [4]. Finally, the mean genetic depth, τD, of the consensus clades sharing either positive or negative responses was calculated (electronic supplementary material, figure S2). To assess the statistical significance of phylogenetic conservation of responses, simulated τD values were calculated by randomizing the responses among the tips 1000 times. The probability of phylogenetic conservation (non-randomness) of the distribution of positive and negative responses was calculated as the fraction 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 n class="Chemical">CO2) and 1079 (P additionpan>) OTUs that were present inpan> three or more experimental locationpan>s. (Note these OTUs were from the pool of nonpan>-rare OTUs inpan> each locationpan>; electronpan>ic supplementary material, table S1.) For each of these widespread OTUs, we averaged the responpan>se values across locationpan>s. This procedure treats the results of each locationpan> equally, regardless of differences inpan> methods and sequencinpan>g effort. Unlike a typical diversity metric (e.g. a metric of richnpan>ess or phylogenetic diversity), the responpan>se ratio parameter that we estimate for each taxonpan> inpan> a locationpan> should not be biased by sequencinpan>g effort, although it will presumably get more accurate with more sequencinpan>g. For each perturbationpan> experiment, we then created a NJ tree of these widespread OTUs and performed the conpan>senTRAIT 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). Notably, this version of SILVA uses the Genome Taxonomy Database (GTDB), which classifies taxonomy based on monophyletic lineages and normalizes taxonomic ranks based on phylogenetic depth [56]. We then calculated the number of OTUs that had a positive or negative response at each phylum, class, order, family or genus level. We performed a two-tailed exact test [57] against the equal distribution of positive and negative responses within each taxonomic group. To compare these results with our previous study about bacterial responses to 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 (PERMANOVA; 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 n class="Chemical">CO2 anpan>d P addition, respectively. The response to the sinpan>gle soil liminpan>g experiment was also phylogenetically patternpan>ed (table 2). As we found previously [14], the consenTRAIT results were robust to the phylogenetic reconstruction method (neighbour-joinpan>inpan>g or maximum-likelihood), as expected from the high correlation between trees for the different locations anpan>d perturbation types (electronic supplementary material, table S2). Notably, the phylogenetic anpan>alysis was more sensitive thanpan> the community anpan>alysis; the locations that displayed overall shifts inpan> bacterial composition were a subset of those that showed significanpan>t phylogenetic conservation. This discrepanpan>cy likely arises because the phylogenetic anpan>alysis gives equal weight to all taxa, whereas the community anpan>alysis primarily considers the most abundanpan>t 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 gene amplicon of 3.4–4.8% among OTUs. These values differed among the perturbation types (one-way ANOVA; p = 0.05, figure 2). This difference was driven by the deeper level of conservation of the elevated n class="Chemical">CO2 response thanpan> 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 anpan>d warminpan>g experiments substanpan>tially overlapped. The soil liminpan>g experiment showed the deepest response of anpan>y 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.)

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 rDNA sequence) than the mean of the individual locations (diamond 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. present at three locationpan>s or more). Coloured nodes and linpan>eages show the conpan>sensus clades inpan> which more than 90% of the descendant OTUs show the same responpan>se directionpan> (blue, positive responpan>se; red, negative responpan>se). Note that all OTUs inpan> the tree responpan>d, but onpan>ly conpan>sensus clades are coloured for clarity. The outer rinpan>g shows the phylum-level taxonpan>omy of OTUs determinpan>ed usinpan>g the RDP classifier trainpan>ed onpan> 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 addition. Similarly, the family Sphingobacteriaceae (phylum Bacteroidetes) responded in opposing directions 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 n class="Chemical">CO2 anpan>d 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 addition. Although one might hypothesize that some shared traits might explain the response to nutrient addition generally, only one family within the order Oligoflexales responded consistently (negatively) to both treatments. This result is consistent with previous studies that found that different taxa (OTUs) from the same community responded to different perturbations (e.g. heat n class="Disease">shock versus cold shock [28] anpan>d warminpan>g versus drought [29]). Our study suggests that this trend might hold true for broader taxonomic groups across a ranpan>ge of locations anpan>d 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 gene amplicon, or approximately the level of a bacterial genus [63]. The one exception was a study that manipulated pH by liming (τD = 0.022), and pH preference appears to be quite deeply conserved relative to other traits [3]. All of these values are relatively shallow compared with traits such as methanogenesis (τD = 0.071), but much more conserved than organic phosphorus acquisition or pan> class="Chemical">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 clone library anpan>alysis [72], responded positively to elevated pan> class="Chemical">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 anpan>d P addition experiments were conducted inpan> grasslanpan>ds. Thus, the degree to which the bacterial responses are context dependent might inpan>crease with the inpan>clusion of additional experiments from a broader ranpan>ge of soil biomes. By contrast, a wider representation of soil communities would inpan>crease the degree of phylogenetic representation across bacterial phylogeny anpan>d improve our ability to predict the responses of particular taxa or clades to particular perturbations.

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 additional evidence that efforts to compile databases of microbial traits [74], classify responses to within experiments [2] and develop new statistical methods [75] make this a tractable goal. A critical next step, however, is to connect shifts in community composition to the functional processes soil bacteria carry out. Many potential functions (as predicted through genomic traits) are phylogenetically conserved [4,74], and new techniques estimating taxon-specific process rates find that these measures are also phylogenetically patterned [10,19]. Finally, recent evidence suggests that phylogeny has a stronger effect than the environment on bacterial processes such as growth rate and n class="Chemical">carbon assimilation [76]. This result suggests that, like their responses inpan> abundanpan>ce, bacterial process rates might not be overly context dependent. In this case, phylogenetic patternpan>s of bacterial responses anpan>d their functional traits could be combinpan>ed to predict how global chanpan>ges will alter ecosystem functioninpan>g, as has been proposed for both microbial anpan>d planpan>t communities [62,77].