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.