Statistical analyses
To summarize the climatic variables (AP, AT, DP, MAXT, MINT, SP, ST and
WP) and soil properties (EC, pH, Soil C, Soil N and SWC), we employed
the “rda” function from the R package “vegan” (v. 2.5.7; Oksanenet al. , 2020) to conduct two principal component analyses (PCAs):
one for climate and one for edaphic factors. For the climate PCA, PC1
was negatively correlated with AP, AT, DP, MINT and WP and positively
correlated with MAXT, SP and ST; in the soil PCA, PC1 was negatively
correlated with pH and positively correlated with EC, Soil C, Soil N and
SWC. We then calculated Spearman rank-order correlations among the
climatic, geographic, plant community, plant defense and soil variables.
Implementing the “diversity” function for each plot, we also
calculated plant species richness (Richness), the Shannon diversity
index (Shannon) and Pielou’s evenness index (Pielou). To summarize plant
community composition, a detrended correspondence analysis (DCA) was
used (implemented using the “decorana” function), given the potential
for arch effects in plant community composition data (Borcard et
al. , 2011).
We evaluated how altitude, latitude and longitude affected CWM
herbivory, intraspecific variability and turnover using a series of
linear mixed-effect models; these were implemented using the “lmer”
function, with “site” as a random effect, and the corresponding
t-values, P-values and marginal R2 values were
calculated. Various climatic, geographical, plant community and soil
factors were also used as independent variables in linear mixed-effect
models to test their effects on CWM herbivory, intraspecific variability
and turnover. We then compared these models with a null model
(i.e. , an intercept-only model) using Akaike’s information
criterion corrected for small sample sizes (AICc);
comparisons were performed using the AICc function in
the “AICcmodavg” package (v. 2.3.1; Mazerolle, 2020). The
log-likelihood (LL) was calculated for each model using the “logLik”
function and AICc-based parameters, including the change
in AICc relative to the top-ranked model
(ΔAICc), AICc weight
(w AICc) and the marginal R2(Burnham et al. , 2011). Linear mixed-effect models (with “site”
as a random effect) were also used to test the effects of “better”
predictors (i.e. , predictors better than the null model as
determined using AICc) on CWM herbivory, intraspecific
variability and turnover.
We then built a piecewise structural equation model (piecewiseSEM) to
test how the climatic, plant community and soil variables affected CWM
herbivory through intraspecific variability and turnover effects. In the
SEM, standardized path coefficients (scaled by their mean and standard
deviation) were calculated, as well as the corresponding significance
(P-values) for each path; the overall fit of both full and final models
was determined using Fisher’s C statistic and AICc using
the “piecewiseSEM” package (v. 1.3.2; Lefcheck 2016).
For population-level herbivory, we tested how altitude, latitude and
longitude affected herbivory in Polygonum macrophyllum ,Saussurea pulchra and Taraxacum mongolicum using linear
mixed-effect models with “site” as a random effect. Linear
mixed-effect models were also used to predict population-level herbivory
using various climatic, geographic, plant community, plant defense and
soil factors; AICc, ΔAICc,w AICc and the marginal R2served as indicators of model goodness of fit. All statistical analyses
were performed in R v4.1.1 (R Development Core Team, 2021).