Bacterial diversity across a highly stratified ecosystem: A salt-wedge Mediterranean estuary

Accepted Manuscript Title: Bacterial diversity across a highly stratified ecosystem: A salt-wedge Mediterranean estuary ˇ Author: M. Korlevi´c L. Supraha Z. Ljubeˇsi´c J. Henderiks I. Cigleneˇcki J. Dautovi´c S. Orli´c PII: DOI: Reference:

S0723-2020(16)30051-0 http://dx.doi.org/doi:10.1016/j.syapm.2016.06.006 SYAPM 25784

To appear in: Received date: Revised date: Accepted date:

1-3-2016 21-6-2016 23-6-2016

ˇ Please cite this article as: M. Korlevi´c, L. Supraha, Z. Ljubeˇsi´c, J. Henderiks, I. Cigleneˇcki, J. Dautovi´c, S. Orli´c, Bacterial diversity across a highly stratified ecosystem: a salt-wedge Mediterranean estuary, Systematic and Applied Microbiology (2016), http://dx.doi.org/10.1016/j.syapm.2016.06.006 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.

1

Bacterial diversity across a highly stratified ecosystem: a salt-wedge Mediterranean

2

estuary

3

M. Korlević1, L. Šupraha2, Z. Ljubešić3, J. Henderiks2, I. Ciglenečki4, J. Dautović4, S.

4

Orlić5,6*

3

7

5

9 10 11 12

6

us

Division for Marine and Environmental Research, Ruđer Bošković Institute, Zagreb, Croatia Division of Material Chemistry, Ruđer Bošković Institute, Zagreb, Croatia

an

4

Division of Biology, Faculty of Science, University of Zagreb, Croatia

Center of Excellence for Science and Technology Integrating Mediterranean Region, Microbial Ecology, Zagreb, Croatia

M

8

Paleobiology, Department of Earth Sciences, Uppsala University, Sweden

cr

2

6

Center for Marine Research, Ruđer Bošković Institute, Rovinj, Croatia

ip t

1

5

*Corresponding author: [email protected]

Ac ce

pt

ed

13

Page 1 of 31

Abstract

14

Highly stratified Mediterranean estuaries are unique environments where the tidal range is

15

low and the tidal currents are almost negligible. The main characteristics of these

16

environments are strong salinity gradients and other environmental parameters. In this study,

17

454 pyrosequencing of the 16S rRNA gene in combination with catalyzed reporter deposition-

18

fluorescence in situ hybridization (CARD-FISH) was used to estimate the bacterial diversity

19

across the Krka estuary in February and July 2013. The comparison of the data derived from

20

these two techniques resulted in a significant but weak positive correlation (R=0.28)

21

indicating a substantial difference in the bacterial community structure, depending on the

22

applied method. The phytoplankton bloom observed in February was identified as one of the

23

main factors shaping the bacterial community structure between the two environmentally

24

contrasting sampling months. Roseobacter, Bacteroidetes and Gammaproteobacteria differed

25

substantially between February and July. Typical freshwater bacterial classes (Actinobacteria

26

and Betaproteobacteria) showed strong vertical distribution patterns depending on the salinity

27

gradient. Cyanobacteria decreased in abundance in February due to competition with

28

phytoplankton, while the SAR11 clade increased its abundance in July as a result of a better

29

adaptation towards more oligotrophic conditions. The results provided the first detailed

30

insight into the bacterial diversity in a highly stratified Mediterranean karstic estuary.

cr

us

an

M

ed

pt Ac ce

31

ip t

13

Page 2 of 31

31 Introduction

33

Estuaries are dynamic ecosystems where the interaction of freshwater and seawater leads to

34

the formation of specific environments strongly influenced by a combination of physical,

35

chemical and biological drivers. The riverine input of terrigenous sediments and organic

36

matter, together with the mixing processes related to tides, represent some of the main factors

37

that characterize these environments. Salt-wedge estuaries are restricted to the areas of low

38

tides where a sharp halocline differentiates a freshwater layer and a marine layer resulting in a

39

strong vertical contrast ecosystem [18]. Bacterial production and biomass in estuarine

40

ecosystems has been well studied [22,23,47,49]. However, studies describing estuarine

41

communities are not common and have been mainly based on 16S rRNA sequencing

42

[10,20,33,36] or FISH approaches [1,5,70]. In recent years, next-generation sequencing has

43

been used successfully to examine biogeographic patterns in coastal areas and estuaries

44

[10,20]. Studies describing the diversity of bacteria in estuarine ecosystems have identified

45

freshwater-specific and seawater-specific bacterial communities [14,15]. In addition, they

46

have shown the formation of particular estuarine bacterial communities depending on the

47

freshwater residence time [15]. Freshwater communities were commonly characterized by the

48

presence of Actinobacteria and Betaproteobacteria, while the marine communities were

49

dominated by Alphaproteobacteria (mainly the SAR11 clade). On the other hand,

50

Bacteroidetes and Gammaproteobacteria did not show a clear relationship with the salinity

51

gradient in estuarine ecosystems [1,5,20,36]. The complex dynamics of environmental

52

parameters in estuarine ecosystems represent a challenge for the study of microbial

53

communities, but they also provide a unique system for the investigation of the community

54

changes that occur in these steep gradients.

55

Highly stratified estuaries are characteristic for the areas of the Mediterranean where the tidal

56

range is low and the tidal currents are almost negligible [18]. The karstic Krka River estuary

57

is situated on the eastern coast of the Adriatic Sea (Croatia). The river is characterized by a

58

low terrigenous input due to its karstic drainage area and a series of travertine barriers, lakes

59

and waterfalls that occur before the 23 km long estuary is reached [43]. The presence of a

60

permanent vertical stratification classifies this estuary as a highly stratified type (salt-wedge

61

estuary) where the depth and thickness of the halocline vary depending on the freshwater

62

input. These features make the Krka estuary suitable for studying the effects of environmental

Ac ce

pt

ed

M

an

us

cr

ip t

32

Page 3 of 31

forcing on the formation of an organic matter layer at the halocline [71]. Phytoplankton

64

community dynamics [8,13,60,66], bacterial abundance, activity and biomass [22,23],

65

together with the dynamics of physical and chemical parameters [12,28,38–40,60,71], have

66

been largely studied in the Krka River estuary, but no systematic estimate of bacterial

67

diversity or dynamics has been carried out to date. The distribution of Synechococcus and

68

Prochlorococcus was only recently determined by Šantić et al. [58] using flow cytometry.

69

Although estuaries are very dynamic and important environments and the microbial

70

communities in many of them have been described, little data is available specifically on

71

bacterial communities in Mediterranean estuaries [64]. Therefore, this current study describes

72

the bacterial diversity along the Krka River estuary in February and July using a combination

73

of 454 pyrosequencing and catalyzed reporter deposition-fluorescence in situ hybridization

74

(CARD-FISH) techniques.

75

Materials and methods

76

Sampling and estimation of environmental parameters

77

Two cruises with the R/V Hidra were conducted in the Krka River estuary in winter (25

78

February) and summer (8 July) of 2013. Samples were collected at three stations located

79

within the Krka River Estuary: E3, E4a and E5, whereas Station AD3, located outside the

80

estuary, was chosen as a reference marine coastal station (Fig. 1). Sampling depths were

81

determined after obtaining the salinity and temperature profile using the Seabird SBE 19 plus

82

CTD probe (SEA-Bird Electronics Inc., USA). The samples for the determination of bacterial

83

community structure, phytoplankton, nutrients and chlorophyll a (Chl a) were sampled with 5

84

L Niskin bottles.

85

The concentrations of nitrate (NO3-), nitrite (NO2-), orthophosphate (PO43-) and orthosilicate

86

(SiO44-) were determined according to Strickland and Parsons [57], while the concentration of

87

ammonia (NH4+)[CJR1] was determined according to Ivančić and Degobbis [32]. Dissolved

88

inorganic nitrogen (DIN) was expressed as the sum of the nitrate, nitrite and ammonia

89

[CJR2]concentrations.

90

carbon (DOC; <0.7 µm) were determined by a high-temperature catalytic oxidation

91

(HTCO)[CJR3] method [16] with a TOC-VCPH-5000 solid-sample total organic carbon

92

(TOC)[CJR4] analyzer (Shimadzu, Japan) coupled to an SSM-5000A solid-sample combustion

93

unit (Shimadzu, Japan), according to previously described protocols [17,45]. For Chl a

Ac ce

pt

ed

M

an

us

cr

ip t

63

The concentrations of particulate (POC; >0.7 µm) and dissolved organic

Page 4 of 31

determination, 1 L of seawater was filtered onto 25 mm GF/F filters (Whatman, UK), stored

95

immediately in liquid nitrogen and afterwards at -80 °C, and the Chl a concentrations were

96

determined by reversed-phase high-performance liquid chromatography (HPLC; Spectra

97

System, Model UV 2000) [4]. The samples for phytoplankton analyses were fixed with

98

formaldehyde (2% v/v final concentration) immediately after collection.

99

Phytoplankton determination

ip t

94

Sub-samples of 50 mL were analysed by a Zeiss Axiovert 200 inverted microscope (Zeiss,

101

Germany) after sedimentation for 24 h [65]. Cells larger than 20 μm were designated as

102

microphytoplankton, and cells between 2-20 μm as nanophytoplankton.

103

454 pyrosequencing

104

Seawater aliquots of 1 L were vacuum-filtered through 0.2 µm Nuclepore™ polycarbonate

105

membrane filters (Whatman, UK) with a peristaltic pump. Filters were stored in 1 mL sucrose

106

buffer (40 mM EDTA, 50 mM Tris-HCl and 0.75 M sucrose), frozen in liquid nitrogen and

107

subsequently stored at -80 °C. The DNA was extracted according to Massana et al. [42]. The

108

bacterial V1-V2 16S rRNA region was amplified using bacterial primers 27Fmod (5’-

109

AGRGTTTGATCMTGGCTCAG-3’) and 519Rmodbio (5’-GTNTTACNGCGGCKGCTG-

110

3’) in four parallel reactions. The primers were modified for 454 pyrosequencing, so that each

111

forward primer contained a gene specific sequence (27Fmod) extended at the 5’-end with a 10

112

bp barcode sequence (specific for each sample) and an adapter sequence A, while the reverse

113

primer contained a gene specific sequence (519Rmodbio) extended at the 5’-end with an

114

adapter sequence B. Each 25 µL PCR reaction contained: 1x Green GoTaq® Flexi Buffer

115

(Promega, USA), 1.5 mM MgCl2, 0.2 mM of each deoxynucleoside triphosphate, 0.15 mg

116

BSA, 0.2 µM of forward and reverse primers, 0.625 U of GoTaq® Flexi DNA Polymerase

117

(Promega, USA) and 20 ng of DNA template. The PCR amplification conditions were: 5 min

118

initial denaturation at 95 ºC, 30 cycles of 40 s denaturation at 94 ºC, 40 s annealing at 53 ºC

119

and 1 min elongation at 70 ºC, finalized by 10 min at 70 ºC. After pooling of the replicate

120

reactions, PCR products were purified using the Wizard® SV Gel and PCR Clean-Up System

121

(Promega, USA) and sent for GS FLX Titanium and GS FLX+ 454 pyrosequencing at

122

Eurofins (Ebersberg, Germany).

123

Sequence analysis

Ac ce

pt

ed

M

an

us

cr

100

Page 5 of 31

Obtained standard flowgram format (SFF) files were extracted using an sff_extract script

125

(available at http://bioinf.comav.upv.es/sff_extract/index.html) applying the sff_extract -c

126

command that allows sequence quality checking. Fasta files were split according to the

127

barcode sequence using mothur [52]. Sequences containing any differences in the barcode or

128

primer sequence were removed in the barcode splitting step. Multifasta files were processed

129

by the SILVAngs pipeline (https://www.arb-silva.de/ngs) [50], as described in Ionescu et al.

130

[31]. Briefly, sequences were aligned against the SILVA SSU rRNA SEED using the SILVA

131

Incremental Aligner (SINA) [48]. Sequences with low alignment quality (50 alignment

132

identity, 40 alignment score reported by SINA) were removed (putative contaminations and

133

artifacts). An additional quality check was carried out by removing all sequences shorter than

134

200 nucleotides, with more than 2% ambiguities or 2% homopolymers. Identical sequences

135

were identified (de-replication) and clustered (operational taxonomic units [OTU]) at 97%

136

sequence identity using cd-hit-est (version 3.1.2; http://www.bioinformatics.org/cd-hit) [41]

137

running in accurate mode and ignoring overhangs. The representative OTU sequence was

138

classified against the SILVA SSU Ref dataset (release 115; http://www.arb-silva.de) using

139

blastn (version 2.2.22+; http://blast.ncbi.nlm.nih.gov/Blast.cgi) with standard settings

140

(Dataset 1) [9]. Sequences obtained in this study have been submitted to the European

141

Nucleotide Archive (ENA) under accession numbers ERS845246 to ERS845301.

142

454 pyrosequencing of 28 samples yielded a total of 336,336 pyrotags (an average of 12,012

143

± 1,736 pyrotags per sample), with the range of the number of pyrotags in individual samples

144

being 9,799-18,378. The average pyrotag length was 464 ± 3 bp with a minimum length of

145

200 bp. On average, 33 ± 21 pyrotags were rejected in the quality control process. The vast

146

majority of pyrotags were successfully taxonomically assigned (11,649 ± 1,684 pyrotags per

147

sample). Each sample contained an average of 1,277 ± 157 OTUs (97% sequence identity), of

148

which 44% were singletons (Table S1). The sequencing effort applied was insufficient to

149

determine the whole bacterial richness, as shown by the rarefaction curves that did not level

150

off even for the samples with the greatest number of pyrotags (Fig. S1).

151

CARD-FISH

152

Water samples were fixed on-board with formaldehyde (1% v/v final concentration) for 24 h

153

at 4 °C. Upon arriving in the laboratory, 80 mL of the water samples were filtered through 0.2

154

µm Isopore™ polycarbonate membrane filters (47 mm diameter; GTTP, Millipore, USA) and

Ac ce

pt

ed

M

an

us

cr

ip t

124

Page 6 of 31

155

stored at -20 °C. CARD-FISH using specific probes (Table 1) was performed according to

156

Pernthaler et al. [46] with a slight modification as described in Korlević et al. [37].

Ac ce

pt

ed

M

an

us

cr

ip t

157

Page 7 of 31

157 Data analyses

159

Observed richness, richness estimators (Chao1 and ACE [abundance-based coverage

160

estimator]) and Shannon’s diversity index were calculated after normalization for sampling

161

effort across samples. A total of 7,843 pyrotags, corresponding to the smallest sampling effort

162

in the dataset, were randomly re-sampled through rarefaction. OTUs that were classified as

163

chloroplasts were not taken into account. Differences in richness or diversity between layers

164

were tested with one-way ANOVA (Systat 12, Systat Software Inc., USA). Normality and

165

homogeneity of variances were tested by Lilliefors and Levene’s tests, respectively.

166

Significant ANOVA results (p<0.05) were followed with post hoc Tukey-HSD multiple

167

comparison tests in order to investigate which of the means were different. Different datasets

168

were constructed in order to estimate the influence of singletons, “No Relative” sequences

169

(sequences that were not successfully taxonomically assigned), and the pooling of sequences

170

at different taxonomic levels (genus-phylum) on bacterial community structure. Datasets were

171

built containing no singletons (OTU-singl.), only taxonomically assigned sequences

172

(OTUannot., without the “No Relative” sequences) and pooled sequences at different

173

taxonomic levels (genus-phylum). In addition, in order to compare 454 pyrosequencing and

174

CARD-FISH, relative abundances of different taxa targeted by the set of probes used in

175

CARD-FISH (Table 1) were extracted from the 454 dataset and compared with the relative

176

abundances of the same taxa detected by CARD-FISH (expressed as a percentage of

177

EUB338I-III signals). Pairwise distance matrices were calculated from the relative abundance

178

data using the Bray-Curtis dissimilarity index [6]. Dissimilarity matrices were compared with

179

Pearson’s product moment correlation coefficient and the significance was determined using

180

the Mantel test followed by the Bonferroni correction. Non-metric multidimensional scaling

181

(NMDS) analysis was performed on a dataset containing relative abundances of each OTU

182

using the Bray-Curtis similarity coefficient in order to compare the samples. All analyses

183

were performed in the R software environment (http://www.r-project.org/) using the vegan

184

package (http://cran.r-project.org/web/packages/vegan/index.html) and custom scripts.

185

To estimate the correlations between community and environmental parameters, Spearman’s

186

correlation coefficients were calculated in order to build the correlation network. Two

187

correlation networks were built: one showing correlations between environmental parameters

188

and bacterial cell counts derived from CARD-FISH, and one showing correlations between

Ac ce

pt

ed

M

an

us

cr

ip t

158

Page 8 of 31

189

environmental parameters and the number of pyrotags belonging to different genera (or

190

different taxa if the genus taxonomic depth could not be achieved) after normalization for

191

sampling. False discovery rates (q values) based on the observed p values were calculated to

192

ensure more stringent criteria. Correlations that matched the criteria of p<0.007, q<0.1 and

193

r>0.65 or r<-0.65 were taken into account. All correlations were performed in R using the

194

vegan

195

Bioconductor

196

(http://www.bioconductor.org/packages/release/bioc/html/qvalue.html).

197

platform Cytoscape was used to visualize the correlation network [54].

198

Results

199

Environmental parameters of the estuary and phytoplankton composition

200

In February, a sharp halocline divided the water column of the estuary forming a strong

201

vertical contrast of physical and chemical conditions (Figs. 2 and S2). Due to the higher

202

freshwater input, the depth of the halocline was deeper than in the July sampling (up to 4 m)

203

and the salinity of the upper water layer ranged from 2.6 to 6.4[CJR5]. The depth of the halocline

204

gradually decreased along the estuary to the reference sea station (AD3) that had no vertical

205

salinity gradient. Surface temperatures inside the estuary ranged from 9.9 °C to 12.4 °C,

206

gradually increasing to 15.0 °C in deeper estuarine layers. Concentrations of total inorganic

207

nitrogen (TIN) and orthosilicate decreased along the estuary, from Station E3 to Station AD3.

208

Concentrations were higher in the surface layer and at the halocline. The highest values of

209

TIN and orthosilicate were detected above the halocline of Station E3 (34.9 and 43.9 µM,

210

respectively). The distribution of orthophosphate did not follow the same pattern. The highest

211

concentration of orthophosphate was detected at the halocline of Station E3 (0.14 µM) and

212

remained highest in the halocline layer. DOC concentrations were uniformly distributed and

213

ranged from 0.94 to 1.38 mg L-1, while POC concentrations were higher in the surface and

214

halocline layer (0.15-0.30 mg L-1) compared to the layer below the halocline or Station AD3

215

(0.05-0.07 mg L-1). A microphytoplankton bloom (2.6 x 105 cells L-1; 2.9 µg L-1 Chl a) was

216

recorded at the surface of Station E3 and it was mostly composed of freshwater diatoms,

217

while at the surface of Station E4a a nanophytoplankton (cryptophytes) bloom was observed

218

(Fig. 6).

219

In July, the vertical gradient of physical and chemical conditions was more pronounced than

220

in February due to strong stratification conditions (Figs. 2 and S2). A sharp halocline layer

Hmisc

(http://cran.r-project.org/web/packages/Hmisc/index.html)

and

ip t

software software

cr

The

Ac ce

pt

ed

M

an

us

package,

Page 9 of 31

occurred at a maximum depth of 2.5 m at Station E3 (minimum salinity 13.8[CJR6]), and

222

gradually decreased along the estuary to the reference sea station (AD3) where the surface

223

salinity was 37.5[CJR7]. A subsurface temperature maximum (23.9 °C) was recorded in the

224

halocline layer of Station E3 probably due to the selective absorption of solar energy by

225

particles along the halocline. The maximum temperature value (24.4 °C) was recorded at the

226

surface of Station E4a, and the minimum (15.7 °C) at the same station at a depth of 30 m. The

227

nutrient concentrations were much lower than in February because of the low precipitation in

228

summer. Maximum concentrations of TIN, orthosilicate and orthophosphate were recorded in

229

the surface layer of Station E3 (10.7, 26.0 and 0.09 µM, respectively). Their concentration

230

decreased downstream, reaching minimum values at Station E5 at a depth of 30 m (1.2, 1.5

231

and 0.02 µM, respectively). In July, the distribution of the DOC and POC concentrations had

232

a similar pattern with higher values at the surface and the halocline. The maximum

233

concentrations of DOC and POC were detected in the halocline layer of Station E3 (1.38 and

234

0.19 mg L-1, respectively). Phytoplankton abundance and composition followed the same

235

trend

236

nanophytoplankton at the surface of Station E3 (1.6 x 105 and 1.3 x 105 cells L-1, respectively;

237

Fig. 6; 2.9 µg L-1 Chl a).

238

In both sampling months it was possible to distinguish three different water layers inside the

239

estuary: the layer above the halocline, the halocline layer and the layer below the halocline.

240

The halocline has been defined as an interval in which ΔS/ΔZ>5 m-1 [39,53,67].

241

Bacterial richness, diversity and community structure

242

Different datasets were compared using Pearson’s correlation coefficient in order to estimate

243

the similarity between community structure at different taxonomic levels, as well as the

244

influence of the removal of singletons or sequences that were not taxonomically assigned. The

245

community structure showed almost no change after removal of singletons or sequences that

246

were not taxonomically assigned (Fig. 3a). In contrast, when the OTU level was compared

247

with higher taxonomic levels (genus-phylum) a drop in the correlation coefficient was

248

observed. Moreover, in order to compare the community structure determined with CARD-

249

FISH and 454 pyrosequencing, relative abundances of different taxonomic groups detected

250

with CARD-FISH were compared with the relative abundances of the same taxonomic groups

251

found in the 454 pyrosequencing data set and with the different taxonomic levels of the 454

252

pyrosequencing data set (Fig. 3b). A weak positive correlation between the CARD-FISH and

nutrients

with

maximum

abundances

of

microphytoplankton

and

M

the

Ac ce

pt

ed

as

an

us

cr

ip t

221

Page 10 of 31

454 pyrosequencing data set was found (R=0.28, p<0.05), indicating that the two methods

254

showed the same direction of community change, but the bacterial structure revealed by the

255

two methods was not the same.

256

Differences in the alpha diversity between layers were described using the total number of

257

bacterial OTUs, Chao1, ACE, and Shannon’s diversity index after a normalization step (Fig

258

S3). In general, small changes in the alpha descriptors between layers and months were

259

detected. In July, significant changes were noticed between layers in the total number of

260

OTUs, Chao1 and ACE estimators and Shannon’s diversity index (one-way ANOVA, p<0.5),

261

while in February no such significant differences were detected. In February, richness values

262

were lower in the halocline layer (Chao1=2,152) compared to the layer below the halocline

263

(Chao1=2,317) and the control Station AD3 (Chao1=2,771) where higher numbers of OTUs

264

were detected. In July, similar numbers of OTUs were observed in the layer above the

265

halocline (Chao1=2,357), the layer below the halocline (Chao1=2,122) and at control Station

266

AD3 (Chao1=2,058) compared to the halocline layer where higher numbers of OTUs were

267

recorded (Chao1=2,576).

268

To identify changes in the bacterial community structure, an NMDS analysis was performed

269

on a table of OTU distributions across samples (Fig. S10). In February, two different

270

communities were identified: the halocline layer community and the community below the

271

halocline layer. In July, three bacterial communities were detected: the halocline layer

272

community, and the community above and below the halocline layer. The grouping of

273

samples in two different communities in February and three communities in July was also

274

supported by ANOSIM (R=0.72, p<0.001). In addition, the differences between the February

275

and July communities were also supported by ANOSIM (R=0.62, p<0.001).

276

Bacterial cell number, diversity and variation

277

In February, a higher average prokaryotic picoplankton cell number was detected in the layer

278

above the halocline (11 x 105 cell mL-1) in comparison with the halocline layer (9.8 x 105 cells

279

mL-1), the layer below the halocline (9.1 x 105 cells mL-1) and Station AD3 (8.4 x 105 cells

280

mL-1; Fig. 6). The maximum abundance (12 x 105 cells mL-1) was detected in the surface

281

layer of Station E4a where the cryptophyte bloom was observed. In July, on average, a higher

282

cell number was characteristic for the layer above the halocline (12 x 105 cell mL-1) and the

283

halocline layer (11 x 105 cells mL-1) compared to the layer below the halocline (8.8 x 105 cells

284

mL-1) and Station AD3 (7.9 x 105 cells mL-1; Fig. 5). The maximum abundance was observed

Ac ce

pt

ed

M

an

us

cr

ip t

253

Page 11 of 31

in the halocline layer of Station E3 (15 x 105 cells mL-1), where a high concentration of

286

organic carbon (POC and DOC) was also observed (Fig. 2). Bacteria dominated the DAPI

287

signals (determined using the EUBI-III probe mix), both in February and July, comprising

288

75% to 94% of the signals (Table S2).

289

The diversity of bacteria was estimated by classifying each reference OTU sequence (using

290

the SILVAngs pipeline, Dataset 1) and by performing the quantitative CARD-FISH analysis

291

using the probes targeting the major taxonomic groups (Table 1). The 454 pyrosequencing

292

was used to obtain an in-depth analysis of diversity, whereas CARD-FISH analysis was used

293

to obtain the cell numbers (absolute and relative).

294

In February and July, the communities were dominated by the alphaproteobacterial clade

295

SAR11 (Figs. 5 and S4). The relative abundance of this oligotrophic clade was higher in the

296

part of the estuary that contained a lower concentration of Chl a, TIN and organic matter

297

(especially POC). Also, a higher abundance was characteristic for July, compared to

298

February. In February, on average, the SAR11 clade comprised a higher proportion of the

299

communities at Station AD3 (41%) and in the layer below the halocline (32%) compared to

300

the halocline layer and the layer above the halocline, where it comprised 22% of the

301

communities. The same distribution was characteristic for July, with a higher contribution at

302

Station AD3 (50%) and the layer below the halocline (48%) compared to the halocline layer

303

(35%) and the layer above the halocline (25%). Another important alphaproteobacterial group

304

was Roseobacter whose abundance was higher, both in February and July, in the layer above

305

the halocline (6% and 9%, respectively) and the halocline layer (14% and 7%, respectively)

306

compared to the other water layers (Figs. 5 and S4). The alphaproteobacterial pyrotag

307

distribution also showed the dominance of SAR11-related pyrotags in both February and July.

308

SAR11 subclade Surface 1 dominated the SAR11-related pyrotags (especially in February),

309

while in June [CJR8]a more diverse distribution of subclades was determined. Pyrotags related to

310

Rhodobacter, Roseibacterium, a Roseobacter clade (NAC11-7 lineage and OCT lineage) and

311

uncultured Rhodobacteraceae had a high contribution that also showed differences in

312

proportion between February and July.

313

Cyanobacteria showed strong change and they were more abundant in July compared to

314

February (Figs. 4 and 5). In February, the maximum abundances were observed at Station

315

AD3 (on average, 17%), while in other layers the relative abundance was <7%. In July, the

316

relative abundance was higher in the halocline layer (16%) and the layer above (11%) and

Ac ce

pt

ed

M

an

us

cr

ip t

285

Page 12 of 31

below (12%) the halocline compared to Station AD3 (10%). The cyanobacterial community

318

was mainly composed of Synechococcus, while Prochlorococcus was a minor constituent, as

319

observed from the distribution of cyanobacterial specific pyrotags (Fig. S5).

320

Betaproteobacteria characterized the riverine layer (above the halocline) with higher relative

321

abundances in February during the higher freshwater input (Figs. 4 and 5). In February and

322

July, the halocline layer contained on average the highest (12% and 7.5%, respectively)

323

abundance of Betaproteobacteria, while in other layers the abundance was <5%. In February,

324

important betaproteobacterial groups, as observed from the pyrotag distribution, were

325

Chlorochromatium and the BAL58 marine group, while in July only BAL58 showed

326

increased pyrotag abundances (Fig. S6). Actinobacteria were also found in the riverine

327

surface waters (the halocline layer and the layer above the halocline) with a pronounced

328

change in abundance between sampling months (Figs. 4 and 5). In February, Actinobacteria

329

comprised, on average, only 3% of the community in the layer above the halocline, while in

330

other layers the relative abundance of this group was <1.5%, as determined by CARD-FISH.

331

In contrast, July was characterized by a high abundance of this group in the halocline layer

332

(8.3%) and the layer above the halocline (14%), while in other layers the average abundance

333

was <1.5%. From the distribution of Actinobacteria-specific pyrotags it could be observed

334

that the high abundance of Actinobacteria in surface waters in July could be attributed to

335

“Candidatus Aquiluna” as it comprised the majority of actinobacterial pyrotags in this water

336

layer (Fig. S7).

337

Bacteroidetes and Gammaproteobacteria (Figs. 4 and 5) showed pronounced differences

338

between February and July. Higher relative abundances were detected in February at Stations

339

E4a and E3 where high Chl a concentrations were also detected. On average, the highest

340

Bacteroidetes relative abundances were detected in February in the halocline layer (21%) and

341

the layer above the halocline (16%) compared to the layer below the halocline (10%) and

342

Station AD3 (3%) where lower values were recorded. In July, lower values were recorded,

343

with minimum values at Station AD3 (7.7%) and in the layer below the halocline (11%)

344

compared to the halocline layer (15%) and the layer above the halocline (13%).

345

Bacteroidetes-specific pyrotags showed a high number of pyrotags related to Owenweeksia,

346

and to marine groups NS4 and NS5 in both February and July (Fig. S8). In February, in

347

addition to these groups, a high number of Flavobacterium- and Formosa-related pyrotags

348

was observed. Gammaproteobacteria showed a similar distribution pattern to Bacteroidetes.

349

On average, the highest relative abundance was detected in February in the halocline layer

Ac ce

pt

ed

M

an

us

cr

ip t

317

Page 13 of 31

(16%) and in the layer below the halocline (17%), while at Station AD3 (13%) and in the

351

layer above the halocline (5.0%) this group comprised a smaller proportion of the prokaryotic

352

picoplankton community. In July, in general, lower values were recorded. The layer above the

353

halocline (13%) and the halocline layer (15%) contained a higher average proportion of

354

Gammaproteobacteria, while in the layer below the halocline (11%) and at Station AD3

355

(7.7%) this proportion was lower. The distribution of gammaproteobacterial pyrotags showed

356

a high number related to the SAR86 clade. In July, in addition to a high number of SAR86

357

related pyrotags, a higher number of Litoricola-related pyrotags was observed (Fig. S9).

358

The relationship between major bacterial groups and environmental factors

359

Association networks using significant correlations (r>0.65 or r<-0.65, p<0.007, q<0.01) were

360

built to examine the interaction between different bacterial groups detected by CARD-FISH

361

(CARD-FISH network) or 454 pyrosequencing (pyrotag network), and physical and chemical

362

environmental parameters (Fig. 6). The pyrotag network was built to examine the interactions

363

at a deeper taxonomic level compared to CARD-FISH, and the environmental parameters

364

were mainly intercorrelated in both networks. Salinity and orthosilicate were at the center of

365

the correlation network emphasizing the estuarial character of the ecosystem. Bacteroidetes

366

were positively correlated with Chl a, POC, DOC and orthosilicate, and negatively correlated

367

with salinity (CARD-FISH network), while groups belonging to this phylum (NS4 and NS5

368

marine groups) were mainly intercorrelated or correlated with other bacterial taxa. Only

369

Owenweeksia was positively correlated with POC, orthosilicate and microphytoplankton.

370

Gammaproteobacteria showed a positive correlation with DOC in the CARD-FISH network,

371

while groups belonging to this class showed a positive correlation with temperature (OM60

372

(NOR5) clade and Litoricola) or positive (SAR86 clade) and negative correlations

373

(Litoricola) with coccolithophorids (<20 µm). Betaproteobacteria were negatively correlated

374

with salinity and positively correlated with TIN, ammonium and nitrate, while no strong

375

significant correlations for different groups that belonged to this class and that accounted for

376

>1% of total pyrotags were found. Actinobacteria showed a negative correlation with

377

orthosilicate (CARD-FISH network), while inside this group only “Candidatus Aquiluna”

378

showed a negative correlation with salinity and a positive correlation with orthosilicate and

379

POC. The SAR11 clade was negatively correlated with Chl a and nitrate (CARD-FISH

380

network), while different SAR11 subclades showed different correlation patterns in the

381

pyrotag network. The Surface 1 subclade was positively correlated with salinity and

382

negatively correlated with orthosilicate, while the Surface 2 subclade was negatively

Ac ce

pt

ed

M

an

us

cr

ip t

350

Page 14 of 31

correlated with TIN, nitrate and ammonium. Roseobacter showed a negative correlation with

384

salinity and positive correlation with DOC, POC, orthosilicate, nitrate and ammonium

385

(CARD-FISH network), while uncultured Rhodobacteraceae showed a negative correlation

386

with salinity and positive correlation with orthosilicate, Chl a, orthophosphate and

387

dinoflagellates (micro fraction) in the pyrotag network. In addition, the pyrotag network

388

showed a positive correlation between the SAR406 clade and salinity, and a negative

389

correlation between SAR116 and Chl a.

390

Discussion

391

Salt-wedge estuaries and coastal regions are contrast ecosystems in which the prokaryotic

392

picoplankton community shows strong spatial and temporal variations, due to high salinity

393

and nutrient gradients. This study describes the spatial bacterial diversity in the Krka River

394

estuary based on the application of two techniques, 454 pyrosequencing of the 16S rRNA

395

gene and CARD-FISH in two environmentally contrasting months. As a result, weak positive

396

correlation between the bacterial community structure determined with CARD-FISH and 454

397

pyrosequencing was observed. The discrepancy between the two methods could arise from

398

DNA extraction efficiency, bias introduced by PCR, clade specific differences in the rRNA

399

operon number, probe specificity or the physiological status of the cells. A similar correlation

400

coefficient between the two methods was found earlier [37].

401

The identified differences between the bacterial communities could be attributed to different

402

environmental conditions in February and July. A number of studies have described season as

403

the main factor shaping bacterial community structure [21,24,25,33,34]. Seasonal change

404

could also be the main factor shaping the bacterial community dynamics in the Krka River

405

estuary but no final conclusion could be drawn from two samplings. Studies describing

406

estuarine bacterial communities have not given uniform results for the main factors

407

influencing the spatial and temporal dynamics [20]. It was also hypothesized that, in the case

408

of small spatial variations, which do not include more than one environment (coastal waters,

409

offshore waters, estuaries etc.), as in the Krka estuary, the dominant component should be

410

seasonal change, while in the case of more environments, the spatial component should be

411

dominant [20]. The probable cause for the formation of different communities in each layer is

412

the sharp halocline on which the formation of a film of organic matter has been described.

413

This film contributes to the stabilization of the density gradient, affects the energy- and mass-

414

transport processes and influences the transformation of particles and solutes [71]. In addition,

Ac ce

pt

ed

M

an

us

cr

ip t

383

Page 15 of 31

the domination of temporal over spatial changes was probably caused by the difference in the

416

hydrographic conditions between February and July. In February, the higher freshwater input

417

caused a higher concentration of nutrients and autotrophic biomass (Chl a and phytoplankton

418

abundances) that subsequently resulted in the accumulation of fresh organic matter on the

419

halocline and a change in the bacterial community structure. The organic matter accumulated

420

on the halocline was produced, at least partially, by the freshwater phytoplankton developing

421

in the freshwater lake located above the estuary (Lake Visovac) [66]. Moreover, a domination

422

of seasonal influences over vertical, spatial changes for the phytoplankton community was

423

observed previously in the Krka River estuary [13].

424

The correlation networks built using the CARD-FISH and 454 pyrosequencing data enabled

425

the influence of environmental parameters on different bacterial groups to be estimated. The

426

interactions were stronger within the environmental parameters and bacterial groups than

427

between these groups, as was shown previously [25]. Interestingly, a stronger interaction

428

within bacterial groups (especially within Bacteroidetes and Gammaproteobacteria) was

429

observed in the network built on the 454 pyrosequencing data. The correlation of different

430

taxonomic groups (at different taxonomic levels of resolution) did not show a clear positive

431

correlation with different phytoplankton groups or Chl a (except for Bacteroidetes, uncultured

432

Rhodobacteraceae and SAR86). Therefore, we hypothesize that, apart from no clear

433

correlation, a phytoplankton bloom in February probably caused the changes in the bacterial

434

community that was most likely blurred by having only two sampling points and by a partial

435

production of organic matter in the freshwater lake above the estuary. In addition, it has been

436

found that the distribution of some bacterial groups depended on the properties of organic

437

matter, which could also explain the lack of correlation between some bacterial groups and

438

the quantity of organic matter in our study [2]. Interestingly, a higher number of bacterial

439

groups was positively correlated with orthosilicate, indicating an indirect link between

440

bacteria and phytoplankton (diatoms). OTU richness in the Krka estuary showed differences

441

between the two samplings, which were particularly notable between the halocline layers in

442

February and July. The lower richness observed in the halocline layer in February could be

443

attributed to the availability of new ecological niches, due to the winter phytoplankton bloom,

444

in which specialist groups can appear, lowering the overall richness [37,61]. A similar pattern

445

of DGGE bands, with the lowest band richness in February and highest in August and

446

October, was identified in Chesapeake Bay [33].

Ac ce

pt

ed

M

an

us

cr

ip t

415

Page 16 of 31

The dynamic of total prokaryotic picoplankton abundance changed according to previously

448

described patterns [22,23]. Bacteria dominated the prokaryotic picoplankton communities in

449

the Krka estuary with the SAR11 clade being the most abundant group, especially at marine

450

Station AD3 and in the deep marine layer where lower concentrations of organic matter were

451

detected (mainly POC) [11,44]. In addition, the contribution of the SAR11 clade to the

452

prokaryotic community was higher in July, when generally lower concentrations of organic

453

matter were detected. The worldwide distribution of this clade [11,44], supported by genomic

454

data, suggests that its members play a major role in the oxidation of low-molecular-weight

455

organic matter [26,51]. Subclade Surface 1 dominated the SAR11 population, which was in

456

contrast to the results from the Baltic Sea [29], while a more diverse SAR11 subpopulation

457

occurred during the summer. In addition, different SAR11 subclades showed distinct

458

correlation patterns, indicating possible specific subclade physiological adaptations [11].

459

Roseobacter was more abundant in riverine waters and in the halocline layers characterized

460

by higher concentrations of Chl a and organic matter [19,62]. Members of this group were

461

found also in other estuaries but their relative abundance in the community was lower

462

[15,33,36]. Also, a high abundance of Roseobacter during phytoplankton biomass increase

463

was observed earlier in accordance with data showing the important role of this group in the

464

transformation of sulfur compounds [7,19,27,62]. Cyanobacteria abundance peaked in July,

465

when nutrient concentrations were low. Due to their higher surface to volume ratio they are

466

considered more competitive in oligotrophic environments compared to other phytoplankton

467

groups. The most common cyanobacterial pyrotags were related to Synechococcus in

468

accordance with previous data for the Krka River estuary derived using flow cytometry [59].

469

Higher actinobacterial and betaproteobacterial abundance in the low-salinity riverine water

470

and a decrease towards the mouth of the estuary and high-salinity waters was also observed in

471

other estuaries studied [5,36,70]. Actinobacteria displayed a strong temporal pattern with

472

higher abundances in July compared to February [3]. In the same month, a high proportion of

473

actinobacterial pyrotags was related to “Candidatus Aquiluna”. The genome of one strain

474

belonging to this genus was sequenced and the proteorhodopsin gene was identified, which

475

indicated a possible adaptive strategy to the more oligotrophic conditions characteristic for

476

July [35]. Betaproteobacteria demonstrated lower differences between the sampling months

477

than Actinobacteria. In surface waters of the Krka estuary a high proportion of

478

betaproteobacterial pyrotags was related to the BAL58 marine group, representing an

Ac ce

pt

ed

M

an

us

cr

ip t

447

Page 17 of 31

oligotrophic marine bacterium that was isolated by dilution to extinction from the Baltic Sea

480

[55], which is often identified in post and phytoplankton blooms [68,69].

481

The presence of Bacteroidetes and Gammaproteobacteria was characteristic for waters with

482

greater concentrations of organic matter, which was also visible in the network analysis. An

483

increase in abundance of these two groups was observed in the North Sea during

484

phytoplankton bloom and postbloom periods [56,61]. Pyrotags related to the marine groups

485

NS3a, NS4 and NS5, as well as Flavobacterium and Formosa, comprised the majority of the

486

Bacteroidetes-specific pyrotags. A high abundance of Formosa during a diatom bloom was

487

described earlier in the coastal North Sea [61]. Other groups identified during the diatom

488

bloom in the North Sea, such as Polaribacter and Ulvibacter, have not been detected in the

489

Krka estuary [61]. The same groups were also not found in a previous study of the South

490

Adriatic bacterial communities, indicating possible biogeographic patterns [37]. SAR86

491

clade-related pyrotags were the most abundant gammaproteobacterial pyrotags in the study

492

area. A high proportion of SAR86-related sequences were also found in the South and North

493

Adriatic [37,63]. Higher proportions of Litoricola-related pyrotags were also identified in the

494

South Adriatic [37]. Genomic studies on a cultivated strain of this genus identified the

495

proteorodopsin gene and the possibility for CO2 fixation that could be an adaptive strategy for

496

the more oligotrophic conditions present in July [30].

497

The present study provided a detailed comparison of microbial communities in a

498

Mediterranean karstic estuary using high throughput sequencing and CARD-FISH.

499

Remarkable differences in this small area were observed with clear shifts in the population

500

structure between the two months studied and along the vertical salinity gradient. However,

501

more time points would be needed to disentangle the specific effects of individual

502

environmental parameters on the community. The phytoplankton bloom observed in February

503

was probably the main source of temporal change in the bacterial structure. Roseobacter,

504

Bacteroidetes and Gammaproteobacteria, known to respond to phytoplankton blooms,

505

showed pronounced temporal differences, while Actinobacteria and Betaproteobacteria,

506

known to be characteristic for freshwater and brackish water, showed strong vertical

507

distribution patterns in relation to the salinity gradient. In addition, Cyanobacteria decreased

508

in abundance in February as a response to competition from phytoplankton, while the SAR11

509

clade displayed an increase in abundance in July as a result of a better adaptation towards

510

more oligotrophic conditions. In conclusion, different bacterial phylogenetic groups

511

responded differently to temporal and spatial changes in the Krka River estuary. However,

Ac ce

pt

ed

M

an

us

cr

ip t

479

Page 18 of 31

future research will be needed in order to identify the key processes present in this

513

environment.

514

Acknowledgements

515

The authors would like to thank the crew and scientific staff on board the R/V Hidra of the

516

Hydrographic Institute of the Republic of Croatia for their help during the fieldwork, as well

517

as to Goran Olujić for nutrient analyses. Also, we would like to thank Bernhard Fuchs from

518

the Max Planck Institute for Marine Microbiology for critical reading of the manuscript.

519

Financial support was provided by the Croatian Science Foundation through the BABAS

520

project (to SO and MK) and partially by the IP_11_2013_1205 SPHERE project (to IC and

521

JD). Further support was provided by the Research Council of Norway (FRIMEDBIO project

522

197823) and the Royal Swedish Academy of Sciences through a grant from the Knut and

523

Alice Wallenberg Foundation (KAW 2009.0287).

524

References [CJR9] [CJR10]

525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552

[1]

[4] [5] [6] [7] [8]

[9] [10]

an

M

ed

pt

[3]

Alonso, C., Gómez-Pereira, P., Ramette, A., Ortega, L., Fuchs, B., Amann, R. (2010) Multilevel analysis of the bacterial diversity along the environmental gradient Río de la Plata–South Atlantic Ocean. Aquat. Microb. Ecol. 61(1), 57–72, Doi: 10.3354/ame01439. Amaral, V., Graeber, D., Calliari, D., Alonso, C. (2016) Strong linkages between DOM optical properties and main clades of aquatic bacteria. Limnol. Oceanogr. 61(3), 906– 918, Doi: 10.1002/lno.10258. Andersson, A.F., Riemann, L., Bertilsson, S. (2010) Pyrosequencing reveals contrasting seasonal dynamics of taxa within Baltic Sea bacterioplankton communities. ISME J. 4(2), 171–181, Doi: 10.1038/ismej.2009.108. Barlow, R., Cummings, D., Gibb, S. (1997) Improved resolution of mono- and divinyl chlorophylls a and b and zeaxanthin and lutein in phytoplankton extracts using reverse phase C-8 HPLC. Mar. Ecol. Prog. Ser. 161, 303–307, Doi: 10.3354/meps161303. Bouvier, T.C., del Giorgio, P.A. (2002) Compositional changes in free-living bacterial communities along a salinity gradient in two temperate estuaries. Limnol. Oceanogr. 47(2), 453–470, Doi: 10.4319/lo.2002.47.2.0453. Bray, J.R., Curtis, J.T. (1957) An ordination of the upland forest communities of Southern Wisconsin. Ecol. Monogr. 27(4), 325-349, Doi: 10.2307/1942268. Buchan, A., Moran, M.A. (2005) Overview of the marine Roseobacter lineage. Appl. Environ. Microbiol. 71(10), 5665–5677, Doi: 10.1128/AEM.71.10.5665. Bužančić, M., Ninčević Gladan, Ž., Marasović, I., Kušpilić, G., Grebec, B., Matijević, S. (2012) Population structure and abundance of phytoplankton in three bays on the eastern Adriatic coast: Šibenik Bay, Kaštela Bay and Mali Ston Bay. Acta Adriat. 53(3), 413–435. Camacho, C., Coulouris, G., Avagyan, V., Ma, N., Papadopoulos, J., Bealer, K., Madden, T.L. (2009) BLAST+: architecture and applications. BMC Bioinformatics 10, 421, Doi: 10.1186/1471-2105-10-421. Campbell, B.J., Kirchman, D.L. (2013) Bacterial diversity, community structure and

Ac ce

[2]

us

cr

ip t

512

Page 19 of 31

[15]

[16]

[19]

[20] [21]

[22] [23] [24]

[25]

pt

[18]

ed

M

[17]

ip t

[14]

cr

[13]

us

[12]

an

[11]

potential growth rates along an estuarine salinity gradient. ISME J. 7, 210–220, Doi: 10.1038/ismej.2012.93. Carlson, C.A., Morris, R., Parsons, R., Treusch, A.H., Giovannoni, S.J., Vergin, K. (2009) Seasonal dynamics of SAR11 populations in the euphotic and mesopelagic zones of the northwestern Sargasso Sea. ISME J. 3(3), 283–295, Doi: 10.1038/ismej.2008.117. Cauwet, G. (1991) Carbon inputs and biogeochemical processes at the halocline in a stratified estuary: Krka River, Yugoslavia. Mar. Chem. 32(2-4), 269–283, Doi: 10.1016/0304-4203(91)90043-V. Cetinić, I., Viličić, D., Burić, Z., Olujić, G. (2006) Phytoplankton seasonality in a highly stratified karstic estuary (Krka, Adriatic Sea). Hydrobiologia 555(1), 31–40, Doi: 10.1007/s10750-005-1103-7. Crump, B.C., Armbrust, E.V., Baross, J.A. (1999) Phylogenetic analysis of particleattached and free-living bacterial communities in the Columbia River, its estuary, and the adjacent coastal ocean. Appl. Envir. Microbiol. 65(7), 3192–3204. Crump, B.C., Hopkinson, C.S., Sogin, M.L., Hobbie, J.E. (2004) Microbial biogeography along an estuarine salinity gradient: combined influences of bacterial growth and residence time. Appl. Environ. Microbiol. 70(3), 1494–1505, Doi: 10.1128/AEM.70.3.1494-1505.2004. Dafner, E.V., Wangersky, P.J. (2002) A brief overview of modern directions in marine DOC studies. Part I: Methodological aspects. J. Environ. Monit. 4(1), 48–54, Doi: 10.1039/b107277n. Dautović, J., Strmečki, S., Pestorić, B., Vojvodić, V., Plavšić, M., Krivokapić, S., Ćosović, B. (2012) Organic matter in the karstic enclosed bay (Boka Kotorska Bay, South Adriatic Sea). Influence of freshwater input. Fresenius Enviromental Bull. 21, 995–1006. Dyer, K.R. (1991) Circulation and mixing in stratified estuaries. Mar. Chem. 32, 111– 120. Eilers, H., Pernthaler, J., Peplies, J., Glöckner, F.O., Gerdts, G., Amann, R. (2001) Isolation of novel pelagic bacteria from the German Bight and their seasonal contributions to surface picoplankton. Appl. Environ. Microbiol. 67(11), 5134–5142, Doi: 10.1128/AEM.67.11.5134-5142.2001. Fortunato, C.S., Herfort, L., Zuber, P., Baptista, A.M., Crump, B.C. (2012) Spatial variability overwhelms seasonal patterns in bacterioplankton communities across a river to ocean gradient. ISME J. 6(3), 554–563, Doi: 10.1038/ismej.2011.135. Fuhrman, J.A., Hewson, I., Schwalbach, M.S., Steele, J.A., Brown, M.V., Naeem, S. (2006) Annually reoccurring bacterial communities are predictable from ocean conditions. Proc. Natl. Acad. Sci. 103(35), 13104–13109, Doi: 10.1073/pnas.0602399103. Fuks, D., Devescovi, M., Precali, R., Krstulović, N., Šolić, M. (1991) Bacterial abundance and activity in the highly stratified estuary of the Krka River. Mar. Chem. 32, 333–346. Fuks, D., Precali, R., Devescovi, M. (1994) Bacterial production in the stratified karstic estuary of the Krka River. Acta Adriat. 34, 21–28. Gilbert, J.A., Field, D., Swift, P., Newbold, L., Oliver, A., Smyth, T., Somerfield, P.J., Huse, S., Joint, I. (2009) The seasonal structure of microbial communities in the Western English Channel. Environ. Microbiol. 11, 3132–3139, Doi: 10.1111/j.14622920.2009.02017.x. Gilbert, J.A., Steele, J.A., Caporaso, J.G., Steinbrück, L., Reeder, J., Temperton, B., Huse, S., McHardy, A.C., Knight, R., Joint, I., Somerfield, P., Fuhrman, J.A., Field, D.

Ac ce

553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602

Page 20 of 31

[31]

[32] [33] [34] [35]

[36] [37] [38] [39] [40] [41]

ip t

cr

us

[30]

an

[29]

M

[28]

ed

[27]

pt

[26]

(2012) Defining seasonal marine microbial community dynamics. ISME J. 6(2), 298– 308, Doi: 10.1038/ismej.2011.107. Giovannoni, S.J., Tripp, H.J., Givan, S., Podar, M., Vergin, K.L., Baptista, D., Bibbs, L., Eads, J., Richardson, T.H., Noordewier, M., Rappé, M.S., Short, J.M., Carrington, J.C., Mathur, E.J. (2005) Genome streamlining in a cosmopolitan oceanic bacterium. Science 309(5738), 1242–1245, Doi: 10.1126/science.1114057. Gonzalez, J.M. (2003) Silicibacter pomeroyi sp. nov. and Roseovarius nubinhibens sp. nov., dimethylsulfoniopropionate-demethylating bacteria from marine environments. Int. J. Syst. Evol. Microbiol. 53(5), 1261–1269, Doi: 10.1099/ijs.0.02491-0. Gržetić, Z., Precali, R., Degobbis, D., Škrivanić, A. (1991) Nutrient enrichment and phytoplankton response in an Adriatic karstic estuary. Mar. Chem. 32, 313–331. Herlemann, D.P.R., Woelk, J., Labrenz, M., Jürgens, K. (2014) Diversity and abundance of “Pelagibacterales” (SAR11) in the Baltic Sea salinity gradient. Syst. Appl. Microbiol. 37(8), 601–604, Doi: 10.1016/j.syapm.2014.09.002. Huggett, M.J., Rappe, M.S. (2012) Genome sequence of strain HIMB30, a novel member of the marine Gammaproteobacteria. J. Bacteriol. 194(3), 732–733, Doi: 10.1128/JB.06506-11. Ionescu, D., Siebert, C., Polerecky, L., Munwes, Y.Y., Lott, C., Häusler, S., BižićIonescu, M., Quast, C., Peplies, J., Glöckner, F.O., Ramette, A., Rödiger, T., Dittmar, T., Oren, A., Geyer, S., Stärk, H.-J., Sauter, M., Licha, T., Laronne, J.B., de Beer, D. (2012) Microbial and chemical characterization of underwater fresh water springs in the Dead Sea. PLoS One 7(6), e38319, Doi: 10.1371/journal.pone.0038319. Ivančić, I., Degobbis, D. (1984) An optimal manual procedure for ammonia analysis in natural waters by the indophenol blue method. Water Res. 18(9), 1143–1147, Doi: 10.1016/0043-1354(84)90230-6. Kan, J., Crump, B.C., Wang, K., Chen, F. (2006) Bacterioplankton community in Chesapeake Bay: Predictable or random assemblages. Limnol. Oceanogr. 51(5), 2157– 2169, Doi: 10.4319/lo.2006.51.5.2157. Kan, J., Suzuki, M.T., Wang, K., Evans, S.E., Chen, F. (2007) High temporal but low spatial heterogeneity of bacterioplankton in the Chesapeake Bay. Appl. Environ. Microbiol. 73(21), 6776–6789, Doi: 10.1128/AEM.00541-07. Kang, I., Lee, K., Yang, S.J., Choi, A., Kang, D., Lee, Y.K., Cho, J.C. (2012) Genome sequence of “Candidatus Aquiluna” sp. strain IMCC13023, a marine member of the Actinobacteria isolated from an arctic fjord. J. Bacteriol. 194(13), 3550–3551, Doi: 10.1128/JB.00580-12. Kirchman, D.L., Dittel, A.I., Malmstrom, R.R., Cottrell, M.T. (2005) Biogeography of major bacterial groups in the Delaware Estuary. Limnol. Oceanogr. 50(5), 1697–1706, Doi: 10.4319/lo.2005.50.5.1697. Korlević, M., Pop Ristova, P., Garić, R., Amann, R., Orlić, S. (2015) Bacterial diversity in the South Adriatic Sea during a strong, deep winter convection year. Appl. Environ. Microbiol. 81(5), 1715–1726, Doi: 10.1128/AEM.03410-14. Legović, T. (1991) Exchange of water in a stratified estuary with an application to Krka (Adriatic Sea). Mar. Chem. 32(2-4), 121–135, Doi: 10.1016/03044203(91)90032-R. Legović, T., Gržetić, Z., Smirčić, A. (1991) Effects of wind on a stratified estuary. Mar. Chem. 32, 153–161. Legović, T., Gržetić, Z., Žutić, V. (1991) Subsurface temperature maximum in a stratified estuary. Mar. Chem. 32, 163–170. Li, W., Godzik, A. (2006) Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics 22(13), 1658–1659, Doi:

Ac ce

603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652

Page 21 of 31

[47] [48] [49] [50]

[51] [52]

[53] [54]

[55] [56] [57]

ip t

cr

us

[46]

an

[45]

M

[44]

ed

[43]

10.1093/bioinformatics/btl158. Massana, R., Murray, A.E., Preston, C.M., Delong, E.F. (1997) Vertical distribution and phylogenetic characterization of marine planktonic Archaea in the Santa Barbara Channel. Appl. Environ. Microbiol. 63(1), 50–56. Moreira-Turcq, P., Martin, J.M., Fleury, A. (1993) Chemical and biological characterization of particles by flow cytometry in the Krka estuary, Croatia. Mar. Chem. 43(1-4), 115–126, Doi: 10.1016/0304-4203(93)90219-E. Morris, R.M., Rappé, M.S., Connon, S.A., Vergin, K.L., Siebold, W.A., Carlson, C.A., Giovannoni, S.J. (2002) SAR11 clade dominates ocean surface bacterioplankton communities. Nature 420(6917), 806–810, Doi: 10.1038/nature01240. Penezić, A., Gašparović, B., Burić, Z., Frka, S. (2010) Distribution of marine lipid classes in salty Rogoznica Lake (Croatia). Estuar. Coast. Shelf Sci. 86(4), 625–636, Doi: 10.1016/j.ecss.2009.11.030. Pernthaler, A., Pernthaler, J., Amann, R. (2002) Fluorescence in situ hybridization and catalyzed reporter deposition for the identification of marine bacteria. Appl. Environ. Microbiol. 68(6), 3094–3101. Preen, K., Kirchman, D.L. (2004) Microbial respiration and production in the Delaware Estuary. Aquat. Microb. Ecol. 37(2), 109–119, Doi: 10.3354/ame037109. Pruesse, E., Peplies, J., Glöckner, F.O. (2012) SINA: accurate high-throughput multiple sequence alignment of ribosomal RNA genes. Bioinformatics 28(14), 1823–1829, Doi: 10.1093/bioinformatics/bts252. Puddu, A., Ferla, R.L., Allegra, A., Bacci, C., Lopez, M., Oliva, F., Pierotti, C. (1998) Seasonal and spatial distribution of bacterial production and biomass along a salinity gradient (northern Adriatic Sea). Hydrobiologia 363(1-3), 271–282. Quast, C., Pruesse, E., Yilmaz, P., Gerken, J., Schweer, T., Yarza, P., Peplies, J., Glöckner, F.O. (2013) The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 41(Database issue), D590– 596, Doi: 10.1093/nar/gks1219. Rappé, M.S., Connon, S.A., Vergin, K.L., Giovanonni, S.J. (2002) Cultivation of the ubiquitous SAR11 marine bacterioplankton clade. Nature 418(August), 630–633. Schloss, P.D., Westcott, S.L., Ryabin, T., Hall, J.R., Hartmann, M., Hollister, E.B., Lesniewski, R.A., Oakley, B.B., Parks, D.H., Robinson, C.J., Sahl, J.W., Stres, B., Thallinger, G.G., Van Horn, D.J., Weber, C.F. (2009) Introducing mothur: opensource, platform-independent, community-supported software for describing and comparing microbial communities. Appl. Environ. Microbiol. 75(23), 7537–7541, Doi: 10.1128/AEM.01541-09. Sempéré, R., Cauwet, G. (1995) Occurrence of organic colloids in the stratified estuary of the Krka River (Croatia). Estuar. Coast. Shelf Sci. 40, 105–114. Shannon, P., Markiel, A., Ozier, O., Baliga, N.S., Wang, J.T., Ramage, D., Amin, N., Schwikowski, B., Ideker, T. (2003) Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13(11), 2498–2504, Doi: 10.1101/gr.1239303. Simu, K., Hagström, Å. (2004) Oligotrophic bacterioplankton with a novel single-cell life strategy. Appl. Environ. Microbiol. 70(4), 2445–2451, Doi: 10.1128/AEM.70.4.2445-2451.2004. Sintes, E., Witte, H., Stodderegger, K., Steiner, P., Herndl, G.J. (2013) Temporal dynamics in the free-living bacterial community composition in the coastal North Sea. FEMS Microbiol. Ecol. 83(2), 413–424, Doi: 10.1111/1574-6941.12003. Strickland, J.D.H., Parsons, T.R. (1972) A practical handbook of seawater analysis. Bull. Fish. Res. Board Can. 167, 310 pp, Doi: 10.1002/iroh.19700550118.

pt

[42]

Ac ce

653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702

Page 22 of 31

[59]

[60]

[65] [66] [67] [68]

[69] [70] [71]

an

M

[64]

ed

[63]

pt

[62]

us

cr

[61]

Šantić, D., Krstulović, N., Šolić, M., Kušpilić, G. (2011) Distribution of Synechococcus and Prochlorococcus in the central Adriatic Sea. Acta Adriat. 52(1), 101–114. Šolić, M., Krstulović, N., Šantić, D., Šestanović, S., Ordulj, M., Bojanić, N., Kušpilić, G. (2015) Structure of microbial communities in phosphorus-limited estuaries along the eastern Adriatic coast. J. Mar. Biol. Assoc. United Kingdom, 1–14, Doi: 10.1017/S0025315415000442. Šupraha, L., Ljubeši, Z., Mihanovi, H., Oluji, G., Mikac, I., Vili, D. (2014) Cryptophyte bloom in a Mediterranean estuary: High abundance of Plagioselmis cf . prolonga in the Krka River estuary (eastern Adriatic Sea). Sci. Mar. 78(3), 329–338. Teeling, H., Fuchs, B.M., Becher, D., Klockow, C., Gardebrecht, A., Bennke, C.M., Kassabgy, M., Huang, S., Mann, A.J., Waldmann, J., Weber, M., Klindworth, A., Otto, A., Lange, J., Bernhardt, J., Reinsch, C., Hecker, M., Peplies, J., Bockelmann, F.D., Callies, U., Gerdts, G., Wichels, A., Wiltshire, K.H., Glöckner, F.O., Schweder, T., Amann, R. (2012) Substrate-controlled succession of marine bacterioplankton populations induced by a phytoplankton bloom. Science 336(6081), 608–611, Doi: 10.1126/science.1218344. Thiele, S., Fuchs, B.M., Ramaiah, N., Amann, R. (2012) Microbial community response during the iron fertilization experiment LOHAFEX. Appl. Environ. Microbiol. 78(24), 8803–8812, Doi: 10.1128/AEM.01814-12. Tinta, T., Vojvoda, J., Mozetič, P., Talaber, I., Vodopivec, M., Malfatti, F., Turk, V. (2014) Bacterial community shift is induced by dynamic environmental parameters in a changing coastal ecosystem (northern Adriatic, northeastern Mediterranean Sea) – a 2year time-series study. Environ. Microbiol. 17(10), 3581–3596, Doi: 10.1111/14622920.12519. Troussellier, M., Schäfer, H., Batailler, N., Bernard, L., Courties, C., Lebaron, P., Muyzer, G., Servais, P., Vives-Rego, J. (2002) Bacterial activity and genetic richness along an estuarine gradient (Rhone River plume, France). Aquat. Microb. Ecol. 28(1), 13–24, Doi: 10.3354/ame028013. Utermöhl, H. (1958) Zur Vervollkommnung der quantitativen PhytoplanktonMethodik. Mitteilung Int. Vereinigung Fuer Theor. Und Amgewandte Limnol. 9, 1–38. Viličić, D., Legović, T., Žutić, V. (1989) Vertical distribution of phytoplankton in a stratified estuary. Aquat. Sci. 51(1), 31–46, Doi: 10.1007/BF00877779. Vukojević, N. (1983) Effect of thermocline on sound velocity in the sea. University of Zagreb, 1983. Wemheuer, B., Güllert, S., Billerbeck, S., Giebel, H.-A., Voget, S., Simon, M., Daniel, R. (2014) Impact of a phytoplankton bloom on the diversity of the active bacterial community in the southern North Sea as revealed by metatranscriptomic approaches. FEMS Microbiol. Ecol. 87(2), 378–389, Doi: 10.1111/1574-6941.12230. Yeo, S.K., Huggett, M.J., Eiler, A., Rappé, M.S. (2013) Coastal bacterioplankton community dynamics in response to a natural disturbance. PLoS One 8(2), e56207, Doi: 10.1371/journal.pone.0056207. Zhang, Y., Jiao, N., Cottrell, M., Kirchman, D. (2006) Contribution of major bacterial groups to bacterial biomass production along a salinity gradient in the South China Sea. Aquat. Microb. Ecol. 43, 233–241, Doi: 10.3354/ame043233. Žutić, V., Legović, T. (1987) A film of organic matter at the fresh-water/sea-water interface of an estuary. Nature 328, 612–614.

ip t

[58]

Ac ce

703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750

Page 23 of 31

Figure Legends

751 752

Fig. 1. Study area and sampling stations inside (E3, E4a and E5) and outside (AD3) the Krka estuary.

753 754 755

Fig. 2. Vertical distributions of temperature (T), salinity (S), chlorophyll a (Chl a), total inorganic nitrogen (TIN), particulate organic carbon (POC) and dissolved organic carbon (DOC) concentrations at stations AD3, E5, E4a and E3 in February and July 2013.

756 757 758 759 760 761 762 763 764 765

Fig. 3. Bacterial community structures were compared after the removal of singletons (OUT[CJR11]-singl.) and “No Relative” sequences (OTUannot.), and after the pooling of sequences at different taxonomic levels (a) using the Pearson’s correlation coefficient. Taxa detected by CARD-FISH were compared with the same taxa from the 454 dataset (labeled as 454, red square) and with the different taxonomic levels of the 454 dataset (b) using the Pearson’s correlation coefficient. The correlation coefficient calculated from the distance matrices resulted in relative abundances of sequences or CARD-FISH counts (expressed as a percentage of the EUBI-III signals). Correlation significances were calculated by the Mantel test with 1,000 matrix permutations. All R values were significant (<0.05) after Bonferroni correction.

766 767

Fig. 4. Taxonomic classification and relative contribution of the most abundant bacterial pyrotags (>2%) in February and July 2013.

768 769 770 771

Fig. 5. Vertical change of total picoplankton cells (DAPI counts), prokaryotic picoplankton community structure determined by CARD-FISH, nanophytoplankton and microphytoplankton cells at stations AD3, E5, E4a and E3 in February and July 2013. A vertical salinity gradient is represented for each station and month.

772 773 774 775 776 777 778 779 780 781 782

Fig. 6. [CJR12]The association network showing significant correlations (Spearman’s correlation coefficient) between environmental parameters and bacterial groups (detected by CARDFISH) (a), and environmental parameters and pyrotag numbers belonging to different genera (or different taxa if genus taxonomic depth could not be achieved) (b). Only significant correlations (r>0.65 or r<-0.65, p<0.007, q<0.01) were plotted. For the 454 data, only taxa that accounted for more than 1% of the total pyrotags are shown. The correlation coefficients (r) are labeled on the lines connecting the nodes. Solid line-positive correlation, dashed linenegative correlation, grey nodes-environmental parameters, blue nodes-bacterial groups, green nodes-phytoplankton groups, T-temperature, S-salinity, TIN-total inorganic nitrogen, NO3-nitrate, NO2-nitrite, NH4-ammonium, PO4-orthophosphate, SiO4-orthosilicate[CJR13], Chl a-chlorophyll a, POC-particulate organic carbon, and DOC-dissolved organic carbon.

cr

us

an

M

ed

pt

Ac ce

783

ip t

750

Page 24 of 31

Table 1. Selected probes and hybridization conditions used for CARD-FISH. Probe

Target organisms

Sequence (5’→3’)

FA (%)a

EUB338

Bacteria

GCTGCCTCCCGTAGGAGT

35

(Amann et al. 1990)

EUB338-II

Bacteria (supplement to EUB338)

GCAGCCACCCGTAGGTGT

35

(Daims et al. 1999)

EUB338-III

Bacteria (supplement to EUB338)

GCTGCCACCCGTAGGTGT

35

(Daims et al. 1999)

NON338

Control

ACTCCTACGGGAGGCAGC

HGC69a

Actinobacteria

TATAGTTACCACCGCCGT

Reference

CF319a

Bacteroidetes

CYA664

ip t

[CJR14]

cr

35

(Wallner, Amann and Beisker 1993) (Roller et al. 1994)

TGGTCCGTGTCTCAGTAC

35

(Manz et al. 1996)

Cyanobacteria

GGAATTCCCTCTGCCCC

35

(Schönhuber et al. 1999)

ALF968

Alphaproteobacteria

GGTAAGGTTCTGCGCGTT

35

(Neef 1997)

ROS537

Roseobacter

35

(Eilers et al. 2001) (Morris et al. 2002)

ed

ATTAGCACAAGTTTCCYCGTGT

SAR11441R(ori)b

(Morris et al. 2002)

TACAGTCATTTTCTTCCCCGAC

pt

SAR11487(mod)b

an

CAACGCTAACCCCCTCC

SAR11-152Rb

SAR11441R(mod)b

us

25

M

783

TACCGTCATTTTCTTCCCCGAC

SAR11 clade

25

(Morris et al. 2002) (Schattenhofer et al. 2009)

SAR11-542Rb

TCCGAACTACGCTAGGTC

(Morris et al. 2002)

SAR11-732Rb

GTCAGTAATGATCCAGAAAGYTG

(Morris et al. 2002)

Ac ce

CGGACCTTCTTATTCGGG

BET42ac

Betaproteobacteria

GCCTTCCCACTTCGTTT

35

(Manz et al. 1992)

GAM42ad

Gammaproteobacteria

GCCTTCCCACATCGTTT

35

(Manz et al. 1992)

784a. Formamide concentration (v/v) in CARD-FISH hybridization buffer. 785b. A mixture of six probes to detect the SAR11 clade, including an unlabelled helper SAR11-487-h3 (5′786 CGGCTGCTGGCACGAAGTTAGC-3′). 787c. Including an unlabelled competitor probe Gam42a (5′-GCCTTCCCACATCGTTT-3′) (Manz et al. 1992).

Page 25 of 31

788d. Including an unlabelled competitor probe Bet42a (5′-GCCTTCCCACTTCGTTT-3′) (Manz et al. 1992). 789

Ac ce

pt

ed

M

an

us

cr

ip t

790

Page 26 of 31

Ac

ce

pt

ed

M

an

us

cr

i

Figure 1

Page 27 of 31

Ac

ce

pt

ed

M

an

us

cr

i

Figure 2

Page 28 of 31

Ac ce p

te

d

M

an

us

cr

ip t

Figure 3

Page 29 of 31

Ac

ce

pt

ed

M

an

us

cr

i

Figure 4

Page 30 of 31

Ac

ce

pt

ed

M

an

us

cr

i

Figure 5

Page 31 of 31