Phytoplankton Species Diversity Patterns and Associated Driving Factors in China’s Jiulong River Estuary: Roles That Nutrients and Nutrient Ratios Play

Understanding diversity patterns and associated driving factors are the critical topics in macroecology and conservation biology. Phytoplankton are highly susceptible to environmental changes in estuaries, particularly eutrophication. This study examined phytoplankton alpha and beta diversity using investigation data in May (springtime), August (summer) and November (autumn) 2009 in China’s Jiulong River estuary, where it was easily polluted because of considerable discharge from a highly dense human population and low self-purification capacity with its limited river basin area, potentially resulting in eutrophication and then influencing phytoplankton diversity. Potential influencing factors were also explored, including dissolved oxygen, salinity, nutrients, nutrient ratios, geographic and hydrologic distance, and so on. The results indicated that Shannon’s index (H’) and Pielou’s index (J) decreased from the estuary’s upper to middle and then increased from middle to lower reaches, Simpson’s (D) observed the opposite trend and species number (S) gradually increased from the estuary’s upper to lower reaches. For beta diversity, all the indices showed a gradual decrease trend from the estuary’s upper to lower reaches, where also, turnover dominated beta diversity for all seasons. It is noteworthy that the significant roles that nutrients and nutrient ratios played in shaping phytoplankton diversity patterns and the nutrient balance were characterized by excess nitrogen (N) and silicon (Si) and limited phosphorus (P), which could potentially cause diatom blooms. Findings also showed that decreasing Si concentrations can help to reduce overall pollution levels as well as the restoration of the estuary’s ecosystem better than just reducing N alone. Accordingly, this study advocates for the protection of the entire estuary system with particular emphasis on its upper reaches. Moreover, greater attention should also be paid to impacts associated with N input and nutrient ratio trade-offs to the prospective watershed management of this estuary. This study provides a practical approach to explore estuarine diversity in a comprehensive way, which can inform effective biodiversity conservation and also be applied to other marine ecosystems to better guide sustainable management and conservation practices.

Understanding diversity patterns and associated driving factors are the critical topics in macroecology and conservation biology. Phytoplankton are highly susceptible to environmental changes in estuaries, particularly eutrophication. This study examined phytoplankton alpha and beta diversity using investigation data in May (springtime), August (summer) and November (autumn) 2009 in China's Jiulong River estuary, where it was easily polluted because of considerable discharge from a highly dense human population and low self-purification capacity with its limited river basin area, potentially resulting in eutrophication and then influencing phytoplankton diversity. Potential influencing factors were also explored, including dissolved oxygen, salinity, nutrients, nutrient ratios, geographic and hydrologic distance, and so on. The results indicated that Shannon's index (H') and Pielou's index (J) decreased from the estuary's upper to middle and then increased from middle to lower reaches, Simpson's (D) observed the opposite trend and species number (S) gradually increased from the estuary's upper to lower reaches. For beta diversity, all the indices showed a gradual decrease trend from the estuary's upper to lower reaches, where also, turnover dominated beta diversity for all seasons. It is noteworthy that the significant roles that nutrients and nutrient ratios played in shaping phytoplankton diversity patterns and the nutrient balance were characterized by excess nitrogen (N) and silicon (Si) and limited phosphorus (P), which could potentially cause diatom blooms. Findings also showed that decreasing Si concentrations can help to reduce overall pollution levels as well as the restoration of the estuary's ecosystem better than just reducing N alone. Accordingly, this study advocates for the protection of the entire estuary system with particular emphasis on its upper reaches. Moreover, greater attention should also be paid to impacts associated with N input and nutrient ratio trade-offs to the prospective watershed management of this estuary. This study

INTRODUCTION
Reducing biodiversity loss has become a global issue on both land and sea (Blowes et al., 2019;Pinsky et al., 2019;Cooper et al., 2020;Trisos et al., 2020). The Earth is presently experiencing an extinction crisis that is largely the result of anthropogenic exploitation, which has caused widespread changes in the global distribution of organisms (Chapin et al., 2000), where roughly one million species are threatened with extinction (Balvanera, 2019). The biodiversity of ocean systems is higher than terrestrial systems, while the rate of loss in marine biodiversity occurs much quicker compared to terrestrial biodiversity caused by climate change and anthropogenic impacts (Blowes et al., 2019;Pinsky et al., 2019;Cooper et al., 2020;Trisos et al., 2020). Furthermore, a decline in marine biodiversity has caused to reduction or loss of some ocean or coastal ecosystem services and functions (Blowes et al., 2019;Pinsky et al., 2019;Cooper et al., 2020;Paul et al., 2020;Trisos et al., 2020).
Estuaries act as transitional zones between brackish and freshwater, where intense material exchanges take place among terrestrial, riverine, and marine ecosystems, subsequently playing a critical role in biogeochemical cycling (Srichandan et al., 2019). High diversity and productivity are critical characteristic of estuarine systems (Dursun and Tas, 2019), which comprise of various ecosystem services and a high natural capital value (Costanza et al., 1997). However, estuaries are susceptible to a variety of natural and anthropogenic disturbances (Gillanders et al., 2011;Atwood et al., 2012;Bucci et al., 2012;O'Brien et al., 2016). It has been estimated that greater than 60% of the world's population and two-thirds of mid-to large-sized cities are distributed within and around estuarine and coastal areas (Ma and Mao, 2014;O'Brien et al., 2016). Over the past half century, intensive anthropogenic disturbances have dramatically increased nutrient inputs into estuaries from inland sources, causing widespread eutrophication (Nixon, 1995;NRC (National Research Council), 2000;Rabalais, 2002;Ho et al., 2019) that ultimately results in the degradation of habitats and a decrease in biodiversity (NRC (National Research Council), 2000;Bricker et al., 2008;Zhao et al., 2010;Clark et al., 2015). Accordingly, nutrient pollution in estuaries and associated effects have garnered increasing attention throughout the world (Liu and Shen, 2001;Li et al., 2014;Paerl et al., 2014;Wallace et al., 2014;Fennel and Testa, 2018;Barletta et al., 2019;Wurtsbaugh et al., 2019).
Being both the most important primary producer and the fundamental component of the food web, phytoplankton are an essential constituent of marine ecosystems (Dursun and Tas, 2019). However, phytoplankton is susceptible to nutrient concentrations (Redfield, 1934;Paerl et al., 2007), whose biomass and patterns will be dramatically altered in the presence of such nutrients. To some extent, the composition, structure, function, and spatiotemporal distribution of phytoplankton indicate the environmental quality status and the biological integrity of aquatic ecosystems (Ekwu and Sikoki, 2006). In macroecology and conservation biology, it is critical to understand phytoplankton spatial diversity patterns (Collen et al., 2013). In doing so, previous studies primarily focused on alpha and gamma diversity (Muylaert et al., 2009;Rochelle-Newall et al., 2011;Guo et al., 2014;Bharathi and Sarma, 2019;Dursun and Tas, 2019), while little attention has been paid to beta diversity. Although both alpha and gamma diversity are indicative of species richness, their spatial scales dramatically differ, where the former refers to in-site diversity (i.e., local species pools) and the latter refers to overall diversity (i.e., regional species pools). Notably differing from alpha and gamma diversity, beta diversity refers to between-site diversity across different habitats, ecosystems, or regions (Whittaker, 1960), which reflect community construction processes as well as biodiversity formation and maintenance (Condit et al., 2002;Kraft et al., 2011;Legendre and De Cáceres, 2013;Li et al., 2016). Given the considerable habitat gradients of estuarine systems (e.g., that associated with salinity levels), beta diversity (used in conjunction with alpha and gamma diversity) and associative patterns provide much more informative results. This approach is key to understanding estuarine response to anthropogenic activities.
The Jiulong River estuary is one of China's most important estuarine systems. It is characterized by an extremely dense human population (13.39 million) that subsists within a limited river basin area (1.47 × 10 4 km 2 ) (Fujian Provincial Bureau of Statistics, 2020), which roughly corresponds to approximately 41% of the areas of the Hudson River basin, Canada, under a population density equivalent to New York City, United States. A highly dense human population in conjunction with highly intensive anthropogenic activities pose great pressure on the Jiulong River estuary, where annual dissolved inorganic nitrogen (DIN) flux is 8.278 × 10 4 t·yr −1 , soluble reactive phosphorus (SRP) flux is 0.195 × 10 4 t·yr −1 , and dissolved silica (DSi) flux is 17.831 × 10 4 t·yr −1 (Wu et al., 2017), which may shape and influence phytoplankton diversity patterns, given that phytoplankton are the most sensitive aquatic community to nutrient pollution. This study used the Jiulong River estuary as the study area to (i) examine spatial patterns of phytoplankton species diversity along estuary to offshore gradients, using both alpha and beta diversity; to (ii) explore potential factors that shape diversity patterns, particularly within an environmental context, and to test whether differences exist between alpha and beta diversity; and to (iii) evaluate the relative roles that nutrient concentrations and nutrient ratios play in shaping diversity patterns.

Study Area
The Jiulong River estuary is situated along China's southeastern coast at the southern point of the East China Sea. This estuary originates from Longyan and Zhangzhou cities, with a total length of 1,923 km and a basin area of 1.47 × 10 4 km 2 (Figure 1; Huang, 2008). Nutrient and pollutant runoff flows into the estuary through three tributaries, namely, the Beixi River, the Xixi River, and the Nanxi River, which together strongly impact the Jiulong River estuary. The estuary is comprised of several ecosystems, such as mangrove, salt marsh, tidal flats, and brackish water, which provide important spawning ground habitats for many marine organisms. In 1988, the Jiulong River Estuary Provincial Mangrove Reserve (with a total area of 420.2 ha) was established to protect its mangrove ecosystem, wetland waterfowl community, and endangered wildlife. The Jiulong River estuary is an also important habitat for the Indo-Pacific humpback dolphin (Sousa chinensis), which the IUCN Red List of Threatened Species has classified as "Vulnerable." The estuary also acts as a migratory "stopover" for thousands of migratory waterfowls in the East Asian-Australasian Flyway.

Methodological Approach
The study approach is provided in Figure 2. First, species diversity indices were selected and calculated, including alpha and beta diversity. Then, a quantile-quantile test (i.e., Q-Q test) was used to verify data normality of alpha diversity, while Ordinary Kriging was used to map distribution. At the same time, the study area was divided into different subregions based on salinity clustering and geographical characteristics, and beta diversity distribution was mapped for pairwise site and subregion comparisons. Finally, Pearson's r was applied to calculate the relationship between environmental variables and diversity indices, after which optimal equations were obtained by means of stepwise regression, which was used to identify critical environmental factors that contribute to spatial patterns in phytoplankton diversity.

Data Sources
The phytoplankton and environmental data used in this study were obtained from Yu et al. (2012). These data were collected at 21 locations along the estuary during three field campaigns, namely, in May (springtime), August (summer), and November (autumn) 2009 (Figure 1). All samples were obtained during low tide. Niskin bottles (5 L) was used to collect water samples, where two water depths were conducted, i.e., surface samples at a water depth below 5 m and bottom samples at a water depth above 5 m. Three duplicate samples were also collected from each location.
Phytoplankton organisms contained in the samples were then transferred to glass sampler bottles and fixed in a 2% glutaraldehyde solution. Referring to traditional phytoplankton identification systems (Jin and Chen, 1960;Hu et al., 1980;Tomas, 1997), species identification and phytoplankton enumeration were conducted using an epifluorescence microscope. The original phytoplankton species data supporting Yu et al. (2012) were provided in the Supplementary Table 2.
The environmental data of 15 variables were extracted from the Figures of Yu et al. (2012), including pH, dissolved oxygen (DO), chemical oxygen demand (COD), temperature (T), salinity (Sal), water depth (WD), Total nitrogen (TN), DIN (i.e., the sum total of ammonium, nitrate, and nitrite), silicate (SiO 4 − -Si), nitrate (NO 3 − -N), nitrite (NO 2 − -N), ammonium (NH 4 + -N), total phosphorus (TP), dissolved inorganic phosphorus (DIP), suspended particulate matter (SPM). A multi-parameter probe (YSI 6600 V2) was used to measure pH and DO. A conductivitytemperature-depth (CTD) system was used to measure T, Sal, and WD. Samples used for nutrient concentration measurements were kept in 2 L plastic bottles at −4 • C and measured within a 12 h period at a laboratory. TN, DIN, SiO 4 − -Si, NO 3 − -N, NO 2 − -N, NH 4 + -N, TP, DIP, SPM, and COD were measured following the National Standard of the Oceanographic Survey of China (GB/T 12763.4-2007). Other environmental variables were also calculated based on the above procedure, including the dissolved inorganic nitrogen (DIC) to phosphorus (P) ratio (N/P), the TN to TP ratio (TN/TP), the silicon (Si) to P ratio (Si/P), and the Si to DIN ratio (Si/N).
Hydrologic data used to calculate the hydrologic distance was extracted from Jing et al. (2011), including hydrologic direction and velocity.

Statistical Analysis
The phytoplankton diversity indices selected for this study cover both alpha and beta diversity. For alpha diversity, four typical and widely used indices were selected, including the Shannon index (H ), Simpson's index (D), Pielou's evenness index (J), and the Species number (S) index. For beta diversity, this study used the  BAS approach (Baselga, 2010(Baselga, , 2012 to calculate and decompose beta diversity, where the Sørensen pairwise dissimilarity index (beta.SOR) was first calculated and then decomposed to a turnover value (beta.SIM) and a nestedness value (beta.SNE). Supplementary Table 1 provides the diversity index equations, which were conducted in the packages "vegan" and "beta.part" of R 4.0.2 (Baselga and Orme, 2012).
Spatial subdivision was conducted to further quantify the diversity distribution. The salinity clustering results showed that five distinctive groups were found according to their geographic Frontiers in Marine Science | www.frontiersin.org location, i.e., tributary (S1-S4, S21), river mouth (S5, S6), southern part of middle estuary (S7, S8, S15), northern and eastern part of middle estuary (S9, S10, S13, S14, S16, S17) and outer estuary (S11, S12, S18-S20) (Figure 3). Based on salinity clustering groups, combined with the consideration of sampling sites number balance, the study area was subdivided into three subregions, including the upper (Area I), middle (Area II), and lower reaches (Area III) of the estuary (Figure 4).
Results from the Q-Q test indicated that alpha diversity data followed normal distribution (see Supplementary Figure 1), suggesting that Ordinary Kriging was applicable in mapping the distribution of alpha and gamma diversity, where a linear function was selected as the semi-variant model used in Ordinary Kriging, given that the linear function was optimal for accuracy of cross-validation evaluation. A pairwise site comparison of beta diversity was mapped using a connective line, and averages were used to represent the beta diversity of each subregion. The "ggplot2" package in R 4.0.2 was used to conduct the Q-Q test. ArcGIS 10.7 was used to conduct Ordinary Kriging and biodiversity mapping.
Pearson's r was used to quantify relationships between environmental and geographic parameters and diversity indices, which was conducted to identify the driving factors of phytoplankton diversity patterns in the Jiulong River estuary. In total, 63 datasets obtained from three field campaigns were used to analyze correlations between environmental variables and alpha diversity. In total, 210 datasets, covering data from three field campaigns that were mixed among any two sites, were used to analyze correlations between beta diversity and environmental and geographic variables, where differences in environmental variables among pairwise site comparisons were used, while geographic distance and hydrologic distance were used to represent geographic variables, which were, respectively, obtained using the point distance tool of proximity analysis and SDM toolbox incorporating the north and east components of mean flow velocity at 6 h during low tide (because of a regular semidiurnal tide in Jiulong River estuary, Liu, 2012) as the cost resistance surface in ArcGIS 10.7. Additionally, based on Pearson's r results, stepwise regression was conducted to obtain optimal equations between variables and diversity indices after data was normalized. This was done to enhance comparisons at the same scale. Forest plots were used to visualize the coefficient strength of key variables in the optimized equation. Pearson's r, normalization, stepwise regression, and forest plots were conducted in packages "corrplot, " "clusterSim, " "car, " "grid, " "magrittr, " "checkmate, " and "forestplot" in R 4.0.2, respectively.

Spatial Patterns of Phytoplankton Alpha Diversity
Supplementary Table 3 and Figure 5 provide the overall spatial variation of phytoplankton alpha diversity indices across three seasons. Shannon index (H') values ranged from 2.23 to 3.51 (2.90 average), where the highest value was measured at S11 (i.e., at the lower reaches of the estuary) and the lowest value was measured at S8 (i.e., at the middle reaches of the estuary). Simpson's index (D) values ranged from 0.15 to 0.71 (0.28 average), where the highest value was measured at S1 (i.e., at the upper reaches of the Beixi River estuary) and the lowest value was measured at S3 (i.e., at the upper reaches of the Xixi River estuary). Pielou's evenness (J) values ranged from 0.55 to 0.84 (0.67 average), where the highest value was measured at S3 (i.e., at the upper reaches of the Xixi River estuary) and the lowest value was measured at S20 (i.e., the lower reaches of the estuary). Species number (S) values ranged from 27 to 59 (45 average), where the highest value was measured at S19 (i.e., the lower reaches of the estuary) and the lowest value was measured at S4 (i.e., the upper reaches of the estuary). Results generally showed a wide range of values in all alpha diversity indices, indicating significant spatial variation among the different sites. A similar gradient was observed in H' and J values, which gradually decreased from upper to middle and then increased from the middle to lower reaches of the estuary (Supplementary Table 3 and Figures 5A,B). D gradually decreased from the upper to lower reaches of the estuary (Supplementary Table 2 and Figure 5C). Conversely, S gradually increased from the upper to lower reaches of the estuary (Supplementary Table 3 and Figure 5D). Among the different subregions, the highest H' and S values were observed in Area III, the highest J value was observed in Area I, and the highest D was observed in Area II.
Season variations of spatial patterns were observed (Supplementary Table 3). The highest H' value was found in autumn (3.43 average), followed by spring (2.79 average) and summer (2.46 average), where spatial patterns varied greatly across different seasons, i.e., H' value in spring decreased from upper to lower reaches of estuary, while the opposite trend was found in autumn. The highest J value was found in autumn (0.75 average), followed by spring (0.70 average) and summer (0.56 average), where their spatial patterns were similar with H'. The highest D value was found in summer (0.42 average), followed by spring (0.24 average) and autumn (0.19 average), where their spatial patterns were nearly opposite to H' and J. The S values in autumn (24 average) and summer (23 average) were higher than that in spring (16 average), with an increasing spatial trend from upper to lower reaches of estuary for all three seasons.

Spatial Patterns of Phytoplankton Beta Diversity
Pairwise beta diversity site comparisons are shown in Figure 6. The beta.SOR values ranged from 0.27 to 0.81 (0.4 average), and this large range indicated significant variation between the different sites. A gradual decreasing trend was observed from the upper to lower reaches of the estuary. The highest beta.SOR value was observed between S3 (the Xixi River) and S21 (the Nanxi River) at the upper reaches of the estuary, characterized by a low hydrologic exchange between these two different tributaries. On the other hand, the lowest beta.SOR value was observed between S10 and S20 at the lower reaches of the estuary, characterized by a small geographic distance and a high hydrologic exchange.
Supplementary Table 3 show the beta diversity of each subregion and associated components. beta.SOR values ranged from 0.34 to 0.51 (0.4 average); beta.SIM values ranged from 0.31 to 0.43 (0.35 average); beta.SNE values ranged from 0.03 to 0.07 (0.05 average). beta.SOR, beta.SIM, and beta.SNE values exhibited the same distribution pattern, namely, a deceasing trend from the upper to lower reaches of the estuary, where the highest value was measured at Area I and the lowest value was measured at Area III. Furthermore, beta.SIM dominated beta.SNE regardless of the area (i.e., the upper, middle, or lower reaches of the estuary), indicating that species turnover was dominant.
The seasonal similarities and differences of beta diversity patterns were provided in the Supplementary Table 3. The same pattern was observed for beta.SOR and beta.SIM, where temporally, their values in summer and spring were slightly higher than that in autumn, and spatially, the highest value was found at Area I and the lowest value was found at Area III, indicating a deceasing trend from the upper to lower reaches of the estuary for all three seasons. Fluctuations in beta.SNE values differed from the other two indices, it showed inconsistent spatial trends across different seasons, where the highest value was found in Area III in spring, while Area II observed the highest value in summer and autumn. It was noted that beta.SIM value was greatly higher than beta.SNE across all subregions and seasons, indicating that species turnover was dominant.
The optimal regression equations between phytoplankton beta diversity indices and environmental and spatial variables were subsequently developed, and their parameters and coefficients are shown in Supplementary

Diversity Patterns Incorporating Alpha and Beta Diversity
Phytoplankton species diversity in the Jiulong River estuary was calculated at two different scales, namely, alpha and beta diversity, where alpha diversity reflects with-in site diversity and beta diversity reflects between-site diversity between sites or subregions. Different diversity scales represent fundamental aspects of biological communities (Socolar et al., 2016), and this is due to the different constructs and mathematical expressions among the various indices.
For overall alpha diversity across three seasons, H' and J values decreased from the upper to middle and then increased from middle to lower reaches of the Jiulong River estuary. D exhibited the opposite trend, and S gradually increased from the upper to lower reaches of the Jiulong River estuary Alpha diversity patterns were observed to be complex, and this was due to their multiple indices, which differ in their theoretical foundation and interpretation (Morris et al., 2014). Moreover, H' was equally sensitive to rare and abundant species; J emphasized evenness; D was sensitive in the presence of abundant species; S was sensitive to rare species (Shannon, 1948;Simpson, 1949;Pielou, 1969;Peet, 1974;Magurran, 1988;Krebs, 1998;Morris et al., 2014). Observed patterns in the different diversity indices indicated that species abundance was relatively evenly distributed in the upper reaches of the estuary, a greater number of abundant species was observed in the tributary estuary, while a greater number of rare species was observed in the lower reaches of the estuary.
For overall beta diversity across three seasons, beta.SOR, beta.SIM, and beta.SNE values exhibited a general gradient, gradually decreasing along with an increase in salinity gradients from the upper to the lower reaches of the estuary. This trend also reflects the strength of environmental gradients under the hypothesis that diversity is greater in heterogeneous environments. In other words, the standard deviation of most environmental variables in the upper reaches of the estuary was higher than corresponding values in the middle or lower reaches, including T, DO, SPM and all nutrients and nutrient ratios except for NO 2 − -N (Supplementary Table 4). Analogously, beta diversity patterns of the three subregions showed that their variations gradually decreased from the upper reaches (Area I) to the middle reaches (Area II) to the lower reaches (Area III) of the estuary (beta.SOR: I−II = 0.146, I−II = 0.024; beta.SIM: I−II = 0.114, I−II = 0.014; beta.SNE: I−II = 0.032, I−II = 0.01). Results indicated that phytoplankton communities in the hyposaline upper reaches of the estuary were the most diverse and heterogeneous in species composition, while hypersaline communities were the least diverse and most homogeneous. Besides, there was not seasonal difference in spatial pattern of beta.SOR and beta.SIM, which was consistent with the decreasing gradient aforementioned from the upper to the lower reaches of the estuary. Slight seasonal variation was found for beta.SNE, but beta.SNE value was largely lower than beta.SIM (see Supplementary Table 3).
Beta diversity decomposition can help clarify processes associated with phytoplankton community construction, which typically decompose into species replacement (or turnover) and richness differences or nestedness (i.e., species gain or loss) (Harrison et al., 1992;Williams, 1996). The relative importance of these two components may help us better understand the underlying mechanisms that shape and maintain community diversity distribution patterns (Baselga, 2010;Podani and Schmera, 2011;Fu et al., 2019). In this study, beta diversity patterns in the Jiulong River estuary were mostly dominated by turnover (spring: 90%, summer: 90%, autumn: 88%), which appears to be a common phenomenon across a wide range of taxa in estuary or marine ecosystems, such as macrobenthos (Alsaffar et al., 2017), aquatic plant communities (Bertuzzi et al., 2018), and marine fish (Liggins et al., 2015).
When alpha and beta diversity was combined, although species richness increased along with gradients (i.e., from the upper to lower reaches of the estuary), phytoplankton composition became more homogeneous as environmental differences decreased (see Supplementary Table 4). With-in and between-site variation are both important to comprehensively FIGURE 8 | Relative contribution of environmental variables and spatial distances to phytoplankton diversity of the overall. explore biodiversity patterns. Since the different scales and the different levels of biodiversity make it relatively complex, a single index cannot reasonably be used to explore diversity patterns. Although a variety of diversity indices have been developed to express various aspects of biodiversity, there is no consensus on which diversity index is best (Socolar et al., 2016). Different indices will yield different diversity patterns, showing different aspects of the same natural phenomenon. In Jiulong River estuary, our study indicated that the species richness of phytoplankton was relatively higher in the lower estuary, while phytoplankton was more diverse and heterogeneous in species composition in the upper estuary, informing flexible conservation strategies in line with different objectives. Thus, it is better to explore diversity patterns hierarchically, incorporating both alpha and beta diversity with various other indices.
Combined with the key factors that influenced alpha and beta diversity, our study found that environmental factors significantly influenced alpha and beta diversity. Moreover, beta diversity was also influenced by spatial distances. Although many studies have shown a strong relationship between spatial distance and beta diversity respective to different taxa in estuarine systems, such as phytoplankton (Josefson and Göke, 2013), benthic invertebrates (Josefson and Göke, 2013), and bivalve mollusks FIGURE 9 | Relative contributions of nutrients, nutrient ratios and spatial distances in beta diversity. (Chust et al., 2013), most of these studies focused on passive dispersal or less mobile species, and the relationship between distance and highly mobile species in estuarine systems is not yet clear. In our study, compared with environmental variables, spatial distance (geographic and hydrologic distance) showed a lower contribution, for example, only about 12% contribution of spatial distances to the beta.SOR (Figure 9), which might be mainly because of the limited spatial scale in Jiulong River estuary. Thus, it was worthy to examine the contribution of spatial distance to phytoplankton beta diversity in a larger scale in further studies.

Roles of Nutrients and Nutrient Ratios in Shaping Phytoplankton Diversity Patterns
Our study found high correlations between phytoplankton diversity indices and environmental variables of nutrients and nutrient ratios in the Jiulong River estuary, indicating that nutrients and nutrient ratios contributed significantly to phytoplankton biodiversity patterns, which is in agreement with findings from previous studies conducted in similar estuarine environments (Schoor et al., 2008;Li et al., 2010;Huang et al., 2013;Chu et al., 2014;Zhang et al., 2015;Anselme et al., 2018). The relative contribution was quantified by calculating the percentage of the absolute coefficient value of all relevant variables that yielded a statistically significant value (p < 0.001) from optimal regression equations (see Figures 9, 10). Results showed that nutrients contributed 25,52,51,40,38, and 88% to S, J, H', beta.SIM, beta.SNE, and beta.SOR, respectively, while nutrients ratio contributed 17, 12, and 25% to S, beta.SNE, and beta.SIM, respectively (Figures 9, 10). In general, nutrients and nutrient ratios contributed significantly to beta diversity and overwhelmingly to alpha diversity, where nutrient ratios influenced S, beta.SIM and beta.SNE and nutrients significantly influenced beta.SOR.
The results indicated that key variables with greatest contribution to beta diversity varied among different seasons, i.e., nutrients in spring, spatial distances in summer and nutrient ratios in autumn, respectively (see Supplementary Table 4). The seasonal difference above might due to the variation of specific variable, such greatest variation of nutrients (TN and TP) found in spring, nutrient ratios (N:P and Si:P) in autumn; by contract, the variation was slight for environmental variables in summer, so spatial distance played a more important role in shaping communities differences between sites (Figure 9).
Previous studies conducted on nutrient ratios primarily focused on N:P and Si:P (Schoor et al., 2008;Li et al., 2010;Huang et al., 2013;Chu et al., 2014;Zhang et al., 2015), while relatively few focused on Si:N (Li et al., 2010;Chu et al., 2014). Justić et al. (1995) proposed that P limitation could potentially occur under conditions where Si:P > 22 and N:P > 22, N limitation could potentially occur under conditions where N:P < 10 and N:Si < 1, and Si limitation could potentially occur under conditions where Si:P < 10 and N:Si > 1. In our study, P limitation could potentially occur under conditions of high Si:P and N:P in the Jiulong River estuary, where the former ranged from 19 to 933 (118 average) and the latter ranged from 16 to 322 (51 average) (Figure 11). Additionally, the fact that both N:P and Si:P exceeded a value of 10 at all sampling sites (Figure 11), the probability of Si and N limitation was low. Our study found that P was a potential limiting factor in the Jiulong River estuary, while the estuary was clearly N-enriched, Si-enriched, and P-limited.
The same nutrient ratio phenomenon was observed in the Yangtze River and Pearl River estuaries, being the two largest and most densely populated and economically developed estuaries in China, where Si:P and N:P values exceeded the recommended threshold (22) (Yin et al., 2001;Li et al., 2010), indicating that similar to the Jiulong River estuary, they were also N-enriched, Si-enriched, and P-limited. However, compared to these two large estuaries, Si:P, N:P, and N:Si values in the Jiulong River estuary were all lower (Figure 12), indicating that the nutrient condition of the Jiulong River estuary was more balanced than these estuaries, even with respect to the general condition along the Chinese coast (Figure 12).
Given that different phytoplankton species differ in their nutrient demands, nutrient concentrations and nutrient ratios may affect phytoplankton diversity (Redfield, 1934;Paerl and Justić, 2013). Additionally, N and P are essential nutrients for all phytoplankton species (Paerl and Justić, 2013), while Si is mostly specific to diatoms (Officer and Ryther, 1980). Taking this into account, we inferred that N played a more important role in phytoplankton diversity patterns than P and Si in the Jiulong River estuary. This is because of a nutrient concentration imbalance in different elements characterized by enriched N, limited P, and specifically Si. We also found that diatoms became more prevalent in the Jiulong River estuary as Si concentrations increased, where diatoms accounted for large proportion of species composition of the phytoplankton community (i.e., 90% in summer and 87% in autumn). Therefore, nutrient concentrations and nutrient ratios may influence phytoplankton diversity patterns. Accordingly, more attention should be paid to trade-offs among nutrient elements instead of just solely focusing on specific nutrient concentrations.

Recommended Jiulong River Estuary Ecosystem Management Practices
Increasingly, beta diversity has been used to guide biodiversity conservation, such as the design of protected areas, the identification of priority areas, and the establishment of ecological corridors (Cody, 1986;Thomas, 1990;Harrison et al., 1992;Pimm and Gittleman, 1992;Nekola and White, 1999;Myers et al., 2000;Wu et al., 2010;Fattorini, 2011). The two beta diversity components reflect two different FIGURE 11 | The range and average of nutrient ratios. Yellow line represented the range values of nutrient ratio in Jiulong River Estuary, and yellow square represented their average. Red dashed line represented the threshold of nutrient ratios to inform potential P, N, and Si limitation (Justić et al., 1995).
FIGURE 12 | The average of nutrient ratios in Yangtze River Estuary (Li et al., 2010), Pearl River Estuary (Yin et al., 2001), Jiulong River Estuary and Chinese coastal waters . Red dashed line represented the threshold of nutrient ratio to inform potential P, N, and Si limitation (Justić et al., 1995). or even opposite processes in community construction (Gaston and Blackburn, 2000;Baselga, 2012); thus, beta diversity decomposition should be incorporated into biodiversity conservation rather than accounting for beta diversity alone (Angeler, 2013). Baselga (2010) suggested that when turnover dominates, the entire area should be protected, while when nestedness dominates, only sites or areas with the highest species richness should be protected to conserve the entire area. On this account, the entire Jiulong River estuary should be protected due to its dominant turnover pattern (accounting for 85%). Additionally, when using alpha and beta diversity, more attention should be paid to conserving diverse species composition of phytoplankton in the upper reaches, and protecting its species abundance in the lower reaches.
This study showed that the Jiulong River estuary is N-enriched, Si-enriched, and P-limited, where an excess in N and Si may lead to eutrophication and impact phytoplankton diversity. Compared to the Yangtze River and Pearl River estuaries, the nutrient status in the Jiulong River estuary was more balanced. Phytoplankton take up extraneous nutrients at a fixed ratio (Redfield ratio: C:N:P = 106:16:1; Redfield, 1934) to meet their own growth and survival requirements, and diatoms specifically require Si in addition to C, N, and P, absorbing Si and N on a 1:1 M basis (Officer and Ryther, 1980). Although there is a strong consensus that N is the primary cause of eutrophication in aquatic marine ecosystems (Paerl, 2009;Cosme et al., 2015), the relatively high Si (compared to N) in the Jiulong River estuary poses a potential algal bloom threat. With an increase in Si concentrations in the Jiulong River estuary, the composition of phytoplankton species richness gradually led to the dominance of diatoms. For example, diatoms accounted for roughly 90% of phytoplankton species richness in our study, which may impact phytoplankton diversity. Therefore, instead of solely focusing on reducing N pollution, it is recommended to also reduce Si concentrations during ecosystem restoration initiatives. This will help restructure the phytoplankton community and prevent algal blooms (via diatoms).

CONCLUSION
Phytoplankton is one of the most susceptible aquatic organisms to environmental change in estuaries, especially regarding eutrophication. In this study, phytoplankton diversity patterns and their potential factors (particularly nutrient concentrations and nutrient ratios) were examined, using both alpha and beta diversity in the Jiulong River estuary, China. Our study indicated that diversity patterns were shown to vary at different scales, in different seasons and in different indices in Jiulong River estuary, where the phytoplankton species richness was higher in autumn and in the lower estuary, while phytoplankton was more diverse and heterogeneous in its composition in the upper estuary, highlighting that it is better to explore diversity patterns hierarchically to inform an effective conservation strategy. Environmental factors had a significant influence on both alpha and beta phytoplankton diversity in Jiulong River estuary, particularly nutrient concentrations and nutrient ratios, characterized by an excess in N and Si but limited P. Our findings provide an important information to enhance biodiversity conservation in Jiulong River estuary in a comprehensive way, and also offer a practical approach and analytically based information to protect estuarine or even other marine ecosystem diversity, which can be applied to other estuaries to better guide sustainable management and conservation practices.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
FG: perform the research. FG, WY, and WH: methodology, data wrangling, and project administration. FG, YW, XL, SA, DZ, WZ, and FK: data analyze. XY, ZM, and ZL: resources and data collection. WY, WH, ZM, and BC: supervision. BC, XY, and WY: funding acquisition. All authors contributed to the conceptualization, writing, review, and editing of the manuscript, and read and agreed to the published version of the manuscript.