Climate drives the spatial distribution of mycorrhizal host plants in terrestrial ecosystems

Mycorrhizal associations have massive impacts on ecosystem functioning, but the mode and magnitude heavily depend on the mycorrhizal type involved. Different types of mycorrhizas are recognized to predominate under different environmental conditions. However, the respective importance of climate and soil characteristics in shaping mycorrhizal global distributions are still poorly understood. We provide a quantitative and comprehensive global analysis of the main climatic and edaphic predictors of the distribution of plants featuring different mycorrhizal types. Estimates on per grid‐cell relative above‐ground biomass of plants holding arbuscular mycorrhiza (AM), ectomycorrhiza (EcM) and ericoid mycorrhiza (ErM) association were related to a set of 39 climatic and edaphic variables. We assessed their relationship by applying a Generalized Additive Models for Location, Scale and Shape (GAMLSS). The best GAMLSS models were able to explain 55%, 41% and 46% of the variance in AM, EcM and ErM distribution, respectively. Temperature‐related factors were the main predictors of distribution patterns for the three different mycorrhizal plant types. AM plants are favoured by warm climates, while EcM plants’ dominance (and to some extent ErM plants too) is favoured by colder climates. Synthesis. The observed lack of importance of soil drivers challenges the predominant view that mycorrhizal plants distribution mainly reflects soil type preferences—as related to its nutrient foraging strategies—of the different mycorrhizal types. Instead, our results highlight climate—and particularly temperature—as the main force shaping the distribution of arbuscular mycorrhiza, ectomycorrhiza and ericoid mycorrhiza host plants at the global scale and suggest that climate change can significantly alter the distribution of mycorrhizal host plants, with a subsequent impact on ecosystem functioning.

According to differences in morphology and plant and fungal taxa, seven major types of mycorrhizas are distinguished (Smith & Read, 2008). Among these types, arbuscular mycorrhiza (AM), ectomycorrhiza (EcM) and ericoid mycorrhiza (ErM) are the most taxonomically and geographically widespread, being present in the majority of terrestrial biomes. It has been estimated that approximately 80% of the Earth's plant species form mycorrhizal associations with AM, EcM and ErM fungi (Brundrett & Tedersoo, 2018).
The majority of plant species is able to form mycorrhizal symbiosis of only one type (Wang & Qiu, 2006), with only a few exceptions in which the same plant species can be colonized by two mycorrhizal fungi types (McGuire et al., 2008).
AM, EcM and ErM associations predominate under distinct edaphic and climatic conditions. This differentiation is presumed to be strongly associated to the different nutrient uptake strategies among AM, EcM and ErM fungi. For example, EcM and ErM fungi are capable of breaking down organic matter through the expression of extracellular lytic enzymes, making these associations more suitable for organic soils (Read, Leake, & Perez-Moreno, 2004). In contrast, AM saprotrophic abilities are less developed, causing AM to mostly rely on inorganic compounds as a source of nutrients, and therefore more prevalent in mineral soils (Smith & Smith, 2011). Based on these insights, Read (1991) and Read and Perez-Moreno (2003) proposed a theoretical model where the abundance of AM, EcM and ErM host plants gradually changes along a latitudinal and altitudinal gradient, driven mainly by the effects of climate on decomposition, which is ultimately reflected in the accumulation of organic C in the soil and the availability of nutrients for plants. According to this model, AM plants dominate in grasslands and tropical forests; EcM trees are abundant in temperate and boreal forests; and, finally, plants featuring ErM associations predominate in heathlands.
Since Read's first approach, only a few attempts have been made to understand quantitatively which environmental drivers explain the distribution of distinct types of mycorrhizal plants. Menzel et al. (2016) focused on AM and analysed the geographical distribution and environmental drivers of AM plants status (obligate, facultative or non-mycorrhizal) on a regional scale (Germany). Bueno et al. (2017) examined how the number of plant species featuring distinct mycorrhizal traits (type and status) varied with different climatic and soil factors at the European scale. Only recently, Steidinger et al. (2019) performed a coarse resolution (1 degree) global analysis on mycorrhizal trees distribution and its environmental drivers although focusing specifically on forest ecosystems. Despite these efforts, the contribution of the different driving forces (e.g. dispersal, climatic factors, edaphic characteristics or evolution) in shaping the biogeography of mycorrhizal vegetation of the entire plethora of plant functional types at global scale and covering all natural biomes and plant growth forms needs better understanding. Moreover, most of the previous studies were based on the number of plant species capable of forming different mycorrhizal associations, without taking the relative abundance of these species in the ecosystems into account.
A quantitative understanding of the relationships between environmental drivers and the relative abundance, in terms of biomass or plant cover, of AM, EcM and ErM host plants is important, because the relative abundance of mycorrhizal types largely underpins ecosystem functioning. Changes in relative abundance of the different mycorrhizal plant types lead to changes in C and nutrient cycling (Phillips et al., 2013;Soudzilovskaia, Van Der Heijden, et al., 2015), soil processes and structure (Rillig & Mummey, 2006), and can even cause deeper modifications in plant community assembly (Van Der Heijden, 2002). In an era of human-induced environmental changes, unravelling the relative importance of soil and climatic factors in shaping the geographical distribution of plant species featuring different mycorrhizal types will lead to better predictions of changes in ecosystem functioning under a future climate.
Here, we present the first quantitative global analysis of the role of climatic and edaphic factors in explaining the distribution patterns of the three main types of mycorrhizal plants that cover all natural biomes and includes all plant growth forms. Our analysis is based on a high-resolution gridded dataset (10 arc-minutes), which includes information about 39 environmental variables and the percentage of above-ground biomass of plant species featuring AM, EcM and ErM mycorrhizal associations. Following Read's hypothesis, we expect a relatively high contribution of soil properties related to organic C content.

| Distribution of biomass fractions of different mycorrhizal associations
Estimates on the relative above-ground biomass of AM, EcM and ErM mycorrhizal associations were obtained from the high-resolu- For the purpose of this paper, non-natural biomes (croplands and urban areas) and bare areas were excluded from the analysis to ensure reliability. This exclusion was performed using the 2015 Land Cover Initiative map developed by the European Space Agency at 300m spatial resolution (https ://www.esa-landc over-cci.org/ ) as a reference. As a result, a total of 270,353 gridded cells were included in the final dataset.

| Climatic and edaphic factors
We assembled a dataset of climatic and edaphic variables that have been proposed to be potential drives of mycorrhizal plant distribution at global scale (Read, 1991;Smith & Read, 2008). In total, our dataset includes information about 39 environmental variables (see Tables S1 and S3). The inclusion of this large number of variables allowed us to evaluate the contribution of temperature, precipitation, seasonality and soil physico-chemical properties to shaping the global distribution of different mycorrhizal plant types.
Climatic variables were obtained from the WorldClim database, Version2 (http://world clim.org/version2; Fick & Hijmans, 2017) at 10 arc-minutes resolution. In total 19 bioclimatic variables were included (see Table S1). These bioclimatic variables are a combination of monthly temperatures and precipitation values. The inclusion of the 19 bioclimatic variables allowed us to determine potential correlations with seasonality or extreme and limiting environmental factors. In addition, Annual Global Potential Evapotranspiration (Global-PET) (https ://cgiar csi.commu nity/categ ory/data/; Zomer et al., 2007;Zomer, Trabucco, Bossio, & Verchot, 2008) was added to the climatic variables due to its ecological relevance. Global-PET was calculated according to the Hargreaves equation (Hargreaves, Hargreaves, & Riley, 1985) which includes mean temperature, daily temperature range and extra-terrestrial radiation.
Data on the main edaphic variables were obtained from the Harmonized World Soil Database (HWSD) (http://dare.iiasa.ac.at/; FAO/IIASA/ISRIC/ISS-CAS/JRC, 2012). We included in total 12 variables (see Table S2) from the soil top layer (0-30 cm), which were scaled up to 10 arc-minutes resolution using the mean of the raster cells as aggregation criterion.
Data on water-holding capacity, Total C, Total nitrogen (N), Total phosphorus (P) and available P is not available in the HWSD database. We considered these variables to have a potential implication on mycorrhizal host plant distribution due to their high ecological relevance, and therefore we prioritized their inclusion.

| Statistical analysis
As climatic variables are highly correlated (Table S3), we applied a principal component analysis (PCA) to alleviate the problematics related to the high degree of collinearity while maintaining a high degree of variance in climate variables. The first two axes (PC1 and PC2) of the PCA explained 79.6% of the total variance in climatic data. PC1 was mainly related to temperature variables; while PC2 incorporated mainly precipitation-related variables ( Figure S1). Soil factors were examined individually due to the low explanatory power of the principal components and difficulties with the ecological interpretation of the PCA axes of the soil variables (see Figure   S2).
Generalized Additive Models for Location, Scale and Shape (GAMLSS) were fitted to relate the percentage of biomass of AM, EcM and ErM plants, respectively, to the soil factors and PC1 and PC2 of the climatic factors using the "gamlss" package. A GAMLSS allows fitting flexible regression and smoothing models and relaxes the assumption of the exponential family distribution for the response variable, replacing it by a general distribution family. Models were fitted using a zero-inflated beta distribution, which is appropriate for modelling proportional data that contain a high proportion of zeros. The smooth functions of each predictor were restricted to a maximum of 3 degrees of freedom, allowing for nonlinearity while detecting only general trends and avoiding overfitting issues. Assuming that different mycorrhizal plant types may vary independently of environmental drivers, EcM, AM and ErM plant distributions were modelled separately. For model simplification, interaction terms were not included.
Model selection was performed by testing competing models that included a set of variables within which each variable explained at least 5% of the data variance, had a Pearson pairwise correlations lower than 0.6 (see Table S4) and variance inflation factors lower than 3. This procedure allowed us to select for sets of non-correlated variables with high explanatory power and to avoid including suppressive variables that would obscure the interpretation of the models. In total, we tested 18 different competing models for AM plant distribution, each of which included eight different variables, six competing models for EcM plant distribution (each including six different variables) and two competing models for ErM plant distribution (each including three variables) (see Tables S5, S6 and S7). For each mycorrhizal plant type, the best model was selected according to the lowest Bayesian information criterion.
After the best models have been selected, a further variable selection was performed. We removed non-significant variables (with p > .05) and variables with low relative importance in the model. We considered that a variable had little explanatory power when the effect of removing the variable did not decrease the Pseudo R 2 (Nagelkerke, 1991) with more than 1%. Finally, degrees of freedom of the smooth terms were reduced to preserve only clearly nonlinear patterns.

The presence of spatial autocorrelation (SAC) in AM, EcM and
ErM final model residuals was tested using Moran's I correlograms with the "sp.correlogram" function in the "spded" package.
Moran's tests confirmed the presence of SAC in the model residuals. The existence of SAC may lead to an overestimation of degrees of freedom and Type I errors may be strongly inflated (Legendre, 1993). The presence of SAC can be alleviated by (a) Including spatial coordinates explicitly in the model as covariates: This can be problematic since they could covary with the environmental variables present in the model (Dormann, 2007;Miller, Franklin, & Aspinall, 2007), which can obscure the interpretation of the relative importance of the predictors. (b) Accounting for SAC in model residuals: There is a wide range of methods available in the mainstream software that allow alleviating SAC in model residuals (Dormann et al., 2007). However, their implementation in the context of a zero-inflated beta distribution is still extremely limited.
This problem is even increased by the large number of data points included (270,353), which makes the computation of the spatial models unfeasible.
Due to these technical limitations, no correction of SAC could be applied to our global high-resolution data. However, filtering the dataset by distances where SAC is significantly reduced as they decrease exponentially with distance (see Figure S3) demonstrated that the presence of SAC does not alter the importance of the predictors in the final models and therefore their interpretation is not biased due to the autocorrelation (more detailed information about the reduced models is provided in the Supporting information). As the main goal of the models is to detect important predictors of mycorrhizal plant distribution and not to serve as a predictive tool, we further discuss the output of the model with the complete dataset.
The final models were validated by 10-fold cross-validation. A difference of less than 10% between the RMSE (root mean squared error) of the final models and cross-validated models was used as a criterion for model validity. Both in AM, EcM and ErM models, the difference was lower than 5%.
Statistical analysis was performed using r 3.5.3 (R Core Team, 2018) and gridded data was processed using ArcGis v10.2.2.

| RE SULTS
The model selection applied to the AM host plant distribution retained in total two different climatic and soil predictors: temperaturerelated factors (PC1), and bulk density. Together, these predictors were able to explain 55% of the variance in AM plant distribution (as indicated by Pseudo-R 2 ). PC1 was, by far, the best single predictor, F I G U R E 1 Predicted relation between AM (a), EcM (b) and ErM (c) relative abundances and the environmental factors maintained in the best models. Each relation was calculated setting the rest of the variables to the mean value. Light coloured shades represent the region within the upper and lower 95% confidence limits. Numbers between brackets in the x-axes correspond to the individual variance explained by each factor in the models  Table 1).
The difference between the sum of Pseudo-R 2 of each variable (0.46) and the Pseudo-R 2 of the final model (0.55) indicates that 9% of the variance explained is shared between the two predictors.
For the relative abundance of EcM plants, the predictors retained by the best model were temperature-related factors (PC1) and base saturation. This set of predictors explained 41% of the total variance (   (Davidson & Janssens, 2006;Doetterl et al., 2015). In line with this, our dataset shows that, at the global scale, the principal components of climatic factors and soil properties are not highly correlated (see EcM and ErM plants being a reflection of their differential ability to take nutrient from organic sources (Read, 1991;Read & Perez-Moreno, 2003) is not supported by our findings. Our results also partially contradict the conclusion drawn by Steidinger et al. (2019), who as well found a strong climatic control over mycorrhizal trees distribution. Steidinger et al. (2019) related the mechanisms explaining this pattern purely to differences in decomposition rates, while they did not find a direct link with soil physicochemical properties. Our results suggest that other mechanisms play a role, as detailed below.

| Environmental predictors of AM plant distribution
Our results clearly highlight the impact of climate (especially temperature) on AM plant distributions. Several studies have reported temperature as an important limiting factor for the growth of AM extraradical mycelium (Gavito, Schweiger, & Jakobsen, 2003;Heinemeyer & Fitter, 2004;Rillig, Wright, Shaw, & Field, 2002). Also, a reduction of intraradical colonization has been commonly reported at temperatures lower than 15°C (Gavito & Azcón-Aguilar, 2012;Hetrick & Bloom, 1984). As an alternative mechanism, Veresoglou (2019) recently proposed that irradiance reduction in higher latitudes contributes to a reduction of AM fungi responsiveness, which may contribute to the detected decline of AM plant abundance in colder climates. In line with these studies, our findings suggest that the physiological restrictions of AM fungi to develop and provide benefits to its plant partner at lower temperatures might be a primarily important driver of AM plant distribution at the global scale, independent of soil properties.
In contrast, soil properties were not relevant in explaining AM abundances (Table 1). Especially surprising is the absence of F I G U R E 2 Predicted global distribution of AM (a), EcM (b) and ErM (c) mycorrhizal host plants and prediction residuals (d-f); here only the 5% of data points with the highest residual values are depicted. Light blue areas denote non-natural biomes, bare areas or regions for which no environmental data was available. Residues are expressed as the difference between predicted and observed AM, EcM, and ErM plant relative abundances. Red points (positive values) indicate zones where the predicted plant relative abundance was overestimated by the model and blue points (negative values) indicate underestimations soil P impacts in the final AM best model, which contradicts the view of AM associations being a key adaptation for P uptake. This view was already challenged by previous research. For instance, Soudzilovskaia, Douma, et al. (2015) reported no significant correlation between P limitation and AM root colonization. Similarly, using a meta-analysis approach, Allison and Goldberg (2002) showed that changes in P availability do not have a consistent effect on mycorrhizal infection at plant community level. These results indicate that, although P availability influences the performance of the plant-fungi relationship at the plant species level (Treseder, 2013), this does not necessarily translate into P availability driving AM distribution patterns at a global scale.
What is clear from these results is that climatic conditions are deeply affecting the global biogeography of AM associations. . Recent research also suggests that the ability of certain AM fungal species to colonize leaf litter may contribute to a higher abundance of this association in organic soils (Bunn, Simpson, Bullington, Lekberg, & Janos, 2019). Incorporating information about specific fungal functional traits and host identities will be key in future studies aimed to better understand AM plant biogeographical patterns.

| Environmental predictors of EcM plant distribution
The relative abundance of EcM plants was mainly explained by temperature-related factors, but showed trends opposite to those of AM. EcM plants showed preferences for moderately cold climates, which is consistent with their greater abundance in Northern temperate and boreal zones (Soudzilovskaia, Van Bodegom, et al., 2019).  Fougère-Danezan, & Bruneau, 2017;Tedersoo, ). These species are suggested to proliferate in nutrient-poor and acidic soils (Campbell, 1996) where specific traits of ectomycorrhizal fungal communities (e.g. the ability to obtain N from organic sources) may give them advantage over AM associations (Alexander & Högberg, 1986;Högberg, 1986

| Environmental predictors of ErM plant distribution
The distribution of ErM plants has been traditionally associated with harsh environments, characterized by nutrient-poor and acidic soils (Read, 1991). This has been related to the ability of ErM fungi to produce hydrolytic and oxidative enzymes (Cairney & Burke, 1998)

| CON CLUDING REMARK S
Our results point at temperature-related factors as the main predictors-instead of soil properties-for the global distribution of the three most abundant mycorrhizal plant types. The observed lack of importance of soil drivers contradicts the traditional view of climate-driven soil properties, such as the rate of organic matter decomposition and nutrient availability as the ultimate mechanisms explaining the latitudinal distribution of mycorrhizal plant types (Phillips et al., 2013;Read & Perez-Moreno, 2003;Smith & Read, 2008;Steidinger et al., 2019).
In contrast, our findings support the role of temperature as a main driving force affecting the global distribution of plant ecological strategies (Moles et al., 2014), and reinforces the view that mycorrhizal type constitutes an important part of these strategies. We suggest that the latitudinal transition between AM, EcM and ErM plants is higher quality soil data and (c) to increase the knowledge of mycorrhizal associations in plant species that have not been investigated yet to extend the analysis beyond the dominant species. This will allow to account for the large heterogeneity of soil properties and to evaluate the importance of smaller-scale processes that cannot be taken into account in this work.