珠江下游浮游细菌群落结构的时空分布特征*

杜宛璘，孙金辉，麦永湛，赖子尼，贾慧娟,3，葛大艳,4，王 超

(1：中国水产科学研究院珠江水产研究所，广州 510380)(2：天津农学院，天津市水产生态及养殖重点实验室，天津 300384)(3：上海海洋大学水产与生命学院，上海 201306)(4：上海海洋大学生态与环境学院，上海 201306)

河流具有高度的环境异质性和微生物群落多样性，是碳氮等元素转化的重要场所. 浮游细菌作为河流生态系统中的重要组成部分在河流生物地球化学循环中发挥着重要作用[1]，其中包括固氮、脱氮、固碳、硫化物氧化[2]和降解水体污染物[3]等多种功能. 已有研究表明，河流浮游细菌群落组成受时间、空间、温度、氮、磷、pH、溶解氧等[4-5]环境因子的影响，在时空分布上呈现显著差异，而细菌群落组成差异又决定了细菌群落在河流生态系统中发挥不同的功能. 因此，阐明浮游细菌群落时空分布及其环境响应机制，有助于深入了解河流内部环境变化规律，揭示细菌群落组成与河流生态系统功能之间的关系. 目前有很多学者对河流生态系统中的浮游细菌群落进行了研究，例如九龙江[4]、密西西比河[6]、长江[7]、东江[8]等，然而迄今为止，有关我国珠江下游浮游细菌群落的研究仍鲜有报道.

珠江是中国南方最大的水系，是中国境内第三长河流. 珠江是沿江地区人们生活、农业、工业用水的主要来源[9]. 珠江下游流经人口密集、经济发达的珠三角城市群和粤港澳大湾区，常年受工农业污染以及居民生活污水排放的影响，珠江下游生态环境遭到严重破坏. 近年研究表明，珠江流域不同地区都存在污染物浓度逐年升高的趋势[10]，且珠江下游河网区域重金属污染程度较高[11]，这极大限制了珠江水资源的可持续利用. 目前有关珠江流域生态环境的相关研究主要涉及土壤[12]、水质[13]、浮游植物[14]和底栖动物[15]，有关浮游细菌群落的研究相对较少.

高通量测序技术因其快速、准确、大批量、并且免于实验室培育等优点[16-17]，已经成为环境微生物学和微生物生物地理学的主要研究技术手段，并被广泛应用于湿地、土壤和生物膜等复杂环境中微生物群落结构的研究[18-20]. 本文采用16S rDNA扩增子测序方法对珠江下游采集的水体样品进行高通量测序，结合多种生信统计方法，探索珠江下游浮游细菌群落动态变化，以期更全面地了解珠江下游水体状况. 我们假设，珠江下游浮游细菌受环境因子影响将在群落结构上呈现明显差异，这种差异有助于确定珠江下游浮游细菌群落与环境特征之间的关系. 因此，本文旨在了解珠江下游浮游细菌群落多样性及分布情况；确定影响浮游细菌群落多样性及分布的关键驱动因子，为珠江流域水生生态环境保护以及珠江水域污染的生态预警，提供科学资料和理论依据.

1 材料与方法

1.1 调查时间和站位

于2016年6月(丰水期)和2016年12月(枯水期)对珠江下游水域进行野外调查. 在该水域共布设16个采样站位，基本覆盖了西江沿线、珠三角河网中部和广州市周边，具体站位如图1所示. 包括封开(FK)、德庆(DQ)、肇庆(ZQ)、青岐(QQ)、左滩(ZT)、外海(WH)、新围(XW)、小榄(XL)、小塘(XT)、北滘(BJ)、榄核(LH)、横沥(HL)、陈村(CC)、市桥(SQ)、莲花山(LHS)、珠江桥(ZJQ). 其中封开、德庆、肇庆、青岐、左滩、外海、新围位于西江沿线，小榄、小塘、北滘、榄核、横沥、陈村和市桥位于珠三角河网中部，珠江桥和莲花山位于广州市周边.

图1 珠江下游调查部位分布Fig.1 Sampling sites in the downstream of Pearl River

1.2 样品采集与处理

使用5 L有机玻璃采水器采集表层(距离水面0.5 m)水样3 L. 其中2 L水样采集后立即转移到无菌取样袋中(TWIRL EM®，LABPLAS INC.，Quebec，Canada)，放于4℃冰箱下冷藏，尽快运回实验室. 取1.5 L水样等分为3份(生物学重复，避免误差)，分别经200 μm尼龙膜(NY-05，Switzerland)初滤去除大型浮游生物和杂质，再通过0.22 μm孔径的聚碳酸酯滤膜(PALL，USA)过滤，并迅速将滤膜转移至10 mL无菌离心管，保存在-80℃条件下用于DNA提取. 其余样品等分为2份，一份经0.45 μm孔径醋酸纤维素滤膜(Whatman，USA)过滤后保留水样；另一份按照《水质样品的保存和管理技术规定》(HJ 493-2009)标准加入固定剂于4℃保存，用于水化学分析.

1.3 环境因子数据收集

1.4 DNA提取、聚合酶链式反应(PCR)和高通量测序

水体样品总DNA采用DNeasy PowerWater®Kit试剂盒(MO BIO Laboratories，United States)提取，按照生产商提供的说明书进行操作. 用分光光度计(NanoDrop ND-2000，USA)检测DNA纯度和浓度，并通过1% 琼脂糖凝胶电泳进行观察. 利用515正向引物(5′-GTGCCAGCMGCCGCGG-3′)和806反向引物(5′-GGACTACHVGGGTWTCTAAT-3′)[21]扩增细菌16S rRNA基因的V4高变区. 在含有10 μL 5× TransStart®FastPfu Fly缓冲液，5 μL dNTPs(2.5 mmol/L)，2 μL正向引物(5 μmol/L)，2 μL反向引物(5 μmol/L)，1 μL TransStart®FastPfu Fly DNA 聚合酶，2.5 μL 模板DNA(25 ng)和27.5 μL ddH2O的50 μL PCR混合物中构建扩增子文库. PCR反应条件为：95℃ 1 min；95℃ 30 s，55℃ 30 s，72℃ 45 s，共27个循环；72℃ 10 min；4℃终止反应. 对每份样品所制备的一式3份的PCR产物等体积混合，进一步通过1%琼脂糖凝胶检测PCR产物纯化效果，最后在Illumina二代测序平台(Illumina Inc，CA，USA)进行测序.

1.5 生信分析

双向测序完成后下机的测序片段(reads)经FLASH(version 1.2.11)[22]组装，并由QIIME(version 1.8.0)[23]拼接. 本文利用UCHIME(version 4.1)[24]识别并去除嵌合序列(chimeric sequences)，由UPARSE(version 7.0.1001)[25]进行可操作分类单元(OTU，operational taxonomic units)聚类(序列一致性设定为97%)，以0.8置信度阈值对OTU进行注释，并将所有叶绿体、真核生物、未知序列以及OTU单体(singleton)去除.

通过RStudio程序包vegan的diversity函数计算香农-威尔指数(Shannon-Wiener index)和Chao 1指数(Chao 1 index). 采用非度量多维尺度(NMDS，nonmetric multidimensional scaling)和相似性分析(ANOSIM，analysis of similarities)方法分析细菌群落组间差异及时空分布特征. 基于群落Bray-Curtis相异度，采用RStudio程序包vegan和grid进行CCA/RDA分析，评估细菌与环境因子之间的关系.

1.6 功能预测分析

使用PICRUSt软件，基于OTU的物种注释和丰度信息，将KEGG数据库中已有基因组的原核生物16S rRNA序列与SILVA[26]数据库中16S rRNA序列进行关连，然后将KEGG数据库已有的原核物种基因组进行序列打断，利用UProC(version 1.2.0)对所有基因组的KO序列进行统计. 最后利用16S的拷贝数对物种数目进行校正，最终实行KEGG预测以及KO丰度统计，从而获得所测基因的功能预测信息. 原始数据上传到NCBI-SRA数据库，登陆号：PRJNA550236.

2 结果与分析

2.1 珠江下游水体理化特征

2.2 珠江下游浮游细菌群落结构组成

测序结果显示，质量过滤后，珠江下游16个站位所采集的48个表层水样品，一共获得4389036个有效reads，并注释268004个OTUs，平均每个样本的OTU数目为3125±1632.

在门分类水平上，检测到的前15个门类，季节特征显示变形菌(Proteobacteria)在丰水期和枯水期均为最优势门类(图2)，相对百分比分别为50.04%和52.62%. 此外，丰水期优势门类的多样性高于枯水期，除了两个季节均占优势的变形菌、放线菌(Actinobacteria)和拟杆菌(Bacteroidetes)外，厚壁菌(Firmicutes)和蓝细菌(Cyanobacteria)在丰水期的个别站位也会形成较大的优势(图2A). 空间特征显示，变形菌在丰水期的空间分布呈现河网中部的相对丰度高于西江沿线和广州市周边(图2A)，而其在枯水期的空间分布则呈相反趋势(图2B). 此外，厚壁菌在丰水期的西江沿线和河网中部的个别站位，相对丰度会超过变形菌(图2A).

在纲分类水平上，β-变形菌(Betaproteobacteria)在丰水期和枯水期均为最优势菌纲(图3)，相对百分比分别为27.89%和25.62%. 此外，丰水期特有的优势菌纲有芽孢杆菌(Bacilli)和蓝细菌(Cyanobacteria)(图3A)，而枯水期特有的优势菌纲有黄杆菌(Flavobacteriia)和ε-变形菌(Epsilonproteobacteria)(图3B). 空间特征显示，丰水期放线菌呈西江沿线>河网中部>广州市周边的趋势，α-变形菌(Alphaproteobacteria)、 γ-变形菌(Gammaproteobacteria)和芽孢杆菌在西江沿线和河网中部个别站位成为最优势菌纲(图3A)，蓝细菌在广州市周边站位相对丰度明显升高；枯水期β-变形菌在西江沿线的相对丰度波动较大，而在河网中部和广州市周边站位较稳定(图3B).

2.3 珠江下游细菌群落结构多样性

Shannon-Wiener 指数和Chao 1指数分析结果显示，珠江下游浮游细菌群落的季节差异不明显，丰水期略高于枯水期(图4).

表1 珠江下游表层水体环境特征

采用基于归一化OTUs丰度的NMDS分析方法，并结合ANOSIM的R值进行β多样性分析. ANOSIM是显示组间差异和组内差异的一个参数. β多样性分析结果表明，珠江下游浮游细菌群落结构在空间尺度上无显著差异(R=0.033，P=0.331)，而在时间尺度上呈现显著性差异(R=0.222，P=0.001). 所有站位被分为4组，丰水期和枯水期可以明显区分开，组1和组3为丰水期站位，组2和组4除丰水期的青岐和市桥两个站位外，均为枯水期站位(图5). 其次每个季节在空间上可以分为两个组群，丰水期的组1站位位于西江三角洲，组3位于北江三角洲；枯水期的组2大部分站位位于广州市周边，组4大部分位于西江沿线和珠三角河网中部.

图2 珠江下游浮游细菌丰水期(A)和枯水期(B)各门的相对丰度(FK.Jun：封开丰水期；FK.Dec：封开枯水期；其他类同)Fig.2 Relative abundance of each phylum of bacterioplankton in wet season(A)and dry season (B)in the downstream of Pearl River

图3 珠江下游浮游细菌丰水期(A)和枯水期(B)各纲的相对丰度(FK.Jun：封开丰水期；FK.Dec：封开枯水期;其他类同)Fig.3 Relative abundance of each class of bacterioplankton in wet season(A) and dry season (B)in the downstream of Pearl River

图4 珠江下游浮游细菌Shannon-Wiener(A)和Chao 1(B)多样性指数Fig.4 Shannon-Wiener(A) and Chao 1(B)indexes of bacterioplankton along the downstream of Pearl River

图5 珠江下游浮游细菌群落NMDS分析(Jun_ZJQ：丰水期珠江桥，其余类同)Fig.5 NMDS analysis of bacterioplankton community in the downstream of Pearl River

2.4 珠江下游细菌群落结构的环境驱动因子

图6 珠江下游浮游细菌群落与环境因子的RDA分析(Jun_FK：丰水期封开，其余类同；绿色站位代表丰水期站位，粉色站位代表枯水期站位)Fig.6 RDA analysis between bacterioplankton community and environment factors in the downstream of Pearl River

图7 珠江下游浮游细菌前15纲与环境因子的RDA分析Fig.7 RDA analysis between the first 15 classes of bacterioplankton community and environment factors in the downstream of Pearl River

2.5 珠江下游浮游细菌群落功能预测分析

功能预测结果(热图)显示，空间分布上大多数样品中，功能谱的相对丰度是相似的. 转运体(Transporters)、ABC转运体 (ABC transporters)、DNA修复和重组蛋白 (DNA repair and recombination proteins)和双组分系统(two-component system)等是珠江下游浮游细菌群落所表达的主要功能. 季节分布上，丰水期所有站位的相对丰度高于枯水期，其中转运体和ABC转运体的相对丰度丰水期明显高于枯水期 (图8).

图8 丰水期(A)和枯水期(B)细菌功能预测相对丰度热图Fig.8 Heat map of bacterial function prediction relative abundance in wet season(A) and dry season(B)

3 讨论

珠江下游浮游细菌群落结构组成显示，在门分类水平上，变形菌门、放线菌门、拟杆菌门、厚壁菌门、疣微菌门和蓝细菌门是珠江下游的优势菌门. 国内外大型江河生态系统的研究结果(表2)与本研究的结果一致. 在纲分类水平上，β-变形菌、放线菌、α-变形菌和γ-变形菌的相对丰度在两个季节均排名前四位，这与国内外多数淡水生态系统的调查结果也是一致的[27-30].

本研究结果显示，浮游细菌群落相似性组成的季节间差异明显，NMDS结果可以将丰水期和枯水期明显区分开(图5). 纲水平的相对丰度图显示(图3)，芽孢杆菌和黄杆菌在两个季节分布差异较明显. RDA分析的结果显示，WT和SD是影响季节间差异的重要因素(图6). Lv等[31]在黄河三角洲，以及Ma等[32]对海河的研究均表明，温度是导致细菌群落结构季节性变化的一个主要原因. 分析认为，由于每种细菌群落都有自己的最适温度范围，且温度可影响细菌的生长速率和新陈代谢[33-35]，丰水期水温高，适宜细菌生长，而枯水期水温相对较低，会抑制细菌的增殖和代谢[36]，从而抑制细菌生长. 芽孢杆菌最适生长温度为37℃[37]，而此次采样的WT范围为18.0～33.8℃，还未达到芽孢杆菌的最适温度，因此随着WT的升高，芽孢杆菌增加. 而黄杆菌的最适生长温度为25℃[38]，低于芽孢杆菌的最适生长温度，因此两者呈现出的结果不同. 除WT的影响外，Jeon等[39]研究发现，厚壁菌门中的芽孢杆菌纲是沙尘细菌的优势菌群. 丰水期，大量降雨将山体泥沙和岸边泥沙卷入河流地表水中，沙尘细菌的引入从而导致芽孢杆菌纲相对丰度显著提高. 而黄杆菌纲属拟杆菌门，拟杆菌门被认为可以降解高分子量有机化合物[40]，枯水期降雨量减少，雨水对河流的稀释度降低，导致人为输入的有机物含量相对增加，因此促进了黄杆菌纲的增加. 然而，相对丰度较高的细菌种群在两个季节的差异不明显. 分析原因在于：(1)两个季节的水温差异不大，平均值相差小于10℃. 寇文伯等[41]对鄱阳湖浮游细菌的研究发现，7月和10月的水温差异不大，浮游细菌群落的群落结构组成变化也不大. (2)丰水期频繁降雨导致径流量增大，水体混合程度高，透明度降低，这有助于大量附着细菌进入表层水体，造成细菌多样性和丰度增加；然而，径流量增大对浮游细菌群落也会产生稀释作用，导致多样性和丰度减小；枯水期径流量小，稀释作用减弱，但是透明度增大对浮游细菌产生的沉降损失作用增强. 因此从整体上看细菌群落组成变化显著，但相对丰度较高的细菌种群较稳定，所占比例并没有发生显著变化.

功能预测分析表明，转运体、ABC转运体、DNA修复和重组蛋白、双组分系统在细菌功能预测中排名前四. 转运体和ABC转运体均为细胞膜上的一种蛋白，参与跨膜运输功能. 其中ABC转运体在细菌中可参与营养物质的摄取或细菌毒素的分泌[52]. DNA修复蛋白可随时纠正和修复损伤，保持DNA结构、功能完整性[53]；DNA重组又可参与DNA修复，并维持遗传的多样性[54]. 双组分系统信号转导机制，可应对多种环境变化、代谢和细胞周期信号，感应 pH、养分、氧化还原状态、渗透压力和抗生素的双组分系统等[55]. Lv等[31]对黄河三角洲功能预测分析表明，低温细菌群落功能受损. 本研究结果同样显示出，丰水期所有站位的相对丰度高于枯水期，其中转运蛋白和ABC转运蛋白的差异明显. 分析原因(1)枯水期温度低，浮游细菌数量低于丰水期，从而导致功能预测丰度降低. (2)转运蛋白是一种蛋白质，ABC转运蛋白是一种ATP酶，温度降低会抑制蛋白质和酶的活性，微生物群落在低温下活性也较低，降低了微生物代谢[56].

表2 国内外基于16S rRNA测序的河流细菌群落组成及影响因子研究结果*

*-表示文献中未进行细菌群落影响因子的研究.

4 结论

1)珠江下游浮游细菌主要菌门为变形菌门、放线菌门、拟杆菌门、厚壁菌门、疣微菌门和蓝细菌门. 优势菌纲为β-变形菌纲、放线菌纲、α-变形菌纲、γ-变形菌纲、芽孢杆菌纲和黄杆菌纲.

2)珠江下游浮游细菌群落结构存在季节和空间变化，多样性指数丰水期高于枯水期，空间分为3个区域：西江沿线、珠三角河网中部和广州市周边.

4)功能预测分析表明，转运体(Transporters)、ABC转运体(ABC transporters)、DNA修复和重组蛋白(DNA repair and recombination proteins)、双组分系统(two-component system)、嘌呤代谢(Purine metabolism)等是珠江下游浮游细菌群落所表达的主要功能，其中转运体和ABC转运体丰水期明显高于枯水期.