基于宏基因組學的酸性森林土壤氨氧化微生物群落特征研究*

唐修峰，秦 華，匡 璐，王欣欣，宋玉翔，高 豪，劉林夢，任 一，單 軍，張煥朝，王保戰，5?

（1. 南京林業大學林學院，南京 210037；2. 土壤與農業可持續發展國家重點實驗室（中國科學院南京土壤研究所），南京 210008；3. 浙江農林大學亞熱帶森林培育國家重點實驗室，杭州 311300；4. 上海美吉生物醫藥科技有限公司，上海 201318；5. 南京農業大學生命科學學院，南京 210095）

硝化作用被認為是全球氮循環的關鍵步驟，是由微生物驅動通過亞硝酸根，將氨氧化成硝酸根的過程[1]。而氮素作為森林生態系統中主要的限制性營養元素[2]，其硝化作用產生的硝酸鹽更易通過淋溶作用導致氮素流失和地下水污染等，此外硝酸根還可通過反硝化作用產生重要的溫室氣體氧化亞氮（N2O）。自1892年Winogradsky[3]發現氨氧化細菌（Ammonia-oxidizing bacteria，AOB）以來，人們一直認為硝化過程由兩類微生物接力完成，分別為AOB催化的氨氧化過程和亞硝酸鹽氧化細菌（Nitrite-oxidizing bacteria，NOB）催化的亞硝酸氧化過程，其中氨氧化過程被認為是硝化作用的限速步驟[4]。2005年，美國科學家David A. Stahl團隊[5]從海洋中分離出第一株氨氧化古菌（Ammonia- oxidizing archaea，AOA），證明除細菌外，古菌也能夠催化地球上的氨氧化過程。2015年，Daims[6]和van Kessel[7]等均發現了具有氨氧化活性的硝化螺旋細菌（Nitrospira），并證明其可將 4NH+氧化成NO3-，稱為完全硝化微生物（Complete ammonia oxidizer，comammox）。Comammox的發現打破了上百年來人們認為硝化過程必須由兩類微生物接力完成的傳統觀念，使人們再次認識到硝化微生物的物種多樣性和生理代謝機制的復雜。

酸性pH和低氨分子濃度是森林土壤的兩個重要特征。一般認為酸性pH選擇了嚴格嗜酸group 1.1a-associated AOA類群[8-10]，且DNA-SIP（穩定性同位素示蹤）實驗證明該類AOA主導了森林土壤氨氧化過程[11]。但也有一些研究發現group 1.1b AOA可存在于酸性土壤[1，12]，且能主導某些酸性農田土壤氨氧化過程[13]，以及group 1.1b AOA在酸性森林土壤中也廣泛存在，且在轉錄組水平發揮了較高的活性[14-15]。盡管關于AOB適應酸性pH及其相對應的低NH3環境的機制目前多為推測[16-17]，但不少森林土壤中確實存在β-proteobacterium（β-變形菌門）中Nitrosospira（亞硝化螺旋菌）屬AOB，而Nitrosomonas（亞硝化單胞菌）AOB較少檢測到[18-20]。最近日本科學家又從pH 3.0左右的酸性茶園土壤中分離到耐酸γ-proteobacterium（γ-變形菌門）AOB菌株[21]，然而森林土壤中是否存在耐酸 γ-proteobacterium AOB仍不清楚。此外最近備受關注的comammox，其在森林土壤中生理生態的研究也剛剛起步，而已報道的comammoxamoA基因PCR引物易產生非特異性擴增，可能導致其基于實時熒光定量PCR（qPCR）的豐度被嚴重高估[22]。土壤總DNA宏基因組測序由于不受PCR引物偏好性等限制，因此在土壤氨氧化微生物amoA基因相對豐度和群落組成的研究中更具優勢。

馬尾松林（Pinus massoniana）是我國南方廣泛分布的主要人工林之一，具有高經濟價值、耐干旱貧瘠等特點[23]，其土壤為典型的酸性森林土壤。本文以浙江省建德市人工馬尾松林土壤為研究對象，綜合利用qPCR、凝膠電泳半定量和土壤總DNA宏基因組測序等手段，研究酸性森林土壤中AOA、AOB和comammox三種氨氧化微生物的豐度和群落組成等，以期加深對酸性森林生態系統中氨氧化微生物潛在生理生態學特征的認識。

1 材料與方法

1.1 供試材料

土壤樣品于2017年4月采自浙江省建德市人工馬尾松林（119°28′E，29°48′N），該地區屬亞熱帶海洋型季風氣候，海拔約200 m，年均氣溫17.4 ℃，年均降水1 600 mm。土壤類型為凝灰巖發揮的紅壤，選擇3個林齡均為55年的馬尾松林（間隔約2 000 m），分別標記為M1、M2和M3。五點法采集0～20 cm土壤，每個馬尾松林采集3個樣品，樣品分成兩份，一份保存于–20 ℃用于提取土壤總DNA，另一份4 ℃保存用于土壤理化性質測定。

1.2 土壤理化性質測定

土壤含水量采用烘干法。pH采用電位法測量，水土比為5∶1。有機碳含量通過濃硫酸-重鉻酸鉀法測定。全氮和堿解氮的含量分別采用半微量凱氏定氮法和堿解擴散法。有效磷測定采用0.05 mol·L–1HCl和0.025 mol·L–1（1/2H2SO4）浸提—鉬銻抗比色法。速效鉀測定采用乙酸銨浸提—火焰光度法。經測定三個森林土壤pH為3.88～4.16，含水量173.1～203.4 g·kg–1，有機碳14.01～16.62 g·kg–1，全 氮 0.76～1.02 g·kg–1，堿 解 氮 123.24～142.81 mg·kg–1，有效磷6.34～10.49 mg·kg–1，速效鉀56.12～60.15 mg·kg–1。

1.3 土壤DNA提取

土壤微生物總DNA提取采用Fast DNA Spin kit for soil試劑盒（MP Biomedicals，美國）。根據試劑盒操作說明提取土壤微生物總DNA，并使用Nanodrop ND-1000 超微量分光光度計（NanoDrop Technologies，Wilmington，DE，美國）和1%的凝膠電泳檢測DNA濃度和質量，將DNA于–20 ℃保藏。

1.4 實時熒光定量PCR和凝膠電泳半定量分析

利用美國伯樂公司CFX96實時熒光定量PCR儀（Bio-Rad Laboratories，Hercules，美國）對人工馬尾松林土壤中AOA、AOB和comammox的amoA基因進行定量分析。反應體系為20 μL：10 μL SYBR Premix Ex Taq（TaKaRa，日本），上下游引物0.5 μmol·L–1，2 μL（約1～10 ng）模板DNA，補充6 μL雙蒸水。qPCR引物和反應條件見表1[22，24-26]。本實驗每個樣品設置三個重復，每個重復包含三個測試平行。將含有AOA、AOB和comammoxamoA基因的單克隆質粒DNA梯度稀釋成101～109的標準濃度作為qPCR標準曲線。通過瓊脂糖凝膠電泳檢查熒光定量PCR擴增產物的特異性。

表 1 定量PCR擴增引物及反應條件 Table 1 Quantitative PCR amplification primers and their reaction conditions

本實驗和已報道研究[27]均發現目前已知的comammox的amoA基因PCR引物易導致非特異性擴增，致使amoA基因豐度高估。因此利用凝膠電泳半定量法對qPCR結果進行數值初步矯正。即通過1.2%凝膠電泳和Bio-Rad公司的Quantity one v4.6.6軟件檢測待測樣品qPCR產生的目標條帶和非目標的相對灰度值，從而計算目標條帶占泳道中總灰度值的百分比，進而用該百分比乘以qPCR檢測結果，即得到樣品中comammox的amoA基因實際拷貝數。同時我們將進一步用宏基因組中AOA、AOB和comammox三種amoA基因相對豐度，以及AOA和AOB的amoA基因的qPCR結果，進一步進行計算驗證（見1.5）。

1.5 宏基因組測序和功能基因拼接及其相對豐度計算

委托上海美吉生物對土壤總DNA進行150 bp雙端測序（Illumina HiSeq 2 000）。每個樣品分別產生約48Gb的原始數據，然后通過SeqPrep（https：//github.com /jstjohn/SeqPrep）去除末端測序接頭（adaptor），然后利用Sickle（https：//github.com/ najoshi/sickle）以默認參數對原始序列進行了質量控制和過濾，三個樣品分別得到約42G bp高質量數據。使用metaSPAdes（version 3.10）宏基因組拼接軟件對測序數據進行de novo從頭組裝，k-mer參數設置為21～53。使用Glimmer 3.0進行開放閱讀框（ORFs）預測。以NCycDB數據庫[28]中AOA和AOB的amoA基因以及Pjevac[22]、Daims[6]和Li[29]等文章中的comammox的amoA基因為參考序列，利用BLASTX與宏基因組拼接片段進行比對，提取Evalue ≤10–5且同源性高于70%的contigs作為amoA的候選序列，并進一步通過NCBI的nr數據庫注釋以及構建amoA進化樹等手段進行驗證去除錯誤注釋的amoA基因序列。

分別采用兩種方法計算宏基因組中AOA、AOB和comammox的amoA基因相對豐度。第一種方法，基于Mapped到amoA基因的HiSeq reads數統計。第二種方法，統計宏基因組中參與三種amoAcontigs拼接所用的HiSeq reads數，reads數高低反映了該類基因的相對豐度高低。即根據amoA基因參考序列，通過BLASTX與原始reads進行比對，統計其中Evalue ≤10–5，同源性≥70%，blast hit區域大于25個氨基酸長度的reads數目，通過比較三種amoA基因mapping的reads總數，確定三者間的相對豐度。

1.6 系統發育進化樹構建

首先通過MOTHUR軟件[30]按照95%序列同源性提取AOA、AOB和comammox的amoA基因的OTUs（operational taxonomic unit）代表性序列，然后利用MEGA 7.0軟件的最大似然法（Maximum likelihood）構建三種氨氧化微生物amoA基因的進化樹。

2 結 果

2.1 AOA、AOB和comammox amoA 基因豐度

qPCR結果顯示酸性馬尾松林土壤中AOA和AOB的amoA基因豐度分別為2.61×106copies·g–1和1.45×106copies·g–1（圖1a）。基于引物Ntsp-amoA- 162F/359R和引物comaA/B-244F/659R的qPCR結果顯示comammoxamoA基因拷貝數分別為4.14× 106copies·g–1和2.64×106copies·g–1。然而，凝膠電泳圖顯示兩組引物的qPCR結果均有明顯的非特異性擴增，導致comammox的amoA基因qPCR結果的高估（圖1a）。通過凝膠電泳和Quantity one v4.6.6對comammox的amoA基因目標條帶進行半定量估算后，得到兩種引物擴增的comammoxamoA基因總拷貝數分別為 1.38×106copies·g–1和 1.47×106copies·g–1（圖1a）。因此，氨氧化微生物功能基因amoA的相對豐度，comammox：AOA為0.53～0.56；comammox：AOB為0.95～1.01（圖1a）。

基于宏基因組Mapped到amoA基因的HiSeq reads數目和參與amoAcontigs拼接的HiSeq reads數目分析發現，三類氨氧化微生物amoA基因相對豐度為：comammox：AOA約為0.55和0.29，comammox：AOB約為0.98和0.59（圖1b）。

2.2 AOA、AOB和comammox的系統發育和群落組成

目前AOA主要包括Nitrosopumilus（group 1.1a）、Nitrosotalea（group 1.1a-associated）、Nitrososphaera（group 1.1b）和Nitrosocaldus（ThAOA）[31]。宏基因組分析發現，馬尾松林AOA其主要類群為group 1.1b，約占AOA總豐度的85.66%～88.07%，其次為group 1.1a-associated，占11.93%～14.34%，未檢測到group 1.1a和Nitrosocaldus類群（圖2）。

AOB 主要分為 β-AOB（Nitrosospira和Nitrosomonas）和γ-AOB（Nitrosococcus，亞硝化球菌屬）[32]。其中Nitrosospira包含cluster 1～4、cluster 0、Nitrosospirasp. Nsp65和Nitrosospirasp. Nsp57/58等分支[33-34]。Nitrosomonas包含N. europaea/N. mobilis、N. communis、N. marina、N. oligotropha、N. cryotolerans、Nitrosomonassp. Nm143和cluster 5等分支[35]。宏基因組分析顯示，Nitrosospira是AOB的主要類群，占AOBamoA基因總豐度的51.71%～63.96%，主要隸屬于cluster 3、0和Nitrosospira-like（圖 3），而Nitrosomonas占 36.04%～48.29%，主要隸屬于Nitrosomonas communis和Nitrosomonas oligotropha（圖 3）。

目前報道的comammox均來源于Nitrospiralineage Ⅱ，且主要分為clade A和clade B兩支[6-7]，近期有研究將comammox clade A進一步分為兩個單系群，clade A.1和clade A.2[36]。土壤宏基因組分析發現，酸性森林土壤中comammox的主要類群為clade B，占comammoxamoA基因總豐度的63.39%～71.39%，此外clade A.1占總序列數的28.61%～36.61%，未檢測到clade A.2序列（圖4）。

3 討 論

馬尾松林土壤屬于典型酸性森林土壤。為了得到AOA、AOB和comammoxamoA基因準確豐度，本研究同時使用了qPCR、凝膠電泳半定量和宏基因組三種方法進行數據比較分析。qPCR和宏基因組均表明AOAamoA基因豐度約為AOBamoA基因的兩 倍（圖1）。而comammoxamoA基因的兩套引物存在明顯的非特異擴增，該結果與Beach和Noguera[27]研究結果一致（圖1a）。通過凝膠電泳半定量對qPCR結果進行矯正以及宏基因組測序分析都發現，comammoxamoA基因拷貝數和AOB的amoA基因處于同一個水平，均為AOAamoA基因拷貝數的50%左右（圖1）。然而，目前comammox的qPCR 定量仍處于起步階段。而凝膠電泳半定量也存在一定缺陷，其反映是qPCR結束時目標電泳條帶和非目標電泳條帶的灰度值之比，無法真正反映qPCR反應過程中特異擴增和非特異擴增釋放的熒光信號強度之比。此外，凝膠電泳半定量目標條帶位置，仍可能包含非特異擴增的序列。而宏基因組技術目前成本較高，數據分析也較繁瑣。因此，仍需進一步發展comammoxamoA基因qPCR更有效的引物，或選擇comammox基因組中其他marker基因設計qPCR引物。

pH和NH3濃度是影響酸性土壤氨氧化微生物群落最重要的兩個環境因子。自首株嚴格嗜酸AOA菌株分離以來[8]，大量基于DNA-SIP實驗證據表明group 1.1a-assocaited AOA主導了很多酸性土壤的氨氧化過程[9，11，37]。然而我們發現group 1.1b AOA也可以在某些酸性土壤中發揮著主導作用[13]。James I. Prosser團隊[38-39]通過全球、區域和田塊三個尺度研究發現group 1.1b AOA也存在于酸性土壤中，且亞枝cluster 11被認為嗜酸生長。本研究發現馬尾松林土壤中group 1.1b AOA占據數量上的優勢，而group 1.1a associated僅占12%（圖2）。此外，酸性土壤中NH3而被大量離子化成NH4+，而目前認為氨氧化微生物的底物為氨分子而非銨離子[40]，因此推測酸性土壤中group 1.1b AOA能夠適應較低氨環境。

酸性土壤中AOB的貢獻長期以來一直存在較大爭議。在AOA發現前，一直認為AOB驅動了酸性土壤的氨氧化過程。但后來發現分離于酸性土壤的Nitrosospira和Nitrosomonas純菌株都不能在實驗室內的酸性pH條件下生長[16]。Prosser等[16-17]推測可能有三種原因：（1）土壤中存在微域中性或堿性pH環境；（2）AOB在土壤中形成生物膜或團聚體抵抗酸性pH；（3）很多AOB菌株具有尿素水解酶活性，能夠通過利用尿素抵抗（酸性土壤中）低氨脅迫。我們認為第三種解釋可能混淆了酸性pH和低氨濃度兩種環境脅迫，通過利用尿素確實可以應對低氨寡營養脅迫，但該AOB菌株能夠在酸性pH條件下生存，還必須同時能夠抵抗低pH脅迫。推測James I. Prosser等當時實驗選擇的Nitrosospira純菌株本身能夠耐受酸性pH，它只需要在此基礎上通過利用尿素應對低氨脅迫即可在酸性pH培養基中生長。這點可以通過Prosser等[16-17]兩篇文章中推測出來，五株來源于酸性土壤的AOB純菌株均可以利用尿素，但只有一株Nitrosospira菌株可以通過利用尿素在酸性pH培養基中生長。因此，利用尿素對以上來源于酸性土壤的五株AOB純菌株來說是適應酸性pH的必要但非充分條件。2017年Hayatsu等[21]分離得到在不添加尿素條件下，能夠在酸性pH中生長的AOB純菌株γ-proteobacterium TAO100。然而目前β-proteobacterium AOB適應酸性土壤的分子生物學機制仍不清楚。。

作為最新發現的comammox，其在地球各個生境中的分布規律備受關注。目前發現comammox在土壤、沉積物、灘涂、植物葉表、地下水、飲用水以及廢水處理系統均有分布[6-7，22，36]。本研究發現酸性森林土壤也存在豐度與AOB（每個AOB細胞平均含有1.5個拷貝的amoA基因）相當的comammox 類群（圖1），其中clade B占comammox總類群的64%，而clade A約占36%。該結果和已報到森林土壤中comammox群落組成基本一致[22，36]。不同的是，Xia等[36]發現森林土壤中clade A主要隸屬于clade A.2分支，而本研究發現馬尾松林土壤中clade A主要隸屬于clade A.1（圖4）。以上數據暗示comammox無論在clade B還是clade A類群中都存在能夠適應酸性土壤環境的物種。Michael Wagner課題組[41]繼2015年發現第一株comammox純菌株后，也明確了已知的comammox純菌株對底物氨具有很高的親和力，其親和力高于所有已知AOB菌株和絕大部分土壤來源的AOA菌株，因此這可以從一定程度解釋comammox如何適應酸性土壤中的低NH3脅迫。然而，comammox如何適應酸性pH脅迫仍不清楚。

4 結 論

本文利用qPCR、凝膠電泳半定量和宏基因組測序等技術研究酸性森林土壤氨氧化微生物的群落結構。qPCR顯示AOA和AOB的amoA基因豐度為2.61×106copies·g–1和 1.45×106copies·g–1。而comammoxamoA基因的擴增引物產生明顯的非特異性擴增，導致qPCR結果的高估。半定量和宏基因組分析發現，氨氧化微生物amoA基因相對豐度為：comammox：AOA=0.55，comammox：AOB=0.98，推算comammoxamoA基因真實豐度為（1.38～1.47）×106copies·g–1。此外，AOA主要類群為group 1.1b占其總豐度的88.07%，而經典嗜酸類群group 1.1a-associated僅占11.93%。AOB主要類群為Nitrosospira（63.96%），而Nitrosomonas為36.04%。comammox主要類群為clade B（63.39%），而clade A僅占36.61%且均隸屬于clade A.1亞枝。

致 謝感謝浙江農林大學孫凱和張云晴等同學在樣品采集和理化性質測試中給予的幫助！