Assessment of Benthic Ecological Quality Status Using Multi-Biotic Indices Based on Macrofaunal Assemblages in a Semi-Enclosed Bay

Semi-enclosed bays have physical and chemical characteristics influenced by both land and sea systems and the quality of the benthic environment is always of great concern. Macrofauna are considered good indicators for evaluating the benthic ecological quality status owing to their biological characteristics. In this study, six biotic indices, namely the Shannon–Wiener diversity index (H′), Abundance-Biomass Comparison (ABC) curve, AZTI’s Marine Biotic Index (AMBI), multivariate-AMBI (M-AMBI), BOPA index, and BENTIX index, were used to evaluate the adaptability of different biological indices in the bioassessment of the benthic environment in a semi-enclosed bay. In the annual environmental assessment of the study area, the average values of the six indices (H′, ABC curve, AMBI, M-AMBI, BOPA, and BENTIX) were 4.494, 0.182, 2.433, 0.791, 0.033, and 3.813, respectively; accordingly, H′, M-AMBI, and BOPA met the high standards whereas the other indices met the good standards, indicating that the whole study bay was slightly disturbed and had good ecological quality. From the perspective of spatial variation, the benthic environment in the middle of the bay was better than that in the north; the environmental problems in the northeast were particularly noteworthy. In terms of temporal patterns, the benthic environment in winter and summer was significantly better than that in spring and autumn, with obvious seasonal differences. The present results indicate that the H′ and ABC curve based on benthic abundance and biomass should be avoided for environmental assessment in mariculture areas. AMBI and M-AMBI should be used with caution when the percentage of unassigned species is high, in which case H′ is the appropriate choice. When there are few unassigned species, M-AMBI is more conducive for accurate evaluation of the benthic environment than AMBI and H′.


INTRODUCTION
Compared with the pelagic ecosystems, benthic ecosystems have high species diversity. The species number of marine benthic invertebrates accounts for 90% in the total marine invertebrates (Xu, 2011). Marine benthic organisms can be divided into macrofauna, meiofauna, and microfauna, according to their body size. As vital components of the benthic food web, they play important roles through complex trophic relationships in the material cycle and energy flow of benthic ecosystems. Macrofaunal communities contribute to promote biogeochemical processes by altering the roughness, original structure, and chemical properties of sediments through activities such as feeding, digging, and tube-building (Graf and Rosenberg, 1997;Snelgrove, 1998). The diversity of macrofaunal assemblages is also an indispensable indicator for evaluating the benthic ecological quality status . Most of macrofauna live on the bottom surface or inside of sediment and usually have weak migration capacities (Pelletier et al., 2010;Ji et al., 2016;Feng et al., 2021). In addition, different species of macrofauna differ in their adaptability to environmental conditions, tolerance and sensitivity to adverse factors such as pollution (Snelgrove, 1998). Macrofauna have often been used to indicate the quality of benthic environments (Pinto et al., 2008;Li et al., 2020) as the relationship between macrofauna communities and the degree of disturbance has been analyzed by previous studies (Zheng et al., 2011;Egres et al., 2019;Yoon et al., 2021).
The negative effects of industry, tourism, maritime transportation, fisheries, and aquaculture on biodiversity are noteworthy in coastal ecosystems (Simboura and Zenetos, 2002). Frequent human activities lead to heavy metal pollution, persistent organic pollutants, eutrophication, seawater acidification, and other issues, changing the structure of marine biological communities and marine biogeochemical cycles and ultimately affecting the structure and function of marine ecosystems (Lv et al., 2016). For example, the reclamation project in Bohai Bay and the expansion of Tianjin Port have permanently occupied coastal tidal flats, weakening the selfpurification capacity and tidal current hydrodynamic force and causing sediment deposition and topographic change, which have a serious impact on the coastal ecosystem (Nie and Tao, 2008;Zhang et al., 2015;Xu et al., 2021). Increasingly frequent human activities in the Yangtze River and adjacent waters have led to substantial increases in the concentrations of nitrogen and phosphorus in seawater, causing marine pollution events such as red tides (Huang et al., 2003;Liu et al., 2013). Many restaurants and entertainment facilities built in the Zhoushan Islands in the East China Sea to develop tourism have discharged untreated sewage into seawater and accelerated the degradation of the marine benthic environment (Zhang et al., 2017). Overfishing in coastal systems, especially bottom trawling, also has a significant impact on benthic communities, resulting in a decrease in community diversity .
In coastal systems, a semi-enclosed bay has wide water area and small waves, so it is an ideal place for development of aquaculture and port shipping. However, owing to the integration of the physical and chemical properties of land and ocean systems (Gao et al., 2021), the ecological environment has potential instability. Environmental problems such as small tidal volume, poor water exchange capacity, and significant influence of human activities have become increasingly noteworthy (Lim et al., 2012;Peng et al., 2013). Jiaozhou Bay is a typical semi-closed bay located on the southern coast of the Shandong Peninsula and is surrounded by densely populated urban areas and coastal industrial areas (Song et al., 2020). The inner bay has wide harbors, deep waters, and calm winds and waves, and the sea is not frozen all year round. It is a natural bay with many rivers flowing into it. The dual effects of frequent human activities and excellent natural conditions make it an ideal area for testing the adaptability of biological indices to semi-enclosed bays.
Shannon-Wiener diversity index (H ), Abundance-Biomass Comparison (ABC) curve, AZTI's Marine Biotic Index (AMBI), multivariate-AMBI (M-AMBI), BOPA index, and BENTIX were used to evaluate the spatial and temporal variation of benthic ecological quality in Jiaozhou Bay. The results of the assessment are also combined with environmental parameters to analyze the environmental quality of Jiaozhou Bay and further explore the impact of anthropogenic activities on the benthic environment. At the same time, taking Jiaozhou Bay as an example, the applicability of different biotic indices in environmental assessment of a semi-enclosed bay was explored and compared. The main advantages and constraints of each indicator are briefly summarized.
Considering the principles of the six indicators, the following assumptions are proposed: (1) The indexes based on biomass have limitations in shellfish farming areas.
(2) M-AMBI is calculated by combining AMBI and H , so it can better explain the status of environmental quality.

Field Sampling
Jiaozhou Bay is a 32 km long and 27 km wide semi-enclosed bay, located in the south of the Shandong Peninsula (35 • 38 -36 • 18 N, 120 • 04 -120 • 23 E). It is a temperate monsoonal climate zone. The tide is a typical semi-diurnal tide. The deposits of Jiaozhou Bay are terrigeneous and organic matter content is 0.38-1.91% which is increasing with time (Wang et al., 2017). In the early 1980s, Jiaozhou Bay was designated as a pilot area for the research and development of marine resources. There are abundant species of macrofauna in the bay and is one of the earliest typical areas to carry out bay ecology research in China (Yu et al., 2006;Yang et al., 2016). With the rapid development of coastal farming and reclamation has resulted in the gradual transformation of most of the coastline from natural to artificial areas (Fu et al., 2018;Quan et al., 2020b;Yin et al., 2021), leading to a decrease in tidal intake and a weakened self-purification ability . In addition, the water quality of Jiaozhou Bay is also affected by industrial and domestic wastewater from many surrounding rivers, including some nutrients, heavy metals, and some organic pollutants, as well as some petroleum and reactive phosphate pollutants from ship operations and fishing vessels (Melet et al., 2020;Gu et al., 2021a). Collectively, all factors put great pressure on the benthic environment of Jiaozhou Bay, causing different degrees of disturbance to benthic organisms and the benthic environment.
Macrofauna were collected from 14 sampling sites (J1-J14) in Jiaozhou Bay during four cruises in February, May, August, and November 2014, representing winter, spring, summer, and autumn cruises, respectively (Figure 1). The sites of each cruise were set consistently, but the J6 site sample was not collected during the spring and autumn cruises. Sites J1-J7 are located near Jiaozhou Bay Bridge, which was built in 2006 and completed in 2011. J2 and J6 are near the Licun River and Dagu River, respectively, which are affected by the impact of the river. J2 and J3 are located in shellfish aquaculture area. J13 and J14 are located at the junction of Jiaozhou Bay and the Yellow Sea. The 14 sites were evenly distributed in most areas of Jiaozhou Bay, reflecting the status of macrofauna and the environment. Sediment samples were collected and washed on deck using a 0.05 m 2 box corer, and four replicate samples were combined into one sample. These samples were subsequently preserved in 75% ethanol until laboratory identification macrofauna analysis.
To characterize the abiotic environment, the water depth, bottom water temperature, salinity, and pH at all stations were measured using a YSI 600XLM Multi-Parameter Water Quality Sonde (YSI Incorporated, United States) in situ. At the same time, surface sediments were sampled for the analysis of sediment water content (W), organic matter content (OM), and median grain size (Md), chlorophyll a content (Chl-a) and pheophorbide content (Pha). They were stored in −20 • C refrigerator until analysis in the laboratory.

Laboratory Analysis
In the laboratory, the sediment samples were stained with Rose Bengal overnight; then, the samples were sieved with a 0.5 mm mesh sieve to collect all macrofaunal organisms. All animals were sorted and identified to the lowest possible taxonomic level using a stereo microscope. Density (ind./m 2 ), and biomass (wet weight, g/m 2 ) were determined using a 0.1 mg precision electric balance (Eleftherious and McIntyre, 2013). The sediment grain size was measured using a laser particle size analyzer (Master Sizer 3000, China), and Chl-a and Pha were determined by F96Pro Fluorescence Spectrophotometry (Shanghai Lengguang Technology Co., Ltd., China).

Statistical Analysis
Microsoft Excel 2019 (Microsoft, United States), PRIMER 6.0 software package (PRIMER-E, United Kingdom), SPSS 22.0 (SPSS, United States), and ArcGIS 10.6 (Esri, United States) were used for data analysis and geographic mapping. The calculation processes of indices are in Supplementary Material.

Evaluation Standard
The European Water Framework Directive (WFD) issued by the European Union in 2000 recommended that the status of ecological environmental quality be divided into five grades; hence, all the evaluation indices in this study eventually divided the status of ecological environmental quality into five grades: high, good, moderate, poor, and bad ( Table 1).

Correlation Analysis
The correlation analysis between the biological index (H , W, AMBI, M-AMBI, BOPA, and BENTIX) and physical and chemical factors was used to test the response of the biological index to human environmental pressure. Pearson correlation analysis was used for correlation significance, which was performed using SPSS 22.0.

Community Structure
A total of 251 species of macrofauna belonging to 10 phyla were identified. The mean macrofaunal abundance and biomass values were 2259.39 ind./m 2 and 501.88 g/m 2 , respectively. Polychetes were the most dominant group, accounting for 41% of the total species and 52% of the total abundance, followed by crustaceans, mollusks, and echinoderms.

Spatial Distribution in Winter
A total of 101 species of macrofauna were recorded during winter cruises. H varied in the range of 3.284-4.406 at 14 sites. The lowest H value was found at site J7, and the highest value occurred at J10. A slight disturbance occurred at J1, J3, J4, J5, J7, J13, and J14. The values of the remaining seven sites were greater than 4.0, and the bioassessment indicated high status. The range of W values generated by ABC curve analysis was 0.220-0.430. The lowest value appeared at site J14 in the east of Jiaozhou Bay, and the highest value appeared at site J11 in the west. As the W value of each site was distributed in the range of 0.15-0.49, each site was slightly disturbed and the environmental assessment result was good.
The range of AMBI values was 1.589-2.585. The lowest value was observed at site J11, and the highest value was at site J14. All the sampling sites were in good status. The range of M-AMBI was 0.680-0.880, with the lowest value at site J1 in the northeastern part of Jiaozhou Bay and the highest value at site J10 in the east. There were seven sites with bioassessment as high status, accounting for 50% of the total sites.
The average values of BOPA in the four seasons were in the range of 0.001-0.036, with the lowest value at site J4 in the northern part and the highest value at site J1. All sites were of a high status.
The highest BENTIX values occurred at site J4 in the north, whereas the lowest values occurred at site J8 in the middle of the bay. Sites J4, J11, and J12 were of a high status, whereas the rest were of a good status.
All sites in the winter cruise were slightly disturbed or undisturbed, and all types of environmental assessments were in good-high status, with relatively uniform evaluation. As shown in Figure 2, H and M-AMBI showed that the environments in central and southern Jiaozhou Bay were better than those in the north. The AMBI and BNETIX bioassessments were similar, and the benthic environments of the J4 site and southwest were better than other regions.

Spatial Distribution in Spring
A total of 141 species were identified during the spring cruise. H varied from 1.688 to 4.855 at 14 sites. The lowest value was found at the J12 site in the south, and the bioassessment result indicated a poor status. The highest value was observed at the J10 site in the east. Sites J2, J3, and J5 were of moderate status, whereas sites J1, J4, J7, and J11 were of good status. Meanwhile, the bioassessment result of J8, J9, J10, J13, and J14 in the central and eastern part indicated a high status.
The range of the W value was 0.069-0.293. The minimum value appeared at the J12 site, and the maximum value was observed at the J1 site in the northeastern direction. Sites J2, J3, and J12 had moderate disturbances. The environmental assessments at the other sites were rated as good.
The range of AMBI values was 2.001-2.660, with the lowest value at site J11 and the highest value at site J2. All sites were rated as having good status. The M-AMBI ranged from 0.510 to 0.910, with the lowest value located at site J12 and the highest at site J10 in the eastern part. Only site J10 was undisturbed and rated at a high status. Site J12, located in the southern part, had the worst environmental rating, which was consistent with the results of H and ABC curve. All sites had an environmental rating of good except for site J12.
The BOPA ranged from 0.002 to 0.045, with the lowest value located at site J2 in the northeast of Jiaozhou Bay and the highest value at site J9. All sites were of a high status.
The range of BENTIX was 2.913-4.439, with the lowest value at the J2 site and the highest value at the J7 site in the western part of the bay. The benthic environment of the J1 and J2 sites was moderately disturbed, and the other sites were evaluated as having a good status.
During the spring cruise, H , ABC curve, AMBI, and M-AMBI showed that some sites were moderately and heavily disturbed, mainly at sites J2, J3, and J5 in the north and J12 in the south of Jiaozhou Bay. The central part is slightly disturbed. As shown in Figure 3, various indices indicated that benthic environment of the central part was better than that of the northeastern part in spring, whereas BOPA indicated that the central benthic environment was worse than other parts.

Spatial Distribution in Summer
A total of 162 species were identified during summer cruises. H ranged from 1.531 to 4.781, with the lowest value at site J2 in the northeast, which was environmentally assessed as poor status, and the highest value at site J8 in the center. J3 and J13 were environmentally assessed as having a good status, and the rest of the sites were assessed as having a high status.
The W values ranged from 0.075 to 0.447. The lowest value was at site J2, and the highest was at site J8. Sites J2 and J12 were rated as moderate, and the rest of the sites were rated as good.
The range of AMBI values was 1.747-2.895, with the lowest value at site J7 and the highest value at site J13, all of which were rated as good. M-AMBI ranged from 0.510 to 0.890. The lowest value was located at site J2 with an environmental rating of moderate, which is consistent with the worst-rated site in the  Frontiers in Marine Science | www.frontiersin.org H , ABC curve, and BENTIX. The central region was assessed as having a high status.
The BOPA range was 0.001-0.063, with the lowest value at site J2 and the highest value at site J13. All the sites were assessed as having a high status. Three sites -J5, J11, and J13 -had an environmental rating of good, whereas the rest of the sites had an environmental rating of high.
The BENTIX range was 2.580-5.053, where site J7 was the least disturbed and had the best environmental rating. Site J2 was moderately disturbed and had a moderate environmental rating. The highest and lowest values are consistent with those of the summer cruise. The J3, J7, J12, and J14 sites had a high environmental rating, and the rest of the sites had a good environmental rating.
As shown in Figure 4, during the summer cruise, different indices generally rated the benthic environment as high. However, various indices showed moderate to heavy disturbance at site J2, with only BOPA showing no disturbance at this site. Some of the sites in the southern part of Jiaozhou Bay were also lightly disturbed, which made the benthic environment worse than that of the central region.

Spatial Distribution in Autumn
A total of 144 species were identified during autumn cruise. H ranged from 3.607 to 4.753, with the lowest value at site J5 and the highest value at site J9. Sites J1, J2, J4, J5, J11, and J14 were assessed as having a good status, whereas the rest were assessed as having a high status.
The W values ranged from 0.045 to 0.349. The lowest value was observed at site J5, and the highest was at site J12. Meanwhile, sites J5 and J14 were rated as moderate, and the rest of the sites were rated as good.
The range of AMBI values was 2.345-3.285, with the lowest value at site J11 and the highest value at site J1. All sites had an environmental rating of good. The M-AMBI ranged from 0.670 to 0.840, with the lowest values located at site J1 and the highest values at sites J2, J8, and J9. The environmental assessment was high, except for sites J1, J3, J4, J5, and J14.
The BOPA range was 0.034-0.101, with the lowest value at site J2 and the highest value at site J1. Sites J2, J3, and J9 were rated as high, and the rest of the sites were rated as good.
The BENTIX range was 2.629-4.390, with the lowest value at site J1 and the highest value at site J11. Seven sites had moderate contamination, namely J1, J2, J3, J4, J5, J8, and J14, whereas the rest of the sites had an environmental rating of good.
The BENTIX value was poor for the autumn cruise, with over 50% of the sites being moderately disturbed, mostly in the northern part of Jiaozhou Bay. As shown in Figure 5, multiple indices showed that site J1 was the most disturbed, followed by site J5. Overall, the environmental assessment of the central and southern parts of Jiaozhou Bay was generally better than that of the northern part.

Season Patterns
Among the four cruises of macrofaunal surveys in 2014, the largest number of species (162 species) were surveyed in  respectively. In terms of biomass, the results of the survey were winter (164.62 g/m 2 ) < spring (524.94 g/m 2 ) < summer (572.08 g/m 2 ) < autumn (814.18 g/m 2 ). Ruditapes philippinarum was the dominant species in the biomass.
As shown in Figure 6, in terms of H , the environmental assessment of the spring cruise was the most unstable, with a wide range of variation. The environmental quality of the autumn cruise was significantly higher than that of the spring cruise. With regard to the ABC curve, the winter cruise had better environmental quality, but there was no significant difference between the W values of the four cruises. AMBI rated the benthic environment as good in all four seasons, but the environmental quality of the autumn cruise was still significantly lower than that of the other cruises. With regard to the M-AMBI value, the environmental quality of the spring cruise was significantly lower than that of the other cruises. Only one site on the spring cruise reached a high standard, which was lower than that of the other cruises. In terms of BOPA, the environmental quality of the autumn cruise was significantly lower than that of the other cruises whereas the environmental quality of all sites during the winter and spring cruises was high; only three sites on the summer cruise were rated as good. For BENTIX, the environmental quality of the autumn cruise was significantly lower than that of the other cruises, with 53.8% of sites rated as moderate.
As shown in Figure 6, H , W, and M-AMBI all showed the worst environmental assessment for the spring cruise. BOPA, AMBI, and BENTIX showed the worst environmental assessment for the autumn cruise. A one-way ANOVA showed that there was a significant difference between the autumn and spring cruises of H . There was no significant difference between the W values of the four cruises. M-AMBI of the autumn cruise was significantly different from that of the other cruises. Among AMBI, BOPA, and BENTIX, the bioassessment of autumn cruises was significantly different from that of the other cruises.

Limitations of the Biotic Indices Based on AZTI's Marine Biotic Index Species List
BOPA, AMBI, M-AMBI, and BENTIX were calculated based on taxon sensitivity classification. BOPA is based on the proportion of polychetes and amphipods, and the other three are referenced to the AMBI species list. The calculation and analysis of the index will be better understood by comparing the composition of EGI-EGV groups in different sites.
Based on the AMBI species list, the distribution of different biological groups at each sampling site was analyzed. As shown in Figure 7, the proportions of the EGII and EGIII groups were similar and high in most sites. In contrast, the  proportion of species at the EGIV and EGV groups was extremely low. A significantly higher proportion of the EGIII group occurred at site J2 during the spring, summer, and autumn cruises. R. philippinarum was the dominant species at this site. The species abundance of unassigned organisms in the spring and summer cruises was relatively high, especially at site J12 (spring). Hemileucon bidentatus was the main unassigned species, and its genus (Hemileucon sp.) was not assigned to the AMBI species list.

Correlation Analysis of Biological Index
The correlation analysis results of each biological index and environmental parameters ( Table 2) showed that H had a significant positive correlation with pH, and W had a significant positive correlation with sediment organic matter content. AMBI was significantly correlated with bottom water salinity, pH, Chl-a, and sediment organic matter content. M-AMBI had a significant positive correlation with pH. BOPA had a highly significant correlation with bottom water salinity. In this study, BENTIX showed no correlation with environmental parameters. Water depth, temperature, Pha content, and median grain size also did not have a significant impact on the biological index.
Correlation analysis results between the six indices ( Table 3) showed that H and W values, M-AMBI, and BENTIX had a significant positive correlation; W had a highly significant positive correlation with M-AMBI, and a significant negative correlation with AMBI. AMBI had a highly significant negative correlation with M-AMBI and BENTIX, and a highly significant positive correlation with BOPA; and BENTIX had a highly significant positive correlation with M-AMBI and BOPA. Overall, M-AMBI had the highest correlation with the other indices. Among them, it had an extremely significant correlation with the four indices and only had no significant correlation with BOPA.

Evaluation of Different Biological Indices
Biological monitoring can comprehensively reflect the interaction between various environmental pollution factors and can continuously monitor ecological health (Wang et al., 2010;Markert et al., 2013). The results of the multiple biotic assessment methods showed that Jiaozhou Bay had a good benthic ecological quality status, which was consistent with the results of other studies (Yang et al., 2021;Yin et al., 2021).
Shellfish aquaculture activities have expanded in the bay since the 1980s, making the region an important shellfish production base in China (Gu et al., 2021b). R. philippinarum has been cultivated in Jiaozhou Bay since 2003 (Fan and Liu, 2016), becoming the dominant species of macrofauna until now (Wang et al., 2011). For nearly 150 years, Jiaozhou Bay area has been reduced by 235.41 km 2 during 1863-2012 (Ma et al., 2014), inevitably causing a sharp decline in water exchangeability in the bay. Industrial and domestic sources along the north coast led to significant increases in N and P in the water (Sun et al., 2011;Gao et al., 2018), but the diffusion and attenuation of these nutrients were impeded by the descending hydrodynamic force and the deterioration of water quality in the semi-enclosed bay . Spatially, Qingdao ports (Qingdao Port and Huangdao Port) are mainly located in the eastern and western waters of Jiaozhou Bay, and pollutants brought by wharf construction and a large number of ships make the concentrations of Cu, Pb, Zn, Hg, and Cr higher in the east of Jiaozhou Bay and the western part of the bay mouth (He et al., 2013). In addition, industrial wastewater mainly flowed into the Bay from river mouth, and quantity of pollutant discharge rapidly from about 7.0 × 10 7 t/a in the 1980s to 1.0 × 10 8 t/a in the 1990s, and then maintained at about 9.8 × 10 7 t/a (He et al., 2013). The estuaries of Dagu River, Baisha River, Licun River, and Haibo River are mainly distributed in the northern and eastern part of Jiaozhou Bay. Among them, COD emissions of Dagu River accounted for about 40.4% of Jiaozhou Bay. The Haibo River located in the east of the bay mouth had the largest discharge of inorganic nitrogen and phosphate, accounting for more than 35% of the total discharge of Jiaozhou Bay, followed by the Dagu River and Licun River (Zhang and Sun, 2007). In this study, the northeastern and southern parts of the bay were moderately disturbed by varying degrees, particularly at sites J2 and J12. Site J2 was in a mariculture area with an extremely high abundance of R. philippinarum, while the abundance of this species was low at other sites, which may cause the poor benthic environmental quality of the site J2 and J12 according to the evaluation. At the same time, it is located in the Licun estuary, and the river carries land-based pollutants into the sea, which may have an impact on the benthic environment. Studies have shown that the Licun and Dagu estuaries were enriched in Cr, Ni, and Cu heavy metals (Xiao et al., 2017), and the benthic environment of the northern waters of Jiaozhou Bay and Licun estuary were more seriously disturbed (Cui et al., 2020;Yang et al., 2021), which is consistent with the results  of this study. The poor benthic environment in the northern areas may also be related to the construction of the Jiaozhou Bay Bridge. Site J12, located in the southern part, showed severe and moderate disturbance in the environmental assessment in spring and summer, which may have been caused by the high abundance of H. bidentatus in spring (2980 ind./m 2 ) and summer (1140 ind./m 2 ), accounting for 78.11% and 29.84% of the total abundance, respectively, which had some impact on the environmental assessment. All the results of the six biotic indices of the central sites (J7, J9, and J10) were in good high status, and the evaluation results were consistent, showing that the benthic environment in the central site was relatively stable, and the living environment of macrofauna was not affected by strong human activities.
Since the 1980s, the intertidal area of Jiaozhou Bay has been decreasing due to the influence of human activities, such as the reclamation of the sea and the development of saltern (Lei et al., 2013). The water exchange capacity of each region has gradually reduced, resulting in serious eutrophication of seawater (Yang et al., 2018;Yuan et al., 2020). At the beginning of the 21st century, with the implementation of a series of protection measures and the abandonment of salt pans and aquaculture ponds (Yang et al., 2018), the water area of Jiaozhou Bay gradually rebounded, while environmental pollution was also reduced (Lei et al., 2013;Gu et al., 2021a;Yang et al., 2021). In recent years, the number of macrobenthic species has generally shown an increasing trend (Yang et al., 2021), and the ecosystem functioning has remained stable. Temporally, the benthic environment was less favorable in spring and autumn than in summer and winter. The higher abundance of R. philippinarum and H. bidentatus at J2, J3, J11, and J12 in the spring cruise was probably the main factor that pulled down the overall benthic environment assessment. In the autumn cruise, the mean abundance of Mediomastus sp. was relatively high, with an average abundance of 507.69 ind./m 2 . At the same time, Mediomastus sp. belonged to the EGIII class of pollution-tolerant species. Therefore, the indices based on the benthic sensitivity classification showed that the environmental assessment of the autumn cruise was poor.

Effects of Environmental Factors on Macrofauna
Changes in habitat directly affect the species composition and community structure of macroinvertebrates. Many studies have shown that the spatial distribution of macrofauna is closely related to water temperature, water depth, hydrodynamic conditions, dissolved oxygen content, and sediment types (Mancinelli et al., 1998;Ysebaert and Herman, 2002;Como and Magni, 2009;Yang et al., 2021). For example, in the survey conducted by Quan et al. (2020a), temperature, nutrients, and sediment types were the environmental factors that mainly affected the distribution characteristics of macrofauna. At the same time, heavy metal pollution can accumulate in the viscera and tissues of benthic organisms, hindering their growth and reproduction (Solà and Prat, 2006). High levels of total phosphorus and nitrogen can lead to the loss of macrobenthic species and an increase in the biomass of tolerant species (Cai et al., 2014). Eutrophication is also associated with ocean surface hypoxia and seawater acidification, which can increase the sensitivity of coastal zone seawater to ocean acidification (Lv et al., 2016;Gao et al., 2020).
The effect of temperature on macrofaunal diversity was also significant in this study, with winter macrofaunal diversity being significantly lower than in other seasons. A number of investigations have concluded that macrofaunal abundance has a highly significant positive correlation with bottom temperature (Shojaei et al., 2015;Neumann et al., 2021), especially in semienclosed bays, which may be related to the drastic changes in bottom temperature (Yang et al., 2021). Multi-biotic indices were significantly or highly significantly correlated with pH, organic matter, sediment salinity, and Chl-a content, suggesting that the marine benthic environment is related to pH, and nutrients, which is similar to other studies (Oleszczuk et al., 2021;Toussaint et al., 2021). In the study by Neumann et al. (2021), bioturbation and temperature had the highest explanatory power to explain the substantial seasonal variability in observed oxygen and nutrients. The total organic matter content is an important environmental variable for anoxic mineralization (Toussaint et al., 2021). Low oxygen conditions can change predator-prey relationships and destroy benthic habitats (Briggs et al., 2017).

Comparisons Among Different Biotic Indices
Because of the complexity of the marine ecosystem and the difference in reference values of different evaluation indices, different evaluation indices usually give different evaluation results in the same region (Medeiros et al., 2012;Wu et al., 2013). In this study, the benthos sensitivity classification indices (AMBI, BOPA, and BENTIX), the benthic biodiversity indices (H and ABC curve), and the composite index (M-AMBI) combined with multiple indices were mutually calibrated. H and ABC curve rely on abundance and biomass relationships and do not consider taxon characteristics. AMBI, M-AMBI, and BENTIX are three types of indices needed to classify the species, which is a tedious and time-consuming operation (dela-Ossa-Carretero et al., 2009). In particular, some endemic species are not included in the classification list; therefore, it is necessary to find the same genus for classification judgment, which has certain subjectivity. For example, in this study, the species of H. bidentatus were highly abundant, but they were not assigned to AMBI and M-AMBI. In addition, AMBI presents some limitations when applied to semi-enclosed systems. The relative frequencies of opportunistic polychetes and amphipods were used for BOPA environmental assessment. Opportunistic polychetes are tolerant to organic-rich sediments, while amphipods are sensitive to changes in the organic matter content (de-la-Ossa-Carretero et al., 2009). The complexity of the classification operation was reduced in this index, but the influence of other organisms on the benthic environment was not considered. Some studies have pointed out that under low disturbance conditions, amphipods are susceptible to other environmental factors, such as temperature, and low disturbance levels cannot cause polychete opportunistic species to become the dominant species (Wang et al., 2017). Therefore, the BOPA is not suitable for evaluating environmental quality when the disturbance is relatively light. BENTIX is considered an ecologically relevant biological index, which is easier to operate than AMBI and M-AMBI and has been successfully applied to areas of organic matter pollution, oil spills, and heavy metal pollution (Simboura and Zenetos, 2002;Zenetos et al., 2004;Simboura et al., 2007). It is more sensitive to the increase in organic matter content in sediments than AMBI (Caglar and Albayrak, 2012), but all opportunistic taxonomic groups (AMBI-EG IV and V) and sensitive taxonomic groups (AMBI-EG III) are given the same proportion in the evaluation; therefore, the ecological state is easily underestimated. At the same time, unassigned organisms remain neglected because of their reliance on the AMBI classification criteria.
In the present study, the abundance of R. philippinarum and H. bidentatus was extremely high at some sites. These sites (J2 and J12) were evaluated as moderately or severely disturbed according to the H and ABC curve, which differed from the results of other indices. It is obviously unreasonable to conclude that the benthic environment was severely disturbed by abundance and biomass alone in mariculture area (Cai et al., 2016), which verifies the first hypothesis. Therefore, it is recommended to avoid using the ABC curve alone for environmental evaluation in aquaculture areas. The disturbance caused by artificial culture for environmental evaluation should be fully considered. M-AMBI had a highly significant correlation with the other indices (expect BOPA) (p < 0.01). Compared with other indices in this study, M-AMBI is more accurate, which supports the second hypothesis. When using AMBI and M-AMBI, the abundance proportion of unassigned species should be fully considered.

CONCLUSION
Long-term monitoring and evaluation play an important role in the healthy development of benthic ecosystems, and the biological index of benthic ecological characteristics is generally not universally adopted. Therefore, the use of a single index for environmental assessment is one-sided and unstable, and the combination of various methods is more likely to obtain the overall environmental pollution in the survey areas. When the abundance of non-polychete opportunistic species and amphipods was high, BENTIX was better than BOPA in evaluating environmental quality. In mariculture areas, the extremely high abundance of a single species brings high biomass, and it is recommended to use the ABC curve with caution. Therefore, in the construction of aquatic ecological evaluation based on macrofauna, we recommend using the biotic indices with caution and, whenever possible, calculating several to reduce potential biases in the quality assessment. The classification of species tolerance should be integrated to weaken the objective difference between the evaluation methods of biological evaluation indices to obtain more comprehensive and accurate results for aquatic ecological evaluation.

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
XLu, JX, ZX, and XLiu contributed to the collection of samples, conception, design, collation, analysis of data, and preparation of the manuscript. XLu collated and analyzed the data and wrote the draft of the manuscript. JX contributed to the data reanalysis and the writing, revision of the manuscript. ZX contributed to the collection and identification of macrofauna and the collection of sediment samples. XLiu contributed to the experiment design, source acquisition, manuscript writing, revision, and approved the submitted version. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We appreciate the help of many members of the Laboratory of Benthos during field sampling, laboratory analysis, and manuscript preparation.