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.