The current survey was carried out from 19 to 23 April 2019, and 30 surface sediments were collected from Houshui Bay and the outer part of Yangpu Bay ( Figure 1 ). The geographic locations of the sampling sites were recorded by using a global positioning system (GPS). Surface sediments were collected with a grab sampler, the grab was carefully opened in a container, and sediments were laid down in their initial position to avoid sediment disturbance (Debenay et al., 2001). The upper 2 cm layer of undisturbed sediment was sampled for laboratory analyses. Two sets of samples were prepared on the boat: one set was for foraminiferal analysis, and the other was for chemical analyses. For the preparation of samples for foraminiferal analysis, approximately 100 cm³ of sediment was placed into plastic bottles containing 4% neutrally buffered formaldehyde solution (Tonghe Biotechnology, Zhejiang, China) and rose Bengal stain (2 g/l). All samples were kept on ice until transport to the laboratory for further processing.

Samples were dried, ground, and reduced to fine powder for the measurement of heavy metals. To detect the concentrations of Cu, Pb, Zn, Cd, Cr, Mn, Mo, and Co, approximately 0.1 g samples were digested by 2 ml of HNO₃ (65%, v/v) and 1 ml of HF (40%, v/v) at 190°C for 48 hours. Then, samples were digested by 1 ml of HNO₃ (65%, v/v) and 0.5 ml of HClO₄ (70%, v/v) and diluted to 100 ml using Milli-Q water. The concentrations were measured with an inductively coupled plasma mass spectrometer (ICP-MS, Agilent 7900, Agilent Technologies Japan, Ltd., Japan). To detect concentrations of Sb, Hg, and As, samples were digested by aqua regia at a controlled temperature (step 1: 120°C for 2 min; step 2: 150°C for 5 min; step 3: 185°C for 40 min). Then, the concentrations were measured by using an atomic fluorescence spectrophotometer (AFS-930, Beijing Jitian Instruments, Beijing, China). Duplicates, blanks, and standard reference materials (GBW07333 for marine sediment) were used for quality assurance and quality control. The measured concentrations of the reference materials were in satisfactory agreement with the certified values. The detection limits for the heavy metals were as follows: Cu (0.6 mg kg⁻¹), Pb (2 mg kg⁻¹), Zn (1 mg kg⁻¹), Cd (0.09 mg kg⁻¹), Cr (2 mg kg⁻¹), Mn (0.4 mg kg⁻¹), Mo (1 mg kg⁻¹), Co: 0.04 mg kg⁻¹, Sb: 0.01 mg kg⁻¹, Hg: 0.002 mg kg⁻¹, and As: 0.01 mg kg⁻¹. Metal concentrations below the detection limit were replaced with half of the detection limit.

For the measurement of total organic carbon (TOC), samples were freeze-dried, homogenized and acidified with diluted HCl (1%) until all carbonates were removed. TOC and total nitrogen (TN) contents were determined through an Elementar Vario Micro Cube elemental analyzer (Elementar Analysensysteme GmbH, Hanau, Germany).

The 16 USEPA priority PAHs were analyzed, including naphthalene (Nap), acenaphthylene (Acy), acenaphthene (Ac), anthracene (Ant), phenanthrene (Phe), fluorene (Flu), fluoranthene (Fla), pyrene (Pyr), benz[a]anthracene (BaA), chrysene (Chr), benzo(a)-pyrene (BaP), benz[b]fluoranthene (BbF), benz[k]fluoranthene (BkF), dibenz[a,h]anthracene (DahA), indeno[1,2,3-cd]pyrene (InP), and benzo(ghi)-perylene (BgP). Approximately 10 g of sample was weighed and dried over anhydrous sodium sulfate and put into a glass tube with 50 μl of decafluorobiphenyl solution. The glass tube was then placed in a Soxhlet extractor, and 100 ml of acetone/n-hexane mixture (1/1 v/v) was added for extraction. The extract was filtered through a glass fiber membrane and cleaned up by an acetone/n-hexane mixture. After adding 5 ml of n-hexane, the eluent was collected and concentrated to 1 ml and further purified by silica gel column chromatography using anhydrous sodium sulfate and a small amount of dichloromethane. The 16 PAHs were detected using a Shimadzu LC-2010A HPLC system (Shimadzu Co., Kyoto, Japan). PAHs were identified by comparing the retention time of the chromatographic peak between the sample and standard solutions. The contents were calculated from the calibration curve and then converted to ng g⁻¹ in the dry weight of the sample. The detection limits of PAHs were 0.3–0.5 ng g⁻¹. PAH concentrations below the detection limit were set at one-half the detection limit.

Samples were washed through a 63-μm sieve to remove mud and debris and then dried at 40°C in an oven. Foraminifera (rose Bengal-stained) were identified under a Nikon SMZ1000 stereomicroscope. The minimum number of individuals was 150 per sample, while for those samples with fewer than 150 individuals, all individuals were counted. Species classification was based on some previous work in the China Seas (Zheng, 1979; Wang et al., 1988; Lei and Li, 2016) and the taxonomy of the World Modern Foraminifera Database (Hayward et al., 2020).

Diversity indices, e.g., foraminiferal density (individuals per 10 cm³), species richness (number of species), Shannon–Weaver index H(S), and Fisher α index, were calculated using PRIMER v5 software (Plymouth Routines In Multivariate Ecological Research; Clarke and Gorley (2001)).

A two-way (simultaneous R- and Q-modes) hierarchical cluster analysis based on the relative abundance of each species was performed by using PCORD version 5.0 (mJm Software). The Sørensen distance (or similarity) coefficients were used for distance measurements, and clusters were defined according to Ward’s merging criterion. Canonical correspondence analysis (CCA) was performed using CANOCO (v. 5.0) software to determine the variation in species composition explained by environmental variables (Šmilauer and Lepš, 2014). To eliminate the effects of stochastic variations, species abundance was ln(abundance + 1)-transformed prior to the CCA analysis (Morey et al., 2005), and environmental variables (heavy metals, TOC, and PAHs) were standardized to a mean of zero and standard deviation of unity.

A total of 11 heavy metals (Cu, Pb, Zn, Cd, Cr, Mn, Mo, Co, Sb, Hg, and As) were measured in 30 sediments from Houshui (17 sites) and Yangpu (13 sites) Bays. Table 1 highlights the ranges and average concentrations of these elements. In Houshui Bay, the average concentrations were Cu, 16.9 mg kg⁻¹; Pb, 28.1 mg kg⁻¹; Zn, 85.3 mg kg⁻¹; Cd, 0.082 mg kg⁻¹; Cr, 63.6 mg kg⁻¹; Mn, 601 mg kg⁻¹; Mo, 32.1 mg kg⁻¹; Co, 13.5 mg kg⁻¹; Sb, 0.32 mg kg⁻¹; Hg, 0.097 mg kg⁻¹; and As, 20 mg kg⁻¹. In Yangpu Bay, the average concentrations were Cu, 7 mg kg⁻¹; Pb, 23.9 mg kg⁻¹; Zn, 59 mg kg⁻¹; Cd, 0.04 mg kg⁻¹; Cr, 39.1 mg kg⁻¹; Mn, 494 mg kg⁻¹; Mo, 22.4 mg kg⁻¹; Co, 10 mg kg⁻¹; Sb, 0.04 mg kg⁻¹; Hg, 0.044 mg kg⁻¹; and As, 15.2 mg kg⁻¹. Obviously, heavy metal concentrations in Houshui Bay were higher than those in Yangpu Bay; in particular, the average concentrations of Cu, Cd, Sb, and Hg in Houshui Bay were more than twice those in Yangpu Bay. Compared with previous investigations of Cu, Pb, Zn, Cd, Cr, Hg, and As in the eastern and western continental shelves of Hainan Island (Hu et al., 2013; Xu et al., 2015), the average concentrations of Cu, Pb, Zn, Cd, Cr, and Hg were at the same levels in Houshui Bay, while they were lower in Yangpu Bay. The As concentration in the two bays was remarkably higher than its levels in the shelves.

The spatial patterns of heavy metals are shown in Figure 2 , and all metals exhibited a similar distribution pattern. In Houshui Bay, higher concentrations were observed in the northwest and southeast areas, and relatively lower concentrations were observed in the central area; in Yangpu Bay, heavy metal concentrations gradually increased from northeast to southwest.

In total, the 30 sediment samples were measured to determine the content of TOC. The TOC content was higher in Houshui Bay (ranging from 0.06% to 0.73%, mean 0.45) and relatively lower in Yangpu Bay (ranging from 0.05% to 0.55%, mean 0.32). The spatial pattern of the TOC content was similar to that of heavy metals: in Houshui Bay, the TOC content reached its highest level in the northwest and southeast areas, and the lowest content was found in the central area; in Yanpu Bay, the TOC content roughly increased from northeast to southwest ( Figure 3 ).

Sixteen PAHs were detected in 19 sediments from the bays of Houshui (10 sites) and Yangpu (9 sites). The concentrations of 18 PAHs (the measured 16 PAHs and two other PAHs) in sediments of the two bays and river and coastal sediments of Hainan Island are summarized in Table 2 . A comparison of the average concentrations of TPAHs in different regions of Hainan Island showed a descending order of Houshui Bay, Yangpu Bay, rivers, and offshore areas. In particular, the concentrations of Ace, BaA, BbF, BkF, and BgP were remarkably higher in the two bays than in rivers and offshore areas.

Spatial distributions of low-molecular weight PAHs (LPAHs, 2- to 3-ring PAHs), high-molecular weight PAHs (HPAHs, 4- to 6-ring PAHs), and TPAHs are presented in Figure 3 , showing a very similar pattern among them. In Houshui Bay, the concentrations of LPAHs, HPAHs, and TPAHs tended to increase from northwest to southeast; in Yangpu Bay, their concentrations were higher in the central part than in the northeast and southwest areas.

Quantitative analyses of benthic foraminifera revealed 127 species belonging to 59 genera and 38 families. Among these species, 112 were found in Houshui Bay (17 sites), and 87 were found in Yangpu Bay (13 sites). Thirty-five species exhibited a relative abundance higher than 3% at least three sites and were thus considered commonly occurring species (see Table 3 ). The foraminiferal communities in both bays were dominated by hyaline-perforate and porcelaneous taxa, which accounted for 47.7% and 50.1% of the average percentages in Yangpu Bay and 48.5% and 37.5% in Houshui Bay, respectively ( Figure 4 ). Table 4 lists the 14 most abundant species (average percentage > 2%) in these two taxa: 7 porcelaneous species (Adelosina colomi, Adelosina longirostra, Massilina laevigata, Quinqueloculina laevigata, Quinqueloculina seminula, Quinqueloculina spp., and Spiroloculina communis) and 7 hyaline species (Ammonia beccarii, Amphistegina exquisite, Amphistegina sp., Elphidium advenum, Elphidium crispum, Elphidium spp., Nonion commune, and Rotalidium annectens). Among these, the dominant species of porcelaneous types were Massilina laevigata, Adelosina laevigata, Quinqueloculina spp., Spiroloculina communis, which occured mainly closer to the coast in southeast of Yangpu Bay and central to outer region of Houshui Bay ( Figure 5 ). The dominant species of hyaline-perforate types were Elphidium crispum, Amphistegina exquisita, Amphistegina sp., Rotalidium annectens, which occured mainly occur in outer of Yangpu Bay. But in Houshui Bay, the dominant hyaline-perforate species were different from those in Yangpu Bay, mainly Elphidium advenum, Rotalidium annectens, and Ammonia beccarii. Elphidium advenum were the most widely distributed and dominant in almost the whole sites, Rotalidium annectens was significantly dominated in the outer bay area, while Ammonia beccarii was also widely distributed but less abundant than the first two in terms of individual stations ( Figure 6 ).

The agglutinated foraminifera were much less diverse and less abundant in the two bays. In Houshui Bay, the agglutinated taxa contribute 2.1% to 28% to the total assemblage at each site, accounting for 14% of the average percentage in relative abundance. In Yangpu Bay, the agglutinated taxa were extremely poor and only 0.1% to 3.5% at each site, accounting for 2.2% of the average percentage in relative abundance. Only 5 species listed in Table 4 were more than 1% of average percentage, including Bigenerina nodosaria, Ammobaculites agglutinans, Textularia spp., Textularia foliacea, Haplophragmoides canariensis, which were significantly dominated in central and outer region of Houshui Bay than in Yangpu Bay ( Figure 7 ).

The average percentages of the dominant species were different in the two bays, and the dominance of the dominant species was stronger in Yangpu Bay than in Houshui Bay ( Figures 5 – 7 ).

The diversity indices showed that species richness, Shannon–Weaver index (H(S)) and Fisher α index were similar in the two bays (3.2 and 14 in Houshui Bay, 3.1 and 12 in Yangpu Bay, respectively) whereas foraminiferal density was much lower in Yangpu Bay (210 per 10cm²) than in Houshui Bay (48 per 10cm²). The spatial patterns of these indices are shown in Figure 8 , and all shared a similar pattern. In Houshui Bay, higher indices were found in the central area, and relatively lower indices occurred in the northwest and southeast areas; in Yangpu Bay, these indices reached their maximum near the coast and decreased with distance away.

All sites and the 35 common species were clustered into two groups by using two-way hierarchical cluster analysis ( Figure 9 ). Sites were roughly classified according to their locations (the bays), which revealed that foraminiferal assemblages within a bay were more similar to each other than to those from different bays. Two species groups (I and II) were divided, each of which was dominated by certain abundant species. For example, R. annectens, Q. laevigata, Q. seminula, A. beccarii, and E. advenum were assigned to group I, and M. laevigata, S. communis, E. crispum, and Elphidium spp. belonged to group II. The agglutinated species Bigenerina nodosaria, Ammobaculites agglutinans, Textularia spp., Textularia foliacea, and Haplophragmoides canariensis were clustered into group I, and most porcelaneous taxa (e.g., Massilina laevigata, Spiroloculina communis) and some hyaline taxa (E. crispum, and Elphidium spp.) were found in group II, while none of agglutinated forms were found in group II. A heat map of the relative abundance and distribution of individual species demonstrated that group I tended to be more abundant in Houshui Bay and that group II was richer in Yangpu Bay ( Figure 9B ).

CCA was applied to reveal the relationship of sites/species to heavy metals (Cu, Pb, Zn, Cd, Cr, Mn, Mo, Co, Sb, Hg, As), TOC, and PAHs. Regarding the different generation processes of HPAHs and LPAHs (Liu et al., 2000), we evaluated their effects separately. To eliminate the adverse influence of rare species (Marchant, 2002), only the 35 common species were used for the CCA analysis. Four axes were identified, and the first two together explained 41% of the variation in foraminiferal assemblages ( Table 5 ). The first axis explained a larger fraction (27.7%) of the total variation, whereas the second axis explained a smaller fraction (13.3%). The variation in the assemblages can be explained by the following environmental variables: Cr (13.8%) > Pb (10%) > Mn (9.9%) > Cd (8.4%) > Mo (8.9%) > Cu (4.8%) > Hg (4.3%) > TOC (3.5%) > Co (4.2%) > Zn (4.1%) > As (3.9%) > HPAHs (3.4%) > LPAHs (3.3%) > Sb (2.2%). Overall, heavy metals had the largest effects on foraminiferal assemblages, followed by TOC, while the effects of PAHs were lower.

Figure 10 shows the positioning of sites and species projected on the CCA biplots, where the arrows represent gradient directions of environmental variables. Along the first axis of the site-environment plot ( Figure 10A ), the sites of Houshui Bay were separated from those of Yangpu Bay. The former appeared on the left side, and the latter appeared on the right side. In the species-environment plot ( Figure 10B ), the 35 common species were divided into two groups, which was consistent with the cluster of species ( Figure 9 ), and the distribution of the two groups were also similarly in two bays ( Figure 10C ). Members of group I were closely located at the upper left, while members of group II were scattered on the right half of the plot.

The relationship of site/species to explanatory variables can be inferred from the projection of the sites/species points onto arrows of the environmental variables (Šmilauer and Lepš, 2014). As shown in Figure 10A , most sites from Houshui Bay were positively related to Cu, Zn, Cd, Cr, Mo, Co, Sb, and Hg, and most sites from Yangpu Bay were positively associated with HPAHs, LPAHs, and TOC. In Figure 10B , a positive relationship existed between species group I and heavy metals such as Cu, Zn, Cd, Cr, Mo, Co, Sb, and Hg and between species group II and HPAHs, LPAHs, and TOC. These two groups had an insignificant correlation with Mn, Pb, and As.

A previous study of a sediment core from the northern South China Sea (near Hainan Island) indicated that concentrations of heavy metals (e.g., As, Pb, Cu, and Mo) are significantly increased by 2–3 times after 1800 cal yr B.P., and human activities contribute largely to the accumulation of heavy metals (Wan et al., 2015). To date, marine culture, port development, shipping, and tourism have already become pillar industries in Hainan (The Annals of Chinese Gulfs Compiling Committee, 1999). The release of heavy metals into the island’s coastal areas can be attributed to industrial activities, agricultural processes, domestic wastes, and vehicle emissions (Hu et al., 2013). A previous survey revealed that industries such as petrochemical products, paper manufacturing, metal smelting, and oil refining surrounding Yangpu Bay discharge industrial wastes, such as Pb, Cr, and other toxic substances, into the nearshore area (He et al., 2014). Natural processes also play an important role in the enrichment of heavy metals in the coastal areas of Hainan Island (Xu et al., 2015). Rock fragments and mineral debris, which originate from weathering of rocks in the central mountain range of the island, are delivered to river networks and ultimately transported to coastal environments (bays, estuaries, and coastal zones).

The source of sedimentary organic matter can be identified from C/N values. Fresh organic matter produced by microalgae has C/N values that are usually between 4 and 10, whereas organic matter from vascular land plants commonly has C/N ratios of 20 and greater (Meyers, 1994). The C/N values in the two bays ranged from 5.2 to 8.1, indicating that organic matter was derived primarily from algae.

PAHs in marine environments can originate from accidental oil spills, offshore production, discharge from routine tanker operations, and municipal and urban runoff (Dachs et al., 1997; Zakaria et al., 2002). Potential sources of PAHs can be inferred from diagnostic ratios of PAH isomers, for example, the proportions of “kinetic” PAHs (less stable) versus their “thermodynamic” isomers (more stable) (Yunker et al., 2002). Three diagnostic ratios were used in this study, including BaA/(BaA+Chr), InP/(InP+BgP), and Fl/(Fl+Py). These ratios varied between 0.26 and 1, between 0.1 and 0.62, and between 0.1 and 0.5, respectively. According to Viñas et al. (2010), a BaA/(BaA+Chr) ratio < 0.2 indicates petroleum, 0.2–0.35 indicates petroleum and combustion (mixed sources), and > 0.35 indicates combustion; InP/(Inp+BgP) < 0.2 petroleum, 0.2–0.5 indicates petroleum combustion, and > 0.5 indicates combustion of coal, grasses and wood; and Fl/(Fl+Py) < 0.4 indicates petroleum, 0.4–0.5 indicates petroleum combustion, and > 0.5 indicates combustion of coal, grasses and wood. As shown in the scatter diagrams ( Figures 11A, B ), PAHs were generated primarily by petroleum and combustion of petroleum products. This result is consistent with the previous investigation of PAHs in seawater and sediment from Hainan Island, which was performed by Xiang et al. (2018). Considering that fishing, shipping, and tourism have developed rapidly in recent years in Hainan, we infer that commercial shipping and tourism lead to enhanced concentrations of PAHs in the sediments.

The present study used CCA to assess the effects of heavy metals, TOC, and PAHs on species distribution and habitats. Sites located in Yangpu and Houshui Bays were separated by the first CCA axis ( Figure 10A ), showing spatial differences in species composition (also indicated by the heat map of species abundance). Two species groups, divided by cluster analysis, were also identified by the ordination of species on CCA ( Figure 10B ). With respect to the species-environment relationship, group I was positively related to heavy metals (Cu, Zn, Cd, Cr, Mo, Co, Sb, and Hg), while group II had a positive correlation with HPAHs, LPAHs, and TOC. These relationships provide a means to explain species’ tolerance to different pollutants.

The tolerance and sensitivity of foraminiferal species to heavy metals have been widely discussed. It is known that Ammonia beccarii and A. tepida, the two most frequently observed nearshore species, are tolerant to heavy metals and widely accepted as indicators for heavy metals worldwide (Armynot du Châtelet et al., 2004; Le Cadre and Debenay, 2006; Frontalini and Coccioni, 2008; Debenay et al., 2009; Donnici et al., 2012; Barras et al., 2014; Li et al., 2015). However, Jorissen et al. (2018) considered A. beccarii to be sensitive to pollution. In this study, A. beccarii was strongly and positively related to heavy metals, while the relationship between A. tepida and heavy metals was relatively weak. The genus Elphidium is also abundant and widespread in marginal marine environments (Sen Gupta, 1999). Ammonia−Elphidium assemblages sometimes dominate in polluted environments [e.g., Vilela et al. (2011)]. But the effect of heavy metals on the Elphidium group is less clearly defined and may be positive, null, or negative (Bergin et al., 2006; Romano et al., 2009; Aloulou et al., 2012). Here, only two members of Elphidium, Elphidium advenum and E. hispidulum, had a positive correlation with heavy metals, whereas other members exhibited negative correlations.

Benthic agglutinated Textulariida appeared to respond positively to heavy metals, as shown in the species-environment biplot ( Figure 10B ). Unfortunately, the interaction between agglutinated foraminifera and heavy metals is rarely reported in the literature. It is known that agglutinated taxa tolerate higher freshwater flux, higher sediment and organic carbon flux, and lower oxygen availability than their calcareous counterparts (Jones, 2013). This study provides a case in which agglutinated foraminifera are somewhat tolerant to heavy metals, and a detailed understanding of their tolerance requires further investigation.

Porcelaneous foraminifera prefer natural, oligotrophic, unpolluted environments and cannot tolerate heavy metals (Elshanawany et al., 2011; Cosentino et al., 2013; Jorissen et al., 2018); this is roughly in accordance with the present study showing that most porcelaneous foraminifera were negatively related to heavy metals (see Figure 10B ). The exception is that two abundant species, Quinqueloculina laevigata and Q. seminula, exhibited a positive correlation with heavy metals, but a previous study in a polluted lagoon reflected a negative correlation (Frontalini et al., 2009). Perhaps the two species are flexible to adapt to the surrounding environment and change their autecological requirements to respond to environmental stress (Zettler et al., 2013; Bouchet et al., 2018).

The benefits and adverse effects of organic matter on benthic foraminifera are well documented (Alve, 1991; Armynot du Châtelet et al., 2009; Armynot du Châtelet et al., 2018): organic matter can serve as a food source for benthic foraminifera, but when its content exceeds the optimum value, the influence becomes adverse, mainly due to oxygen depletion under nutrient-enriched conditions. Thus, whether organic matter is beneficial or harmful is highly dependent on its content. In addition to the amount, the kind and quality of organic matter modifies foraminiferal assemblages (Murray, 2006). The calculated C/N ratios suggested that algae supplied organic matter to sediments, while this labile organic matter was bioavailable to benthic foraminifera (Martins et al., 2015). This could explain the positive response of porcelaneous taxa to the increase in TOC content. Furthermore, the TOC content was not high enough to cause pH lowering, which is particularly harmful for porcelaneous taxa (Li et al., 2014).

Much attention has been given to the influence of PAHs on benthic foraminiferal assemblages [e.g., Armynot du Châtelet et al. (2004); Bergamin et al. (2009); Romano et al. (2009); Schintu et al. (2015); Machain-Castillo et al. (2019), and Vilela et al. (2011)]; however, it is still hard to evaluate in a quantitative manner. Machain-Castillo et al. (2019) showed that very high PAH concentrations, e.g., above 5000 ng g⁻¹ (a level far higher than that in the study area), could cause large amounts of damage to foraminiferal assemblages. The problem is how to evaluate their effect at moderate concentrations. In this study, the impact of PAHs (LPAHs and HPAHs) was not as strong as that of heavy metals, probably because the PAH concentration was not sufficiently high. According to the species-environment relationship, porcelaneous taxa appeared to benefit from PAHs, while they were harmful to other species. The reason was yet unclear.

What thickness of the surface sediment layer is required to obtain a reliable dataset of benthic foraminiferal assemblages for biomonitoring studies? In many studies, researchers have chosen thicknesses of 1, 2, 4, and 5 cm [e.g., Frontalini and Coccioni (2008); Bergamin et al. (2009); Cosentino et al. (2013), and Martins et al. (2010)]. Barras et al. (2014) compared foraminiferal assemblages in the 0–1 cm layer and in the 0–4 cm layer and found no significant difference in community parameters and dominant species between the two intervals, although the latter’s density was much higher. However, the FOBIMO group strongly recommends the first centimeter of sediment for biomonitoring purposes because if an interval is too thick, surface-dwelling foraminifera will be diluted by a large number of dead individuals from deeper layers (Schönfeld et al., 2012). In this study, we chose the top 2 cm of sediments for foraminiferal analysis, while a more accurate investigation should be performed in the future.

This work emphasized the effects of pollution on benthic foraminiferal assemblages, including heavy metals, TOC, and PAHs. It should be mentioned that various natural factors, e.g., salinity, temperature, oxygen, nutrients, substrate, and food supply, also influence the distribution of benthic foraminifera (Murray, 2006). The influence of natural factors is not evaluated in this work and should be considered for future work. To measure the association between foraminiferal assemblages and environmental variables, we chose CCA as the primary analysis method. In CCA, environmental variables are used to describe the most important community gradients, but environmental variables are often measured for reasons of convenience rather than the “best” variables from the viewpoint of biology (McCune, 1997). Furthermore, environmental variables can be noisy or irrelevant and are sensitive to CCA (McCune, 1997).