Material and methods

Study Site

Dongting Lake (111°40’–113°10′ E, 28°38’–29°45’ N) is China’s second largest freshwater lake. Due to land reclamation, the whole lake has been divided into three sub-lakes (East, South, and West Dongting Lake), which are hydrologically connected through main river channels. Our study site, the West Dongting Lake (111°57′ - 112°17′ E, 28°47′ - 29°07′ N) is the most upstream section (Fig. 1). Water from the Yuan River, the Li River and two escaping channels of Yangtze River flows through West Dongting into South Dongting, and then into East Dongting, finally discharges into the Yangtze (Fig. 1). The prevailing climate in the area is continental subtropical monsoon, with four distinct seasons featuring a warm and humid summer and a cold and dry winter. The annual rainfall ranges from 1,200 to 1,700 mm, of which up to 60% occurs in summer (Guan et al., 2014).
The West Dongting Lake (WDTL) is a seasonally dynamic system, with large inter- and intra-annual variations in both lake area and water level. In a typical year, the water level increases from the end of March and reaches the highest level in July or August. The water level begins to decrease in September and the lowest level occurs from December to March in the next year. During the summer wet season, water level reaches an average level of 34.30 m (long-term mean, based 1965-2014 records) at Nanzue hydrological station (Fig. 1), and the majority of the area is under water. During the winter dry season, the lake water level decreases to 28.03 m (long-term mean, based 1965-2014 records). As water levels recede, the Lake divides into four broad habitat types: open water, Carex sedge, Phragmites communis marsh andPopulus nigra plantation. While the first two types of habitat are in relatively natural conditions, the other two, especially thePopulus nigra plantation, are subjected to great degree of artificial modifications.
From late 1980s, plantation of Populus nigra , an introduced fast-growing tree for wood pulp, started in WDTL. To ensure young trees survive the summer flooding, elevated ridges were built with sediments excavated from the nearby lakebed, leaving a network of artificial ditches with depth varies from 1.0 to 2.5 m (Li et al, 2020). The plantation creates a novel habitat, which has reduced lateral hydrological connectivity, flattened hydrograph (i.e. the average peak water level is greatly reduced), and stable water level during the low water period. In addition, the artificial ridges present a significant barrier for fish movement even during high-water seasons. Another notable change is that the excavation destroyed submerged plants and sedges, resulting in a habitat with less structural complexity.

Fish survey

During the period of 2017 – 2019, we surveyed benthic fish at 12 locations across the WDTL, of which 6 were from ditches at Populus nigra plantations, 3 were from Carex sedges and 3 were from open waters (Fig. 1). The year 2017 is considered as a flooding year with maximum water level reaching 36.57 m at Xiaohezuo (Fig. 2), which was nearly 3 m higher than the long-term mean maximum water level, while 2018 and 2019 is normal with maximum water level of 33.63 m. For each year, we collected monthly samples from April to October, which covers a completed cycle of water level rise and recession (Fig. 2). In April, the water level is normally high enough to inundate most of theCarex sedges and Phragmites communis marshes. Fish move to the newly flooded areas to breed and forage, taking average of the abundant food sources (King et al. 2003). October is near the end of water recession phase, and water level is high enough to maintain the hydrological connection between channels and Carex sedges.
Fish samples were collected using long fyke nets (length: 30 m, width: 20 cm, mesh size: 10 mm). The net is designed by local fisherman, and is efficient for benthic fish sampling in areas with complex habitat structure (e.g. macrophyte with different height, irregular topography). For each sampling events, the net was set at 9:00 -10:00 am and the catch was collected at 9:00 – 10:00 the next day. All fish were transferred to laboratory for identification (to species level) and counting.

Environmental variables

Before setting the net, temperature, PH and conductivity (SPC) were measured in site using a YSI multiprobe (YSI, Yellow Springs, Ohio, USA), clearance was also measured using a 20 cm white Secchi disc. We also collected instantaneous water grab samples to measure total phosphorus (TP), total nitrogen (TN), total suspended solids (TSS), and chlorophyll a (chla). Vertically integrated water samples (sub-samples were taken from the surface, 0.5 m, 1 m and 1.5 m) were collected and preserved following USEPA methods (Kopp et al. 1983). Water samples were analyzed in laboratory pursuant to APHA protocols (APHA. 2005).

Assessing ecological stochasticity

We calculated both spatial and temporal ecological stochasticity to explore the mechanisms shaping β-diversity patterns in WDTL. We used spatial ecological stochasticity (Ning et al. 2019) to compare the deterministic and stochastic processes shaping benthic fish communities at different habitat types. To quantify the spatial ecological stochasticity, we aggerated all sample sites for each survey occasions within a habitat type to form a metacommunity, assuming that fish can move freely between sampling locations. We used temporal stochasticity to assess the relative importance of deterministic and stochastic processes on community dynamics (see below). To calculate temporal stochasticity, for each site, the monthly community matrices in every sampling year were considered as a metacommunity. Following Chase et al. (2009), we use null models to determine the ecological stochasticity in shaping the benthic fish community. Briefly, the community dissimilarity (we used the abundance-weighted Bray-Curtis index) is compared with the null expectations to quantify ecological stochasticity under different situations (Zhou et al. 2014). In this study, the null communities were generated by randomizing the observed community structure 1,000 times based on the SIM2 algorithm in Gotelli (2000), i.e. species occurrence totals are fixed but all sites within a habitat type are equiprobable. The normalized stochasticity ratio (NST), was then developed with 50% as the threshold between more deterministic (<50%) and more stochastic (>50%) assembly. NST ranges from 0 to 1, where the minimum value is 0 when extremely deterministic assembly drives communities completely the same or totally different, and the maximum value is 1 when completely stochastic assembly makes communities the same as null expectation. The detailed mathematic formulas can be found in Ning et al. (2019).

Quantify the magnitude and variability of temporal fish community dynamics

We calculate four metrics to describe patterns of the fish community change: total temporal β-diversity, stability, synchrony, and variation ratio. These metrics were calculated for each year at every sampling site using the abundance-based community matrix.

Temporal β-diversity

We calculated the total β-diversity (BDTotal) as the total variance in the community matrix (Legendre & De Caceres, 2013) using the abundance based Sørensen index with the beta.div function in the “adespatial” package (Dray et al. 2016).

Stability

Community temporal stability was calculated as µ/σ, using the time series of fish catch data collected for each sampling sites, where µ is the average aggregate abundances of a site across all months and σ is the temporal standard deviation in the fish catch of a site (Tilman 1999).

Variance ratio

Variance ratio (VR) characterizes species covariance patterns in a community (Schluter 1984), which compares the variance of the community (C ) as a whole relative to the sum of the individual species abundance (Si ) variances:
\(VR=\ \frac{\sigma^{2}(C)}{\sum_{i=1}^{N}{\sigma^{2}(S_{i})}}\) (1)
and
\(\sigma^{2}(C)=\sum_{i=1}^{N}{\sigma^{2}\left(S_{i}\right)+2\times\sum_{j<i}^{N}{\sigma^{2}(S_{i}{,S}_{j})}}\),N is the number of species in a community.
If species vary independently (i.e. they do not covary), then the variance ratio will be close to 1. A variance ratio <1 indicates predominately negative species covariance, whereas a variance ratio >1 indicates that species generally positively covary (Schluter 1984; Hallett et al. 2014). We tested if the observed variation ratio significantly differs from 1 using null model approach. By randomly selecting different starting points of time series for each species, we generated 10,000 null communities in which species abundances vary independently but within-species auto correlation is maintained (Hallett et al. 2014). The significance is confirmed if the observed variation ratio falls out of the 1st and 3rd percentile of the null distribution.
The calculation of community stability and variation ratio, and significance tests were conducted using the “codyn” package (Hallett et al. 2016) in R (R Core Team 2019).

Synchrony

Community synchrony provides another measure of species covariance. In study, we defined community synchrony as the degree to which the abundance of fish species in a community rise and fall together through time. Using the community time series of abundance, we calculated community synchrony using the metric of Loreau and de Mazancourt (2008):
\(\varphi=\frac{\delta^{2}[\sum_{i=1}^{N}{S_{i}(t)}]}{{(\sum_{i=1}^{N}{\delta[S_{i}\left(t\right)]})}^{2}}\)(2)
where \(\delta^{2}\) is the variation operator; and\(S_{i}\left(t\right)\) is the abundance of species i at sampling time t . The numerator represents the temporal variance of the aggregate community-level abundance; and denominator is the sum of the population variances squared. \(\varphi\) ranges from 0 (perfect asynchrony) to 1 (perfect synchrony) (Loreau and de Mazancourt 2008).
We tested the significance of the observed community synchrony by comparing it with the distribution of synchrony of 10,000 random communities. The random communities were generated using a Monte Carlo randomization procedure, which shuffles the columns of the community matrix independently. Community synchrony calculation and test were conducted using the “synchrony” package (Gouhier and Guichard 2014).

Statistical analysis

To determine if and how community temporal dynamics is affected by environmental factors, we use generalized linear mixture models (GLMM) or generalized additive mixed models (GAMM), to count the possible non-linear relationship (Wood 2006). Starting from the simplest null model with no explanatory variable, we increased the model complexity by adding sequentially independent variables (i.e. TP, TN, TN:TP, Depth, habitat type and sampling year) as either simple or smooth terms. Because the relatively small sample size, we restricted the number of independent variables to three in any single model. For each community dynamic metric, a total of 126 models were fitted. We select and report the best model using leave-one-out cross-validation (loo), which is a fully Bayesian model selection procedure, similar to the widely applicable weighted Akaike’s information criterion (Vehtari et al., 2017). Moreover, we also investigated the relationships between stability and diversity, VR and synchrony in the same way to explore the effects of biological interactions on community stability (Hallett et al. 2014). To explore the effects of species diversity on community stability, we aggregated all monthly fish catches for each site and calculated the inverse Simpson’s diversity index as community diversity.
The posterior distributions of model parameters were estimated using Markov chain Monte Carlo (MCMC) methods with NUTS sampler implemented in the R package “brms” (version 2.1.0, Bürkner, 2017) by constructing four chains of 100,000 steps, including 50,000-step warm-up periods, so a total of 200,000 post-warm up draws (i.e. (100,000-50,000)×4=200,000) are retained to estimate posterior distributions. All modeling procedures are conducted in a Bayesian framework using the programme Stan (Carpenter et al. 2017) implemented in the R package “brms” (Bürkner 2017).
We used the same Bayesian approach to explore the temporal and spatial patterns of NST, and to link NST with environmental factors.