Geochemical Control of PAHs by Inflowing River Water to West Nanao Bay, Japan, and Its Influences on Ecological Risk: Small-Scale Changes Observed under Near-Background Conditions at an Enclosed Bay

Polycyclic aromatic hydrocarbons (PAHs), even at low concentrations, have been shown to trigger changes in life cycles and provoke abnormal behaviors in numerous marine organisms. From May 2019 to September 2020, particulate and dissolved PAH concentrations were analyzed on the surface water of West Nanao Bay, Japan, to determinate their levels, emission sources, environmental pathways, and ecological risks at this remote but semi-enclosed bay. The 14 targeted PAHs were analyzed by HPLC-fluorescence detector. Mean total PAH concentrations were lower than 20.0 ng L−1 for most samples. Based on fluoranthene (Flu) to pyrene (Pyr) ([Flu]/[Flu + Pyr]) and benzo[a]anthracene (BaA) to chrysene (Chr) ([BaA]/[BaA + Chr]) isomeric ratios and a varimax rotated PCA, it was established that biomass combustion was the principal source in the particulate phase and that liquid fossil fuel combustion was the principal source in the dissolved phase. From salinity and turbidity distribution, riverine discharges were determined to be the major and continuous transportation pathway of particulate PAHs. It was observed that rain events had a role in the transport of dissolved PAHs. The risk quotients (RQ∑14 PAHs (NCs): 0–84.53) indicated that PAHs represented a very low to low acute environmental risk. The results of this study will contribute to filling the paradigm gap of ecotoxicological studies in remote areas, working as a booster for future in-lab studies of non-lethal implications of endocrine disruptors such as PAHs.


Introduction
Polycyclic aromatic hydrocarbons (PAHs), released to the environment mainly from the incomplete combustion of fossil fuels and biomass as well as from the leakage of petroleum derivates [1,2], are hazardous pollutants with adverse carcinogenic and/or mutagenic potential [3]. PAHs, even at low concentrations, have been shown to trigger changes in life cycles and to provoke abnormal behaviors in numerous marine organisms [4][5][6][7][8].
As mainly a hormone-regulated process, fish reproduction is very sensitive to physical and chemical stressors, such as PAHs [9]. Nicolas (1999) summarized the direct effects of PAHs in increasing vitellogenesis in various species of fish [4]. Recent studies have continued developing in this area [10][11][12][13][14][15]. Vitellogenesis in most species of fish occurs only once in the reproductive cycle. Consequently, even in non-lethal concentrations, PAHs have the potential to affect the whole reproductive cycle of the organism, leading to a gradual decrease in the population [4,[11][12][13][14]. Additionally, PAHs can provoke "behavioral toxicity" syndromes categorized by change in activeness and responsiveness [16]. The

Study Area and Water Sampling
West Nanao Bay is the innermost segment of Nanao Bay, isolated from Toyama Bay's waters by Noto Island. For this study, West Nanao Bay was divided into three subareas: west, central, and east areas ( Figure 1). The Kumaki, Otsu, and Ninomiya rivers flow into the west area of the bay, generating an estuary-like environment with shallow water (2-3 m). Fresh water discharges result in salinity that oscillates between 30 and 32 PSU. The seabed in this area is characterized by fine-grained (mud) sediments rich in organic matter and abundant in seagrass (Zostera marina) during the warm months. In the east area of the bay, the lack of seagrass, a sandier seafloor, and deeper waters (>10 m) are characteristic [24]. Salinity in the east area is usually over 32 PSU. Geologically, the central area presents the characteristic of being a transition of the previously described west and east subareas.
A set of 10 L surface seawater samples were collected from 15 points in West Nanao Bay ( Figure 1) every three months from May 2019 to September 2020. Additional samples (10 L) were taken in August 2020 and June 2021 at Kumaki, Otsu, and Ninomiya River mouths (salinity 0-15 PSU) and at the nearest fishing port to sampling point 11 (37.114475 • N, 136.883191 • E). Samples were collected and stored in polished stainless well-closed containers. In addition to the sample from the fishing port, no oil slicks were visible on the water surface during the surveys. Water quality, such as water temperature, salinity, and turbidity, were measured on-site with a direct-reading comprehensive water quality meter

Reagents and Chemicals
Extraction recoveries were quantitatively determined by an internal standard mix (ISTD) confirmed by two deuterated PAHs in phenanthrene (Phe-d10) and pyrene (Pyr-d10). Supelco EPA 610 PAH Mix was used as the reference standard for the target PAHs. All reagents were at least analytical reagent grade. Specific suppliers and purity grades were previously reported in Mundo et al. (2020) [28].  [28][29][30]. Briefly, particulate phase PAHs and dissolved phase PAHs were simultaneously separated by tandem filtration through 0.5 μm glass-fiber filters (GC50, Advantec, Hyogo, Japan) and C18 (octadecyl) disks

Reagents and Chemicals
Extraction recoveries were quantitatively determined by an internal standard mix (ISTD) confirmed by two deuterated PAHs in phenanthrene (Phe-d 10 ) and pyrene (Pyr-d 10 ). Supelco EPA 610 PAH Mix was used as the reference standard for the target PAHs. All reagents were at least analytical reagent grade. Specific suppliers and purity grades were previously reported in Mundo et al. (2020) [28].  [28][29][30]. Briefly, particulate phase PAHs and dissolved phase PAHs were simultaneously separated by tandem filtration through 0.5 µm glass-fiber filters (GC50, Advantec, Hyogo, Japan) and C18 (octadecyl) disks (Empore, St. Paul, MN, USA) at a flow rate of <200 mL min −1 ( Figure S1). After being air-dried, the glass fiber filters were dried over silica gel in the dark. Filters were stored at −20 • C until extraction with 100 mL dichloromethane. Dimethyl sulfoxide (DMSO; 100 µL) was added before concentration. Samples were finally reconstituted to 1 mL with acetonitrile. After the resulting solution was filtered with a membrane disk (HLC-DISK3, pore size 0.45 µm, Kanto Chemical Co., Tokyo, Japan), PAHs were separated and characterized by HPLC-fluorescence analysis. HPLC-fluorescence analysis of PAHs is described in detail in the Supplementary Material.

Quality Control
Method recoveries, including measurement sensitivities for HPLC assay of each sample, were estimated using internal standards of deuterated analogues of phenanthrene (Phe-d 10 ) and pyrene (Pyr-d 10 ), with which samples were spiked at the ultrasonic extraction stage for particulate PAHs and after seawater samples were filtered by C18 disks for dissolved PAHs. Overall recovery rates for dissolved PAHs' internal standards were 80.3 ± 13.3% for Phe-d 10 and 88.6 ± 14.2% for Pyr-d 10 . Recovery rates for particulate PAHs' internal standards were 95.3 ± 10.1% for Phe-d 10 and 96.6 ± 16.2% for Pyr-d 10 . To calculate the loss due to extraction procedures, recoveries of Phe-d 10 and Pry-d 10 were utilized for 3 rings-and 4-6 rings PAHs, respectively.
Background derived from pretreatment and measurements of PAHs were measured several times during the current research, and an average of all measurements were subtracted from sample concentrations. For each background measurement, blank tests were performed with 100 mL of ultra-pure water. The background concentrations in 10 L samples were Σ 14 PAH part = 0.21 ng L −1 and Σ 14 PAH diss = 0.27 ng L −1 .

PAH Source Estimation
PAH compositions vary among petrogenic, pyrogenic, biogenic, and diagenetic sources, although all of them can be found naturally in many areas of the world after specific events such as petroleum seeps and forest fires. Nowadays, petrogenic and pyrogenic events are mainly anthropogenic, releasing great amounts of PAHs to the environment and expanding their distribution ranges even to remote areas such as the Arctic and the Antarctic [31][32][33]. Their emissions origins, such as petrogenic and pyrogenic, are commonly estimated by using molecular indices based on ratios of PAH concentrations, according to the relative thermodynamic stabilities of each pair [34][35][36].
A [Flu]/[Flu + Pyr] ratio less than 0.4 indicates the predominance of petrogenic sources, while a range between 0.4 and 0.5 suggests the dominance of fuel combustion processes, and values greater than this suggest grass, coal, and wood combustion.
Principal component analysis (PCA) was additionally performed to clarify the discrepancies that arose from the diagnostic ratios method [38][39][40][41]. A varimax rotation and normalization was performed so the obtained data matrix would be comparable with those in previous studies. The rotation is commonly used to maximize the sum of the variances of the squared loadings within principal components and so to simplify the expression in terms of fewer PAHs each. The new data set obtained after varimax rotation contained fewer non-zero weights, which made the results easier to interpret and compare among studies.

Ecological Risk Assessment
The potential ecological risk of a given PAH was estimated by its risk quotients (RQ) for aquatic organisms [42][43][44][45]: where [PAH] i and QV i are the total concentration of the detected PAH in each water sample and the corresponding quality value, respectively. Two sets of quality values-the negligible concentrations (NCs) and the maximum permissible concentrations (MPCs)were calculated by Kalf et al. [43]. Therefore, for any given PAH, the values of RQ NCs(i) and RQ MPCs(i) can be defined as follows: where RQ NCs(i) and RQ MPCs(i) represent the quality value of the NCs and MPCs for the given PAH, respectively. Based on RQ NCs and RQ MPCs values, risk assignment parameters were determined [44]. Quality values for ten of the sixteen USEPA priority PAHs were originally reported [43]. For the remaining six, toxic equivalent factors (TEFs) were utilized to define the missing quality factors [46]. PAHs with similar TEFs should have similar negligible concentrations (NCs) and maximum permissible concentrations (MPCs).
To determine the overall risk of the 14 studied PAHs, the RQ values for ∑ 14 PAHs was calculated as defined by Cao et al. (2010) [45]:

Oceanographical Conditions
Salinity and turbidity average values at each subarea are summarized for each sampling season in Table 1. In Figures S2 and S3, the surface distribution of salinity and turbidity of the five sampling seasons are respectively presented. The west area, representing the riverine discharges, presented the lowest salinity in all five sampling seasons (Table 1, Figure 2). Turbidity presented an exact opposite trend, decreasing when moving away from the west area (Table 1, Figure 2). Tsushima warm current salinity varied between 31.6 and 35 PSU throughout the year, with salinity below 32 PSU in summer to autumn and above 34 PSU in winter to spring [47]. In West Nanao Bay, salinity was certainly lowest in summer (30.5-32.5), but it was never above 33.5 PSU (<32.5 in August, February, and September and <33.5 in May, November). As observed in Table 1, salinity in each subarea of the bay was lowest in summer and winter due to rainfall and snowfall, supported by the summer and winter monsoons [30]. The salinity trends observed in the present study support the observations of Yamazi, who found that oceanographical conditions in West Nanao Bay were largely controlled by riverine outflows and North Nanao Bay's seawater and were not directly affected by the Tsushima Warm Current [24].   The observed levels of PAHs at West Nanao Bay are comparable with other remote areas in Singapore, the northern part of Japan, and the Russian Far East, summarized elsewhere [28]. The 3-and 4-ring PAHs were the dominant species in dissolved PAHs, of which the sum of Fle, Phe, Flu, and Pyr occupied more than 80% of Σ14PAHdiss in all samples. For particulate PAHs, 4-, 5-and 6-ring PAHs were the major components, yet a more homogeneous distribution was observed in which each one of Flu, Pyr, BaA, Chr, BbF, BkF, BaP, DBA, BPe, and IDP occupied ±10% of Σ14PAHpart in all samples. The 5-and 6-ring PAHs were distributed mainly in the particulate phase, with almost negligible amounts in the dissolved phase. The difference in composition between particulate and dissolved PAHs is explained by the greater water solubility of PAHs with decreasing molecular size, due to the differences in water-octanol partition coefficients [29,34,48,49]. Examined by geographical subareas, the highest Σ14PAHs in the west area, decreasing towards the lowest Σ14PAHs in the east  The observed levels of PAHs at West Nanao Bay are comparable with other remote areas in Singapore, the northern part of Japan, and the Russian Far East, summarized elsewhere [28]. The 3-and 4-ring PAHs were the dominant species in dissolved PAHs, of which the sum of Fle, Phe, Flu, and Pyr occupied more than 80% of Σ 14 PAH diss in all samples. For particulate PAHs, 4-, 5-and 6-ring PAHs were the major components, yet a more homogeneous distribution was observed in which each one of Flu, Pyr, BaA, Chr, BbF, BkF, BaP, DBA, BPe, and IDP occupied ±10% of Σ 14 PAH part in all samples. The 5-and 6-ring PAHs were distributed mainly in the particulate phase, with almost negligible amounts in the dissolved phase. The difference in composition between particulate and dissolved PAHs is explained by the greater water solubility of PAHs with decreasing molecular size, due to the differences in water-octanol partition coefficients [29,34,48,49]. Examined by geographical subareas, the highest Σ 14 PAHs in the west area, decreasing towards the lowest Σ 14 PAHs in the east area, was also the general trend in August 2019, when abnormally high Σ 14 PAHs levels were detected in the central area of the bay (Figure 2).

PAH Levels
The all-points average Σ 14 PAHs was 11.01 ng L −1 , of which 43.9% was composed by the particulate phase, and 56.1% was the dissolved phase. In the west area, the all-points mean Σ 14 PAHs was 13.52 ng L −1 , of which the particulate phase occupied 58.1%, and the dissolved phase occupied 41.9%. In the central area, the annual mean Σ 14 PAHs was 11.42 ng L −1 , with 22.0% particulate phase and 78.0% dissolved phase. In the east area, the annual mean Σ 14 PAHs was 6.13 ng L −1 , with 30.2% particulate phase and 69.8% dissolved phase. Particulate PAHs were dominant in the west area, while in the central and east area, dissolved PAHs dominated. Mean Σ 14 PAH part was also highest in the west area, decreasing towards the central and east areas, although the mean Σ 14 PAH diss did not always follow such a pattern (Table 1, Figure 2). Analyzed by sampling survey, the particulate phase percentage of Σ 14 PAHs varied between 44% and 64% in all surveys except for August, when particulate PAHs occupied no more than 9.3%. area, was also the general trend in August 2019, when abnormally high Σ14PAHs levels were detected in the central area of the bay (Figure 2). The all-points average Σ14PAHs was 11.01 ng L −1 , of which 43.9% was composed by the particulate phase, and 56.1% was the dissolved phase. In the west area, the all-points mean Σ14PAHs was 13.52 ng L −1 , of which the particulate phase occupied 58.1%, and the dissolved phase occupied 41.9%. In the central area, the annual mean Σ14PAHs was 11.42   Figure 4. Diagnostic ratio averages are presented separately by sampling season (    It must be noted that some isomer pairs, such as in BaA-Chr, with considerable differences in photo-degradation rates can affect the reliability of the original source appointment method. Additionally, since the method is based on statistical analysis, the PAH sources directly interfering with the studied area may not be properly reflected within the set ranges.

Principal Component Analysis
PCA is one of the most common statistical methods for source apportionment of PAHs. PCA was separately performed for particulate and dissolved phases, with 59 samples and 14 PAHs each. The component's loading matrix is shown in Table 3. For each dissolved and particulate sample, two principal components were selected, so that 91.7% and 80.1% of the total variability would be explained in the dissolved and particulate phase, respectively. The first component (PC1) in the particulate phase is highly loaded by 4-6 ring PAHs (Pyr, Flu, BbF, BPe, BaA, Chr, BkF, Ant, and IDP). PC1 is responsible for 71.0% of the total variance. This distribution is characteristic of pyrogenic events fueled by biomass [39]. The second component (PC2) of particulate phase is responsible for 9.2% of total variance and is loaded by 3-ring PAHs (Ace, Fle) and one 5-ring PAH (DBA). The loading by low molecular weight PAHs is an indication of petrogenic sources. The PC1 for the dissolved phase is responsible for 84.5% of the total variance. The marked presence of 3-ring PAHs (Ace, Fle, Phe, Ant) would indicate a certain petrogenic character, but the presence of higher molecular weight PAHs is characteristic for fossil fuel-powered pyrogenic events [50,51]. PC2 for the dissolved phase, accounting for 7.2% of variability, is dominated by IDP plus Chr to a certain degree. Chr, according to Khuman et al. (2018), acts as a marker for coal combustion [38]. Additionally, a relatively (expressed as [IDP]/[IDP + BPe] ratio) high concentration of IDP would indicate pyrogenicity, coal burning in particular [36].

Environmental Pathways
The mean values of Σ 14 PAH part by area (west, central, east) are presented in Table 1. Σ 14 PAH part is highest in the west area, decreasing towards the central and east areas ( Figure 2). The west area also presents the highest mean turbidities and lowest mean salinities (Table 1, Figure 2, Figures S2 and S3). An inverse correlation was observed between turbidity and salinity. The r 2 values for the five seasonal sampling expeditions were 0.976, 0.626, 0.945, 0.393, and 0.864, respectively. The lowest r 2 values found in February can be attributed to snowfalls brought by winter monsoon. Less erosion and lower temperatures, limiting biological activity, would lower suspended solids concentrations, while melting snow would decrease salinity in a smoother gradient. To evaluate the correlation, turbidity versus Σ 14 PAH part means by area were plotted, giving correlation coefficients (r 2 ) of 0.999, 0.960, 0.947, 0.985, and 0.989 for each seasonal sampling expedition, respectively. Similarly, the anti-correlation of salinity versus Σ 14 PAH part was evaluated by the r 2 coefficients for their plots. The r 2 coefficients were 0.979, 0.427, 0.999, 0.516, and 0.927, corresponding to each sampling expedition. May, November, and September presented stronger anti-correlations between salinity and Σ 14 PAH part than August and February due to the overall lower salinity in all the areas of the bay due to summer and winter monsoon rainfalls.
The mean values of Σ 14 PAH diss by area (west, central, east) are presented in Table 1. Mean Σ 14 PAH diss variability among sampling areas did not always follow the west > central > east tendency shown by Σ 14 PAH part . The r 2 coefficient sets for the turbidity-Σ 14 PAH diss and salinity-Σ 14 PAH diss plots were (0.118, 0.393, 0.924, 0.0986, 0.0062) and (0.0373, 0.0004, 0.998, 0.877, 0.086), corresponding to each sampling expedition. With longer retention times, compared with suspended solids, once inside the bay basin, the Σ 14 PAH diss correlation with turbidity and salinity is less apparent [52,53].
August 2019 presented seven points with abnormally high dissolved PAH levels (>10 ng L −1 ) in the west and central areas. To solve the incongruity, seawater samples around the affected area were re-taken in the summer seasons that followed, together with water samples from the three river mouths (salinity 0-15 PSU) during both non-rain and rain events (rainfall: 52.5 mm day −1 ). Finally, the closest fishing port was also sampled.
Because the fishing port had relatively high Σ 14 PAHs (24.80 ng L −1 ) and Σ 14 PAH diss over 5 times Σ 14 PAH part , the dissolved phase had a petrogenic character; therefore, the fishing port was discarded as the main pollution source. Σ 14 PAHs diss from river water samples under normal conditions were all of a similar order of magnitude (12-21 ng L −1 ) and higher than most samples in the bay (Table S2). Thus, they can be considered as constant input concentrations during the summer season, but they were not necessarily responsible for the values recorded in August 2019. The Σ 14 PAHs in river water samples during rain events were an order of magnitude higher than average Σ 14 PAH diss in West Nanao Bay. Together with the overall higher fluxes, rain events brought considerably great amounts of PAHs into the bay system. Daily rainfall data are presented in Figure S4, where it can be seen how nearly 200 mm of rainfall was recorded during one week before sampling in August 2019. With Σ 14 PAH diss > 40 ng L −1 in all river water samples during rain events and the prolonged number of days of rainfall prior to sampling, the main pollution source in August 2019 is thus explained as being the result of runoff, followed by the respective water mixing-dilution. As demonstrated by Nakada et al., Kumaki and Ninomiya river plumes converge in the area between sampling points 11 and 3 during rain events (rainfall: 24 mm day −1 ) [54]. Despite no conclusive studies being found regarding PAH concentration variations within rain-induced runoff plumes, considerably higher PAH concentrations should be expected at headwaters after the initial "wash-out" of drieddeposited PAHs at the watershed. The Σ 14 PAHs diss found at points 11 and 3 were as high as "river water during rain event" samples. Furthermore, during the rain events, riverine plumes transported greater amounts of eroded particles. The greater particle densities and sizes accelerated scavenging, explaining the low particulate phase percentage of Σ 14 PAHs in August (only 9.3%) and the poorer anti-correlation between salinity and Σ 14 PAH part .

Ecological Risk Assessment
Results of the assessment of ecological risk in surface seawater of West Nanao Bay are presented in Table 4. Some sampled points in the east area, such as B1 and 8, had low RQ ∑14PAHs(NCs) values, even equal to 0 in some seasons, representing very low to negligible risk. In contrast, many sampled points on the west side present RQ ∑14PAHs(NCs) values greater than 20 more than once. Sampling point C3, located at the Otsu River mouth, had RQ ∑14PAHs(NCs) values of 20.68, 18.79, 53.80, 17.86, and 84.53 in the five sampling surveys, becoming consistently the most polluted point in the bay. Although sampling point 12 at the Ninommiya river mouth also presented high RQ ∑14PAHs(NCs) values, sampling point 2 at the Kumaki river mouth was one of the less contaminated points in the west area of the bay. Thus, not all river mouths are necessarily high-risk areas. Differences in RQ ∑14PAHs(NCs) values reflect the differences in land use of the three watersheds. RQ ∑14PAHs(MPCs) was less than 1 for all sampled points in all sampled seasons.  (28.37). RQ MPCs , on the other hand, were less than 1 for all PAHs in all samples.
Finally, it is known that many polar PAH metabolites-such as 4-hydroxybenz[a] anthracene; 3-, 4-, and 10-hydroxybenz[a]anthracenes; and 2-hydroxychrysene-have TEF values up to 10 times greater than original PAHs, in which case, ecological risks would be sub-estimated [19,23,55]. Although these and many other polar PAHs were not considered, the method proposed by Cao et al. provides a simple way to assess risk based solely on the concentrations and toxic equivalent factors of the 16 US EPA-designated PAHs, providing valuable information on the corresponding monitored sites.

Conclusions
From May 2019 to September 2020, particulate and dissolved PAH average concentrations were repeatedly analyzed in a four-seasons-long monitoring survey performed on surface waters of West Nanao Bay, Japan. Mean total PAH concentrations were lower than 20.0 ng L −1 for most samples. Based on [Flu]/[Flu + Pyr] and [BaA]/[BaA + Chr] isomeric ratios and a varimax rotated PCA, it was established that biomass combustion was the principal source in the particulate phase, and liquid fossil fuel combustion was the principal source in the dissolved phase. From salinity and turbidity distribution, riverine runoffs were determined to be a major and continuing transportation pathway of particulate PAHs towards the bay. River water samples suggest that dissolved PAHs had the same input routes as particulate PAHs but within rain-dominated pulses. The risk quotients (RQ ∑14PAHs(NCs) : 0-84.53) indicated that PAHs represented a very low to low environmental risk. BaA, BaP, and BPe were the three PAHs that presented single-component RQ NCs > 1 in the greater number of samples. The most polluted sample and the one with the highest ecological risk was consistently sampling point C3, located in the west area of the bay, where direct riverine runoffs are present.
This work will directly contribute to studies of the chronic effects of PAHs on environments with low contamination levels (e.g., wildlife protected coastal areas) within Japan and overseas. As reliable field data, paradigms such as negligible concentrations could continue to encroach bodies of water, justifying future research in the sublethal toxic effects of PAHs. Finally, the understanding of PAHs' environmental pathways in remote areas will contribute to improve remediation initiatives at the local level.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/ijerph181910310/s1, Figure S1: Illustration of tandem filtration system adaptation for simultaneously separate particulate PAHs and perform solid phase extraction of dissolved PAHs. Figure S2: Silinity seasonal changes in surface water of West Nanao Bay, Japan; August 2019-September 2020. * different scale, Figure S3: Turbidity (FTU) seasonal changes in surface water of West Nanao Bay, Japan; August 2019-September 2020. * different scale. Figure S4: Daily rainfall (mm) before sampling surveys at observation station Nishiyachi (Kumaki river upstream). Pluvial data obtained from Ishikawa Prefectural Civil Engineering Department River Division. Table S1: Concentrations of the targeted 14 PAHs in the particulate phase (Σ14PAHpart) and dissolved phase (Σ14PAHdiss) at the 15 seasonally sampled points in West Nanao Bay during 2019-2020, Table S2: Concentrations of the targeted 14 PAHs in the particulate phase (Σ14PAHpart) and dissolved phase (Σ14PAHdiss) at end members in West Nanao Bay, Table S3: PAHs ecological risk assessment guideline, Table S4: Toxic equivalent factor, negligible concentrations, and maximum permitted concentrations for USEPA 16 priority PAHs.