Anthropogenic pollutant-driven geographical distribution of mesozooplankton communities in estuarine areas of the Bohai Sea, China

Mesozooplankton communities in marine ecosystems are mainly influenced by both anthropogenic pollutants (e.g. nutrients and heavy metals) and natural variables (e.g. temperature, salinity and geographic distance). To achieve a deeper understanding of the effects of anthropogenic pollutants on mesozooplankton communities, we analyzed the community structure of mesozooplankton from 91 stations representing five typical estuarine regions in the Bohai Sea and assessed the relative importance of anthropogenic pollutants and natural variables by using multiple statistical approaches. Cd was identified as the leading pollutant for observed community variation among the five regions, followed by NH4-N and COD. Redundancy analysis (RDA) model demonstrated that mesozooplankton communities were largely determined by both anthropogenic pollutants and natural variables, and the indicator species of mesozooplankton also varied when responding to different factors. Variance partitioning analysis showed both anthropogenic pollutants and natural variables posed significant influences (ANOVA, P < 0.05) on the mesozooplankton community structure, but the explanatory power of anthropogenic pollutants overrode the natural variables. These observations highlighted the importance of anthropogenic pollutants in the shifts of zooplankton structures among different regions. Our results obtained in this study provided new insights into the mechanism of the influence of anthropogenic pollutants on mesozooplankton communities in estuarine areas.

Zooplankton inhabiting coastal ecosystems are key components in the food webs, acting an important trophic link between the marine primary producers and higher trophic levels 3,19 . The crucial role in marine food webs indicates that decline in biodiversity of zooplankton could decrease the survival rates of higher trophic organisms such as fish 20 , and eventually pose far-reaching consequences for the marine food webs 21 . Moreover, zooplankton assemblages consist of species highly sensitive to environmental pollutants such as nutrients and heavy metals, and their variation could be attributed to different levels of pollution in marine systems, such as estuarine systems 19,22,23 and coastal systems 24,25 , also making them good indicators to explore ecological effects of anthropogenic pollutants. Available studies mainly focused on ecological effects of environmental variables and hydrological processes on zooplankton community structure 26,27 . However, natural variables and anthropogenic pollutants may act synergistically, additively or antagonistically [28][29][30][31][32][33] , suggesting that it is essential to distinguish respective effects of these two categories of factors.
The Bohai Sea, which is a semi-closed inner sea located in the northeast China, has been impacted by serious pollution caused by anthropogenic activities. The coastal regions support various industries derived from extremely rapid urbanization and industrialization in the past decades, and such rapid development along the coastal regions has caused severe pollution problems to Bohai Sea, especially for the coastal estuaries. The pollution levels varied among the coastal estuaries of Bohai Sea, especially for heavy metal pollutions 3,[34][35][36][37] . For example, Jinzhou Bay, which is surrounded by one of the old industrial bases in China, was found with the highest concentrations of Cd, the pollution levels gradually decreased along Jinzhou Bay, Luanhekou Estuary, Laizhou Bay, Shuangtaizi Estuary and Bohai Bay 36 . Thus, the environmental gradients along the coastal bays or estuaries provide premise conditions for the exploration of ecological response of zooplankton under anthropogenic stress.
In this study, we aimed to: (1) characterize the community structure of mesozooplankton in the typical Bohai estuarine regions, (2) identify the key factors responsible for local and regional community structure variation, and further (3) explore the respective ecological roles of anthropogenic pollutants and natural variables on community structure. Results environmental features. Generally, the average value of all environmental factors showed significant difference (P < 0.05) among the five regions ( Fig. 1; Table 1). Tests based on nonmetric multidimensional scaling (NMDS) and the analysis of similarity (ANOSIM) also provided similar results (Stress = 0.14, Global R = 0.409, P = 0.001; Fig. 2A), supporting significant environmental gradients among the five regions.
Specifically, the concentrations of As, Hg, Pb and Cd were the highest in Jinzhou Bay among the five regions and significantly differed from those in other regions (Table 1). While, the concentration of Zn in Laizhou Bay was significantly higher than that in the rest four regions. The concentration of Cu in Jinzhou Bay, Luanhe Estuary and Bohai Bay was significantly higher than that in Laizhou Bay and Yellow River Estuary. One organic pollutant indicator, COD, in Bohai Bay and Laizhou Bay were significantly higher than that in Yellow River Estuary, Luanhe Estuary and Jinzhou Bay. Luanhe Estuary and Bohai Bay showed significantly higher salinity than that in Jinzhou Bay, Laizhou Estuary and Yellow River Estuary. The concentration of NH 4 -N in Yellow River Estuary and Bohai Bay was significantly higher than that in the rest three regions. While, the concentration of NO 3 -N in Laizhou Bay was highest and significant difference was observed between Laizhou Bay and the other regions (Table 1).  www.nature.com/scientificreports www.nature.com/scientificreports/ Zooplankton community structure. A total of 43 mesozooplankton species were identified in the five regions across all 91 sampling stations ( Table 2). The mesozooplankton communities were dominated by Paracalanus parvus (26.94%), Acartia bifilosa (25.30%), Oithona similis (11.50%) and A. pacifica (9.21%). When the four diversity index (Species richness, Shannon-Wiener index, Pielous' evenness and Simpson diversity index) within each region were subjected to the non-parametric Mann-Whitney U test, we found significant difference between regions (P < 0.05 for most pairs, Table 2). In addition, a coincident result was also observed based on statistical percentiles, as there appeared different median values and 25 th and 75 th percentiles, such as for Shannon-Wiener index and Pielous' evenness (Fig. 3).
Explanatory variables for observed community structure. A total of 17 environmental factors and 18 spatial vectors transformed from geographic coordinate values for sampling stations were considered as possible drivers for structuring observed zooplankton community. Based on forward selection, 13 key factors showed significant influence on the observed community structure (P < 0.05, Table S2), including nine anthropogenic pollutants (Cd, NH 4 -N, COD, As, Hg, DO, NO 2 -N, Cu and NO 3 -N) and four natural variables (two spatial variables (V1and V2), salinity and temperature). The selected factors were used for the construction of parsimonious RDA model, which was globally significant (P = 0.001) with 49.52% of adjusted R 2 , and the first two axis (RDA1 and RDA2) explained 23.96% and 13.21% of total variations of mesozooplankton communities, respectively (Fig. 4). For anthropogenic pollutants, Cd was the largest contributor in affecting zooplankton community structure, followed by NH 4 -N, COD, As, Hg, DO, NO 2 -N, Cu and NO 3 -N (Fig. 4, Table S2). For example, one high abundance species O. similis showed positive correlation with Cd (Fig. 4, Table 2). For natural variables, V2 was the leading contributor of variation in zooplankton community structure, followed by salinity, V1 and temperature (Fig. 4, Table S2). Specifically, O. similis, P. crassirostris and A. longirostris were negatively correlated with V2, whereas P. parvus was positively correlated with V1 (Fig. 4). In addition, we detected different influence of anthropogenic pollutants on mesozooplankton. For example, the abundance of O. similis, P. crassirostris and Acanthomysis longirostris were mainly affected by Cd and As, whereas both heavy metals posed no obvious influence on A. bifilosa (Fig. 4). Indeed, the abundance of A. bifilosa was mainly affected by COD, whereas COD had no obvious influence on P. parvus (Fig. 4). Spearson correlation analyses also confirmed this pattern: the abundance of O. similis was positively correlated with the concentration of Cd (rho = 0.336, P = 0.01, Fig. S1) and A. bifilosa was positively correlated with COD (rho = 0.301, P = 0.04, Fig. S1).
To explore the relative roles of anthropogenic pollutants and natural variables in structuring zooplankton community, variance partitioning was performed for explanatory variables based on nine anthropogenic pollutants (Cd, NH 4 -N, COD, As, Hg, DO, NO 2 -N, Cu and NO 3 -N) and four natural variables (V1, V2, salinity and temperature). The results showed that the shared explained percentage between anthropogenic pollutants and natural variables was 25.7%. Anthropogenic pollutants alone explained 13.9% of the total variation of community structure when excluding the influence of natural variables. Conversely, natural variables alone explained 9.8% of the total variations of community structure when removing the anthropogenic pollutant influence (Fig. 5).

Discussion
The present study showed that the marine water quality varied largely among the five regions, ranging from very poor for both Jinzhou Bay and Bohai Bay, due to sewage pollution, to good in farming locations such as Luanhe Estuary. Those observations highlighted the differences in their surrounding land-use. A high level of heavy metal pollution was noted in Jinzhou Bay which was characterized with a high level of Cd, As and Hg pollutions in this investigation, mainly due to wastewater discharge by surrounding mining industries in the past century 36 . The organic pollutant concentrations in marine water columns were found to be generally high for Bohai Bay, with high COD being found in the Bohai Bay sampling stations because of sewage discharge from the big cities (Beijing and Tianjin) 38,39 . Whereas the opposite water quality states was observed in Yellow River Estuary, as shown in Fig. 4. Our observations are in agreement with recent studies which have demonstrated that the differences in land-use largely contribute to changes in pollutant variables in aquatic systems 29,36,40 . Hence, the understanding of anthropogenic pollutions induced by different land-use would provide insights into the management of aquatic resource. www.nature.com/scientificreports www.nature.com/scientificreports/ Our results clearly showed that both anthropogenic pollutants and natural variables play a crucial role in structuring variation of mesozooplankton communities in the Bohai Sea, the influence of anthropogenic pollutants on mesozooplankton communities overrode that of natural variables. Our RDA model showed that the  Table 3. Results of the similarity percentage (SIMPER) analysis for each region in Bohai Sea.   www.nature.com/scientificreports www.nature.com/scientificreports/ variations of mesozooplankton communities were in line with the changes in anthropogenic pollutants and natural variables among the five regions, with the effects of anthropogenic pollutants and natural variables being integrated into overall resultant mesozooplankton communities, which was in agreement with recent studies 28,29 . Some species may be sensitive to anthropogenic pollutants, some species may be influenced by natural variables and some other species can be affected by both factors. Those findings were consistent with previous studies in other aquatic systems such as lake, river and marine systems, they have demonstrated that the species' responses to environmental factors varied largely 19,22,23,33,[41][42][43][44][45] . Thus, different mesozooplankton species responded differently to anthropogenic pollutants and natural variables in the current study, suggesting some species of mesozooplankton communities at different sampling locations can be used as reliable ecological indicators to reflect local environmental pollutions.  47,48 . Oithona was characterized with flexible diet, shorter life cycle and higher reproduction rate compared to other mesozooplankton 48,49 , those physiology features could partly explain its higher tolerance to anthropogenic pollutions.
In addition, A. bifilosa showed positive response to COD concentration (A. bifilosa-COD, rho = 0.301, P = 0.04), indicating the variation of A. bifilosa abundance may mirror the fluctuation of COD in marine coastal regions. Indeed, COD contents can reflect the concentrations of anthropogenic organic matters, which can provide aquatic animals with foods and often promote the proliferation and growth of some species of zooplankton communities 50 . Moreover, the abundance of P. parvus was negatively correlated with NH 4 -N concentration (Fig. 4), which was also selected by forward selection as an anthropogenic pollutant with significant influences on mesozooplankton communities. High concentration of NH 4 -N can decrease the rates of growth and survival for copepod species 51 . Hence, it is expected that NH 4 -N toxicity can cause similar negative influence on P. parvus. However, other factors such as high trophic predators may affect the determinate of bioindicator species, and species may show different tolerant abilities to toxic pollutants under different environments 29 , more studies based on different trophic levels and toxicological validations under laboratory conditions should be performed to verify those findings in the field.
In the current study, Cd was identified as the leading anthropogenic pollutant structuring zooplankton communities in coastal regions of Bohai Sea. Whereas inorganic nitrogen was considered as the main driver for observed zooplankton structure in Tangshan Bay 52 . Moreover, total phosphate, NH 4 -N and Mg 2+ were identified as the leading factors in Chaobai river, Beiyun river and Fuyang river of China, respectively 41,43,44 . Our results combined with other findings revealed that different leading anthropogenic pollutants were recognized in different ecosystems, even in different regions of the same ecosystems. Those findings indicated that the leading anthropogenic pollutants for observed structure of zooplankton community may be regional-specific, suggesting more works should be performed to clarify driving mechanism of zooplankton structure in regional scale.
Zooplankton community diversity within each region is traditionally used to reflect local environmental conditions 53,54 . Our study showed that Bohai Bay and Jinzhou Bay were both heavily polluted by anthropogenic pollutants. However, two Bays showed large variations in diversity index in the current study. Bohai Bay was mainly suffered from organic pollution, which may pose toxic effects on some species but also may provide more food resources to support more species 55,56 , thus resulted in high species diversity. Jinzhou Bay was mainly polluted www.nature.com/scientificreports www.nature.com/scientificreports/ by heavy metals, which may only have negative effects on aquatic species, leading to decline of mesozooplankton diversity. Those observations suggested zooplankton community diversity may be not an ideal index to mirror local anthropogenic pollutions, whereas maybe an ideal index to reflect local organic pollutions. Alternatively, the shifts of mesozooplankton community composition in major groups between regions, such as the variation of bioindicator species, may be good indexes to represent local anthropogenic pollutions, as have been revealed at other communities such as benthic communities 29,57,58 .
Traditionally, natural variables including hydrological processes, temperature and salinity are often considered as the major determinants of community structures, and other variables including pollutions played less importance [30][31][32] . For example, the identification of V2 as leading natural factor indicated that spatial variables such as ocean current may play a crucial role in driving spread of some species, including O. similis, P. crassirostris and A. longirostris, and eventually affect local community structures. However, variance partitioning in our study showed that the explained percentage of community variation by pollutants alone was larger (13.9%) than that by natural variables (9.8%) alone, indicating that anthropogenic pollutants contributed larger than natural variables to the variation of zooplankton community in the coastal regions of Bohai Sea. In the open seas, natural variables such as ocean current and water temperature may play major roles in structuring zooplankton structure 59,60 . Whereas in semi-closed seas such as Bohai Sea, both higher pollution levels 36,46 and weaker ocean currents 61 may lead to the observed distribution patterns of mesozooplankton communities in the present study. Thus, the relative roles of anthropogenic pollutants and natural variables in shaping zooplankton structure may largely depend on the relative strengths of both factors. In addition, about quarter of total variations of zooplankton community structure was simultaneously explained by anthropogenic pollutants and natural variables, suggesting strong colinearity between two types of factors exist in the coastal regions of Bohai Sea. Strikingly, P. parvus showed positive correlation with V1, dispersal alone cannot explain this observation. One reasonable reason may be that V1 such as ocean current-dispersal largely shaped unmeasured variables or biological processes, which may be actual factor affecting P. parvus. Those observations indicated natural variables such as salinity and hydrological processes in this study may affect degradation and dispersal of pollutants 6,46 , highlighting the necessity of excluding the influence of natural variables when exploring the ecological effects of anthropogenic pollutants on zooplankton composition in the field. This idea is especially applicable for marine systems at coastal waters because of complexities of hydrological processes and obvious gradients of natural variables.
In conclusion, our study clearly showed that the mesozooplankton communities among the five regions varied significantly along the environmental gradients. Multiple analyses identified that both anthropogenic pollutants and natural variables were major factors driving mesozooplankton communities in the coastal marine system. Cd was identified the leading anthropogenic pollutants factor structuring mesozooplankton community, followed by Hg, COD, NH 4 -N, As, Zn, NO 2 -N. The species responses to those environmental factors varied largely and mainly depended on organism taxa, suggesting some species can be used as potential bioindicators of environmental pollutants. Further analyses showed that anthropogenic pollutants still played a major role with significant influence on the mesozooplankton community even after removing the natural variable influence, highlighting the necessity of considering negative effects of anthropogenic pollutants on coastal ecosystems in environmental management and monitoring programs. Methodologically, our results emphasized the importance of excluding influence of natural variables including hydrological processes, temperature and salinity when exploring the ecological effects of anthropogenic pollutants on plankton community structure, especially at coastal waters. However, this study was only performed on mesozooplankton that were adequately identified based on morphological features, other zooplankton such as microzooplankton have not been test on this issues. More works on different trophic levels should be carried out using feasible molecular-based methods such as metabarcoding-based identification approach 62 . Material and Methods study region and sampling stations. This study focused on five important estuarine areas of the Bohai Sea. The sampling stations mainly distribute in the shallow coastal areas, the water depth of sampling stations is between 2.5 m and 17.0 m. The water column is generally mixed homogeneously due to strong tidal mixing. Neither thermocline nor halocline was observed in the sampling stations of this study because the summertime stratification of the water column mainly occurs in the deep basins (25~35 m depth) in the central Bohai Sea 63 . A total of 91 sample stations were set up over the coastal area of Bohai Sea (Fig. 1), including Jinzhou Bay (12 stations), Luanhe Estuary (24 stations), Bohai Bay (20 stations), Yellow River Estuary (16 stations) and Laizhou Bay (19 stations). Jinzhou Bay, located in northwest of Liaodong Gulf, is a semi-closed shallow water area. Six rivers including Lianshan River, Wuli River, Lao River, Cishan River, Zhouliu River and Tashan River flow into Jinzhou Bay. It is famous as an old industrial base, and become one of the most polluted coastal area in China. Luanhekou Estuary is located on the northwest coast of Bohai Sea with water depths less than 20 m. Freshwater and sediment discharges have decreased greatly since the 1980s due to large dams and reservoirs built along the Luanhe River. Bohai Bay is located on the west of Bohai Sea, near the city of Tianjin and Beijing. Bohai Bay is a typical semi-enclosed coastal area and has limited water exchange with the ocean. Large quantities of industrial and domestic wastewater discharges flow into Bohai Bay from rivers of Beijing-Tianjin. The western coast of Bohai Bay locates the Tianjin Ports, the 10th largest port in the cargo throughout in the world. Yellow River Estuary is located in the southwestern part of Bohai Sea, the end of the second largest river (Yellow River) of the world in terms of sediment load. Yellow River Estuary is characterized with high concentration of Ammonia nitrogen 64 . Laizhou Bay is located on the southern part of Bohai Sea, accounting for up to 10% of the total area. It's a semi-closed shallow area with average water depth less than 10 m. There are more than a dozen of rivers running into the Laizhou Bay, among which Yellow River, Xiaoqinghe River and Wei River are the most important. All samples were collected in the August, 2015 (Table 1).
www.nature.com/scientificreports www.nature.com/scientificreports/ Zooplankton sampling and enumeration. Mesozooplankton samples were quantitatively collected in each sampling station. Specifically, we firstly measured water depth for each sampling station and collected mesozooplankton samples using a plankton net (505 μm mesh size, 50 cm mouth diameter) by towing vertically from 2 m above the bottom to the surface with a speed of 0.5-0.8 m/s. The filtered water volume (m 3 ) was measured using the rope length multiplied by the mouth area (0.2 m 2 ). The samples were collected and preserved immediately in 5% formaldehyde. In the laboratory, all individuals (zooplankton larvae were not included) were identified into species and enumerated. The abundance (ind./m 3 ) of each species was calculated as the number of individuals divided by the filtered water volume. In cases when the mean is presented, the standard deviation was provided (mean ± SD). environmental variable sampling and analysis. Surface seawater samples were collected with a 5 L Niskin bottles from 0.5 m below the surface at each station. The seawater salinity and temperature was measured in situ with a multiparameter sensor YSI6600, and pH values were determined with a pH meter. The seawater for dissolved oxygen (DO) analysis was collected with a tube reaching the bottom of bottle until the water overflowed. Suspended matter samples were filtered through pre-weighted Whatman GF/F fiber filters (25 mm). The samples for metal determination were filtered immediately through Whatman GF/F fiber filter (0.45 mm), and then acidified with 10% HNO 3 , placed in an ice box and transported to the laboratory. Concentrations of NO 3 -N, NO 2 -N, NH 4 -N and PO 4 -P in seawater were determined according to the methods described by Grasshoff et al. 65 . DO was determined using the Winkler titration method. Chlorophyll-a (Chl-a) was determined by filtering 100-200 mL seawater onto GF/F fiber filter by a cascading filtering device under low vacuum pressure. After extraction with 90% acetone, Chl-a was determined by a Turner Design fluorometer (TD Trilogy). The concentrations of dissolved heavy metals were determined using the inductively coupled plasma mass spectrometry (ICP-MS, Thermo X series) for Cd, Pb, Zn and Cu, while the content of Hg and As was determined using the atomic fluorescence spectrometer (AFS-920). spatial variables. Besides environmental variables, the variations on the spatial distribution of aquatic communities are traditionally correlated with geographical distances between sampling stations 43,44,66 . The spatial distances were generated based on Cartesian coordinates and Euclidian distance matrix, which were transformed from longitude and latitude among the sampling stations. In detail, the longitude and latitude were converted to Cartesian coordinates using the geoXY function available in the SoDa packages in R software v.3.4.1 67 . Then, an Euclidian distance matrix on this Cartesian coordinates was computed using the dist function and PCNM (Principal Coordinates of Neighbor Matrices) analysis (permutations = 1000) was performed on this matrix using PCNM function implemented in the PCNM package. The method of PCNM 68,69 can effectively model spatial structure in biological communities among sampling stations 70 and has been increasingly used in various groups including bacteria and phytoplankton 71,72 . In this study, we attempted to apply the method of PCNM to mesozooplankton in order to understand the effects of spatial variables on mesozooplankton community. The number of PCNM variables formed is always dependent on the number of sampling stations and their spatial relations. At last, a total of 18 PCNM vectors (V1-V18) showing positive spatial autocorrelation were formed and used as spatial variables for subsequent redundancy analysis (RDA) and forward selection. In detail, the first PCNM vectors indicate spatial relations among sampling stations at a large scale (e.g. between sampling stations across regions) and the last PCNM vectors represent spatial relations among a small scale (e.g. between sampling stations in the same region). statistics analysis. In order to separate the effects of anthropogenic pollutants on mesozooplankton communities from natural environmental factors, the 17 environmental variables together with 18 spatial variables were reclassified into two groups: natural variables (temperature, salinity, and spatial variables) and anthropogenic pollutants (COD, suspend matters, DO, Chl-a, pH, PO 4 , NO 2 -N, NO 3 -N, NH 4 -N, As, Hg, Cu, Pb, Cd, and Zn). The average value and standard deviation for each environmental variable and study location were calculated. One-way ANOVAs were used to compare means of environmental variables among study locations, after testing for homogeneity of variances (Levene's test, P < 0.05) and normality of distribution (Shapiro-Wilk test, P < 0.05) using Paleontological Statistics (PAST) version 3.01 73 . Significant ANOVAs (P < 0.05) were followed by Tukey HSD post hoc analysis to identify differences between study locations using PAST version 3.01.
Before statistical analyses, all measured environmental factors (except for pH) and mesozooplankton data were log10 (x + 1) transformed to improve normality. To characterize distribution patterns of zooplankton, the composition and abundance of zooplankton were analyzed using non-parametric multivariate methods implemented in PRIMER 5.0 74 . The abundance of zooplankton between regions was compared using nonmetric multidimensional scaling (NMDS) and the analysis of similarity (ANOSIM), which is based on Bray-Curtis distance and rank dissimilarity. The major species driving distribution patterns of zooplankton assemblages at both inter regions and intra regions were identified using similarity percentage analysis (SIMPER) with a cutoff of 90% contributions. The NMDS, ANOSIM and SIMPER analyses were performed using PRIMER 5.0 74 To recognize the major factors responsible for observed zooplankton community structure, we performed the linear ordination method of RDA, which was chosen mainly based on a preliminary detrended correspondence analysis (DCA) on zooplankton community. The DCA showed that the longest length of gradient (3.03) was shorter than four, indicating that the majority of taxa showed a linear response to explanatory variables 75 . To avoid multicollinearity problems and construct parsimonious RDA model, which has been proved to have greater predictive power for the relationship between zooplankton communities and explanatory variables 76 , we conducted forward selection to select significant explanatory variables including environmental factors and spatial variables using the forward.sel function (ANOVAS; 1000 permutations) in packfor package in R, which simultaneously taken account for significance (P < 0.05) and adjusted R 2 of the global RDA model with all available (2019) 9:9668 | https://doi.org/10.1038/s41598-019-46047-5 www.nature.com/scientificreports www.nature.com/scientificreports/ explanatory variables 77 . To verify the correlations obtained from RDA analysis, additional Spearman correlation analysis was also performed.
To evaluate the ecological effects of anthropogenic pollutons on mesozooplankton structures, variance partitioning and partial redundancy analysis (pRDA) were performed to estimate explained percentage of the significant anthropogenic pollutions and natural variables selected by forward selection. Variance analyses (ANOVAS; 1000 permutations) were performed to test the significance of RDA and pRDA. Those analyses including RDA, pRDA, ANOVA and DCA analyses were computed using vegan package in R software.