Statistical analysis
We fitted generalized linear mixed-effects models (GLMMs) for each mycorrhizal type to test for the effects of road disturbance (hypothesis 1), temperature and elevation, and their interactions with road disturbance (hypothesis 2) on the respective mycorrhizal type covers (2822 data points from 894 individual plots). We used beta regressions with the glmmTMB package (Brooks et al. 2017) in R (R Core Team 2021) after transformation of the response variable (i.e. proportion data) to avoid extreme values of 0 and 1: (response variable value * (number of observations – 1) + 0.5) / number of observations (Cribari-Neto and Zeileis 2010). The explanatory variables used were: (i) distance to the road as a proxy for disturbance, as a three-level factor for each plot (0 to 2 m from the road, 2 to 52 m from the road and 52 to a 102 m from the road); (ii) mean annual soil temperature; (iii) elevation; as well as (iv) the two-way interactions between these two factors. The elevation values used were relative to each region’s elevation gradient obtained by scaling elevation individually for each region using the scale function in base R (Becker 2021) resulting in gradients between -1 and 1 for all regions, with the lowest elevation of each gradient being the valley where the roads start, i.e. the point at which the elevation gradient starts and not sea level. We chose this because we were interested in the elevational distance relative to the bottom of the gradient and not in the absolute elevation of a region, as the latter is not easily comparable. The initial model contained both elevation and temperature, as elevation can serve as a measure of non-climate driven and more local gradients while temperature takes into account the differences between regions as well as large-scale climate-driven trends within a region. However, after testing for multicollinearity using the VIF (Variable Inflation Factor) through the vif function in R (Fox and Weisberg 2018), elevation and temperature were found to be too strongly correlated (VIF value of 5.812 for elevation). We consequently omitted the effect of elevation from the final model. However, to make sure that temperature and elevation patterns did indeed behave similarly we also ran the model selection strategy in parallel with elevation instead of temperature. The random intercept term of transect nested in road nested in region was added to consider the hierarchical nature of our design, as well as a random intercept term for year of observation to consider repeated surveys, and a random slope term for plots. Candidate models with all possible combinations of fixed effects were then derived from the complete model and compared using AICc (Akaike Information Criterion, corrected for small sample sizes). Only models with a ΔAICc of less than 2 units compared to the best candidate model were retained (Zuur et al. 2009). We then applied a variation partitioning approach to the selected models using the performance package in R to determine the proportion of variation in mycorrhizal type proportions explained by disturbance, mean annual soil temperature and elevation (Lüdecke et al. 2021).
To investigate the variation between regions, we used a partial pooling approach (Harrison et al. 2018). This was done using models derived from the aforementioned beta regressions for each mycorrhizal types with the same explanatory variables that were retained after model selection, i.e. temperature and disturbance, but with additional random slope terms for both of these terms. The addition of these random effect allows for the intercept and slope associated with each region and each variable to deviate while still capturing the overall trends from the larger dataset. We then extracted the coefficients associated with each variable for every region, allowing for comparisons of trends between these regions and the global models.
A similar modeling approach was then used to investigate how the proportion of plants associated with the different mycorrhizal types correlated with the proportion of non-native plant cover (hypothesis 3). We ran separate tests for the percentage of total non-native plant cover in the roadside plots (0-2 m from the road) and in the neighboring vegetation in the furthest adjacent plots (52-102 m from the road). The 2-52 m plots were left out for this particular analysis, because we know from previous studies (Clavel et al., 2020; McDougall et al., 2018) that using the 2-52 m adjacent plot could be misleading for non-native species as it in some cases still included roadside vegetation when the roadside was more than 2 m wide. As the presence of non-native species is linked to changes in the local balance of mycorrhizal association types, we used the percentage cover of vegetation associated to each given mycorrhizal type amongst native species only as a predictor instead of the proportion of total vegetation cover. The models included native cover percentage of a certain mycorrhizal type and mean annual soil temperature as well as their interaction as explanatory variables and either roadside plot or adjacent plot non-native plant cover percentage as a response variable in order to distinguish between non-native species simply benefitting from the road disturbance and more established non-native species present in the surrounding vegetation. As previously, a random intercept of road nested into region was included to account for the survey’s hierarchical design. Model selection was then performed by comparing AICc values as described above.