Global projections of the soil microbiome in the Anthropocene.

AIM: Soil microbes are essential for maintenance of life-supporting ecosystem services, but projections of how these microbes will be affected by global change scenarios are lacking. Therefore, our aim was to provide projections of future soil microbial distribution using several scenarios of global change. LOCATION: Global. TIME PERIOD: 1950-2090. MAJOR TAXA STUDIED: Bacteria and fungi.
METHODS: We used a global database of soil microbial communities across six continents to estimate past and future trends of the soil microbiome. To do so, we used structural equation models to include the direct and indirect effects of changes in climate and land use in our predictions, using current climate (temperature and precipitation) and land-use projections between 1950 and 2090.
RESULTS: Local bacterial richness will increase in all scenarios of change in climate and land use considered, although this increase will be followed by a generalized community homogenization process affecting > 85% of terrestrial ecosystems. Changes in the relative abundance of functional genes associated with the increases in bacterial richness are also expected. Based on an ecological cluster analysis, our results suggest that phylotypes such as Geodermatophilus spp. (typical desert bacteria), Mycobacterium sp. (which are known to include important human pathogens), Streptomyces mirabilis (major producers of antibiotic resistance genes) or potential fungal soil-borne plant pathogens belonging to Ascomycota fungi (Venturia spp., Devriesia spp.) will become more abundant in their communities. MAIN
CONCLUSIONS: Our results provide evidence that climate change has a stronger influence on soil microbial communities than change in land use (often including deforestation and agricultural expansion), although most of the effects of climate are indirect, through other environmental variables (e.g., changes in soil pH). The same was found for microbial functions such as the prevalence of phosphate transport genes. We provide reliable predictions about the changes in the global distribution of microbial communities, showing an increase in alpha diversity and a homogenization of soil microbial communities in the Anthropocene.

Introduction

Global assessments continue to provide strong evidence that humans are causing unprecedented loss of biodiversity (Díaz et al., 2019). However, existing information is strongly biased towards a few groups of vertebrates and plants (Díaz et al., 2019; Eisenhauer et al., 2019), and much less is known about the potential changes in belowground communities under global change scenarios. With hundreds of thousands of taxa per gram of soil, bacteria and fungi are the most diverse groups of soil-dwelling organisms across the globe (Delgado-Baquerizo & Eldridge, 2019; Maron et al., 2018; Tedersoo et al., 2014). Soil microbial communities are largely an unseen majority, which conducts a wide range of ecosystem functions (Delgado-Baquerizo & Eldridge, 2019; Maron et al., 2018; Tedersoo et al., 2014) (e.g., decomposition, pathogenesis, nutrient cycling; Fierer et al., 2012; Schwarz et al., 2017; Stürmer et al., 2018), with implications for human well-being (e.g., human health, food production; Wall et al., 2015) and ecosystem sustainability (Bahram et al., 2018; Tedersoo et al., 2014). Recent studies have provided crucial advances in the ecological preferences and biogeographical patterns of soil microbial communities globally (Bahram et al., 2018; Delgado-Baquerizo & Eldridge, 2019; Delgado-Baquerizo, Oliverio, et al., 2018; Ramirez et al., 2019). Such efforts have culminated in the first global (Bahram et al., 2018; Delgado-Baquerizo, Oliverio, et al., 2018; Tedersoo et al., 2014; van den Hoogen et al., 2019) and national (Karimi, Prévost-Bouré, et al., 2018) atlases for soil communities, whereas atlases for plants and other organisms have been available for decades (Tedersoo et al., 2014). Unlike the situation for plants, birds or mammals (Di Marco et al., 2019; Newbold et al., 2015; Powers & Jetz, 2019), only recently have we started to understand how drivers of global change affect soil microbial communities (Delgado-Baquerizo & Eldridge, 2019), which is a crucial step towards forecasting their fate accurately in an increasingly human-dominated world.

Biodiversity projections have advanced our knowledge of the potential changes in biodiversity under scenarios of future global change and enable decision-makers to adjust policies in line with potential future outcomes (Di Marco et al., 2019; Dornelas et al., 2014; Newbold et al., 2015; Thuiller et al., 2019). However, most of these efforts have focused on aboveground communities, with little being done to understand the response of the global soil microbiome to contrasting scenarios of global change (Eisenhauer et al., 2019). Forecasted increases in aridity and temperature (Huang et al., 2017) will result in important changes in the major ecological drivers controlling belowground communities (e.g., reductions in plant cover and increases in pH; Slessarev et al., 2016). Such changes can, in theory, influence the trajectories of soil biodiversity and promote or reduce entire groups of soil taxa in future scenarios of global change (Delgado-Baquerizo & Eldridge, 2019; Fierer, 2017; Tedersoo et al., 2014). Increasing our knowledge about the future trajectories of soil biodiversity could be of great utility for land managers and decision-makers seeking to understand the sensitivity of soil biodiversity and ecosystem functions to global change and how this could affect society (Delgado-Baquerizo et al., 2020). Therefore, the aim of the present study was to produce the first integrated assessment of future projections for global bacterial communities under different scenarios of climate and land-use change.

Methods

Data collection and site characteristics

We used a global dataset (compiled from Delgado-Baquerizo, Oliverio, et al., 2018; Egidi et al., 2019) containing information on amplicon sequencing for bacteria (16S ribosomal RNA) and fungi [ribosomal internal transcribed spacer region (ITS)] for 231 locations world-wide. These sites include a wide range of ecosystem types and climatic regions across 18 countries and six continents. The mean annual precipitation and temperature in these locations ranged from 67 to 3,085 mm and from -11.4 to 26.5 °C, respectively. Soil sample collection took place between 2003 and 2015. At each site, a composite soil sample (top c. 7.5 cm depth) was collected under the most common vegetation. For all soil samples, the pH, texture and total organic carbon (soil C) concentration were measured using standard laboratory methods (Anderson & Ingram, 1993; Kettler et al., 2001). In brief, soil pH ranged from 4.04 to 9.21, soil C from .23 to 34.77%, and fine texture fraction (percentage clay + percentage silt) from 1.40 to 92.00%.

We obtained information on the mean annual temperature and precipitation for all sampling locations from the WorldClim database (www.worldclim.org), which has a 1 km resolution (Hijmans et al., 2005), and elevation data from the Global Multi-resolution Terrain Elevation Data 2010 (GMTED2010). We used an estimation of the fractional amount of green vegetation (hereafter, “vegetation cover”) as a proxy for plant coverage and plant primary productivity (Filipponi et al., 2018). This index provides a measure of the percentage of soil covered by green vegetation across Earth’s landscapes for a given composite period. Vegetation cover data were obtained from MODIS moderate-resolution images (Filipponi et al., 2018). We calculated the monthly average value for this variable for the period between 2001 and 2015 (c. 1 km resolution) to match the period of time for which these data are available and until the soil sampling was conducted.

Mapping and prediction datasets

For the temporal analyses, we used temporally explicit datasets coming from the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) (Hempel et al., 2013) and the Land-Use Model Intercomparison Project (LUMIP) (Lawrence et al., 2016) activities from the Intergovernmental Panel for Climate Change (IPCC), respectively. This selection followed the protocol laid out by Kim et al. (2018).

In terms of climate datasets, we used a bias-corrected historical and future (for three Representative Concentration Pathways: RCP2.6, RCP6.0 and RCP8.5) ISMIP2a dataset (Hempel et al., 2013) with forcing data of the IPSL-CM5A-LR Earth System Model (Dufresne et al., 2013). The monthly means of the daily temperatures (mean, minimum and maximum) and daily total precipitation > 1 mm were calculated for the available time period. For the purpose of this study, we selected the projection period between 1950 and 2100, in line with Kim et al. (2018). To avoid outliers and to make this dataset as comparable as possible with the original model-fitting datasets, we calculated 20-year climatologies using a moving window centred in each year step, resulting in a projection period between 1950 and 2090. The dataset created was used as climate input for all model runs. For each shared socio-economic pathway (SSP) × RCP combination, we used two different general circulation models (GCMs; gfdl-esm2m and noresm1-m) (Dufresne et al., 2013).

In the case of the land-use projections, we built on the dataset provided by the land-use Harmonized v.2.0 project (http://luh.umd.edu/) (Hurtt et al., 2011). This dataset was produced in the context of the World Climate Research Program Coupled Model Intercomparison Project 6 (CMIP6) (O’Neill et al., 2014; Popp et al., 2017; Riahi et al., 2017) and contains a harmonized set of land-use scenarios that are consistent between historical reconstructions and future projections. These modelled projections reproduce annual land-use reconstructions for historical land-use forcing (covering the period from 850 to 2015) and for different integrated assessment models (IAMs) and SSPs (from 2015 to 2100) at 0.25° resolution. These pathways represent a range of plausible futures based on different socio-economic challenges for mitigation of climate change (low in SSP1 and SSP4; high in SSP5) and potential challenges for adaptation (low in SSP1 and SSP5; high in SSP4). Full descriptions of the scenarios are given by Popp et al. (2017). We selected the combinations SSP1 × RCP2.6 (high climate-change mitigation and adaptation, with expected low CO2 concentrations; sustainability pathway), SSP4 × RCP6.0 (high climate-change mitigation, with challenges for adaptation and expected high CO2 concentrations; inequality pathway) and SSP5 × RCP8.5 (low climate-change mitigation, with high adaptation and expected high CO2 concentrations; fossil-fuelled development pathway). The land-use projections were harmonized for CMIP6 using an updated version of the land-use harmonization methodology (Hurtt et al., 2011), which was developed and used widely to support future biodiversity projections (Newbold et al., 2015; Powers & Jetz, 2019).

Correlation network and identification of dominant microbial phylotypes

We used correlation network analyses to identify ecological clusters formed by bacteria and taxa belonging to three major groups of fungi (mycorrhizal fungi, saprobes and pathogens) strongly co-occurring with each other. We included only relatively common taxa in our network analyses (present in > 25% of locations). A total of 1,916 phylotypes were included in our correlation network. To build the co-occurrence network, we first calculated pairwise Spearman’s rank correlations (ρ) between all common phylotypes. We focused exclusively on positive correlations, because they provide information on microbial phylotypes that might respond in a similar manner to environmental conditions and they are expected to co-occur with each other. We considered a co-occurrence to be robust if the Spearman’s correlation coefficient (ρ) was > .65 and p < .00001 (Barberán et al., 2012). This network was visualized with the interactive platform Gephi (Bastian et al., 2009). Finally, we used default parameters from the interactive platform Gephi to identify ecological clusters of soil bacterial and fungal taxa strongly co-occurring with each other (Bastian et al., 2009). We then computed the relative abundance of each ecological cluster by averaging the standardized relative abundances (ranging from zero to one) of the taxa that belong to each ecological cluster. By standardizing our data, we ruled out any effect of comparing the diversity from different soil groups (i.e., bacteria and three groups of fungi).

Functional predictions

We calculated functional predictions based on the representative genomes from our database (Ortiz-Álvarez et al., 2018). We used this representative genomes approach for the estimation of genomic and metabolic potential and obtained functional gene information from the Kyoto Encyclopedia of Genes and Genomes database (www.genome.jp/kegg/) using the Integrated Microbial Genomes & Microbiomes (IMG/M) system (https://img.jgi.doe.gov). We included in our analyses only those genomes that matched > 97% of a reference genome and were more than c. 90% complete. In our study, this information is used as a proof of concept; however, given the nature of our data, further studies are encouraged to confirm our observations based on alternative functional predictions. We calculated weighted communities, combining the functional data and the relative abundance matrix of genomic matches. We assessed the relative amount of genes associated with: (a) high-efficiency inorganic phosphate (P) transport (averaged standardized abundance of K02036, K02037 and K02038; pstA, pstB and pstC; hereafter, “phosphate transport genes”) associated with the capability of microbial communities to take up and mobilize inorganic phosphate [the abundance of this group of genes in our database is negatively correlated with total soil P (ρ = -.228; p < .001; n = 231) and positively correlated with the activity of phosphatase measured in a subset of our soil samples (ρ = .294; p = .011; n = 74)]; (b) enzyme activities (averaged standardized abundance of K01176, K01183, K01187, K01188 and K01205; alpha-amylase, chitinase, alpha-glucosidase, beta-glucosidase and alpha-N-acetylglucosaminidase genes); (c) carbon fixation [averaged standardized abundance of genes prkB, K00855 and rbcS, K01602; this group of genes was positively related to the relative abundance of soil organic carbon (ρ = .273; p < .001; n = 231)]; and (d) nitrogen fixation (averaged standardized abundance of genes anfG, K00531; nifD, K02586;, glnA, K01915; and nifK, K02591). This group of genes was positively related to the concentrations of soil total nitrogen (N) (ρ = .241; p > .001; n = 231).

Development of a meta-model of the global soil microbiome

Structural equation models use multiple structural equations to depict and model multivariate relationships (Grace, 2006). These can be used to represent and test indirect (cf. mediation) and direct relationships, partial contributions of correlated explanatory variables, and alternative hypotheses (Eisenhauer et al., 2015). This methodological approach has made its way into several scientific fields, including ecology, being used to understand complex systems and/or mechanistic relationships (Eisenhauer et al., 2015; Grace et al., 2016). The construction of a structural equation model depends on what is known or suspected about the variables describing the system being studied (Eisenhauer et al., 2015; Grace, 2008; Grace et al., 2016). Although for many soil organisms (e.g., protists, mites and Collembola) information may be lacking (Cameron et al., 2018), today there is an extensive pool of literature informing the drivers and range of conditions that favour the soil microbiome (Bahram et al., 2018; Delgado-Baquerizo & Eldridge, 2019; Delgado-Baquerizo, Reith, et al., 2018; Griffiths et al., 2016; Maestre et al., 2015; Ramirez et al., 2018; Tedersoo et al., 2014). We used this knowledge to develop and justify a meta-model describing the global relationships associated with the soil microbiome (Supporting Information Appendix S1). Attention was given to the potential direct and indirect effects of different drivers to account for mediator effects (Grace et al., 2012), particularly in the case of soil-driven effects, climate regulation and the impact of land use on the soil microbiome (Jansson & Hofmockel, 2020). For our final model, vegetation type, elevation, latitude, precipitation, temperature, soil carbon content and texture were selected as independent variables (i.e., predictors). Moreover, pH, the percentage of green vegetation (owing to their importance to determine general patterns) and the different diversity (i.e., bacterial richness and community dissimilarity) and functional (i.e., phosphate transport genes) indicators were treated as response variables and, thus, used for prediction (Supporting Information Appendix S1).

Once established, the assumptions of the meta-model were initially verified using a variance partitioning of the dataset to ensure that they depicted the main variance explained (Supporting Information Appendix S2). Variance partitioning was also used initially to preselect the biodiversity groups, functions and indicators that could be used in the proceeding analysis. For this, an initial cut-off value of 40% variance explained was used as a preselection value. From this analysis, we selected bacterial richness, community dissimilarity, the relative abundance of phosphate transport genes and the ecological clusters formed from the correlation network analysis previously described.

Model fit and structural equation modelling

Following the overall model description (Supporting Information Appendix S1), the initial variable pool was tested for normality and linearity, and selected variables were transformed accordingly (e.g., pH, elevation and texture were log10-transformed). The initial tests validated the assumptions of linearity across the different models although, in the case of the ecological clusters, our interpretation of these should be taken more conservatively. After the initial diagnostics, we proceeded with the fitting of four different structural equation models (i.e., for bacterial richness, community dissimilarity, phosphate transport genes and the ecological clusters). Before performing a model fit, we did a conditional independence test to account for possible missing links within the structural equation model formulation. To do this, we applied both the dagitty and piecewise packages in R (Lefcheck, 2016; Textor et al., 2016), and all potential missing links were identified and evaluated. We used the lavaan package to obtain a global estimation for the different structural equation models and used robust maximum likelihood parameter estimates (MLM) to account for potential bias in the input data and obtain the path coefficients of each model (Rosseel, 2012). The model fitting was done with independent variables measured in the field to avoid potential circularity associated with the use of modelled environmental data to fit causal relationships (see Supporting Information Appendix S1). We also took into account the potential influence of the available number of samples for modelling (Supporting Information Appendix S3).

A full description of the resulting structural equation models and their fit measures can be found in the Supporting Information (Appendix S3). Given that pH was found to have a strong effect on the distribution and temporal evolution of the soil microbiome, we performed an extra analysis to illustrate the determinants of pH using a second dataset (Batjes et al., 2017; Orgiazzi et al., 2018). A full description of the dataset and the subsequent analysis is given in the Supporting Information (Appendix S5).

Predictions in space and time

Using the parameter estimates obtained from the structural equation model fitting, we formulated two linear equations to predict the values of percentage vegetation and pH. Both the equations and the unstandardized coefficients are given in the Supporting Information (Appendix S3). In a similar way, we set up linear equations for bacterial richness, community dissimilarity, abundance of phosphate transport genes and each ecological cluster. We used the predicted values for percentage vegetation and pH in the overall prediction model to account for the indirect effects of climate and land-use change on the soil microbiome. The model runs were calculated with a yearly time step according to the available climate and vegetation data. All temporal changes (both backcasting and forecasting predictions) were calculated using 2015 as a baseline to which all predictions were compared. Although other microbial attributes were included in this global database (Delgado-Baquerizo, Oliverio, et al., 2018), we excluded them from further analyses, because our environmental models did not predict > 40% of their variance (for a detailed explanation of the methods, see Supporting Information Appendix S1 and S2).

In order to assess the accuracy of these predictions, we determined how much the parameter space of the predictors of the original dataset differed from the ones used to make predictions. For this, we used the Mahalanobis distance of any multidimensional point (here considering the seven dimensions given by the seven exogenous variables) to the centre of the known distribution (calculated based on the 231 data points from the original dataset) and the distance of any multidimensional point to the convex hull formed by the 231 data points used to estimate the model. Further explanation of the method used can be found in the Supporting Information (Appendix S4).

Accuracy of projections and distribution of outliers

One of the difficulties of predicting response variables is the extrapolation of results to areas where the environmental conditions of the initial dataset (used to estimate the structural equation models) are not covered. Attempting to use predictive equations to estimate areas that differ considerably from the original dataset can lead to unreliable predictions (Brooks et al., 1988). Although this is still a current topic of research interest in multidimensional statistical analysis (Ebert et al., 2014), we propose to use the Mahalanobis distance to identify outliers with regard to the centre of the observed distribution given by the 231 points present in our dataset.

The Mahalanobis distance of a point to a dataset of points in multidimensional space is calculated using the vector of means of the dataset and its covariance matrix (Figure 1). This distance is often used to detect outliers in point cloud distributions that are assumed to follow a multivariate normal distribution (Jackson & Chen, 2004; Rousseeuw & van Zomeren, 1990). When each of the variables is normally distributed, the Mahalanobis distance follows a chi-square distribution with d degrees of freedom, where d is the dimension of the multidimensional space (d = 7 in our case). To inform the quality of our predictions, we computed the Mahalanobis distance of each input point to the dataset of 231 points. The Supporting Information (Appendix S4, Figure S4-1) uses a colour gradient to indicate the quantile of the chi-square distribution with seven degrees of freedom that each point belongs to (Mallavan et al., 2010). Although this distance is an informative measure of how close a new data point is to the distribution of points used to estimate the model, we used a second measure to assess whether a new data point is an extrapolation or an interpolation. We used the resulting outlier identification to mask our outputs and provide more reliable predictions.

Figure 1

Quantiles of the chi-square distribution with seven degrees of freedom calculated for the conditions in 2015 based on the Mahalanobis multidimensional distance (corresponding to the seven dimensions determined by the exogenous variables). Using this method, we can determine that 60.6% of the terrestrial systems fall under the .975 quantile, with the rest being considered as outliers to the original distribution of 231 data points

Independent validation

We used data from an independent global survey (Delgado-Baquerizo et al., 2019) that includes 80 locations across the globe to validate further the spatial distributions of bacterial richness across samples. This dataset includes soils from contrasting vegetation (forests, grasslands and shrublands) and climates (arid, temperate, tropical and continental) from five continents and eight countries. Note that any comparisons between this independent dataset and our dataset need to be considered carefully, given methodological differences in the primer sets used (here 16S 341F/805R vs. 16S 515F/806R). Soils are comparable in terms of sampling design and soil depth. We calculated bacterial richness and community dissimilarity using a similar approach as for the 16S 341F/805R data (for further details, see Appendix S1).

Results and Discussion

We estimated the distribution of multiple components of the soil microbiome, providing global projections for soil organisms in the Anthropocene similar to those existing for plant and aboveground animals (Newbold et al., 2015; Powers & Jetz, 2019). Within these, we found soil pH and vegetation cover to be the major direct ecological predictors of the distribution of the soil microbiome, particularly for bacterial richness and community dissimilarity (Supporting Information Appendix S3; Figure S3-1). The spatial distribution of both bacterial richness and the ecological clusters follows the ones obtained in current literature using different methods and datasets (Figure 2; Bahram et al., 2018; Delgado-Baquerizo & Eldridge, 2019; Delgado-Baquerizo, Oliverio, et al., 2018; Ramirez et al., 2018). Our results also show that globally, community dissimilarity is negatively correlated with bacterial richness (Spearman’s ρ = .904, p < .001). This suggests that the low local bacterial richness shown in several regions of the globe (Delgado-Baquerizo & Eldridge, 2019; Delgado-Baquerizo, Oliverio, et al., 2018) (e.g., in the tropics or in high latitudes) might be related to higher local differentiation; that is, although locally we are predicting lower (alpha) diversity, the compositions of these communities are expected to be highly differentiated from each other, resulting in higher regional (gamma) diversity (Figure 2).

Figure 2

Global distribution of: (a) bacterial richness, (b) community dissimilarity, and (c–f) the ecological clusters A (c), B (d), C (e) and D (f). All maps are for the baseline of 2015 and are classified using a quantile distribution. Colours correspond to each of the quantiles, from low (first quantile) to high (fourth quantile) values

Using the spatial predictions for current time and cross-validating them with this new dataset, we obtained significant positive correlations for bacterial richness (Spearman’s ρ = .236, p = .035, n = 80) and community dissimilarity ((Spearman’s ρ = .241, p = .031, n = 80). Given the completely different methods associated with each independent dataset [training (Delgado-Baquerizo, Oliverio, et al., 2018) and validation (Delgado-Baquerizo et al., 2019)], the relationships found provide further support for the validity of our approach to estimate and predict the global soil microbiome. Nevertheless, given the spatial and temporal resolution of our predictions, our models and maps should be used only to explore broad patterns and trends of the soil microbiome and not for fine-scale (pixel-level) analyses, given the lower reliability of predictions in some regions (Figure 1).

In terms of the ecological clusters, cluster A defines a highly diverse ecological cluster composed of both fungi and bacteria and is associated with areas with high temperature variation, low precipitation and alkaline soils with low carbon concentrations. It is prevalent mostly in drylands, including deserts, and mountain grasslands and shrublands in Africa, Central Asia, Australia and South America. This ecological cluster is mainly affected positively by vegetation structure and pH, whereas latitude and vegetation cover correspond to the main negative effects (see Supporting Information Appendix S3, Table S3-3). Cluster B, which includes 10 bacterial taxonomic groups (see Supporting Information Appendix S2, Figure S2-2), occupies a similar spatial range to cluster A, although it extends further north, including mediterranean and some temperate to dry systems with low vegetation cover and high pH. Cluster A has a higher prevalence in South America, whereas cluster B expands further, to Central and North America. The latter is mainly affected (positively) by pH, climate (directly by temperature and indirectly by precipitation) and latitude, with vegetation cover having the only negative effect.

In comparison to the previous two clusters, the ecological clusters C and D are related mainly to more humid regions. Cluster D reflects the lowest diversity of the four clusters (six taxonomic groups of bacteria) and is concentrated mainly in temperate forests and grasslands. This cluster is characterized by having positive effects from both pH and vegetation cover, which makes it more prevalent in areas with higher vegetation cover combined with more alkaline soils. At the same time, it shows negative effects from both temperature and vegetation structure (see Supporting Information Appendix S3). Finally, cluster C is associated mainly with soils with lower pH, with climate and land use having only indirect effects. Overall, cluster C is distributed across humid systems with low pH and high vegetation cover. The overall model fit for bacterial richness returned an R 2 of .449 (p = .666, rmsea (root mean square error of approximation) = .000, cfi (comparative fit index) = 1.000). In the case of community dissimilarity this was .702 (p = .581, rmsea = .000, cfi = 1.000) and for phosphate transport genes .541 (p = .299, rmsea = .046, cfi = .996). Finally, for the ecological clusters we obtained an R 2 for A, B, C and D of .677, .681, .666 and .416, respectively (p = .148, rmsea = .046, cfi = .996; for further information, see Supporting Information Appendix S3, Tables S3-1 and S3-2).

Our projections indicate major changes in bacterial diversity and a global-scale homogenization process for the soil microbiome. Overall, we found that the local richness of bacteria will increase in future conditions [especially for scenarios SSP4 (regional inequality) and SSP5 (fossil-fuelled development)]. This increase of bacterial richness with land-use intensification is in line with other studies carried out at national scales (Karimi, Terrat, et al., 2018). However, the community composition of bacteria is expected to undergo a strong and generalized global homogenization processes (i.e., a decrease in community dissimilarity, with local communities becoming more similar to each other) affecting > 85% of terrestrial ecosystems assessed (Figure 3; Supporting Information Appendix S4, Figure S4-2). This homogenization trend is consistent with previous studies warning about a rapid global homogenization of marine and aboveground biodiversity (Dornelas et al., 2014; Magurran et al., 2015).

Figure 3

Global projections of the soil microbiome show large-scale anthropogenic homogenization. Relative changes (gains and losses) in the historical patterns (a,b) and future projections [scenario agreement (2015–2090); e,f] for bacterial richness (left) and community dissimilarity of bacteria (right). (c,d) Projection lines correspond to the global median value with interquartile range (IQR) for the historical (1910–2015; in grey) and future projections (IQR bars for each scenario). For the shared socio-economic pathways (SSP), SSP5 corresponds to a fossil-fuelled development scenario, SSP4 to an inequality scenario and SSP1 to a sustainability scenario (Supporting Information Appendix S4). All projections represent combined land-use and climate effects

Future ecosystems are, therefore, expected to share a larger number of bacterial phylotypes at the local scale, making several bacterial taxa potentially more abundant in soil communities under global change scenarios (Karimi, Terrat, et al., 2018). This result contrasts with the often-reported negative associations between global change drivers and the local biodiversity of aboveground communities (e.g., plants and mammals; Di Marco et al., 2019; Newbold et al., 2015; Powers & Jetz, 2019). Despite enhanced local bacterial richness, future communities will be based on a more homogeneous species pool in the majority of terrestrial systems, reducing the dissimilarity of bacterial communities across locations (Figure 3d). The observed bacterial richness and community dissimilarity trends were associated with increases in soil pH linked to important shifts in precipitation and temperature, and reductions in vegetation cover that are negatively associated with soil pH. Such environmental dynamics will select for species that benefit from higher soil pH.

These changes in soil biodiversity might also result in important changes in soil ecosystem functions (Bardgett & van der Putten, 2014) related to the spread of specific functional genes (Delgado-Baquerizo et al., 2020). The lower bacterial dissimilarity across locations, in addition to the often-reported strong link between microbial taxonomic and functional community composition (Allison & Martiny, 2008; Fierer et al., 2012, 2013; Torsvik & Øvreås, 2002), suggest that cross-site functional redundancy might increase in the near future, assuming that the links between functionality and taxonomy remain constant through time. Thus, similar bacterial taxa with similar functional capabilities will live in soils across the globe, reducing specialization and, potentially, the capacity of ecosystems to adapt to new environmental realities. Nevertheless, we acknowledge that the assumption of constant functional links might be flawed given that, while microbial communities change, their functional links might also result in adaptations to the new climatic and land-use conditions. Although this aspect deserves further study, we also observed combined gains in bacterial richness and community dissimilarity in 7.3% of the globe, considering historical trends (green areas in Supporting Information Appendix S4, Figure S4-4c). These areas result from less dramatic change in land cover and more positive effects coming from climate change that resulted in net gains for bacterial richness and a predicted gain in community dissimilarity. Therefore, the microbial communities in these areas are expected to become more unique from their counterparts.

Although limited in functional scope, our projections also show changes in the relative abundance of functional genes associated with the increases in bacterial richness. For example, our global projections support an increase in the relative abundance of phosphate transport genes (positively associated with activity of phosphatase) in many regions across the globe, particularly after 1970 with an increase in land-cover change, especially deforestation (Supporting Information Appendix S4, Figure S4-5). Such an increase in the abundance of genes associated with phosphate transport indicates a larger capacity to mobilize and take up soil phosphorus, which is one of the most important factors limiting plant and microbial production (Peñuelas et al., 2013). Our initial statistical threshold to exclude variables from being predicted (see Section 2.5) limits our ability to go further in the description of future functional consequences. Nevertheless, these first insights into shifts in microbial community composition call for further work on the functional consequences of changes in the soil microbiome under global change scenarios.

Our projections further provide evidence for the existence of microbial taxa that will prevail under global change scenarios. These results are based on the relative abundance of ecological clusters of co-occurring bacterial and fungal taxa within a global correlation network (for a description of these clusters, see Supporting Information Appendix S2; for a list of taxa included in each cluster, see Supporting Information Appendix S8). We found two main ecological clusters that we predict to have important gains (cluster B) or losses (cluster C) under current and, even more so, future conditions (Figure 2a,b). The expected reductions in the relative abundance of cluster C will lead to a decrease in one of the most dominant taxa found in soil (Bradyrhizobium spp.) and in functionally important taxa, such as Methylocystaceae (including methanotrophs), both of which are known to prefer low-pH environments (Delgado-Baquerizo, Oliverio, et al., 2018).

In contrast to clusters B and C, the trends of the relative abundance of clusters A and D changed relative to the past (1910–2015) from losses to substantial increases (cluster A; the most diverse cluster; Supporting Information Appendix S2) or from gains to decreases (cluster D; the least diverse cluster; Supporting Information Appendix S2) under global change scenarios (Figure 2b,c). Our projections also indicate that the trends identified for the historical period will become more abrupt in the future, resulting in stronger shifts in community composition towards communities driven by an expansion of drylands (clusters A and B). These future communitiestend to be less specialized and to aggregate a wider group of more common microbial taxa. In fact, the relative losses projected for the communities associated with more humid environments (cluster C) show generalized losses independently of the scenario considered, with exception of Southeast Asia.

These results suggest that taxa included in cluster A will become more common in the future. This cluster includes phylotypes such as Geodermatophilus spp. (typical desert bacteria; Delgado-Baquerizo, Oliverio, et al., 2018), Mycobacterium spp. (which are known to include important human pathogens; Delgado-Baquerizo et al., 2020), Streptomyces mirabilis (major producers of antibiotic resistance genes) or potential fungal soil-borne plant pathogens belonging to Ascomycota fungi (Venturia spp. and Devriesia spp.). Our results also indicate that the relative abundance of cluster D might have important local reductions (with a high model agreement across all scenarios; Figure 4c), which include consistent declines in particular bacterial taxa (e.g., family Gaiellaceae). For most of the taxa included in this cluster, we do not even know their genus, which makes it difficult to evaluate the potential consequences of such a reduction in their relative abundance. Although, in this analysis, we did not consider the relationships between biodiversity and the magnitude or combination of multiple soil ecosystem functions, this projected community shift is expected to have strong effects on the multifunctionality of soils across the globe, with potential implications for human health and wellbeing (Evans & Wallenstein, 2014; Kardol et al., 2010).

Figure 4

Projected changes in ecological clusters (including both bacteria and fungi) show divergent trends when considering past and future conditions.

(a) Distinction between clusters with relative net gains in relative abundance and clusters with relative losses. Only clusters with the largest net changes are represented within each pixel. (b) Future projections of the different ecological clusters according to the three shared socio-economic pathways considered. (c) Scenario agreement (2015–2090) in terms of gains and losses for clusters A and D. For the shared socio-economic pathways (SSPs), SSP5 corresponds to a fossil-fuelled development scenario, SSP4 to an inequality scenario and SSP1 to a sustainability scenario (Supporting Information Appendix S4). All projections represent combined land-use and climate effects

Conclusions and Implications

Using structural equation models, we were able to differentiate between the potential interactions between climate (e.g., increases in temperature) and land-use effects (e.g., changes from forest to grasslands) on the soil microbiome. We also assessed how these two drivers affect the soil microbiome, both directly and indirectly, through changes in vegetation cover and pH. Across the different scenarios considered, we found that climate change has a greater effect than land-use intensification on both bacterial richness and community dissimilarity. In the case of bacterial richness, these changes reflect only the indirect effects of climate change (Supporting Information Appendix S3). The increase in local bacterial species richness results only from the direct and indirect effects of land use and climate. However, other factors, such as competitive exclusion or effects of microhabitat, could prevent some of the species from co-existing, introducing a potential overestimation in the “increase” trends of local species richness.

For some scenarios (for SSP1, see Supporting Information Appendix S4), the trend obtained from the projections changed from positive (in the case of bacterial richness) to negative and vice versa, when considering only land use as opposed to both drivers of global change (land use and climate change). Such an effect is likely to be associated with the indirect effects of climate and land-use change on vegetation cover (e.g. vegetation cover decreases owing to deforestation and increases in aridity) and pH (e.g., pH increases owing to an increase in temperature and decrease in precipitation). These results suggest that although researchers often focus on changes between types of land use when assessing impacts on soil communities, at a macroecological scale the trends associated with these microbial communities are mostly driven indirectly by climate (because most changes occur by indirect effects through vegetation cover and pH). Therefore, models not considering the indirect effects of global change drivers might be ineffective in predicting changes in the soil microbiome (Singh et al., 2010).

Our results provide global projections for belowground communities under global change scenarios and forecast major changes in microbial diversity (local increases in bacterial diversity and a global homogenization of the soil microbiome) driven by changes in climate and land use. These findings contrast with current global projections of aboveground biodiversity declines, but they do not necessarily provide a more positive view of the future of nature. Although increases in local microbial diversity might seem positive, these hide strong reductions in community complexity (Kortz & Magurran, 2019) in the majority of terrestrial systems, with implications for ecosystem functioning. This belowground trend is now mostly missing from policy documents and nature conservation assessments, which diminishes the capacity of policy-makers to make informed decisions about the future of soil ecosystems. Our global projections of the soil microbiome provide a baseline for future discussion on the importance, trajectories and conservation needs of soil organisms. While informing the biodiversity decline debate, by showing the spatial and temporal trends of an important group of soil organisms, our results also contribute to closing important gaps identified in current global assessments and conventions (e.g., IPCC, IPBES and UNCCD) and lay the foundation for inclusion of soil organisms in future assessments of the response of ecosystem to drivers of global change.

Supplementary Material

Additional supporting information may be found online in the Supporting Information section.