Adaptive shifts of bacterioplankton communities in response to nitrogen enrichment in a highly polluted river

Accepted Manuscript Adaptive shifts of bacterioplankton communities in response to nitrogen enrichment in a highly polluted river Yuzhan Yang, Yangchun Gao, Xuena Huang, Ping Ni, Yueni Wu, Ye Deng, Aibin Zhan PII:

S0269-7491(18)31257-0

DOI:

https://doi.org/10.1016/j.envpol.2018.11.002

Reference:

ENPO 11829

To appear in:

Environmental Pollution

Received Date: 23 March 2018 Revised Date:

25 October 2018

Accepted Date: 1 November 2018

Please cite this article as: Yang, Y., Gao, Y., Huang, X., Ni, P., Wu, Y., Deng, Y., Zhan, A., Adaptive shifts of bacterioplankton communities in response to nitrogen enrichment in a highly polluted river, Environmental Pollution (2018), doi: https://doi.org/10.1016/j.envpol.2018.11.002. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ED M AN

ACCEPTED MANUSCRIPT 1

Adaptive shifts of bacterioplankton communities in response to

2

nitrogen enrichment in a highly polluted river

3

Yuzhan Yanga, Yangchun Gaoa,b, Xuena Huanga,b, Ping Nia,b, Yueni Wua,b, Ye Denga,b,

5

Aibin Zhana,b,*

6

a

7

Shuangqing Road, Haidian District, Beijing 100085, China

8

b

9

Beijing 100049, China

RI PT

4

SC

Research Center for Eco-Environmental Sciences, Chinese Academy of Sciences, 18

M AN U

University of Chinese Academy of Sciences, 19A Yuquan Road, Shijingshan District,

Conflict of interest

11

Declarations of interest: none.

12

Correspondence

13

Dr. Aibin Zhan, Research Center for Eco-Environmental Sciences, Chinese Academy

14

of Sciences, 18 Shuangqing Road, Haidian District, Beijing 100085, China. Tel & Fax:

15

+86-10-62849882. E-mail: [email protected]

EP AC C

16

TE D

10

1

ACCEPTED MANUSCRIPT

ABSTRACT

18

Anthropogenic activity-mediated nutrient pollution, especially nitrogen enrichment,

19

poses one of the major threats to river ecosystems. However, it remains unclear how

20

and to which extent that it affects aquatic microbial communities, especially in

21

heavily polluted rivers. In this study, a significant environmental gradient, especially

22

nitrogen gradient, was observed along a wastewater receiving river, the North Canal

23

River (NCR). The pollution was highest, moderate, and lowest in the up-, middle, and

24

downstream, respectively. The community composition of bacterioplankton

25

transitioned from being Betaproteobacteria-dominated upstream to

26

Gammaproteobacteria-dominated downstream. Copiotrophic groups, such as

27

Polynucleobacter (Betaproteobacteria) and Hydrogenophaga (Betaproteobacteria),

28

were dominant in the upstream. Despite this variation, we identified a core

29

microbiome of 32 operational taxonomic units (OTUs), which were dominated by

30

Betaproteobacteria and Gammaproteobacteria taxa. Multiple statistical analyses

31

indicated that total nitrogen (TN) was the most important factor driving the adaptive

32

shifts of community structure. Analyses of co-occurrence networks showed that the

33

complexity of networks was disrupted in the up- and middle streams, while enhanced

34

in the downstream. Our findings here suggested that microbial interactions were

35

reduced in response to the aggravation of nutrient pollution. Similar to these changes,

36

we observed significant dissimilarity of composition of functional groups, with

37

highest abundance of nitrogen metabolism members under the highest level of

38

nitrogen enrichment. Further analyses indicated that most of these functional groups

AC C

EP

TE D

M AN U

SC

RI PT

17

2

ACCEPTED MANUSCRIPT belonged to Betaproteobacteria, suggesting the potential coupling of community

40

composition and functional diversity. In summary, adaptive shifts of bacterioplankton

41

community composition, as well as altered species interactions, occurred in response

42

to nutrient pollution in highly polluted water bodies.

43

Capsule

44

Bacterioplankton communities showed adaptive shifts in response to nutrient

45

pollution in the highly disturbed river ecosystems.

46

Keywords

47

Microbial community, adaptive shifts, molecular ecological networks, heavy nitrogen

48

pollution, microbial interactions

AC C

EP

TE D

M AN U

SC

RI PT

39

3

ACCEPTED MANUSCRIPT

Introduction

50

In the past half century, anthropogenic activities have markedly altered the

51

biogeochemical recycles of nitrogen, phosphorous and other nutrients in terrestrial,

52

atmospheric, and aquatic ecosystems, leading to severe ecological consequences

53

(Camargo and Alonso, 2006; He et al., 2011). For example, in disturbed aquatic

54

ecosystems, nutrient enrichment has caused significantly diverse ecological effects,

55

including toxic algal blooms, oxygen depletion, and biodiversity loss (Dudgeon et al.,

56

2006; Woodward et al., 2012). Being a crucial component of aquatic ecosystems,

57

aquatic community plays important roles in ecosystem functioning while being

58

directly or indirectly impacted by nutrient pollution. However, it remains largely

59

unexploited on ecological mechanisms that affect the assembly of aquatic

60

communities by rapid changes of complex environmental factors.

TE D

M AN U

SC

RI PT

49

In aquatic ecosystems and their food webs, microbes play important roles in

62

nutrient cycles, energy flow, as well as pollutant degradation (Azam and Malfatti,

63

2007). This significance has attracted tremendous interests to investigate patterns of

64

their composition and distribution, as well as associated mechanisms (Read et al.,

65

2015; Jeffries et al., 2016; Niño-García et al., 2016; Wang et al., 2017). Pioneering

66

results on riverine bacteria have confirmed the non-random distribution of community

67

composition, with inconsistent predictors being identified for different rivers (Jones,

68

2012; Zinger et al., 2012; Zeglin, 2015). In less polluted rivers, water

69

physicochemical parameters are frequently identified, including water temperature

70

(Crump and Hobbie, 2005), water residence time (Read et al., 2015), and transparency

AC C

EP

61

4

ACCEPTED MANUSCRIPT (Liu et al., 2013). While in polluted rivers, nutrients (nitrogen, phosphorous and

72

others) and hypoxia are usually identified as most important influential factors

73

(Ibekwe et al., 2016; Jeffries et al., 2016; Ma et al., 2016; Wang et al., 2017). In

74

addition, copiotrophic groups are often abundant in nutrient-enriched river segments

75

(Jeffries et al., 2016; Ma et al., 2016). Therefore spatial distribution and diversity of

76

bacterial communities could reflect their ecological functions, as well as riverine

77

ecosystem functioning. However, current conclusions are inconsistent on the factors

78

responsible for community assembly and controversial on the relationship between

79

community structure and functioning (Zinger et al., 2012; Freimann et al., 2013;

80

Zeglin, 2015). Besides, abiotic factors alone are frequently incomplete to predict the

81

dynamics of community assembly, implying the necessity to incorporate biotic factors,

82

such as species-species interactions (Louca et al., 2016). Thus, further investigations

83

on these themes could improve our understanding of the response and adaptation of

84

microbes to nutrient enrichment.

SC

M AN U

TE D

Microbial community, which is often described as a black box, is not just the

EP

85

RI PT

71

physical assembly of diverse species (Faust and Raes, 2012; Cardona et al., 2016). In

87

a given habitat, diverse microbial assemblages interact and form complex networks

88

through various types of interactions, such as competition and/or mutualism (Faust

89

and Raes, 2012; Ghoul and Mitri, 2016). Consequently, microbial interactions can

90

influence the performance and functioning of bacterial communities, even the

91

functioning of an entire ecosystem (Deng et al., 2016; Morriën et al., 2017; Li et al.,

92

2017). For instance, laboratory experiments have successfully proved that microbial

AC C

86

5

ACCEPTED MANUSCRIPT interactions contributed to rapid successions of microbial communities (Datta et al.,

94

2016). The complexity of species interactions, such as network size and average

95

clustering coefficient, is supposed to provide useful indices for evaluating the stability

96

of microbial communities (Faust and Raes, 2012). Stressful state such as lower

97

precipitation led to unstable and relatively simple interaction networks in desert

98

ecosystem (Wang et al., 2018). Similarly, heavy eutrophication of coastal ecosystems

99

witnessed disrupted inter-specific interactions (Dai et al., 2017). However, to the best

SC

RI PT

93

of our knowledge, few studies analyzed to what extent that the complexity of

101

microbial interactions is altered under nutrient enrichment in riverine ecosystems.

102

Although analysis of microbial interactions is still at its infancy, it could provide

103

alternative investigations on the response and adaptation of bacterial communities to

104

nutrient pollution.

TE D

105

M AN U

100

The North Canal River (NCR), which runs through megacities such as Beijing and Tianjin, provides a good model system to investigate adaptive shifts of microbial

107

communities to nutrient enrichment. Diverse organic and inorganic pollutants have

108

been released into NCR, including pharmaceuticals and personal care products

109

(PPCPs), nutrient pollutants, heavy metals, and pesticides (Heeb et al., 2012). Nutrient

110

pollution, especially the enrichment of nitrogen and phosphorous, has intensified

111

massively in the past four decades (Pernet-Coudrier et al., 2012; Shan et al., 2012).

112

According to previous research, water quality of NCR displays significant spatial and

113

seasonal heterogeneity (Ma et al., 2016; Yang et al., 2018), which is largely owing to

114

different sources of pollution. The upstream of NCR is mainly influenced by effluents

AC C

EP

106

6

ACCEPTED MANUSCRIPT of wastewater treatment plants (WWTPs) scattering in Beijing (Guo et al., 2012). The

116

middle stream is dominated by agricultural areas, with 90% of irrigation water

117

deriving from the discharge of upstream (Pernet-Coudrier et al., 2012). Water quality

118

is therefore impacted by both upstream residuals and agricultural diffuse pollution.

119

The downstream mainly receives discharges from upstream rivers and industrial

120

wastewater from Tianjin (Heeb et al., 2012). The potential environmental gradient can

121

exert selective pressures on microbial communities. For instance, pioneering studies

122

have found that spatial distribution and community composition of bacterioplankton

123

in the upstream of NCR was mainly influenced by total phosphorous (Yu et al., 2012),

124

while by total nitrogen in the downstream (Ma et al., 2016). However, none of these

125

studies about NCR investigated community-level functional changes of microbial

126

communities under severe pollution, nor did they provide investigations on changes of

127

microbial interactions caused by nutrient pollution.

SC

M AN U

TE D

128

RI PT

115

In this study, we sampled bacterioplankton from NCR to investigate the response of microbial communities and inter-species interactions to heavy nutrient pollution.

130

With a significant environmental gradient (mainly the nitrogen gradient), we aim to

131

address: (i) the influence of eutrophication (especially nitrogen enrichment) on the

132

spatial patterns of bacterioplankton communities, (ii) the alternation of bacterial

133

functions in response to environmental changes, (iii) the relationship between the

134

strength of microbial interactions (i.e., the complexity of networks) and eutrophication

135

(especially nitrogen enrichment).

AC C

EP

129

136 7

ACCEPTED MANUSCRIPT

Materials and methods

138

Sample collection and physicochemical measurements

139

Water samples were collected from 39 sites of NCR (Figure 1) in March 2017. Water

140

samples (0.5 L) for bacterioplankton analysis were filtered through 0.22 µm mixed

141

cellulose ester membranes (Millipore, Bedford, MA) to collect microbial cells.

142

Membranes were stored at -80 ºC until further treatments. Coordinates of each site

143

were recorded using a Garmin Handheld GPS navigator.

SC

RI PT

137

For each site, water physical parameters including temperature (T), electric

145

conductivity (EC), pH, oxidation-reduction potential (ORP) and total dissolved solid

146

(TDS) were measured using multi-parameter water quality sonde (MYRON Company,

147

USA). Dissolved oxygen (DO) and concentration of chlorophyll a (Chl_a) were

148

measured using portable dissolved oxygen meter (HACH Company, USA) and

149

Handheld Fluorometer (Turner Designs, USA), respectively. The concentration of

150

total nitrogen (TN), nitrate (NO3-N), ammonia (NH4-N), total phosphorous (TP),

151

chemical oxygen demand (COD), and soluble reactive phosphorous (SRP) was

152

measured in the laboratory using methods by Xiong et al. (2016). The concentration

153

of total organic carbon (TOC) was measured using Shimadzu-TC (Shimadzu, Japan).

154

We also measured the concentration of potassium (K), calcium (Ca), sodium (Na), and

155

magnesium (Mg) using inductively coupled plasma optical emission spectrometry,

156

and the concentrations of chromium (Cr), iron (Fe), manganese (Mn), nickel (Ni),

157

copper (Cu), zinc (Zn), arsenic (As), and lead (Pb) using inductively coupled plasma

158

mass spectrometry (see all environmental variables in supplementary Figure S1).

AC C

EP

TE D

M AN U

144

8

ACCEPTED MANUSCRIPT DNA extraction, PCR amplification, and high-throughput sequencing

160

Total genomic DNA was extracted from membranes using the CTAB extraction

161

protocol modified from Huang et al. (2009). DNA extracts were used as templates to

162

amplify the hypervariable V4 region of 16S rRNA with the primer pair of 515F

163

(GTGCCAGCMGCCGCGGTAA) and 806R (GGACTACHVGGGTWTCTAAT)

164

(Caporaso et al., 2011). The primers for each sample were modified with an addition

165

of 12 nucleotide tag at the 5’-end to distinguish all samples. A total of three replicates

166

were performed for each sample to avoid biased amplification (Zhan et al., 2014;

167

Yang et al., 2016). The amplification process was 4 min of denaturation at 94 ºC, 35

168

cycles of 45 s at 94 ºC, 60 s at 50 ºC, and 90 s at 72 ºC, followed by a final elongation

169

for 10 min at 72 ºC. Three replicates of each sample were then pooled together and

170

purified using the Sangon column PCR product purification kit (Sangon Biotech,

171

Shanghai, China). The high-throughput sequencing was conducted using the Illumina

172

Miseq PE250 Platform.

173

Data processing and functional annotations

174

Raw sequence data was de-multiplexed based on unique tags. Paired-end reads were

175

merged before trimming tags and primers. We discarded low-quality sequences (i)

176

containing ‘N’ (undetermined nucleotides), (ii) with quality score less than 20 (< Q20),

177

and (iii) with length shorter than 200 bp (Yang et al., 2016). All filtered sequences

178

were then clustered into operational taxonomic units (OTUs) at the threshold of 97%

179

similarity using UPARSE pipeline without discarding singletons (Edgar, 2013). Based

180

on the Ribosomal Database Project (RDP), taxonomic assignment of representative

AC C

EP

TE D

M AN U

SC

RI PT

159

9

ACCEPTED MANUSCRIPT sequences was obtained by searching against the SILVA_128 database at the 85%

182

confidence (Wang et al., 2007). Sequences classified as unknown, archaea, and

183

mitochondria were removed from subsequent analyses. OTU table was normalized by

184

rarifying down to the smallest sequence number. To further assess ecological

185

implication of community shifts, FAPROTAX (Functional Annotation of Prokaryotic

186

Taxa) was used to annotate OTUs to get functional groups (Louca et al., 2016). Using

187

the current literature on cultured strains, FAPROTAX maps microbial taxa to

188

established metabolic or other ecologically relevant functions. Therefore, resampled

189

OTU table was then converted into putative functional profiles, based on taxa

190

identified in each sample. All these analyses were conducted using the in-house

191

pipeline Galaxy (http://mem.rcees.ac.cn:8080).

192

Ecological and statistical analysis

193

Before statistical analyses, all environmental variables, except pH, were log10(x+1)

194

transformed to improve normality. Changes of environmental variables along the river

195

were illustrated with a heatmap and Principal Component Analysis (PCA). Based on

196

the clustering result, the extent of dissimilarity of environmental variables in different

197

river sections was tested using multi-response permutation procedure (MRPP),

198

analysis of similarity (ANOSIM), and analysis of distance matrices (ADONIS). All

199

these analyses were performed using R ‘vegan’ package (Oksanen et al., 2013).

200

AC C

EP

TE D

M AN U

SC

RI PT

181

Rarefaction curves were generated to test the sequencing depth. Resampled OTU

201

table was used to calculate α diversity indices, including the number of observed

202

OTUs, Chao1 value and Shannon index. The difference of α diversity in different 10

ACCEPTED MANUSCRIPT river sections was tested using Mann-Whitney U test. Non-metric multidimensional

204

scaling (NMDS) was used to display the difference of OTU-level community

205

composition along the river. Based on the clustering result, dissimilarity of microbial

206

community composition was tested using MRPP, ANOSIM, and ADONIS. In addition,

207

OTUs appearing in all sites were clustered as the core microbiome (Staley et al.,

208

2014). To investigate the influences of environmental variables on community

209

structure, forward selection was performed to select relatively more important

210

variables (Blanchet et al., 2008). Then a parsimonious redundancy model (RDA) was

211

constructed with those selected variables. The relationship between selected OTUs

212

and selected variables were examined with Spearman’s rank correlation coefficient.

213

All these analyses were performed using R ‘vegan’ package (Oksanen et al., 2013),

214

except correlation analysis and Mann-Whitney U test which were conducted using

215

SPSS 18.0.

TE D

M AN U

SC

RI PT

203

Changes of relative abundance of functional profiles along the river were

217

examined using NMDS. Significance of dissimilarity of functional profiles among

218

different group pairs was tested using MRPP, ANOSIM, and ADONIS. The difference

219

of relative abundance of each functional group in different sections was examined

220

using Mann-Whitney U test. To investigate the influence of environmental variables

221

on functional profiles, forward selection was performed to select relatively more

222

important variables (Blanchet et al., 2008). Then a parsimonious redundancy model

223

(RDA) was constructed with those selected variables. The relationship between

224

selected environmental variables and relative abundance of functional profiles was

AC C

EP

216

11

ACCEPTED MANUSCRIPT tested using Spearman’s rank correlation coefficient. All these analyses were

226

performed using R ‘vegan’ package (Oksanen et al., 2013), except Mann-Whitney U

227

test which was performed using SPSS 18.0.

228

Species interaction networks of bacterioplankton and associations with

229

environmental changes

230

To investigate the relationship between environmental changes and microbial

231

interactions, a full phylogenetic molecular ecological network (pMEN) of all sites was

232

constructed using 16S rRNA sequencing data (Zhou et al., 2010; Deng et al., 2012).

233

According to the instructions of pMEN construction, only OTUs that appeared in at

234

least 20/39 samples and with the relative abundance > 0.1% were kept for network

235

construction (Zhou et al., 2010; Deng et al., 2012). The relationships between

236

microbial network topology and environmental variables were examined using Mantel

237

and partial Mantel tests (Deng et al., 2012). Firstly, the OTU significance (GS) was

238

calculated and defined as the square of the correlation coefficient of OTU abundance

239

profile with environmental traits. Secondly, GS was correlated with OTU connectivity.

240

To display changes of species interactions along the river, individual networks were

241

constructed for each section. Student’s t-test was used to test the difference of network

242

properties, including average connectivity (avgK), modularity, geodesic distance (GD),

243

and clustering coefficient (avgCC). Sub-networks of top OTUs from core microbiome

244

were constructed to understand the succession of inter-species interactions along the

245

river. Network construction was conducted via Galaxy pipeline

246

(http://ieg4.rccc.ou.edu/mena/). The networks were visualized using Cytoscape 3.3.0

AC C

EP

TE D

M AN U

SC

RI PT

225

12

ACCEPTED MANUSCRIPT (Shannon et al., 2003).

248

Results

249

Identification of a nitrogen gradient along the North Canal River

250

PCA plot of environmental variables showed that all the sites were grouped into three

251

distinctive clusters, mainly responding to up- (Section I), middle (Section II) and

252

down streams (Section III) of NCR (Figure 2a). Dissimilarity tests based on the

253

Euclidean distance showed that the clustering of samples was significant (p < 0.05 for

254

all section pairs, supplementary Table S1). In addition, environmental variables varied

255

along the tested river (supplementary Figure S1). For example, the concentrations of

256

nutrients (i.e., TN, NH4-N, NO3-N, TP, and SRP) decreased from upstream towards

257

downstream. Particularly, concentration of TN showed a 100-fold variation, with

258

highest value of 80.642 mg/L (58.129~134.144 mg/L) in Section I, moderate value of

259

29.389 mg/L (5.454~62.015 mg/L) in Section II, and lowest value of 3.621 mg/L

260

(1.043~12.459 mg/L) in Section III. All these analyses showed the existence of a

261

nitrogen gradient along NCR.

262

Variation of α diversity and microbial community composition

263

The high-throughput sequencing generated a total number of 1,349,123 sequences,

264

and the number of sequences per sample ranged from 22,738 to 44,063. The filtered

265

sequences were clustered into 4,567 OTUs. Rarefaction curves indicated that most

266

samples reached or almost reached the saturation stage, suggesting that the majority of

267

biodiversity was recovered based on current sequencing depth (supplementary Figure

268

S2). Three indices of α diversity were calculated, including Chao1 value, the observed

AC C

EP

TE D

M AN U

SC

RI PT

247

13

ACCEPTED MANUSCRIPT number of OTUs, and Shannon index. The former two estimated species richness,

270

while the latter one considered both species richness and evenness. Alpha diversity

271

showed a slightly decreasing trend from upstream towards downstream

272

(supplementary Table S2). Significant difference of α diversity was detected between

273

Sections I and III (two extremes of the nitrogen gradient), including the Shannon

274

index (Section I > Section III; p = 0.000), Chao1 value (Section I > Section III; p =

275

0.015), and number of observed OTUs (Section I > Section III; p = 0.001). The

276

observed number of OTUs was significantly higher in Sections I than II (p = 0.035).

277

For all three sections, microbial community composition varied at different

278

taxonomic levels (supplementary Figure S3). At the phylum level, the dominant phyla

279

of the whole river were Proteobacteria (mostly classes of Gammaproteobacteria and

280

Betaproteobacteria), Bacteroidetes, Actinobacteria, and Firmicutes (Figure 1).

281

However, the relative abundance of each phylum/class was different in the three

282

sections. For example, the relative abundance of Gammaproteobacteria in Section I

283

(18.38%) was substantially lower than that in Section II (37.49%, p = 0.012) and

284

Section III (43.65%, p = 0.000); while the relative abundance of Betaproteobacteria

285

in Section I (35.73%) was much higher than that in Section II (19.56%, p = 0.000) and

286

Section III (10.82%, p = 0.000). We also detected obvious difference of relative

287

abundance at the family level (supplementary Figure S3b). For example, the relative

288

abundance of Comamonadaceae was highest in Section I (19.97%), moderate in

289

Section II (10.02%), and lowest in Section III (3.50%); while the relative abundance

290

of Moraxellaceae was lowest in Section I (8.94%), moderate in Section II (19.12%),

AC C

EP

TE D

M AN U

SC

RI PT

269

14

ACCEPTED MANUSCRIPT and highest in Section III (41.25%).

292

Significant dissimilarity of community β diversity and the correlation with

293

environmental variables

294

Significant difference of community composition between section pairs was identified,

295

with plot of NMDS showing that most samples belonging to the same section were

296

clustered together (Figure 2b). In spite of some deviations, further results of

297

OTU-level community composition analyses confirmed the significant dissimilarity

298

among three sections (p < 0.001 between section pairs; supplementary Table S1).

299

Despite this variation, 1,003 OTUs (21.97% of total OTUs) were shared by three

300

sections (supplementary Figure S4). The percentage of unique OTUs was higher in

301

Section I (47.35%, 1584/3345) than that in Section II (27.06%, 645/2398) and Section

302

III (28.25%, 417/1476). In addition, we found 32 common OTUs (0.7% of total OTUs)

303

in all sites and formed the core microbiome of NCR (Table S3). These OTUs were

304

dominated by Gammaproteobacteria (10 OTUs), Betaproteobacteria (10 OTUs), and

305

Actinobacteria (6 OTUs). They contributed to more than 30% of total sequences for

306

34/39 samples, and more than 50% for 9/39 samples. This was mainly due to the

307

overlap of core OTUs with dominant OTUs. For example, three core OTUs were the

308

top three dominant ones, including OTU_1 (Gammaproteobacteria / Acinetobacter),

309

OTU_2 (Betaproteobacteria / Polynucleobacter), and OTU_3 (Betaproteobacteria /

310

Unclassified).

AC C

EP

TE D

M AN U

SC

RI PT

291

311

Total nitrogen (TN), chlorophyll a (Chl_a), total dissolved solid (TDS),

312

magnesium (Mg), water temperature (T), pH, soluble reactive phosphorous (SRP), 15

ACCEPTED MANUSCRIPT and dissolved oxygen (DO) were selected as crucial environmental variables that

314

structured bacterioplankton biogeography (supplementary Table S4). All these

315

variables explained 40.9% of the total variation of bacterial geographical distribution,

316

with TN alone explaining the proportion of 20.2%. The plot of RDA indicated the

317

separate clusters of samples under different conditions, especially according to the

318

concentration of TN (Figure 3). In addition, TDS and Mg separated Sections I and II

319

from Section III; while Chl_a separated Section I from Sections II and III (Figure 3).

320

In addition, these variables exerted discrepant influence on different taxa. For

321

example, the relative abundance of 16 high fitness OTUs was significantly correlated

322

with TN (Figure 3; Figure S5). The Gammaproteobacteria OTUs were generally

323

negatively affected by TN; while the Betaproteobacteria OTUs were generally

324

positively affected by TN (Figure 3). We also found 17 out of 32 core OTUs were

325

significantly correlated with TN (Table S3), with discrepant responses of different

326

taxa.

327

Distinct co-occurrence networks of bacterioplankton in three sections

328

A full phylogenetic molecular ecological network (pMEN) of all sites was constructed

329

(supplementary Figure S6). The nodes with max degree (OUT_2) and max

330

betweenness (OUT_3) all belonged to Betaproteobacteria, suggesting the importance

331

of this class in the full network. We used the trait-based OTU significance measure to

332

determine the relationships between network interactions and environmental variables.

333

Mantel test revealed significant correlations between connectivity and the OTU

334

significance of all environmental variables in all detected OTUs (p = 0.001) or several

AC C

EP

TE D

M AN U

SC

RI PT

313

16

ACCEPTED MANUSCRIPT phylogenetic groups such as Actinobacteria (p = 0.035), Betaproteobacteria (p =

336

0.001), and Gammaproteobacteria (p = 0.012) along the river (supplementary Table

337

S6). Significant correlations were also observed between connectivity and the OTU

338

significance of selected nutrient variables in all detected OTUs (p = 0.001) and

339

several phylogenetic groups such as Actinobacteria (p = 0.027), Betaproteobacteria (p

340

= 0.001), and Gammaproteobacteria (p = 0.002) (supplementary Table S6). These

341

results suggest that species interactions were dramatically shifted by the

342

environmental gradient along the gradient.

SC

M AN U

343

RI PT

335

In order to determine the succession of species interactions under different environments, individual network of each river section was constructed

345

(supplementary Figure S7). The network size was similar among three networks.

346

However, the network of Section III (1134) contained much more links than Section I

347

(287) and Section II (260). In addition, significant difference was identified for

348

network properties, including average connectivity, path length, and clustering

349

coefficient (supplementary Table S6). All of the three indices measured the

350

connectivity status of one node to other nodes. The higher the avgK and avgCC, and

351

the lower GD indicated the more complex interactions of nodes. The network of

352

Section III was most complex with largest number of links, highest connectivity,

353

smallest path length, as well as highest clustering coefficient. However, in all of the

354

three networks, the majority nodes belonged to Betaproteobacteria,

355

Gammaproteobacteria, and Flavobacteriia (Figure S7).

AC C

EP

TE D

344

17

ACCEPTED MANUSCRIPT 356

In order to show succession of species interactions at a higher resolution, we constructed sub-networks of six top Betaproteobacteria OTUs, which were

358

significantly correlated with environmental factors (especially TN) and appeared in all

359

three networks (Figure 5). All of these OTUs had much more complex interactions in

360

Section III than corresponding OTUs in Section I and Section II (Figure 5).

361

Significant changes of functional groups among three river sections

362

Using FAPROTAX (Louca et al., 2016), 33.41% of observed OTUs could be assigned

363

to at least one functional group, and a total of 69 groups were identified. Dominant

364

functional groups were related to chemoheterotrophy, animal parasites, human

365

pathogens and nitrogen metabolism (supplementary Figure S8). NMDS plot displayed

366

separate clusters of functional groups according to river section (Figure 5a).

367

Dissimilarity tests revealed that the relative abundance of functional groups was

368

significantly different among three sections (p < 0.01 for all section pairs,

369

supplementary Table S1). TN was the most important factor that led to changes of

370

functional groups (supplementary Table S3b). The relative abundance of most

371

nitrogen metabolic OTUs was significantly different among section pairs (Figure 5b).

372

Most importantly, the relative abundance of nitrogen metabolism groups (9/12) in

373

Section I was higher than that in the other two sections. The correlation analysis

374

confirmed that relative abundance of most of these functional groups was significantly

375

correlated with the concentration of TN (supplementary Table S7). By associating the

376

functional groups to bacterial classes, we found Betaproteobacteria contributed a

377

largest fraction to most nitrogen cycle groups (Figure S9). Therefore, the higher

AC C

EP

TE D

M AN U

SC

RI PT

357

18

ACCEPTED MANUSCRIPT abundances of nitrogen cycle groups might result from the higher abundance of

379

Betaproteobacteria in Section I (Figure 1).

AC C

EP

TE D

M AN U

SC

RI PT

378

19

ACCEPTED MANUSCRIPT

Discussion

381

In this study, we explored the impacts of nutrient pollution (especially nitrogen

382

enrichment) on microbial community composition, function, and microbial

383

interactions in a highly eutrophic river - the North Canal River (NCR). NCR was

384

separated into three sections with different levels of nutrient pollution, thus providing

385

a promising model to answer questions associated with impacts of nutrient pollution

386

on microbial communities. In response to environmental variation, community

387

characteristics were significantly changed in different river sections. Total nitrogen

388

was identified as the most important factor driving these changes. Altogether, the

389

results obtained in this study indicate that nutrient pollution exerts a selective pressure

390

on microbial communities, leading to adaptive shifts and inter-specific interactions

391

along the gradient observed.

392

Factors impacting microbial communities in the highly polluted NCR

393

In highly disturbed streams and rivers, nutrient pollution caused by excess loadings of

394

nitrogen and phosphorous has been widely observed (Carpenter et al., 1998; Camargo

395

and Alonso, 2006; Smith and Schindler, 2009). Nutrient pollution could lead to

396

cascading changes of water physicochemical parameters and aquatic biodiversity,

397

which complicates the interpretation of aquatic community dynamics (Dudgeon et al.,

398

2006; Smith and Schindler, 2009). In this study, we found significant changes of

399

abiotic variables, forming the significant environmental gradient along the tested river

400

(Figure 1). In respect to variation of community composition, several dominant

401

environmental factors were identified, including TN, Chl_a, TDS, Mg, T, pH, SRP,

AC C

EP

TE D

M AN U

SC

RI PT

380

20

ACCEPTED MANUSCRIPT and DO. Consequently, the combined effects of different variables have influenced the

403

microbial communities and their geographical distributions. However, TN alone could

404

explain more than half of the total variation explained by all environmental variables.

405

This could be contributed to >100-fold span of TN concentration, ranging from

406

134.14 mg/L in the upstream to 1.03 mg/L in the downstream. In addition, the

407

concentration of TN in the upstream was much higher compared to other polluted

408

rivers (Ibekwe et al., 2016; Wang et al., 2017). Consequently, NCR in our study

409

provided an ideal system to test the influence of nitrogen gradient on microbial

410

communities in the wild.

SC

M AN U

411

RI PT

402

Diverse studies of nitrogen accumulation in the field or nitrogen addition experiments have confirmed the influence of nitrogen on microbial community

413

composition, functional diversity, or microbial biomass (Wang et al., 2016; Ma et al.,

414

2016; Zeng et al., 2016). In aquatic ecosystems, nitrogen pollution was supposed to

415

impact ecological communities through at least two mechanisms: water acidification

416

and/or oxygen depletion (Howarth et al., 2000; Camargo and Alonso, 2006).

417

Following these findings, TN might indirectly function through altering water pH and

418

DO in our study, which could further impact growth and metabolism of microbial

419

species. For example, the heavily polluted upstream (Section I) held lower pH and DO

420

concentration, while the lightly polluted downstream (Section III) held higher pH and

421

DO concentration (Figure S1). In addition, both pH and DO were selected as

422

important factors to drive the variation of microbial community structure, with

423

smaller explanation fractions than TN (Table S3). Therefore, distinct patterns of

AC C

EP

TE D

412

21

ACCEPTED MANUSCRIPT community structure and distribution might result from both the selection of

425

environmental pressures and specific adaptation strategies of microbial species.

426

Adaptive shifts of microbial community composition to nutrient pollution

427

Our study found that microbial community composition was significantly changed in

428

three sections of NCR with TN being the most influential factor. At the phylum level,

429

the most dominant composition of bacterial communities was Proteobacteria (> 50%).

430

However, the relative abundance of Betaproteobacteria was highest in Section I and

431

decreased towards downstream; while Gammaproteobacteria showed the opposite

432

trend (Figure 1). This was consistent with previous studies, affirming the capability of

433

Betaproteobacteria to tolerate massively polluted water (Brümmer et al., 2000; Hahn

434

et al., 2012; Ma et al., 2016). Community changes at the phylum or class levels could

435

be explained by the composition variation at lower taxonomic ranks (Figure S3). For

436

example, Betaproteobacteria in Section I was dominated by Comamonadaceae at the

437

family level (Figure S3). Comamonadaceae was abundant in activated sludge during

438

wastewater treatment and contained diverse denitrifying members (Khan et al., 2002).

439

Therefore, Comamonadaceae could potentially indicated the disturbance of

440

wastewater in the upstream, along which many wastewater treatment plants (WWTPs)

441

were scattered (Heeb et al., 2012). In addition, this family also contained fast-growing

442

and copiotrophic groups (Newton et al., 2013), further indicating that high nutrient

443

concentration may serve as one important driver in the microbial assemblages.

444

Betaproteobacteria in Section I was dominated by Polynucleobacter at the genus

445

level (Figure S3). Polynucleobacter has been widely found in polluted rivers and lives

AC C

EP

TE D

M AN U

SC

RI PT

424

22

ACCEPTED MANUSCRIPT as chemoorganotrophs by utilizing low-molecular-weight substrates derived from

447

photooxidation of humic substances (Hahn et al., 2012; Ma et al., 2016). Despite the

448

significant changes of community composition, 32 OTUs were found to appear in all

449

sites and formed the core microbiome (supplementary Table S4). These OTUs mostly

450

belonged to Betaproteobacteria (i.e. Polynucleobacter) and Gammaproteobacteria

451

(i.e. Pseudomonas). Most of these OTUs also showed adaptive shifts to nutrient

452

pollution with significant correlation between their relative abundance and TN. As

453

these OTUs contributed to a large fraction of microbial abundance, they were greatly

454

responsible for mediating community-level response to nutrient pollution. Microbial

455

community, especially those core groups, could serve as good indicators of water

456

quality and anthropogenic disturbance (Staley et al., 2014). Understanding the process

457

and mechanism of changes of their relative abundance may provide further insights

458

into ecological consequences of nutrient pollution on aquatic microorganisms.

TE D

M AN U

SC

RI PT

446

In our study, α diversity of microbial communities slightly decreased from the

460

upstream towards downstream (supplementary Table S2). This trend has frequently

461

been observed in other polluted river ecosystems (Savio et al., 2015; Wang et al.,

462

2017). Higher biodiversity is assumed to safeguard stronger tolerance to

463

environmental pressures, while the loss of biodiversity would impair ecosystem

464

functioning (Cottingham et al., 2001; Matias et al., 2013). Therefore, the relatively

465

higher α diversity in heavily polluted rivers may contribute to microbial adaptation to

466

aggravation of nutrient pollution. With high concentration of nutrients, Section I

467

might provide diverse and abundant substrates, which favor diverse species. The

AC C

EP

459

23

ACCEPTED MANUSCRIPT nursery effect of headwaters could be one important source of those diverse species

469

(Besemer et al., 2013). For example, microbes derived from riparian soils or

470

groundwater could be flushed into the river and contributed to the increase of

471

biodiversity (Savio et al., 2015; Crump et al., 2012). When they move downstream

472

passively or actively, microbes with weaker tolerance might be filtered out, which

473

also partially explained the lower biodiversity in middle and downstream (Liu et al.,

474

2013; Savio et al., 2015). Consequently, loss of α diversity of bacterioplankton could

475

be an alternative indicator of the influence of nutrient pollution (Savio et al., 2015;

476

Wang et al., 2017).

477

Adaptive shifts of microbial interactions to nutrient pollution

478

Apart from environmental variables, biotic variables such as predation or

479

species-species interactions could also shape the distribution and structure of bacterial

480

communities, and consequently regulate the response of communities to disturbance

481

(Shade et al., 2012). Since it is difficult to directly explore the relationship among

482

diverse species, conceptual frameworks provide powerful tools (Ponomarova and

483

Patil, 2015; Cardona et al., 2016). For example, random matrix theory-based

484

(RMT-based) network construction provides a useful pathway to characterize

485

microbial interactions with high throughput sequencing data (Zhou et al., 2010; Deng

486

et al., 2012). This method has been successfully used in studies in various ecosystems,

487

especially soil microorganisms (Deng et al., 2016; Li et al., 2017). However, few

488

studies have used this method to explore adaptation of microbial interactions to

489

nutrient pollution, especially in river ecosystems (but see Dai et al., 2017; Hu et al.,

AC C

EP

TE D

M AN U

SC

RI PT

468

24

ACCEPTED MANUSCRIPT 2017). In this study, network analyses indicated that nutrient pollution impacted not

491

only bacterial community composition and function, but also species-species

492

interactions. Most importantly, the complexity of individual network was highest in

493

Section III under the moderate eutrophication, while lowest in Section I under the

494

heavy eutrophication. Sub-networks of OTUs from the core microbiome were also

495

relatively more complex in Section III. Based on previous studies, a simple network

496

structure may indicate an unstable and vulnerable microbial community when various

497

disturbances occur (Faust and Raes, 2012; Dai et al., 2017; Wang et al., 2018).

498

Therefore, the heavy eutrophication occurring in the upstream in this study disrupted

499

the species-species interactions, while the moderate eutrophication enhanced the

500

species-species interactions. In addition, central nodes (OTUs) in networks

501

overlapped with taxa of core OTUs (such as Polynucleobacter/Betaproteobacteria),

502

further suggesting the ecological importance of the core microbiome. Consequently,

503

microbial interactions also showed adaptive responses to nutrient pollution. However,

504

more solid experimental evidence is still needed to confirm the direct linkages

505

between nutrient pollution and species-species interactions.

506

Adaptive shifts of microbial functional diversity to nutrient pollution

507

Extensive studies explored the relationship between community structure and

508

functional diversity (Zeglin, 2015; Louca et al., 2016; Maynard et al., 2017). However,

509

it remains unclear to which extent that community functions would be dependent on

510

community composition (Zeglin, 2015; Louca et al., 2016; Maynard et al., 2017).

511

Some studies insisted that functional redundancy could contribute to the weaker

AC C

EP

TE D

M AN U

SC

RI PT

490

25

ACCEPTED MANUSCRIPT diversification of community function than community composition (Staley et al.,

513

2014); while others found composition-dependent changes of community function

514

(Woodward et al., 2012; Freimann et al., 2013). Using FAPROTAX (Louca et al.,

515

2016), we obtained functional annotations of OTUs and revealed significant

516

dissimilarity of functional group profiles in three sections, with TN being the most

517

influential environmental variable (Figure 5). Most importantly, most of the nitrogen

518

metabolism-related groups showed highest abundance in Section I, which was heavily

519

polluted by nitrogen enrichment. With Betaproteobacteria contributed dominantly to

520

those nitrogen-metabolism groups, our results suggest the potential coupling of

521

community composition and functions. This phenomenon was consistent with other

522

studies in polluted rivers or streams where coupling of community composition and

523

functions was detected (Jeffries et al., 2016). Therefore, bacterial community showed

524

functional adaptation to nutrient pollution and conceived potentials to attenuate

525

negative impacts of environmental changes (Li et al., 2017). In addition, many

526

functional groups were found to be human or animal pathogens (supplementary

527

Figure S5), suggesting aquatic microorganisms could also be the sources of diseases

528

and impose threats/risks to human health (Ibekwe et al., 2016; Li et al., 2017). Given

529

the significance of microorganisms’ participation in fundamental biogeochemical

530

cycles, more efforts are deserved to assess the interactions between water pollution

531

and dynamics of aquatic microbial community (Woodward et al., 2012). In addition,

532

seasonal and yearly monitoring is necessary in future studies as both items would

533

change temporally, as well as their relationship (Fortunato et al., 2012; Gonzalez et al.,

AC C

EP

TE D

M AN U

SC

RI PT

512

26

ACCEPTED MANUSCRIPT 534

2012).

AC C

EP

TE D

M AN U

SC

RI PT

535

27

ACCEPTED MANUSCRIPT

Conclusions

537

This study showed that bacterial community composition, potential interactions, and

538

community function were all significantly affected by nutrient pollution, especially

539

nitrogen enrichment in the highly polluted river. The shifts of microbial community

540

composition and functional groups were caused by several dominant factors, with

541

total nitrogen (TN) being the most important one. This could be derived

542

from >100-fold changes of TN concentration along the whole river. Bacterial groups

543

with strong tolerance to nutrient pollution were enriched in river sections with heavy

544

eutrophication. Meanwhile, functional prediction revealed highest abundance of

545

nitrogen metabolic groups in the heavily polluted river section. The relative

546

abundance of these functional groups were significantly correlated with TN.

547

Consequently, microbial communities present in certain areas were tolerant to

548

distinctive levels of nutrient pollution (especially nitrogen enrichment). However, the

549

co-occurrence networks were dramatically simple and vulnerable under heavy nutrient

550

pollution. This hints that the potential interactions of microbial communities were

551

impacted by nutrient pollution, which would further impair ecosystem functioning.

552

Therefore, species interactions should also be taken into consideration when studying

553

potential influence of water pollution on biological communities.

AC C

EP

TE D

M AN U

SC

RI PT

536

554

28

ACCEPTED MANUSCRIPT

Acknowledgements

556

This work was supported by the National Natural Science Foundation of China [grant

557

nos.: 31572228; 31800419], the Innovation in Cross-functional Team Program of the

558

Chinese Academy of Sciences [grant no.: No.: 2015], Chinese Academy of Science

559

[grant no.: ZDRW-ZS-2016-5], the State Key Joint Laboratory of Environment

560

Simulation and Pollution Control [RCEES, Chinese Academy of Sciences; grant no.:

561

15K01ESPCR], and the Water Pollution Control and Treatment Special Project [grant

562

no.: 2015ZX07201-008-09].

563

Author contributions

564

A.Z. and Y.Y. designed the experiment; Y.Y., Y.G., and X.H. collected samples in the

565

field; Y.Y., P.N., Y.W., and Y.D. conducted the laboratory analysis of water

566

physiochemical variables; Y.Y. analyzed the data, mad the figures and tables, and

567

wrote the manuscript. All authors contributed to the revision of this manuscript.

568

Data accessibility

569

The high-throughput sequencing data will be deposited into NCBI Sequence Read

570

Archive (SRA) upon acceptation.

AC C

EP

TE D

M AN U

SC

RI PT

555

29

ACCEPTED MANUSCRIPT 571

References

572

Azam F, Malfatti F. (2007). Microbial structuring of marine ecosystems. Nat Rev

574

Microbiol 5:782. Besemer K, Singer G, Quince C, Bertuzzo E, Sloan W, Battin TJ. (2013). Headwaters

RI PT

573

575

are critical reservoirs of microbial diversity for fluvial networks. Proc R Soc B

576

280:20131760.

579

SC

578

Blanchet FG, Legendre P, Borcard D. (2008). Forward selection of explanatory variables. Ecology 89:2623-2632.

M AN U

577

Brümmer IHM, Fehr W, Wagner-Döbler I. (2000). Biofilm community structure in

580

polluted rivers: abundance of dominant phylogenetic groups over a complete

581

annual cycle. Appl Environ Microb 66:3078-3082.

Camargo JA, Alonso Á. (2006). Ecological and toxicological effects of inorganic

TE D

582 583

nitrogen pollution in aquatic ecosystems: a global assessment. Environ Int

584

32:831-849.

587 588 589 590 591 592

EP

586

Canfield DE, Glazer AN, Falkowski PG. (2010). The evolution and future of Earth’s nitrogen cycle. Science 330:192-196.

AC C

585

Caporaso JG, Lauber CL, Walters WA, Berg-Lyons D, Lozupone CA, Turnbaugh PJ, et al. (2011). Global patterns of 16S rRNA diversity at a depth of millions of sequences per sample. Proc Natl Acad Sci USA 108:4516-4522.

Cardona C, Weisenhorn P, Henry C, Gilbert JA. (2016). Network-based metabolic analysis and microbial community modeling. Curr Opin Microbiol 31:124-131. Carpenter SR, Caraco NF, Correll DL, Howarth RW, Sharpley AN, Smith VH. (1998). 30

ACCEPTED MANUSCRIPT 593

Nonpoint pollution of surface waters with phosphorus and nitrogen. Ecol Appl

594

8:559-568. Chafee M, Fernàndez-Guerra A, Buttigieg PL, Gerdts G, Eren AM, Teeling H, et al.

596

(2017). Recurrent patterns of microdiversity in a temperate coastal marine

597

environment. ISME J 12:237.

600

temporal variability of ecological systems. Ecol Lett 4:72-85.

SC

599

Cottingham KL, Brown BL, Lennon JT. (2001). Biodiversity may regulate the

Crump BC, Amaral-Zettler LA, Kling GW. (2012). Microbial diversity in arctic

M AN U

598

RI PT

595

601

freshwaters is structured by inoculation of microbes from soils. ISME J

602

6:1629–1639.

Dai W, Zhang J, Tu Q, Deng Y, Qiu Q, Xiong J. (2017). Bacterioplankton assembly

604

and interspecies interaction indicating increasing coastal eutrophication.

605

Chemosphere 177:317-325.

606

TE D

603

Datta MS, Sliwerska E, Gore J, Polz MF, Cordero OX. (2016). Microbial interactions lead to rapid micro-scale successions on model marine particles. Nat Commun

608

7:11965.

610 611

AC C

609

EP

607

Deng Y, Jiang YH, Yang Y, He Z, Luo F, Zhou J. (2012). Molecular ecological network analyses. Bioinformatics 13:113.

Deng Y, Zhang P, Qin Y, Tu Q, Yang Y, He Z. (2016). Network succession reveals the

612

importance of competition in response to emulsified vegetable oil amendment for

613

uranium bioremediation. Environ Microbiol 18:205-218.

614

Dudgeon D, Arthington AH, Gessner MO, Kawabata ZI, Knowler DJ, Lévêque C, et 31

ACCEPTED MANUSCRIPT 615

al. (2006). Freshwater biodiversity: importance, threats, status and conservation

616

challenges. Bio Rev 81:163-182.

619 620 621

amplicon reads. Nat Methods 10:996-998.

RI PT

618

Edgar RC. (2013). UPARSE: highly accurate OTU sequences from microbial

Faust K, Raes J. (2012). Microbial interactions: from networks to models. Nat Rev Microbiol 10:538.

Fortunato CS, Herfort L, Zuber P, Baptista AM, Crump BC. (2012). Spatial variability

SC

617

overwhelms seasonal patterns in bacterioplankton communities across a river to

623

ocean gradient. ISME J 6:554.

624

M AN U

622

Freimann R, Bürgmann H, Findlay SE, Robinson CT. (2013). Bacterial structures and ecosystem functions in glaciated floodplains: contemporary states and potential

626

future shifts. ISME J 7:2361-2373.

629 630 631 632 633 634

Trends Microbiol 24:833-845.

Gonzalez A, King A, Robeson II MS, Song S, Shade A, Metcalf JL, et al. (2012).

EP

628

Ghoul M, Mitri S. (2016). The ecology and evolution of microbial competition.

Characterizing microbial communities through space and time. Curr Opin

AC C

627

TE D

625

Biotech 23:431-436.

Guo J, Jing H, Li J, Li J. (2012). Surface water quality of Beiyun River basin and the analysis of acting factors for the recent ten years. Environ Sci 33:1511-1518.

Hahn MW, Scheuerl T, Jezberová J, Koll U, Jezbera J, Šimek K, et al. (2012). The

635

passive yet successful way of planktonic life: genomic and experimental analysis

636

of the ecology of a free-living Polynucleobacter population. PLoS One 7:e32772. 32

ACCEPTED MANUSCRIPT 637

He B, Kanae S, Oki T, Hirabayashi Y, Yamashiki Y, Takara K. (2011). Assessment of

638

global nitrogen pollution in rivers using an integrated biogeochemical modeling

639

framework. Water Res 45:2573-2586. Heeb F, Singer H, Pernet-Coudrier B, Qi W, Liu H, Longrée P, et al. (2012). Organic

RI PT

640

micropollutants in rivers downstream of the megacity Beijing: sources and mass

642

fluxes in a large-scale wastewater irrigation system. Environ Sci Technol

643

46:8680-8688.

645 646

Howarth RW, Anderson DB, Cloern JE, Elfring C, Hopkinson CS, Lapointe B, et al.

M AN U

644

SC

641

(2000). Issues in ecology: Nutrient pollution of coastal rivers, bays, and seas. Hu A, Ju F, Hou L, Li J, Yang X, Wang H, et al. (2017). Strong impact of anthropogenic contamination on the co-occurrence patterns of a riverine

648

microbial community. Environ Microbiol 19:4993-5009.

649

TE D

647

Huang WE, Ferguson A, Singer AC, Lawson K, Thompson IP, Kalin RM, et al. (2009). Resolving genetic functions within microbial populations: in situ analyses using

651

rRNA and mRNA stable isotope probing coupled with single-cell

652

Raman-fluorescence in situ hybridization. Appl Environ Microb 75:234-241.

654 655

AC C

653

EP

650

Ibekwe AM, Ma J, Murinda SE. (2016). Bacterial community composition and structure in an Urban River impacted by different pollutant sources. Sci Total Environ 566:1176-1185.

656

Jeffries TC, Schmitz Fontes ML, Harrison DP, Van-Dongen-Vogels V, Eyre BD, Ralph

657

PJ, et al. (2016). Bacterioplankton dynamics within a large anthropogenically

658

impacted urban estuary. Front Microbiol 6:1438. 33

ACCEPTED MANUSCRIPT 659 660 661

Jones SE, Cadkin TA, Newton RJ, McMahon KD. (2012). Spatial and temporal scales of aquatic bacterial beta diversity. Front Microbio 3:318. Khan ST, Horiba Y, Yamamoto M, Hiraishi A. (2002). Members of the family Comamonadaceae as primary poly (3-hydroxybutyrate-co-3-hydroxyvalerate) -

663

degrading denitrifiers in activated sludge as revealed by a polyphasic approach.

664

Appl Environ Microbiol 68:3206-3214.

Li X, Meng D, Li J, Yin H, Liu H, Liu X, et al. (2017). Response of soil microbial

SC

665

RI PT

662

communities and microbial interactions to long-term heavy metal contamination.

667

Environ Pollut 231:908-917.

668

M AN U

666

Lima-Mendez G, Faust K, Henry N, Decelle J, Colin S, Carcillo F, et al. (2015). Determinants of community structure in the global plankton interactome. Science,

670

348:1262073.

TE D

669

Liu L, Yang J, Yu X, Chen G, Yu Z. (2013). Patterns in the composition of microbial

672

communities from a subtropical river: effects of environmental, spatial and

673

temporal factors. PLoS One 8:e81232.

675 676 677 678 679 680

Louca S, Parfrey LW, Doebeli M. (2016). Decoupling function and taxonomy in the

AC C

674

EP

671

global ocean microbiome. Science 353:1272-1277.

Ma L, Mao G, Liu J, Gao G, Zou C, Bartlam MG, et al. (2016). Spatial-temporal changes of bacterioplankton community along an exhorheic river. Front Microbiol 7:250. Matias M, Combe M, Barbera C, Mouquet N. (2013). Ecological strategies shape the insurance potential of biodiversity. Front Microbiol 3:432. 34

ACCEPTED MANUSCRIPT 681

Maynard DS, Crowther TW, Bradford MA. (2017). Competitive network determines

682

the direction of the diversity–function relationship. Proc Natl Acad Sci USA

683

201712211. Morriën E, Hannula SE, Snoek LB, Helmsing NR, Zweers H, De Hollander M, et al.

RI PT

684 685

(2017). Soil networks become more connected and take up more carbon as nature

686

restoration progresses. Nat Commun 8:14349.

Newton RJ, Huse SM, Morrison HG, Peake CS, Sogin ML, McLellan SL. (2013).

SC

687

Shifts in the microbial community composition of gulf coast beaches following

689

beach oiling. PLoS One 8:e74265.

690

M AN U

688

Niño-García JP, Ruiz-González C, del Giorgio PA. (2016). Interactions between hydrology and water chemistry shape bacterioplankton biogeography across

692

boreal freshwater networks. ISME J 10:1755.

695 696 697

‘vegan’. Community ecology package, version, 2(9). Paerl HW, Dyble J, Moisander PH, Noble RT, Piehler MF, Pinckney JL, et al. (2003).

EP

694

Oksanen J, Blanchet FG, Kindt R, Legendre P, Minchin PR, et al. (2013). Package

Microbial indicators of aquatic ecosystem change: current applications to

AC C

693

TE D

691

eutrophication studies. FEMS Microbiol Ecol 46:233-246.

698

Pernet-Coudrier B, Qi W, Liu H, Müller B, Berg M. (2012). Sources and pathways of

699

nutrients in the semi-arid region of Beijing–Tianjin, China. Environ Sci Technol

700 701 702

46:5294-5301. Ponomarova O, Patil KR. (2015). Metabolic interactions in microbial communities: untangling the Gordian knot. Curr Opin Microbiol 27:37-44. 35

ACCEPTED MANUSCRIPT 703

Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. (2012). The SILVA

704

ribosomal RNA gene database project: improved data processing and web-based

705

tools. Nucl Acids Res 41:D590-D596.

709 710

RI PT

708

Catchment-scale biogeography of riverine bacterioplankton. ISME J 9:516-526. Savio D, Sinclair L, Ijaz UZ, Parajka J, Reischer GH, Stadler P, et al. (2015). Bacterial diversity along a 2600 km river continuum. Environ Microbiol 17:4994-5007.

SC

707

Read DS, Gweon HS, Bowes MJ, Newbold LK, Field D, Bailey MJ, et al. (2015).

Shade A, Peter H, Allison SD, Baho DL, Berga M, Bürgmann H, et al. (2012).

M AN U

706

711

Fundamentals of microbial community resistance and resilience. Frontiers

712

Microbiol 3:417.

713

Shan B, Jian Y, Tang W, Zhang H. (2012). Temporal and spatial variation of nitrogen and phosphorous and eutrophication assessment in downstream river network

715

areas of North Canal river watershed. Environ Sci 33:352-358.

TE D

714

Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. (2003).

717

Cytoscape: a software environment for integrated models of biomolecular

718

interaction networks. Genome Res 13:2498-2504.

720 721

AC C

719

EP

716

Smith VH, Schindler DW. (2009). Eutrophication science: where do we go from here?. Trends Ecol Evol 24:201-207.

Staley C, Gould TJ, Wang P, Phillips J, Cotner JB, Sadowsky MJ. (2014). Core

722

functional traits of bacterial communities in the Upper Mississippi River show

723

limited variation in response to land cover. Front Microbiol 5:414.

724

Sun MY, Dafforn KA, Brown MV, Johnston EL. (2012). Bacterial communities are 36

ACCEPTED MANUSCRIPT 725 726

sensitive indicators of contaminant stress. Mar Pollut Bull 64:1029-1038. Wang J, Pan F, Soininen J, Heino J, Shen J. (2016). Nutrient enrichment modifies temperature-biodiversity relationships in large-scale field experiments. Nature

728

Commun 7:13960.

RI PT

727

Wang P, Wang X. Wang C, Miao L, Hou J, Yuan Q. (2017). Shift in bacterioplankton

730

diversity and structure: Influence of anthropogenic disturbances along the

731

Yarlung Tsangpo River on the Tibetan Plateau, China. Sci Rep 7:12529.

Wang Q, Garrity GM, Tiedje JM, Cole JR. (2007). Naive Bayesian classifier for rapid

M AN U

732

SC

729

733

assignment of rRNA sequences into the new bacterial taxonomy. Appl Environ

734

Microbiol 73:5261–5267.

Wang S, Wang X, Han X, Deng Y, Fortin MJ. (2018). Higher precipitation strengthens

736

the microbial interactions in semi-arid grassland soils. Global Ecol Biogeogr

737

00:1-11.

738

TE D

735

Woodward G, Gessner MO, Giller PS, Gulis V, Hladyz S, Lecerf A, et al. (2012). Continental-scale effects of nutrient pollution on stream ecosystem functioning.

740

Science 336:1438-1440.

742 743 744

AC C

741

EP

739

Xiong W, Li J, Chen Y, Shan B, Wang W, Zhan A. (2016). Determinants of community structure of zooplankton in heavily polluted river ecosystems. Sci Rep 6:22043.

Yang Y, Deng Y, Cao L. (2016). Characterising the interspecific variations and

745

convergence of gut microbiota in Anseriformes herbivores at wintering areas. Sci

746

Rep 6:32655. 37

ACCEPTED MANUSCRIPT 747

Yang Y, Ni P, Gao Y, Xiong W, Zhao Y, Zhan A. (2018). Geographical distribution of zooplankton biodiversity in highly polluted running water ecosystems: Validation

749

of fine-scale species sorting hypothesis. Ecol Evol 8:4830-4840.

750

Yu Y, Wang X, Zhang P. (2012). Spatial distribution of planktonic bacterial

RI PT

748

751

community and its relationship to water quality in Beiyun River. Asian J Ecotox

752

7:337-344.

755

SC

754

Zeglin LH. (2015). Stream microbial diversity in response to environmental changes: review and synthesis of existing research. Front Microbiol 6:454.

M AN U

753

Zeng J, Liu X, Song L, Lin X, Zhang H, Shen C, Chu H. (2016). Nitrogen fertilization

756

directly affects soil bacterial diversity and indirectly affects bacterial community

757

composition. Soil Biol Biochem 92:41-49.

Zhan A, Bailey SA, Heath DD, Macisaac HJ. (2014). Performance comparison of

TE D

758 759

genetic markers for high-throughput sequencing-based biodiversity assessment in

760

complex communities. Mol Ecol Resour 14:1049-1059.

763 764 765

EP

762

Zinger L, Gobet A, Pommier T. (2012). Two decades of describing the unseen majority of aquatic microbial diversity. Mol Ecol 21:1878-1896.

AC C

761

Zhou J, Deng Y, Luo F, He Z, Tu Q, Zhi X. (2010). Functional molecular ecological networks. MBio 1:e00169-10.

38

ACCEPTED MANUSCRIPT Figure 1. Study area and sampling sites. A total of 39 sampling sites were chosen

767

along the North Canal River (NCR). The abundance of dominant phyla in each river

768

section was shown in pie charts. The region in the left bottom is the Haihe River

769

Basin where NCR is located.

M AN U

SC

RI PT

766

EP AC C

771

TE D

770

39

ACCEPTED MANUSCRIPT Figure 2. Clustering results of environmental variables and bacterial community

773

composition along the North Canal River (NCR). (a) Principal component analysis

774

(PCA) plot of all the 25 environmental variables based on Euclidean distance; (b)

775

Non-metric multidimensional scaling (NMDS) plot of bacterial community

776

composition based on Bray-curtis distance at the operational taxonomic unit (OTU)

777

level.

EP

779

AC C

778

TE D

M AN U

SC

RI PT

772

40

ACCEPTED MANUSCRIPT Figure 3. The ordination plot of redundancy analysis (RDA) of bacterial community

781

composition at all sites. Operational taxonomic units (OTUs) weakly associated with

782

the first two axes (with fitness < 60%) were omitted for clarity. The size of all sites

783

(dots in the figure) wase scaled according to the concentration of TN. Black arrows

784

represent selected environmental variables, and blue arrows represent OTUs.

EP

786

AC C

785

TE D

M AN U

SC

RI PT

780

41

ACCEPTED MANUSCRIPT Figure 4. The dynamic network connections of top six operational taxonomic units

788

(OTUs) and their first neighbors. (a) Network interactions of the top six OTUs in

789

Section I; (b) Network interactions of the corresponding OTUs in Section II; (c)

790

Network interactions of the corresponding OTUs in Section III. Blue and red links

791

represented positive and negative interactions, respectively.

M AN U

SC

RI PT

787

792

AC C

EP

TE D

793

42

ACCEPTED MANUSCRIPT Figure 5. Variation of relative abundance of functional groups. (a) Non-metric

795

multidimensional scaling (NMDS) plot of all the 69 categories based on bray-curtis

796

distance; (b) Comparison of dominant nitrogen relevant groups. a, b, and c show the

797

significant difference (p < 0.05) between Sections I and II, I and III, and II and III,

798

repectively.

EP AC C

799

TE D

M AN U

SC

RI PT

794

43

ACCEPTED MANUSCRIPT Highlights •

Nitrogen pollution was crucial to structure microbial communities.

•

Microbial communities showed diverse changes in response to nitrogen

•

RI PT

pollution. Microbial interactions were potentially changed in response to nitrogen pollution.

SC

Potential nitrogen degradation microbes were enriched under severe pollution.

AC C

EP

TE D

M AN U

•