Assessment of the Ecological Status of Rongjiang Estuary (China) Under Human Pressure, Using Biotic Indices Based on Benthic Macroinvertebrates

Rongjiang River, the second largest river system in Guangdong Province, flows through the main urban areas of Jieyang and Shantou cities before reaching the South China Sea. Human activities in the surrounding area pose significant threats to this aquatic ecosystem. The ecological status (ES) of the benthic ecosystem of the Rongjiang River estuary has not yet been conducted using indices based on the macrobenthic fauna, which is important for evaluating environmental health. Here, we used four biotic indices (the AZTI’s Marine Biotic Index (AMBI) and Multivariate AZTI’s Marine Biotic Index (M-AMBI), and taxonomic distinctness indices (average taxonomic distinctness Δ+ and variations in taxonomic distinctness Λ+) to appraise the current ES of benthic communities in Rongjiang estuary. Samples were taken from 11 sampling sites located in six general regions: western aquaculture zones, aquaculture zones, Hanjiang River water channel, Shantou City, Shantou Port, and near the ocean. The benthic ecosystem of this estuary is greatly disturbed: the ES of the aquaculture zones and the sites near the Hanjiang River water channel, Shantou City, and close to the ocean was poorer compared with that of other areas; ES was also poorer in winter than in summer. Generalized linear models revealed that Shannon-Wiener index was negatively correlated with dissolved inorganic nitrogen (p < 0.01), M-AMBI was negatively correlated with temperature and dissolved inorganic nitrogen (p < 0.05 and p < 0.01, respectively), and Λ+ was negatively correlated with pH (p < 0.05). The AMBI, M-AMBI, Δ+, and Λ+ indices were suitable for assessing the ES of the benthic ecosystem in an anthropogenically disturbed estuary.


INTRODUCTION
With their complex hydrodynamics, estuaries have been important areas for human settlement, and they also represent some of the most productive natural ecosystems on Earth (Costanza et al., 1997;Lotze et al., 2006;van der Linden et al., 2012). Their health is key to the development of human society, but anthropogenic and natural stressors have increased estuarine ecosystem deterioration and biodiversity loss (Barbier et al., 2011;Gao et al., 2019;Hillman et al., 2020;Spreitzenbarth and Jeffs, 2021). For example, many anthropogenic pollutants such as dissolved inorganic nitrogen (DIN) can lead to eutrophication, which decrease the survival of sensitive species, reduces species diversity, and simplifies food webs (Alexander et al., 2017;Gao et al., 2019;Ji et al., 2020). In addition, high variability in the natural environment (e.g. hydrodynamics) can also affect the structure of biological communities because only a few highly dominant species are capable of thriving in such environments (Elliott and Quintino 2007;Haase et al., 2012;Van der Linden et al., 2017).
Macrobenthic invertebrates are an important component of aquatic ecosystems, and they have been commonly used to assess the ecological status (ES) of the benthos (Díaz et al., 2004;Pinto et al., 2009). The population dynamics of benthic taxa are sensitive to stresses such as human disturbance because benthic taxa show high site fidelity, have long life spans, and often reside at the sediment-water interface where many pollutants concentrate (Dauer et al., 1993;Dauvin et al., 2007;Dauvin et al., 2012). Selection of the most suitable macrobenthos indices for use in marine ecological assessments remains a major challenge given the high diversity of biotic indices that have been developed. The AZTI's Marine Biotic Index (AMBI) (Borja et al., 2000) and Multivariate AZTI's Marine Biotic Index (M-AMBI) (Muxika et al., 2007) are determined by classifying species into one of five ecological groups (EGs), and these indices were shown to be the most effective for evaluating marine ecological quality among 35 benthic indices examined (Patrícioet al., 2009;Borja et al., 2015). The other two indices of taxonomic distinctness, average taxonomic distinctness Δ + , and variations in taxonomic distinctness Λ + (Clarke and Warwick 1998;1999;Warwick and Clarke 1995) are based on diversity and related attributes and have also been shown to be useful for terrestrial and marine ecological assessments (Hu and Zhang 2016;Arbi et al., 2017). All of the aforementioned indices have been used in the seas around China (Hu and Zhang 2016;Arbi et al., 2017;Qiu et al., 2018;Li et al., 2021). However, these indices might not be effective in some circumstances (Borja et al., 2000;Salas et al., 2006;Muxika et al., 2007;Patrícioet al., 2009;Li et al., 2021). For example, the robustness of AMBI and M-AMBI is reduced if only 1-3 taxa or individuals can be identified in a sample, or more than 20% of taxa in a sample are not assigned to an EG (Borja and Muxika, 2005;Muxika et al., 2007). In addition, too much interference information for the taxonomic diversity (TD) indices-Δ + and Λ + is added by tolerant taxa, such as some tolerant Oligochaeta and Mollusca, because these species contribute the greatest path lengths, which lead to equal or even higher scores in disturbed sites than in undisturbed sites (Hu and Zhang 2016). Given the complex nature of biological systems, robust ES evaluations require the use of several biotic indices (Salas et al., 2006;Wetzel et al., 2012;Martinez-Haro et al., 2015).
Rongjiang River, the second largest river system in Guangdong Province, flows through the urban environment between Jieyang and Shantou cities before entering the South China Sea (Wang et al., 2016). Aquaculture ponds on this river between these two cities cover approximately 21,297 ha (Gu et al., 2018). The aquaculture, industrial pollution, agriculture, domestic sewage, and other human activities surrounding this river and estuary all pose significant threats to the aquatic ecosystem (Gu et al., 2018;JYEPB Jieyang Environmental Protection Bureau 2018;STEPB 2018). Furthermore, the estuary ecosystem of this river is affected by other areas, through the water channel between this river and the Lianjiang River and Hanjiang River. The ES of the benthic ecosystem of the Rongjiang River estuary has not yet been fully evaluated because no ES assessment of the Rongjiang River estuary has been conducted using indices based on the macrobenthic invertebrate fauna, which is important for evaluating environmental health.
The aim of this study was to 1) evaluate the ES of the benthic ecosystem of Rongjiang estuary using the AMBI, M-AMBI, Δ + , and Λ + indices; 2) assess the suitability of the AMBI, M-AMBI, Δ + , and Λ + indices for evaluating the ES of an anthropogenically disturbed estuary system; and 3) explore the environmental factors affecting the ES of the benthic ecosystem.
Frontiers in Environmental Science | www.frontiersin.org November 2021 | Volume 9 | Article 728196 pond aquaculture area is more than 2000 ha. Shantou Port is one of China's key ports in the Belt and Road Initiative, which exchanges cargo with 268 other ports in 57 countries (Li, 2013). Rongjiang estuary has a mild, wet subtropical climate with prevailing northeasterly and southwesterly monsoons from November to May and April to September, respectively. Thus, the area experiences wet and dry (70 and 30% freshwater discharge, respectively) seasons (Chen, 2010). Sampling sites were located within six general regions: western aquaculture zones (S1 and S2), aquaculture zones (S3-S6), Hanjiang River water channel (S7), Shantou City (S8 and S9), Shantou Port (S10), and near the ocean (S11).

Sampling and Measurements
Samples were taken from 11 sampling sites (S1-S11) from the upper to lower reaches of Rongjiang River (Figure 1). The macrobenthos was sampled in June 2017 (wet season) and January 2018 (dry season). Benthic invertebrates were sampled using a grab sampler in an area of 0.1 m 2 . Three replicates were taken at each site to ensure the reliability of the data, and the grab sample was retained if it was more than 75% full. Samples were gently washed and sieved through 0.5 mm mesh; they were then preserved in 4% borax-buffered formalin before identification and counts.
We measured 18 environmental factors associated with various ecological factors. These environmental factors were grouped into two classes−natural variables and anthropogenic pollutants−each with similar attributes. The natural variables were the temperature, pH, salinity, and dissolved oxygen (DO) of the bottom water and sediment grain size. The anthropogenic pollutants were chemical oxygen demand (COD Mn ), suspended matter concentration (SS), chlorophyll a (Chla), dissolved inorganic phosphorus (DIP), dissolved inorganic nitrogen (DIN) of the bottom water, metal concentrations (As, Pb, Zn, Cd, Cr, Cu and Mn) and the soil organic matter (SOM) content of the sediment.
The temperature, pH, salinity, and DO of the bottom water were measured in the field using a multiparameter water quality monitor (YSI 6600, Brannum Lane Yellow Springs, United States). COD Mn was measured by the permanganate method, Briefly, NaOH was added to water samples, followed by KMnO 4 digestion solution. The sample were heated for 10 min in a water bath at 96-98°C. After digestion, the sample was cooled immediately in an icebox to quench the reaction. The remaining permanganate in the sample solution was determined by titration. SS samples were filtered through pre-weighed Whatman GF/F fiber filters (25 mm). Chla was determined by filtering 100-200 ml of seawater onto GF/F fiber filter by a cascading filtering device under low vacuum pressure. After extraction with 90% acetone, Chla was determined by a Turner Design fluorometer (TD Trilogy) (GB/T17378.4-2007; GB/T12763.6-2007) (General Administration of Quality Supervision, 2008a;2008b). Water column nutrient variables were analyzed following the methods of Kirkwood et al. (1996). Measurement of the SOM content, metal concentrations and grain size was based on the methods of Gu et al. (2014Gu et al. ( , 2018.

Requested Potential Ecological Risk Index and Mean Effects Range-Median Quotient Based on the Metal Concentrations in Sediments
Heavy metals always occur in sediments as complex mixtures; therefore, the potential ecological RI and mean ERM quotient (MERMQ) method were used to determine the comprehensive biological effect of combined toxicant groups (Hakanson, 1980;Long, 2006). The equation (Eq. 1) is as follows: where C i f is the contamination factor, C i surface is the measured concentration of the examined component (i) in surface sediment, C i n is the standard preindustrial reference level, RI is the requested potential ecological risk index, Er i is the potential ecological risk factor for the metal i, Tr i is the "toxic-response" factor for the given metal i, C i is the measured concentration of the examined component (i) in sediment, ERM i is the ERM for metal i, and n is the number of metals.

Benthic Indicators and Their Relationships to Environmental Factors
Species richness (S), Shannon-Wiener index (H′), AMBI, and M-AMBI were calculated using AMBI index software (version 5.0, http://ambi.azti.es), and the assignment of the identified species to the five EGs was based on the species list (2017 June) in the AMBI index software (Borja and Muxika, 2005). Reference conditions play a crucial role in calculations of M-AMBI, and these vary naturally with location, water type, and habitat. Following Borja et al. (2008), the highest H' and S values increased by 15% used as the reference conditions in our study. Threshold levels and ES for AMBI and M-AMBI are shown in Table 1 (Borja et al., 2000;Borja and Tunberg, 2011).
Δ + and Λ + indices were calculated by PRIMER 6 (v. 6.1.12) software using six taxonomic levels (species, genus, family, order, class, and phylum). Pollution status can determined by the location of Δ + in a funnel plot. Disturbed sites are usually located far from the 95% confidence interval, whereas undisturbed sites are the 95% confidence interval. In some situations, Δ + results were complemented by those of Λ + reflecting taxonomic unevenness.

Statistical Analysis
A bivariate correlation analysis based on Spearman's correlation coefficient of all environmental variables was conducted to filter the variables; when two variables were highly correlated within the same variable class (r > 0.75, p < 0.05), one was discarded. A heatmap was made to visualize the spatial−temporal distribution patterns of standardized environmental factor data. Generalized linear models (GLMs) were used to select the predictor variables that could significantly explain the ES of the benthic ecosystem. A step-wise multiple regression was then used to identify the best-fit models (variable selection and model evaluation were both based on Akaike information criterion (AIC) values). The bivariate correlation analysis was performed by PRIMER 6 with PERMANOVA + add-on (v. 6.1.12), and the heatmap of environmental factors and GLMs were conducted in R 3.6.2. Sampling maps were plotted using Surfer 14.0 software.

Environmental Conditions
DO was significantly correlated with temperature and salinity (r −0.84, p < 0.01; r 0.80, p < 0.01, respectively), and SOM was significantly correlated with grain size and MERMQ (r 0.79, p < 0.05; r 0.87, p < 0.05, respectively). DO and SOM were thus omitted from subsequent analyses.
The temperature, Chla, RI, and SS in the bottom water were higher in summer than in winter; pH, salinity, DIN, DIP, and MERMQ showed the opposite pattern. COD Mn in the aquaculture zones was lower in summer than in winter; in all other areas, the opposite pattern was observed. Temperature decreased from the upper river reaches to a site near the Hanjiang River water channel and then increased; pH showed a similar pattern to temperature in summer but the opposite pattern in winter. Salinity increased from the upper to the lower river reaches; the site near the Hanjiang River water channel was the only exception to this pattern. Chla and RI were higher in the aquaculture zones and areas near Shantou City. The patterns of DIN and MERMQ varied in different periods. In summer, DIN and MERMQ were lower in the aquaculture zones; in winter, DIN and MERMQ were higher in the aquaculture zones and areas near Shantou City. COD Mn was higher in the aquaculture zones than in the other areas. SS was higher in the western aquaculture zones and lower reaches; DIP was lower in the western aquaculture zones. Surface sediments were generally dominated by silts, and the silt, sand and clay contents being 88.76%-99.13%, 0.84%-14.24% and 0-0.04%, respectively (Supplementary Appendix SA).

Macrobenthic Faunal Composition
A total of 22 species from 6 phyla, 9 classes, 16 orders, 19 families, and 21 genera were identified. In both summer and winter, the greatest number of species detected belonged to the phylum Annelida (8 species), followed by Mollusca and Arthropoda (5 species each). However, Annelida species only comprised 9.4% of all specimens sampled; Mollusca and Arthropoda species comprised 57.0 and 32.3% of all specimens sampled, respectively. We identified 13 species from these phyla in both summer and winter. Most species belonged to Annelida in winter and Arthropoda and Mollusca in summer (Supplementary Appendix SB).
AMBI and M-AMBI One (7.7%) species was not initially assigned to an EG in either summer or winter, and this species remained unassigned in both seasons after assignment (Supplementary Appendix SB).
Four sites (S1, S3, S10 and S11) could not be sampled in summer because of inclement weather. The mean AMBI for summer ranged from 0.75 to 5.99: the western aquaculture zones (S2) were slightly disturbed (good ES). The aquaculture zones (S4-S6) with extensive pond aquaculture and multiple water channels from the Rongjiang and Lianjiang Rivers were slightly to heavily disturbed (poor ES). The site near the Hanjiang River water channel (S7) was slightly disturbed, and areas near Shantou City (S8 and S9) were moderately disturbed and undisturbed, respectively. In winter, the mean AMBI ranged from1.50 to 5.25, and the ES of each region was as follows: the western aquaculture zones (S1 and S2) were moderately disturbed (moderate ES); the aquaculture zones (S3-S6) were slightly disturbed (good ES); S8 (near Shantou City) and S10 (near Shantou Port) were slightly disturbed (good ES); S7 (near the Hanjiang River water channel) was heavily disturbed (poor ES); and both S9 and S11 were moderately disturbed (moderate ES) (Figures 2A,B).
According to the M-AMBI, in summer, the ES of the western aquaculture zones (S2) was high, and that of the aquaculture zones S4, S5, and S6 was good, high, and bad, respectively. The ES of the site near the Hanjiang River water channel (S7) and areas near Shantou City (S8) was poor, and the ES of S9 was good. The ES was poorer in winter than in summer. The ES of S1 and S2 in the western aquaculture zones was good and moderate, respectively; the ES of aquaculture zones (S3, S5, and S6) was poor; the ES of (S4 and S7) was poor; the ES of S10 was high; and the ES of S8, S9, and S11 was poor, moderate, and bad, respectively ( Figures 2C,D).
Taxonomic Distinctness Indices: Δ + and Λ + Confidence plots for summer revealed that most Δ + values were within the 95% confidence interval, and S4 and S8 were close to the upper and lower limit of the 95% confidence intervals, respectively; S7 fell outside the 95% confidence limits. In winter, more sites (S2, S3, S5, S6, S8, and S11) fell outside the lower 95% confidence limit; S9 was close to the lower limit; and S1 and S4 were close to the upper 95% confidence limit. In some cases, the Δ + results were complemented by those of Λ + , reflecting taxonomic unevenness. The Λ + plot differed from the Δ + plot because as S4, S6, S7, and S8 in summer were near the lower 95% confidence limit (Figure 3).  Frontiers in Environmental Science | www.frontiersin.org November 2021 | Volume 9 | Article 728196 5

ES of the Benthic Ecosystem and its Relationship With Environmental Factors
Estuaries are transitional zones between rivers and oceans and are characterized by complex hydrodynamics and chemical conditions (Van der Linden et al., 2017); they are also some of the most anthropogenically disturbed aquatic systems and are sensitive to environmental changes (Wetzel et al., 2012;Gao et al., 2019;Hillman et al., 2020). High variability in this physicochemical has a substantial effect on biological communities, and only a few highly dominant species are capable of thriving in such environments (Elliott and Quintino 2007). Hydrodynamics imposes strong selective forces on the macrobenthos through direct physical stress (Warwick and Uncles 1980); top-down control (Leonard et al., 1998), regulation of larval or juvenile transport (Haase et al., 2012); and control of marine waters from entering the estuary, creating salinity gradients (Wetz et al., 2013). Human activities also affect estuarine environments and estuarine dynamics (Dutertre et al., 2013;Wetz et al., 2013;Gao et al., 2019;Hillman et al., 2020).
Our results showed that the ES of the Rongjiang estuary was higher in summer than in winter, and both anthropogenic pollutants (DIN) and natural variables (temperature of bottom water) play a crucial role in determining the ES of this area. The growth and development of benthic communities are known to be significantly affected by water temperature (Day et al., 1989;Vaquer-Sunyer and Duarte 2011;Spreitzenbarth and Jeffs, 2021). The reproductive activity of species (e.g., Limnodrilus hoffmeisteri and Branchiura sowerbyi) were high, and the proportion peaks during the winter (Li et al., 2012). Additionally, the study area has a mild, wet subtropical climate with northeasterly (winter) and south-westerly (summer) monsoons. Given that winter is characterized by low minimum rainfall and river flow (Chen 2010), the capacity of the Rongjiang River to dilute incoming substances is almost six times lower in winter than in summer, and this facilitates the accumulation of pollutants; this also explains why the concentrations of DIN were three times higher in winter than in summer (Huang et al., 2015). Consequently, the ES was poorer in winter. In winter, the dominance of larger and longer-lived species such as the polychaete Micronephthys oligobranchia decreased, and small and pollution-resistant species with short life spans (e.g., the polychaete Heteromastus filiformis and the oligochaetes L. hoffmeisteri and B. sowerbyi) which were assigned to EG IV and EG V, became dominant (Zhen et al., 1985;Li et al., 2012;Lin et al., 2013). This is a common feature of disturbance in estuaries that have experienced the loss of large macrofaunal animals (Pearson and Rosenberg 1978;Hillman et al., 2020).
Overall, the ES of the aquaculture zones, the site near the Hanjiang river water channel, areas near Shantou City on the Rongjiang River, and the site closest to the ocean was poor or bad. The pond aquaculture surrounding the Rongjiang River results in the rapid deterioration of the aquatic ecosystems (Gu et al., 2018) because the discharge of aquaculture wastewater contributes to nutrient enrichment in this area (He et al., 2005). The low ES of the site near the Hanjiang River water channel and in areas near Shantou City on the Rongjiang River were also due to man-made wastewater discharge (DIN concentration were still high although slightly lower than the aquaculture zones) and all samples met the class III criteria based on the Chinese Sea Water Quality Standard (0.4 mg/L), which is used for harbors (SEPA 1997). High concentrations of nutrients result in eutrophication; in eutrophic regions, the survival of sensitive benthic species is restricted, and they are gradually replaced by tolerant species. The accumulation of nutrients can provide abundant food sources for tolerant species, which allows them to reach high densities, and this increase the unevenness of benthic assemblages, decrease species diversity, and simplifying food webs (Alexander et al., 2017;Ji et al., 2020). Our results were consistent with this observation: species diversity was lower in areas with higher DIN concentration, where species (P. laevis) tolerant of excess organic matter enrichment and second-order opportunistic species (Borja 2000) represented more than 50% of taxa. In addition to eutrophication, the poor ES of the lower Rongjiang estuary was driven by intense hydrodynamic disturbance, which can limit the initial colonization of species (Haase et al., 2012;Van der Linden et al., 2017). This is consistent with the hypothesis proposed by previous studies, that hydrodynamics is a stress to most macrofaunal species in estuaries ( Van der Linden et al., 2012;Van der Linden et al., 2017).

Value of Using Multiple Indices to Assess the ES of Benthic Ecosystems and Recommendations for Environmental Protection
The composition and distribution of macrobenthic communities in estuaries are shaped by high variability in the physicochemical environment (Day et al., 1989;Leonard et al., 1998;Vaquer-Sunyer and Duarte 2011;Haase et al., 2012;Dutertre et al., 2013). These communities are characterized by the high dominance of a few species (Elliott and Quintino 2007). A single ecological indicator is thus often of limited use in ES assessments. For instance, the robustness of results based on AMBI and M-AMBI is reduced if only 1-3 taxa or individuals can be identified in a sample, or more than 20% of taxa in a sample are not assigned to an EG (Borja and Muxika 2005;Muxika et al., 2007). In addition, in case of the taxonomic diversity indices (TD) indices-average taxonomic distinctness Δ + and variations in taxonomic distinctness Λ + , too much interference information is added by tolerant taxa, such as some tolerant Oligochaeta and Mollusca, because these species contribute the greatest path lengths, which lead to equal or even higher scores in disturbed sites than in undisturbed sites (Hu and Zhang 2016). In this study, 66.7% of taxa were assigned to EG II at S4 in winter, and the ES of S4 was classified by AMBI and M-AMBI as good; however, only two taxa were identified. Healthy ecosystems are not dominated by a single functional group or by a few species (Peng, 2013). Classification of the ES by AMBI and M-AMBI indices was relatively accurate. At S7 in winter, more than 90% of taxa were assigned to EG III, EG IV, and EG V, and the ES of S7 was classified as good by the TD index. Thus, the use of multiple biotic indices can provide a more effective evaluation of the ES of certain areas, and especially in estuaries, where physicochemical environments are highly variable. We found that the AMBI, M-AMBI, Δ + , and Λ + indices were suitable for assessing the ES of benthic ecosystems in aquatic areas experiencing anthropogenic disturbance, and that the use of these indices may enhance our understanding of the effects of environmental variables on the distribution of species.
These results indicate that some measures need to be implemented, such as strict national and local standards to control the discharge of terrestrial pollutants and sewage into Rongjiang, Lianjiang, and Hanjiang Rivers. There is also a need to improve the arrangement of cages, implement controls on the scale of breeding, and improve the breeding structure in multi-species operating models to promote the healthy, stable, and sustainable development of aquaculture.

CONCLUSION
The main conclusions of this study are outlined below: 1) The ES of the benthic ecosystem in Rongjiang estuary has been seriously affected by anthropogenic activities. Several measures need to be implemented to prevent the loss of benthic habitat.
2) The ES of Rongjiang estuary differed in summer and winter.
3) The use of multiple ecological indices permits more robust assessments of the ES of benthic ecosystems in anthropogenically disturbed estuary systems.