点源炼油废水处理系统微生物多样性及代谢功能研究

宋佳宇， 秦 榕， 杜显元， 吴慧君， 梁立宝

1.石油石化污染物控制与处理国家重点实验室， 北京 102206

2.中国石油集团安全环保技术研究院有限公司， 北京 102206

3.中国石油大学(北京)化学工程与环境学院， 北京 102206

4.中国石油新疆油田公司百口泉采油厂， 新疆 克拉玛依 834011

炼油废水是一种高浓度有机废水，具有水量大、污染物种类多、生物降解难度大等特点[1-3]. 目前，石化企业普遍采用隔油—浮选—生化“老三套”工艺处理此类难降解点源废水[4]，但保障其出水水质稳定达标GB 31570—2015《石油炼制工业污染物排放标准》尚有困难. 此外，原油重质化、劣质化情况加剧，导致点源炼油废水处理系统下游生物处理单元冲击风险加大，污水处理厂无法实现长周期稳定运行[5].

作为生物处理单元的保护屏障，上游一级处理单元对保障废水处理系统稳定高效运行具有重要意义[6-7]. 目前，已报道的研究多集中于对炼油废水处理系统生物处理单元活性污泥群落结构变化的评价[8-12]，如石油类物质对活性污泥沉降性能影响及其干扰处理效果评价[13]，以及硫化物对活性污泥中亚硝酸盐氧化细菌(NOB)和氨氧化细菌(AOB)的抑制作用[14]等. 关于炼油废水处理系统沿程微生物群落结构变化及代谢功能的研究却鲜见报道.

目前，国内外研究者主要采用荧光原位杂交技术(FISH)、变性梯度凝胶电泳技术(DGGE)、末端限制性片段长度多态性技术(T-RFLP)、克隆文库(Clone library)和高通量测序等技术[15-18]开展微生物群落结构研究，其中，高通量测序技术因具有测序通量高、准确、全面等优势[19]，已成为微生物生态结构和功能解析的强大技术支持手段.

综上，该文以某石化污水处理厂点源炼油废水为研究对象，开展废水微生物多样性及群落结构研究，明确点源炼油废水处理系统(简称“系统”)微生物群落动态行为，阐明环境因子与废水微生物群落的响应关系，预测其代谢功能，以期为现行废水处理工艺综合诊断及优化奠定理论基础.

1 材料与方法

1.1 废水处理系统处理效果及采样信息

某石化企业污水处理厂运行稳定，COD、TOC (总有机碳)、TN、NH3-N、石油类等水质排放指标达到GB 31570—2015《石油炼制工业污染物排放标准》排放标准. 该污水处理厂一级处理单元各出水口设置采样点，采样点布设如图1所示，水样依采样顺序用英文字母编号(分别为A、B、C、D、E). 每个样品采集3个平行样，共计45个样品. 样品采集后，0～4 ℃冷藏运输至实验室，48 h内完成水质指标检测分析. 各单元主要污染物去除贡献率如表1所示. 由表1可知，该系统沿程主要污染物的去除效果显著，去除率均在90%以上.

表1 各单元主要污染物去除贡献情况

图1 点源炼油废水处理系统采样点分布

1.2 水质指标及测试方法

水质指标分为无机指标、有机指标和金属指标3类(见表2)，其中，无机指标包括pH、DO、ORP (氧化还原电位)、电导率、TN、NH3-N、TP ；有机指标包括COD、TOC、石油类；金属指标包括Cd、Cu、Ni、Ba、Cr、Mn、Fe、Mg、Se、As、Pb、Al、Ca、K、Zn、Na. 分析测试方法参照《水和废水监测分析方法》(第四版)[20].

表2 水质指标及测试方法

1.3 高通量测序方法

1.3.1样本前处理

采用无菌微孔膜(0.22 μm，美国密理博公司)及无油式真空泵(Rocker 300，台湾洛科仪器有限公司)抽滤1 L水样，所得滤膜即为微生物菌群样品，立即放入-80 ℃冰箱保存.

1.3.2DNA抽提、PCR扩增及Pair End文库构建

根据DNeasy®96 PowerSoil®Pro QIAcube®HT Kit 说明书进行微生物群落总DNA抽提，使用超微量分光光度计(NanoDrop2000，美国赛默飞世尔科技公司)测定DNA浓度和纯度(1%的琼脂糖凝胶). 确定DNA浓度及纯度后，微生物多样性分析通过引物338F (5′-ACTCCTACGGGAGGCAGCAG-3′)和806R(5′-GG ACTACHVGGGTWTCTAAT-3′)对16S rRNA基因V3～V4可变区进行PCR扩增；宏基因组测序通过非接触式超声波破碎仪(Covaris M 220，美国Covaris公司)将DNA片段化，筛选约400 bp的片段，用于构建PE文库.

1.3.3上机测序

微生物多样性采用Miseq PE300(Illumina，美国Illumina公司)测序平台进行测序；宏基因组采用NovaSeq 6000 (Illumina，美国Illumina公司)测序平台测序. 上述测序平台由上海美吉生物医药科技有限公司提供.

2 结果与讨论

2.1 系统沿程微生物多样性变化情况

Alpha多样性指数[21]是用于评价微生物群落多样性和丰富度的综合指标，其中，Shannon-Wiener和Simpson指数反映群落多样性[22]，Ace和Chao指数反映群落丰富度. Shannon-Wiener指数越大，说明群落多样性越高；与之相反，Simpson指数越大，说明群落多样性越低.

测序分析结果如表3所示. 由表3可见，所有样品Coverage均在0.998以上，表明测序样品微生物覆盖率高，满足菌群多样性分析需求. 沿程微生物多样性分析表明，Shannon-Wiener指数基本呈逐渐上升趋势，表现为A

表3 水样Alpha多样性分析

2.2 系统沿程微生物群落结构变化情况

为明确系统沿程微生物群落结构变化情况，开展各单元出水样品组间相似性分析及群落组成变化分析. 由相似性分析(analysis of similarities, ANOSIM)结果可知，不同处理单元出水群落结构具有显著性差异(见图2, Bray-Curtis ANOSIM=0.388 1,P<0.01)，表明处理工艺的不同对微生物群落结构影响较大.

注： Ave表示组间距离排名的平均值.

系统各单元出水门水平微生物菌群结构分布情况如图3所示. 进水中变形菌门(Proteobacteria)相对丰度最高(90.55%)，与Shen等[24]研究结果一致. 变形菌门属革兰氏阴性异养细菌，是石油化工废水中最常见的微生物之一[25-26]，具有高效降解石油烃化合物的能力. 该研究中，Proteobacteria在整个处理过程中亦占主导地位.

图3 门水平上微生物群落结构组成

系统沿程各单元出水优势菌组成均发生变化，其中，调节池、涡凹气浮池、斜板气浮池出水中Proteobacteria相对丰度分别降至54.90%、62.59%、52.89%，Epsilonbacteraeota相对丰度分别增至17.03%、29.60%、41.53%，厚壁菌门(Firmicutes)相对丰度分别增至26.13%、6.02%、2.33%. 已有研究[27-29]表明，Epsilonbacteraeota和厚壁菌门具有深度降解废水中难降解有机物的能力. A/O池出水中Proteobacteria相对丰度降至39.69%，髌骨细菌门(Patescibacteria)和拟杆菌门(Bacteroidetes)相对丰度分别增至23.95%、15.65%. 其中，髌骨细菌门为难降解有机废水处理工艺中常见优势菌[30-31]；拟杆菌门具有较强的营养物质代谢能力[32]，为污水脱氮的重要参与者. 生物单元对系统TN和NH3-N去除贡献率(分别为72.202%和28.750%)均较高，表明拟杆菌门在该单元脱氮过程中发挥重要作用.

2.3 微生物群落多样性与环境因子响应分析

微生物群落多样性与环境因子密切相关[33-35]. 进行环境因子分析前，使用VIF(variance inflation factor analysis，方差膨胀因子分析)对3类环境因子指标(14种)进行环境因子筛选，保留多重共线性较小的环境因子进行后续研究. 筛选后的环境因子(VIF<10)分别为COD(2.954)、TN(4.983)、TP(1.769)、pH(1.491)、NH3-N(6.555)、Na(2.878)和Ca(2.339)(见表4).

表4 筛选后环境因子的VIF值

通过db-RDA分析(基于Bray-Curtis距离)揭示了样品与筛选后环境因子的相关性(见图4). 在db-RDA 二维图中，带箭头的线段表示环境因子，样点到数量型环境因子的箭头投影长度反映该环境因子对样品群落组成的影响. 结果表明，筛选后的环境因子均与A/O出水呈负相关，COD、NH3-N、TN浓度均与进水微生物群落结构呈正相关. 其中，COD(r2=0.599,P<0.05)、NH3-N(r2=0.574,P<0.05)浓度是影响废水中微生物群落变化的最显著因素.

图4 环境因子与微生物群落基于db-RDA分析的相关性分析结果

为明确炼油废水处理效果与微生物群落结构的响应关系，选取显著相关的环境因子(COD、NH3-N)去除率与门水平优势物种(相对丰度前10)进行相关性分析(见图5). 结果表明，NH3-N去除率与Bacteroidetes (r=0.588,P<0.05)呈显著正相关，COD去除率与Patescibacteria (r=0.530,P<0.05)和Bacteroidetes (r=0.706,P<0.01)均呈显著正相关，这与陈重军等[36-37]研究结果一致.

注： *表示0.01≤P≤0.05，** 表示0.001≤P≤0.01.

2.4 微生物代谢功能预测分析

基于宏基因组测序结果，采用PICRUST软件与KEGG基因数据库[38-39](Kyoto encyclopedia of genes and genomes, http://www.genome.jp/kegg)和通路数据库进行对比分析，预测各样品功能基因组成情况. KEGG通路数据库主要包括6类功能基因，分别为新陈代谢、遗传信息处理、环境信息处理、细胞过程、人类疾病和机体系统. 该研究中，各样品中相对丰度最高的功能基因为新陈代谢功能基因，占比均高于60%(A、B、C、D、E样品分别为68.37%、63.67%、63.20%、63.97%、62.59%). 其中，代表有机物去除的新陈代谢功能基因为能量代谢和外源化合物生物降解[40]，因此将其作为该研究的重点.

从能量代谢角度分析，沿程各单元能量代谢基因相对丰度大小表现为A>C>B>D>E(见图6)，推测由于易降解物质逐级耗尽导致后续单元能量代谢基因相对丰度降低. 从外源化合物生物降解角度分析，一级处理工艺相对丰度波动较小，均在3.00%以下；进入生物处理单元后，其相对丰度迅速提升至15.38%，其中，注释到石油烃降解代谢通路的功能基因相对丰度显著增加，如苯甲酸降解基因(相对丰度为3.15%)、二甲苯降解基因(相对丰度为1.18%)、硝基甲苯降解(相对丰度为2.15%)等. 另外，生物处理单元对系统各主要污染物去除贡献率较大(见表1)，进一步证实炼油废水一级处理工艺无法完全去除废水中难降解有毒有害组分，有毒物质的长期累积可能增加生物处理单元冲击风险. 据此，亟待开展点源炼油废水处理系统沿程废水生物毒性削减负荷核算及关键致毒物质识别工作.

图6 KEGG外源生物降解和能量代谢

3 结论

a) 系统沿程微生物多样性Shannon-Wiener指数和丰富度Ace指数基本呈逐渐上升趋势，表现为进水<涡凹气浮池<斜板气浮池<调节池

b) ANOSIM分析揭示了系统沿程微生物群落结构差异显著，各单元出水优势菌组成均发生变化. 其中，调节池、涡凹气浮池、斜板气浮池出水中Proteobacteria占比分别由90.55%降至54.90%、62.59%、52.89%，Epsilonbacteraeota占比分别增至17.03%、29.60%、41.53%，厚壁菌门(Firmicutes)占比分别增至26.13%、6.02%、2.33%；A/O池出水中Proteobacteria占比降至39.69%，髌骨细菌门(Patescibacteria)和拟杆菌门(Bacteroidetes)占比分别增至23.95%、15.65%.

c) 环境因子与微生物群落相关性分析表明，COD、NH3-N在系统微生物群落形成过程中起重要作用(P<0.05)，其中，NH3-N去除率与Bacteroidetes呈显著正相关(P<0.05)，COD去除率与Patescibacteria(P<0.05)、Bacteroidetes(P<0.01)均呈显著正相关.

d) KEGG功能注释分析可知，生物处理后外源化合物生物降解相关基因占比迅速增加，其中，与石油烃降解代谢通路相关的功能基因丰度增加尤为明显.