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.