Nutrients and sea surface temperature drive harmful algal blooms in China’s coastal waters over the past decades

Eutrophication under climate change is well known to affect the community, productivity, and distribution of phytoplankton. However, the specific drivers of harmful algal blooms (HABs) in coastal China are not fully understood. Using HAB observed data since 1981, we quantified the distribution changes in HABs and estimated key environmental drivers (e.g. nutrients, sea surface temperature (SST), and precipitation) of HABs in China’s near seas. After 1981, the geographic range of HABs significantly expanded; moreover, the annual impacting period (AIP) of four of China’s near seas increased from 259% to 1090%. We found that rising total nitrogen (TN) or SST dominated the increase in the AIP in each near sea. Compared to the major contribution of TN to AIP in the other three near seas, TN has relatively weaker impacts than SST on AIP in the South China Sea (SCS). The peak of AIP in the SCS is highly correlated with extremely high SST. The significant contributions of climate change to HABs underscore the growing urgency to strictly control watershed nutrient input to mitigate marine eutrophication.


Introduction
Harmful algal blooms (HABs) are ecological phenomena caused by blooms of microscopic algae, phytoplankton, or macroalgae in seawater or freshwater, which can endanger public health and aquatic ecosystems (Anderson et al 2012).HABs threaten marine ecosystems, as they can damage or even be lethal to marine organisms and humans by producing toxins or depleting oxygen (Berdalet et al 2016, Plaas and Paerl 2021, Wang et al 2021a).In recent decades, the frequency, area, and distribution of HABs have increased significantly throughout China (Su 2001, Yu andLiu 2016).
Anthropogenic eutrophication caused by human activity, such as industrial expansion and agricultural practices, has led to drastic changes in nutrient structure along the coasts of China (Li et al 2014, Wang et al 2021b), which can further intensify the proliferation of HABs or even alter the structures of phytoplankton community (Glibert et al 2006, Li et al 2014, Glibert 2019).Meanwhile, specific physiological adaptive strategies of some phytoplankton (e.g.dinoflagellates) have permitted them to thrive by exploiting nutrients under high total nitrogen to total phosphorus ratio (TN:TP) conditions (Heisler et al 2008, Glibert et al 2014).Under complex hydrodynamic conditions, the increase and transformation of HABs along each regional coast responded differently to increased nutrient loads (Glibert et al 2018, Yu et al 2020).Of particular concern is the question of the extent to which changes in HAB intensity and distribution are contributed by nutrients to China's near seas.
In addition, climate change can also shape the spatial extent, duration, and toxicity of HABs (Thomas et al 2012, Gobler et al 2017).Some studies have shown that increasing sea surface temperature (SST) has facilitated potential growth rates and duration of bloom seasons of some toxic HABs (e.g.Alexandrium fundyense) (Moore et al 2015, Gobler et al 2017) by controlling phytoplankton metabolic processes and nutrient structure (Wells et al 2015(Wells et al , 2020)).Increased stratification as a result of warming can favour phytoplankton (e.g.dinoflagellates), due to their ability to migrate vertically to access nutrients (Wells et al 2015).Meanwhile, precipitation can alter terrestrial nutrient loads and concentrations by driving runoff, rates of fresh water inflow and outflow, influencing HAB changes (Wells et al 2015, Sinha et al 2017, Zhou et al 2019).
Consequently, in view of the nonnegligible risk caused by nutrient input and climate change (Hallegraeff 2010, Gobler 2019, Zeng et al 2019), it is imperative to identify the overall effect of such changes affecting HAB variability to control this severe marine disaster.The annual duration of HABs is an important indicator that can reflect the distribution and intensity changes of HABs over a long timescale (Gobler et al 2017, Gobler 2019).However, the overall effects of anthropogenic nutrient loading and climate change on the annual duration of HABs in China are limited.Therefore, based on HAB observed data since the 1980s in China, as well as the observed or reanalysis environmental variables, we aim (1) to evaluate the temporal and spatial distribution changes in HABs and (2) to quantify the dominant drivers of the annual impacting period (AIP) of HABs by constructing a statistical framework to further understand the driving mechanism of long-term HAB changes in China.
To include all potential riverine nutrient inputs to each coastal seas, we delineated seven watersheds (figure S1) and calculated the total discharge by summing the average riverine discharge for each watershed.The annual total nitrogen (TN, kgN km −2 yr −1 ) data was obtained from the PANGAEA website (https://doi.pangaea.de/10.1594/PANGAEA.892940)(Xu et al 2018).The watershed data were obtained from Resource and Environment Science and Data Center (www.resdc.cn/Default.aspx).Consequently, the TN value for each coastal sea represents the summation of average discharge from the outlets of all connected watersheds (Sinha et al 2017).Specifically, TN in the Haihe River watershed, Songhua and Liaohe River watershed, and Yellow River watershed was calculated for the BS.TN for the YS was obtained from Huaihe River watershed.For the ECS, TN in the Yangtze River watershed and Southeast watershed was calculated.For the SCS, TN in the Pearl River watershed was calculated.

Statistical analysis of HABs
The AIP is the number of HAB days occurring every year in near seas and can simply quantify change of HAB intensity over a long timescale.The AIP is an integer and discrete variable; it typically follows a Poisson distribution (Zhou et al 2021).Due to Bayesian Information Criterion (BIC) was used to select the combination of candidate variables with the best explanatory power for each of China's near seas (Zhou et al 2021).The BIC comprehensively considers both the number of variables and goodness of fit in candidate combinations and tends to select simpler combinations with higher precision (Anderson et al 1998): where y denotes the AIP in each near sea, AIP represents the estimated AIP of GLM for a given combination involving candidate variables, k is the number of variables in the combination, n is the number of data points in observed AIP, and T represents the transposition of the matrix.Here, eight candidate variables in each near sea are considered in the model (table 1).Additionally, only the combinations with at most three variables were selected as the optimal model to avoid overfitting of the AIP.The interactive effects of the variables in the optimal models can represent the complex environmental mechanism of HAB variability.
Once the variables are identified for the best model, the GLM was used to evaluate the relationship between the AIP and environmental variables of each near sea.
where β 0 represents the constant, X i is the variables in GLM, and β i denotes the regression coefficient of X i .GLM analysis was performed separately to identify the optimal predictor variables for explaining the AIP in each sea.
In this paper, extreme HABs are defined as those with long durations (>4 d) and large areas (>100 km 2 ).Only 10% of all HABs in the four near seas meet the criteria for extreme HABs.

Geographic extent variation of HABs
The geographic range of HABs has expanded significantly in China's four near seas over the past four decades (figure 1).Between 1981 and 1991, HABs were concentrated on regional coasts with intensive maricultural industries, such as coastal Shenzhen the Yangtze River Estuary, and the Dalian coastal sea (figure 1(a)).After 1991, HABs started expanding to the surrounding coastal waters.For example, HABs expanded from Laizhou Bay to coastal Hebei and the Liaoning, as well as the central BS (figure 1(b)).From 2001 to 2010, HABs rapidly expanded and covered most of China's coast, including the South of Yellow Sea (YS), the whole East China Sea (ECS), and the Guangxi coast of the South China Sea (SCS) (figure 1(c)).In particular, the ECS experienced over 500 HAB events, with high frequency (30 events) in the Yangtze River Estuary.After 2011, the frequency and area of HABs decreased significantly (−392 events and −117 319 km 2 , respectively), compared to 2001-2010.However, HABs expanded to the south of Hainan coast and to the west of coastal Guangxi in the SCS (figure 1(d)).
With the expansion of the geographic range of HABs in the four coastal seas of China, the dominant species of HABs also transformed to more toxic and harmful ones.Prior to 2001, Dinoflagellata

Changes in the AIP of HABs
For our analysis, the annual AIP of HABs increased significantly (p < 0.01) in the four near seas from 1981 to 2018 (figure 3).AIP increased at a rate of 2.2-4.   in the SCS, the results suggested that the fluctuation of nonextreme HABs has more important impacts on the expansion of the AIP (r = 0.94) of all HABs.

The key drivers of AIP changes in four of China's seas
For the nutrient-related variables, TN was selected as the predominant driver of interannual variability of the AIP in all the four China's near seas, indicating that increased nitrogen loading is a key factor for HAB proliferation at the multi-decadal scale (Wang et al 2021b) (table 2).Compared to the major contribution of TN on AIP in other three near seas (i.e. the BS, the YS, and the ECS), TN has relatively weaker impacts on AIP in the SCS.
For the climate change-related variables, SST a4-9 dominated the expansion of the AIP in the SCS, causing a prolonged blooming season and rising growth rates especially for the dominant phytoplankton groups of Ochrophyta and Dinoflagellata (Gobler et al 2017, Griffith andGobler 2020)

Discussion
The significant role of anomalous SST on HABs of China's near seas at the multi-decadal scale was identified through our model.Comparing to studies in the Chesapeake Bay, North Atlantic, and North Pacific oceans (Anderson et al 2010, Cai et al 2016, Zhou et al 2021), the fact that SST a4-9 is selected over SST su or SST asu for AIP models suggests that the proliferation of some algae (e.g.dinoflagellates) is more promoted by the duration of the warm season rather than higher SST su (Gobler et al 2017, Wells et al 2020).The results confirm that increasing SST significantly expand spatial extent and frequency of HABs in mid and high latitude coastal waters in the 21st century (Griffith et al 2019, Dai et al 2023).Although reduced nutrient loadings after 2010 lead to decreased area and frequency of HABs (Zeng et al 2019), the AIP increased (except the ECS) and almost peaked with abnormally high SST a4-9 since 2010 (figures 4-7).Reduced application of fertilizer as implemented in other regions (e.g.Europe, US, and Japan) (Dai et al 2023) would not effectively control the HAB frequency, since the benefits form nutrient control might be counterbalanced by climate warming in China coasts.Especially  et al 2018).This work highlighted that SST a4-9 could dominantly drive a peak of AIP, which will have profound implications for designing marine governance, especially in the SCS, as climate warming is dramatically transforming marine ecosystems.
Our results highlight the mixed effects of precipitation on AIP changes along the coasts of China.The contributions were dependent on the aggregation period (annual or seasonal) of Pre and its patterns (total or extreme) (Ho and Michalak 2020).Consistent with the load driving hypothesis (i.SCS has the highest SST, causing the longest duration for HAB growth.Specifically, the SCS is the largest marginal sea in the western tropical Pacific Ocean, with a maximum depth reaching 5000 m, and its continental slope is relatively steep (Shaw and Chao 1994).Ocean circulation and mixing induced by monsoon and Kuroshio intrusion (Gan et al 2006) cause significant water exchange that predominantly regulates the phytoplankton and nutrient dynamics in the SCS (Shaw and Chao 1994).Additionally, this highly dynamic water system constrained the growth and proliferation of HABs to some extent, especially extreme HABs (figure 3(d)), as the algae could be more easily advected from nearshore to offshore region that with relatively lower nutrient supply and higher zooplankton grazing pressure (Feng et al 2022).Conversely, the ocean currents are conducive to proliferation of dinoflagellate in BS, YS, ECS, Japan and South Korea coasts by boring more nutrients from the continental shelf (Griffith et al 2019, Zhou et al 2019).
By investigating statistical relationship of the AIP and environmental variables, the study further explains how the interactions of nutrient and climate changes drive the interannual variability of AIP along the coasts of China.This study improves our understanding of the anthropogenic perturbation co-stressor mechanism on the long-term HAB changes.However, due to observational data limitations, the complex, dynamic, and nonlinear biological responses of different HAB species to environmental factors were not considered in this work.In addition, the exploration of the environmental driving mechanisms on HAB changes for different algal species is limited.An additional consideration in the assessment of how co-stressors of anthropogenic perturbations affect HABs that has yet to be studied is the extent to which phytoplankton may acclimate their growth response to changing co-stressors.With gradually increasing marine monitoring and data sharing programs, future modelling studies could focus on evaluating the biological responses of different phytoplankton species to the overall effects of eutrophication and climate warming to further quantify the driving mechanism of HAB variability.

Conclusions
Overall, the AIP of HABs have expanded by 259%-1090% in China's four near seas since 1981.The GLM results demonstrate that the SST dominantly increases AIP in the SCS, while TN predominantly promote the expansion of AIP in the other three seas.Our analysis revealed that water quality management strategies including nutrient loading control should consider the large impacts of climate change on proliferation of HABs.

Figure 2 .
Figure 2. Changes in the dominant species of HABs in four of China's near seas since 1981.The frequency (bars) and area (red lines) changes in the dominant species of HABs during (a) 1981-2000 and (b) 2001-2018.

Figure 3 .
Figure 3. in AIP of HABs in the four near seas from 1981 to 2018.AIP of all HABs and extreme HABs in the (a) Bohai Sea, (b) Yellow Sea, (c) East China Sea, and (d) South China Sea.r is the correlation coefficient of the AIP of all HABs and extreme HABs in each sea.
8 d per year for these coasts, with the highest increase occurring in the YS, indicating that AIP increased by 2-6 months during 1981-2016.In the YS, Chlorophyta with long duration (>50 d) caused the AIP to increase rapidly to above 65 d after 2008 (figure 3(b)).In 2016, two extreme HABs that lasted over 120 d caused the AIP peaked at 226 d.Additionally, the increased AIP in the BS, peaked at 130 d in 2013, due to Ochrophyta that lasted 99 d (figure 3(a)).The AIP in the ECS increased to above 50 d starting in 2002, but it fell back after reaching its peak (162 d) in 2007 (figure 3(c)).However, For the SCS, the AIP has increased since 1981, with a peak (92 d) occurring in 1987, followed by a decrease.Subsequently, the AIP again increased starting in 1995, and exceeded 100 d during 2014-2016 (figure 3(d)).These observations evidenced extreme HABs with more than 60 d seem to play a vital role in the peak of the AIP.Except for the SCS, AIP of extreme HABs was highly correlated with that of all HABs (r > 0.85) in the other three seas (figures 3(a)-(c)), implying that duration changes of extreme HABs can effectively drive the AIP changes of all HABs.Furthermore, the extreme HABs that lasted over 90 d directly exacerbated AIP peak after 2000, particularly in the YS.However, a low correlation (r = 0.54) between the AIP of extreme HABs and that of all HABs was reported in the SCS (figure 3(d)).Compared to the low correlation of extreme HABs on the AIP of all HABs

Figure 4 .
Figure 4. Observed AIP, GLM simulated AIP, and individual variables included in the GLM of the Bohai Sea during 1981-2016.
Figure 5. Observed AIP, GLM simulated AIP, and individual variables included in the GLM of the Yellow Sea during 1981-2016.

Figure 6 .
Figure 6.Observed AIP, GLM simulated AIP, and individual variables included in the GLM of the East China Sea during 1981-2016.
Figure 7. Observed AIP, GLM simulated AIP, and individual variables included in the GLM of the South China Sea during 1981-2016.

Table 1 .
Candidate variables used to explain the interannual variability of the AIP of HABs. the watersheds of each near sea (10 3 kg N km −2 yr −1 ) DIN Dissolved inorganic nitrogen concentration of certain sea areas (mg l −1 ) DIP Dissolved inorganic phosphorus concentration of certain sea areas (mg l −1 ) WD, VWD, UWD Average wind speed, meridional wind speed, zonal wind speed over warm season of each near sea (m s −1 ) back of HABs and environmental variables through the logarithm link function(Li et al 2014, Zhang et al  2016).Therefore, to explore the responses of the AIP to nutrients and climatic variables (e.g.SST, wind, and Pre), we formulated a GLM to quantify driving factors of the AIP interannual variability in each of China's near seas.

Table 2 .
The estimated coefficient and associated standard of each variable in the selected GLM in China's near seas during 1981-2016.All variables were normalized for use in the multiple linear regression model by subtracting the 1981-2016 mean for that year and dividing by the standard deviation.