Results

Overall fish diversity

We recorded a total of 29 benthic species during 252 surveys over the three-year study. Raw sample species richness ranged from 0 to 13 and average richness at sites was 4.2. There were two occasions that we recorded no fish and these events were discarded in the analysis. The inter- and intra- annual differences in richness were not significant (p = 0.411 and 0.542 for monthly and annually richness, t-test), however, the richness difference between habitat types was significant with more species recorded at natural sites (Carex sedges and open waters) than at plantation sites (t-test, p = 0.002 and <0.001 forCarex sedges and open waters, respectively). For abundance, although there was no significant inter-annual difference (p = 0.452, t-test), the catch was significantly higher in May than other months (p = 0.028, t-test), coincident with the breeding season. In contrast to richness, fish abundance at Carex sedge sites was significantly higher than that at plantation sites (P <0.001, t-test) and open water sites (P = 0.044, t-test), which had comparable abundance (p = 0.200, t-test).

Temporal community variation and its drivers

Synchrony, variance ratio and stability varied greatly across habitat types and between sampling years (Table 1). The total β-diversity showed great spatial and temporal variations (Table 1). Specifically, the temporal variation of benthic fish communities was significantly lower in the flooding year (i.e. 2017) than in normal years (i.e. 2018 and 2019) (p < 0.01, t-test). In plantations, the total temporal β-diversity was significantly lower than open waters andCarex sedges (p < 0.01, t-test). Although the difference in temporal β-diversity was significant (p = 0.04, t-test) in the flooding year, it became non-significant in normal years. The best model fitted for temporal β-diversity included water total phosphorus concentration (TP) and habitat type. This model explained 68% of the variations of β-diversity across the sampling sites (Table 2). The estimated model coefficient indicated that temporal β-diversity decreased significantly among the TP gradients, and was significantly lower in plantations than in the natural habitat of Carex sedges and open waters.
The observed synchrony of fish abundance ranged from 0.09 to 0.90, and in 10 occasions out of the total of 36, the synchrony was significantly greater than null expectations. These results suggested that most taxa abundances were changing in the same direction, especially in the modified habitat (Table 1).
We did not detect any spatial or temporal patterns in synchrony variations (Table 2). The best model describing the synchrony variation included a single environmental variable (i.e. TN, Table 2), which explained 38% of the variation in the spatial and temporal variation of community synchrony (Table 2). Total nitrogen had a significantly (at 0.05 level) negative effects on synchrony (the lower 2.5% and upper 97.5% intervals do not include zero). Other tested environmental variables such as water depth, TP, and water clearance had no statistically significant relationship to species synchrony.
The calculated VR ranged from 0.30 – 2.83 with the majority cases greater than one (Table 1), of which seven were significantly higher than the null VR. Similar to species synchrony, these results indicated that the benthic fish community largely positively covaried, especially at the plantation sites. As for variation ratio, there was no clear spatial or temporal patterns. In addition, no model with environmental variables was better than the null based on the loo model selection procedure.
Community stability had a weakly negative relationship with TN. The relationship was not significantly at 0.05 level (Table 2), but was at 0.1 level (the lower 5% and upper 95% interval of the coefficient for TN was -0.12 to -0.01, data not shown). Like synchrony and VR, other tested environmental variables had no significant impacts on community stability.

Biotic mechanisms of community temporal dynamics

Species synchrony and the variance ratio were significant predictors for community stability (53% and 66% of the total variation was explained by the VR and synchrony model, respectively, Table 3). In comparison, both α and β diversity models had considerably lower explanatory power with a mean Bayes-R2 of 0.33 (0.30, 0.54) and 0.23 (0.06, 0.40) (Table 3). This lower performance was also reflected in the relatively larger confidence interval of the response curve, especially when abundance stability was high (Fig. 3). While the modelled relationship with VR and β diversity was linear, stability had non-linear relationships with synchrony and α diversity (Table 3 and Fig. 3). Community stability decreased with both VR and synchrony, however, the marginal effect curve of the stability – α diversity was hump-shaped unimodal, suggesting lower stability at both the extreme ends. In contrast, the relationship with β-diversity was linear and community abundance stability increased significantly with β-diversity (Table 3, Fig. 3).
Habitat type was also included in these models. The estimated model coefficients suggested that the plantation sites had significantly higher community stability than the sedge sites (the 95% interval did not include zero for all models, Table 3). However, the difference between plantation and open water was not significant when the effects of species variance rate, synchrony or α-diversity were accounted for (zero was within the bracket of the 95% confidence interval, Table 3).

Ecological Stochasticity and its effects on community dynamics

There were no inter- and intra- annual patterns in spatial ecological stochasticity in benthic fish community. However, there were strong spatial patterns: the modified sites had significantly lower stochasticity than the more natural sites (Fig. 4). The modelled coefficient indicated that ecological stochasticity was the predominant process at natural sites (mean ecological stochasticity greater than 0.50, Fig. 4), whereas deterministic processes were prevailing at the modified sites (mean ecological stochasticity was 0.37, Fig. 4). In addition, none of the tested environmental variables had a significant effect on ecological stochasticity and the best regression model included only habitat types as explanation variable. The model had a relatively high Bayse-R 2 of 0.42 (95% confidence interval was 0.20 - 0.57).
The temporal ecological stochasticity was a strong predictor of all investigated metrics of community temporal dynamics (the mean Bayes-R 2 and its 95% confidence interval were 0.41 [0.22, 0.55], 0.30 [0.08, 0.50], 0.59 [0.42, 0.69] and 0.67 [0.55, 076] for stability, variance rate, synchrony, and total temporal β-diversity respectively, Table 4). While the normalized ecological stochasticity had significantly positive effects on species variance rate and synchrony, it had significantly negative effects on community stability (Fig. 5). The effects of stochasticity on temporal β-diversity was highly non-linear: the β-diversity increased with stochasticity at low level but stayed more or less stationary when the NST reached ~ 0.4.
The community stability was significantly higher in the plantation habitat than in the sedge habitat while the difference between plantation and open water was not significant (Table 4). Moreover, community stability was significantly lower in the flooding year of 2017 than in the normal water years of 2018 and 2019. For species synchrony and variance rate, the difference between habitat types was not significant (Table 4). Regarding the yearly variation, species synchrony was significantly lower in flooding year of 2017, and the between-year difference for species variance rate was not significant (Table 4).