基于高通量测序的宁波三江干流微生物多样性分析

张晓镇，陆志波，季树勋，夏珊珊，王 颖

（1. 同济大学环境科学与工程学院，上海 200092；2. 宁波市水利发展规划研究中心，浙江 宁波 315000）

河流是淡水生态系统的重要组成部分，河道生态系统中的浮游微生物是河道生物多样性的重要组成部分，且与河流中溶解氧、pH、硝态氮、总磷和重金属等众多环境因子息息相关[1-3]，做好河流微生物多样性研究是河流健康评价、水系连通与生态修复、供水安全保障等工作的重要基础。

近 30 年来，各国形成了侧重点不同、各具特色的评价方法，如，美国环境保护局（US EPA）提出的快速生物监测协议（Rapid Bio-assessment Protocols，RBP）[4]、欧盟的《水框架法令》(WFD)[5]等都已在河流水生态系统健康管理中得到应用。2020 年世界自然基金会联合我国8 家科研院所共同发布《长江生命力报告2020》，提出长江生命力指数。这些评价方法均将水生态作为重要维度。

微生物作为有机化合物的主要分解和矿化者[6-7]，在食物链和食物网中所处的营养级比大型水生动植物更低[8]，水体微生物网络的复杂性一定程度上能反应水体中大型生物的网络复杂性与社区对扰动的抵抗力和恢复力[9]。微生物群落结构、丰度等也可以作为有机物分解、反硝化功能的指标[10]。

近一个世纪以来，地表水微生物安全主要是通过检测粪便污染的细菌指标来确定[11]，但除致病菌外，在地表水微生物监测中污染指示微生物更能反应水体受污染的现状[12]。有研究者[13]提出可以创建可靠、标准化的参考数据库，作为指示健康微生物群落的基线指标。通过相关性计算可以分析环境因子与物种分布间的相关关系[14]。已有研究发现微生物多样性指数与总磷和磷酸盐含量呈极显著负相关关系[15]。SOUFFREAU et al[16]也发现总磷和氨氮对湖泊的微生物群落组成有显著影响。

微生物多样性研究的传统方法是微生物分离纯化培养，存在如培养周期长、分类学鉴定不准确等弊端[17]。分子生物学技术可以高效、准确地从遗传物质水平上研究环境中微生物的群落结构[18]。高通量测序技术是近年来微生物群落分析的常用技术[19-21]，具备测序快速、通量高等优势，且可对多样品的多个可变区同时测序[22]。通过高通量测序技术研究三江干流水体微生物群落结构及其多样性，可以为构建三江流域生态系统管理方案提供科学依据。

宁波市的余姚江、奉化江和甬江合称“三江”。三江口由于其独特的自然地理条件，成为宁波市的重要标志[23-24]。三江河道的治理方案主要为水利工程建设[20]，需要对河道水生态进行长期综合监测以提供数据支撑，极有必要对其中微生物群落结构与多样性展开全面研究。当前对于宁波三江流域微生物群落结构与多样性的研究较为缺乏，本文开展的研究说明通过高通量测序的方法开展河流的生态指标监测具有可行性。

1 材料与方法

1.1 样品采集

此次三江干流巡航测量范围为余姚江（牟山闸/永思桥～蜀山大闸～姚江大闸～三江口）、奉化江（三江口～澄浪堰）、甬江（三江口～入海口）。航测总长度约105 km。以下称三江干流均为上述范围。巡航测量时间为2021 年5 月11～14 日。测量方式为使用航测船只沿余姚江、奉化江和甬江深泓航道线单向巡航一次。途中在三江干流12 处点位停船进行水体取样测量。测量断面及取样点位置分布，见表1 和图1。

表1 三江干流航测断面及取样点位置分布

图1 宁波三江干流采样点分布

取样过滤装置含真空泵或蠕动泵、六联过滤器和直径为47 mm、孔径0.45 μm 的混合纤维素酯膜。

1.2 生物信息提取与分析

DNA 提取和测序在Personal Biotechnology Co. Ltd.（中国上海）进行。按照制造商的说明，通过DNeasy PowerSoil Kit（QIAGEN，Hilden，Germany）提取水样中总微生物基因组 DNA。通过琼脂糖凝胶电泳和 NanoDrop ND-1 000 分光光度计 (Thermo Fisher Scientific, Waltham, MA, USA) 评估提取的DNA 的质量和数量。

使用引物对515F_926R 和18SV4F_18SV4R来分别对 16S rRNA 基因、18S rRNA 基因的可变区进行PCR 扩增。

采用Illumina HiSeq X-10（Illumina，San Diego，CA，USA）平台对群落DNA 片段进行双端（Pairedend）测序。使用DADA2 方法进行去引物，质量过滤，去噪(denoise)，拼接和去嵌合体等步骤。用QIIME2 (2 019.4)进行物种分类学注释。将有效数据在97%的相似水平上划分为不同的运算分类单位(operational taxonomic unit, OTU), 将OTU 代表序列与相应的微生物数据库比对, 利用Greengene数据等数据库对物种进行注释, 从而得到每个样本所含的物种信息。通过对OTUs 进行Alpha 多样性指数、丰度等分析, 并对物种注释在各个分类水平上进行群落结构的统计分析, 得到细菌群落的组成信息, 再通过Beta 多样性分析, 系统分析三江干流不同流域样本间的微生物群落结构差异。

使用Alpha 多样性来衡量局域均匀生境下的物种数目。α 多样性使用Hill’s number family 多样性计算[25]，主要包括Chao 指数；Faith’s pd 指数；Simpson 指数；Shannon 指数；Good’s coverage 指数；Pielou’s evenness 指数；Observed species 指数。

Chao1 是基于丰度的估计量，计算见式（1）[26]：

式中，Schao1为预测丰富度，Sobs为观测到的物种数量，n1为单序列OTUs 数目，n2为双序列OTUs 数目。

Faith’s pd 指数反映了谱系多样性，等于生命发育树（The Tree of Life Web）的分支数量[27]。

Simpson 指数是衡量个体在分类时的集中程度，见式（2）[25]：

式中，Sobs为观测到的OTUs 数量，ni为第i个OUT 对应的个体数目，N为群落中所含的总个体数目（所有 OTU 的丰度之和）。

Shannon 指数量化了预测从数据集中随机抽取的个体物种身份的不确定性，见式（3）[28]：

式中，Sobs为观测到的OTUs 数量，ni为第i个OUT对应的个体数目，N为群落中所含的总个体数目。

Good’s coverage 指数用于评估测序结果代表的community 物种的总数[29]，见式（4）：

式中，F是单序列 OTU 的数量，N是群落中所含的总个体数目。

Pielou’s evenness 指数用于描述物种中的个体的相对丰富度或所占比例[30]，见式（5）：

式中，H’为Shannon 指数的值，H’max为H’的最大值。

Observed species 指数为群落中丰度＞0 的物种数之和，可以反映物种种类信息。

Beta 多样性是指沿着环境梯度变化的不同群落之间，物种组成的相异性或物种沿环境梯度的更替速率，因此也被称为生境间多样性(betweenhabitat diversity)。β多样性分析所使用的Braycurtis 距离基于OTUs 的计数统计，比较2 个群落微生物的组成差异，见式（6）[31]：

式（6）中，SA,i和SB,i表示第i个OTU 分别在A 群落和B 群落中的计数。min 表示取两者最小值。

2 结果与讨论

2.1 三江干流微生物群落组成

通过Silva 数据库对所得数据进行比对，根据序列物种分类学注释的结果对物种组成进行分类，展示相对丰度最高的10 个分类，其余合并显示为others。统计样本的物种注释结果中门(Phylum) 、纲(Class) 、目(Order) 、科(Family) 和属(Genus)水平上含有的分类单元的数量。

对于原核微生物，在门水平上，变形菌门(Proteobacteria)的平均丰度最高（图2a），在C1、C3、C4 和PC2 点位处相对丰度均超过60%，但在靠近入海口的C13 点位处拟杆菌门(Firmicutes)的相对丰度最高，达到42.8%。需要注意在C4、C5 和C7 等采样水域厚壁菌门(Firmicutes) 表现出较高丰度，说明已经存在污水排放导致的DO 含量下降，存在厌氧环境[32]。在纲水平上（图2b），γ-变形菌纲(Gammaproteobacteria) 的平均丰度最高，其中在C1-C5、PC1、PC2 的7 个点位处相对丰度均超过50%，在临近入海口的C13 处γ-变形菌纲的丰度明显下降，拟杆菌纲(Bacteroidia)、梭菌纲(Clostridia)的相对丰度呈现上升趋势。在门水平与纲水平上，优势菌种与其他研究中的结果一致[33-34]。在目水平上（图2c ），在余姚江河段，假单胞菌目(Pseudomonadales)平均相对丰度较高，达到46.4%，在奉化江和甬江河段异常球菌目(Deinococcales)相对丰度可达到20%以上，在C9 处达到最高值42.7%。在科水平上（图2d），莫拉菌科(Moraxellaceae)在余姚江河段各采样点处都表现出较高的丰度，在C2-C5 点位处丰度均高于40%，尤其在C4 处达到54.5%，异常球菌科(Deinococcaceae)在甬江流域也表现出较高丰度，在C9 处达到最高值42.7%。在属水平上，不动杆菌属(Acinetobacter)的平均丰度最高，异常球菌属(Deinococcus)在甬江流域也表现出了较高的丰度，最高在C9 支流交汇处，达到42.7%。最丰富的细菌种类包括奇异球菌(Deinococcus grandis)(9.44%)、鲁氏不动杆菌(Acinetobacter lwoffii) (1.18%)、波西米亚不动杆菌(Acinetobacter bohemicus) (1.09%)等。

图2 原核生物群落结构组成

对于真核微生物，在门水平上，绿藻门(Chlorophyta)的平均丰度最高（图3a），在C2、C6 和PC2 点位处相对丰度均超过40%，在C6、PC1、PC2 和C12 等采样点处子囊菌门(Ascomycota) 也表现出较高的相对丰度，在PC1 点位处子囊菌门(Ascomycota)的相对丰度最高，达到65.5%。在纲水平上（图3b），绿藻纲 (Chlorophyceae)的平均丰度最高，为25.6%，其次为引藻纲(Cryptomonadales)(12.7%)、甲藻纲(Dinophyceae)(7.8%)等分类。在临近入海口的C12、C13 处绿藻纲的丰度明显下降，拟杆菌纲(Bacteroidia)、梭菌纲(Clostridia)的相对丰度呈现上升趋势。在目水平上（图3c），在C1～C6， 环 藻 目 (Sphaeropleales)、 海 链 藻 目(Thalassiosirales)相对丰度较高，达到46.4%，在C7、C10 处，隐藻目(Cryptomonadales)相对丰度分别达到85.0%、69.1%，表现出显著优势。在科水平上（图3d），隐鞭藻科(Cryptomonadaceae) 平均丰度达到13.2%，在C7、C10 处表现出显著的高丰度，尤其在C7 处达到85.0%，但在C1～C5 各点位处其丰度并不显著。在属水平上，隐藻属(Cryptomonas)的平均丰度最高，主要在C7、C10 处表现出显著的高丰度，小环藻属(Cyclotella)在C1 处相对丰度达到46.6%。最丰富的真核微生物种类包括弯隐藻(Cryptomonas curvata) (6.85%)、寄生索囊藻(Choricystis parasitica)(0.79%)和马索隐藻(Cryptomonas marssonii) (0.62%)等。

图3 真核生物群落结构组成

2.2 微生物群落Alpha 多样性分析

以Chao1 和Observed species 指数表征丰富度，以Shannon 和Simpson 指数表征多样性，以Faith’s PD 指数表征基于进化的多样性，以Pielou’s evenness 指数表征均匀度，以Good’s coverage 指数表征覆盖度，微生物群落Alpha 多样性分析结果，见图4。

图4a 可知，对于原核微生物，YYJ 样品的Chao1 和Observed species 指数最高， FHJ 和YJ 样品的Chao1 和Observed species 指数值差异不大，说明余姚江流域的群落丰富度最高，Chao1 值在3 000 以上。YYJ 样品的Shannon 指数值相对较高，但3 组样品整体差异不显著，Simpson 指数基本稳定在0.98 以上，说明三江干流流域各采样点原核微生物群落多样性较丰富，稳定性较好，其中余姚江流域的原核微生物多样性相对更高。YYJ 流域样品的Faith’s PD 指数较小，而YJ 流域样品Faith’s PD指数较大，FHJ 居中，说明三江干流自上游至下游原核微生物的亲缘关系逐渐复杂化，基于进化的多样性较高。YYJ 流域样品的Pielou’s evenness 指数较高，FHJ 其次，YJ 最小，说明3 个流域的物种均匀度自上游向下游逐渐减低。Good’s coverage 指数值均在0.97 以上，说明测序对群落中物种的覆盖度较高，未被检测出的物种所占的比例很少。与澜沧江[35]、科希河[36]和博斯腾湖流域[37]等流域相比，宁波三江流域原核微生物群落丰富度与多样性更高。

图4b 可知，对于真核微生物，FHJ 样品的Chao1和Observed species 指数最高， YJ 其次，YYJ 样品最低，说明奉化江采样监测流域的群落丰富度最高。FHJ 样品的Shannon 指数值相对较高，67%的采样点的Simpson 指数在0.9 以上，3 组样品Simpson 指数值无显著差异，说明三江干流流域各采样点真核微生物群落多样性较丰富，稳定性较好，其中FHJ 流域的真核微生物多样性相对更高。YYJ 流域样品的Faith’s PD 指数较小，而YJ 流域样品Faith’s PD 指数较大，FHJ 居中，说明三江干流自上游至下游真核微生物的亲缘关系与原核微生物的发展一致，也呈现逐渐复杂化的趋势，基于进化的多样性较高。FHJ 流域样品的Pielou’s evenness 指数较高，YJ 其次，YYJ 最小，说明3 个流域的物种均匀度自上游向下游逐渐升高。Good’s coverage 指数值均在0.999 以上，说明测序对群落中物种的覆盖度较高，未被检测出的物种所占的比例极少。与九龙江[38]、西安渭河、灞河和浐河[39]等流域相比，宁波三江流域真核微生物丰富度更低，多样性相近。

图4 Alpha 多样性分析

2.3 微生物群落Beta 多样性分析

根据Bray-curtis 距离算法，对物种进行非度量多维尺度(non-metric multidi-mensional scaling，NMDS)分析(图5)，发现YYJ 的样品在维度NMDS1 和维度NMDS2 能够很好地与FHJ、YJ 区分开，说明余姚江流域样品的原核、真核微生物种类与FHJ、YJ 有较为明显的区分性。其中Stress 值分别为0.057 1、0.107，均＜0.2，很好地证明了数据的可靠性。

图5 基于相似性＞95%的OTU 水平上非度量多维度分析

采用层次聚类(hierarchical clustering)的分析方法，以等级树的形式展示样本间的相似度，通过聚类树的分枝长度衡量聚类效果的好坏。对Braycurtis 距离矩阵采用UPGMA 算法（即聚类方法为average）在门水平上进行聚类分析，见图6。可以看出，三江干流的YYJ、FHJ 和YJ 的3 部分采集的样品中的原核微生物在不同河段呈现出鲜明的聚类效果，说明三江干流自上游向下游原核微生物在门水平上有比较明显的变化，不同流域的原核微生物相似度较高，入海口C13 处则展示出完全不同的群落组成。位置相近的采样点的真核微生物的聚类也比较相似，如干流上游方位的C1～C5 点位处样本之间的相似度很高，但与相对处于下游方向的点位之间的相似度较低。在支流汇入后的样品的群落构成会出现比较明显的变化，如C9 采样点与其他点位的相似度就相对较低。

图6 门水平上的层次聚类分析

2.4 物种组成热图

根据测序结果和注释OTU 数据进行样本间丰度和相似性聚类，构建宁波三江干流微生物属层次物种分类Heatmap 图谱，以颜色梯度变化表现丰度值的高低，见图7。使用平均丰度前30 位的属的丰度数据绘制热图。样本按照物种组成数据的欧式距离(euclidean distance)进行UPGMA聚类。

图7 三江干流微生物属层次物种分类Heatmap

原核微生物在属上的Heatmap 图谱显示，YYJ 流域样品中微小杆菌属(Exiguobacterium) 的丰度最高，次为马赛菌属(Massilia)、黄杆菌属(Flavobacterium)、芽孢杆菌属(Bacillus)和不动杆菌属(Acinetobacter)等；FHJ 流域样品中丰度最高的属为新鞘脂菌属(Novosphingobium)、砂单胞菌属(Arenimonas)和短波单胞菌属 (Brevundimonas)等；YJ 流域样品中分支杆菌属(Mycobacterium)、鞘脂菌属(Sphingobium)和异常球菌属(Deinococcus)等丰度更高。原核乳杆菌属(Lactobacillus)、鼠杆菌属(Muribaculaceae)等只在入海口C13 处测得，该处组成与其余点位有明显的区别。

真核微生物在属水平上的Heatmap 图谱显示，YYJ 流域样品中绿藻属(Pseudomuriella)、麦可属(Mychonastes) 和海链藻属(Thalassiosira)的丰度最高， FHJ 样品中囊藻属(Tetracystis)、 扁藻属(Tetraselmis)和直链藻属(Aulacoseira)表现出较高的丰度，YJ 样品中尖尾藻属(Teleaulax)、隐藻属(Cryptomonas)相对丰度较高。入海口C13 采样点处测得铜绿微囊藻属(Hererocapsa)、小二仙草属(Haloragis)、剑水蚤属(Mesocyclops)和赤潮异弯藻属(Heterosigma)等仅在此点位显示出较高丰度的真核生物属。

3 结论与展望

3.1 结论

通过高通量测序，发现三江干流优势菌门主要有变形菌门、放线菌门、异常球菌-栖热菌门和拟杆菌门，真核微生物优势门包括子囊菌门、甲藻门、绿藻门和硅藻门等，这些均为典型淡水微生物门类，与其他水域微生物群落组成检测结构相似[40-42]。如刘峰等[43]发现汾河入黄口水体中变形菌门、放线菌门和厚壁菌门处于优势地位，在湖泊水体中变形菌门、蓝细菌门、拟杆菌门、放线菌门和疣微菌门也是最主要的5 个原核生物门类[44]，这些典型的淡水菌门在淡水生态系统中的广泛分布已经得到研究证实[45-46]。在三江干流下游异常球菌-栖热菌门丰度的增加可能意味着该河段受到较多的热污染影响[45]。在西安渭河、沣河、灞河和浐河的研究中，绿藻和真菌界子囊菌门也是在各类群中的优势物种[39]，与宁波三江流域真核微生物优势门一致，个别点位处（如C4、C9 等处）甲藻门丰度较高可能与温度的升高或pH 的降低有关[46]，与原核微生物的优势菌门中异常球菌-栖热菌门的出现共同说明宁波三江干流受到热污染的可能性很高。

同时通过Beta 多样性分析也注意到，不同流域微生物群落组成表现出明显的区分性，如异常球菌-栖热菌门、子囊菌门在三江干流的中下游水域丰度更高，而在上游的采样点C1～C5 处的样品中丰度较低。厌氧性厚壁菌的出现说明存在工业或生活废水的排入，已经造成水生态环境的破坏[47]。入海口处水域的微生物组成也与其他采样点处明显不同。由物种组成热图发现，微生物组成在相邻采样点间相似性较高，YYJ、FHJ 和YJ 3 组样品中也有不同的优势物种，可以作为流域水环境状况的指示生物。

Alpha 多样性分析反映了宁波三江干流流域微生物多样性的综合指标，自三江干流上游向下游进行采样检测，各采样点水域处的原核微生物丰富度指数和多样性指数均较高，YYJ 样品与FHJ 样品、YJ 样品相比较，原核微生物丰富度差异不大，YYJ流域样品多样性最高。真核微生物方面，三江干流流域整体原核微生物丰富度与其他流域研究结果相比偏低，多样性维持在较高的水平，其中FHJ 样品的物种丰富度和多样性相对更高。自上游至下游微生物的亲缘关系呈现逐渐复杂化的趋势。

3.2 展望

本实验结果仅就宁波三江干流流域微生物多样性进行表征，水体微生物多样性受到多种环境因子的影响，一方面水体微生物多样性与环境因子之间的关系，不同水域所处的地理位置对于水体微生物群落的影响都值得进一步探索。

同时，水体的微生物作为重要的环境指示生物也可以作为水生态环境现状的指示，通过研究流域微生物群落结构多样性与群落结构，可以为构建完善的监测与管理体系提供理论支持。在常态化生态监测工作中微生物多样性的监测与分析工作也具备可行性。通过对微生物群落结构多样性与其驱动因子的研究，可以开发宁波三江干流流域特殊功能菌群，为区域水环境治理提供针对性支持。