Seasonal differences in the relationship between biodiversity and ecosystem functioning in an overexploited shelf sea ecosystem

Biodiversity loss results in reduced ecosystem functioning. However, the relationship between biodiversity and ecosystem functioning (the BEF relationship) in overexploited ecosystems remains controversial largely because of the lack of validation in nature. Therefore, in this study, through four seasonal trawl surveys of demersal fish assemblages off the shelf sea of Shandong Province, China, we explored BEF relationships and their seasonal variability. Furthermore, we assessed the relative importance of the environment and biodiversity on biomass variations.


| INTRODUC TI ON
Unprecedented biodiversity loss resulting from both anthropogenic activities (Ceballos et al., 2017;Maxwell et al., 2016) and global climate change (Hillebrand et al., 2018;Teagle & Smale, 2018) has been spreading worldwide, and is known to be detrimental to ecosystem functioning and services Worm et al., 2006). To date, biodiversity loss has been reported in both terrestrial and aquatic ecosystems, and involves many organisms, such as birds, ants, and marine fish (Luypaert et al., 2020;Philpott et al., 2008). Such loss threatens ecosystem resilience, functioning, and services, thereby affecting the protein supply of over a billion people (MEA, 2005;Truchy et al., 2015). Therefore, understanding the effects of biodiversity loss on ecosystem functioning is not only crucial for ecosystem conservation and management but also for human well-being (Naeem et al., 2016;Yasuhara et al., 2016). In this context, the relationship between biodiversity and ecosystem functioning (the BEF relationship) has garnered increasing attention (Baldrighi et al., 2017;Cardinale et al., 2006;Teixeira et al., 2019;Tilman et al., 2014).
Although an increasing number of recent studies have assessed the significance of biodiversity for ecosystem functioning, the mechanisms driving the BEF relationships in natural ecosystems remain controversial (Baldrighi et al., 2017;Séguin et al., 2014).
Previous studies have already reported positive, negative, and irrelevant BEF relationships (Daam et al., 2019). Inconsistent patterns of BEF relationships have even been noted for the same biodiversity indices Maureaud et al., 2019). For instance, species richness, one of the most frequently measured diversity indices in BEF studies, often positively affects ecosystem functioning (Gamfeldt et al., 2014;Loreau & Hector, 2001). However, its negative effects have also been described in other studies (Mora et al., 2014). Furthermore, Maureaud et al. (2019) observed that ecosystem functioning (represented by fish biomass) was related to species evenness and environmental conditions but not to species richness. Overall, we have a limited understanding of the relationship between various biodiversity indices and ecosystem functioning in overexploited shelf sea ecosystems (Strong et al., 2016). Specifically, it remains unclear whether BEF relationships are stable across different seasons because marine ecosystems represent seasonally dynamic environments (Bellino et al., 2019). This knowledge gap regarding BEF relationships has impeded efforts to conserve and manage these overexploited ecosystems (Luypaert et al., 2020).
Therefore, for the sustainable development of the fisheries industry, it is imperative to understand the effects of biodiversity loss on shelf sea ecosystem functioning.
The shelf sea of Shandong Province is located in the central and northern Yellow Sea of China (35°-38°30′N and 118°-124°E) ( Figure 1). This region is particularly important as a key contributor to ecosystem services and functions, in terms of both economic and ecological values. It has rich fishery resources and serves as spawning and feeding grounds for many important commercial species (SPOFD, 2010). In different seasons, each fishing ground shows different water temperatures and other environmental variables, and each spawning ground is characterized by a unique ecology (Wan & Jiang, 2000). Consequently, some species exhibit seasonal variations (Shan et al., 2012). Over the past several decades, the ecosystem structure in this area has been greatly altered as a result of changes in the species composition and a decline in the biomass of commercially important species, due to intensive fishing (Jin et al., 2006;Pang et al., 2018;Tang, 2009). The typical dominant species, including the small yellow croaker (Larimichthys polyactis) and largehead hairtail (Trichiurus lepturus), have been replaced by fish of lower trophic levels, such as anchovy (Engraulis japonicus) and chub mackerel (Scomber japonicus) (Li et al., 2015). This has significantly decreased the biodiversity of the region (Liu & Ning, 2011). Furthermore, under the dual pressures of climate change and anthropogenic activities, it remains unclear whether biodiversity is closely linked to ecosystem functioning and whether this relationship shows seasonal variability. Therefore, this study aims to fill this knowledge gap.
To investigate the BEF relationships, seasonal fishery resource surveys were conducted in the shelf sea of Shandong Province.
We sought to address the following research questions: (i) Which "complementarity" and "selection" mechanisms maintain higher biomass? Presumably, high biomass is often driven by "complementarity" and "selection" mechanisms, which respectively reflect positive species richness-biomass and negative species evenness-biomass relationships (Loreau & Hector, 2001). Therefore, we hypothesized that species richness or species evenness, or both, are significantly correlated to biomass. Functional diversity (FD) quantifies the range and value of organism traits that affect their performance and ecosystem functioning (Violle et al., 2007). In general, higher volumes of niche space occupation (i.e., functional richness [FRic]) and functional redundancy lead to higher ecosystem functioning (Mouillot et al., 2013). Therefore, we hypothesized that positive FRic-biomass relationships may be stable, as stated by Yasuhara et al. (2016).
(ii) Do BEF relationships follow consistent patterns across seasons? Both abiotic and biotic factors shape BEF relationships, and BEF relationships to ecosystem-based management, at least in overexploited shelf sea ecosystems.

K E Y W O R D S
BEF relationship, demersal fish community, environmental factors, overexploited shelf sea ecosystem, seasonal dynamics, taxonomic and functional diversity environmental conditions are closely linked to these relationships (Brose & Hillebrand, 2016;Duncan et al., 2015). Therefore, we hypothesized that with obvious seasonality in species composition or environmental conditions, BEF relationships would follow different patterns among seasons. (iii) What is the relative importance of biodiversity and the environment in explaining the changes in ecosystem functioning? We hypothesized that both the environment and biodiversity play important roles in maintaining ecosystem functioning, as suggested by Belley and Snelgrove (2016).
To the best of our knowledge, this study is the first to consider seasonal shifts in BEF relationships in demersal fish communities in an overexploited shelf sea ecosystem. Our findings provide an important reference for including biodiversity conservation in ecosystem-based management (EBM). trawl survey stations were surveyed in autumn, winter, spring, and summer, respectively. An otter trawl (width, 15 m; mesh size, 46 mm) was towed for 1 h at a speed of 3 knots at each survey station (Zhang et al., 2020). All catches were sorted onboard and identified to the species level in the laboratory. The individuals of each species were F I G U R E 1 Map of the shelf sea of Shandong Province in China. Black circles indicate survey stations counted and their weights were measured. Only demersal fish species were included in the analysis, and all invertebrates and strictly pelagic species were excluded because of their low selectivity with trawling. The biomass (total weight) of the demersal species per station was measured as a proxy for ecosystem functioning. All demersal fish species were further categorized into three thermal groups:

| Trawl surveys and functional trait assessment
warm-, temperate-, and cold-water groups, according to their optimal temperature (Ma et al., 2019).
Nine functional traits were measured to calculate the functional distance between the demersal fish according to their biology, taxonomy, trophic ecology, and life history. Six of these traits were continuous numerical variables, including mean temperature preference, growth, food consumption-to-biomass ratio (Q/B), trophic level, generation time, and length at first maturity. The remaining traits (body shape, swimming mode, and reproductive guild) were categorical traits ( Table 1). These traits are largely uncorrelated to one another (Trindade-Santos et al., 2020). The values of all functional traits were obtained from FishBase (Froese & Pauly, 2019).

| Environmental variables
Two environmental variables, sea bottom temperature (SBT) and sea surface chlorophyll-a concentration (Chl-a) were selected as biotic and abiotic factors, respectively. The SBT dataset was derived from four-dimensional variational ocean reanalysis outputs of the Western North Pacific (FORA-WNP30) for the periods during which the fish data were collected (i.e., autumn and winter in 2016, and spring and summer in 2017) (http://metad ata.diasjp.net/dmm/doc/ FORA_WNP30_JAMST EC_MRI-DIAS-en.html) (Usui et al., 2017).

| Biodiversity indices
In this study, we focused on three biodiversity indices for taxonomic and functional diversity each (TD and FD, respectively). TD was measured based on the species richness, Simpson's evenness, and variance index, representing three facets (richness, evenness, and divergence) of the fish community (Gotelli et al., 2011;Magurran, 2004;Maureaud et al., 2019). The TD indices were computed using the "vegan" package in R (Oksanen et al., 2020).
Functional diversity was measured in a multidimensional functional space, based on the functional traits of species and community composition . First, we generated a trait distance matrix between species based on the Gower distance (Gower, 1971), which assigns equal weights to a mixture of different types of variables (e.g., categorical and continuous). We then generated principal coordinate analysis (PCoA) trait vectors based on the distance matrix. The first three PCoA axes were retained to preserve the quality of the functional space (Maire et al., 2015). Gower distance analysis and PCoA were calculated using the "gowdis" function (Laliberté et al., 2015) and the "pcoa" function (Villéger et al., 2008) in the R package, respectively.
The FD indices were measured using three complementary components: FRic, functional evenness (FEev), and functional divergence (FDiv). The FRic of a community was measured as the volume within the minimum convex hull that contains all species (Villéger et al., 2008).
FEev describes both the regularity of the spacing of species within a multidimensional functional space and the evenness of abundance distribution in functional trait space (Villéger et al., 2008). FDiv is the level of complexity in the combination of species traits and considers abundance distribution within the volume of the functional trait space occupied by the species (Mason et al., 2005). FD was computed using the "dbFD" function in the R package (Laliberté et al., 2015).

| Statistical analyses
We conducted trawl surveys to explore seasonal differences in fish community structure. Using PRIMER 5.0, we generated a Bray-Curtis similarity matrix for fish assemblages between any two seasons. We then used one-way crossed analysis of similarity (ANOSIM) to test the differences in fish community structure among different seasons based on a similarity matrix; values of p < .05 were considered significant. Similarity percentage analysis (SIMPER) was used to identify the contribution of each species to the seasonal differences in demersal fish community structure. One-way analysis of variance (ANOVA) was used to test the effects of seasonal factors on biodiversity indices and biomass.
To investigate the BEF relationships in the shelf sea of Shandong Province, we applied linear models to identify the relationships between the biodiversity indices and biomass. Pearson's correlation analysis was used to test these BEF relationships and values of p < .05 were considered significant. Before analysis, the fish data and biodiversity indices were log 10 (X + 1) transformed to meet the assumptions of multivariate normality and to moderate the effects of outliers. SPSS 22.0 (SPSS Inc.) was used to perform these statistical analyses.
To answer our third question, we first used generalized additive models (GAMs; Hastie & Tibshirani, 1990) to analyse the relationships between environmental variables and biomass (and biodiversity indices). To avoid collinearity and for consistency, the biomass and biodiversity indices for each season were modelled as a smoothing function of a single environmental variable as follows: where C i is the response variable representing the standardized biomass and biodiversity indices of the demersal fish community for each survey station in each season; predictor variable V i refers to environmental variables; s() is the smoothing function; α is the intercept of the model; ɛ is the random error term of the normal distribution with zero mean and finite variance; and N is the number of survey stations in each season. To avoid overfitting, the effective degrees of freedom (edf; representing the level of nonlinearity) were restricted to a maximum of 4.
Furthermore, variation partitioning combined with redundancy analysis (RDA) were used to disentangle the drivers of biomass patterns (Cai et al., 2019). The explanatory variables in the RDA models were selected using a forward selection method that avoids overestimation of Type I error and explained variance (Blanchet et al., 2008). Two stopping rules, including p > .05, and the adjusted R 2 value of a reduced model exceeding that of the global model, were used to conduct forward selection. Thereafter, we used variation partitioning to uncover the pure and shared effects of three explanatory variable groups (TD, FD, and environmental variables) on biomass (Peres-Neto et al., 2006). RDA analyses were implemented using the "ordiR2step" function in the R package (Oksanen et al., 2020). Variation partitioning analyses were conducted using the function "varpart" in the R package (Oksanen et al., 2020).

| Fish species composition
In this study, a total of 106 species were collected, of which 27 were cold-water species and 28 were warm-water species (Table S1). The species differed in terms of their seasonal preferences ( Figure 2).
Variations in the relative abundances of pholis fangi, Indian perch (Jaydia lineata), branded goby, spiny red gurnard (Chelidonichthys spinosus), and fat greenling contributed the most to the differences in fish communities between spring and autumn (Table S3).
Furthermore, variations in the relative abundances of pholis fangi, pinkgray goby, and fat greenling contributed the most to the differences in fish communities between spring and winter (Table   S4). Variations in the relative abundances of pholis fangi, branded goby, verticalstriped cardinalfish, spiny red gurnard, and pinkgray goby contributed the most to the differences in fish communities between summer and autumn (Table S5). Variations in the relative abundances of pholis fangi, pinkgray goby, branded goby, and striped seasnail (Liparis liparis) contributed the most to the differences in fish communities between summer and winter (Table   S6). Finally, variations in the relative abundances of pinkgray goby, Indian perch, branded goby, spiny red gurnard, and pholis fangi contributed the most to the differences in fish communities between autumn and winter (Table S7).

| Seasonal variations in biodiversity and biomass
One-way ANOVA revealed that the total fish biomass, Simpson's evenness, and all three FD indices (FRic, FEev, and FDiv) Table 2). Biomass values at each station were highest in summer, followed by spring, and lowest in autumn. FRic was highest in winter, when species richness was the lowest. FEev was highest in winter and lowest in autumn, whereas Simpson's evenness was highest in autumn and lowest in spring. FDiv was higher in spring and winter than in summer and autumn, whereas the variance was higher in summer and spring than in autumn and winter. Overall, the patterns of change in the TD and FD indices were inconsistent among the seasons.

| BEF relationships
There were significant relationships between two TD indices (species richness and variance) and biomass across all seasons ( Figure 3).

| Effects of environmental variables on BEF relationships
The GAM significantly fitted curves of both biomass and biodiversity indices, in relation to environmental variables, in at least one season ( Figure 5). Specifically, biomass was found to decline significantly F I G U R E 2 Changes in species composition by weight and number based on the bottom-trawl survey data obtained during spring, summer, autumn, and winter with Chl-a in autumn and winter. SBT was found to significantly impact species richness and FRic in winter; FDiv in autumn; the Simpson's evenness and variance index in spring, summer, and autumn; and FEev in all seasons. Chl-a had a significant impact on species richness and FDiv in spring, autumn, and winter. It also affected Simpson's evenness in all seasons, the variance index in spring and winter, FEev in spring and summer, and FRic and FDiv in autumn and winter.
Variation partitioning analysis showed that 65%, 85%, 64%, and 69% of the variations in biomass were explained by TD indices, FD indices, and environmental variables in spring, summer, autumn, and winter, respectively ( Figure 6). In all three types of scenarios, the observed variation was best explained by TD, with 49% in spring, 78% in summer, 50% in autumn, and 54% in winter, respectively. This was followed by environmental variables in spring (31%), summer (29%), and winter (52%), and by FD in autumn (41%). The least variation was explained by FD in spring (20%), summer (20%), and winter (41%), while none was explained by environmental variables in autumn. Of the unique effects, TD was clearly the most important one in spring (23%), summer (50%), and autumn (24%), while environmental variables uniquely explained the most variation in winter (13%). Other unique fractions only explained a small amount of the total variation, less than 10%, with the exception of FD (15%). Clearly, 7%, 13%, 26%, and 34% variations were explained by all facets of the scenarios in spring, summer, autumn, and winter, respectively. Moreover, there was a larger amount of shared effects among the TD metrics and environmental variables than that among the FD metrics and environmental variables in spring, summer, and winter ( Figure 6).

| Seasonality of species composition
Overall, the species composition was dynamic because of the nonstatic abiotic environment, which is typical of marine ecosystems.
Non-static environments affect the life cycle of marine organisms, as evidenced by obvious seasonal patterns in terms of growth, reproduction, and abundance (Maynou & Cartes, 2000).  in the survey in that season is low and vice versa. For instance, inshore hagfish (Eptatretus burgeri) move to deeper waters (50 m) from shallow areas (6-10 m) between mid-July and mid-October (Ichikawa et al., 2000). Moreover, the head croaker (Collichthys lucidus) often migrates to resource-rich coastal waters for breeding during spring and winter (Huang et al., 2010). Similar to what was seen in those reports, we only collected the head croaker during spring, autumn, and winter, as these fish likely migrated out of the survey area during summer when the species richness was the lowest.
In addition, seasonal species composition was reflected in variations in fishing effort, which was intensive from September to April of the following year because fishing activity was forbidden from May to August (fishing moratorium). The selective removal of large predatory fish from the ecosystem by fishing may lead to declines in predation mortality of smaller organisms, such as small pelagic fish, through species interactions. In this way, the species composition becomes altered with community miniaturization trends, lower trophic levels, and earlier maturation (Engelhard et al., 2015;Heino & Godø, 2002).

| Are changes in TD and FD correlated?
We noted that the seasonal patterns of FD were uncoupled from

| BEF relationships in a shelf sea ecosystem
Four of the six biodiversity indices were significantly related to the biomass in at least one season. These results support the consideration of biodiversity loss in the EBM. Species richness was positively related to total biomass in all seasons, which is consistent with the reports of Gamfeldt et al. (2014) and Yasuhara et al. (2016).
This trend remained consistent across seasons with no significant relationship between Simpson's evenness and biomass, indicating that high biomass levels in the study area were mainly driven by "complementarity" mechanisms. In keeping with our hypothesis, the positive FRic-biomass relationship remained stable across seasons, corroborating the findings reported by Yasuhara et al. (2016).
The relationship between species evenness and biomass is complex because species evenness is sensitive to species interactions, such as competition and predation (Maureaud et al., 2019).
In European seas, which represent an overexploited ecosystem, in winter fish assemblages. In contrast, other reports have suggested a positive relationship between species evenness and biomass (Hodapp et al., 2015). Contrary to their findings, we noted a nonsignificant relationship between Simpson's evenness and biomass in our study area in all seasons. Two key mechanisms drove the high biomass by evenness. At stations with lower evenness and higher biomass, a few dominant species can outcompete many other species to exploit the available resources more efficiently, as proposed by Hodapp et al. (2015). At these stations, we found that some small fish, such as pholis fangi, Indian perch, and pinkgray goby, were more abundant (i.e., low evenness), and served as prey for large predatory fish, such as yellow goosefish and striped seasnails, which ultimately contributes to high biomass. However, at other stations with higher evenness, greater species richness coupled with greater species evenness within a community may contribute to the higher biomass, as observed by Hillebrand et al. (2008).
Although BEF relationships have previously been explored in marine ecosystems (Leduc et al., 2013;Maureaud et al., 2019;Strong et al., 2016), this study offers novel insights into the seasonality of these relationships in demersal fish communities in an overexploited shelf sea ecosystem. Although only one of the six biodiversity indices (i.e., FEev) exhibited seasonal differences in the patterns of its correlations with biomass, our results underscore the need to conduct seasonal scientific surveys. Meanwhile, SBT and Chl-a may contribute to the seasonal differences in BEF relationships because environmental factors are known to shape BEF relationships (Duncan et al., 2015). However, species composition clearly differed across seasons, resulting in seasonal changes in biodiversity indices, further complicating the explanation of biodiversity relative to biomass dynamics. Fluctuations in biodiversity (species richness, community composition, and trait dominance within species) can markedly alter BEF relationships (Brose & Hillebrand, 2016). Fishing is an important contributor to the decline in the biomass of large commercial fish because of the increased risk of depleting spawning stocks and recruitment failure (Edgar et al., 2018). Longterm overfishing throughout the study region has reduced the spawning stocks of some fish, such as largehead hairtail and small yellow croaker (FAO, 2016). Possible changes in previous demersal fish communities due to overexploitation that occurred before our analysis cannot be overruled (Dulvy et al., 2008). Accordingly, the seasonal BEF relationships observed in this study may be the result of past anthropogenic disturbances. Thus, whilst the findings of this study appear to be extensively applicable, they are limited to demersal fish communities in overexploited shelf sea ecosystems.
However, at least in the area examined, we advocate considering the seasonal factors affecting BEF relationships in conservation policymaking and ecosystem management.

| Relative importance of biodiversity and the environment to biomass
According to previous studies, environmental variables may play a more important role than biodiversity indices in explaining the changes in biomass (Yasuhara et al., 2016). However, some studies have found that biodiversity is a major driver of biomass (Balvanera et al., 2014;Cardinale et al., 2012). In this study, both biodiversity and environmental variables were found to significantly affect biomass, which is consistent with the reports of Teichert et al. (2018).
Furthermore, we noted that the relative importance of biodiver-  (Table S8).
This resulted in zero variation being explained by the environment.
Moreover, although SBT was retained in the RDA model in winter with a high percentage of explained variation (47.1%), and Chl-a was F I G U R E 6 Results of variation partitioning, illustrating the contributions of different factors to biomass in spring, summer, autumn, and winter. All fractions are based on adjusted R 2 values shown as percentages of the total variation. The components higher than 5% are in bold font. Negative values are shown as "×." Total adjusted R 2 and residuals (Res) are shown in the lower left and right corners, respectively. Env, environmental variables; FD, functional diversity; TD, taxonomic diversity retained in the RDA model in spring and summer (Table S8), we detected no significance between their relationship with biomass in the respective seasons ( Figure 5). Therefore, more attention should be paid to the interactions among the different environmental variables. In addition, the interactions among biodiversity indices should not be neglected.
Overall, our findings offer valuable insights into realizing sustainable fishery development through the integration of biodiversity loss, environmental changes, and fishing efforts into predictive models to estimate changes in ecosystem functioning. For instance, combining FD with a trait-based model across different spatiotemporal scales may reveal future trends in ecosystem functioning (Holzwarth et al., 2015). Therefore, further studies should combine habitat models to explore BEF relationships and predict future changes in ecosystem functioning in diverse habitats. In addition, information on the BEF pattern and the role of biodiversity in ecosystem management benefits biodiversity conservation. Notably, some studies have demonstrated that EBM can be used to design more effective operational management measures and policies to protect biodiversity . Overall, more effort should be made in research to better understand the BEF relationship to support effective decision-making across spatiotemporal scales and further promote the transition from EBM as an academic concept to practical applications. This may then fortify efforts to better protect biodiversity and, simultaneously, achieve sustainable ecosystem development.

| CON CLUS IONS
This study revealed significant BEF relationships in an overexploited marine ecosystem. Notably, FEev exhibited different patterns of BEF relationships across seasons. Both environmental variables and biodiversity indices affected biomass, and the responses of biomass varied among the seasons. These results indicate that seasonality is an important factor in the response of ecosystem functioning to biodiversity loss and environmental change. Our findings highlight the need to consider seasonal factors and to apply BEF relationships to ecosystem management. Therefore, we recommend the conservation and assessment of ecosystem functioning based on seasonal biodiversity and environmental variables. In the future, we advocate further research on the BEF relationship in other ecosystems and under different contexts (e.g., spatiotemporal scales and trophic levels).

ACK N OWLED G EM ENTS
This work was partially supported by the National Natural Science Foundation of China (NSFC) (Grant No. 41930534 and 41861134037).

CO N FLI C T O F I NTE R E S T
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

PEER R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13521.

DATA AVA I L A B I L I T Y S TAT E M E N T
The