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