Responses of Functional Traits of Macrobenthic Communities to Human Activities in Daya Bay (A Subtropical Semi-Enclosed Bay), China

The biological trait analysis (BTA) is regarded as a promising approach to unravel how ecosystem functions respond to human-induced disturbances. This study considered the four sampling locations associated with different human activities in Daya Bay, that is, the domestic and industrial sewage discharge area (SED), mariculture area (MRC), nuclear power plants thermal discharge area (NTD), and an area with relatively low human disturbance as a reference (REF). Thirty modalities of nine traits were selected in BTA. Our results showed a clear shift in the functional structure of macrobenthic communities between the sampling locations, except for the case between NTD and REF. The trait composition in the communities did not highlight any seasonal patterns. Bioturbation, longevity, tolerance, body size, feeding habit, and environmental position were the key traits to characterize the functional structure of macrobenthic communities and demonstrated predictable responses along the environmental gradients. Water depth, DO, Chl-a, NH4 +, and petroleum contaminants in sediments were the main variables influencing the trait composition. In addition, the taxonomic index (H′) and functional diversity index (Rao’s Q) showed clear differences among the sampling locations. Although there were no significant differences between NTD and REF in terms of the trait composition and functional diversity, a potential function loss in NTD still can be detected through the integrated analysis with taxonomic diversity. We suggest that the traits (except for fragility, larval development, and living habits) selected and the diversity indices (H′ and Rao’s Q) could serve as promising indicators of ecological conditions in Daya Bay.


INTRODUCTION
Macrobenthic fauna are key components in marine ecosystems and drive important processes such as sediment reworking, bio-irrigation, nutrient cycling, and organic matter decomposition (Widdicombe et al., 2004;Lindqvist 2014;Bonaglia et al., 2017). As their limited migration ability, relatively long life span, and sensitivity to environmental changes, benthic macrofauna are usually regarded as effective indicators for ecosystem health assessment (Snelgrove 1998;Patricio et al., 2009;Borja et al., 2010). Traditionally, biodiversity and how it is related with environmental change and ecosystem health have often been measured and quantified through taxonomy-based approaches, which include taxonomic diversity indices and multivariate methods, based on species richness, abundance, and biomass (Vandewalle et al., 2010;Strong et al., 2015). But these approaches do not take ecosystem functions into account and ignore that different species may perform similar functional roles in communities and ecosystems (Frid et al., 2000;Dolédec and Bonada 2013). Furthermore, ecosystem processes and services are primarily determined by the functional identity, rather than by taxonomic identity (Bellwood et al., 2012;Naeem et al., 2012). Since the beginning of this century, the trait-based approach, characterizing organisms by a suite of functional identities (functional traits), has received an increasing recognition in the marine community ecology (Bremner et al., 2003;Gagic et al., 2015;Barton et al., 2016;Beauchard et al., 2017;Degen et al., 2018).
A representative trait-based approach is biological trait analysis (BTA) which uses a set of biological traits of organisms to depict the variation patterns of the functional structure of communities along environmental gradients (Bremner et al., 2006). Generally, the functional structure of a community consists of two components, which can be quantified by different indices (Dias et al., 2013). The first component, the "community-weighted mean trait value" (CWM), can be calculated for each trait as the average trait value in a community weighted by the relative abundance/biomass of the species carrying each value (Garnier et al., 2004). This index is usually used to identify the dominant traits in a community and can be particularly useful to understand the response of the biotic communities to the changes in disturbances and environmental gradients due to the environmental selection for certain traits (Vandewalle et al., 2010;Ricotta and Moretti 2011). The other component, the functional diversity (FD), describes the extent of trait differences among coexisting species within a community (Petchey and Gaston 2006). A general assumption is that the stressed communities would lead to a decrease in FD , thereby decreasing the diversity in resource use strategies (Strong et al., 2015). FD can be further decomposed into three complementary subcomponents, that is, functional richness, functional evenness, and functional divergence (Mason et al., 2005). An increasing number of indices have been proposed or are being developed for the quantification of one or more of these subcomponents, for example, FRic, FEve, FDiv, FDis, FR, and Rao's Q. FRic (functional richness index), FEve (functional evenness index), and FDiv (functional divergence index), respectively, measure the three subcomponents of functional diversity (Mason et al., 2005;Villeger et al., 2008). Functional dispersion index (FDis) measures the mean distance in the multidimensional trait space of individual species to the centroid of all species (Laliberte and Legendre 2010). FR (functional redundancy index) corresponds to the average number of species sharing similar sets of traits (Mouillot et al., 2014). Rao's Q (Rao's quadratic entropy) incorporates the pairwise functional differences of species weight by their relative abundance (Schmera et al., 2017). Among these indices, CWM and Rao's Q (Rao's quadratic entropy) are most frequently used in previous studies (Ricotta and Moretti, 2011).
Coastal areas host more than half of the planet's human population, depending on the ecosystems and the services they provide (Barbier et al., 2008;Halpern et al., 2008). Most coastal ecosystems have been somewhat altered by human activities at all spatial scales (Lotze et al., 2006;Halpern et al., 2008;Halpern et al., 2015;Cloern et al., 2016). Benthic communities are extremely susceptible to human disturbances, which include domestic and industrial sewage and extensive mariculture, resulting in eutrophication, oxygen deficiency, and contaminant enrichment in sediments (Villnas et al., 2011;Franzo et al., 2016). This situation could become worse, given the weak hydrodynamic force in enclosed and semi-enclosed bays; therefore, the concentration of pollutants in sediments often reaches up to the orders of magnitude higher than in the overlying water (Trannum et al., 2004;Cardellicchio et al., 2007). Recently, more and more nuclear power plants (NPPs) have been or are being constructed for the surge in energy demand in China (Ye et al., 2017). Thermal discharge and the use of radionuclides pose substantial risks to aquatic organisms living nearby after the operation of NPP (Shiah et al., 2006;Teixeira et al., 2009;Thinova and Trojek 2009). BTA has been applied in the marine macrobenthic community to assess human impacts such as bottom trawling (Tillin et al., 2006;Bolam et al., 2014), pollution (Oug et al., 2012;Gusmao et al., 2016;van der Linden et al., 2016;Nasi et al., 2018;Egres et al., 2019;Nunes de Souza et al., 2021), dredging, and dredged material disposal (Cooper et al., 2008;Bolam et al., 2016). Most of the research studies focus on Europe, South America, and Arctic region; the application of BTA on the macrobenthic community is relatively scarce in Asia. Above all, previous studies only consider the impact of singular human activities, while the comparisons on the effects between different human activities are rarely involved.
Daya Bay is embedded in the key economic development zone in south China. It is surrounded by petrochemical industries, mariculture, NPP, and urban development (Wu et al., 2009;Liu et al., 2018), which makes it an ideal paradigm for studying the response of macrobenthic communities to multiple human activities. Numerous previous studies demonstrated that human activities had severely changed the biotic and abiotic environment in Daya Bay (Gao and Chen 2008;Wu et al., 2009;Wang et al., 2011;Ma et al., 2014;Arbi et al., 2017). The mariculture area in Daya Bay has expanded from 440 ha in 1988 to 14,000 ha in 2005, a nearly 600-fold growth over the past decades (Wu et al., 2009). The annual average Chl-a concentration increased by 1.9 mg/m 3 , and harmful algal blooms increased significantly (Song et al., 2009). The cooling sea water from NPP has been discharged at a rate of 95 m 3 /s and with a temperature of 65°C since 2002, and sea water +1°C area can cover approximately 56 km 2 (Yu et al., 2010). The thermal discharges decreased the abundance and diversity of phytoplankton and zooplankton and consequently declined the fishery production (Hao et al., 2016). With the rapid development of the petrochemical industries in Daya Bay, around 1,150 m 3 /h sewage, containing COD, petroleum hydrocarbons, heavy metals, sulfide, ammonia, and so on, is directly discharged into the bay (Xu et al., 2014).
In this study, BTA was used to illustrate the functional responses of macrobenthic communities to human activities. We targeted three major human activities (domestic and industrial sewage discharge, mariculture, and NPP thermal discharge) in Daya Bay and aimed to 1) evaluate potential differences in the functional traits of macrobenthic communities influenced by different types of human activities; 2) identify the most influential environmental variables impacting the functional trait patterns; and 3) discern the most responsive traits to environment stress.

Study Area
Daya Bay is a semi-enclosed bay located in the northern edge of the South China Sea, with a subtropical marine climate. It covers an area of 600 km 2 and has an average water depth of 10 m. The minimal sea surface temperature occurs in winter (average 17.3°C) and the maximum in summer (average 29.3°C). It is mainly influenced by the East Asian Monsoon system, whereby southwest winds prevail in summer and northeast winds in winter. The average annual rainfall is 1984 mm, 80% of which happens in spring and summer (wet season) and 20% in autumn and winter (dry season) (Wang et al., 2008). Tidal current is the dominant current component in Daya Bay which is of the irregular semi-diurnal type. Despite the low velocity (2-6 cm/ s), there also exist two types of seasonal circulation: cyclonic circulation in summer and anticyclonic circulation in winter (Wu et al., 2007). There are mariculture areas (shellfish farming that has lasted over the past decades) in the northeast coast, petrochemical industries (operated since 2001) in the northwest coast, and two nuclear power plants (operated since 1994 and 2002, respectively) in the western coast. Aotou, located in the northwest coast, is the most densely populated town in Daya Bay, with more than 60,000 residents. Over the past decades, the ecological environment of Daya Bay has been severely influenced by the aforementioned types of human activities (Wang et al., 2008;Chen et al., 2010;Qu et al., 2018;Ke et al., 2019).

Field Sampling and Sample Processing
Three stations were set in the four sampling locations associated with different human activities, respectively, resulting in a total of 12 stations (Figure 1). These locations are domestic and industrial sewage discharge area (SED), mariculture area (MRC), sea water +1°C area induced by NPP thermal discharge (NTD), and an area in the mouth of the bay with relatively low human disturbance as a reference (REF). Sampling surveys were conducted seasonally (November in 2017, January, April, and July in 2018, representing autumn, winter, spring, and summer, respectively). At each station, four replicate sediment samples were collected using a 0.05-m 2 Van Veen grab and washed through a 0.5-mm mesh sieve, and the residues were transferred to the sample containers with 5% formalin buffer in situ for faunal analysis. In the laboratory, animals were identified to the lowest possible taxon and enumerated under a dissecting microscope. Taxonomic data from the four replicates were pooled to obtain faunal composition at each sampling station.
The biological and physicochemical samplings were carried out simultaneously. A set of human activity-related environmental variables were selected. One extra sediment sample from the surface up to 10 cm depth was obtained at each station for the analyses of grain size, total organic carbon (TOC), petroleum contaminants (PCs), and heavy metals including As and Hg, which were mainly influenced by the anthropogenic input (Zhao et al., 2016). A total of 1,000 ml water sample was collected from the bottom layer per station for measurements of phosphate (PO 4 3-), ammonia (NH 4 + ), and chlorophyll a (Chl-a). The water samples were stored in PE bottles with special fixatives. All samples were preserved under 4°C prior to analyses. Water depth (Dep) was measured by a handheld depth sounder (Hondex PS-7). Temperature (T), salinity (Sal), and dissolved oxygen (DO) were measured in situ using a handheld multiparameter meter (YSI Pro Plus). In the laboratory, the sediment grain size was determined using a granulometer (Mastersizer 2000) capable of measuring particle Frontiers in Environmental Science | www.frontiersin.org November 2021 | Volume 9 | Article 766580 sizes from 0.02 to 2,000 μm. Shepard's sediment classification method was used to describe the types of sediment (Shepard, 1954). The analytical methods of remaining variables refer to the Marine Monitoring Specification (GB 17378.4-2007;GB 17378.5-2007).

Biological Trait Analysis
Nine biological traits subdivided into 30 modalities were used to describe the functional properties of macrobenthic communities in Daya Bay ( Table 1). The traits included morphological (body size and fragility), life history (larval development and longevity), and behavioral (living habit, feeding habit, bioturbation, tolerance, and environmental position) characteristics. The selection of these traits and their implications of ecosystem functions were based on the recent reviews (Beauchard et al., 2017;Degen et al., 2018;Degen and Faulwetter 2019;Lam-Gordillo et al., 2020). Information regarding the biological traits were gathered from literature (Oug et al., 2012;Queirós et al., 2013;Jumars et al., 2015), online databases (http://www.marlin.ac.uk/biotic/; http://www. marinespecies.org/; https://www.univie.ac.at/arctictraits/; http:// www.marinespecies.org/traits/index.php; http://polytraits. lifewatchgreece.eu/polytraits), and local specialists. When little or no information is available for a taxon, we referred to the information from the same genus or family. A "fuzzy coding" procedure was used to code for each taxon according to their affinity for the modalities of each trait (Chevenet et al., 1994;Bremner et al., 2003). Fuzzy-coding procedure was used to account for multiple trait modalities of an individual taxon, that is, a scoring range of 0-3 was adopted, where 0 expresses no affinity for the given trait category, 1 or 2 express partial affinity, and 3 expresses exclusive affinity. In present study, to avoid bias among different traits, the affinity scores for each trait were standardized so that the sum was equal to 1 (Darr et al., 2014).

Data Analysis
The principal component analysis (PCA) was used to investigate the variation in environmental variables (Z-score normalized) among the sampling locations. The quantification of the functional structure of communities was proceeded after obtaining the taxa matrix (table L, taxa by stations) and the fuzzy coded trait matrix (table Q, traits by taxa). CWM was calculated for each of the 30 trait modalities. The trait modality  (Bremner et al., 2006). In FCA, the correlation ratio (cr) indicates the contribution of each trait to the overall variance, whereas traits with a CR value greater than 0.1 can be considered as the traits with the greatest influence on variance between stations . Then the CWM table was calculated by Bray-Curtis measure to construct the resemblance matrix. Based on this resemblance matrix, the permutational multivariate analysis of variance (PERMANOVA main test, with 9,999 permutations) was performed to test whether the trait composition of macrobenthic communities showed significant differences among the sampling locations/seasons. The similarity percentage (SIMPER) analysis based on the CWM table and taxa matrix with Bray-Curtis measure was performed to determine the dominant traits and species/taxonomic unit characterizing the dissimilarity among the communities.
FD was estimated by Rao's quadratic entropy (Rao's Q) because it incorporates functional richness, functional evenness, and functional divergence . To compare the relationship between functional diversity and taxonomic diversity, the classical taxonomic diversity index, that is, the Shannon-Wiener diversity index (H′, based on log2), was also calculated.
In order to unravel the trait-environment relationships, the BIOENV analysis was used to assess which environmental variables correlate best with the CWM values. The RLQ and fourth-corner methods were further applied to assess which traits response to environmental gradients, based on the following three tables, that is, R table (environmental variables by stations), L table, and Q table (Dray et al., 2014).
Non-parametric Kruskal-Wallis tests and Wilcoxon tests were performed to test whether the median values of environmental variables and indices (functional and taxonomic) showed significant differences between the sampling locations and between seasons. Kruskal-Wallis tests, Wilcoxon tests, PCA, FCA, and the calculation of functional and taxonomic diversity indices were run in R-3.5.3 (R Core Team, 2019) using the packages "stats," "factoextra," "vegan" (Oksanen et al., 2019), "ade4" (Dray and Dufour, 2007), and "FD" (Laliberté et al., 2014). PERMANOVA main test, SIMPER, and BIOENV were performed in PRIMER 7 (with PERMANOVA+).

Environmental Variables
For environmental variables, notable differences among the sampling locations and among seasons were observed (Supplementary Table S1). Compared with the stations in SED and MRC, the stations in NTD and REF had higher values of water depth and salinity, whereas they showed a lower content of PCs and Hg in sediments (p < 0.05). SED showed a higher concentration of Chl-a than the other sampling locations (p < 0.05) and a higher content of TOC in sediments than MRC. Salinity, temperature, DO, Chl-a, PO 4 3− , and NH 4 + differed significantly between seasons (p < 0.05). The highest value of salinity occurred in autumn. The highest content of PO 4 3− and NH 4 + occurred in spring, while the highest concentration of Chl-a occurred in summer. The sediment types in the four sampling locations were all silty clay.
The first two principal components of the PCA accounted for 51.0-56.0% of the total variation in environmental conditions over the studied stations ( Figure 2). The seasonal variation of the environmental conditions of the sampling locations was clear. But, in general, SEDs were associated with higher concentration of Chl-a and higher contents of petroleum PCs and Hg in sediments, NTD had greater percentages of silt and clay, REF was characterized by greater water depth. MRC had relatively large intra-group variations and did not show concordant characteristics among seasons.

Functional Structure of Macrobenthic Communities
A total of 163 taxa, including 80 polychaetes, 50 molluscs, 25 crustaceans, 3 echinoderms, 2 sipunculid, 1 platyhelminth, 1 nemertean and 1 chordate, were identified during the four sampling surveys in Daya Bay. Taxa matrix and fuzzy coded trait matrix are shown in Supplementary Tables S2, S3. FCA depicted the patterns in the trait composition for the stations grouped by the sampling locations and seasons ( Figures 3A,B). The first two FCA axes accounted for 56.4% of the variance in the trait composition between the stations, 38.1% along axis 1, and 18.3% along axis 2. The traits, that is, bioturbation, longevity, tolerance, body size, and feeding habit were significantly related to axis 1, and environmental position was strongly correlated with both axes ( Table 2).
For the trait composition of macrobenthic communities, a significant difference among the sampling locations was confirmed by the PERMANOVA main test (Pseudo-F 9.587, p < 0.001). However, no significant difference among seasons was detected (Pseudo-F 1.329, p 0.215). Therefore, the comparisons of the trait composition between seasons will not be carried out in the following analyses. Pair-wise tests showed that there were significant differences between the sampling locations (t 2.031-4.455, p < 0.01), excluding the case between NTD and REF (t 1.268, p 0.152).
The average values of CWM for each sampling period were plotted to visualize the functional composition of macrobenthic communities in the sampling locations (Figure 4). The SIMPER test performed on CWM values showed that the macrobenthic communities of the four sampling locations were all characterized by the following trait modalities (Table 3), that is, fragile (F1), pelagic larval development (LD1), short longevity (A1), borrower (LH3), deposit feeder (FH1), and high tolerance (T3). Nevertheless, SED and MRC exhibited a higher percentage of Frontiers in Environmental Science | www.frontiersin.org November 2021 | Volume 9 | Article 766580 5 small body size (S1) than NTD and REF, and MRC, NTD, and REF had higher percentages of small-medium body size (S2) than SED; NTD, and REF had higher percentages of medium longevity (A2) and biodiffusor (B2) than SED and MRC, while SED was dominated by upward conveyor (B3) and MRC was dominated by surficial modifier (B1); REF had a higher percentage of intermediate tolerance (T2) than the other sampling locations; SED had a higher percentage of epifauna (EP2) than the other sampling locations. The SIMPER test performed on taxa matrix showed that ( Table 4) the macrobenthic communities of SED were represented by two Polychaeta species (Polydora sp. and Prionospio polybranchiata); MRC were represented by one Nemertea species and six Polychaeta species (Cirratulidae, Goniada maculata, Mediomastus chinensis, Notomastus latericeus, Paraprionospio cristata, and Sigambra hanaokai); NTDs were represented by one Nemertea species and five Polychaeta species (Aglaophamus lyrochaeta, Cirratulidae, Cossura dimorpha, M. chinensis, and S. hanaokai); REFs were represented by one Nemertea species, nine Polychaeta species (A. dibranchis, A. lyrochaeta, Cirratulidae, Glycera sp., Laonice cirrata, Listriolobus brevirostris, Lumbrineris sp., M. chinensis, and S. hanaokai), one Crustacea species (Photis hawaiensis), and two Echinodermata species (Amphiodia (Amphispina) obtecta and Protankyra bidentata).
There were significant differences in Shannon-Wiener diversity (H′) and functional diversity (Rao's Q) among the sampling locations (p < 0.01). The highest H′ and Rao's Q were both observed in REF, and the lowest values both occurred in SED ( Figures 5A,B). The H′ was significantly higher in REF than that in NTD (p 0.013), but there was no significant difference in Rao's Q between them ( Figures 5A,B). The H′ was slightly higher in MRC than that in NTD, and Rao's Q was slightly higher in NTD than that in MRC ( Figures 5A,B). However, there were no significant differences between MRC and NTD for both indices. Additionally, Rao's Q responded rapidly as H′ increased in SED, MRC, and NTD, but for REF with relatively Frontiers in Environmental Science | www.frontiersin.org November 2021 | Volume 9 | Article 766580 6 higher species diversity, the variation of Rao's Q was small ( Figure 5C).

Relationships Between Traits and Environmental Variables
According to the BIOENV analysis, functional trait pattern was best correlated with the combination of Dep, DO, Chl-a, NH 4 + , and PCs (ρ 0.470). The combination of RLQ and fourth-corner analysis revealed 14 out of the 30 trait modalities were significantly correlated (p < 0.05) with the environmental variables ( Figure 6). Specifically, the frequency of small body size (S1) was positively associated with the contents of NH 4 + , As, Hg, and PCs, and negatively associated with water depth, salinity, and the percentage of clay. Small-medium (S2) was negatively associated with the contents of NH 4 + and Hg. Pelagic larval development (LD1) was more abundant in higher percentage of clay sediments and benthic larval development (LD2) had the opposite trend. Animals with short longevity (A1) were more abundant in PCs/Hg enriched sediments but less abundant in greater water depth, salinity, and clay percentage environments. The frequency of medium longevity (A2) was positively influenced by water salinity, and negatively influenced by PCs, silt, Hg, and NH 4 + . Filter/suspension feeder (FH2) was positively correlated with NH 4 + , and predator (FH4) was positively correlated with water salinity. Upward conveyor (B3) and downward conveyor (B4) were more abundant in DO, Chl-a, and PCs enriched environments, while the surficial modifier (B1) was less abundant in DO enriched environments, and downward conveyor (B4) was less abundant in PO 4 3--enriched environments. Animals with high tolerance (T3) were scarce in deeper waters. Epifauna (EP2) was more abundant in NH 4 + -, As-, and Hg-enriched environments and in fauna (EP1) had the opposite trend.

Effect of Anthropogenic and Natural Drivers on Environmental Conditions in Daya Bay
Previous studies have confirmed the significant impacts of human activities on Daya Bay environments, for example, petrochemical sewage discharge (Xu et al., 2014;Liu et al., 2018), mariculture (Qi et al., 2019), and thermal discharge (Yu et al., 2010;Hao et al., 2016). In the present study, water depth was associated with environment stress. When water depth increased, the content of nutrients (PO 4 3-, NH 4 + ) in water columns and contaminants (PCs, Hg) in sediments decreased, and vice versa. It can be partly attributed to the human activities such as sewage discharge and mariculture that increase the chemical pollution load located in shallow waters. However, the relatively low Chl-a and TOC contents in MRC may be due to its higher sea water renewal and the dilution and dispersion rates of contaminants. MRC is located in the eastern part of Daya Bay, which is the main channel for water exchange with the outer sea (Wu et al., 2007). Due to the  diurnal fluctuations in water temperature, the heating effect of thermal discharge was not reflected from our in situ surveys. Furthermore, the precipitation and the East Asian Monsoon also have an important influence on water quality in Daya Bay. We observed that the higher content of PO 4 3− , NH 4 + , and Chl-a occurred in wet season (spring and summer). When the East Asian Monsoon prevails, the nutrients are washed out from the land into the bay, causing differences in the water quality between the wet season and dry season (Wu et al., 2009). Additionally, the eastern Guangdong coastal upwelling promotes primary productivity explaining the high Chl-a concentration in summer (Gu et al., 2012). As a result, the environmental conditions of Daya Bay were determined by both anthropogenic and natural factors.

Functional Response of Macrobenthic Communities to Human Activities
Assessments of benthic community responses to human activities are an important task when considering the significance of the impacts on the sustainability of marine ecosystem functions. Trait distribution emerges through evolution and natural selection and is mediated by the environment, biological interactions, and anthropogenic drivers (Barton et al., 2016). Given that the sensitivity of traits to changes in habitat disturbance, the order of extinction of traits matters for assessing the decline in ecosystem functioning (Boersma et al., 2016;Kenny et al., 2018). For instance, organic pollution affects mainly organisms sensitive to oxygen depletion, whereas acidification affects those sensitive calcium depletion (Bolam et al., 2017).
We explored the functional structure of macrobenthic communities in Daya Bay subject to the potential impacts from different human activities, in particular, assessing the potential effects of contaminants and other environmental variables on macrofaunal functional identity and diversity. The FCA carried out on trait composition highlighted the differences in the sampling locations associated with different human activities (see Figure 3A), a result corroborated also by the PERMANOVA main test (Pseudo-F 9.587, p < 0.001). However, no significant differences were observed between NTD and REF (t 1.268, p 0.152). We found that bioturbation, longevity, tolerance, body size, feeding habit, and environmental position, were key traits in characterizing the functional structure of macrobenthic communities (See  Table 1 for the meaning of the trait codes. Frontiers in Environmental Science | www.frontiersin.org November 2021 | Volume 9 | Article 766580 Table 2). Small body size (<10 mm), short longevity (<2 years), and high tolerance animals were more abundant in SED and MRC (higher proportions in SED), which were represented by the species of Spionidae, that is, Polydora sp., P. polybranchiata, and P. cristata. These are typical features of opportunistic species within nutrient-enriched or stressed habitats (Oug et al., 2012;Nasi et al., 2018). Although fragile animals were dominant in Daya Bay, this was more evident in SED and MRC, indicating they were not suitable for animals with shells, for example, mollusks and crustaceans. The dominant bioturbation mode in SED was upward conveyor (B2). Upward conveyors are efficient at moving deep horizon particles to the sediment surface    Table 1 for the meaning of the trait codes.
Frontiers in Environmental Science | www.frontiersin.org November 2021 | Volume 9 | Article 766580 10 (Kristensen et al., 2012). Surficial modifiers (B1) in general with filter/suspension feeding habit, a greater amount of food provided by mussel biodeposition can explain why they were more abundant in MRC (Lacoste et al., 2019). The activities of surficial modifiers are restricted to the surficial layer (0-2 cm) of the sedimentary profile (Queirós et al., 2013); therefore, the limited bioturbation potential could lead to the deterioration of sediment quality in MRC. Less infauna (EP1) found in SED was supposedly related to the anoxic sediments there. Comparatively, small-medium body size (S2), medium longevity (A2), and intermediate tolerant (T2) species were more abundant in NTD and REF (higher proportions in REF), indicating their relatively good environmental status. A twodecade study demonstrated that the body size of the macrobenthic fauna became smaller and biodiversity was increasingly simplified in Daya Bay (Wang et al., 2008). It is true, from our findings, that a small Polychaeta species Sigambra hanaokai was formerly rare in Daya Bay and has recently become ubiquitous.
On the other hand, we found a similar functional structure between seasons (Pseudo-F 1.329, p 0.215), indicating low temporal variability for the trait composition in macrobenthic communities. While the taxonomic composition was quite different between seasons (Pseudo-F 1.682, p 0.006). Most benthic organisms reproduce and settle at specific times of the year (Reiss and Kroncke, 2005), which makes it difficult to compare works conducted in different seasons. Yet the seasonal stability of functional structure seems to make up for this deficiency.
The taxonomic and functional diversity indices showed clear differences between different human activities ( Figures  5A,B). The relatively lower taxonomic and functional diversity were found in the disturbed areas (i.e., SED, MRC, and NTD). It is consistent with previous studies which noted that the inner bay of the Daya Bay disturbed by mariculture and sewage pollutions had lower taxonomic richness and diversity (Arbi et al., 2017), and the biomass and species richness near nuclear power plants had dropped dramatically relative to the results of pre-operational surveys (Wang et al., 2008). Compared with MRC, NTD had relatively lower taxonomic diversity but higher functional diversity. This result suggested that thermal discharge might cause a decline in taxonomic diversity but not affect the functional diversity, if the reduced species contribute little to function in the community. Thermal discharge have been shown to affect the abundance, richness, and feeding habit of freshwater macroinvertebrate (Worthington et al., 2015;Han et al., 2017). However, most previous studies have demonstrated that the influence of thermal discharge on marine macrobenthic community was insignificant (Lardicci et al., 1999;Kwon et al., 2017;Jung et al., 2018;Lin et al., 2018), possibly benefitting the huge bulk of seawater. Kwon et al. (2017) assumed that if thermal discharge could exert an influence on the macrobenthic community depending on the season, summer may be a bad one. Lin et al. (2018) admitted that the elevated seawater temperature at the bottom might not be high enough to significantly affect the macrobenthos, but the stratification and elevated temperature could affect the phytoplankton at the bottom. In fact, the shift of pelagic food web in thermal discharge areas has been identified (Chen et al., 2015;Hao et al., 2016). The relatively low proportion of deposit feeder (FH1) in NTD might be a response to this shift. Deposit feeders ingest sediments and use organic matter and microbial organisms in sediments as food, such as benthic diatoms. However, the effects of thermal discharge on benthic diatoms have been demonstrated (Ingleton and McMinn, 2012). Despite the higher level of taxonomic diversity, REF showed little increase in functional diversity over NTD. In other words, there were a larger number of functionally similar species in REF. This indicates that this community potentially has a higher functional redundancy than the other communities. Functional redundancy is a policy that insures against changes in the functioning of the ecosystem due to the loss of species (Dolédec and Bonada, 2013). Although NTP and REF are similar in functional structure, a potential decrease in functional diversity in NTP can be expected when environmental degradation continues. Therefore, the combination of taxonomic and functional diversity can be a reliable method for ecological risk assessment.
In the present study, the trait composition was primarily influenced by water depth, DO, Chl-a, NH 4 + , and petroleum contaminants (PCs) in sediments. Generally, water depth and sediment type are considered as the main natural factors that obviously influence the distribution of macrobenthic fauna (Bolam et al., 2017;Kenny et al., 2018;Liu et al., 2019). Although natural factors showed influence on the contaminants, given the consistency of the sediment type (silty clay) and the small water depth range (4-16 m) in Daya Bay, we assumed that the humaninduced environmental changes were the main factors affecting macrobenthic communities. Combined with the correlations between the environmental variables and trait modalities ( Figure 6), we suggested that the traits "body size," "longevity," "feeding habit," "bioturbation," "tolerance," "environmental position," and the diversity indices H′ and Rao's Q could be served as the promising indicators of how macrobenthic fauna react to the stress caused by human activities in Daya Bay.

Ecosystem Functioning Effects
The present study on macrobenthic communities in Daya Bay showed that the BTA method detected specific functional features correlated with human disturbances, for instance, the change from small body size organisms to larger body size organisms along pollution gradients. Of particular interest is the demonstration of the impact of human activities on body size, feeding habit, and bioturbation. All these traits express the aspects of biological activities that are important for sediment reworking, and nutrient uptake and turnover. It is a first step toward understanding how human disturbances can affect ecosystem functioning. However, the accurate assessment of the magnitude of biogeochemical processes loss in response to functional structure shifts still requires further studies (Gusmao et al., 2016).

Cautions and Limitations
The performance of the BTA relies on the selected traits and metrics used during the analytical process (Bremner, 2008). In the present study, a priori selection of traits was based on a trade-off between their efficacy for describing ecological functioning and the effort required to gather information on the selected traits. Some traits, such as reproductive type and mobility, are especially important to analyze recolonization when disturbances are absent. At present, it seems unrealistic to incorporate all possible functional features that could be of importance. Darr et al. (2014) suggested that biomassbased BTA represent an informative tool to describe the functional structure of macrobenthic community. In contrast, Gusmao et al. (2016) considered that abundance is a more suitable metric than biomass. In fact, the discrepancy may depend on the selected traits. Because some traits (e.g., larval development) are linked to abundance and others (e.g., bioturbation) are directly related to biomass. However, how to effectively integrate abundance and biomass into the corresponding traits is still a challenging question.

CONCLUSION
In the present study, BTA was used to discern shifts in the functional structure of macrobenthic communities associated with different human activities. The effects of contaminants and other environmental variables on macrofaunal functional identity and diversity were observed. However, the seasonality of trait expressions in communities was not significant. Bioturbation, longevity, tolerance, body size, feeding habit, and environmental position were the traits most affected by the human disturbances. And water depth, DO, Chl-a, NH 4 + , and petroleum contaminants (PCs) in sediments were the main factors related to the trait composition. Small body size, short longevity, and high tolerance species were more abundant in the contaminated areas (SED and MRC). The taxonomic and functional diversity indices revealed obvious differences among the sampling locations. Although NTD and REF were similar in trait composition, a potential function loss in NTD can be detected through the integrated analysis with taxonomic diversity.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
LC conceived of the study and obtained funding. LC, YR, XC, XZ, SF, and HH were responsible for field and laboratory work. LC conducted the identification of benthic macrofauna and created many of the manuscript elements. YR led the data analysis and writing effort, with significant contributions to these by all authors.