Three common submerged macrophytes and two common aquatic animals, Hydrilla verticillata, Vallisneria natans, Potamogeton maackianus, Cyprinus carpio, and Cipangopaludina cahayensis were chosen to evaluate the effects of species diversity within and across trophic levels on EMF. H.verticillata, V. natans, and P. maackianus are common macrophytes in the middle-lower reach of the Yangtze River, and they often can grow together and form a stable aquatic plant community. They were collected from Shahu Lake, a eutrophic lake near Hubei University, the proximity facilitating transport. The plants were then washed with running tap water to remove particles and other adjacent organisms. Cyprinus carpio that were 7.78 ± 0.80 cm (n = 24) in length and 4.93 ± 1.95 g (n = 24) in weight, and Cipangopaludina cahayensis were 5.68 ± 0.40 g (n = 12) in weight, were brought from Chongqing Yongchuan Yuxi aquaculture base. The protocols involving animals were reviewed and approved by the Animal Ethics and Welfare Committee of Hubei University. Before the experiment, the three macrophytes and the two aquatic animals were cultivated in aquariums with tap water in the greenhouse; the temperature varied from 12.5 °C to 30.4 °C and the average temperature was 23.1 ± 6.1 °C during the experiment.

A species-addition experiment, consisting of twelve 110 L polyethylene buckets (diameter = 0.4 m, height = 0.65 m), was carried out from August 20 to November 20, 2021. The buckets were placed in the experimental greenhouse of Hubei University and filled with tap water which had total nitrogen (TN) and total phosphorus (TP) concentrations of 1.8 and 0.02 mg/L, respectively, and a water depth of 47 cm. The buckets were randomly divided into four treatments, each replicated three times, including H. verticillata monoculture (Mon), H. verticillata + V. natans mixed culture (Mix_two), H. verticillata + V. natans + P. maackianus mixed culture (Mix_three), and H. verticillata + V. natans + P. maackianus + C. carpio + C. cahayensis mixed culture (Mix_five). During the experiment, the water temperature varied from 10.7 °C to 28.6 °C, the average water temperature was 21.3 ± 6.0 °C, and the pH of the water varied from 8.3 to 9.9.

Before the experiment, H. verticillata, V. natans, and P. maackianus were washed thoroughly with tap water followed by deionized water, and the lateral branches or stolons of the plants were removed. Then, complete clean plants of V. natans (10 cm length) and apices of H. verticillata and P. maackianus cut at 10 cm below the top were prepared and planted in small plastic pots (length 18 cm and width 12 cm), which were filled with 5 cm fine sand. Fourteen plants were planted in each plastic pot. The density of H. verticillata was 14 in monoculture, the densities of H. verticillata and V. natans were eight and six in Mix_two, and the densities of H. verticillata, V. natans, and P. maackianus were eight, three, and three in Mix_three in each pot, respectively. In total, 168 young seedlings of similar size were transplanted into 12 small plastic pots, which included three H. verticillata monoculture pots, three H. verticillata + V. natans mixed culture pots, and six H. verticillata + V. natans + P. maackianus mixed culture pots. Then, 12 small plastic pots were placed randomly into the 12 polyethylene buckets. At the same time, three polyethylene buckets of the H. verticillata + V. natans + P. maackianus mixed cultures and two C. carpio and three C. cahayensis were randomly added to each polyethylene bucket. We used a line to lift and lower the pots, which allowed gentle and nondestructive sampling. During the experiment, deionized water was replenished every 2 days to compensate for evaporation losses and assure that the water levels were similar in all 12 polyethylene buckets. Also, dead C. carpio were replaced with live ones of approximately similar size.

At the beginning of the experiment, the height/length and fresh weight of all 168 macrophytes were measured with a ruler and an analysis balance. Macrophytes were measured on days 14, 28, 42, 56, and 70 after the start of the experiment. This was done by pulling up the tray on which the plant pots were standing, making handling and measuring easy and non-destructive. The height and length of each branch were measured with a ruler. After each measurement (except the last time), we put the macrophytes back into the corresponding polyethylene bucket. At the end of the experiment, all macrophytes were harvested and weighed, and the biomass of phytoplankton in the water column was proxied by total chlorophyll. Chlorophyll content was measured according to the method described by Jampeetong and Brix (2009). Filamentous macroalgae growing on the macrophytes or at the water surface were collected by hand, and epiphyton growing on the bucket wall and on the sand surface was collected by brush and then weighed.

The macrophyte biomasses of different sampling times were estimated by the fitted linear relationship between plant height and biomass using the data of height and biomass measured at the beginning of the experiment and at 70 days. The linear fitting equations for the three tested macrophytes were as follows:

where x is the plant height in cm and y is the plant biomass in g.

The water temperature, dissolved oxygen content and pH, total nitrogen (TN), and total phosphorus (TP) concentration of each polyethylene bucket were measured at the beginning of the experiment and on days 14, 28, 42, 56, and 70. The water temperature, dissolved oxygen content, and pH of the water bodies were measured at 30 cm depth between 10:30 and 12:00 using a YSI (YSI 6600V2). The irradiance was measured by an illuminance meter (MODEL ZDS-10W-2D). TP content was measured using a colorimeter, an AutoAnalyzer (Bran+Luebbe GmbH, Inc., Germany), following sulfuric acid/hydrogen peroxide digestion and the ammonium molybdate ascorbic acid method. TN content was measured using an IL-500N nitrogen analyzer (Hach Company, Loveland, USA). During the experiment, different volumes of NH₄CL and KH₂PO₄ were added every 14 days to maintain 1.8 mg/L TN and 0.02 mg/L TP concentration in the water column in each bucket.

Three ecological functions, including total net primary productivity (the total biomass of macrophyte, phytoplankton, filamentous macroalgae, and epiphyton) and cumulative variation of TN and TP, at the end of the experiment in each polyethylene bucket, were used to quantify EMF.

The cumulative variation of TN and TP was calculated as:

where C_{TN ( TP )} represents the cumulative variation of TN or TP, C _i represents the TN or TP content of the water measured at different sampling periods, and C ₀ represents the value of the TN or TP content at the beginning of the experiment.

According to Maestre et al. (2012), we used the mean value method to quantify EMF because it simply and intuitively reflects the ability of communities to maintain EMF (Jing et al., 2015). We calculated the average of standardized data of the total net primary productivity and TN and TP cumulative variation. The multifunctionality index (M) for each polyethylene bucket was the average of the three evaluated functions M can be calculated according to the method used by Xu et al. (2016):

where F represents the number of functions being measured, f_i represents the measures of function i, r_i represents a mathematical function that converts f_i to a positive value, and g represents standardizing of all index values (the top 5% of the measured values of each function are averaged as the maximum value of the function, and the ratio of each measured value to the maximum value is calculated).

The mean value method cannot distinguish between one function being provided at a high level and another being provided at a low level vs. two functions being provided at an intermediate level (Byrnes et al., 2014). As such, we calculated MF _t using 20%, 30%, 40%, 50%, 60%, and 70% of our threshold values. Thus, we supplemented the mean value method with a threshold analysis method for further analysis of the relationship between species richness and EMF.

Community temporal stability of macrophyte productivity was calculated for each treatment using the ratio of μ to σ:

where μ represents mean net macrophyte primary productivity and σ represents net macrophyte primary productivity standard deviation.

Before statistical analysis, all data were tested for normality and were subjected to Levene’s tests. Non-normal data were transformed (log10) to obtain normality. Repeated-measures ANOVA was used to test the effects of species richness, time, and their interaction on the total primary productivity, C_{TN ( TP ),} and temporal stability of macrophyte productivity. A one-way ANOVA was applied to determine the statistical significance (P< 0.05) of the differences in the total primary productivity and the productivity of macrophytes, filamentous macroalgae and epiphyton, and phytoplankton at the end of the experiment, as well as C_{TN ( TP ),} EMF, and temporal stability of macrophyte productivity among different treatments. The relationships between species richness and individual EFs, EMF, and temporal stability of macrophyte productivity were evaluated using simple OLS regressions within and across trophic levels. Variation partitioning analysis was used to analyze the degree of explanation of three indicator systems (environmental factors including water temperature, dissolved oxygen, pH, irradiance, time, and species richness) for the variation in EMF within and across trophic levels.

All the statistical analyses were performed with SPSS 20.0 for Windows (IBM Inc., Chicago, IL, USA) and R version 4.2.1 (https://www.r-project.org/), and figures were generated in Origin 2021 for Windows (OriginLab, Inc., USA).

For the three tested individual EFs (total primary productivity and the cumulative variation of TN and TP) and the temporal stability of macrophyte productivity, repeated-measures ANOVA showed that experimental treatments only significantly influenced total primary productivity (P< 0.001, Table 1 , Figure 1 ). Time significantly influenced total primary productivity and the cumulative variation of TN ( Table 1 , Figure 2 ). The interactions of experimental treatments and time affect the three tested individual EFs ( Table 1 ).

The response patterns of the total primary productivity and the macrophyte productivity for experimental treatments were similar. At the end of the experiment, the highest total primary productivity and macrophyte productivity both occurred in the Mix_two treatment, but significant differences in total primary productivity and macrophyte productivity were only found between the Mon and Mix_two treatments (P< 0.05, Figures 1A, B ). For other primary producers, the productivity of phytoplankton showed no significant difference among the four experimental treatments, and the productivity of filamentous macroalgae and epiphyton under the Mix_three treatment was significantly higher than Mon and Mix_two treatments (P< 0.05, Figures 1C, D ).

For the cumulative variation of TN and TP, the values of the two tested individual EFs in the Mix_five treatment were significantly lower than that under other experimental treatments at the end of the experiment (P< 0.05). However, no significant differences were found in the amount of cumulative variation of the two aforementioned individual EFs among the three experimental treatments within a trophic level ( Figures 2A, B ). For the temporal stability of macrophyte productivity, there was no significant difference among the four experimental treatments ( Figure 2C ).

A one-way ANOVA showed that trophic complexity had a significant effect on EMF (P< 0.05, Figure 3 ). At the end of the experiment, the value of EMF across trophic levels was significantly lower than that within a trophic level (P< 0.05), while there was no significant difference among different macrophyte species richness treatments within a trophic level ( Figure 3 ).

For individual EFs, simple OLS regression analyses showed that species richness positively correlated with the total primary productivity within and across the trophic levels ( Figure 4A ). Across the trophic levels, there was a significant negative relationship between species richness and the cumulative variation of TP and EMF (P< 0.05, Figures 4C, D ). However, there were no significant correlations between species richness and the cumulative variation of TN and temporal stability of macrophyte productivity within and across the trophic levels ( Figures 4B–E ).

A threshold analysis showed that species richness positively affected the number of functions exceeding the threshold at values of 40% and 50%, and the relationships between species richness and the number of functions became flatter at higher and lower threshold values of 70% and 30%, respectively. By contrast, when the threshold value was 20%, species richness negatively affected the number of functions exceeded ( Figure 5 ).

Variation partitioning analyses showed that the number of variances of EMF explained by environmental factors, time, and species richness were similar among, within, and across the trophic levels, and the total explanations for EMF within and across trophic levels both were 99.60%, although the independent effects of the three driving factors were weak ( Figure 6 ). Variation partitioning showed that the shared effects of environmental factors and time accounted for a high percentage, with their combined effects for EMF within and across trophic levels 94.10% and 89.00%, respectively ( Figure 6 ). Furthermore, our results showed that time has the strongest independent effects on EMF within and across trophic levels among the three tested factors.

For individual EFs, as found in several other studies (Tilman et al., 1997; Engelhardt and Ritchie, 2001; Isbell et al., 2011; Chang et al., 2022), we showed that species richness–productivity relationships are positive both within and across trophic species diversity levels. Previous studies showed that the two main mechanisms responsible for the positive species richness–productivity relationship are the niche complementarity effect hypothesis and the selection effect hypothesis (Loreau and Hector, 2001). In this study, our results indicated that positive richness–productivity relationships may be largely due to the selection effect. Some studies demonstrated that higher productivity from greater species richness sometimes arises from the greater chance that superior competitors with high productivity are present in species mixtures (Tilman et al., 1997; Engelhardt and Ritchie, 2001). Our results indicated that H. verticillata became dominant due to environmental selection and competition and provided higher productivity within a community. Similar life-forms of submerged macrophyte species involved in our study may reduce the niche complementarity effect on system productivity. This may be because complementarity among submerged macrophyte species, especially for species of a similar type, might be rather weak, resulting in species diversity having no significant effect on vegetative production (Zhang et al., 2019). Thus, it is important to note that the types of species used in such experiments may lead to different outcomes (Zhang et al., 2019). Besides macrophyte productivity, our results showed that there were higher filamentous macroalgae and epiphyton productivity in communities with greater species richness. The most plausible explanation may be through a niche complementarity effect because various macrophytes can access different and more overall resources (Vanderstukken et al., 2011; Zhang et al., 2019). Therefore, our results suggest that the underlying mechanisms of species richness–productivity relationships were primary producer taxonomic group-specific.

Different from the positive correlation between species richness and productivity, our findings showed that species richness did not significantly relate to the cumulative variation of TN and TP within a trophic level. Previous studies demonstrated that BEF relationships for different EFs are not always congruent, and there are many different BEF relationships with positive, negative, nonlinear, and non-significant correlations (Mittelbach et al., 2001; Hooper et al., 2005; Duffy, 2009; Jaillard et al., 2014; Brose and Hillebrand, 2016). The observed non-significant correlations between species richness and the cumulative variation of TN and TP can be explained by other important processes besides nutrient uptake by plants and algae, including nutrient recycling by microorganisms or nutrients retained by sediment or surface structures, playing vital roles in nitrogen and phosphorus removal (Engelhardt and Ritchie, 2001; Zhang et al., 2019). We found a significant negative relationship between species richness and the cumulative variation of TP across trophic levels, but there were no significant correlations between species richness and the cumulative variation of TN. The main reason may be due to the different removal mechanisms of TN and TP in the water column, e.g., TN loss is mainly via denitrification by microorganisms; however, TP loss is mainly by plant absorption (Kyambadde et al., 2004; Paerl et al., 2011; Ma et al., 2017; Xing et al., 2020). Another reason may be that consumer excretion and decomposition of dead bodies release nutrients to the water column, resulting in increased phosphorus concentrations in the water column (Northcote, 1988). Our results further suggest that consumers play an important role in the EFs of freshwater systems (Duffy, 2002; Thébault and Loreau, 2003; Duffy et al., 2005) and that multitrophic species diversity has a greater impact on EMF than macrophyte diversity (Lefcheck et al., 2015).

Quite different from most research results that biodiversity enhances EMF (Hector and Bagchi, 2007; Lefcheck et al., 2015; Perkins et al., 2015; Moi et al., 2021), our results showed a significant negative correlation between species richness and EMF across trophic levels. Early research about the relationship between biodiversity and EF showed similar negative species richness–EMF relationships and suggested that greater diversity caused by poor competitors results in reduced productivity (Loreau, 1998; Lasky et al., 2014). In our study, the negative species richness–EMF relationships were not due to primary productivity because there were non-significant differences in total primary productivity. Two reasons account for this phenomenon. First, it may be due to selected individual EFs and the functions used to calculate EMF (Edlinger et al., 2020; Garland et al., 2021; Jing and He, 2021), because different individual EFs would influence EMF through different trade-offs or synergies (Zavaleta et al., 2010; Byrnes et al., 2014). Second, the average approach used to calculate EMF in our study may affect the objective weight and general weight of the three selected individual EFs, which can lead to higher contribution rates of TN and TP cumulative variation for EMF (Byrnes et al., 2014). The effects of nutrients released via consumer excretion and the decomposition of dead bodies may be additive (or synergistic). Thus, our results emphasize the importance of selecting individual EFs, as well as the fact that the calculation of EMF can result in different outcomes for BEMF relationships.

Increased diversity facilitates the temporal stability of community productivity (Cottingham et al., 2001; Valone and Hoffman, 2003; Tilman et al., 2006). Although we hypothesized that there were positive diversity–ecosystem stability relationships within and across trophic levels, the lack of consistent differences in temporal stability of productivity among the four experimental treatments did not support this hypothesis. In agreement with some other studies (Jiang and Pu, 2009; Sasaki and Lauenroth, 2011; Ma et al., 2017), our results also showed that increased macrophyte species richness does not affect the temporal stability of macrophyte productivity within a trophic level. The most probable explanation is that productivity was controlled by dominant macrophyte species, rather than species richness, which results in no effect of species richness on macrophyte productivity (Sasaki and Lauenroth, 2011; Ma et al., 2017). Another possible reason may be lower species richness and lower trophic level complexity, as the relationship between species richness and temporal stability of community productivity was critically dependent (Jiang and Pu, 2009). Our study uses low species richness and fewer trophic levels and, therefore, almost certainly underestimates the levels of diversity needed for EMF.

Our results showed that the links between EMF and biodiversity within and across trophic levels are mediated by the combination of environmental factors and time, more so than by species richness. Many studies have proved that environmental factors and time can affect the interactions between the various elements of an ecosystem which results in effects on the EMF (Bai et al., 2022). Changing environmental factors can affect EF both directly via altering rates of ecosystem processes (Spaak et al., 2017) and indirectly via altering species dynamics and interactions (Dornelas et al., 2014; Blowes et al., 2019); environmental factors could also interact with biodiversity in regulating EF (Eisenhauer et al., 2019; Benkwitt et al., 2020). For example, individual responses to environmental variations cause changes in population abundance, generating interactions between populations, which in turn affect community structure and composition that extend to entire ecosystems (Hofmann and Todgham, 2010). Some studies indicated the effects of environmental factors on EMF may change over time (Tilman et al., 2001; Xu et al., 2021; Guo et al., 2022; Hong et al., 2022), and the response rate of EF to environmental change depends on the resilience of ecosystem functions over time (Wu et al., 2018; Qu et al., 2019). Long-term experiments are revealing that the slope of the BEF relationship can vary with environmental context across sites and through time (Isbell et al., 2011; Qiu and Cardinale, 2020). Thus, some researchers strongly emphasize that the dimension of time should be taken into consideration in experiments and observational studies which can predict whether, when, and to what extent ecosystems will be affected by altered conditions (Edlinger et al., 2020).

We found a weak influence of species richness on EMF within and across trophic levels. This also may be due to selected individual EFs and the diversity of the taxonomic group considered. In our study, we only took macrophytes, snails, and fish as representatives for biodiversity. Other groups of organisms, such as microbes, can drive many ecosystem processes, especially nutrient cycling, but our study does not directly address the functional role played by these and other important groups. Obviously, our results underestimated the effects of some organisms on EMF, especially for microbes which can significantly influence the nitrogen cycle in our study. Therefore, our results further highlight the significance of selecting individual EFs and considering the diversity of the taxonomic groups for BEMF relationship research (Vaughn, 2010; Zhang et al., 2019).