近年来, 人为来源活性氮的过量排放严重干扰了河口生态系统的氮循环过程, 由此引发了水体富营养化等系列环境问题.因此, 河口近岸生态系统氮循环过程及其功能微生物类群成为广泛关注的研究热点.硝化过程(将氨氮(NH$_{4}^{+})$氧化为亚硝氮(NO$_{2}^{-})$, 再将NO$_{2}^{-}$氧化为硝氮(NO$_{3}^{-})$)作为氮循环的重要环节, 在氮的生物地球化学循环中起着重要作用^[1-2].硝化微生物是参与硝化反应、促进环境中氮元素迁移转化的重要微生物.硝化微生物包括氨氧化细菌(ammonia oxidizing bacteria, AOB)、氨氧化古菌(ammoniaoxidizing archaea, AOA)、亚硝酸盐氧化菌(nitrite oxidizing bacteria, NOB)以及新发现的能进行NH$_{4}^{+}$氧化和NO$_{2}^{-}$氧化两个过程的全程氨氧化微生物(complete ammonia oxidizers, comammox)^[3-4].

目前, 有关自然环境中AOB与AOA的菌群多样性、空间分布差异等方面的研究较多^[5-7], 而NOB的研究相对滞后.一方面, 因为亚硝酸盐氧化不是硝化过程的限速步骤, 所以亚硝酸盐氧化过程的关注度较低; 另一方面, 因为NOB Nitrospira属难培养, 发现新种系的可能较小, 所以NOB的相关研究也较少^[8].而comammox的发现更晚于AOB、AOA和NOB类群, 并迅速成为氮循环过程的研究热点.首次报道的三株comammox(Candidatus Nitrospira inopinata, Candidatus Nitrospira nitrosa和Candidatus Nitrospira nitrificans)均发现于工业系统^[3-4], 之后在饮用水系统、活性污泥、森林土壤和湖泊沉积物等生态系统都检测到comammox的存在^{[3-4, 9-16]}.

河口近岸区域是地球上最多产的生态系统之一, 是生物地球化学循环的重要组成部分, 是氮循环研究的热点区域^[17-18].目前, 在河口潮滩生态系统中对硝化功能微生物还缺乏深入认识, 尤其是新发现的comammox类群, 其生态位及与其他硝化功能菌的相互作用还不甚清楚.受潮汐变化、盐水入侵以及脉冲式营养盐输入等海陆交互作用的影响, 河口沉积环境微生物群落在系统分类和功能代谢上较为复杂^[19].富集培养是微生物探究的重要方法, 可以放大研究目标, 简化研究体系^[20].为获得河口潮滩生态系统具有较高硝化能力的微生物群体, 提供后续硝化微生物代谢和影响机制探究的实验基础, 简化河口潮滩复杂环境因子的影响, 本研究进行了硝化微生物的富集培养.在此基础上, 结合宏基因组和宏转录组技术, 从DNA和mRNA两个层面对硝化富集培养物进行了种群结构和基因表达活性的分析.研究工作深化了对河口潮滩氮素生物地球化学循环的理论认识, 对于深入开展河口环境硝化过程的分子生态学研究有重要借鉴.

为了避免复杂环境因子的干扰, 更具体地分析河口潮滩地区的硝化过程, 选择氨氧化(硝化过程的限速步骤)细菌的最适生长条件, 进行了硝化微生物的富集培养.将4.0 L工作容积的生物反应器湿热灭菌后, 接种50 g新鲜的长江口潮滩沉积物, 并配备pH、温度、溶解氧和液面高度监测的探针.设置培养基的加入流速1.4 mL/min, 搅拌速度为每分钟250转.通过加入100 g/L的KHCO$_{3}$溶液维持反应体系的pH=7.49$\pm $0.09.控制反应器的温度$T= (26.9\pm 0.8)$℃, 反应体系溶解氧DO=(94.3$\pm $3.6)%(相对于水体饱和溶解氧).当反应体系液面超过4.0 L, 液体经膜过滤后泵出.培养基由纯水配置, 包含0.147 g/L CaCl$_{2}\cdot $2H$_{2}$O, 1.72 g/L NaCl, 0.58 g/L MgSO$_{4}\cdot $7H$_{2}$O, 0.054 g/L KH$_{2}$PO$_{4}$, 1 mL/L微量元素溶液1和2.培养基的NH$_{4}$Cl浓度由9 mg/L逐渐增加到36 mg/L.每升微量溶液1含有5 g EDTA和5 g FeSO$_{4}$.每升微量溶液2含有15 g EDTA, 0.1 g ZnSO$_{4}\cdot $7H$_{2}$O, 0.002 g CoCl$_{2}\cdot $6H$_{2}$O, 0.2 g MnCl$_{2}\cdot $4H$_{2}$O, 0.02 g CuSO$_{4}\cdot $5H$_{2}$O, 0.1 g Na$_{2}$MoO$_{4}\cdot $2H$_{2}$O, 0.05 g H$_{3}$BO$_{4}$, 0.1 g NiCl$\cdot $6H$_{2}$O.培养基及反应器中NH$_{4}^{+}$、NO$_{2}^{-}$和NO$_{3}^{-}$浓度使用Skalar营养盐自动分析仪(SAN plus)分析, 其中NH$_{4}^{+}$采用靛酚蓝法, NO$_{2}^{-}$和NO$_{3}^{-}$采用重氮偶氮法^[21].

为分析硝化菌群$\beta $-AOB($\beta $-proteobacterial ammonia-oxidizing bacteria)和AOA在培养过程中的丰度差异, 提取不同培养时期(培养天数分别为0, 5, 10, 15, 20, 25, 30和35 d)反应器富集物的总DNA(PowerSoil®DNA Isolation Kit(MO BIO, USA)试剂盒), 使用ABI Prism 7500 real-time PCR system(Applied Biosystems, Canada)进行定量. 20 $\mu $L反应体系包含10 $\mu $L 2$\times $TransStart Top Green qPCR SuperMix (Transgene Biotech, Beijing, China), 7.5 $\mu $L双蒸水, 0.4 $\mu $L Passive Reference Dye Ⅱ (50$\times)$, 0.4 $\mu $L上下游引物(10 mmol L$^{-1})$, 1$\mu $L DNA模板($1\sim 20$ ng).扩增程序为: 94 ℃预变性5 min, 94${^\circ}$变性5 s, 退火温度15 s, 72 ℃延伸34 s并采集信号, 共40个循环. $\beta $-AOB定量的上下游引物为amoA1F (5$'$-GGGGTTTCTACTGGTGGT-3$'$)和amoA2R(5$'$-CCCCTCKGSAAAGCCTTCTTC-3$'$), 退火温度为57 ℃^[22]. AOA定量实验的引物为Arch-amoAF (5'-STAATGGTCTGGCTTAGACG-3')和Arch-amoAR(5'-GCGGCCA- TCCATCTGTATGT-3'), 退火温度为56℃^[23].携带有细菌及古菌amoA基因片段的质粒采用Plasmid Mini Preparation试剂盒(Tiangen)从本研究构建的大肠杆菌中抽提出来, 用于标准曲线的构建.质粒DNA浓度用Nanodrop-2000(Thermo)超微量分光光度计测定.本研究中构建的$\beta $-AOBamoA基因标准曲线的浓度范围是42.1$\sim $4.21$\times $10$^{5}$ copies/ng DNA, AOB amoA基因标准曲线的浓度范围是22.4$\sim $2.24$\times $10$^{5 }$ copies/ng DNA.临界循环次数(threshold cycle, CT)与标准质粒拷贝数的对数值呈现显著的线性相关($R^{2}> 0.99$), AOB与AOA amoA基因的扩增效率大于90 %.不加DNA模板的体系作为空白, 所有试验均3个平行.

经一个月的富集培养后, 使用PowerSoil®DNA Isolation Kit (MO BIO, USA)提取3份反应器富集物的总DNA.将3份DNA样品混匀, 使用Covaris S220 Focused Ultrasonicator将总DNA随机打断成350 bp的片段.通过PCR反应对这些片段进行末端修复, 连接poly-A尾并添加接头构建测序文库, 基于Illumina Hiseq平台, 进行150-bp双末端测序.得到的原始数据进入信息分析流程.首先去除原始测序reads的接头, 去除threshold quality$\le $38, 长度$\le $40 nt, 含N碱基$\ge $10的3种reads.使用MEGAHIT^[24]组装软件对高质量的有效数据进行组装(参数: -presets meta-large), 组装片段从N碱基连接处打断, 保留长度在500 bp以上的片段^[25-26].采用MetaGeneMark^{[25, 27]}进行开放阅读框(open reading frames, ORFs)的预测(参数:采用默认参数), 并从预测结果出发, 删除长度小于100 nt^[28]的片段.采用CD-HIT^[26]软件进行去冗余, 以identity 95%, coverage 90%的标准聚类, 选取最长的序列为代表序列(参数: -c 0.95, -G 0, -aS 0.9, -g 1, -d 0^[29]).采用SoapAligner^[30]将样品的有效数据比对至初始基因集, 计算每个基因在样品中比对上的reads数目(比对参数: -m 200, -x 400, identity $\ge \;95 \%$^[30]), 过滤掉在样品中支持reads数目$\le $2的基因^[31], 获得用于后续分析的非冗余基因集(Unigene catalogue).如计算得到基因$k$在样品中的丰度信息$G_{k}$, 计算公式^[32]为

上式中, $r$为比对上基因$k$的reads数目, $L$为基因的长度.使用DIAMOND软件^[33]将Unigenes与NCBI的NR数据库比对(BLASTP, e-value$\le $10$^{-5})$.选取e-value$\le $10倍最小e-value^[34]的比对结果, 采用MEGAN软件, 使用LCA(Lowest Common Ancestor)算法进行系统分类注释^[35], 获得各个样品在各个分类层级(界门纲目科属种)上的丰度信息.本研究中测序reads在NCBI数据库中的序列号为SRP122254.

经一个月的富集培养后, 使用EZNA®Soil RNA kit (Omega Bio-tek, Norcross, GA, USA)试剂提取三份反应器富集物的RNA.使用Ribo-Zero Magnetic Kit (Epicenter, Charlotte, NC, USA)去除rRNA, 混合三份样品, 使用TruSeq$^{\rm TM}$ RNA Sample Prep Kit (Illumina, San Diego, CA, USA)试剂生成单链cDNA, 使用cBot Truseq PE Cluster Kit v3-cBot-HS试剂产生DNA簇, 在Illumina Hiseq平台测序.使用软件Seqprep(https://github.com/jstjohn/SeqPrep)对序列3$'$端和5$'$端进行质量剪切, 使用软件Sickle(https://github.com/najoshi/sickle)去除剪切后长度小于50 bp、平均质量值低于20以及含N碱基的reads.使用软件SortMeRNA(http://bioinfo.lifl.fr/RNA/sortmerna/)将reads比对SILVA SSU (16S/18S)与SILVALSU (23S/28S)数据库, 去除rRNA reads^[36].获得高质量reads后, 使用Trinity(http://trinityrnaseq.github.io/)将所有测序读段进行从头组装, 挑取长度$\ge $300bp的转录本, 利用TransGeneScan软件(http://sourceforge.net/projects/transgenescan/)对序列进行ORF预测^[37].用CD-HIT软件(http://www.bioinformatics.org/cd-hit/)将基因序列进行聚类(参数: 95% identity、90% coverage), 每个类取最长的基因作为代表序列, 构建非冗余基因集^[31].基因表达水平的比较采用FPKM(Fragments Per Kilobase of exon model per Million mapped reads)值进行衡量, 即每一百万条序列中, 每个基因以一千个碱基为单位, 比对上的reads条数^[38], 使用RSEM(http://deweylab.biostat.wisc.edu/rsem/)计算^[39].将基因集序列与eggNOG数据库进行比对(e-value$\le $10$^{-5})$, 获得基因对应的COG分类(Clusters of orthologous groups of proteins, 直系同源序列聚类).测序reads在NCBI数据库中的序列号为SRP122263.

一个月的培养过程中, 反应器体系各理化因子测量值变化趋势如图 1.培养基的NH$_{4}^{+}$浓度从0.5 mmol$\cdot$L$^{-1}$逐步上升至2 mmol$\cdot$L$^{-1}$.反应器中, NH$_{4}^{+}$浓度逐渐降低, 大约在20 d以后维持在0 mmol$\cdot$L$^{-1}$水平; NO$_{2}^{-}$起始为0 mmol$\cdot$L$^{-1}$, 一周左右上升至最高值, 随后在25 d左右降至0 mmol$\cdot$L$^{-1}$; NO$_{3}^{-}$浓度持续增高, 并在培养基中NH$_{4}^{+}$浓度上升时, 同步升高.反应器体系的pH在设置为7.45的条件下, 先升高, 最后逐步稳定为7.45, 推测由于氨氧化反应产生H$^{+}$增多^[4], pH呈下降趋势, 由碱性溶液调整, 维持pH稳定.液体中的溶解氧逐步降低, 可能为硝化过程消耗的O$_{2}$增加, 在30 d左右氧气的消耗量与外界的补充量达到平衡.整体表现为反应器中富集物硝化能力的逐渐增强.

图 1(Fig. 1)

图 1 培养基及反应器中NH$_{4}^{+}$、NO$_{2}^{-}$、NO$_{3}^{-}$浓度及pH、溶解氧变化趋势 Fig.1 The variation trend of NH$_{4}^{+}$, NO$_{2}^{-}$, NO$_{3}^{-}$ concentration and pH, dissolved oxygen in medium and enrichment culture

然而, $\beta $-AOB和AOA菌群的定量分析结果(见图 2)显示, 富集物中二者的丰度都呈现先升高后降低的变化趋势. $\beta $-AOB amoA的基因丰度从起始的2.1$\times $10$^{3}$ copies/ng DNA, 在培养20 d左右上升至最高水平, 为4.1$\times $10$^{4}$ copies/ng DNA, 之后在35 d左右降低至1.2$\times $10$^{3}$copies/ng DNA. AOA amoA的基因丰度从沉积物中的42.8 copies/ng DNA, 于培养20 d左右上升至最高水平, 为151.3 copies/ng DNA, 之后在35 d左右降低至48.7 copies/ng DNA.

图 2(Fig. 2)

注:图中左侧纵坐标表示每ng DNA中AOAamoA的基因拷贝数, 右侧是每ng DNA中$\beta$-AOBamoA的基因拷贝数 图 2 生物反应器中氨氧化微生物丰度 Fig.2 The abundance of ammonia oxidizing microorganisms in the bioreactor system

氨氧化微生物的丰度变化与反应器中硝化能力增强的趋势并不一致, 暗示反应体系中可能存在进行氨氧化作用的非传统微生物.与此同时, 新近comammox的发现指明反应体系中可能存在该类群微生物.

对富集物总DNA测序获得20 Gb数据, 构建的非冗余基因集包含762, 516个ORFs, 其平均长度为650.40 bp.统计获得菌群丰度结果如图 3, 图中显示了硝化类群($\beta $-AOB, AOA, Nitrospira)和其他微生物在种水平丰度由高到低排前5位的菌群.注释到的硝化类群占总测序reads的34.7%.

图 3(Fig. 3)

图 3 生物反应器中微生物类群丰度 Fig.3 The abundance of microorganisms in the bioreactor system

图 4显示硝化微生物群体的具体组成, 图中内圈扇形表示富集物在种水平上注释到的硝化菌群的丰度, 外圈表示富集物中4类硝化菌群comammox, NOB, $\beta $-AOB, AOA的丰度.硝化菌群中呈主导优势的微生物包括2种严格的NOB(占硝化菌群的50%左右)和3种comammox(占48%左右).注释到$\beta $-AOB有9种, 虽然种类相对较多, 但仅占硝化菌群的2%;注释到的AOA有3个类群, 丰度最低, 占0.1%的极小比例.

图 4(Fig. 4)

注:各类群丰度为该类群占硝化菌群的比例.图中内圈表示富集物在种水平上注释到的主要硝化菌群的丰度(右侧图例), 外圈表示富集物中4类硝化菌群的丰度(左侧图例) 图 4 生物反应器中硝化菌群丰度 Fig.4 The abundance of nitrifiers in the bioreactor system

宏基因组结果从种群结构角度揭示了富集的硝化微生物中comammox和NOB的优势地位. AOB、AOA和comammox三类菌群的竞争关系源于它们都通过对底物NH$_{4}^{+}$的氧化获得能量来源^[3-4], comammox的高丰度暗示了本研究培养条件下, comammox具有更强的竞争力.

宏基因组结果呈现了富集物的种群结构, 揭示了非传统氨氧化微生物——comammox的主导地位.结果暗示在长江河口潮滩沉积物中, 存在comammox菌群.硝化微生物的注释结果中, $\beta $-AOB和AOA类群的丰度低(约为2.1%)而微生物种类多(12种), NOB和comammox类群的丰度高(约为98%)而细菌种类少(5种).一方面可能由于富集培养的选择过程, 菌群多样性发生变化^[20], 注释结果表现了菌群的多样性差异; 另一方面, 由于研究难度等差异, 氨氧化微生物的相关研究更深入^[8], 而NOB和comammox的研究相对滞后, 所以数据库中氨氧化微生物的相关信息更丰富, NOB和comammox的信息较少, 注释结果呈现为comammox和NOB的低多样性.

生物反应器富集培养过程中, 设定培养条件(pH, 溶氧, 温度)为氨氧化(硝化过程的限速步骤)细菌的最适生长条件, 使培养基中NH$_{4}^{+}$浓度逐步升高, 以达到富集硝化菌群的目的.宏基因组数据验证了硝化菌群的富集成果, 表明了富集培养物较高的氨氧化(由comammox和氨氧化微生物承担)和亚硝酸盐氧化(由comammox和NOB承担)潜力.结合富集过程中培养条件的变化、氨氧化微生物降低的趋势、注释到comammox和NOB的高丰度, 硝化菌群结构的变化可能是培养基NH$_{4}^{+}$浓度的变化引起的.

对富集物总RNA反转录测序后获得14 Gb数据, 预测得到43 417个ORFs.比对eggNOG数据库将功能基因分为25个功能大类, 富集培养物4类硝化菌群的基因表达量在25个功能大类分布情况的统计结果如图 5. Comammox的代谢中, 翻译后修饰、蛋白质转换相关的基因表达量最高.能量生产及转化, 细胞壁、细胞膜合成, 信号传导以及翻译、核糖体构建、生物合成的基因表达次之, 暗示了富集物中comammox非常高的代谢活性和细菌增殖趋势. NOB基因表达量的高低分布情况与comammox相近, 但略低于comammox. AOA代谢中, 防御机制相关的基因表达量最高, 且远高于其他3个菌群, 翻译、核糖体构建、生物合成等代谢次之.防御机制相关基因表达量高的现象可能是生物反应器选择了AOB的富集培养条件(与富集培养AOA的pH、培养基成分等条件不同^[40-41]), 不适合的培养环境刺激产生了防御机制.氨氧化细菌在翻译, 核糖体构建, 生物合成和无机离子传输、代谢方面基因表达量较高.氨氧化微生物类群细胞壁、细胞膜合成代谢的基因表达量非常低, 暗示富集物中氨氧化微生物较小的增殖空间.

图 5(Fig. 5)

注:图中NOB指严格亚硝酸盐氧化的NOB, COM指全程氨氧化微生物, AOA指氨氧化古菌, $\beta $-AOB指$\beta $-变形菌门的氨氧化细菌 图 5 生物反应器中硝化菌群表达量热图 Fig.5 Heatmap of gene expression of NOB (strictly nitrite oxidizing Nitrospira), COM (comammox), AOA (ammonia oxidizing archaea), and $\beta $-AOB ($\beta $-proteobacterial ammonia oxidizing)

4类菌群预测性的和未知功能的基因表达量都较高, comammox的最高, 暗示菌群代谢中较多的知识空白.总体上, 4类菌群中, comammox基因表达量略高于NOB, 二者远高于氨氧化微生物.宏基因组数据中comammox与NOB的高丰度优势与宏转录组中二者的高表达量优势一致.结果表明本研究培养条件下, comammox和NOB具有较高的代谢活性, 与comammox在极贫营养条件下更有竞争优势的特点^[42-43]一致. Nitrospira在实验室难培养^{[8, 44]}的特点暗示本研究的富集培养对于Nitropira属和comammox的研究具有重要意义.

本研究采集河口潮滩沉积物, 进行硝化微生物的富集培养.宏基因组结果揭示了硝化微生物的种群结构, 主要包括4类硝化菌群——comammox, NOB, $\beta $-AOB和AOA.富集物催化氨氧化反应的主导微生物并不是传统的AOB和AOA, 而是最新发现的comammox.宏转录组结果进一步证实了高丰度的comammox具有高的基因表达量, 暗示了其较高的代谢活性.本研究所运用的富集培养的方法, 有利于后续推进对氮循环微生物驱动机制的探索.研究结果也初步呈现了硝化菌群之间的竞争关系, 对于探讨自然环境中的硝化过程具有重要意义.后续研究可以增加反应器数量, 在设置多个培养条件的基础上进一步探究菌群变化特征、菌群对某些特殊环境因子的响应以及响应机制.