Data analysis
Multiple generalized linear models (GLMs), using R package lme4with the “glm” function (Bates et al.,
2015), were employed to analyze effects of light intensity, nutrient
loading, snail biomass and their interactions on primary producers
(including P. crispus shoot biomass, shoot number and internode
length, V. spinulosa total biomass, periphyton Chl a and
phytoplankton Chl a content, C, N, P content and their ratios). Light
intensity and nutrient loading were treated as categorical variables,
and snail biomass as a continuous variable. Variance inflation factors
were calculated to test for collinearity problems that might arise from
including interaction terms in the same model, and no collinearity
problems were detected in all the models. Due to the non-linear response
of P. crispus shoot biomass and shoot number to snail biomass
under full light, the quadratic of snail biomass (poly function) was
used as a predicted variable to analyze these effects, and the nonlinear
responses were further confirmed by generalized additive models using
“gam” function (Table S1). Since no significant effects of nutrient
loading were observed on the growth of macrophytes and periphyton (Table
1), nutrient loading and the interaction with other factors were
excluded and analyzed again in the GLMs (Table S2) and plotted. Multiple
linear mixed-effects models were applied to test effects of the
treatments on water temperature, dissolved oxygen, conductivity, TN and
TP, with sampling date as a random factor (Table S3). Estimated marginal
means were compared after each linear model test to compare the
difference between the treatments, using package emmeans. QQplot
and residual plot were used to visually assess the normality of data. If
the data were not normally distributed, data were transformed (data
transformation is added in Table 1 and S2). One sample of P.
crispus (in the E1S1 treatment) was missing after harvesting. This
resulted in 39 P. crispussamples being available for the
analysis of the three plant elemental compositions
(plant C, N, and P content), and
three plant stoichiometric traits
(C: N, C:P and N: P ratios).
A structural equation model (SEM) was constructed to summarize the
effects of shading, snail biomass, and nutrient loading on biomass of
the macrophytes. This allowed assessing the complete graphical network
of the interactions and relationships, with the directions of paths in
the SEM diagram indicating causal influences
(Rosseel, 2012). Quadratic of snail
biomass was used to predict the response of P. crispus biomass.
Three indices of model fit were used with conventional significance
thresholds to assess the overall fit of the SEM, with the
χ2 p -value (p > 0.05), the
standardized root mean squared residual (SRMR ≤ 0.08), and the
comparative fit index (CFI ≥ 0.95) (Hu &
Bentler, 1999). All SEM procedures were conducted with thelavaan (version 0.6-3) package in R
(Rosseel, 2012). All of the analyses were
performed using R software (ver. 4.0.1; R Development Core Team, 2021).