Statistical analysis
Firstly, the effects of N addition, drought, delaying snowmelt and their interactions on plant community features (plant biomass, plant richness, plant foliar N), environmental factors (plant height, CVheight and soil moisture) and the diversity and abundance of total arthropod and each functional group were analyzed using repeated-measures ANOVA, with plot and year as random factors. The species diversity of arthropod communities was estimated using cumulative effective diversity index, calculated aseH , where H is the Shannon-Wiener index of arthropod community. Compared with Shannon-Wiener index, cumulative effective diversity could reduce errors by controlling the difference in species richness caused by differences in sampled insect individuals (Ricklefs &Miller, 2000).
Next, we used Mantel test to analyze the correlation of two paired arthropod functional groups and build interaction networks. Plant community, environmental features and each arthropod functional group are represented as a matrix. Dependent parameters mentioned above has been partitioned into 5 parts: (1) plant biomass; (2) the presence/absence of plant species in the matrix (coded as 1 and 0 respectively); (3) N content of plant species under all the treatments; (4) soil moisture, plant height and CVheight displayed as a matrix; (5) each arthropod functional group. Bray-Curtis dissimilarity was calculated for each plant community character, environmental features and arthropod functional groups across all plots of a same treatment. Then, Mantel tests were used in turn for all pairs of functional groups from different trophic levels and from the same trophic level, respectively. We also calculated the relationship between environmental features and each arthropod functional group by Mantel tests, and the link between the two groups existed whenever r (correlation index) was significant (p < 0.05). In the condition of vertical interactions, the indirect biotic effects of non-adjacent trophic levels are mediated by correlated functional groups, so we did not analyse the relationship between non-adjacent trophic levels. Considering the modest effect of delaying snowmelt in the peak of growing season, the network webs were plotted for ambient, N addition, drought and their interaction for two years.
To characterize the interaction network, we calculated three indices, including connectance (CFG), interaction diversity (ID) and interaction strength (IS) (Rzanny & Voigt, 2012). These indices can be used to predict community stability and ecosystem functioning and were found to be sensitive to multiple global change factors in interaction networks (Tylianakis et al., 2008; Tylianakis et al., 2010; Ebeling et al., 2011). IS and ID were calculated by the r value of above-mentioned Mantel tests. CFG and ID were explained as different aspects of ecosystem complexity, which is the diversity of interactions among different functional groups. We took omnivores and detritivores as herbivores when calculating CFG, ID and IS. The validity of so doing lied in the following facts: (1) Most of omnivorous beetles feed on plant stems, roots and litters and the rest omnivores, mainly ants, were significantly correlated with plant chewers (r 2=0.46, p <0.01); (2) Detritivore species (mainly Diptera) feed mainly at larval stage on plant litter, stage, which was significantly correlated with plant communities (biomass: r 2 = 0.52,p <0.01; richness: r 2 = 0.28,p <0.01; N content: r 2 = -0.17, 0.05<p <0.1).
CFG was calculated as the number of realized links between two arthropod functional groups in each pair (L) divided by the number of potential interactions (LMax).
CFG = L / LMax
ID was calculated using the Shannon index: wherepi is the proportion of interaction i to the total sum of n observed interactions (i.e. each significant Mantel test was divided by the total number of all significant Mantel tests):
ID = \(\sum_{i=0}^{n}p\)ilog2 pi
Mean IS is the arithmetic mean of number of all significant interactions:
IS = \(\frac{1}{n}\sum_{i=0}^{n}p\)i
In addition, the IS of environmental factors with all the arthropod functional groups was assessed and Duncan’s test was used to analyze the significance of difference of IS between all the experimental treatments including ambient.
Finally, we explored the driving factors of the effect of plant communities on arthropod network in β-diversity level. The variation of plant richness (including species turnover and nest) and biomass (including balanced-variation and abundance-gradients components) could strongly impact arthropod networks. We followed the additive taxonomic β-diversity partitioning method on plant absence-presence matrix (Baselga, 2010) and plant biomass matrix (Baselga & Chao, 2017). In plant biomass matrix, the total β-diversity (βbc) has been partitioned into βbal(balanced-variation) and βgra (abundance-gradients components). In plant absence-presence matrix, the total β-diversity (βsor) has been partitioned into βsim(species turnover) and βsne (nestedness). To compare the direction and magnitude of responses to experimental treatments, we calculated the response ratio [ln (RR)] and the 95% confidence interval of plant β-diversity under N addition, drought and their interactions.
All statistical analyses were performed using the statistical software R 1.2.5033 by ‘vegan’ (version 2.5-6), ‘betapart’ (version 1.5-4) package.