Spatial turnover of multiple ecosystem functions is more associated with plant than soil microbial β -diversity

. Biodiversity — both above-and belowground — in ﬂ uences multiple functions in terrestrial ecosystems. Yet, it is unclear whether differences in above-and belowground species composition ( β - diversity) are associated with differences in multiple ecosystem functions (e.g., spatial turnover in ecosystem function). Here, we partitioned the contributions of above-and belowground β -diversity and abiotic factors (geographic distance, differences in environments) on the spatial turnover of multiple grassland ecosystem functions. We compiled a dataset of plant and soil microbial communities and six indicators of grassland ecosystem functions (i.e., plant aboveground live biomass, plant nitrogen [N], plant phosphorus [P], root biomass, soil total N, and soil extractable P) from 18 grassland sites on four continents contributing to the Nutrient Network experiment. We used Mantel tests and structural equation models to disentangle the relationship between above-and belowground β -diversity and spatial turnover in grassland ecosystem functions. We found that the effects of abiotic factors on the spatial turnover of ecosystem functions were largely indirect through their in ﬂ uences on above-and belowground β -diversity, and that spatial turnover of ecosystem function was more strongly associated with plant β -diversity than with soil microbial β - diversity. These results indicate that changes in above-and belowground species composition are one mechanism that interacts with environmental change to determine variability in multiple ecosystem functions across spatial scales. As grasslands face global threats from shrub encroachment, conversion to agriculture, or are lost to development, the functions and services they provide will more strongly converge with increased aboveground community homogenization than with soil microbial community homogenization.


INTRODUCTION
Biodiversity influences a variety of ecosystem functions such as primary productivity, litter decomposition, and carbon and nutrient cycling in terrestrial ecosystems (Wardle et al. 2011, Naeem et al. 2012).Since the 1990s, numerous experimental studies have documented that some dimensions of biodiversity (typically localscale plant species richness, hereafter α-diversity) promote plant primary productivity (Cardinale et al. 2006, Schmid et al. 2009, Duffy et al. 2017).However, most early studies focus on how biodiversity loss influences a single ecosystem function, such as a nutrient pool size or the rate of an ecosystem function (but see Hector andBagchi 2007, Gamfeldt et al. 2008).More recent work demonstrates that biodiversity positively affects an ecosystem's ability to support multiple ecosystem functions (i.e., ecosystem multifunctionality; Hector and Bagchi 2007, Maestre et al. 2012, Manning et al. 2018).When species are lost from a community, the capacity of that community to perform multiple ecosystem functions may be diminished.Connecting biodiversity to multifunctionality therefore enhances our ability to understand the consequences of biodiversity change for sustaining overall functioning in realworld ecosystems (Gamfeldt et al. 2008, Eisenhauer et al. 2016, Mori et al. 2018, Manning et al. 2019).
Most experimental studies of biodiversity and ecosystem function research have been conducted at relatively small spatial scales (van der Plas 2019) and have focused on how changes in α-diversity affect individual ecosystem functions or an aggregated measure of multiple ecosystem functions, that is, α-multifunctionality (van der Plas et al. 2016Plas et al. , H ölting et al. 2019)).While it is still controversial whether species are being lost at the spatial scale at which these empirical studies are conducted (Cardinale et al. 2018), it is clear that species composition is changing (Wardle et al. 2011, Bannar-Martin et al. 2018, Blowes et al. 2019).For example, β-diversity, spatial turnover in species composition among communities (Anderson et al. 2011), increases with species heterogenization in some areas, whereas it decreases as communities become homogenized in other areas (Arroyo-Rodríguez et al. 2013).Clearly, differences in species composition (i.e., β-diversity) can lead to differences in ecosystem function (Pasari et al. 2013, Grman et al. 2018, Hautier et al. 2018, Winfree et al. 2018): A pine plantation and a tropical rainforest have very different functions.But, it remains unclear whether differences in function that arise in these divergent ecosystems are the product of differences in aboveground or belowground community structure (Talbot et al. 2014, Gravuer et al. 2020).There is evidence that plant (van der Plas et al. 2016, Hautier et al. 2018), soil fungal (Mori et al. 2018), and soil multi-trophic (Martinez-Almoyna et al. 2019) β-diversity are broadly linked to spatial turnover in multiple ecosystem functions (i.e., β-multifunctionality, the pairwise differences in multiple ecosystem functions (H ölting et al. 2019, Peters et al. 2019, van der Plas et al. 2016)) from local to regional scales.But again, it remains unclear whether spatial turnover in multiple ecosystem functions arises more strongly from differences in aboveground community structure (aboveground β-diversity) or differences in belowground community structure (belowground β-diversity).
In this study, we examined the relative importance of spatial differences in above-and belowground species composition (β-diversity) on the spatial turnover of multiple ecosystem functions in grasslands.Specifically, we predicted that (1) plant and soil microbial β-diversity jointly determine the spatial turnover of multiple ecosystem functions, while (2) abiotic factors mediate the relationship between β-diversity and spatial turnover in multiple ecosystem functions.To test these hypotheses, we used data from 18 grassland sites on four continents (Appendix S1: Table S1, Fig. S1).At each site, we collected data on plant and soil microbial communities as well as six indicators of ecosystem function that are associated with biological productivity and nutrient pools in both above-and belowground compartments of grassland ecosystems: plant aboveground biomass, plant nitrogen (N), plant phosphorus (P), root biomass, soil total N, and soil extractable P (Garland et al. 2021).

Data collection
We collected data from a subset of 18 sites from the Nutrient Network presented in Prober et al. (2015; Appendix S1: Fig. S1, Table S1).The Network is a globally coordinated experiment to investigate the impacts of human activities (e.g., nutrient fertilization) on grassland structure and function (Borer et al. 2014).In this study, we used a total of 47 unfertilized control plots (1 m × 1 m) of 18 Nutrient Network sites, which were selected based on the availability of data characterizing the soil microbial communities and plant N and P. The unfertilized plots were established between 2007 and 2011, and included a single plot in three sites, two plots in three sites, three plots in 11 sites, and five plots in one site (Appendix S1: Table S1).The dataset covers 11 ecosystem types (alpine grassland, annual grassland, mesic grassland, montane grassland, old field, pasture, savanna, semiarid grassland, shortgrass prairie, shrub steppe, and tallgrass prairie) across four continents (Africa, Australasia, Europe, and North America).The diversity of study sites represents a wide range of climatic and environmental conditions across grasslands worldwide with mean annual precipitation ranging from 262 to 1898 mm, mean annual temperature ranging from 0.3 to 18.4°C, soil pH ranging from 4.6 to 8.3, soil total N ranging from 0.03 to 1.38%, soil extractable P ranging from 1 to 253 ppm, and plant productivity ranging from 23 to 1022 gÁm −2 Áyr −1 .
Climate data were obtained from the bioclimatic dataset of WorldClim (version 1.4; data source http://www.worldclim.org;Hijmans et al. 2005).We used the 30 arc-second resolution (~1 km) climate data to estimate the spatial variation in climate for the 18 sampling sites.Mean annual temperature (BIO1) and precipitation (BIO12) between 1950 and 2000 were extracted for further analysis.
During the growing season, plant aboveground biomass (gÁm −2 ) was collected in 2011 and 2012 from two 0.1-m 2 strips within an unfertilized 1-m 2 plot.Aboveground plant biomass was separated into dead and live and dried at 60°C to a constant mass.Samples for the measurement of plant aboveground N and P were collected between 2010 and 2014 from the same unfertilized plots as aboveground plant biomass and sorted by functional group (annual/perennial grasses, other graminoids, forbs, and legumes).Content of plant aboveground N and P was measured following the methods presented by Anderson et al. (2018).The content of plant aboveground N and P and the live biomass of plant functional groups were then used to calculate plant aboveground N (%) and P (%).The percent coverage of all vascular plant species was visually estimated to the nearest 1% within the unfertilized 1-m 2 plot in the field.Maximum percent coverage was assigned for those species that have two growth peaks within a growing season in 2011 or 2012.
Soil and root samples (five soil cores 2.5 cm in diameter × 10 cm in depth) were collected from the 0.1-m 2 strips the same as the plant aboveground biomass harvest plots.Roots were washed and dried at 40°C to a constant mass (Cleland et al. 2019).Soil pH was measured using a pH probe (Fisher Scientific Waltham, Massachusetts, USA) in a ratio of 2:5 dry soil to deionized water.Soil total N (%) and extractable P (ppm) were measured using the same methods as plant N and P. Soil DNA was extracted by the MoBio PowerSoil kit (Mo Bio Laboratories, California, USA) and amplified using the V4 region of the 16S rRNA gene with the 515f/806r primer set for bacteria and the first internal transcribed spacer (ITS1) region of the rRNA gene with ITS1-F/ITS2 primer set for fungi.Samples were sequenced on the Illumina HiSeq and MiSeq instruments at the University of Colorado Next Generation Sequencing Facility.The molecular and bioinformatics analyses were processed following the protocol as described by Prober et al. (2015).All the experimental protocols are provided on the website of Nutrient Network (https://nutnet.org/methods).

Quantifying geographic and environmental distance
We calculated great circle (geographic) distance (km) among sites using the longitude and latitude of the geographic locations.We calculated environmental distance using Euclidean distance (i.e., the square root of the sum of squared differences between environmental variables of two plots).Mean annual temperature and mean annual precipitation were considered as key indicators of climatic differences among sites (Prober et al. 2015).We also used a third variable (soil pH) to account for environmental differences among sites (Fierer and Jackson 2006).We computed climatic and soil pH distance matrices separately.Prior to the calculation of environmental distance, we standardized environmental variables by taking their rangerelevant standardization for each variable as follows: where X represents the values of the environmental variable in the original raw dataset, X min and X max represent the minimum and maximum values of the environmental variable, respectively.This standardization approach has the advantage that it allows the environmental variables to have different means and standard deviations but with equal ranges (Grace et al. 2018).

Quantifying community β-diversity
We quantified community β-diversity using two approaches.Firstly, we calculated three indices of β-diversity for plant, soil bacterial, and fungal communities: Sorensen index, Horn index, and Morisita-Horn index.The three indices are a function of order q which determines the sensitivity of a diversity index to rare and abundant species (Jost 2007).Specifically, Sorensen index is a β-diversity index of q order zero and uses species presence/absence data (i.e., weights are equal for both rare and abundant species).Horn index (q order one) is weighted in proportion to species frequency in a community and Morisita-Horn index (q order two) assigns abundant species more weight than rare species.
Secondly, we decomposed the Sorensen index of β-diversity (i.e., total β-diversity) into two complementary indices, that is, richness difference and species replacement.Richness difference refers to the differences in the number of species while species replacement refers to the pure species turnover among communities (Podani andSchmera 2011, Legendre 2014).The two components of β-diversity measure the relative richness difference and species replacement for any pair of communities that can be linked to the functions and processes of ecosystems (Marini et al. 2013, Legendre 2014).We used species presence/absence data following the methods proposed by Podani and Schmera (2011).The beta.div.compfunction (Legendre 2014) was used to compute richness difference and species replacement.

Quantifying spatial turnover in multiple ecosystem functions
We calculated a multivariate index of spatial turnover of multiple ecosystem functions using Euclidean distance (Martinez-Almoyna et al. 2019, Peters et al. 2019).We used Euclidean distance because it is metric and can be derived by tracing back to each individual variable based on the triangle inequality (Goslee 2010; see sensitivity analysis in Appendix S1: Fig. S2).Aboveground plant live biomass, plant N, plant P, root biomass, soil total N, and soil extractable P were selected as the indicators of ecosystem functions (Spehn et al. 2005, Allan et al. 2013, Jing et al. 2015, Delgado-Baquerizo et al. 2016, Hautier et al. 2018, Garland et al. 2021).To quantify whether spatial turnover in multiple ecosystem functions is different from the null expectation, we developed a null model to calculate a standardized effect size for each pairwise Euclidean distance (L2 SES ) as follows: where L2 obs is the observed pairwise Euclidean distance, L2 sim is the simulated pairwise Euclidean distance, μ(L2 sim ) is the mean of the simulated Euclidean distance, and σ(L2 sim ) is the standard deviation of the simulated Euclidean distance.
To calculate L2 SES , we used the matrix of the six ecosystem functions (functions are in columns and plots are in rows).We generated a randomized matrix of ecosystem functions by taking samples without replacement.This v www.esajournals.orgrandomization process was to constrict our simulations to the levels of ecosystem functions within a sample which assumed that the values of an ecosystem function were randomly drawn from the distribution of this specific set of six ecosystem functions.We repeated this 1000 times.The mean and standard deviation of the simulated Euclidean distance were then calculated.We classified each pairwise Euclidean distance into three groups.That is, |L2 SES | > 1.96 represents that the spatial turnover of multiple ecosystem functions is significantly greater (group I) or less (group II) than the null expectation.|L2 SES | < 1.96 represents that the spatial turnover of multiple ecosystem functions is not different from the null expectation (group III).Furthermore, if L2 SES ≥ 1.96, the spatial turnover of multiple ecosystem functions is large enough to detect a significant difference between two communities.If L2 SES < 1.96, the observed value of Euclidean distance is in the range expected by chance or is less than the null expectation, thereby we considered the spatial turnover of multiple ecosystem functions is not large enough to detect a significant difference between two communities.

Statistical analyses
We conducted sensitivity analysis using simple Mantel tests to select the indices of β-diversity and the number of ecosystem functions.Firstly, we inspected the bivariate associations between β-diversity and the spatial turnover of multiple ecosystem functions.Five indices of β-diversity (Sorensen index, the replacement and richness difference of Sorensen index, Horn index, and Morisita-Horn index) were separately assessed for each organismal type (plant, bacteria, and fungi).We calculated Spearman correlation coefficients (rho) to account for the non-linear bivariate associations.The replacement and richness difference had lower correlations with the spatial turnover of multiple ecosystem functions than Sorensen index; Horn index; and Morisita-Horn index were not significantly different from Sorensen index in most cases (Appendix S1: Table S2).Therefore, we selected the Sorensen index for further analyses.Secondly, we determined how the bivariate associations between plant β-diversity and the spatial turnover of multiple ecosystem functions change with the number of functions.
We calculated all combinations of the six indicators of ecosystem functions and calculated the Spearman correlation coefficients using simple Mantel tests.After inspection, the bivariate correlations increased with the number of ecosystem functions (Appendix S1: Fig. S2; see Jing et al. 2020 for more information about the relationships between biodiversity effects and the number of ecosystem functions considered).We therefore used the multivariate index of spatial turnover in multiple ecosystem functions combining all six indicators of ecosystem functions.
To examine whether geographic, climatic, and pH distances were associated with above-and belowground β-diversity, we conducted simple Mantel tests.In addition, we used partial Mantel tests controlling for the influences of other abiotic factors to examine the bivariate associations between geographic and environmental distances with β-diversity.We used the same methods examining whether abiotic factors and plant and soil microbial β-diversity were associated with spatial turnover in multiple ecosystem functions.Specifically, simple Mantel tests were used to assess the bivariate associations and partial Mantel tests were used to assess whether abiotic factors influence the bivariate associations between β-diversity and spatial turnover in multiple ecosystem functions.
To overcome the drawbacks of correlative inference from Mantel tests and to tease apart the relative influences of geographic distance, environmental distance and β-diversity of plant, soil bacteria, and soil fungi on spatial turnover in multiple ecosystem functions, we used structural equation models (SEMs).A prior framework for testing the hypothesis was built upon previous work to address the relative importance of past and current environmental change on aboveand belowground community assemblages and ecosystem functions (Barnes et al. 2016, Delgado-Baquerizo et al. 2016, Martinez-Almoyna et al. 2019, Prober et al. 2015).We specifically assumed that (1) both geographic distance and environmental distance affect above-and belowground β-diversity.That is, the legacies of historical events and contemporary environmental conditions determine the biogeography of plant and soil microbial community assemblages (Martiny et al. 2006) diversity (Appendix S1: Fig. S3).We tested two SEMs, one tested the hypothesis that plant βdiversity influences soil microbial β-diversity and the other tested the opposite directionality (that soil microbial β-diversity influences plant βdiversity).We excluded soil pH distance from the SEMs because the influence of soil pH distance on β-diversity and spatial turnover in multiple ecosystem functions was substantially smaller than geographic and climatic distances.
In addition, we included two residual correlations between geographic distance and climatic distance and between soil bacterial and soil fungal β-diversity.The significance tests for the path coefficients were estimated by a bootstrap procedure with 9999 random samples.Non-significant paths were removed sequentially from the candidate models.As a result, we kept only path coefficients that were statistically significant.The goodness-of-fit was assessed using several indices including χ 2 tests (P > 0.05), comparative fit index (CFI > 0.90), root mean square error of approximation (RMSEA < 0.10), and standardized root mean square residual (SRMR < 0.10; Grace 2020, Rosseel 2012).All statistical analyses were carried out in R version 3.6.1 (R Development Core Team 2019).All Mantel tests were performed with 9999 permutations.We extracted the climate data using the raster package (version 2.8-4; Hijmans 2018).We calculated the geographic distance using rdist.earthfunction in the fields package (version 8.4-1.6;Furrer et al. 2015), the Euclidean distance using distance function in the ecodist package (version 2.0.7;Goslee and Urban 2007), and the Sorensen index, Horn index, and Morisita-Horn index using sim.

RESULTS
We first compared the biogeography of aboveand belowground taxa.We determined the relative importance of geographic distance vs. differences in climate and soil pH on plant, soil bacterial, and soil fungal community assemblages.Specifically, we found that geographic distance was more positively associated with plant and soil microbial β-diversity than was environmental distance (Fig. 1a).The result was maintained even if we controlled for the influences of climatic distance, differences in soil pH or both (Fig. 1b).Differences in climate were a strong predictor of plant and soil microbial βdiversity (Fig. 1a).However, when geographic distance was controlled, plant β-diversity was not significantly associated with climatic distance (Fig. 1b).Differences in soil pH were more positively associated with soil bacterial β-diversity (Spearman correlation coefficient, hereafter, ρ = 0.34) than with soil fungal (ρ = 0.15) or plant β-diversity (ρ = 0.07; Fig. 1).
Next, we assessed whether the geographic and environmental distances meditate the relationship between β-diversity and spatial turnover in multiple ecosystem functions.Plant β-diversity (ρ = 0.34) was a stronger predictor of the spatial Points in red represent spatial turnover in multiple ecosystem functions that are significantly higher than null expectation, in blue lower than null expectation, and in gray no difference from null expectation.(c) Partial Mantel tests for the bivariate associations between β-diversity and spatial turnover in ecosystem function given geographic and environmental turnover of multiple ecosystem functions than were geographic (ρ = 0.23), climatic (ρ = 0.15), or soil pH (ρ = 0.09) distances (Fig. 2a, b).For belowground communities, we found that soil bacterial (ρ = 0.22) and soil fungal (ρ = 0.23) βdiversity were as important, or even more important, than geographic and environmental distances (Fig. 2a, b) in accounting for variability in the spatial turnover of multiple ecosystem functions.
Most interestingly, our null models showed that, for any two communities, whether they were close to one another or far apart (or with either a low or high environmental difference), spatial turnover in ecosystem function could be high or low (i.e., the red points in Fig. 2a).This pattern arises because the highest non-null values of spatial turnover in multiple ecosystem functions were widely distributed across geographic and environmental distances.However, we showed that large differences in multiple ecosystem functions were always associated with high β-diversity between communities (red points in Fig. 2b illustrating that the non-null highest values of differences in multiple ecosystem functions were associated with that of βdiversity).In addition, the relationships between plant and soil microbial β-diversity and spatial turnover in multiple ecosystem functions remained positive and strong when we controlled for the influences of geography, climate, soil pH, and any combinations of the three factors (Fig. 2c).Our results were maintained when we examined the spatial differences of the single ecosystem functions (Appendix S1: Table S3).However, we found that the effects of abiotic factors on the relationship between β-diversity and spatial turnover in multiple ecosystem functions depended on the identity of ecosystem functions considered.For example, when abiotic factors were controlled, plant β-diversity was positively associated with plant aboveground biomass, plant P, and soil extractable P. Soil bacterial βdiversity was positively associated with root biomass and soil total N, while soil fungal βdiversity was only significantly associated with plant aboveground biomass (Appendix S1: Table S3).
Finally, we used SEMs to assess the direct and indirect effects of geographic and environmental distances on spatial turnover in ecosystem function through changes in above-and belowground β-diversity.SEMs indicated that 14% of variance in spatial turnover in multiple ecosystem functions was explained by plant and soil fungal β-diversity (Fig. 3).Plant β-diversity (standardized path coefficient, hereafter β std = 0.28) was more directly and positively associated with spatial turnover in multiple ecosystem functions than was soil fungal β-diversity (β std = 0.12; Fig. 3; Appendix S1: Tables S4 and S5).Meanwhile, plant β-diversity was strongly associated with soil fungal β-diversity (β std = 0.47 in Fig. 3a; β std = 0.62 in Fig. 3b), and soil fungal β-diversity was strongly associated with soil bacterial βdiversity (β std = 0.59 in Fig. 3a; β std = 0.63 in Fig. 3b).The models did not detect significant direct effects of geographic or climatic distances on spatial turnover in multiple ecosystem functions (Fig. 3).

DISCUSSION
We found clear evidence that greater differences in the spatial turnover of multiple ecosystem functions were associated with greater differences in plant and soil microbial β-diversity, but not geographic or environmental distances.The bivariate associations between above-and belowground β-diversity and spatial turnover in multiple ecosystem functions were retained when we controlled for geographic and environmental distances.Moreover, spatial turnover in multiple ecosystem functions was strongly associated with plant β-diversity and was weakly associated with soil fungal β-diversity.We found indirect effects of soil bacterial β-diversity through soil fungal β-diversity on spatial distances.Xs are one of the groups of β-diversity (i.e., plant, bacterial, and fungi).Points represent the Spearman correlation coefficients (rho) and error bars 95% confidence intervals.Details of the statistical summary are given in Appendix S1: Table S3.
(Fig. 2. Continued) v www.esajournals.orgturnover in multiple ecosystem functions.Our results support our first hypothesis that plant and soil microbial β-diversity determine the spatial turnover of multiple ecosystem functions.In addition, although our results support our second hypothesis that abiotic factors mediate the relationship between β-diversity and the spatial turnover of multiple ecosystem functions, we found that the impacts of geographic and environmental distances on spatial turnover in multiple ecosystem functions were mainly through changes in above-and belowground species composition.These novel findings highlight that the legacies of historical events vs. contemporary environmental factors should be incorporated for a more nuanced understanding of the causes of changes in above-and belowground species composition and their consequences for multiple ecosystem functions in naturally assembled communities.

Mechanisms explaining above-and belowground β-diversity effects on spatial turnover in multiple ecosystem functions
Two potential mechanisms explain why βdiversity, particularly of plants, is a stronger predictor of spatial turnover in multiple ecosystem functions than is soil microbial β-diversity.First, plant species differ in abundance (Soliveres et al. 2016a), dominance (Lohbeck et al. 2016), functional traits (Lefcheck and Duffy 2015), and evolutionary history (Cadotte et al. 2017), each of which influences the magnitude of ecosystem functions.The plant communities examined here span several continents, capturing large species pools that are likely to be influenced by a wide range of evolutionary and ecological processes, leading to wider trait dispersion in the plant communities among sites (Seabloom et al. 2015).In turn, the differences in functional traits among plant species are major factors driving spatial differences in multifunctionality (Roger et al. 2016, Gross et al. 2017, Blesh 2018, Grman et al. 2018, Villnäs et al. 2018).We found changes in species abundance and dominance in plant species composition did not influence the bivariate association between plant β-diversity and the spatial turnover of multiple ecosystem functions (see more detailed results in Appendix S1: Table S2).Our results are supported by two empirical studies, which found minor effects of species abundance on spatial variation in multifunctionality in grasslands (Grman et al. 2018) and in subtropical coniferous forests (Li et al. 2021).Our findings therefore support the idea that differences in functional traits and evolutionary history are likely the main mechanisms explaining plant βdiversity effects on the spatial turnover of multiple ecosystem functions.
A second explanation for the pattern we observed is that differences in plant communities are highly associated with differences in composition of other trophic levels (e.g., soil microorganisms in this study), which promote multiple ecosystem functions (Soliveres et al. 2016b, Schuldt et al. 2018, Martinez-Almoyna et al. 2019, Anujan et al. 2021).Previous work from the Nutrient Network using the same grassland sites as this study found that decreases in plant βdiversity were associated with the decreases in soil microbial β-diversity (Prober et al. 2015) and the extent of multifunctionality (αmultifunctionality; Hautier et al. 2018).These findings are in broad agreement with the results presented here.Thus, we propose that taxonomic covariance, or linkages among plant and soil microbial communities, contribute to differences in multiple ecosystem functions.For instance, our work showed that soil fungal β-diversity indirectly affected the spatial turnover of multiple ecosystem functions through changes in plant β-diversity or vice versa (Fig. 3).Although we detected no direct effects of soil bacterial βdiversity on the spatial turnover of multiple ecosystem functions, soil bacterial β-diversity was significantly associated with fungal βdiversity, which indicates that bacterial βdiversity could indirectly affect the spatial turnover of multiple ecosystem functions through changes in soil fungal β-diversity.
(Fig. 3. Continued) (SRMR).***bootstrap P < 0.001 and **P < 0.01.Note the highly significant goodness of fit for the second model (P = 0.002) suggests that the first model is likely the better model of the two.Details of the statistical summary are given in Appendix S1: Tables S4 and S5.

Indirect effects of abiotic factors on spatial turnover in multiple ecosystem functions
The weak associations between geographic or environmental distances and spatial turnover in multiple ecosystem functions probably indicates that other key abiotic factors are important but were not measured in this study and thus cannot be explored (e.g., soil parent material, soil types, soil heterogeneity, and vegetation structure (Grman et al. 2018, Hu et al. 2020, Song et al. 2020)).However, in our study, the most likely explanation for the weak associations is that abiotic factors have indirect rather than direct effects on spatial turnover in multiple ecosystem functions, perhaps through their effects on community composition (Allan et al. 2015, Grman et al. 2018).We used SEMs to test this idea and found evidence that geographic and climatic distances were not directly associated with the spatial turnover of multiple ecosystem functions (Fig. 3).Our results are partially supported by an elevational study, which found that climate had indirect effects on spatial turnover in multiple ecosystem functions, while soil properties had both direct and indirect effects on spatial turnover in multiple ecosystem functions (Martinez-Almoyna et al. 2019).Despite our inability to detect direct effects of abiotic factors, no studies that we are aware of have examined the abiotic effects of geographic and environmental distances combining the biotic effects of above-and belowground β-diversity at continental scales.
Our findings highlight the potential roles of legacies of historical events and contemporary environmental factors on above-and belowground community assemblages.That is, plant and soil microbial assemblages differed in different locations and were not randomly distributed among the grassland sites.Although abiotic factors had strong effects on plant and soil microbial β-diversity, soil bacterial and fungal β-diversity are more driven by geographic distance and differences in climate and soil pH than was plant βdiversity.This suggests that the spatial distributions of soil microorganisms are more driven by past (e.g., dispersal limitations and past climatic conditions) and contemporary environmental conditions (e.g., climate and soil pH; Ladau et al. 2018, Nottingham et al. 2018) than is aboveground plants.However, our findings indicate that as grasslands face global threats from shrub encroachment, conversion to agriculture, or are lost to development, the functions and services they provide will be more convergent with an increase in aboveground community homogenization than with soil microbial community homogenization.In addition, we found soil bacterial β-diversity was more driven by climate and soil pH than was soil fungal β-diversity.These results suggest that different organismal types differ in the processes of shaping community composition assemblages.Our work therefore supports a unified research framework to study the mechanisms underlying the biogeography of all life using the knowledge of macroecology (Shade et al. 2018).Furthermore, we found that different organismal types in above-and belowground had different influences on individual ecosystem functions.This finding highlights the importance of biodiversity in maintaining multiple ecosystem functions and reinforces the role of spatial differences in above-and belowground species compositions as a mechanism of shaping spatial differences in multiple ecosystem functions.

CONCLUSIONS
Our work differs from classical experiments in the study of biodiversity and ecosystem function relationships (Fukami et al. 2001, Loreau and Hector 2001, Hector et al. 2009, Jochum et al. 2020).We are not asking whether biodiversity is related to an increase or decrease in a single ecosystem function, nor are we asking whether biodiversity is related to multiple ecosystem functions (Byrnes et al. 2014).Instead, we examine what controls similarity in multiple ecosystem functions.Although there is evidence that soil fungal (Mori et al. 2018) and soil multitrophic (Martinez-Almoyna et al. 2019) βdiversity are broadly linked to similarity in multiple ecosystem functions, our work provides evidence that changes in aboveground plant species composition are more associated with similarities in multiple ecosystem functions in grasslands around the world than with changes in soil microbial composition.Importantly, the effects of abiotic factors are largely indirect.These indirect effects suggest that the impacts of past and ongoing climatic and edaphic change on ecosystem functions are mediated through their strong v www.esajournals.orginfluences on plant and soil microbial communities.Those indirect effects, coupled with the ongoing homogenization in above-and belowground community composition, have substantial effects on the diversity of ecosystem functions performed by grasslands.Therefore, our work is consistent with the idea that spatial variability in species composition contributes to multiple ecosystem functions (Mori et al. 2018, Winfree et al. 2018) and may be one mechanism that interacts with environment and soil properties across spatial scales.

Fig. 2 .
Fig. 2. Relationships between abiotic and biotic factors with spatial turnover in ecosystem function.Associations of spatial turnover in ecosystem function with abiotic factors (geographic, climatic, and pH distances) (a) and β-diversity (b).Spearman correlation coefficients (rho) and the lower and upper 95% confidence intervals of simple Mantel tests are shown.Lines represent the trends of the bivariate associations.Points in red represent spatial turnover in multiple ecosystem functions that are significantly higher than null expectation, in blue lower than null expectation, and in gray no difference from null expectation.(c) Partial Mantel tests for the bivariate associations between β-diversity and spatial turnover in ecosystem function given geographic and environmental