Bacterioplankton community responses and the potential ecological thresholds along disturbance gradients

Science of the Total Environment 696 (2019) 134015

Contents lists available at ScienceDirect

Science of the Total Environment journal homepage: www.elsevier.com/locate/scitotenv

Bacterioplankton community responses and the potential ecological thresholds along disturbance gradients Lixia Xuan a,b, Zheliang Sheng b, Jiaqi Lu a,b, Qiongfen Qiu b, Jiong Chen a,b, Jinbo Xiong a,b,⁎ a b

State Key Laboratory for Quality and Safety of Agro-products, Ningbo University, Ningbo 315211, China School of Marine Sciences, Ningbo University, Ningbo 315211, China

H I G H L I G H T S

G R A P H I C A L

A B S T R A C T

• Roles of local variables in governing bacterioplankton assembly are quantified. • Nonlinear responses of bacterioplankton to disturbance intensity and frequency. • The significant tipping points afford a warning line for coastal management. • A few pollution-discriminatory taxa accurately diagnose coastal pollution level.

a r t i c l e

i n f o

Article history: Received 8 June 2019 Received in revised form 24 July 2019 Accepted 19 August 2019 Available online 20 August 2019 Editor: Julian Blasco Keywords: Coastal pollution Bacterioplankton community Non-linear responses Tipping points Distance-decay OPI-discriminatory taxa

a b s t r a c t Increasing intensity and frequency of coastal pollutions are the trajectory to be expected due to anthropogenic pressures. However, it is still unclear how and to what extent bacterioplankton communities respond to the two factors, despite the functional importance of bacterioplankton in biogeochemical cycles. In this study, significant organic pollution index (OPI) and offshore distance gradients, as respective proxies of disturbance intensity and disturbance frequency, were detected in a regional scale across the East China Sea. A multiple regression on matrices (MRM) revealed that the biogeography of bacterioplankton community depended on spatial scale, which was governed by local characters. Bacterioplankton community compositions (BCCs) were primarily governed by the conjointly direct (−0.28) and indirect (−0.48) effects of OPI, while offshore distance contributed a large indirectly effect (0.52). A SEGMENTED analysis depicted non-linear responses of BCCs to increasing disturbance intensity and disturbance frequency, as evidenced by significant tipping points. This was also true for the dominant bacterial phyla. Notably, we screened 30 OPI-discriminatory taxa that could quantitatively diagnose coastal OPI levels, with an overall 79.3% accuracy. Collectively, the buffer capacity of bacterioplankton communities to increasing disturbance intensity and disturbance frequency is limited, of which the significant tipping points afford a warning line for coastal management. In addition, coastal pollution level can be accurately diagnosed by a few OPI-discriminatory taxa. © 2019 Elsevier B.V. All rights reserved.

Abbreviations: OPI, organic pollution index; BCCs, bacterioplankton community compositions; DDR, distance-decay relationship; CAP, constrained analysis of principal coordinates; MRM, multiple regression on matrices; PLS-PM, partial least squares path modeling; ANOSIM, analysis of similarity; perMANOVA, permutational multivariate analysis of variance; GoF, Goodness-of-Fit. ⁎ Corresponding author at: State Key Laboratory for Quality and Safety of Agro-products, Ningbo University, Ningbo 315211, China. E-mail address: [email protected] (J. Xiong).

https://doi.org/10.1016/j.scitotenv.2019.134015 0048-9697/© 2019 Elsevier B.V. All rights reserved.

2

L. Xuan et al. / Science of the Total Environment 696 (2019) 134015

1. Introduction Coastal zones are among the most productive and diverse ecosystems globally, where are suffering a continuous supply of contaminants from rapid urbanization, industrial effluents and agricultural runoff (Nogales et al., 2011; Qiu et al., 2019). Microorganisms contribute indispensible roles in mediating costal biogeochemical recycles, thus being at the base of coast trophic networks (Cao et al., 2017; Isabwe et al., 2018). However, it has been shown that bacterioplankton community is resistant (Allison and Martiny, 2008), resilient (Shade et al., 2011) or sensitive (Xiong et al., 2015) to coastal pollutions. These inconsistencies could result from conditional dependency, disturbance intensity and/or disturbance frequency (Berga et al., 2012; Shade et al., 2012; Xiong et al., 2014; Dai et al., 2017), which in turn make an uncertainty of the coastal ecosystem stability and services. Indeed, it has been proposed that incorporating microbial responses into ecosystem models could improve the predictions to a changing environment (Zhou et al., 2012). Under these scenarios, it is a fundamental prerequisite to understanding how bacterioplankton communities respond to increasing disturbance intensity and disturbance frequency, a trajectory to be expected due to anthropogenic activities (Hu et al., 2017a; Valerie et al., 2018). Bacterioplankton appears to disperse globally and rapidly in water currents (Gibbons et al., 2013). However, increasing evidence has depicted a divergence in microbial community similarity with an increasing geographic distance (Wang et al., 2015; Dai et al., 2017; Meyer et al., 2018). It is now become clear that the distance-decay relationship (DDR) of bacterioplankton community compositions (BCCs) reflects the combined effects of both past events (i.e., spatial distance imposed dispersal limitation, physical barrier) and contemporary environmental variations (i.e., local environmental conditions) (Stegen et al., 2013; Goldmann et al., 2016; Wang et al., 2019). Nevertheless, it is still controversial on the relative importance of geographic distance in the biogeography of BCCs, with reports of negligible (Gibbons et al., 2013), moderate (Wang et al., 2017) or strong role (Dai et al., 2017). These divergences could be attributed to the differences in spatial scale (Martiny et al., 2011), taxonomic resolution (Storch and Šizling, 2008). However, these assertions are rejected by a recent metaanalysis, showing that adjusting for differences in spatial scale, altering taxonomic resolution, or excluding inactive species, insufficiently affects the turnover rate of microbial biogeography (Meyer et al., 2018). Given the functional importance of bacterioplankton, exploring the biogeographic pattern of BCCs could predict the function potentials in ocean (Larsen et al., 2015), which in turn facilitates designing effective costal management and conservation strategies. There is evidence that the degree of changes in BCCs mirrors the intensity of pollution, which is dose-dependent in both experimental and field data (Ager et al., 2010; Washburn et al., 2016; Dai et al., 2017). However, when resource spectrum is sufficiently extreme, a critical transition–a tipping point–occur, which imposes the divergence of BCCs in a way that outlast the trajectory (Moore, 2018; Ratajczak et al., 2018). In other words, disturbance regimes generate predictable but non-linear pattern in BCCs. Consistently, intermediate disturbance hypothesis posits that biodiversity always exhibits an inverted V shape along pollutant gradients (Yuan et al., 2016; Guitet et al., 2018). The tipping point, in turn, affords a descriptive statistic for evaluating the buffer capacity of bacterioplankton communities to an increasing coastal eutrophication. For this reason, understanding when an ecological regime shifts from one state to an alternative state once breached a disruptive threshold (that is, tipping point) is therefore of considerable interest, since it will afford a biological integrity for coastal regulatory authorities. Traditional geochemical analysis has provided important information on evaluating water pollution status, i.e., water quality index (Tyagi et al., 2013), while it does not afford a viewpoint upon biological exposure. By contrast, microbial community could systematically

indicate and record both biotic and abiotic (including unmeasured and unobserved factors) pressures, even after the contaminants have been fully degraded (Smith et al., 2015). In parallel, microorganisms are highly susceptible to environmental disturbances, which can either augment or buffer environmental change (Vanacker et al., 2016). Accordingly, increasing studies have identified microbial indicators, whose abundances are significantly associated with the levels of local pollution. These findings provide a direct clue on the application of bioindicators for diagnosing coastal pollution levels. Indeed, relevant studies have produced long lists of indicative taxa (Borja, 2018), or just simply discriminate the disturbed areas from reference ones (Fodelianakis et al., 2014; Nemati et al., 2018; Xiong et al., 2015). It should be noted that natural and anthropogenic pollutions are generally continuous gradients, rather than discrete, thereby hampering the application of these indicators in practical. In this context, current assessment methods present an important gap regarding the quantitatively predicting pollution levels in situ. It is now recognized that both disturbance intensity and disturbance frequency govern the BCCs (Shade et al., 2012; Danneyrolles et al., 2019), while there is little discussions on how bacterioplankton communities respond to the two counterparts. For this reason, here we recruited organic pollution index (OPI, an index is more comprehensive than a single water variable and clearly reflects water quality) and offshore distance as the proxies of disturbance intensity and disturbance frequency, respectively. The underlying hypotheses are that pollution severity reflects disturbance intensity, while remote sites are less frequently affected by anthropogenic activities. It should be noted that the two indices are complementary, rather than exclusively. To this end, we chose data sets where alone focused on environmental factors that shaping the regional biogeography of BCCs in the East China Sea (Wang et al., 2015). We integrated ecological approaches to explore the responses of BCCs to disturbances with the following concerns: i. Which factors drive the DDR of BCCs? It has been exemplified that costal BCCs and their metabolic potentials can be predicted by a few environmental parameters (Larsen et al., 2015). Thus, information on the factors governing BCCs would add the prediction of coastal services to pollution. To determine this, we used a multiple regression on matrices (MRM) to identify the relative importances of contributing factors for DDR. Further, we employed a partial least squares path modeling (PLS-PM) to evaluate the interactive effects of disturbance intensity and disturbance frequency on BCCs. ii. Are there potentially ecological threshold-inducing shifts in triggering non-liner responses of BCCs along environmental gradients? It is known that the buffer capacity of BCCs to costal pollution is limited, rather than resistance (Dai et al., 2017). Thus, we applied SEGMENTED model to identify the tipping points that cause collapse divergences in BCCs. iii. Whether bio-indicators quantitatively predict coastal pollution? There is evidence that the abundances of some sensitive taxa are closely associated with the levels of organic pollution (Washburn et al., 2016). In this case, we used a randomForest model to identify OPI-discriminatory taxa, which subsequently facilitate independent variables for diagnosing coastal pollution.

2. Materials and methods 2.1. Experimental design and sample collection Surface seawater samples were collected from the East China Sea during a routine-monitoring project of Marine Environmental Monitoring Center of Ningbo (Fig. S1). To negate the effects of seasonality and tide on BCCs, the samples were collected within two weeks in summer (15–28 August 2013) during low tide period. Seawater samples were taken with consistent sampling producers, including 0.5 m depth,

L. Xuan et al. / Science of the Total Environment 696 (2019) 134015

0.5 L water sample filtered onto a 0.2 μm polycarbonate membrane (Millipore, USA), and storage in a tank filled with liquid nitrogen. The sampling region included eight representative locations, where are suffering from excess nutrient pollutants and anthropogenic disturbances, with the exception of Yushan Reserve. Specifically, Yushan Reserve is a marine ecological special protection zone (anthropogenic activities are restricted), where is oligotrophic. Sanmen Bay and Xiangshan Bay are semi-enclosed bays, where are intensively farmed for decades. Hangzhou Bay is located downstream of the Qiantang River and south of the Yangtze River estuary, where is annually subjected to huge volumes of freshwater and sediment discharges. Additionally, anthropogenic inputs have largely increased in the boundary of island-chain and Zhoushan (Ren et al., 2010; Dai et al., 2017). Shipu is a fishing port, while Jiushan is an inshore fishery. The two locations are intensively affected by anthropogenic activities (Fig. S1). In total, 82 samples were enrolled in the analysis, as well as corresponding biogeochemical variables from our previous work (Wang et al., 2015). The disturbance intensities of coastal zones were indexed by OPI CODi DINi DIPi DOi that is defined as ΟΡΙ ¼ þ þ − where COD is chemCODs DINs DIPs DOs ical oxygen demand (mg/L), DIN is dissolved inorganic nitrogen (the sum level of ammonia, nitrate and nitrite, mg/L), DIP is dissolved inorganic phosphorus (mg/L), and DO is dissolved oxygen (mg/L). CODs = 2 mg/L, DINs = 0.2 mg/L, DIPs = 0.015 mg/L and DOs = 6 mg/L are the standard concentrations following the Sea Water Quality Standard of China, GB 3097–1997 (SEPA, 1998). Pollution levels were categorized as follows: water quality is uncontaminated when OPI b 1 and is beginning to be contaminated when 1 ≤ OPI ≤ 2; it is moderately polluted when 2 b OPI ≤ 4, and heavily polluted when OPI N 4. The offshore distances were served as proxy of disturbance frequency, which were calculated by geographic information system (GIS) in ArcGIS. Specifically, offshore distance is calculated by the distances between the sampling site and its nearest shore, with the exception of samples from Hangzhou Bay where the water current is clear. Instead, offshore distances were the distance between sampling site and the mouth of QianTang River (Fig. S1).

2.2. DNA extraction, PCR amplification, sequencing and data processing DNA was extracted using a Power Soil DNA Isolation Kit (Mo Bio, Carlsbad, USA). To minimize reaction-level PCR bias, each sample was amplified in triplicate with 10 ng purified DNA as template for amplifying the V4 region of bacterial 16S rDNA gene. Equimolar amounts of purified amplicons from every sample were pooled and sequenced on a single run using the MiSeq 2 × 250 bp platform (Illumina, San Diego, CA) (Wang et al., 2015). The paired reads were joined with FLASH using default setting (Magoč and Salzberg, 2011). Then, the merged reads were quality filtered and analyzed with Qiime2 (Caporaso et al., 2010). In brief, reads were truncated at any site of more than three sequential bases receiving a Phred quality score (Q) b 20 and any read containing ambiguous base calls was deleted. The remaining sequences were chimera detected and removed using UPARSE (Edgar, 2013). Bacterial sequences were clustered into operational taxonomic units (OTUs) at 97% cutoff. The representative sequences (most abundant read of each OTU) were classified taxonomically against the Greengenes database using PyNAST (De Santis et al., 2006). The OTUs that were not affiliated with bacterial domain, as well as Archaea, Chloroplast, and b0.005% of the total bacterial reads were removed as proposed elsewhere (Bokulich et al., 2013). To correct unequal sequencing efforts and avoid biases associated with repeatedly subsampling the smallest library, the OTU table was randomly rarefied the subset of 11,200 sequences per sample (approximately 80% of the lowest sequencing depth) prior to further analysis. The raw reads were openly accessible under the accession number DRA002865 (Wang et al., 2015).

3

2.3. Statistical analysis The following analyses were performed in R 3.4.1, unless otherwise stated (R Core Team, 2015). A constrained analysis of principal coordinates (CAP) was used to depict the association between BCCs and the selected environmental variables (Legendre and Anderson, 1999). Analysis of similarity (ANOSIM) and permutational multivariate analysis of variance (perMANOVA) were performed to test the significance of BCCs among sampling locations based on Bray–Curtis dissimilarity metrics using ‘vegan’ package (Oksanen et al., 2015). Pairwise geographic distances between sites were calculated based on geographic coordinates. The spatial turnover rate of BCCs was the slope of distance– decay linear regression, which was plotted as a logarithmic community similarity against a logarithmic distance. A MRM model was used to quantify the relative importances of multiple factors in contributing to the DDR using ‘ecodist’ package (Legendre et al., 1994). Partial regression coefficients of an MRM model afforded a measure of the change rate in community similarity for variables of interest when other variables were held constant (Martiny et al., 2011). To quantitatively describe the relationships between environmental variables and bacterial communities, we used a PLS-PM model using ‘plspm’ package (Sanchez, 2013). This approach provides a framework for analyzing multiple relationships between a set of blocks of variables (Sanchez, 2013). We ran PLS-PM to validate the estimates of path coefficients and the coefficients of determination using 1000 bootstraps. Path coefficients (i.e., direct effects) are the strength and direction of the linear relationships between variables, while indirect effects represent the multiplied path coefficients between a predictor and a response variable, adding the product of all possible paths except the direct effect (Barberán et al., 2014). The fit of model was evaluated by the Goodnessof-Fit (GoF) statistic. We implemented the SEGMENTED analysis to identify tipping point of environmental variables as a threshold by integrating linear and generalized linear regression models with the ‘segmented’ package (Muggeo, 2008). The tipping point chosen by SEGMENTED is the mean combined with its standard error of the potential abrupt change, which was observed within the range of a concerned variable (Vanacker et al., 2016). To evaluate if the tipping point is relevant, Davies test was applied to test whether a difference in slope was significant (Davies, 1977). To screen the OPI-discriminatory taxa (OTUs level), relative abundances of OTUs were fitted against their corresponding OPI values using the ‘randomForest’ package (Liaw and Wiener, 2002; Cutler et al., 2007). A smaller model was constructed by including only the most important variables by ranking the order of importance (Strobl et al., 2007). The probability of pollution level was stratified by using the profiles of OPI-discriminatory taxa and corresponding weight coefficients as independent variables. A consistency between observed and predicted OPI category was termed a correct classification, otherwise was termed as false prediction. 3. Results 3.1. Biogeography of BCCs along environmental gradients After quality control, there were a total of 2,617,300 high-quality bacterial sequences (accounting for 97.9% of the raw sequences) and 14,061–56,520 sequences per sample (mean = 31,918 ± 8830, standard deviation), of which 98.4% clean reads was assigned at the phylum level. By contrast, only 77.9% clean reads can be classified at the family level. The predominant bacterioplankton phyla were Bacteroidetes (mean relative abundance, 37.5%) and Proteobacteria (29.6%), followed by Planctomycetes (15.1%), Acidobacteria (5.2%) and Cyanobacteria (3.8%) across the enrolled 82 samples (Fig. S2). In general, relative abundances of the dominant phyla were significantly associated with OPI or offshore distances (Fig. S2). The CAP biplot depicted clear differences in

4

L. Xuan et al. / Science of the Total Environment 696 (2019) 134015

BCCs according to sampling location, OPI level and offshore distance (Fig. 1). Sampling location effects were further supported by ANOSIM (Global r = 0.673, P = 0.001) and perMANOVA (Global F = 12.23, P = 0.001) (Table S1). Strikingly, OPI gradients impose directional changes in BCCs along the first two axes of CAP, which significantly fitted the ordination of bacterioplankton communities (Fig. 1B; Table S2). Consistently, disturbance intensity (OPI) and disturbance frequency (offshore distance) respectively constrained 14.5% (P b 0.001) and 13.3% (P b 0.001) variation in BCCs, as well as their interaction (10.1%, P = 0.001) (Table S3).

3.2. Spatial turnover of bacterioplankton communities Across the study region (ranged from 1.47 to 225.9 km between any two sites), there was a significant and negative (r = −0.609, P b 0.001) DDR, with a turnover rate of −0.070 (Fig. 2). The distance–decay slope within locations was significantly shallower than the overall slope (slope = −0.049, r = −675, P b 0.001), and steeper among locations than the overall slope (slope = −0.072, r = −0.478, P b 0.001) (Fig. 2). Additionally, for every location, the spatial turnover rates were varied markedly, ranged from −0.052 to −0.218, or a random biogeography (Fig. S3). There is no significant correlation between the slope DDR and corresponding averaged OPI (Fig. S3), or alpha diversity (Fig. S4). Notably, BI and JS had similar OPI values (Fig. S5). However, BI exhibited a DDR pattern with a turnover rate of −0.116, while the later was randomly distributed, suggesting that local variables affected the BCCs. A partial Mantel test revealed that environmental distance (ρ = 0.353, P b 0.001) alone affected BCCs when geographical factors were controlled over all spatial scales, and vise versa (ρ = 0.467, P b 0.001) (Table S4). Thus, the both factors independently affected BCCs. Consistently, a forward selection procedure identified that OPI, offshore distance, SP, salinity, pH, and geographic distance were the key variables that were significantly (all P b 0.01) associated with variations in BBCs (data not shown). We further used a MRM model to qualify the relative importances of these contextual environmental variables in governing the DDR pattern (Table 1). Overall, the selected six variables were significant (P = 0.001) with proportion R2 = 0.753 (P = 0.001). For individual variable, seawater pH contributed the largest partial regression coefficient of 0.570–0.952 (P = 0.001), followed by OPI (0.162–0.169, P b 0.01), SP (0.093–0.136, P b 0.01) and geographic distance (0.097–0.113, P b 0.01) (Table 1). Attending to these results, the importances of a given environmental parameter were substantially varied

Fig. 2. Distance–decay curves for the bacterioplankton communities. The solid blue line denotes the regressions across the enrolled spatial scales. The solid black or red lines denote separate regressions within locations or among locations. The turnover rate is estimated using a least-squares linear regression (ln–ln space approach) fit between pair-wise similarities and geographic distances. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

among different spatial scales and sampling locations, which might reflect differences in the underlying mechanisms.

3.3. Relationships between BCCs and environmental variables To better integrate the interrelationships between BCCs and contextual environmental variables, we constructed a PLS-PM model (Fig. 3). The GoF value for the best model represented here was 0.501. OPI exerted the largest contribution (−0.76, imposed by the conjoint direct (−0.28) and indirect (−0.48) effects) on the variation in BCCs. Water OPI was significantly and negatively correlated with salinity (−0.94), which reinforced that anthropogenic discharge and river-borne (i.e., Qiantang River and Yangzi River, also affected by anthropogenic discharge) runoff of freshwater led to eutrophication in the East China Sea. By contrast, pH directly regulated (0.54) the bacterioplankton community. It should be noted that offshore distance inadequately explained the variation in bacterioplankton community (directly effect of −0.06), while it contributed a large indirectly effect (0.52) on BCCs (Fig. 3 and Table S5).

Fig. 1. Constrained analysis of principal coordinates (CAP) of bacterioplankton community structures (BCCs) coded by sampling location (A) and OPI category (B) using Bray–Curtis distance. Arrows indicate direction and magnitude of the variables (OPI and offshore distance) governing BCCs.

L. Xuan et al. / Science of the Total Environment 696 (2019) 134015

5

Table 1 Relative importance of environmental factors contributing to the correlation by multiple regression on matrices (MRM) analysis.

All (R2 = 0.753, P = 0.001) HZ a (R2 = 0.558, P = 0.001) ZS (R2 = 0.841, P = 0.001) SM (R2 = 0.822, P = 0.069) SP (R2 = 0.568, P = 0.219) JS (R2 = 0.612, P = 0.336) XS (R2 = 0.591, P = 0.074) BI (R2 = 0.842, P = 0.001) YS (R2 = 0.105, P = 0.851)

Log (geo)

Log (Od)

Log (OPI + 1)

Log (SP)

Log (Chl a)

pH

Salinity

0.113 (0.001) b 0.002 (0.975) 0.097 (0.003) −0.028 (0.487) 0.178 (0.110) −0.038 (0.327) −0.091 (0.367) 0.045 (0.375) −0.073 (0.374)

−0.022 (0.013) −0.004 (0.956) −0.001 (0.952) −0.048 (0.275) −0.344 (0.046) 0.058 (0.698) −0.229 (0.158) 0.141 (0.394) −0.237 (0.720)

0.162 (0.001) −0.007 (0.977) 0.169 (0.006) 0.069 (0.478) −2.120 (0.724) −0.033 (0.845) −0.183 (0.741) −0.037 (0.674) −0.043 (0.633)

0.113 (0.001) 0.132 (0.004) 0.093 (0.001) 0.191 (0.475) 2.406 (0.049) 0.124 (0.071) 0.058 (0.567) −0.064 (0.595) −0.064 (0.562)

0.003 (0.873) −0.028 (0.383) 0.096 (0.001) 0.351 (0.004) 0.710 (0.602) 0.081 (0.532) 0.355 (0.083) 0.130 (0.004) 0.058 (0.355)

0.570 (0.001) 0.444 (0.155) 0.290 (0.003) 0.393 (0.134) −4.257 (0.094) 0.396 (0.530) 0.573 (0.576) 0.952 (0.001) 1.222 (0.272)

0.003 (0.006) 0.008 (0.007) −0.004 (0.226) 0.094 (0.096) −0.010 (0.803) −0.021 (0.807) 0.190 (0.062) −0.008 (0.685) −0.216 (0.378)

a HZ: Hangzhou Bay, ZS: Zhoushan islands, SM: Sanmen Bay, SP: Shipu, JS: Jiushan, Xiangshan Bay, BI: the boundary of island-chain, YS: Yushan Reserve, Log (geo): Log (geographic distance), Log (Od): Log (Offshore distance); the same as follows. b Values in bracket were P values, of which significant effects were bolder at P b 0.05 level.

3.4. Non-linear responses of bacterioplankton to environmental gradients To exploit the potential ecological thresholds, the responses of bacterial community to OPI, offshore distance, salinity and suspended particle, were evaluated by using the scatter plot and SEGMENTED method. Results showed that OPI, offshore distance and salinity driving the BCCs at threshold levels of 4.61, 46.5 and 33.3‰, respectively. Davies test revealed the significances of the three tipping points (Fig. 4 and Table S6). However, there was no significant tipping point for bacterial diversity in response to increasing OPI and offshore distance (Fig. S6). For the five dominant phyla, their responses to a given environmental gradient were divergent, as well as the significance (Fig. S7). For example, Cyanobacteria and Proteobacteria exhibited a hump-backed trend, while Bacteroides showed an opposite pattern as OPI increased. In addition, there were non-linear and significant responses of Proteobacteria and Planctomycetes phyla along OPI gradient, with tipping points of 4.08 and 14.3 (P b 0.001 in both cases). In contrast, no significant tipping point of OPI was observed for Bacteroidetes and Cyanobacteria (Fig. S7 and Table S6).

were indicative of pollution level. Using the Random Forests machine learning algorithm, the relative abundances of OTUs were regressed against the OPI of each sample. The regression explained 83.8% of the variance related to disturbance intensity. Accordingly, 30 top OPIdiscriminatory taxa were selected according to their feature importances. These taxa were mainly affiliated with Bacteroidetes and Gammaproteobacteria, which divergently responded to the increasing OPI gradient. Notably, relative abundances of overall taxa were significantly (P b 0.05) correlated with OPI levels (Fig. 5A). For this reason, profiles of the 30 OPI-discriminatory taxa were employed to diagnose the OPI level. Notably, the predicted OPI values fitted well (r = 0.975, P b 0.001) with these of the measured OPI in situ (Fig. 5B). In addition, there was a high diagnostic performance, with an overall accuracy of 79.3%. The false diagnoses were mainly attributed to samples within moderate OPI category. However, the prediction accuracies were pretty good in high OPI (90.7% accuracy) and uncontaminated sites (100% accuracy) (Fig. 5C and Table 2). These findings underscore the potential of OPI-discriminatory taxa as biomarkers in diagnosing costal pollution. 4. Discussion

3.5. Organic pollution-discriminatory lineages of coastal water quality

4.1. Factors determine the spatial pattern of BCCs

Given that the BCCs were directionally remodelled along OPI gradient (Fig. 1), an open question is whether certain bacterial lineages

An increasing severity and frequency of coastal pollution is the trajectory to be expected due to anthropogenic activities (Nogales et al., 2011; Xiong et al., 2015). Given that bacterioplankton affords a variety of important ecosystem services, such as pollutant degradation and nutrient cycling, exploring the responses of BCCs could facilitate an initial step to predict the coastal services to a changing environment. As expected, both disturbance intensity and disturbance frequency were the major driving force in governing BCCs (Fig. 1B), which cumulatively contributed 37.9% variation in BCCs (Table S3). The large unexplained portion could be attributed to unmeasured biotic factors, such as protists predation and phytoplankton community that are known to be key variables in governing BBCs (Fortunato and Crump, 2011; Xiong et al., 2016). Notably, OPI was the subset of variables that best explained bacterial community dissimilarities, followed by water salinity and offshore distance (Fig. 3 and Table S5). This makes sense because the three variables define the differences in the environments contributing the divergences of BCCs observed here. OPI and salinity were varied among the estuary, bay and coastal locations, due to different magnitudes of river runoff, tide and anthropogenic discharge. For example, Hangzhou Bay is the estuary of Qiantang River (Fig. S1), where river and ocean communities mix along both salinity and OPI gradients. Ample evidence has shown that salinity is a determinative factor in governing the separation of BCCs (Laque et al., 2010; Xiong et al., 2016). OPI has also been reported to be a driving factor in shaping BCCs (Liu et al., 2011). Although offshore distance was part of the nature in sampling, it can be considered a proxy for many biological and physical factors that vary in multidimension but were not directly measured, including primary production. Consistently, there was a signification association (P b

Fig. 3. Directed graph of the partial least squares path model (PLS-PM). Each box represents an observed variable. The loading for nonmetric multidimensional scaling (NMDS) scores of bacterioplankton communities that create the latent variables are shown in the dashed rectangle. Path coefficients are reflected in the arrow width, with blue and red indicating positive and negative effects, respectively. Dashed arrows show that coefficients insignificantly differ from 0 (P N 0.05). For the best model represented here, the GoF value was 0.501. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

6

L. Xuan et al. / Science of the Total Environment 696 (2019) 134015

Fig. 4. Tipping points for bacterioplankton community responded to key environmental factors. The solid and dashed lines represent the generalized linear model with its confidence interval at 95%. The horizontal solid line represents the standard error of the tipping point.

0.001) between offshore distance and chlorophyll a. In particular, the major effect of offshore distance on BBCs was attributed to the indirect effects (0.52), which reinforce the roles of some unmeasured biotic factors. Consistent with this assertion, the levels of chlorophyll a were significantly (Mantel r = 0.231, P = 0.001) correlated the variations in BCCs. This result is similar to the pattern seen in other coastal zones, showing that bacterial production and chlorophyll a were the top factors in explaining variability of BCCs among locations (Fortunato and Crump, 2011). However, it is worthy to note that disturbance intensity and disturbance frequency in structuring BCCs do not completely erase the effects of other variables. Indeed, there was a clear separation of BCCs according to sampling locations (Fig. 1A and Table S1), indicating the occurrence of local adaptation. This explanation is supported by community data, where locations of SM and XS had similar OPI values, while the BCCs were significantly distinct (Table S1). Thus, bacterioplankton communities within each location could be further partitioned at finer spatial scales (see our discussions below), which are governed by local variables. Although ample evidence has shown DDR of microbial communities, there are very few studies that attempt to disentangle the underlying drivers of this pattern. There was a significant distance–decay slope at the overall scale (pairwise against each other) (Fig. 2), though the relative contribution of geographic distance is relative low (Table 1). This could be due to spatially structured environmental gradients (Meyer et al., 2018; Wang et al., 2019). Consistently, our previous study shows that pure spatial factors solely contribute 9.4% variation in BCCs (Wang et al., 2015). However, a partial Mantel test supported the spatial distance effect when conditioned the role of water variables (r = 0.467, P b 0.001). It has been reported that the biogeography of microbial community is governed by different drivers at different spatial scales (Martiny et al., 2011), which is what we observed here (Fig. 2). The turnover rate at smaller spatial scale (within locations) was lower than that spanning a larger range of spatial scale (among locations) (Fig. 2). One possible explanation for this pattern is that

bacterioplankton communities are dispersal limited, rather than random distribution (Gibbons et al., 2013). Accordingly, more new taxa tend to be detected a larger spatial scale, then a steeper slope should result. Alternatively, bacterioplankton communities exhibit narrower environmental tolerances (Xiong et al., 2016), resulting in heterogeneity of communities under a given set of locally environmental conditions. Across all locations, a significant and large proportion of the variation is exerted by geographic distance and environmental variables (R2 = 0.753, P = 0.001, Table 1). The MRM analysis illustrated that environmental selection exerted a major role in governing those BCCs, showing that water properties except salinity exhibited larger regression coefficients than spatial distance, with the largest contribution from water pH (Table 1). Indeed, water pH has been shown to be the key factor assembling BCCs (Ren et al., 2015). Notably, the slope, significance and driving factor of DDRs were substantially varied among locations (Fig. S3 and Table 2). Possible reasons could include a narrow distance scale embraced, limited sampling sizes for each location, and/ or other unaccounted factors, i.e., biotic variables. For example, at Hangzhou Bay, the spatial distribution of bacterioplankton communities was governed by the gradients of OPI and SP along water flow. In contrast, insignificant DDR in BCCs was observed in the Yushan Reserve, where the bacterioplankton communities were not affected by any measured abiotic factors. A plausible explanation is that the good water quality in this reserved location exerts little selection on bacteria, resulting in stochastic births, death and migration of a given species (Hanson et al., 2012). Additionally, an increasing dispersal rate could counteract compositional differentiation imposed by selection and eliminate the DDR (Zinger et al., 2014), such as Shipu. This location is a fishing port, where is frequently and intensively affected by anthropogenic activities. Therefore, homogenized abiotic factors induced by those activities could contribute more variation in BCCs than spatial distance, which offsets distance effects. Consistently, offshore distance (as proxy of disturbance frequency) imposed a significant effect on the BCCs at Shipu (Table 1). Indeed, it has been shown that the biogeographic patterns

L. Xuan et al. / Science of the Total Environment 696 (2019) 134015

7

Fig. 5. OPI-discriminatory taxa for defining coastal pollution level. A. The top 30 OPI-discriminatory bacterial taxa and their correlation with OPI. Species were shown in descending order of their importances to the accuracy of the model. Importance was determined based on the percentage increase in mean-squared error of microbiota OPI prediction when the relative abundance values of each taxon were randomly permuted. B. The consistency (r = 0.984, P b 0.001) between observed and diagnosed (based on the profiles of OPI-discriminatory taxa) OPI level. C. The predicted probabilities of coastal OPI category based on profiles of the OPI-discriminatory taxa. Inconsistency between predicted and observed category was termed as false diagnose (red cycle) with a cutoff of 50%. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

of microbial communities have differed drivers at different locations (Gao et al., 2019; Martiny et al., 2011). We infer that structuring physical factors in zones, such as geographic distance, disturbance intensity and pH, define the boundaries of regional distinct microbial habitats.

But within each coastal location, variability in microbial communities is explained by local features or unaccounted biotic factors.

Table 2 Predicting the accuracy of pollution stratification by using the abundance of seawater pollution indicator as an independent variable.

A breakpoint occurs where the system responds rapidly to a relatively small change in nutrient enrichment, which indicates the potential ecological processes (Taylor et al., 2014; Lamentowicz et al., 2019). The SEGMENTED analysis depicted significantly non-linear responses of BCCs to increasing disturbance intensity and disturbance frequency (Fig. 4). This pattern reinforces again that the buffer capacity of bacterioplankton communities to disturbance is limited, rather than resistant (Xiong et al., 2015). It has been proposed that tipping points should not be considered as thresholds above or below which we will

True\predicted

OPI N 4

2 b OPI ≤ 4

1 ≤ OPI ≤ 2

OPI b 1

Accuracy (%)

OPI N 4 2 b OPI ≤ 4 1 ≤ OPI ≤ 2 OPI b 1 Overall

39 7 0 0

4 17 2 0

0 1 0 0

0 0 3 9

90.7 68 0 100 79.3

4.2. Non-linear responses of BCCs along gradients of environmental variables

8

L. Xuan et al. / Science of the Total Environment 696 (2019) 134015

have a certain number of species recover, instead as a point from which a regime devastatingly shifts (Vanacker et al., 2016). On this bias, and given the functional importance of bacterioplankton communities in ecosystem services, tipping points could afford a warning line for coastal management. For example, when the OPI level exceeded 4.61 (Fig. 4), it might cause an irreversible disruption of BCCs. This corresponds to the notion that ecosystems with alternate states harbor hysteresis, wherein a return to original pressure fails to recover the transition (Harris, 1999). The response pattern and sensitivity to increasing OPI and disturbance frequency were substantially varied among the dominant phyla. For example, Cyanobacteria and Proteobacteria imposed a humpbacked pattern as OPI increased, with the highest relative abundances detected at intermediate levels of OPI, while Bacteroides exhibited the opposing trend (Fig. 4). Bacteroides populations are capable of particulate organic matter mineralization (Fernándezgómez et al., 2013), which is consistent with their high relative abundance in oligotrophic environments, such as Yushan reserve (Fig. S1). Accordingly, the relative abundance of Bacteroides decreased as OPI increased, with a relative low value of tipping point (Fig. 4). By contrast, fast-growing and copiotrophic Proteobacteria and Planctomycetes populations exhibited significant and higher tipping points of 4.08 and 14.3, respectively. For a given variable, the tipping points and patterns of these dominant phyla were disparate and divergent (Fig. 4), which indicate different sensitivity and ecology of bacterioplankton lineages in response to disturbance. Thus, additional works are required to explore the functional changes in BCCS, which would provide further insights into ecological consequences of coastal pollution.

4.3. OPI-discriminatory taxa are indicative of coastal water quality Microbial communities are generally the first sentinels to external perturbation (Hu et al., 2017b; Yang et al., 2019). As have stated above, BCCs and the dominant phyla exhibited non-linear responses to increasing OPI (Fig. 4), which arise the question whether sensitive bio-indicators could predict the degrees of pollution. To minimize the list of indicators, we screened the top 30 OPI-discriminatory taxa according to their feature importances. These taxa mostly belonged to Bacteroides (i.e., Flavobacteriaceae) and Gammaproteobacteria (were diverse at lower taxonomic level). Despite non-linear responses of BCCs to increasing disturbance intensity, relative abundances of all the 30 OPI-discriminatory taxa were significantly associated with OPI (Fig. 5A). Thus, the OPI-discriminatory taxa exhibited adaptive shifts to increasing pollution. Intriguingly, the OPI-discriminatory Flavobacteriaceae species were consistently and negatively associated with OPI. This pattern is expected, as marine Flavobacteriaceae species have a higher proportion of genes involved in CO2 fixation, as well as preference for consuming polymers (Fernándezgómez et al., 2013). These features in turn facilitate their competitive advantage in nutrient-poor, ocean surface. In contrast, Gammaproteobacteria members divergently responded to increasing OPI, which is concordant with the physiologically diverse characters in this group. For example, clade OM60 was negatively affected by OPI, while Thiohalorhabdaceae and Halomonadaceae species exhibited the opposite trend (Fig. 5A). Consistently, global ocean survey datasets support that OM60 has a clear preference for coastal marine waters (where OPI level is higher), whereas there were fewer from open-ocean surface waters (Yan et al., 2009). On the other hand, Halomonadaceae family consistently enriched at eutrophic location (Xiong et al., 2015). Thus, these OPIdiscriminatory taxa could serve as bio-indicators of water quality and anthropogenic disturbance. It is worthy to note that most of the OPIdiscriminatory taxa could not be classified up to the bacterial genus level (Fig. 5A). This could be attributed to an insufficient coverage of current Greengenes database. Alternatively, there are novel bacterial species in the East China Sea.

Notably, there was a high consistency between predicted and detected level of OPI (Fig. 5B). In addition, profiles of the 30 OPIdiscriminatory taxa can quantitatively and accurately diagnose the OPI categories (an overall 79.3% accuracy, Table 2). The major inconsistencies were found for the samples from beginning to be contaminated sites (Fig. 5C). A possible explanation is the low sampling size (n = 5) in this OPI category. Alternatively, the low pollution is sufficient for bacterioplankton growth, thereby resulting in random birth. That is, the abundance of a given OPI-discriminatory taxon is decoupled from its ecological function (Zhou and Ning, 2017). As a consequence, the profiles of the OPI-discriminatory taxa do not mirror local pollution level. By contrast, the predictive accuracies were 100% for uncontaminated, and 90.7% for heavily polluted sites (Fig. 5C). This makes sense, given that extreme pollutions exert strong filtering on bacterial taxa (Yang et al., 2019), thus resulting in tight associations between the abundances of pollution-discriminatory taxa and the surrounding pollution level. Also, oligotrophic sites may select specialists that are capable of survival in nutrient limited condition, such as high abundance of Flavobacteriaceae species in Yushan Reserve. Collectively, there are several important implications. (i) The buffer capacity of bacterioplankton community is limited. That is, incrementally increasing disturbance intensity and frequency could push beyond a boundary (i.e., tipping point) (Moore, 2018; Ratajczak et al., 2018). These non-linear, abrupt changes afford early warning of regime shifts, which in turn guild coastal management. (ii) Despite non-linear shifts in the whole BCCs, changes in relative abundances of the OPI-discriminatory taxa are gradual. This gradient pattern affords an objective index for diagnosing pollution status by tracking the abundances of a few OPI-discriminatory taxa. However, whether the screened OPI-discriminatory taxa also hold true across coastal zones is not yet known, thus further works are required to validate this pattern. 5. Conclusions The findings obtained here provide comprehensive insights into how increasing disturbance intensity and frequency affect the spatial distribution of bacterioplankton communities. Overall, BCCs are primarily shaped by OPI, offshore distance, suspended particle, pH and salinity. Differences among these factors help to explain the bacterioplankton provincialism. However, at each location, BCCs and DDRs are governed by different local variables (Table 1). Non-linear responses of BCCs to increasing disturbance intensity and frequency reinforces that their buffer capacity is limited. In this regard, the significant tipping points afford a warning line for coastal management. In particular, using the profiles of a few OPI-discriminatory taxa, coastal pollution level can be accurately diagnosed. This is the few attempts to quantitatively diagnose the coastal pollution, though additional survey is required to validate this model. Declaration of competing interest The authors declare no conflict of interest. Acknowledgements This work was supported by the Natural Science Fund for Distinguished Young Scholars of Zhejiang Province (LR19C030001), the National Natural Science Foundation of China (31872693), the Technology Innovation Team of Ningbo (2015C110018), and the K.C. Wong Magna Fund in Ningbo University. Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi. org/10.1016/j.scitotenv.2019.134015.

L. Xuan et al. / Science of the Total Environment 696 (2019) 134015

References Ager, D., Evans, S., Li, H., Lilley, A., Gast, C.v.d., 2010. Anthropogenic disturbance affects the structure of bacterial communities. Environ. Microbiol. 12, 670–678. Allison, S.D., Martiny, J.B., 2008. Resistance, resilience, and redundancy in microbial communities. Proc. Nat. Acad. Sci. U. S. A. 105, 11512–11519. Barberán, A., Ramirez, K.S., Leff, J.W., Bradford, M.A., Wall, D.H., Noah, F., 2014. Why are some microbes more ubiquitous than others? Predicting the habitat breadth of soil bacteria. Ecol. Lett. 17, 794–802. Berga, M., Székely, A., Langenheder, S., 2012. Effects of disturbance intensity and frequency on bacterial community composition and function. PLoS One 7, e36959. Bokulich, N.A., Subramanian, S., Faith, J.J., Gevers, D., Gordon, J.I., Knight, R., Mills, D.A., Caporaso, J.G., 2013. Quality-filtering vastly improves diversity estimates from illumina amplicon sequencing. Nat. Methods 10, 57–59. Borja, Á., 2018. Testing the efficiency of a bacterial community-based index (microgAMBI) to assess distinct impact sources in six locations around the world. Ecol. Indic. 85, 594–602. Cao, X., Wang, J., Liao, J., Gao, Z., Jiang, D., Sun, J., Zhao, L., Huang, Y., Luan, S., 2017. Bacterioplankton community responses to key environmental variables in plateau freshwater lake ecosystems: a structural equation modeling and change point analysis. Sci. Total Environ. 580, 457–467. Caporaso, J.G., Kuczynski, J., Stombaugh, J., Bittinger, K., Bushman, F.D., Costello, E.K., Fierer, N., Peña, A.G., Goodrich, J.K., Gordon, J.I., Huttley, G.A., Kelley, S.T., Knights, D., Koenig, J.E., Ley, R.E., Lozupone, C.A., McDonald, D., Muegge, B.D., Pirrung, M., Reeder, J., Sevinsky, J.R., Turnbaugh, P.J., Walters, W.A., Widmann, J., Yatsunenko, T., Zaneveld, Z., Knight, R., 2010. Qiime allows analysis of high-throughput community sequencing data. Nat. Methods 7, 335–336. Cutler, D., Edwards, T.J., Beard, K., Cutler, A., Hess, K., Gibson, J., Lawler, J., 2007. Random forests for classification in ecology. Ecology 88, 2783–2792. Dai, W., Zhang, J., Tu, Q., Ye, D., Qiu, Q., Xiong, J., 2017. Bacterioplankton assembly and interspecies interaction indicating increasing coastal eutrophication. Chemosphere 177, 317–325. Danneyrolles, V., Dupuis, S., Fortin, G., Leroyer, M., de Römer, A., Terrail, R., Vellend, M., Boucher, Y., Laflamme, J., Bergeron, Y., Arseneault, D., 2019. Stronger influence of anthropogenic disturbance than climate change on century-scale compositional changes in northern forests. Nat. Commun. 10, 1265. Davies, R., 1977. Hypothesis testing when a nuisance parameter is present only under the alternatives. Biometrika 64, 247–254. De Santis, T.Z., Hugenholtz, P., Keller, K., Brodie, E.L., Larsen, N., Piceno, Y.M., Phan, R., Andersen, G.L., 2006. NAST: a multiple sequence alignment server for comparative analysis of 16S rRNA genes. Nucleic Acids Res. 34, 394–399. Edgar, R.C., 2013. UPARSE: highly accurate otu sequences from microbial amplicon reads. Nat. Methods 10, 996–998. Fernándezgómez, B., Richter, M., Schüler, M., Pinhassi, J., Acinas, S.G., González, J.M., Pedrósalió, C., 2013. Ecology of marine Bacteroidetes: a comparative genomics approach. ISME J 7, 1026–1037. Fodelianakis, S., PapageoSrgiou, N., Pitta, P., Kasapidis, P., Karakassis, I., Ladoukakis, E.D., 2014. The pattern of change in the abundances of specific bacterioplankton groups is consistent across different nutrient-enriched habitats in Crete. Appl. Environ. Microbiol. 80, 3784–3792. Fortunato, C.S., Crump, B.C., 2011. Bacterioplankton community variation across river to ocean environmental gradients. Microb. Ecol. 62, 374–382. Gao, Q., Yang, Y., Feng, J., Tian, R., Guo, X., Ning, D., Hale, L., Wang, M., Cheng, J., Wu, L., Zhao, M., Zhao, J., Wu, L., Qin, Y., Qi, Q., Liang, Y., Sun, B., Chu, H., Zhou, J., 2019. The spatial scale dependence of diazotrophic and bacterial community assembly in paddy soil. Glob. Ecol. Biogeogr. 28, 1093–1105. Gibbons, S.M., Caporaso, J.G., Pirrung, M., Field, D., Knight, R., Gilbert, J.A., 2013. Evidence for a persistent microbial seed bank throughout the global ocean. Proc. Natl. Acad. Sci. U. S. A. 110, 4651–4655. Goldmann, K., Schröter, K., Pena, R., Schöning, I., Schrumpf, M., Buscot, F., Polle, A., Wubet, T., 2016. Divergent habitat filtering of root and soil fungal communities in temperate beech forests. Sci. Rep. 6, 31439. Guitet, S., Sabatier, D., Brunaux, O., Couteron, P., Denis, T., Freycon, V., Gonzalez, S., Hérault, B., Jaouen, G., Molino, J.F., Pélissier, R., Richard-Hansen, C., Vincent, G., 2018. Disturbance regimes drive the diversity of regional floristic pools across guianan rainforest landscapes. Sci. Rep. 8, 3872. Hanson, C., Fuhrman, J., Horner-Devine, M.C., Martiny, J.B., 2012. Beyond biogeographic patterns: processes shaping the microbial landscape. Nat. Rev. Microbiol. 10, 497–506. Harris, G.P., 1999. Comparison of the biogeochemistry of lakes and estuaries: ecosystem processes, functional groups, hysteresis effects and interactions between macroand microbiology. Mar. Freshw. Res. 50, 791–811. Hu, A., Ju, F., Hou, L., Li, J., Yang, X., Wang, H., Mulla, S.I., Sun, Q., Bürgmann, H., Yu, C., 2017a. Strong impact of anthropogenic contamination on the co-occurrence patterns of a riverine microbial community. Environ. Microbiol. 19, 4993–5009. Hu, A., Wang, H., Yang, X., Hou, L., Li, J., Li, S., Yu, C., 2017b. Seasonal and spatial variations of prokaryoplankton communities in a salinity-influenced watershed, China. FEMS Microbiol. Ecol. 93 (fix093). Isabwe, A., Yang, J.R., Wang, Y., Liu, L., Chen, H., Yang, J., 2018. Community assembly processes underlying phytoplankton and bacterioplankton across a hydrologic change in a human-impacted river. Sci. Total Environ. 630, 658–667. Lamentowicz, M., Gałka, M., Marcisz, K., Słowiński, M., Kajukało-Drygalska, K., Dayras, M.D., Jassey, V.J., 2019. Unveiling tipping points in long-term ecological records from sphagnum-dominated peatlands. Biol. Lett. 15, 20190043.

9

Laque, T., Farjalla, V.F., Rosado, A.S., Esteves, F.A., 2010. Spatiotemporal variation of bacterial community composition and possible controlling factors in tropical shallow lagoons. Microb. Ecol. 59, 819–829. Larsen, P.E., Scott, N., Post, A.F., Field, D., Knight, R., Hamada, Y., Gilbert, J.A., 2015. Satellite remote sensing data can be used to model marine microbial metabolite turnover. ISME J 9, 166–179. Legendre, P., Anderson, M., 1999. Distance-based redundancy analysis: testing multispecies responses in multifactorial ecological experiments. Ecol. Monogr. 69, 1–24. Legendre, P., Lapointe, F., Casgrain, P., 1994. Modeling brain evolution from behavior: a permutational regression approach. Evolution 48, 1487–1499. Liaw, A., Wiener, M., 2002. Classification and regression by randomForest. R News, pp. 18–22. Liu, S., Lou, S., Kuang, C., Huang, W., Chen, W., Zhang, J., Zhong, G., 2011. Water quality assessment by pollution-index method in the coastal waters of Hebei Province in western Bohai Sea, China. Mar. Pollut. Bulletin. 62, 2220–2229. Magoc, T., Salzberg, S.L., 2011. FLASH: fast length adjustment of short reads to improve genome assemblies. Bioinformatics 27, 2957–2963. Martiny, J.B.H., Eisen, J.A., Kevin, P., Allison, S.D., Horner-Devine, M.C., 2011. Drivers of bacterial beta-diversity depend on spatial scale. Proc. Natl. Acad. Sci. U. S. A. 108, 7850–7854. Meyer, K., Memiaghe, H., Korte, L., Kenfack, D., Alonso, A., Bohannan, B., 2018. Why do microbes exhibit weak biogeographic patterns? ISME J 12, 1404–1413. Moore, J.C., 2018. Predicting tipping points in complex environmental systems. Proc. Natl. Acad. Sci. U. S. A. 115, 635–636. Muggeo, V.M.R., 2008. Segmented: An R package to fit regression models with brokenline relationships. R News. vol. 8, pp. 20–25. Nemati, H., Shokri, M.R., Ramezanpour, Z., Ebrahimi Pour, G.H., Muxika, I., Borja, Á., 2018. Sensitivity of indicators matters when using aggregation methods to assess marine environmental status. Mar. Pollut. Bulletin. 128, 234–239. Nogales, B., Lanfranconi, M.P., Piña-Villalonga, J.M., Bosch, R., 2011. Anthropogenic perturbations in marine microbial communities. FEMS Microbiol. Rev. 35, 275–298. Oksanen, J.F., Blanchet, G., Kindt, R., Legendre, P., Minchin, P.R., O'Hara, R.B., Stevens, M.H.H., Oksanen, M.J., 2015. Vegan: community ecology package. R package version 2.3–0. Qiu, Z., Kennen, J.G., Giri, S., Walter, T., Kang, Y., Zhang, Z., 2019. Reassessing the relationship between landscape alteration and aquatic ecosystem degradation from a hydrologically sensitive area perspective. Sci. Total Environ. 650, 2850–2862. R Core Team, 2015. A Language and Environment for Statistical Computing. The R Foundation for Statistical Computing, Vienna, Austria http://www.R-project.org/. Ratajczak, Z., Carpenter, S.R., Ives, A.R., Kucharik, C.J., Ramiadantsoa, T., Stegner, M.A., Williams, J.W., Zhang, J., Turner, M.G., 2018. Abrupt change in ecological systems: inference and diagnosis. Trends Ecol. Evol. 33, 513–526. Ren, J.L., Zhang, J., Li, D.D., Cheng, Y., Liu, S.M., 2010. Behavior of dissolved inorganic arsenic in the Yellow Sea and East China Sea. Deep-Sea Res. PT II. 57, 1035–1046. Ren, L., Jeppesen, E., He, D., Wang, J., Liboriussen, L., Xing, P., Wu, Q., 2015. pH influences the importance of niche-related and neutral processes in lacustrine bacterioplankton assembly. Appl. Environ. Microbiol. 287, 3104–3114. Sanchez, G., 2013. PLS Path Modeling with R. Trowchez Editions, Berkeley. SEPA (State Environmental Protection and Administration of China), 1998. State Oceanic Administration of China. Sea water quality standard GB 3097-1997 (in Chinese). Shade, A., Read, J.S., Welkie, D.G., Kratz, T.K., Wu, C.H., McMahon, K.D., 2011. Resistance, resilience and recovery: aquatic bacterial dynamics after water column disturbance. Environ. Microbiol. 13, 2752–2767. Shade, A., Peter, H., Allison, S.D., Baho, D.L., Berga, M., Bürgmann, H., Huber, D.H., Langenheder, S., Lennon, J.T., Martiny, J.B.H., Matulich, K.L., Schmidt, T.M., Handelsman, J., 2012. Fundamentals of microbial community resistance and resilience. Front. Microbiol. 3, 417. Smith, M.B., Rocha, A.M., Smillie, C.S., Olesen, S.W., Paradis, C., Wu, L., Campbell, J.H., Fortney, J.L., Mehlhorn, T.L., Lowe, K.A., Earles, J.E., Phillips, J., Techtmann, S.M., Joyner, D.C., Elias, D.A., Bailey, K.L., Hurt Jr., R.A, Preheim, S.P., Sanders, M.C., Yang, J., Mueller, M.A., Brooks, S., Watson, D.B., Zhang, P., He, Z., Dubinsky, E.A., Adams, P.D., Arkin, A.P., Fields, M.W., Zhou, J., Alm, E.J., Hazen, T.C., 2015. Natural bacterial communities serve as quantitative geochemical biosensors. MBio 6 e00326-00315. Stegen, J.C., Lin, X., Fredrickson, J.K., Chen, X., Kennedy, D.W., Murray, C.J., Rockhold, M.L., Konopka, A., 2013. Quantifying community assembly processes and identifying features that impose them. ISME J 7, 2069–2079. Storch, D., Šizling, A.L., 2008. The concept of taxon invariance in ecology: do diversity patterns vary with changes in taxonomic resolution? Folia. Geobotanica. 43, 329–344. Strobl, C., Boulesteix, A.L., Zeileis, A., Hothorn, T., 2007. Bias in random forest variable importance measures. BMC Bioinformatics 8, 1–10. Taylor, J.M., King, R.S., Pease, A.A., Winemiller, K.O., 2014. Nonlinear response of stream ecosystem structure to low-level phosphorus enrichment. Freshw. Biol. 59, 969–984. Tyagi, S., Sharma, B., Singh, P., Dobhal, R., 2013. Water quality assessment in terms of water quality index. Am. J. Water Resour. 1, 34–38. Valerie, D.A., Icoquih, Z.P., Jazmín, B.l., Augusto, C.P.H., Bruno, C.M., Marcos, G.L., Niza, G.T., Maribel, H.R., Luis, E.E., Valeria, S., 2018. Understanding the mechanisms behind the response to environmental perturbation in microbial mats: a metagenomicnetwork based approach. Front. Microbiol. 9, 2606. Vanacker, M., Wezel, A., Arthaud, F., Guérin, M., Robin, J., 2016. Determination of tipping points for aquatic plants and water quality parameters in fish pond systems: a multiyear approach. Ecol. Indic. 64, 39–48. Wang, K., Ye, X., Chen, H., Zhao, Q., Hu, C., He, J., Qian, Y., Xiong, J., Zhu, J., Zhang, D., 2015. Bacterial biogeography in the coastal waters of northern Zhejiang, East China Sea is highly controlled by spatially structured environmental gradients. Environ. Microbiol. 17, 3898–3913.

10

L. Xuan et al. / Science of the Total Environment 696 (2019) 134015

Wang, P., Zhao, J., Xiao, H., Yang, W., Yu, X., 2019. Bacterial community composition shaped by water chemistry and geographic distance in an anthropogenically disturbed river. Sci. Total Environ. 655, 61–69. Wang, X., Lü, X., Yao, J., Wang, Z., Deng, Y., Cheng, W., Zhou, J., Han, X., 2017. Habitatspecific patterns and drivers of bacterial β-diversity in China's drylands. ISME J 11, 1345–1358. Washburn, T., Rhodes, A., Montagna, P., 2016. Benthic taxa as potential indicators of a deep-sea oil spill. Ecol. Indic. 71, 587–597. Xiong, J., Sun, H., Peng, F., Zhang, H., Xue, X., Gibbons, S.M., Gilbert, J.A., Chu, H.Y., 2014. Characterizing changes in soil bacterial community structure in response to shortterm warming. FEMS Microbiol. Ecol. 89, 281–292. Xiong, J., Chen, H., Hu, C., Ye, X., Kong, D., Zhang, D., 2015. Evidence of bacterioplankton community adaptation in response to long-term mariculture disturbance. Sci. Rep. 5, 15274. Xiong, J., Xiong, S., Qian, P., Zhang, D., Liu, L., Fei, Y., 2016. Thermal discharge-created increasing temperatures alter the bacterioplankton composition and functional redundancy. AMB Express 6, 68.

Yan, S., Fuchs, B.M., Lenk, S., Harder, J., Wulf, J., Jiao, N.Z., Arnann, R., 2009. Biogeography and phylogeny of the NOR5/OM60 clade of Gammaproteobacteria. Syst. Appl. Microbiol. 32, 124–139. Yang, Y., Gao, Y., Huang, X., Ni, P., Wu, Y., Deng, Y., Zhan, A., 2019. Adaptive shifts of bacterioplankton communities in response to nitrogen enrichment in a highly polluted river. Environ. Pollut. 245, 290–299. Yuan, Z., Jiao, F., Li, Y., Kallenbach, R.L., 2016. Anthropogenic disturbances are key to maintaining the biodiversity of grasslands. Sci. Rep. 6, 22132. Zhou, J., Ning, D., 2017. Stochastic community assembly: does it matter in microbial ecology? Microbiol. Mol. Biol. Rev. 81, e00002–e00017. Zhou, J., Xue, K., Xie, J., Deng, Y., Luo, Y., 2012. Microbial mediation of carbon-cycle feedbacks to climate warming. Nat. Clim. Chang. 2, 106–110. Zinger, L., Boetius, A., Ramette, A., 2014. Bacterial taxa-area and distance-decay relationships in marine environments. Mol. Ecol. 23, 954–964.