处理酯化废水EGSB反应器中病毒群落动态变化*

张 堃，谢乔光，徐小云，韩晓鸽，郭雄飞，杜 娟，程学庆

（1.广东轻工职业技术学院， 广东 广州 510300；2.中国电器科学研究院股份有限公司,广东 广州 510300）

地球生态系统中病毒无处不在，几乎可以感染所有的细胞生命体，通过裂解宿主细胞、用辅助代谢基因重编程宿主代谢或对宿主基因进行水平基因转移等来影响宿主功能特性［1-3］。如病毒通过感染裂解宿主细胞，使得胞内营养物质被释放到环境中影响土壤碳循环［4］；通过影响宿主介导的代谢通路溶解环境中的有机碳和总氮［5］；有机肥的长期使用增加了细菌与病毒间的相互作用，提高了两者间功能基因的交换频率［6］。总之，病毒会对宿主微生物的组成和代谢功能造成压力并最终影响元素的地球化学循环。

废水处理设施（Wastewater treatment plants，WWTP）采用独特工程微生物生态系统的活性淤泥进行城市生活污水或工业废水处理，是现代生物技术的重要应用之一［7］。活性淤泥中涵盖多种功能性微生物群体，以细菌群体为主进行污水中营养物质（氮和磷等）或有机物的降解，而细菌群落结构的组成对系统稳定性及淤泥沉降等具有重要的影响［8］。以往研究证明WWTP中含有高浓度的生物量以及营养物质，是病毒良好的栖息地，其中病毒数目可能是自然水生生态系统中含量的 10 到 1 000 倍［9］。Li等人发现 WWTP 的活性淤泥中病毒的存在比较广泛，涵盖一些致病性病毒［10］。Wang等人对2处WWTP的长期研究表明地理隔离的污水处理系统中存在一群核心病毒，对污水处理可能存在未知但稳定的影响［11］。生存条件如温度、pH值和营养状况等的动态变化会直接影响环境中病毒的活性，促进病毒群落的变异［12］。目前为止，对于WWTP连续运行过程中微生物群落组成和功能的研究很多，但对于其中病毒类群的组成和动态变化的研究较少。

由于病毒变异是高度动态的过程［13］，短期研究可能无法全面反映在WWTP的复杂环境中病毒类群的动态变化。本次研究以WWTP中EGSB厌氧反应器中的污泥为研究对象，同时为更好地研究病毒的空间分布差异，我们每月从EGSB罐体的上层、中层及下层分别取样，连续12个月共计获得36个样品，结合宏基因组高通量测序技术对样品内病毒的组成和丰度进行研究，主要包括：（1）鉴定WWTP的EGSB体系内病毒群落组成和多样性；（2）挖掘不同月份、罐体不同部位间病毒的动态变化规律；（3）鉴定系统内病毒-宿主关联性及预测病毒潜在功能组成。本研究加深了对EGSB反应器内病毒群落结构及功能变异模式的理解。

1 材料和方法

1.1 污水处理系统及样品采集

本研究以位于广东省的一处工业规模的生产性膨胀颗粒污泥床（Expanded granular sludge beds，EGSB）反应器污水为研究对象。实验取样从2020年1月持续到12月，每个月15号分别从反应器的上层、中层和底层采样管处用无菌针筒吸取污水污泥样本，并干冰运输立即转移到实验室中-40℃保存。

1.2 DNA提取及宏基因组高通量测序

采用QIAGEN强力土壤基因组DNA提取试剂盒（天根生化科技（北京）有限公司，型号2745288）对所有收集到的36个污泥样品进行总DNA提取，具体提取方法参考试剂盒提供的方法说明书。采用琼脂糖凝胶电泳检验提取到的总DNA的质量，采用Qubit3.0对样品DNA进行准确定量。将符合标准的DNA送至广东美格基因科技有限公司（广州，中国）进行Illumina NovaSeq 6000双端测序。平均每个样品原始测序数据量约 12 Gb，总数据量 440 Gb。

1.3 病毒序列分析

利用 Trimmomatic（Version 0.39）软件对下机数据进行质控；采用Megahit（Version 1.2.5-beta）软件进行数据组装；CD-HIT（Version 4.6.7）进行序列去冗余（相似性＞90%，覆盖度＞80%）。VIBRANT（Version 1.2.1）对筛选保留到的长度5 Kb以上的序列进行病毒鉴定；采用CheckV（Version 0.8.1）对预测到的病毒信息进行质量评估，去除质量评级中未明确的序列。采用Prodigal（Version 2.6.3）对鉴定获得的病毒序列进行ORF预测，并结合DRAMv（Version 1.2.4）进行功能注释。采用VPF-tools对病毒序列进行物种归类及宿主预测，其中对病毒物种注释域水平上隶属度（Membership ratio）与置信度（Confidence score）阈值设定0.2以上；科水平二者阈值设定为0.3；属水平二者阈值设定为0.25；而对病毒-不同分类层级宿主（域、门和科）预测中隶属度（Membership ratio）与置信度（Confidence score）阈值均设定 0.2 以上［14］。采用Bowtie 2计算病毒在各样品中的相对丰度。采用 QIIME2（Version 2020.8）进行病毒群落α和β多样性分析。

1.4 统计分析

样品间共存和特有病毒类型花瓣图采用R语言包“Venn diagram”进行分析绘制；采用R（Version 4.2.1）对不同处理组间病毒群落显著性差异统计分析，P＜0.05判定为组间存在显著性差异；从NCBI refseq下载已知的噬菌体terL氨基酸序列，与本次注释到的terL基于最大似然法（Maximum likelihood）采用 MEGA（Version X）构建蛋白树，并采用iTOL在线网页进行树图美化；采用本地Python编制的脚本（Local script）进行病毒与宿主关联桑基图绘制等。

2 结果与讨论

2.1 病毒序列鉴定及丰度动态变化

36个样品中共获得3 360条非冗余病毒序列，仅7.5%的序列为溶原性病毒（Prophage），一般认为溶原性病毒对宿主的代谢以及环境元素循环的影响远低于烈性噬菌体［15］。样品中Chao1指数在 1 503.60 到 3 046.80 间；Shannon 指数在8.20到9.58间，说明持续运行过程中污泥内病毒群落存在较高的动态变化。其次，对组间α多样性的显著性差异分析发现Chao1指数在连续月份间呈现规律性变动趋势：多样性先从一二月份显著升高至三四月份后急剧降低，并在十月份开始再次发生显著升高至十二月；而Shannon指数则无明显规律（见表1）。由于病毒丰度不均匀，导致病毒丰富度最高的三四月份整体多样性并不最高，可能存在一些优势病毒类型影响样品中病毒整体多样性的体现。反应器不同部位的病毒群落均不存在显著性差异（P＜0.05），可能差异明显的采样月份是导致反应器不同部位不存在显著差异的主要原因。

表1 组间Chao1和Shannon指数差异显著性分析

采用Bray-Curtis PCoA法进行不同分组间病毒类群整体离散性分析（见图1）。样品按照采样月份能各自聚类，表明相应月份内的病毒整体组成上较为相似；而在采样部位主导下，分组内离散性较高，不存在明显聚类模式。有研究表明AS体系中微生物群落结构出现明显季节性演替，春季和冬季采集的样品细菌群落更相似［16］。有趣的是，本体系中病毒群落在春季和冬季也出现同样的规律，可能源于病毒主要依靠侵染的宿主存活，其更倾向于直接受到病毒群落结构的影响［17］。

图1 基于Bray-Curtis PCoA的不同分组下样品散点聚类分析

2.2 病毒物种组成及遗传多样性

所有样品中长尾噬菌体科（Siphoviridae）（35.91%±7.84%）、短尾病毒科（Podoviridae）（13.35%±4.54%） 和肌 病毒 科（Myoviridae）（10.44%±2.54%）相对丰度最高，均属于噬菌 体 目（Caudovirales）（ 见 表2）。 与 以 往 对WWTP中病毒类群的研究结果一致［18］。样品间病 毒 属 以 T4virus（4.74%±2.41%）、P22virus（2.94% ±1.45%）、Lambdavirus（2.63%±0.68%）及P12024virus（2.21%±0.40%）等为主，其中丰T4virus属的成员被认为具有裂解性但不会产生毒性蛋白［19］。

表2 病毒科层级物种组成信息统计表

噬菌体目（Caudovirales）在不同科中末端酶大亚基（terminase large subunit，terL）具有一定保守性，利用terL基因的系统发育分析，可以了解噬菌体的进化地位和遗传距离［20］。本次提取到227个terL氨基酸序列中存在部分与参考蛋白差异过大，并在系统遗传树中分支为一些独立的簇（见图2），表明本研究体系中未知噬菌体科存在较高多样性和新颖性，与近期Shi等人的研究一致［21］。

图2 基于terL氨基酸序列的噬菌体科水平系统发育树。红色的树枝以及红色的圆点代表未能准确定位的噬菌体蛋白。

2.3 病毒-宿主关联预测

体系中鉴定到 2 785 条噬菌体（Phage或Prophage）和3条噬真菌体（Mycophage）。由于病毒的复杂侵袭能力，还存在532条病毒序列可以同时侵染细菌、真菌或古菌。Gulino等人对纽约14处污水处理厂的研究发现活性淤泥中噬菌体的比例高达90%［22］， 与本研究结果一致。

环境条件影响噬菌体对宿主的选择范围［23］，本体系中噬菌体主要侵袭的宿主细菌包括变形杆菌门（Proteobacteria）、放线杆菌门（Actinobacteria）及厚壁杆菌门（Firmicutes）等；主要宿主属包括分支杆菌属（Mycobacterium）、假单胞杆菌属（Pseudomonas）及芽孢杆菌属（Bacillus） 等（ 见 图3）。 其 中， 短 尾 病毒科（Podoviridae）预测到的宿主大部分属于变形杆菌门（Proteobacteria）；长尾噬菌体科（Siphoviridae）的宿主则以放线杆菌门（Actinobacteria）的分支杆菌属（Mycobacterium）为主。噬菌体-宿主感染的模式会对潜在的菌群进化和元素循环产生影响，而噬菌体-宿主的感染可以跨越多种门和属的范围［24］。

图3 病毒科与宿主属或门间的潜在关联桑基图

2.4 病毒促氮代谢功能潜力

废水可能是环境中活性氮的最大来源之一，会导致水体和生态系统富营养化；而WWTP中的工程微生物可促进氮代谢，经过包括反硝化、硝化反应及氨同化作用等途径促进活性氮的转变［25］。本次共检测到7条噬菌体含有氮相关功能基因（见图4），其中ctg1051569为沙粒病毒属（Omegavirus），被检测含有亚硝酸盐还原酶基因（nitrite reductase，nirK），可以在反硝化通路中促进亚硝酸盐还原为氧化一氮［26］。而ctg467926（Cd119virus）含有羟胺还原酶基因（hydroxylamine reductase，hcp），可以将羟胺还原为铵根离子进入厌氧氨氧化途径或氨同化作用途径［27］。其余5条病毒序列被预测到含有氨同化作用过程中的天冬酰胺合成酶（asparagine synthase，ASN），研究发现该酶可能会提高宿主耐铵性［28］。

图4 7条含有氮代谢相关基因的病毒结构基因组序列图

此外，ctg1051569（Omegavirus）在10月份和12月份时丰度较高，提示可能该月份下反硝化作用细菌丰度可能较高；ctg467926（Cd119virus）在8月份样品中最高，提示该月份下厌氧氨氧化途径或氨同化作用途径参与的宿主可能丰度较高；拥有天冬酰胺合成酶（asnB）的病毒在月份间分布较为分散，提示各个月份下促进氨同化作用的病毒丰度差异较小；整体上十月和十一月拥有氮代谢相关AMG功能的病毒丰度最高，说明这两个月份下可能菌群代谢更为活跃（见图5）。

图5 7条关键病毒在12个月份间的丰度聚类热图

3 结论

（1）连续运行12个月的WWTP活性淤泥内病毒群落在月份间存在较高的动态变化，Chao1指数出现有规律的动态变化，即呈先升高后降低并再次升高的趋势。

（2）各自月份下的病毒群落相似性较高，但采样罐部位的影响不显著；春天和冬天体系下的样品中病毒群落组成更为接近。

（3）淤泥体系中烈性噬菌体占比最高，且存在一些高度新颖且多样化的病毒科。

（4）体系中含有氮代谢相关AMG的噬菌体，且在十及十一月份相对丰度更高。