Sixty-five grid-based sampling stations were established and divided into three sectors (inner, middle, and outer) according to their relative geographic locations and environmental conditions (Rao et al., 2020a; Rao et al., 2021). Sampling surveys were conducted in the autumn (November 2017), winter (January 2018), spring (April 2018), and summer (July 2018). However, in this study, we only selected data from 36 sampling stations, namely, 12 sampling stations in the inner bay, 12 sampling stations in the middle of the bay, and 12 sampling stations outside the bay ( Figure 1 ).

Benthic macrofauna were collected using a 0.05 m² van Veen grab at each station. The sediments were washed through a 0.5 mm mesh sieve, and the residues were transferred to sample containers with 5% formalin buffer in situ for further identification. In the laboratory, benthic macrofauna were identified to the lowest possible taxon and enumerated under a dissecting microscope. They were then weighed using an electronic balance (0.1 mg, FA1204) after blotting surface water off with clean absorbing paper (Rao et al., 2020a; Rao et al., 2021).

There are many benthic animals in the Daya Bay subtidal zone (Jiang et al., 1990a; Jiang et al. 1990b; Du et al., 2008a; Du et al., 2008b; Du et al., 2009; Du et al., 2011; Rao et al., 2020a; Rao et al., 2021). In this study, ten species of benthic animals were considered based on the following three points: (1) according to the relative importance index (IRI) value (Zhang et al., 2017), based on four surveys in 2017 and 2018 that found Aglaophamus dibranchis, P. undulatus, L. brevirostris, A. laevis and P. cristata to be the top five species for average IRI values, and Sigambra hanaokai and T. lata to be the first seven and the first 20, respectively; (2) potential disaster-causing animals in the nuclear power coldsource system, including A. dibranchis, P. cristata, P. undulatus, L. brevirostris, Apionsoma (Apionsoma) trichocephalus, and P. bidentata (Cai et al., 2022b); (3) historically recorded dominant species, including T. scabra, P. undulatus, L. brevirostris, and A. laevis (Jiang et al., 1990a; Jiang et al., 1990b; Du et al., 2008a; Du et al., 2008b; Du et al., 2009; Du et al., 2011; Yuan et al., 2017; Rao et al., 2020a; Rao et al., 2021).

Water depth was determined using portable bathymetry (SPEEDTECH SM-5A, USA). Salinity and dissolved oxygen (DO) were measured using a portable water quality analyzer (WTW Multi 3430, Germany). Oils were extracted with n-hexane and measured using UV spectrophotometry procedures (Ehrhardt and Burns, 1993). The metal concentrations of Cr (chromium), As (arsenic), and Pb (lead) were determined using inductively coupled plasma mass spectrometry (Agilent 7700x, Agilent Technologies, USA). Quality assurance was performed using the procedure described by Wang et al. (2019). The samples for analyzing total nitrogen (TN) and total phosphorus (TP) were measured using the method of Varol and Şen (2012). Total nitrogen was measured by converting all nitrogen forms to nitrate via alkaline persulfate oxidation and subsequent analysis of nitrate was performed using spectrophotometric procedures. The total phosphorus was determined spectrophotometrically using the ascorbic acid method after persulfate digestion. The silicate content was measured using standard silicon molybdenum blue spectrophotometric procedures. Samples were filtered with 0.45 μm polycarbonate filters and analyzed later using the method of Dai et al. (2008). The samples for analyzing suspended solids (SS) were collected using 0.45 μm polycarbonate filters and were subsequently measured gravimetrically (Wang et al., 2015).

Variance analysis (ANOVA) and correlation analyses were performed using SPSS v25 software. ANOVA was used to determine whether there were significant differences in the ten species of benthic macrofauna and the 12 water environmental parameters in different seasons and regions. Correlation analysis was used to determine whether the densities and biomasses of the common benthic macrofauna were significantly correlated with the 12 water environment parameters. BIO-ENV analyses were performed using PRIMER v7 (Anderson et al., 2008). Similarities in benthic macrofauna between each pair of sites were determined using the Bray-Curtis similarity measure based on the fourth root transformed abundance data. Nonmetric multidimensional scaling (NMDS) ordination based on Bray-Curtis similarity was performed to explore the seasonal and site variation of the macrofaunal community. BIO-ENV analyses were used to examine the major environmental factors affecting the ten species of benthic macrofauna.

From the spatial distributions of ten species of benthic macrofauna in the subtidal zone near the Daya Bay Nuclear Power Plant, the densities of A. trichocephalus, A. laevis, and P. bidentata increased from inside the bay to outside the bay. The densities of P. cristata, T. lata, and T. scabra decreased from inside to outside the bay. The densities of A. dibranchis, and P. undulatus were higher in the middle bay than in the inner bay and outside the bay. The densities of S. hanaokai and L. brevirostris were lower in the middle bay than in the inner bay and outside the bay ( Figure 2 ). In addition, the density distribution of L. brevirostris was inconsistent with the biomass distribution, while the density and biomass distributions of the other nine macrobenthic species were consistent ( Figure 3 ).

Two-way ANOVA results indicated that the densities of A. trichocephalus, A. laevis, P. bidentata, P. cristata, A. dibranchis, and P. undulatus showed significant regional variation ( Table 1 ). The biomasses of A. trichocephalus, A. laevis, P. cristata, P. undulatus, and L. brevirostris showed significant regional variation ( Table 2 ). The densities of A. laevis, P. cristata, A. dibranchis, P. undulatus, and T. lata showed significant seasonal variation ( Table 1 ). The biomass of A. laevis, A. dibranchis, P. undulatus, and T. lata showed significant seasonal variation ( Table 2 ).

From the spatial distributions of 12 water environmental parameters in the subtidal zone near the Daya Bay Nuclear Power Plant, water depth, salinity, DO, suspended matter, Cr, and Pb increased from inside the bay to outside the bay. Oils, TP, and silicate decreased from the inside to the outside of the bay. The pH was higher in the central bay than in the inner bay and outside the bay. The TN was lower in the central bay than in the inner bay and outside the bay ( Figure 4 ).

A two-way ANOVA indicated that water depth, salinity, Cr, silicate, and pH showed significant regional variation. Salinity, DO, suspended matter, Cr, Pb, silicate, oils, TP, pH, and TN showed significant seasonal variation ( Table 3 ).

In terms of water environmental parameters, water depth was significantly correlated with the densities of six species of benthic macrofauna. There was a significant positive correlation between the densities of A. trichocephalus, A. laevis, P. bidentata, and A. dibranchis and water depth. There was a significant negative correlation between the densities of P. cristata and T. lata and water depth. Regarding benthic macrofaunal density, A. laevis and T. lata were significantly correlated with the five water environmental parameters. There was a significant positive correlation between the density of A. laevis and water depth, salinity, and Pb. There was a significant negative correlation between the densities of A. laevis, silicate, and TP. There was a significant positive correlation between the densities of T. lata, silicate, and TP. There was a significant negative correlation between the density of T. lata and water depth, salinity, and pH ( Table 4 ).

In terms of environmental parameters, water depth was significantly correlated with the biomass of six species of benthic macrofauna. There was a significant positive correlation between the biomass of A. trichocephalus, A. laevis, and P. bidentata and water depth. There was a significant negative correlation between the biomass of P. cristata, T. lata, S. hanaokai and water depth. Among the benthic macrofaunal biomass, A. laevis and T. lata were significantly correlated with the five water environment parameters. There was a significant positive correlation between the biomass of A. laevis and water depth, salinity, As, and Pb. A significant negative correlation was observed between the biomass of A. laevis and silicate. There was a significant positive correlation between the biomass of T. lata and silicate and TP. There was a significant negative correlation between the biomass of T. lata and water depth, salinity, and pH ( Table 5 ).

The BIO-ENV analysis showed that water depth was the first factor affecting 10 macrobenthic species, with a correlation coefficient of 0.191. The second factor affecting 10 macrobenthic species was salinity, with a correlation coefficient of water depth + salinity of 0.215.

This study involved three polychaetes: P. cristata, A. dibranchis, and S. hanaokai. Three species of polychaetes were the dominant species of macrozoobenthic communities in Daya Bay (Jiang et al., 1990a; Jiang et al., 1990b; Du et al., 2008a; Du et al., 2008b; Du et al., 2009; Du et al., 2011; Yuan et al., 2017; Cai et al., 2022b), and they were also the dominant species along the coast of China (Cai and Zheng, 1994; Huang et al., 2012; Chen et al., 2020; Rao et al., 2020b). Sigambra hanaokai is the dominant species in the macrobenthic community in the southern Yellow Sea (Xu et al., 2021).

P. cristata (It was previously identified as Paraprionospio pinnata in China and Korean, therefore, the P. pinnata reported in Chinese waters should be P. cristata) collected in Jiaozhou Bay and the Yellow Sea may be P. inaequibranchia and P. coora (Zhou et al., 2008). P. cristata has been widely reported in Korean waters (Yokoyama and Choi, 2010). P. cristata is considered an opportunistic species (Ji et al., 2022). In coastal waters along the Arabian Sea, Cossura coasta, the dominant species during the pre-monsoon period, was replaced by the surface deposit feeder P. pinnata during the monsoon and post-monsoon periods (Rehitha et al., 2019). Paraprionospio pinnata is only located in low-oxygen habitats in the Humboldt upwelling ecosystem (Fajardo et al., 2018) and it is the dominant species of the macrobenthic community in winter in the cage culture area of Daya Bay (Huang et al., 2005). The results of this study showed that the density of P. cristata was significantly negatively correlated with water depth and suspended solids, and significantly positively correlated with oil, whereas the biomass of P. cristata was significantly negatively correlated with water depth, salinity, and pH, and significantly positively correlated with oils and TP.

The relationship between S. hanaokai and the environmental factors in Chinese sea waters has been rarely reported, only a correlation analysis showed that there were significant negative correlations between the density of S. hanaokai and the content of organic carbon and nitrogen in the intertidal zone of Xiamen Crocodile Island (Cai et al., 2022c).

In winter and spring, the grappler method risk index (GMRI) of A. dibranchis is more than 50%, the risk of blocking the nuclear cold source system is at a medium-risk level. In summer and autumn, the GMRI of A. dibranchis is less than 50%, the risk of blocking the nuclear cold source system is at a medium-risk level. The risks of blocking the nuclear cold source system P. cristata and S. hanaokai are both at a low-risk level in Daya Bay in four seasons (Cai et al., 2022b).

Apionsoma trichocephalus was not the dominant species in previous surveys in Daya Bay (Jiang et al., 1990a; Jiang et al., 1990b; Du et al., 2008a; Du et al., 2008b; Du et al., 2009; Yuan et al., 2017), but it poses a potential risk of blocking nuclear cooling source systems (Cai et al., 2022b). Although the average density of A. trichocephala was low in all four seasons in Daya Bay (below 4 ind./m²), it was the dominant species in the subtidal zone of the East China Sea during all four seasons (Shou et al., 2018) and in the subtidal zone of the Chinese islands (Huang et al., 2012). However, there are few reports on the relationship between A. trichocephala and environmental factors. Our study found that A. trichocephala showed a significant positive correlation with water depth and a significant negative correlation with silicate.

L. brevirostris is one of the dominant species in macrozoobenthic communities in Daya Bay (Jiang et al., 1990a; Jiang et al., 1990b; Du et al., 2011; Yuan et al., 2017; Cai et al., 2022b). The results of this study showed that the density of L. brevirostris was significantly positively correlated with TN and TP. The density of L. brevirostris was high in sea areas with high organic matter content in Daya Bay (Jiang et al., 1990b).

Although the individual weight of L. brevirostris is large, its density is low and the distribution range is small; therefore, the risk of blocking the nuclear power cold source system is also low. Except for the GMRI in summer, which is over 50%, the risk of blocking is less than 50% in the other three seasons in Daya Bay (Cai et al., 2022b).

P. undulatus and T. scabra are the dominant species in macrozoobenthic communities in Daya Bay (Jiang et al., 1990a; Jiang et al., 1990b; Du et al., 2009; Du et al., 2011; Yuan et al., 2017; Cai et al., 2022a). The undulated surf clam, Paphia undulata (the new revised name in WoRMS is P. undulatus), is commercially cultured on the southern coast of China (Zhang et al., 2022a). Therefore, there have been many studies on the relationship between P. undulatus and environmental factors. P. undulatus showed a significant positive correlation with the DO, clay, and chlorophyll-a content in the Beibu Gulf (Ye et al., 2010). Mud substrates with ≥ 40% water content in the temperature range of 20–30°C and salinity range of 20–40 psu were appropriate for P. undulatus burrowing and may be appropriate for its culture (Zhang et al., 2022a). P. undulatus is not only a cultured species but also a common dominant species in the macrozoobenthic community along the southeast coast of China (Ye et al., 2010; Du et al., 2011; Shu et al., 2015; Yuan et al., 2017; Rao et al., 2021; Cai et al., 2022a).

T. scabra is not an economically important species, but it is also a common dominant species in the macrozoobenthic community along the coast of Guangdong (Du et al., 2011; Li et al., 2016). T. lata is neither an economically important nor a dominant species on the southeast coast of China (Jiang et al., 1990a; Jiang et al., 1990b; Zhang et al., 2007; Du et al., 2008a; Cai et al., 2022a). Why do heavy metals have less impact on the Daya Bay bivalves? That is, there is no significant correlation between bivalves and heavy metals. This is because the sediment quality radically improved after the late-2000s, and heavy metals in nearshore sediments of Daya Bay were all closely related to the import of anthropogenic and/or terrestrial material, whereas those offshore were likely to be related to the joint influence of anthropogenic activities and natural processes (Du et al., 2008; Qu et al., 2018).

The GMRI of P. undulatus in all four seasons was less than 50%; therefore, the risk of blocking the nuclear power cold source system was low in Daya Bay (Cai et al., 2022b).

Amphioplus (Lymanella) laevis has been a dominant species in Daya Bay for over 30 years (Jiang et al., 1990a; Jiang et al., 1990b; Zhang et al., 2007; Du et al., 2008a; Du et al., 2008b; Du et al., 2009; Huang et al., 2012; Zhang et al., 2017; Cai et al., 2022b). It is also a common dominant species on the coast and in the Gulf of China (Cai and Zheng, 1994; Chen et al., 2020). The genus Amphioplus is also dominant in some open-sea areas. Amphioplus sinicus was not only the dominant species of macrozoobenthic communities in Hailing Bay, western Guangdong, but also outside the bay, in addition to being the annual dominant species (Li et al., 2018). A. laevis is the dominant species in subtidal amphioxus habitats in Xiamen Bay (Chen et al., 2020; Rao et al., 2020b). BIO-ENV analysis revealed significant seasonal variations in environmental factors affecting community structure (including A. laevis) in semi-enclosed waters in Bohai Bay, China (Shi et al., 2022). This study found that the density of A. laevis showed a significantly positive correlation with water depth, salinity, and Pb and a significantly negative correlation with silicate and TP.

P. bidentata has not been the dominant species in previous surveys in Daya Bay (Jiang et al., 1990a; Jiang et al., 1990b; Du et al., 2008a; Du et al., 2008b; Du et al., 2009; Zhang et al., 2017; Yuan et al., 2017; Rao et al., 2020a), but it was a common species in the sea area near the nuclear power plant (Cheng et al., 2018; Li et al., 2020; Zhang et al., 2022b), It was a dominant species on the northern coast of Zhejiang Province (Yan et al., 2020) and in the surrounding waters of Qinshan Island (Mao et al., 2022). This study found that the density and biomass of P. bidentata were significantly positively correlated with water depth and DO.

The GMRI of P. bidentata in all four seasons was less than 50%, so the risk of blocking the nuclear power cold source system was low in Daya Bay (Cai et al., 2022b).

The main sources of pollution near the Daya Bay nuclear power plant were the warm drainage from the power plant and the shallow aquaculture area of Dapeng Ao, but the warm drainage from the power plant had less effect on benthic macrofauna in nearby waters (Zhang et al., 2007). From 1982 to 2004, the ecological environment of Daya Bay changed from 237 species in 1987 to 194 species in 1997 (Wang et al., 2004), and the mean biomass and species of benthic animals near power plants ranged from 317.9 g/m² in 1991 to 45.24 g/m² in 2004, and from 250 species in 1991 to 177 species in 2004 (Wang et al., 2008). The correlation analysis between the macrobenthos community and environmental factors indicated that secondary productivity was significantly affected by the content of inorganic nitrogen, phosphorus, dissolved oxygen of seawater, and organic carbon content of sediment (Liu et al., 2018). Changes in the ecological environment in Daya Bay are related to sewage discharge and mariculture. Small body size, short longevity, and high tolerance species are more abundant in industrial sewage discharge areas and mariculture areas (Rao et al., 2021). We believe that to accurately identify the main environmental factors affected, long-term monitoring and sufficient comparable data are needed.