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).