基于宏基因组学的三峡库区重庆主城段水体浮游微生物群落的组成和功能分析

罗 芳，鲁伦慧，李 哲，付 川，闫 彬

（1.重庆三峡学院 环境与化学工程学院，重庆 404100；2.中国科学院 重庆绿色智能技术研究院，重庆 400714；3.中国科学院水库水环境重点实验室，重庆 400714）

微生物在水生生态系统中种类繁多且密度大，影响着淡水生态系统功能的发挥和人类健康[1]。因此，淡水生态系统中的微生物成为研究的热点，主要研究内容包括浮游微生物群落组成结构及其与环境因子之间的相互作用关系等。三峡大坝的建设改变了大坝上下游河流的水文和水生生态[2]，水环境的物理、化学和生物特性均发生了一定程度的变化[3]。水环境改变会影响水库系统中的水生微生物，尤其是对环境因子变化敏感的浮游细菌[4]。浮游生物对环境刺激表现出复杂和敏感的反应，可以作为生态系统发生变化的指标[5]。

三峡库区是一种特殊水库生态系统[6]，人为扰动及三峡工程运行对库区环境均有直接影响。三峡库区重庆主城段是经济和人口增长迅速的区域之一[7]，研究该区域的浮游微生物群落组成结构及功能，可以进一步理解水体浮游微对修建三峡大坝环境变化的响应，及其水库系统中发挥的生物地球化学循环作用[3]。

宏基因组这一概念在1998年第一次提出，它是指生境中微生物群落总的基因组信息[8]。传统的微生物独立培养只能获得1%的微生物[9]，宏基因组学避开了微生物分离培养，直接从微生物群落获得基因组信息[10]。宏基因组学的诞生推动了对浮游细菌结构[11]、海洋环境中原核生物[12-13]和真核生物[14-15]群落组成、河流和湖泊生态系统浮游生物群落[16-17]的研究。目前，学者们利用16S扩增子测序技术[4]和宏基因组测序技术[18]揭示了浮游细菌在三峡水库生态系统的多样性。然而关于三峡水库浮游微生物群落基因组研究，包括基因和功能的研究明显不足。

本研究主要采用HiSeq 4000（Illumina，美国）高通量测序技术对三峡库区重庆主城三处水域环境样本浮游微生物（主要针对浮游细菌和浮游植物）进行分析，以获取微生物群落物种分类信息（组成结构）、功能基因信息（基因组成）以及基因功能（COG function，KEGG pathway），为理解其地球化学循环机制奠定基础。同时通过对重庆主城段三处断面水域水环境不同理化因子的比较分析，为水体水环境保护提供基础数据。

1 材料与方法

1.1 研究区域与样品采集

2015年5月12日于春季水华前期在三峡库区重庆主城段共设置3个断面，分别位于重庆市上游的朱沱、下游的木洞，以及嘉陵江支流的北碚。朱沱镇地处长江上游北岸，位于重庆市永川南端。朱沱断面（105°51′26″E，29°01′27″N）属于亚热带季风气候，年平均温度为17.7 ℃。木洞断面（106°32′01″E，29°25′18″N） 位于巴南区北部，北临长江。北碚断面（106°21′09″E，29°54′25″N）属亚热带季风性湿润气候，多年平均温度为18.2 ℃。各断面按重庆主城上游（长江干流）、重庆主城下游（长江干流）及嘉陵江支流划分。在各断面沿河道中泓线设置3个采样点（图1），每个采样点取10 L水，采集表层（0.5 m深度）水体并分析其水质情况。现场将每个断面3个采样点的5 L表层水样充分混合，取1 L 水通过1.2 μm 滤膜（Whatman，美国）后再由0.22 μm 混合纤维素酯滤膜（Millipore，美国）过滤，滤膜放置液氮罐内保存运输，运回实验室后放置于-80 ℃超低温冰箱，进行样品DNA 的提取。剩余水样放置于冷藏保温箱运回实验室后，立即开展理化指标的测试。

图1 三峡库区重庆主城段采样断面分布图Fig.1 The distribution of sampling sites in the Chongqing urban section of Three Gorges Reservoir

1.2 水体基本理化指标测定

现场测试分析的指标包括：pH 值、电导率（Cond）、溶解氧（dissolved oxygen，DO）、水温（T）、气温、大气压和总碱度（total alkalinity，TA）。其中pH 值、电导率、DO 值、水温由YSI ProODO 测定，现场气温、大气压由手持式数字大气压测定，TA 由GDYS-103SI 总碱度测定仪测定。水体的有机碳（total organic carbon，TOC）、总碳（total carbon，TC）、水溶性有机碳（water-soluble orgnic carbon，WSOC） 由 vario TOC 分析仪 （elementar，GER）测定。硝态氮（N-NO3）、亚硝态氮（N-NO2）、氨氮（N-NH4）、总磷（total phosphorus，TP）、溶解性正磷酸盐（soluble orthophosphates，SRP）、总氮（total nitrogen，TN）、叶绿素a（Chl a）（估算水体浮游植物生物量水平与营养状态的重要指标），测定方法参考文献[19]。

1.3 宏基因组测序

1.3.1 全基因组测序分析

总DNA 的提取使用D5525-01Water DNA Kit（Omega，美国），提取方法参见试剂盒使用说明书。将提取的DNA 样本送到上海美吉生物医药科技有限公司进行质检。DNA完整性检测使用1%琼脂糖凝胶电泳；DNA 纯度检测使用Nanodrop2000（赛默飞世尔科技，美国）；DNA 的浓度检测使用TBS-380。送检的DNA 样本经检测满足：有明显主带，且清晰完整，无降解；OD260/280在1.7～2.2之间（为1.89～1.90）。将检测合格的DNA样品构建短输入（500 bp）文库，使用HiSeq 4000平台（Illumina，美国）进行测序。有效数据（Clean Data）的获取方法：使用Seqprep（https://github.com/jstjohn/SeqPrep），对序列3′ 端和5′ 端进行质量剪切，获得不含N碱基的较长序列片断（reads）；使用Sickle（https://github.com/najoshi/sickle）保留的高质量的双末端和单端reads，即舍去剪切后的部分reads（长度＜50 bp以及平均质量值＜20）。

1.3.2 基因拼接和预测

使用软件SOAPdenovo （http://soap.genomics.org.cn/, Version 1.06），对高质量的reads 进行拼接和组装，然后根据长度为k的核苷酸序列（kmer）之间的重叠关系，构建德布鲁因图（De-Brujin graphs），获取重叠群（contigs），把reads 映射（mapping）到contig，根据reads 之间的双末端（pair-end）关系把contigs 拼接成多个contig 的连接物（scaffolds）。将kmer 值设定为 39～47，在scaffolds的gap处，将scaffolds打断成新的contigs，优选contig（长度≥500 bp）。contig 的 ORF 预测使用 MetaGene（http://metagene.cb.k.u-tokyo.ac.jp/）。并将其基因序列（核酸长度≥100 bp）翻译为氨基酸序列。

1.3.3 非冗余基因集的构建与基因注释

CD-HIT软件（http://www.bioinformatics.org/cdhit/）将所有样品预测出来的基因序列按照95%一致性、90%覆盖率进行聚类，选取各类中长度最长的序列进行非冗余基因集的构建。使用SOAPaligner 软件 （http://soap.genomics.org.cn/），按照95%一致性将每个样品中经优化的reads 与非冗余基因集对比，然后统计样品中的基因丰度信息。

选用 BLASTP（BLAST Version 2.2.28+，http://blast.ncbi.nlm.nih.gov/Blast.cgi）软件将得到的样品中的基因集序列分别与非冗余蛋白质的氨基酸序列（Amino acid sequence of non-redundant protein，NR）、同源聚类基因群（homologous cluster gene group，eggNOG）和京都基因和基因组（Kyoto Encyclopedia of Genes and Genomes， KEGG） 数据库进行比对，即可分别获得物种、COG 功能和KEGG 功能注释。比对参数e-value 均设置为 1×10-5。

2 结果

2.1 水环境基本参数

各水域断面的水环境理化参数值见表1。采样当天，水温（T）在各水域断面变化很小。Cond，TP，TN，TC，Chl a在各采样点间变化较大。与朱沱、木洞采样点的理化参数值相比，北碚水域水环境中Cond，N-NO2，TN，Chl a含量最高，Chl a含量为2.09 μg/L，pH值，TP，TC，TOC以及DOC含量最小，pH 小于7，呈酸性。朱沱水体SRP、N-NO3和TC均是三者中最高的。木洞点位水体TP含量为0.092 mg/L，但是Chl a 含量却是最小的。朱沱和木洞水体的大部分理化因子数据比较接近，可能是由于都位于长江干流段，有机质含量较高，水体环境比较类似。

参考《地表水环境质量标准》（GB 3838—2002）对三峡库区重庆段（北碚、朱沱、木洞三个断面）的TP、TN和DO进行评价。在3个断面水体中，TP含量为0.023～0.092 mg/L，只有北碚水体符合Ⅱ类标准，朱沱、木洞水体TP 含量均小于0.1 mg/L，符合Ⅳ类水质标准；TN 含量均大于1.5 mg/L，符合Ⅴ类水质标准。TP 和TN 的浓度远远超过国际公认的水体富营养化TN、TP含量的阈值 （TP 和 TN 浓度分别为 0.02 mg/L 和 0.2 mg/L）时，从营养盐单因子考虑，超过藻类疯长的“水华”现象的水体富营养化初始浓度限值[20]，是春季水华的主要成因。DO 含量为7.84～8.01 mg/L，符合Ⅰ类水标准。可以初步判断，水体氮污染比较严重，TN，TP 是影响重庆段水体水质的主要因素。

表1 各断面的理化因子数据Tab.1 Summary of the physiochemical characteristics for each of the sample sites

2.2 三峡库区重庆主城段水体浮游微生物基因集

宏基因组学测序共产生1.23 亿条原始读长（每个样品0.39 亿～0.45 亿），3 组样品（北碚、朱沱和木洞）的原始序列长度分别为5.82 Gb，6.07 Gb和6.75 Gb。单个样品经数据质控后的Clean reads的序列长度占原始序列长度的93%左右。浮游微生物宏基因组共获得230 740条组装序列（contigs），单个样品contigs 总序列长度分别为86.03 Mb，80.16 Mb 和 62.99 Mb （N50 为 944～1 046 bp）。这些组装序列共产生398 263个（每个样品11.23万～14.69 万个）开放阅读框序列（ORFs），ORFs 在3 组样品中的长度分别为76.76 Mb，72.76 Mb 和57.24 Mb（每个样品平均长度509.81～522.96 bp），每个样品的最长ORF 序列长度为13.59 kb 左右，最短ORF 序列为100 bp。去冗余前所有样品的基因数为39.83万条（每个样品11.23万～14.69万条），所有基因的总序列长度为0.21 Gb。非冗余基因集基因数为28.46万条，非冗余基因集基因的总序列长度为0.15 Gb，平均序列长度为536.47 bp。

2.3 基于宏基因组测序的物种与功能注释

2.3.1 物种分类学注释

根据宏基因组序列分析，共得到5域6界78门139 纲 314 目 571 科 1 727 属 6 771 种 。 在 域 水 平上（表2），3组样品水体微生物类别中相对丰度最大的均是细菌（Bacteria）和病毒（Virus），真核生物（Eukaryota）和古生菌（Archaea）含量较低。所有样品门分类学水平上的物种丰度条形图见图2。将各样品物种在不同分类中相对丰度小于0.1%的门类合并到“Others”组，一共注释到8 个门类。可以看出放线菌门、变形菌门和拟杆菌门在北碚、朱沱、木洞水体相对丰度均较高，依次为44.39%，24.22% 和 9.32%； 36.28%， 26.59% 和 17.85%；44.45%，28.09%和11.43%。蓝细菌门在北碚水体含量略高于拟杆菌门，为9.71%，而在朱沱和木洞水体含量仅为0.51%和0.49%。

表2 所测样品在域水平上物种的相对丰度Tab.2 The relative percent contribution of planktonic mi crobe community at the domain level for the different sample

图2 所有水域环境样本浮游微生物门分类水平物种分类条形图Fig.2 The relative percent contribution of planktonic microbe community at the phyla level for the different sample

将基因集与NR数据库进行比对，选取属水平上丰度前50 的属绘制热谱图（图3）。norank（属未定）在北碚、朱沱和木洞断面含量分别高达49.47%，41.72%和47.62%。除norank（属未定）外，含量较为丰富的为蓝细菌门中的颤藻属和微鞘藻属，未分类放线菌门的菌属（Actinobacte-ria_unclassified）、沉积岩杆菌属 （Ilumatobacter）、ß-变形菌门的Limnohabitans菌属，土壤杆菌属（Sediminibacterium）、乳酸球菌属 （Lactococcus）、噬冷菌属（Algoriphagus）和黄杆菌属（Flavobacterium）。值得注意的是，微鞘藻属和颤藻属只出现在北碚断面的嘉陵江表层水体，其余两处几乎没有。

图3 属水平上不同断面浮游微生物群落的热谱图（相对含量前50个属类别）Fig.3 Heatmap based on relative abundances of the top 50 genera of planktonic microbial communities among the different cross sections

2.3.2 COG功能注释

将基因产物与eggNOG 数据库①http://eggnog.embl.de/（evolutionary genealogy of genes: Non-supervised Orthologous Groups）对比获得对应的直系同源序列聚类（COG），得到所有样本COG功能分类统计图（图4）。在3组样品宏基因组序列中注释到未知功能（Function unknown，S）的基因序列最多，北碚断面水域环境样本浮游微生物基因序列数目已经接近8×105。其次为注释到氨基酸转运和代谢功能（Amino acid transport and metabolism，E）的基因序列，在每个样品中所占的比例分别为11.34%，11.10% 和11.63%。然后为注释到翻译、核糖体结构和生物合成功能（Translation, ribosomal structure and biogenesis，J）的基因序列，分别占总体的8.68%，8.27%和8.72%。注释到能量生成和转化功能（Energy production and conversion，C）的基因序列数排在第四位，基因序列的相对丰度在北碚、朱沱和木洞断面水域分别为8.42%，8.23%和8.72%，与翻译、核糖体结构和生物合成功能相对丰度比较接近。

2.3.3 KEGG功能注释

KEGG数据库是一个整合基因组、化学和系统功能信息的数据库。本研究以KEGG途径数据库为依据，在对3组样品一共99 472条基因序列进行代谢通路分析时，有54 340条基因注释到KEGG数据库中。统计被注释的基因发现，只得到了3 505个K编号。再对注释到的216个通路分析可得3组样品注释序列较多的代谢通路为ABC 转运（ko02010）、嘧啶代谢（ko00240）和嘌呤代谢（ko0230），其中嘌呤代谢通路相对丰度最大，分别为北碚4.95%、朱沱5.10%和木洞5.01%。将注释得到的代谢通路中与碳、氮、硫循环相关的部分代谢路径绘制成一个热图，见图5。3 个断面水体中的浮游微生物参与碳元素循环相关的甲烷代谢（ko00680）、原核生物的碳固定（ko00720）、三羧酸循环（ko00020）的代谢路径较多。微生物参与氮素循环相关的氮代谢（ko00910）途径在朱沱和木洞两处断面水体样本的相对丰度分别为1.13%和1.14%，略低于北碚断面（1.23%）。与硫循环相关的代谢（ko00920）在北碚、朱沱、木洞3个断面水体样本中相对丰度偏小，分别为0.41%、0.45%和0.42%。

图5 基于碳、氮、硫代谢的不同断面水体KEGG 路径热谱图Fig.5 Heatmap showing the differences among three sampling sections based on the KEGG pathway in carbon, nitrogen and sulphur metabolism

2.4 三峡库区重庆段水体浮游微生物组成分析

2.4.1 多样本物种、功能和基因数目比较

用Venn 图来统计3 组样品的浮游微生物群落的物种、功能和基因数目的组成情况（图6）。从图中可以看出3组样品共有基因数高达170 606条，朱沱和木洞水体样本共有的基因数目最多，达到211 902 条。北碚、朱沱和木洞水体样本独有的基因数分别为18 724 条、5 583 条和4 165 条，相比之下北碚水体浮游微生物独有的基因数最多。对于属水平上的物种数目，三组样本共有的属达到1 599 个，而木洞独有2 个属，根据统计得到的组间成分表，获得这两个属是紫球藻属（Porphyridium）和噬菌体属（Tunalikevirus）。3组样本中共有的COG 功能数目和KEGG功能数目分别为8 125个和207个，独有的KEGG 功能数目较少。

图6 多样本比较Venn图Fig.6 Venn diagram of multi-sample comparison

值得关注的是北碚水体的蓝细菌门（Cyanobacteria），属于原核生物，其相对丰度为9.71%，这一值明显高于朱沱水体的0.51%和木洞水体的0.49%，但其绝对含量，仍需相关的定量手段来加以佐证。另外“Others”这一类别的物种相对丰度占比较小，说明在门分类水平上所选取物种可以很好地反映整体样品的丰度。此外，在淡水生态系统扮演初级生产者角色的浮游植物，即浮游藻类，在物质、能量等方面发挥着极大的作用[21]。对于真核藻类，在研究区域的表层水体中，仅发现了少量的硅藻门（Bacillariophyta），分别占北碚、朱沱和木洞水体总物种的0.02%，0.01%和0.01%。

2.4.2 样本的非度量多维标度（NMDS）分析

使用布雷-柯蒂斯（Bray-Curtis）距离计算样本两两间的距离，分析样本间的物种或功能丰度的分布差异。将样本间的物种或功能距离矩阵用热图（图7）表示，分别可得样本间物种组成矩阵[图7 （a）]、基因矩阵[图7 （b）]、COG 功能的矩阵[图7 （c）]和 KEGG 功能矩阵热图[图7 （d）]。结果显示：朱沱和木洞水体在物种、基因、COG功能上的差异是最小的，北碚和朱沱水体之间差异最大。在KEGG功能上，北碚和朱沱水体的差异最大，朱沱和木洞水体具有相似性。NMDS排序图在朱沱和木洞的物种组成相似性与浮游微生物在门和属分类水平上的聚类图具有一致性。

图7 样本的非度量多维标度排序图Fig.7 Non-metric multidimensional scaling (NMDS)ordination diagram of samples

2.5 三峡库区重庆段理化因子与浮游微生物类群间的关系

选取门分类水平上丰度靠前的12 个物种和环境水文因子进行Pearson 相关性分析（表3）。由表3 可知，放线菌门与大部分环境因子呈现高度相关性。其中与N-NO2和N-NH4显著正相关（r=1，P＜0.05），与N-NO3呈显著负相关（r=-1，P＜0.05），与DOC 呈极显著负相关（P＜0.01）。绿弯菌门（Chloroflexi）与pH 值呈显著正相关（P＜0.05）。蓝细菌门与T呈极显著正相关（P＜0.01），而与TC呈极显著负相关（P＜0.01）。异常球菌-栖热菌门（Deinococcus-Thermus）与DO呈显著负相关（P＜0.05）。芽单胞菌门（Gemmatimonadetes）与Chl a呈显著负相关（P＜0.05）。变形菌门与TP呈显著负相关（P＜0.05）。

表3 浮游微生物优势门与环境因子的Pearson相关性分析相关系数表Tab.3 Pearson correlationship analysis on dominant phyla of planktonic microbial communities with environmental factors

3 讨论

3.1 三峡库区重庆段水体浮游微生物优势菌属群落结构分析

通过研究属水平上的物种丰度（图3）发现，在朱沱断面，拟杆菌门的黄杆菌属和β-变形菌门的Limnohabitans菌属相对丰度最高。已有研究发现黄杆菌属在沣河和鄱阳湖微生物群落中大量存在，黄杆菌属主要靠分解有机物生存，且具有溶藻活性[21]。此外，黄杆菌属的某些菌种具有致病作用，主要表现为对鱼的致病作用，需要重视对该菌属的防护[22]。已有研究表明Limnohabitans菌属在水体浮游细菌群落中扮演着重要的角色，因为它们在藻类衍生的基质上有很高的基质吸收和生长速率[23]。在北碚断面水体发现大量的颤藻属和微鞘藻属，它们都属于蓝细菌门。高营养浓度导致北碚断面水体Chl a 偏高，这与在该断面水体中蓝细菌的生物量较高具有一致性。三峡库区重庆段水体的浮游微生物物种丰度分布遵循生态学最古老、最普遍的规律之一，即水体中稀有物种占多数，常见物种占少数[24]。

3.2 三峡库区重庆主城段水环境浮游微生物的功能分析

在COG注释（图4）中，绝大多数的基因序列归类于代谢类下的氨基酸转运，代谢亚类和能量产生与转化亚类，信息存储及处理类的翻译、核糖体结构和生物合成亚类。根据KEGG 路径热图（图5）分析，在北碚断面水体样本中的浮游微生物参与能量代谢相关的光合作用途径多于另外两处断面水体，这与得到的北碚断面表层水体Chl a含量明显高于另外两处断面水体的结果具有一致性。统计可得，北碚断面水体基因注释到与参与氮循环相关代谢途径的数目高于另外两处断面水体，且所得到北碚断面水体硝酸盐还原酶高于另外两处断面水体，这与表1中北碚断面水体亚硝酸盐氮含量高于另外两处断面水体的结果相符。此外，根据KEGG数据库得到注释序列最多的代谢通路与膜运输和核苷酸代谢路径有关。通过NMDS排序图（图7）可知朱沱和木洞水体中浮游微生物具有的COG功能和KEGG功能类别是相似的。KEGG路径的热图中包含的芳香族化合物多环芳烃的降解，可能涉及微生物降解过程。有研究表明，黄杆菌属能以低分子量多环芳烃作为唯一碳源和能源，降解和矿化多环芳烃[25]。KEGG路径热图中包含甲烷代谢和硫代谢路径，本研究还未能分辨参与的甲烷代谢过程为产甲烷过程、甲烷氧化过程还是同时包括这两个过程，如若探讨甲烷微生物过程机制，需对相应的功能基因做进一步分析。KEGG 路径热图中涉及了氮代谢过程，参与硝化、反硝化的菌种的确定和功能分析是接下来需要进一步分析的内容。KEGG热图中的原核生物碳固定路径以及光合作用路径相对丰度较高，可能是蓝藻参与反应。

3.3 环境因子对三峡库区重庆主城段水环境浮游微生物群落的影响

为了确定水体理化因子对浮游微生物（主要是浮游细菌）群落组成的影响，本研究进行了Pearson 相关性分析 （表3）。T，TC，DOC 是影响三峡库区重庆段浮游微生物群落的主要环境因子。因为微生物都有一个利于生存的最适温度，温度往往是潜在的限制性因素[26-27]。在相关性分析结果中，蓝细菌门与温度具有明显的相关性，这与之前的研究[28]相符。在确定对蓝细菌门生长的影响因素时应结合理化及水文学等参数综合考虑。相关性分析发现芽单胞菌门与Chl a 呈显著负相关（P＜0.05），Chl a与浮游植物的生物量关系密切[29]。这可能是由于浮游植物的某物种抑制了芽单胞菌门的生长。

4 结论

（1）基于宏基因组学技术对重庆段进行分析，得到的物种注释表明：在门水平上，重庆段水体（北碚、朱沱和木洞）四大优势菌门为放线菌门、变形菌门、蓝细菌门和拟杆菌门。古菌的优势类群为广古菌门。在属水平上，除属未定的微生物以外，北碚水体的颤藻属和微鞘藻含量最高；朱沱水体含量最高的是变形菌门的Limnohabitans菌属；木洞水体含量最高的为Limnohabitans菌属，其次为沉积岩杆菌属。COG功能注释表明：3个采样点水体的微生物解释得最多的功能就是氨基酸转运和代谢。

（2）Pearson 相关分析结果表明：浮游微生物群落丰度与环境因子有着很强的关联性，其中T，TC和DOC的影响程度最大。

（3）在本研究中，水体的营养状态均是通过Chl a，TN和TP、浮游植物群落结构进行单因子评价。从Chl a 单因子评价，研究水体均属于贫营养状态。从TN和TP的指标来看，各调查水域均处于中营养至重富营养状态。如果从浮游植物群落进行评价，比如硅藻门，可以初步判断调查水体水质属于中营养型，因为中营养型湖泊通常是甲藻、隐藻、硅藻占优势[30]。但是考虑到水库生态系统本身的复杂性和野外环境的不确定因素等，还需要通过结合三峡水库重庆段水体的评价标准来评价水体的贫富营养状态。