污水处理厂排放废水对受纳河段沉积微生物群落结构的影响——以太平河为例

田雨露,张柏桓,鄢仕伟,孙昊田,宋进喜,郭家骅

(1.陕西省河流湿地生态与环境重点实验室,陕西 渭南 714099;2.西北大学 城市与环境学院/陕西省地表系统与环境承载力重点实验室,陕西 西安 710127)

随着我国城市化进程的加快与人口的飞速增长,污水处理厂的日处理量与负荷也逐渐增大[1]。这导致污水处理厂出水被直接排入水环境中,从而增加了河流水环境容量负荷[2]。尽管污水处理厂会将污水经过一定处理后达标排放,但是污水处理厂的废水排放依然会对受纳水生态系统的结构和功能产生影响。

沉积物是城市河流系统的重要组成成分, 是流域陆源污染物的“汇”, 也是水体各种污染物的“源”和“汇”[3]。 相较于浮游微生物, 沉积微生物受水体流动的影响较小, 其群落的结构也相对稳定, 能够较好地反映出环境因子对沉积微生物种群的长期影响[4]。 近些年,随着高通量测序技术的不断发展, 人们可以快捷地获取环境样品中的16S rRNA等基因序列, 这些基因序列信息通过与数据库中的已知信息进行比对, 可揭示不同环境介质中微生物群落的特点[5]。 Xie等基于DNA元条形码技术研究原生生物群落演变特征, 可为生态压力源的监测和评估提供有效支撑[6]; Yang等基于生态基因组技术对沉积物中有毒压力源进行生态风险评估, 发现此技术可高效监测环境微生物群落的组成变化[7]; Guo等基于高通量测序技术研究藻类生长, 得到了绿藻群落结构对抗生素的响应机制[8]; Xie等使用eDNA技术研究石油污染下沉积物中微生物群落结构的变化, 发现此技术可有效地评估人为因素对环境的影响[9]; 杨琴等使用16S rRNA 技术揭示了陕北地区石油污染和未污染土壤中细菌群落结构的差异[10]。

目前,关于污水处理厂废水排放对河流沉积物中微生物群落影响的相关研究较少,杨思航等探究了不同气候类型对污水处理厂排放废水的影响,发现污水处理厂排放废水中微生物群落对温度、pH和电导率的响应最为显著[11];赵琛研究了受纳河流微生物群落对污水处理厂排放废水的响应,发现水体中的微生物群落存在空间性差异[12];Wang等对比了不同污染源对太湖沉积物中微生物群落的影响,发现污水处理厂排放污染对其影响最为显著[13]。由此可见,当前研究多集中于揭示微生物群落多样性,而缺少微生物之间相互关系的网络推断分析和菌群的功能潜能预测分析。因此,在特定的时空变化和环境过程的驱动下,微生物之间的相关性尚不清楚,无法准确获取微生物群落在河流生态系统中的生态功能。本研究利用16Sr RNA测序技术,基于应用生物信息学分析平台(QIIME 2),并结合多种统计学分析,研究太平河在西安市第六污水处理厂和第七污水处理厂(以下简称“六污”和“七污”)排水的影响下,河流沉积微生物群落的组成多样性、微生物群落的关联网络和功能潜能预测。在此基础上,揭示显著影响河流沉积微生物群落的环境因子。本研究结合微生物群落、河流水质的理化性质与环境因子,可更为全面地探究污水处理厂排放废水对河流造成的复合效应,同时可为评价西安市其他污水处理厂对受纳河流水生态系统的影响提供科学依据和数据支撑。

1 材料与方法

1.1 研究区概况和采样点布设

西安市位于关中平原中部,北濒渭河,南依秦岭,属于暖温带半湿润大陆性季风气候,春季风多干燥,夏季炎热多雨,秋季天高气爽,冬季寒冷干燥,四季分明但春秋较短,降水分布不均。西安市内地形平坦,河流众多,自古以来素有“八水绕长安”之称。其东西走向的河流有泾河和渭河,南北走向的河流有涝河、沣河、皂河、浐河、灞河等。作为西安市内最大的过境河流,渭河发源于甘肃省渭源县鸟鼠山,于陕西潼关汇入黄河,是黄河最大支流[11,14]。其支流皂河发源于长安区杜曲街道,河流长度为35 km,依次流过长安区、雁塔区、高新区、未央区,于渭河城市运动公园处汇入渭河,是一条主要接纳污水的排污渠。

太平河是皂河的主要支流,全长约26 km,始于西安市第七污水处理厂,流经长安区、高新区、雁塔区、莲湖区、未央区,于未央区八兴滩村附近汇入皂河[15]。太平河是一条人工河渠,主要接纳河流沿岸污水处理厂排水、工业排水、雨水,具有接纳污水,防洪蓄水等功能。本研究以太平河在七污和六污排水口附近的河段为研究对象,根据《水质采样方案设计技术规定(HJ 495—2009)》进行采样点的布设。在七污排水口处(S1)及中下游(S2、S3为中游,S4、S5为下游)共布设5个采样点,在六污排水口上中下游(S6、S7为上游,S8、S9为中游,S10、S11为下游)共布设6个采样点,S1到S11的位置分布如图1所示。

图1 太平河采样点布设图

1.2 河流沉积物和水样的采集与保存

沉积物样品:本次研究于2021年1月13日在全部11个采样点采集沉积物样品,并使用GPS记录了每一个采样点的经纬度坐标。在河道中使用采泥斗和不锈钢铲来采集表层5 cm的沉积物样品,在1 m2的区域随机抽取3个子样本进行混合,作为一个样品,将混合之后的沉积物装入已清洗干净并贴有编号标签的25 mL离心管中。运回实验室之后,于-80 ℃低温保存。24 h内将样品低温(干冰)运送至派森诺生物公司进行高通量测序和微生物群落的分析。

水样:在沉积物采样点处,同时采集距沉积物表面5～10 cm水样,将其装入已清洗干净并贴有编号标签的500 mL采样瓶中。运回实验室之后,于-20 ℃低温保存。1周内完成水样化学分析。

1.3 水样的化学分析

在采样现场,使用川步tp101工业温度计测量水样温度(t,单位:℃);使用YD-1H笔式盐度计测量水样的盐度(salinity,单位:10-6);使用JBP-607A溶解氧分析仪来测量水样的溶解氧含量(DO,单位:mg/L);使用pH-100型pH计来测量水样的pH值;使用ORPscan10笔式ORP计测量水样的氧化还原电位(ORP,单位:mV);使用力辰科技笔式电导率仪/TDS笔测量水样的电导率(conductivity,单位:S/m)和溶解性总固体(TDS,单位:10-6)。

1.4 沉积物样品中微生物测序

本文采用Illuminm测序[16]方法,对沉积物样品进行基因组DNA抽提,抽提完成后,利用0.8%琼脂糖凝胶电泳检测抽提的基因组DNA,利用紫外分光光度计对DNA进行定量分析。PCR采用NEB Q5 DNA高保真聚合酶,选用细菌16S rDNA V3-V4区特异性引物,正向引物为F:ACTCCTACGGGCAGCA,反向引物为R:TCGGACTACHVGGGTWTCTAAT。PCR反应过程见表1。利用Illumina公司的TruSeq Nano DNA LT Library Prep Kit进行文库构建。对合格的文库,在4 MiSeq机器上利用MiSeq Reagent Kit V3(600cycles)进行2×250 bp的双端测序。

表1 PCR 扩增过程

1.5 生物信息学分析

1.5.1 物种分类学注释

对于原始测序数据,首先进行初步筛查,对问题样本进行重测、补测,随后将初步筛查完成后的样本进行文库和样本划分,并去除barcode序列,按照QIIME2 dada2分析流程对序列进行去噪处理,获得ASVs。采用QIIME2的classify-sklearn算法对每个ASV特征序列进行物种分类学水平注释;采用QIIME2对ASVs特征序列进行抽平处理,便于后续分析。

1.5.2 分类学组成分析

采用QIIME2和自编pre脚本统计各样本的分类单元数,并进行分类学组成分析,采用R语言(ggraph、ggplot、ggtree和phyloseq包)对抽平后的ASVs数据进行统计,可以实现各样本在门、纲、目、科、属、种6个分类水平上组成分布的可视化。本研究基于门水平进行分类学组成分析,并以柱状图的形式展现结果。

1.5.3 Alpha多样性分析

采用QIIME2和R语言(ggplot包)计算Alpha多样性指数,使用R语言脚本绘制箱线图,直观地呈现出不同样本间Alpha多样性指数的差异。本研究以7种不同的Alpha指数来表征微生物群落不同的指标。多样性指数标签下的数字为Kruskal-Wallis 检验的p值,当p≤0.05 时,代表其具有较高的统计学差异可信度[17]。

使用R语言的vegan包,绘制稀疏曲线与物种累积曲线。通过绘制稀疏曲线的方式来探究Alpha多样性指数与抽平深度的关系,绘制物种累积曲线以探究样本数量是否能够有效反映群落丰富度。

1.5.4 Beta多样性分析

采用QIIME2和R语言(ape、vegan、ggtree包)计算距离矩阵,进行主座标(PcoA)分析,并将结果绘制成散点图,以表征Beta多样性。

1.5.5 物种差异与标志物种分析

采用Python LEfSe包进行LEfSe (LDA effect size)分析,通过分类学分支图的形式,对所有分类水平同时进行差异分析,以快速找出引起组内差异的物种。使用基于多组学数据的机器学习方法,进行随机森林分析,通过热图与柱状图组合的形式呈现引起组建差异的标志物种。

以上生物信息学分析均在派森诺基因云平台上进行(www.genescloud.cn)。

1.6 统计学分析

1.6.1 关联网络分析

将相关阈值(r)和卡方检验值(p)构建相关性矩阵,使用随机矩阵理论确定阈值,再采用igraph[18]构建关联网络数据,并绘制丰度大小前15的物种模块关联网络图。使用指示关键物种的拓扑指数Zi值(within-module connectivity,节点与其所属模块内其他节点的连接情况)和Pi值(among-module connectivity,节点与其他模块内节点的连接情况)进行分析,将2个值可视化,便于寻找关键物种。

1.6.2 功能潜能预测

使用PICRUSt2软件[19]对微生物群落的代谢通路(pathways)进行比较分析。通过柱状图的形式展现太平河上下游的微生物群落功能差异。再使用分层的样本代谢通路丰度表进行通路的物种组成分析,探究承担差异功能的菌种。

1.6.3 RDA冗余分析

RDA 能够有效地对多个指标进行统计学检验, 并确定对水质环境变化具有最大解释能力的最小变量[20]。本研究使用R语言软件包vegan对 ASVs 的丰度数据(门水平)和水样的水质指标数据进行线性回归整合(见表2),以三序图的形式呈现微生物群落与环境因子之间的关系。

表2 样本测序量统计表

2 结果与讨论

2.1 微生物的测序结果

本次测序一共得到2 403 313对序列,去除低质量序列之后得到2 161 929对序列,去噪处理后得到有效序列2 072 120对,经过拼接和去除嵌合体后得到1 443 218对高质量序列,再将序列中样本丰度为1的单元去除,得到1 399 712对序列(见表2)。本次测序共得到23个门、51个纲、113个目、181个科、276个属和86个种,经过处理得到的ASVs中,物种分类学注释能够注释到的域、门、纲、目、科、属、种的占比分别为1%、0%、4%、3%、27%、59%和6%(见图2)。

图2 物种分类学注释

2.2 微生物群落多样性分析

2.2.1 微生物群落的物种组成

图3为基于门水平的微生物群落的物种组成,由图可知,太平河中优势菌种为变形菌门(proteobacteria)。变形菌种类繁多,在污水处理厂的活性污泥中有较高丰度[21],可以作为河流中的优势菌种。在S1点发现蓝细菌门(cyanobacteria),且丰度很高,这是由于S1点距离七污排污口很近,废水排出导致此处溶解氧和营养物质含量较高(溶解氧8.8 mg/L,总磷1.92 mg/L,均为所有样本中最高含量,见表3),适合蓝细菌的生长[22]。

表3 水样的水质指标数据

图3 基于门水平的分类学组成分析

2.2.2 微生物群落的Alpha多样性

图4 Alpha多样性指数检验

图5为物种稀疏曲线,由图可知,每一个样本的曲线都在抽平深度30 000～40 000 bp趋于平缓。图6为物种累积曲线,由图可知,随着样本数量的增加,曲线逐渐趋于平缓。由此可见,本次测序深度和样本数量都足够反映微生物群落的丰富度变化情况。

图5 稀疏曲线

图6 物种累积曲线

2.2.3 微生物群落的Beta多样性

图7为PCoA分析结果,由图可知,在置信度较高的PCo1维度上(置信度35.3%,大于PCo2维度的14.8%),受七污影响的河段与受六污影响的河段微生物群落之间有较大的差异,以七污排污口和中游与六污中下游的差异最为明显。

图7 PCoA分析结果

2.3 物种差异与标志物种分析

图8为LEfSe分析结果,可以看出,受七污影响引起差异的物种为蓝细菌门、α-变形菌门(alphaproteobacteria)和变形菌门;受六污影响的河段引起差异的物种为酸杆菌门(acidobacteria)、厚壁菌门(firmicutes)和放线菌门(actinobacteria)。

图8 LEfSe分析结果

图9为随机森林分析结果,可以看出,引起样本之间差异最显著的物种为蓝细菌门和厚壁菌门。蓝细菌门在七污口处尤为显著,是因为七污口处溶解氧和营养元素含量较高;而厚壁菌门在六污影响河段中最为显著,六污以改良型AAO工艺为主,陈重军等研究发现厚壁菌常存在于厌氧氨氧化反应器的污泥中[25],厚壁菌随着六污工艺的废水排放进入太平河中。

图9 随机森林分析结果

2.4 关联网络分析

图10为关联网络绘制图,由图可知,蓝细菌门与大部分的变形菌门呈正相关关系,与小部分呈负相关关系;拟杆菌门与蓝细菌门呈较强的负相关关系;衣原体门与绝大部分变形菌门呈负相关关系;小部分变形菌门模块与其他变形菌门模块呈负相关关系。图11将Zi值和Pi值可视化,得到太平河微生物群落中的关键物种为变形菌门中的α-变形菌纲。

图10 关联网络绘制图

图11 纲水平Zi-Pi图

2.5 功能潜能预测

图12为太平河上下游代谢通路的差异,通过对比 KEGG 数据库得到太平河上游显著代谢通路为L-赖氨酸发酵醋酸盐和丁酸盐(P163-PWY);下游显著代谢通路为dTDP-N-乙酰托莫胺生物合成(PWY-7315)。图13为代谢通路物种组成分析的结果,由图可知,承担太平河上游微生物群落功能的主要物种为生氧光细菌(oxyphotobacteria);承担下游微生物群落功能的主要物种为放线菌(actinobacteria)。

图12 代谢通路差异分析

(a) P163-PWY 的物种组成 (b) PWY-7315 的物种组

2.6 环境因子与微生物群落的相关性

图14(a)为采样点微生物群落与环境变量之间的关系,受七污排水影响的河段与受六污排水影响的河段在轴1上具有显著差异,其中:受七污排水影响的河段对总磷、溶解氧和总固体等因子的响应最显著;受六污排水影响的河段对氨氮、氧化还原电位、温度等因子的响应较为显著。从物种角度来看,蓝细菌门与总磷和溶解氧呈正相关,而与盐度和温度呈负相关;厚壁菌门与温度、氧化还原电位、氨氮等呈正相关,与硝酸盐氮、溶解氧等呈负相关。表4为不同环境因子的解释率,其中温度(30.2%)、溶解氧(25.3%)和盐度(9.8%)为解释率最高的环境因子,本次采样时间在冬季,排污口排出的废水温度会随着河流流动而迅速降低,所以温度因子的影响最为显著。图14(b)为它们与微生物群落之间的关系,由图可知,总固体(2.5%)、总氮(1.4%)、CODMn(0.8%)等环境因子对于整个太平河微生物群落的影响不显著。

(a) 采样点微生物群落与环境因子的关系 (b) 微生物群落与起显著影响环境因子的关系

3 结论

本次测序共得到1 443 218对高质量序列,有多半的ASVs可以注释到属以上,证明了本次测序的准确性。在物种组成上,太平河沉积微生物群落优势种为变形菌,随着采样点与排污口距离的增加,微生物群落的丰富度逐渐增加,说明了污水处理厂排放的废水对河流沉积微生物群落存在消极的影响,引起差异的标志物种为蓝细菌门和厚壁菌门。承担群落功能的关键物种为变形菌门中的α-变形菌纲。太平河上游与下游的差异功能是L-赖氨酸发酵醋酸盐和丁酸盐和dTDP-N-乙酰托莫胺生物合成,它们分别是由生氧光细菌和放线菌来承担的。温度、盐度和溶解氧是对菌群影响程度最大的环境因子,其中,蓝细菌门分别和总磷与盐度呈正相关与负相关,厚壁菌门则分别和温度与溶解氧呈正相关与负相关。

综上所述,污水处理厂排放废水会改变河流沉积微生物群落结构,而河流沉积微生物群落结构的变化可以用来反映污水处理厂废水对河流水质的影响,进而可以表征受纳河流生态系统的健康情况。