PD-Anammox快速启动过程中氮代谢功能基因的表达特征

刘 宁，何 阳，王德朋，杨 庆，邱晨晨，罗 琦，，黄开龙，，张徐祥

（1.南京大学金陵学院，江苏南京 210032；2.南京大学环境学院，污染控制与资源化研究国家重点实验室，江苏南京 210023；3.南京江岛环境科技研究院有限公司，江苏南京 210019）

短程反硝化-厌氧氨氧化（PD-Anammox）耦合工艺通常以氨态氮和硝酸盐氮作为基质富集培养厌氧氨氧化菌和反硝化菌，从而实现耦合工艺的启动〔1〕。受厌氧氨氧化菌倍增周期长以及耦合工艺中有机物对厌氧氨氧化细菌富集的影响〔2〕，PD-Anammox耦合工艺的常规启动过程十分漫长〔1〕。近年有研究发现，通过调控进水基质比例，将耦合工艺初始进水氮源由硝酸盐改为亚硝酸盐，待脱氮效率稳定后逐步将亚硝酸盐氮转换为硝酸盐氮，可以有效减轻进水基质（特别是进水中的有机物）对厌氧氨氧化菌的影响，从而实现耦合工艺的快速启动〔3〕。需要指出的是，进水基质的变化可显著改变微生物的代谢功能〔4〕。在长期驯化过程中，厌氧氨氧化系统中的微生物群落结构及其功能均呈现出明显的演替过程〔5-8〕。目前，关于PD-Anammox耦合工艺中微生物菌群演替及其氮代谢功能基因表达随基质改变的瞬时响应机制的研究较少，而相关机理的解析对耦合工艺的快速启动和稳定运行具有重要意义。

笔者拟通过切换进水的氮源基质类型（亚硝酸盐氮完全转化为硝酸盐氮）实现耦合工艺的快速启动，重点揭示启动过程中微生物功能菌群的变化规律，并基于宏转录组学技术深入解析基质改变时功能菌群的瞬时响应机制，阐明耦合工艺中脱氮效能恢复的微生物学机理，为实现PD-Anammox耦合工艺的快速启动提供技术支撑和理论依据。

1 材料与方法

1.1 实验装置与接种污泥

实验装置为膨胀污泥床反应器（EGSB）（见图1），有机玻璃材质，内径5 cm，有效容积1 L。运行过程中，用恒温槽（HX-08，无锡沃信仪器有限公司）保持反应器内水温为（30±1.5）℃，水力停留时间（HRT）为8 h，同时设置回流比为10∶1。接种污泥取自实验室运行的厌氧氨氧化-反硝化耦合反应器〔9〕，污泥接种量为200 mL，SS为11.3 g∕L。

图1 EGSB耦合反应器Fig.1 EGSB coupled reactor

1.2 PD-Anammox耦合工艺的快速启动

采用EGSB反应器及模拟污水，通过改变进水基质比例实现PD-Anammox耦合工艺的快速启动。该过程分为3个阶段（阶段Ⅰ进水氮素为氨氮和亚硝酸盐氮，阶段Ⅱ进水氮素为氨氮、亚硝酸盐氮和硝酸盐氮的混合液，阶段Ⅲ进水氮素为氨氮和硝酸盐氮），运行时间为77 d，具体运行情况见表1。耦合工艺运行过程中严格控制反应器内p H在7.0～8.0，模拟污水进水碳源的配制参照文献〔10〕。为保证微生物正常生长，在1 L进水中添加1 mL微量元素储备液。运行期间，每48 h采集1次出水样品（3个平行），并用0.45μm滤头过滤，统一保存在4℃冰箱中，用于水质检测。污泥样品每次采集3个平行样品，当天提取污泥样品中的RNA送往测序公司，同时在收集的污泥样品中加入等体积的无水乙醇，置于-20℃冰箱中用于后续提取样品总DNA。

表1 PD-Anammox耦合工艺启动过程中进水水质及污泥采样信息Table 1 Influent quality and sludge sampling information during the start-up of PD-Anammox coupled system

1.3 分析方法

氨氮采用纳氏试剂分光光度法测定；亚硝态氮采用N-（1-萘基）-乙二胺分光光度法测定；硝态氮采用紫外分光光度法测定；总氮采用碱性过硫酸钾消解紫外分光光度法测定；COD采用重铬酸盐法〔11〕测定。

PD-Anammox耦合工艺中总氮去除的贡献率按式（1）〔12〕计算，全程反硝化贡献率按式（2）计算。

式中：（NH4+-N）inf——进水氨氮，mg∕L；

（NH4+-N）eff——出水氨氮，mg∕L；

TNinf——进水总氮，mg∕L；

TNeff——出水总氮，mg∕L。

1.4 高通量测序与生物信息学分析

1.4.1 16SrRNA基因高通量测序与数据分析

采集的污泥样品用FastDNA®Spin Kit土壤试剂盒（MP Biomedicals，美国）提取总DNA，将DNA样品送往上海美吉生物医药科技有限公司，采用ABIGeneAmp®9700型PCR仪扩增16SrRNA基因的序列V3-V4区域〔13〕。PCR扩增结果经纯化和建库后上机测序。原始数据通过qiime2软件合并数据，再通过质量筛选、模块降噪后生成特征表和代表序列，最终获得物种分类信息。

1.4.2 宏基因组测序与数据分析

对DNA样品进行宏基因组测序，每个样品的测序深度为6 GB。DNA数据处理流程如图2所示。

图2 DNA数据处理流程Fig.2 Flow of DNA data processing

1.4.3 宏转录组数据分析

污泥样品中的RNA采用全式金试剂盒（TransZol Up Plus RNA Kit）提取，将RNA样品置于干冰中寄送至上海美吉生物医药科技有限公司进行宏转录组测序，每个样品的测序深度为6 GB。RNA原始数据处理流程如图3所示。

图3 RNA原始数据处理流程Fig.3 Flow of raw RNA data processing

1.4.4 生物统计学分析

挑选LDA绝对值＞4.0物种进行统计分析，在线 分 析 网 址http：∕∕huttenhower.sph.harvard.edu∕galaxy∕。利用T-检验对水质指标进行差异性分析，再用R语言EasyAovWlxPlot程序包对菌群丰度和氮代谢基因进行方差分析，将符合正态分布、方差齐性＞0.05的平行数据用Tukey检验进行多重比较分析，不同字母表示具有显著性差异（p＜0.05）。

2 结果与讨论

2.1 PD-Anammox耦合工艺的快速启动

PD-Anammox耦合工艺快速启动过程中，出水氮素浓度及总氮去除率变化情况如图4所示。

图4 耦合反应器出水氮素及其去除率变化Fig.4 Variations of effluent nitrogen and their removal rates in the coupled system

阶段Ⅰ中，进水中氨氮为30 mg∕L，亚硝酸盐氮为30 mg∕L，硝酸盐氮为0。从开始至运行20 d稳定后，总氮去除率从最初的（73.7±2.4）%升至（90.9±1.6）%，相应的总氮负荷去除速率由0.190 kg∕（m3·d）增至0.241 kg∕（m3·d）。出水中亚硝酸盐氮和氨氮能被完全去除，存在少量硝酸盐氮累积（硝酸盐氮平均值2.9 mg∕L），此时反应器中Anammox的总氮去除贡献率为100%。

阶段Ⅱ中，进水亚硝酸盐氮降至16 mg∕L，硝酸盐氮增至16 mg∕L，相应的碳源质量浓度由20 mg∕L提升至40 mg∕L。运行前2天总氮去除率降至（77.4±0.6）%，出水中的硝酸盐氮为10.5 mg∕L，表明进水改变初期耦合工艺的反硝化活性较弱，导致硝酸盐累积。运行20 d后（阶段Ⅱ末期），去除率回升至（89.4±0.7）%，总氮去除速率为0.224 kg∕（m3·d），此时反应器中PD-Anammox总氮去除贡献率为100%。

阶段Ⅲ中，进水硝酸盐氮提高至33 mg∕L（COD为75 mg∕L）。进水第1天总氮去除率突降至66.4%，出水中氨氮、硝态氮和亚硝态氮分别为8.6、11.0、0.9 mg∕L，第2天总氮去除率上升至70.6%，第4天快速上升至76.6%；经30余天的驯化，总氮去除率回升至（89.1±0.5）%，出水中的氨氮、硝酸盐氮和亚硝酸盐氮分别为2.0、3.3、1.7 mg∕L，总氮负荷去除速率升至0.227 kg∕（m3·d）。此时PD-Anammox总氮去除贡献率为97.1%，全程反硝化贡献率为2.9%。综上，氮源基质的改变（亚硝酸盐氮完全转化为硝酸盐氮）在短期内显著降低了耦合工艺的脱氮效能，但长期运行后PD-Anammox耦合工艺与全程反硝化可达到新的平衡，从而恢复系统的脱氮效果。

2.2 氮源冲击过程微生物群落结构变化特征

通过16SrRNA基因高通量测序技术获得耦合工艺中微生物群落结构受氮源改变瞬时冲击（阶段Ⅱ切换至阶段Ⅲ）时的响应特征，结果如图5所示。

图5 瞬时冲击下耦合反应器门水平微生物的群落结构变化Fig.5 Variations of microbial community structure（phylum level）in coupled system under transient impact

在阶段Ⅱ稳定期，变形菌门（Proteobacteria，38.0%～63.6%）、浮霉菌门（Planctomycetes，10.7%～22.5%）和绿弯菌门（Chloroflexi，8.7%～13.2%）是耦合工艺中的主要菌门。变形菌门包含广泛的反硝化菌〔14〕，当进水硝酸盐和碳源增加2 h（S1）后，变形菌门相对丰度从38.0%（S0）显著升高至59.3%（p＜0.05），并 在 进 水1 d后（S2）丰 度 达 到 最 大 值（63.6%），但随后其丰度在第3天（S3）时显著降低至42.0%，与阶段Ⅲ稳定阶段（S4，42.2%）相比并无显著差异（p＞0.05）。进水硝酸盐氮含量增加可诱导变形菌门丰度升高，表明进水硝态氮的冲击短时间内诱导了耦合系统中反硝化菌群的快速增殖〔5〕。浮霉菌门丰度在进水2 h后降至11.9%，并在1 d后达到最低丰度（10.7%）；3 d后（S3）其丰度升至16.7%，在稳定期丰度为12.1%。进水改变瞬间，碳源浓度增加对厌氧氨氧化菌产生了抑制，导致浮霉菌门丰度显著下降〔15〕。随着厌氧氨氧化菌对进水条件的适应，其丰度逐渐上升，但稳定期后碳源升高导致异养菌滋生过多，造成厌氧氨氧化菌相对丰度的降低。与此同时，异养菌Calditrichaeota门的丰度在稳定运行期间显著高于氮源冲击前。该菌门可利用细胞裂解的碳源和蛋白质进行富集生长〔16〕，是运行后期浮霉菌门相对丰度下降的主要原因。

为进一步了解耦合工艺快速启动过程中微生物群落结构的变化情况，在属水平对微生物进行LEfSe分析，以获得氮源冲击下的差异微生物信息，如图6、图7所示。

如图7所示，S0～S4样品中丰度差异最大的菌属分别是陶厄氏菌（Thauera）、动胶菌（Zoogloea）、球衣菌（Sphaerotilus）、Denitratisoma和JdFR_76。从进化树分支（图6）可知，γ-变形菌纲（Gammaproteobacteria）包含的5个特异性菌属分别为陶厄氏菌（Thauera）、动胶菌（Zoogloea）、球衣菌（Sphaerotilus）、Denitratisoma和固氮螺菌（Azospira），其中陶厄氏菌（Thauera）〔17〕、动 胶 菌（Zoogloea）〔18〕和Denitratisoma〔19〕都属于反硝化菌，表明γ-变形菌亚门内反硝化菌差异物种在瞬时冲击后的富集显著。此外，Calorithrix（2.4%～5.8%）是耦合工艺启动后的差异菌，其隶属的Calditrichaeota门丰度也在进水冲击后显著上升。Calorithrix属中包含一种新型细菌Calorithrix insularissp〔20〕，能 异 化 硝 酸 盐 还 原 为 铵（DNRA），推测在阶段Ⅲ中部分硝酸盐被异化还原为铵盐。

图6 瞬时冲击下耦合反应器微生物群落结构LEfSe分析Fig.6 LEfSe analysis of microbial community structure in coupled system under transient impact

图7 瞬时冲击下耦合反应器微生物群落结构LDA分析Fig.7 LDA analysis of microbial community structure in coupled system under transient impact

为进一步确定脱氮功能菌在氮源冲击下的响应特征，考察冲击后的丰度变化情况，如图8所示。

图8 瞬时冲击下耦合反应器脱氮功能菌群的变化情况Fig.8 Variations of nitrogen removal functional bacteria in coupled system under transient impact

随着进水中硝酸盐和碳源含量的增加，耦合工艺中的厌氧氨氧化菌Candidatus Brocadia丰度由21.9%（S0）降为11.2%（S1），1 d后（S2）达到谷值（9.5%），其原因在于有机物浓度突升对厌氧氨氧化菌产生一定抑制作用〔21〕。3 d后（S3）Candidatus Brocadia丰度回升至16.1%。该趋势与总氮去除率的变化趋势一致。在阶段Ⅲ稳定运行时期（S4），厌氧氨氧化菌丰度相对于氮源冲击第3天后（S3 16.1%）降至11.0%，原因可能与长期运行下异养菌（如Calorithrix）的滋生有关。

相对丰度较高的反硝化菌分别为Denitratisoma（2.3%～7.6%）、动 胶 菌（2.2%～6.4%）和 陶 厄 氏 菌（2.3%～5.2%）。其中内源性反硝化菌Denitratisoma的相对丰度在瞬时冲击2 h后下降为2.3%，并在1 d（3.4%）和3 d（7.6%）后迅速上升，在阶段Ⅲ稳定运行时期（S4），反硝化菌Denitratisoma丰度由7.6%降至3.9%。Denitratisoma主要利用细胞死亡裂解的有机碳源富集生长〔22〕。受到氮源冲击后，部分细胞死亡裂解促进了Denitratisoma的生长，该菌是厌氧氨氧化菌重要的共生细菌〔23-24〕。与此同时，陶厄氏菌是多种PD-Anammox耦合工艺中的优势短程反硝化菌〔8，25-26〕，其丰度在瞬时冲击2 h后也呈现先下降（2.3%）后上升（3.5%）的趋势，表明Denitratisoma和陶厄氏菌是氮源冲击后耦合工艺恢复阶段的关键反硝化菌。

2.3 氮源冲击过程对耦合工艺氮代谢基因差异表达的影响

利用宏转录组学技术分析氮代谢功能基因表达特征，揭示耦合工艺受氮源冲击后的瞬时与长期响应，以及氮源改变前后的基因变化规律，如图9所示。

图9 瞬时冲击下耦合反应器氮代谢功能基因的表达量（TPM）Fig.9 Expression of nitrogen metabolic genes in coupled system under transient impact

图9中，当进水硝酸盐和碳源升高2 h后（S1），厌氧氨氧化反应的主要功能基因hdh∕hzo和hzsA∕B∕C表达丰度分别从阶段Ⅱ运行末期（S0）的38 652和54 505 TPM显著下降至5 537和7 687 TPM（p＜0.05），1 d后（S2）开始显著上升（p＜0.05），同时反应器中的Candidatus Brocadia菌呈现明显富集趋势。随着反应时间的延长，耦合工艺脱氮效能逐渐恢复，hdh∕hzo和hzsA∕B∕C基 因 表 达 量 无 显 著 变 化（p＞0.05），说明厌氧氨氧化反应活性自第3天开始趋于稳定。阶段Ⅱ运行末期（S0）的hdh∕hzo和hzsA∕B∕C基因表达活性下降，主要是由于进水中的亚硝酸盐氮减半，硝酸盐难以被厌氧氨氧化菌直接利用，且碳源升高对厌氧氨氧化菌有抑制作用〔27〕。此外，稳定时期（S4）hdh∕hzo和hzsA∕B∕C基因表达丰度显著低于进水改变前（S0），主要原因在于进水中硝酸盐和碳源浓度的升高分别促进了反硝化细菌和异养菌的生长与活性〔28〕。

硝酸盐还原基因napA∕B∕C和narB∕C、亚硝酸盐还原基因nirS和一氧化二氮还原酶基因nosZ的表达量在进水硝酸盐和碳源升高2 h后（S1）下降明显，氮源冲击1 d后（S2），反硝化相关功能基因napA∕B∕C、nirK、nirS及nosZ表达量显著上升（p＜0.05），S2阶段反硝化相关功能基因表达量升高，其反硝化活性增强，系统中的有机物被反硝化菌利用，进而减少了有机物对厌氧氨氧化菌活性的抑制。另外，narG∕H基因表达受到氮源冲击后（S1～S3）表达丰度没有显著差异，有研究显示narG基因丰度受NO3--N浓度变化的影响不显著〔29-30〕。随着反应时间的延长，硝酸盐还原基因中除napA∕B∕C基因表达量显著下降，narB∕C∕G∕H∕I∕J∕Y∕Z基因表达量均呈显著上升趋势，而nosZ基因表达量无显著变化。表明在长期适应过程中，细胞周质中的硝酸盐还原活动减弱，细胞内硝酸盐还原活动相对上升，而一氧化氮还原为氮气的氮代谢过程并未加强。与进水改变前（S0）比较，稳定时期（S4）的硝酸盐还原酶基因narG∕H、narI∕J、narY∕Z和异化硝酸盐还原为铵的关键基因nrfA∕B∕C∕D均呈显著上升趋势，铵化的相关功能菌（Calorithrix）的丰度（5.8%）与氮源冲击前的丰度（2.5%）相比也明显上升，表明在稳定阶段Calorithrix在短程反硝化过程中起到重要作用，narG基因已被证实在短程反硝化中发挥关键作用〔31〕，高浓度硝酸盐和碳源能促进异化硝酸盐为铵的反应活性〔32〕。

2.4 氮源冲击下耦合工艺氮代谢通路的变化趋势

基于宏转录组学结果，构建了氮源冲击条件下耦合工艺中氮代谢活性通路图，如图10所示。

图10 瞬时冲击下氮代谢通路Fig.10 Nitrogen metabolic pathways under transient shock

由基因表达活性可知，活性较高的氮代谢途径为厌氧氨氧化作用和反硝化作用。氮源瞬时冲击（2 h）会对上述2种氮代谢途径产生抑制，其中厌氧氨氧化功能基因hdh∕hzo和hzsA∕B∕C，反硝化功能基因napA∕B∕C、norB∕C和nosZ均显著下降（p＜0.05），总氮去除率降至最低（66.4%）。氮源冲击1 d后，厌氧氨氧化功能基因hdh∕hzo和hzsA∕B∕C，反硝化功能基因napA∕B∕C和nosZ表达量逐步上升（p＜0.05），3 d后其基因表达量升至峰值，运行末期厌氧氨氧化功能基因hdh∕hzo和hzsA∕B∕C表达量丰度降低且保持稳定，反硝化作用功能基因napA∕B∕C和nosZ表达量丰度升高且趋于稳定。

3 结论

（1）硝酸盐氮冲击在短期内显著降低了耦合工艺的脱氮效能，但长期运行后系统中的PDAnammox与全程反硝化反应可达到新的平衡，进而恢复系统的脱氮效能。

（2）进水硝酸盐氮含量升高导致变形菌门丰度升高，碳源浓度的增加导致浮霉菌门丰度显著下降。γ-变形菌门内的不同反硝化菌种在受到氮源瞬时冲击后富集情况差异显著，Calorithrix为耦合工艺启动后最主要的差异菌属。Candidatus Brocadia为氮源冲击后耦合工艺恢复效能的关键厌氧氨氧化功能菌，Denitratisoma和陶厄氏菌为氮源冲击后耦合工艺恢复效能的关键反硝化功能菌。

（3）氮源冲击会抑制耦合工艺中的厌氧氨氧化和反硝化氮代谢途径。其中，厌氧氨氧化功能基因hdh∕hzo和hzsA∕B∕C，硝 酸 盐 还 原 基 因napA∕B∕C和narB∕C，亚硝酸盐还原基因nirS、一氧化二氮还原酶基因nosZ表达量均显著下降；氮源冲击引起的抑制作用在1 d后开始减弱，hdh∕hzo和hzsA∕B∕C，napA∕B∕C和nosZ基因表达量显著上升；系统恢复脱氮效能稳定运行后，与冲击前相比，hdh∕hzo和hzsA∕B∕C基因表达丰度降低，nosZ基因表达量升高。