Statistical analyses
Contrast between urban and rural soil, plant traits and AMF variables: To test for differences between urban and rural environments in soil and plant traits (question 1), and AMF variables (AMF-colonization rates, spore density, and diversity; question 2), we ran independent one-way ANOVAs to test for the effect of the environment factor (three levels: RS, DUS, OUS). Variables such as soil properties (e.g. N concentration, pH), spore density, richness, and diversity were pooled at the site level; spore diversity was calculated using the Shannon-Wiener diversity index. On the other hand, because replication for variables such as plant morphometrics and AMF-colonization rates were at the plant level, we implemented the lmer function from package lmer4 in R (Bates et al. , 2015) to perform linear mixed models considering environment and site as fixed and random factors, respectively. The significance of the fixed effect was evaluated using a Type II Wald’s test running the Anova function from the car package in R (Fox & Weisberg, 2019), assessing the significance of the random effect with a likelihood-ratio test (LRT). To discard potential spatial autocorrelation, Mantel’s tests were used to assess the correlation between geographic distance and Euclidean distance for each variable recorded using the package ade4 in R (Chessel et al. , 2004); however, in none of the considered dependent variables spatial autocorrelation was detected. To test for changes in community composition based on spore morphospecies between environments, we conducted a permutational analysis of variance (PERMANOVA) with 999 permutations using the adonis function from the vegan package in R (Oksanen et al. , 2020). Finally, a Principal Component Analysis (PCA; prcomp , R), based on a correlation matrix, was used to test for differences in the multidimensional variation of soil properties and plant attributes between each environment (question 1). A one-way ANOVA tested theenvironment (RS, DUS, OUS) effect on PC1. This PC1 summarized the variation on soil nutrient concentration (N, P, K, Na, and Ca) and pH. Before PCA, we equalized P and Ca variances by dividing raw values by 100.
Soil and AMF associations: To explore the soil properties that may predict AMF-colonization rates on roots of R. nudiflora(question 3), we performed pairwise correlations between AMF-colonization rate and soil attributes at the global level (i.e. using data from all environments) as well as for each environment separately. As a next step, we used the Z-scores of PC1 (see previous section) to explore the association between PC1 and AMF-colonization rate (using pooled data). However, because PCA visual inspection and previous statistical analyses on the first principal component (PC1; see previous section) indicated strong differences in Z-scores between urban and rural environments we ran an independent linear model (i.e. AMF-colonization rate ~ PC1) within each environment.
Structural equation modeling: Because analyses based on principal components obscure the relationship between soil attributes mediating the levels of AMF-colonization rate, we implemented Structural Equation Modelling (SEM) using lavaan package in R (Rossel, 2012), and tested for contrasts between urban (DUS, OUS) and rural (RS) environments using a multigroup test. We used SEM approach to disentangle the direct and indirect effects that soil nutrients and pH have on AMF-colonization rate (question 3). We designed an initial causal model (Fig. S2) based on previous knowledge of soil properties and nutrient interactions (Fig. S2). We considered Ca as one of the main drivers of pH, which in turn can affect K, P, and N availability (Osman, 2013). However, the two later macronutrients (P and N) have been reported to negatively affect K availability (Osman, 2013; but see Dibb & Thompson, 1985). In turn, it is expected that N and P will have important and negative effects on the AMF-colonization rate (Salvioli di Fossalunga & Novero, 2019; Chen et al. , 2020). Finally, Na can also increase the benefit of the association with AMF in the context of salinity stress (Evelin et al. , 2009), and, in soils with low K availability, AMF can play an important role to improve plant K nutrition (Garcia & Zimmermann, 2014). The goodness-of-fit of the initial causal and alternative models (Fig. S2) was assessed with an LRT to test the null hypothesis that the predicted covariance matrix of the model is not different from the observed (Iriondo et al. , 2003). A significant χ2 indicates that the model does not fit the observed covariance matrix. Because a good-fit model may result from an inadequate statistical power (Mitchell, 1992), we calculated additional indices of goodness-of-fit, such as the Comparative Fit Index (CFI; cut-off good fit ≥ 0.95), the Root Square Mean Error of Approximation (RMSEA; cut-off good fit ≤ 0.06), and the Standardized Root Mean Square Residual (SRMR; cut-off good fit ≤ 0.08), which are insensitive to sample size (Hooper et al., 2008). All these indices were estimated using the fitMeasures function in the lavaan package in R (Rossel, 2012), and were used in the process of model selection.
Once we obtained a fitted and general causal model, we ran a multigroup test to assess whether path coefficients contrast between DUS, OUS, and RS environments. The multigroup test imposes cross-group constraints on the path model regression coefficients, and then compared with an LRT, the constrained and unconstrained models using the lavTestLRTfunction (lavaan package in R; Rossel, 2012). In particular, we used the function lavTestScore, a Lagrange Multiplier Test (LMT), as a guide to identify a set of constrained path coefficients that, if released, would result in a significantly better model (i.e. a lower goodness-of-fit χ2 ; Bentler, 1989).
All statistical analyses were performed using R version 4.0.3 (R Core Team, 2020). Furthermore, all analyses, were evaluated with transformed variables when needed to meet the normality assumption.