汉江流域上游污水处理厂中微生物多样性及PICRUSt功能预测分析

2022-02-26蒙小俊王秋利龚晓松

安康学院学报 2022年1期

蒙小俊，王秋利，龚晓松

（安康学院旅游与资源环境学院，陕西安康 725000）

随着工农业迅速发展和人们生活水平的提高，城市污水排放量迅速增加，水污染状况日趋严重[1]。WWTPs是控制水污染，保证人类可持续发展的重要基础设施，活性污泥是WWTPs生物处理工艺的功能主体，污泥菌群多样性以及优势菌群之间的相互作用决定了WWTPs生物处理功能的稳定性，群落的改变常常影响WWTPs的处理效率与出水水质[2]3699[3-4]。深入分析污水处理中微生物的群落结构和功能，对控制和提高污水处理的稳定性具有重要作用[5]，可为污水处理系统精准调控和优化设计提供理论依据。

已有研究利用纯培养从污水处理系统中分离出多株功能菌，但WWTPs中微生物数量多达1018个，个体微小、多样性丰富和可培养性低等特征导致传统培养法对微生物开展研究十分困难[6-7]。近年来，随生物技术的发展，高通量测序技术被广泛用于WWTPs中微生物群落结构的研究。Zheng等[8]利用高通量分析技术对外源性进水群落和环境条件如何影响WWTPs膜生物反应器中的活性污泥群落进行了分析，发现进水COD与Pseudomonas呈正相关，进水群落无法为活性污泥提供硝化菌，但能接种反硝化菌，同时增强了活性污泥群落的生物多样性和抵抗外部干扰的能力。李国令等[9]利用高通量测序技术对OAO和AO工艺处理城镇生活污水的微生物群落特征进行了分析，发现两种工艺间微生物群落差异性较大，OAO工艺的微生物种群丰度和多样性均大于AO工艺；AO工艺优势菌属为Pseudomonas、 Thiothrix、 Dechloromonas， OAO 工艺则为 Thermomonas、 Dechloromonas、 Rhodobacter。韩文杰等[10]5037对长三角地区5个采用MBBR泥膜复合工艺污水厂低温季节的微生物群落变化进行了研究，发现各污水厂相对丰度较高的菌属主要有 Nitrospira、Mycobacterium、Defluviicoccus、Hyphomicrobium和Macellibacteroides等，悬浮载体的投加极大程度上强化了优势硝化菌属Nitrospira的富集。张晓红等[11]1896以京津冀区域内典型WWTPs活性污泥为研究对象，利用高通量测序对5个WWTPs活性污泥的微生物种群结构特征进行了详细解析，发现5个WWTPs活性污泥种群结构具有一定差异；不同污水厂，影响其活性污泥群落结构组成的环境因素不同，特殊的进水水质会对污泥菌群组成和生物多样性产生影响。WWTPs活性污泥微生物的多样性和丰度指数受进水水质的影响较为明显，群落组成还受到工艺类型和地理位置等因素的影响，不同工艺、不同水质和不同地理位置的WWTPs微生物多样性存在差异[12]41。

地理位置特殊的汉江流域上游WWTPs以活性污泥为主体的微生物多样性研究鲜见报道。发源于陕南秦巴山地的汉江流域是丹江口水库的主要水源地，其在确保南水北调中线工程水质，促进区域经济发展中具有重要的战略地位[13]，但汉江流域上游支流污染严重、生态环境质量差、潜在威胁多[14]，污水排放是最主要的威胁之一。2017年以来对汉江流域上游污水处理设施进行提标改造，很大程度上缓解了污水对水源涵养地水环境的污染。为保证水源涵养地水质，确保南水北调中线的稳定供水，需要加强该区域WWTPs处理设施的优化和监控。同时，利用高通量分析WWTPs微生物的研究主要关注的重点为群落结构。PICRUSt（Phylogenetic investigation of communities by reconstruction of unobserved states）功能预测分析已在土壤、沉积物、水体和动物胃液等环境中进行应用[15]422[16-18]。孙峰等[15]421对丹江口库区库滨带植被土壤细菌群落多样性及PICRUSt功能预测分析表明，库滨带植物根际细菌主要涉及次生产物代谢的生物合成、转录、多糖生物合成和代谢、细胞生长和死亡等38个子功能。基于高通量测序的PICRUSt功能预测分析在WWTPs微生物研究中分析相对较少，因此对微生物进行功能多样性预测分析，可为提高WWTPs的运行效果提供理论支撑，为污染物去除系统提供优化运行的技术支持。

本研究利用高通量测序技术对采自汉江流域上游5个不同WWTPs好氧单元中活性污泥的微生物群落组成和多样性进行研究，结合PICRUSt分析预测活性污泥微生物功能，从微生物群落和功能角度分析不同样品的差异，以期为汉江流域上游WWTPs的稳定运行、精准调控和优化及丹江口水源涵养地水环境保护提供参考依据。

1 材料与方法

1.1 主要试剂与仪器

E.Z.N.ATMMag-Bind Soil DNAKit， OMEGA；Qubit3.0 DNA检测试剂盒，Life；2×HieffRobust PCR Master Mix，Yeasen；Hieff NGSTMDNA Selection Beads，Yeasen。Pico-21台式离心机，Thermo Fisher；DYCZ-21电泳槽，北京市六一仪器厂；FR-1000凝胶成像系统，上海复日科技有限公司；Q32866 Qubit3.0荧光计，Invitrogen；PCR仪，北京东胜创新生物科技有限公司。

1.2 样品采集

污泥样品分别采自汉江流域上游5个不同的WWTPs好氧单元，每个污泥样品采集3次，采样时间为2020年10月23日—10月29日，所有采集的污泥样品置于聚乙烯密封袋中并排尽空气密封，立即用冰盒储存带回实验室-80℃保存，用于DNA提取和分子生物学实验。样品编号分别为O1、O2、O3、O4和O5，各WWTPs工艺及主要进出水平均水质指标如表1所示。

表1 污水处理厂水质指标mg/L

1.3 微生物总DNA提取、PCR扩增和Illumina测序

污泥样品微生物总DNA提取严格按照OMEGA试剂盒E.Z.N.ATMMag-Bind Soil DNA Kit试验步骤进行。PCR进行两轮扩增。第一轮扩增，利用Qubit3.0 DNA检测试剂盒对基因组DNA精确定量后，确定PCR反应加入的模板DNA量，PCR扩增引物341F：CCTACGGGNGGCWGCAG和805R：GACTACHVGGGTATCTAATCC，扩增区域为V3～V4[19]。反应条件94℃3 min；94℃30 s，45℃20 s，65℃30 s，5个循环；94℃20 s，55℃20 s，72℃30 s，20个循环；72℃5min，10℃。30μLPCR反应体系，2×HieffRobust PCR Master Mix 15 μL，Bar-PCR primer F 1 μL， Primer R 1 μL，模板 PCR 产物10～20ng，ddH2O9～12μ L。第二轮扩增，引入Illumina桥式PCR兼容引物。反应条件95℃3 min；94℃20s，55℃20s，72℃30s，5个循环；72℃5min，10℃。30 μL PCR反应体系，2×HieffRobust PCR Master Mix 15μ L，Primer F 1 μ L，Index-PCR Primer R1μ L，第一轮 PCR 产物 20～30ng，ddH2O9～12μ L。2%琼脂糖凝胶电泳检测文库大小，为得到均匀和高质量的测序数据，使用Qubit3.0荧光定量仪进行文库浓度测定。利用Illumina MiSeq 2000进行测序，委托生工生物工程（上海）股份有限公司完成。将得到的微生物16S rRNA基因高通量测序结果提交至NCBI（登录号为PRJNA698055）。

1.4 序列数据分析

所有测序原始数据利用Prinseq软件经过质控过滤后得到高质量序列，采用Usearch软件进行Operational taxonomic units（OTUs）聚类分析，并利用RDP classifier（OTU物种注释，默认情况下16S扩增片段使用RDP classifier比对RDP数据库（http://rdp.cme.msu.edu/misc/resources.jsp））贝叶斯算法在97%相似水平上对OTU进行分类[11]1899，运行软件Mothur（http://www.mothur.org）/中的summary.Single命令，计算生物相关指数，包括丰富度指数（Chao1/ACE指数）、覆盖率指数（Good’s coverage）和多样性指数（Shannon/Simpson指数）。利用R软件基于主坐标分析（Principal co-ordinates analysis，PCoA）进行β多样性分析。使用Origin、Excel等软件进行微生物的α多样性和β多样性的相关图绘制。

1.5 PICRUSt功能预测分析

PICRUSt是一款菌群功能和代谢的预测软件，可以通过16S rRNA基因序列预测对应细菌的代谢功能谱功能。采用PICRUSt软件进行分析，利用Quantitative Insights Into Microbial Ecology（QⅡME）获得的colsed OTU Table分别与Clusters of Orthologs Groups（COG）数据库和Kyoto Encyclopedia of Genes and Genomes（KEGG）数据库进行比对，获得不同的数据库功能预测信息，具体分析步骤基于在线分析平台（http://picrust.github.io/picrust/）[20]816。

2 结果与分析

2.1 微生物多样性分析

所有污泥样品原始序列经过质控后α多样性指数如表2所示，有效序列在54499～77420之间，样品 OTUs在 1550～2064之间，覆盖率均达到0.99，表明未被检出的微生物种类较少，得到的α多样性指数可靠性较高。微生物组成丰富度稀疏曲线图1表明，随着测序深度增加，OTUs（a）和香农指数均逐渐增加并趋于饱和，说明所用样品的序列数（b）能够充分展示样品中的微生物多样性，与序列高覆盖率（均达到0.99）反映的结果一致。图2显示,尽管每个污泥样品的OTUs较多，但共有的OTUs仅为652，样品 O1、O2、O3、O4和O5独有的OTUs分别为65、83、269、187和47，说明样品O3、O4微生物与其他样品相比差异较大；5个样品香农指数分别为5.84、5.74、6.30、6.13和5.82，辛普森指数分别为0.0069、0.0150、0.0067、0.0066和0.0068，Chao1指数分别为1959、1924、2291、2190和1784，Ace指数分别为1902、1888、2213、2145和1761（见表2）。不难看出，样品O3香农指数、Chao1指数和Ace指数均较高。香农指数越高，表明群落的多样性越高，物种分布越均匀；Chao1指数越大，OTU数目越多，说明该样本物种数比较多；Ace指数越大，表明群落的丰富度越高；Simpson指数值越大，说明群落多样性越低[10]5039，表明O3的多样性相对较高。

表2 微生物α多样性指数

图1 微生物组成丰富度稀疏曲线

图2 微生物组成OUT韦恩图

聚类和PCoA β多样性分析均能描述不同样品间群落的差异，若两个样品距离较近，表示这两个样品的物种组成较相似。图3微生物组成聚类分析（a）和PCoA分析（b）表明，样品O2和O3聚类在一起且距离较近，样品O1和样品O4虽聚类在一起但分布在不同的象限内，样品O5未与其他单个样品聚类且单独分布在一个象限内，说明样品O2和样品O3微生物组成相似度较高，样品O5微生物组成与其他样品相似度较低。

图3 微生物组成聚类分析和PCoA分析

2.2 微生物结构门水平上的组成和丰度

为进一步了解不同污泥样品微生物的组成和丰度，对质控后的序列利用RDP classifier软件进行生物分类，图4（每种门的丰度至少在一个样品中不小于1%）为微生物群落结构在门水平上的组成。优势门主要为 Proteobacteria、Bacteroidetes、Planctomycetes、Chloroflexi和 Acidobacteria，这 5类门类丰度总和在61.72%～78.25%之间（见图4），此研究结果与Ma的研究结果类似[21]。Proteobacteria和Bacteroidetes为WWTPs污泥中常见的优势门[22][23]97，样品O1、O2、O3、O4和O5中Proteobacteria的丰度分别为26.514%、33.864%、36.008%、21.510%和31.790%，Bacteroidetes的丰度分别为15.621%、12.164%、10.071%、6.779%和33.379%。Bacteroidetes和Chloroflexi普遍存在于WWTPs中并参与有机物的降解[24]，Chloroflexi在五个污泥样品中的丰度为3.123%～13.856%。Bacteroidetes和Chloroflexi存在于WWTPs中并对有机物进行降解，有利于脱氮自养菌的生长，从而增强氮素的去除；缺氧条件下反硝化菌以有机物为电子供体将NO3—或者NO2—转化成N2，若碳源不足将影响总氮的去除，工业上常添加甲醇、葡萄糖、乙酸和乙酸钠等作为碳源以提高总氮的去除效率。研究表明，分别添加葡萄糖和乙酸钠为单一碳源时，Acidobacteria和Proteobacteria含量会上升[25]。值得一提的是，Planctomycetes在样品 O1、O2、O3、O4和O5的丰度分别为7.245%、13.866%、13.416%、8.639%和7.327%时，可能存在厌氧氨氧化菌（Anaerobic ammonium oxidation bacteria，AAOB）。王蕾等的研究表明WWTPs好氧污泥中存在一定数量的AAOB（8.63×105～1.90×107拷贝数/g污泥）[26]。以AAOB为驱动的Anammox反应因无需外加碳源、运行成本低、耗氧量低和污泥产率低等优点，成为迄今为止最简捷的污水生物脱氮途径，如何实现主流Anammox并加以工业化应用是未来污水处理领域最热门的研究热点和焦点之一。Actinobacteria在5个污泥样品中的丰度范围为1.295%～3.826%，McLellan等发现Actinobacteria为污水处理厂进水中的优势细菌[27]，推测Actinobacteria可能是外源微生物，群落结构能在一定程度上反映WWTPs的进水性质及所在城市的地理环境特征[28]。Nitrospirae在样品O1、O2、O3、O4和O5中的丰度分别为1.599%、0.237%、1.934%、2.479%和0.310%，表明5个污泥样品具有良好的硝化性能，这与张冰等研究得出的Nitrospirae在城市活性污泥群落中的相对丰度在0.27%～1.61%基本相同[2]3701。Nitrospirae为化能自养菌，是WWTPs中典型的硝化细菌，大量异养菌的存在将与Nitrospirae竞争O2，不利于氮素的硝化，若好氧生物池存在大量COD，大量异养菌会增值，因此应尽可能保持好氧进水较低的COD浓度。除此之外，Verrucomicrobia、Firmicutes和Ignavibacteriae等在5个污泥样品的丰度均较高，未分类的门类丰度在6.455%～21.621%之间。

图4 微生物群落结构在门水平上的组成

2.3 微生物结构属水平上的组成和丰度

微生物属水平上的组成如表3所示（Nitrosomonas、Pseudomonas、Thauera、Zoogloea和Unclassified Saprospiraceae除外，每种属的丰度至少在一个样品中不小于1%）。样品O1、O2、O3、O4和O5中未分类但可以归类到相应属的丰度分别为45.227%、41.094%、43.707%、43.074%和54.260%，推测这些可能是相应属的新种，除Dokdonella、Nitrosomonas、Thauera在样品O2中的丰度和Zoogloea在样品O3中的丰度为零外，所有属在5个样品中均有检出，优势属被所有污泥样品共享，但丰度不同，完全未分类的丰度在6.455%～21.621%之间，这完全不同于贺赟等对北京市5座不同处理工艺污水处理厂污泥生物群落结构的分析结果[29]，而与高晨晨等的研究结果相似[12]40，这可能主要与各系统污水水质、工艺类型和运行参数有关，水质是一个关键的决定因素[30]。

表3 微生物群落结构在属水平上的组成 %

续表

研究表明，Dechloromonas是氨氧化菌，同时具备反硝化除磷的能力[31]，Saprospiraceae可参与脱氮除磷[32]；Anaerolineaceae 是反硝化菌[33]；Sphingobacteriales和 Xanthomonadaceae可进行脱氮[34]；Comamonadaceae是构成活性污泥菌胶团的主要微生物[12]40。常见的亚硝化氧化菌Nitrosomonas在污泥中的丰度在0.042%～0.196%之间，其他参与氮素反应的微生物如 Pseudomonas、Thauera、Terrimonas和Thiobacillus等均有检出，Terrimonas菌能发生好氧反硝化[35]，Pseudomonas能发生好氧反硝化且具有异养硝化的能力[36]，T hiobacillus和

Thauera具有反硝化的功能[37]，Thauera能将硝酸盐还原成亚硝酸盐。构成活性污泥菌胶团的微生物为Flavobacterium、Acinetobacter和 Zoogloea[12]40，丰度在5个污泥样品中的范围分别为0.019%～4.184%，0.086%～1.978%和 0.033%～0.253%，同时Flavobacterium能将周围环境中的有机磷转化为无机磷[23]97，Zoogloea是兼性好氧细菌，能降解有机物，可脱氮除磷[38]，Acinetobacter可进行异养硝化[39]。尽管大部分属为非优势属，丰度在5个样品中的范围为17.440%～28.623%，但污水处理系统功能的稳定性一定程度上需要多样化群落结构间的相互协同，生物处理单元是一个复杂的微生物群体，污水中污染物的有效去除是多种菌群共同作用的结果。上述结果表明，WWTPs中功能微生物多数具有脱氮除磷的功能，脱氮除磷是WWTPs关注的重点，若污水处理不达标或工艺运行失稳，大量氮磷流入汉江将对供水水质造成严重威胁。加强进出水水质监测是关键。脱氮除磷微生物在泥龄（SRT）上存在矛盾，一定条件下，SRT越长，有机负荷越低，硝化细菌越易成为优势菌，硝化效果越好；SRT越短，有机负荷越高，除磷效果越好。为兼顾脱氮除磷，SRT应控制在脱氮菌的最小值和除磷菌的最大值之间。同步脱氮除磷需COD充足，若WWTPs进水COD浓度较低，应缩短进水初沉池的水力停留时间或取消初沉池，若需要，可依据实际情况添加葡萄糖、乙酸和乙酸钠等进行碳源补充，同时可考虑对工艺进行优化分段进水。低温是WWTPs生物脱氮除磷效率降低的重要原因，可适当提高污水水力停留时间和DO含量，增大污泥回流比以及提高污泥SRT，以应对低温的不利影响。生物单元中添加填料形成的MBBR，可提高生物量、处理负荷，强化处理效果，实现原池提量，A2/O、SBR和氧化沟工艺均能较易改造成MBBR，在好氧池原位原池进行，对工艺流程影响小。总体原则为先优化运行，后工程措施；先内部碳源，后外加碳源。

2.4 微生物群落功能预测分析

为获得不同污泥样品微生物的功能，采用PICRUSt软件进行菌群功能预测。利用PICRUSt预测微生物的功能，采用weighted nearest sequenced taxon index（weighted NSTI）指数来评估每个样品中的序列与已知数据库中测定基因的加权平均距离[20]819。NSTI指数的数值越小，说明预测样品中的微生物组成与已知的数据库匹配度越高[15]425，不同样品在COG和KEGG数据库预测的NSTI指数在0.218～0.263之间，表明污泥样品微生物功能预测匹配度较高。COG和KEGG数据库预测的结果分别见图5、图6和图7。

图5 COG数据库对不同样品预测功能基因分析

图6 KEGG数据库对不同样品预测功能基因二级功能分析

图7 KEGG数据库对不同样品预测功能基因三级功能分析

图5显示，能量产生与转化（Energy production and conversion）、氨基酸运输和代谢（Amino acid transport and metabolism）、碳水化合物运输和代谢（Carbohydrate transport and metabolism）、辅酶转运和代谢（Coenzyme transport and metabolism）、脂质运输和代谢（Lipid transport and metabolism）、翻译、核糖体结构和生物合成（Translation，ribosomal structure and biogenesis）、转tion）、复制，重组和修复（Replication，recombination and repair）、细胞壁/细胞膜/膜结构的生物合成（Cell wall/membrane/envelope biogenesis）、翻译后修饰、蛋白转换与分子伴侣（Post-translational modification，protein turnover，and chaperones）、无机离子转运与代谢（Inorganic ion transport and metabolism）以及信号转导机制（Signal transduction mechanisms）十二大类为主要功能基因家族，占比之和范围为66.618%～67.518%。样品O1、O2、O3、O4和O5中的预测功能基因占比相近，氨基酸运输和代谢占比为7.312%～7.786%，产能与转化占比为5.899%～6.290%、无机离子转运与代谢占比为4.786%～5.191%、碳水化合物运输和代谢占比为5.746%～6.748%、信号转导机制占比为6.380%～6.812%、细胞壁和细胞膜占比为6.464%～6.783%，这些研究结果与彭永臻等的研究结果基本一致[23]99。未能很好注释的基本预测功能家族（General function prediction only）和未知功能基因家族（Function unknown）占比分别为12.387%～12.757%和7.504%～8.022%。

图6和图7分别为KEGG数据库对不同样品预测功能基因的二级功能分析和三级功能分析。图6显示，氨基酸代谢（Amino Acid Metabolism）、膜转运（Membrane Transport）、碳水化合物代谢（Carbohydrate Metabolism）、复制和修复（Replication and Repair）、能量代谢（Energy Metabolism）、维他命及辅因子代谢（metabolism of vitamins）、翻译（Translation）和次生产物代谢的生物合成（Biosynthesis of Other Secondary Metabolites）等 20个子功能预测基因占比在5个污泥样品中都约为90%，其中氨基酸代谢、膜转运、碳水化合物代谢、复制和修复和能量代谢为预测的优势功能基因，占比分别为10.343%～10.975%、9.574%～11.337%、10.047%～10.721%、7.001%～7.299%和 5.895%～6.302%，这些检测结果与COG数据库对不同样品预测功能基因分析得出的结论基本相似。图7显示，除丁酸酯代谢、TCA循环和氮代谢的预测基因外，其他所有功能预测基因丰度在每个样品中均高于1%，丰度总和在43.027%～44.938%之间，大部分预测功能基因与氨基酸及其派生物和蛋白质的新陈代谢、遗传信息和物质代谢等有关，这与COG数据库预测得出的结论相吻合。如，检测到主要的碳水化合物的代谢途径糖酵解（Glycolysis/Gluconeogenesis）和三羧酸循环（TCA cycle）的预测功能基因占比分别为1.051%～1.178%和0.856%～0.922%（图7），结果与Yang等的研究结果基本一致[40]，碳水化合物是典型高效的活性污泥中微生物能量产生和生物合成的过程。与氮素代谢的基因在样品O1、O2、O3、O4和O5中的占比分别为0.718%、0.748%、0.741%、0.675%和0.703%；氧化磷酸化（Oxidative phosphorylation）是物质在微生物体内氧化时释放的能量并通过呼吸链供给ADP与无机磷酸合成ATP的偶联反应，基因丰度在1.509%～1.653%之间；氮磷的去除是活性污泥中微生物最重要的代谢活动之一。甲烷代谢的预测基因在1.033%～1.220%之间，表明有部分菌能在厌氧条件下发酵产生甲烷，产生可再次利用的能源。