三峽消落帶土壤古菌群落組成與代謝功能預測

葉 飛，吳曉冰，姜 毅，吳勝軍，王 雨*

（1.中國科學院 重慶綠色智能技術研究院，重慶 400714；2.中國科學院水庫水環(huán)境重點實驗室，重慶 400714）

長期以來，古菌被認為僅存在于強堿、強酸或者高鹽等極端環(huán)境中。當與極端環(huán)境古菌高度相似的“泉古菌”在海洋浮游生物中被發(fā)現(xiàn)后[1-2]，一系列后續(xù)研究無不證明古菌是地球生物圈中一類豐富且多樣的微生物，廣泛地分布于地球的各類環(huán)境中，改變了人們對古菌的傳統(tǒng)認識[3-7]。數(shù)量龐大的古菌在生物地球化學過程中的重要作用遠超預期。

在碳循環(huán)中，產(chǎn)甲烷古菌產(chǎn)生的甲烷占全球甲烷排放量的70%[8]。除此之外，古菌還積極參與硫循環(huán)，能將二硫化碳（CS2）快速轉(zhuǎn)化為硫化氫（H2S）和二氧化碳（CO2）[9]。進行化能自養(yǎng)的氨氧化古菌能將氨（NH3）氧化為亞硝酸鹽（NO2-）[10]。氨氧化古菌的發(fā)現(xiàn)徹底改變了氨氧化過程由細菌驅(qū)動的傳統(tǒng)認識，它們在自然界中的龐大數(shù)量意味著其在自然界物質(zhì)轉(zhuǎn)化尤其是全球碳氮循環(huán)方面扮演著重要角色[4，7]。

三峽水庫實行夏落冬漲的調(diào)度運行方案，夏季汛期（6—9月）保持145 m 最低水位，汛期末（10月）開始蓄水，冬季（11月）達到175 m最高水位，之后水位逐漸消落，次年5月底降至最低水位，水位變幅大且同自然洪枯規(guī)律相反。季節(jié)性水位漲落使庫區(qū)淹沒土地周期性露出水面，形成了面積達348.9 km2的消落帶。三峽庫區(qū)消落帶由于其所處的特殊空間位置，同時受到人類活動和水位波動的干擾[11-12]，導致消落帶土壤環(huán)境比較復雜，進而影響消落帶微生物群落的分布。然而目前對消落帶古菌群落組成的空間分布規(guī)律還缺乏深刻的認識。此外，了解古菌群落的潛在功能對于更好地認識其在生物地球化學循環(huán)中的角色是至關重要的，有助于理解長江流域氮磷營養(yǎng)鹽循環(huán)過程，對流域水資源保護有重要意義。近年來越來越多的研究采用宏基因組測序的方法來揭示環(huán)境中微生物群落功能[13-14]。作為一種比較經(jīng)濟的生物信息學方法，PICRUSt（phylogenetic investigation of communities by reconstruction of unobserved states）可以通過16S rRNA基因數(shù)據(jù)推斷微生物群落的功能屬性，具有較好的預測結(jié)果，已廣泛應用于人體或哺乳動物[15]、土壤[16]、海水和沉積物[17]等生境研究，是了解古菌群落代謝功能的有效途徑。本研究以不同高程消落帶土壤為研究對象，分析其古菌群落組成沿高程梯度的分布規(guī)律，并采用基于16S rRNA 高通量測序的PICRUSt 功能預測分析消落帶土壤古菌群落的代謝功能特征。

1 材料與方法

1.1 樣品的采集與保存

本研究選取的白家溪消落帶位于長江一級支流澎溪河流域，處于三峽庫區(qū)腹地的重慶市開州區(qū)境內(nèi)（圖1）。該區(qū)域?qū)儆谥衼啛釒駶櫦撅L氣候，年平均氣溫18.5 ℃，四季分明，年均降雨量大約1 200 mm，雨季長、雨日多、降水豐富。土壤類型以發(fā)育自紫色砂巖的紫色土為主。

圖1 研究區(qū)和采樣點示意圖Fig.1 Map of the study region and sampling sites

以海拔150 m 為最低點，沿消落帶150～175 m海拔高程選擇一個垂直樣帶，每隔5 m 設置采樣點，共6個采樣點。在消落帶水位上升初期的2013年10月10日和2014年9月23日，以及水位下降末期的2014年4月1日和2015年5月6日共進行兩年4次土壤樣品的采集。水位線以上的高程樣點，用不銹鋼土鉆采集；水位線以下用彼得森采樣器采集。每個樣點采集3 個重復樣，充分混勻成一個混合樣，共形成24 個混合樣。樣品采集后放置于無菌密封袋，置于冰盒中保存運輸至實驗室。其中一份立即進行理化指標分析，另一份于-20 ℃儲存，用于DNA提取和后續(xù)分子生物學實驗。

1.2 土壤理化性質(zhì)分析

土壤pH采用pH計（FE20，Mettler Toledo，美國）測定，水土比為1∶5。含水率（Moisture）用烘干質(zhì)量測定法。土壤NH4+和NO3-的含量用2 mol/L KCl 浸提后，采用連續(xù)流動注射分析儀（FIA Star 5000，F(xiàn)OSS Tecator，瑞典）測定。土壤有機質(zhì)（origanic matter，OM）采用燒失量法測定?？偺迹╰otal carbon，TC）、總氮（total nitrogen，TN）和總硫（total sulfur，TS） 用元素分析儀 （Vario EL cube，Elementar，德國）測定[18]。

1.3 DNA的提取和Illumina高通量測序

采用Power Soil?DNA Isolation Kit（Mobio，美國）試劑盒提取土壤樣品中的總DNA。用超微量紫外分光光度計NanoVue Plus Spectrophotometer（GE Healthcare，英國）確定DNA濃度。用濃度為1%（w/v）的瓊脂糖凝膠電泳檢驗DNA的質(zhì)量。提取的DNA溶液保存于-20 ℃冰箱中待用。

針對古菌16S rRNA 基因的V3-V5 高變區(qū)，使用帶標簽序列（barcode）的引物Arch344F（5′-ACGGGGYGCAGCAGGCGCGA-3′） 和 Arch915R（5′-GTGCTCCCCCGCCAATTCCT-3′）[19]進行聚合酶鏈反應（polymerase chain reaction，PCR 擴增。擴增采用TransGen AP221-02（TransGen，中國）的 20 μL 體系，具體為：5 μL 5 × FastPfu buffer，2 μL dNTPs（脫氧核糖核酸，2.5 mmol/L），0.3 μL FastPfu DNA 聚合酶（2.5 U），正、反向引物（5 μmol/L）各0.8 μL，0.2 μL BSA（牛血清白蛋白，20 μg/μL）（Takara，中國），DNA模板（10 ng），雙蒸水補足至20 μL。PCR 擴增程序為：95 ℃預變性3 min，然后95 ℃變性30 s，55 ℃退火30 s，72 ℃延伸45 s，共32 個循環(huán)，最后再72 ℃延伸10 min。擴增所用的PCR 儀為ABI GeneAmp?9700。PCR 產(chǎn)物使用AxyPrepDNA 凝膠回收試劑盒（AXYGEN，美國）純化后用2%瓊脂糖凝膠電泳進行檢測。參照電泳初步定量結(jié)果，將PCR 產(chǎn)物用QuantiFluor?-ST 藍色熒光定量系統(tǒng)（Promega，美國）進行檢測定量，之后按照每個樣本的測序量要求，進行相應比例的混合后雙端（pair-end，PE）測序（2 × 300）。本研究的Illumina 高通量測序委托上海美吉生物醫(yī)藥科技有限公司①http://www.Majorbio.com在Illumina Miseq PE300 平臺上完成。在一個淹水周期內(nèi)，160 m 樣點的淹水和出露時間與165 m 樣點非常接近，因此，在進行高通量測序及后續(xù)相關分析過程中將160 m高程采集的樣品排除在外。

1.4 高通量測序序列處理

MiSeq測序得到的是雙端序列數(shù)據(jù)，首先根據(jù)雙末端片段（PE reads）之間的重疊（overlap）關系，將成對的片段（reads）拼接（merge）成一條序列，同時對reads 的質(zhì)量和拼接的效果進行質(zhì)控過濾，根據(jù)序列首尾兩端的barcode 和引物序列區(qū)分樣品得到有效序列，并校正序列方向，進行數(shù)據(jù)優(yōu)化。優(yōu)化后的高質(zhì)量序列使用Usearch 軟件（vsesion 7.0，http://drive5.com/uparse/）在97%的相似度下劃分獨立操作單元（Operational Taxonomic Units，OTUs）。為了得到每個OTU對應的物種分類信息，采用核糖體數(shù)據(jù)庫項目分類器（RDP classifier）貝葉斯算法對OTU 代表序列進行分類學分析，分析比對使用SILVA rRNA 數(shù)據(jù)庫（Release 128，http://www.arb-silva.de）。為了便于比較各樣品間的差異性，根據(jù)測序樣品中的最小序列數(shù)將所有樣品的序列稀釋到同一測序深度（每個樣品21 068條序列）。

1.5 PICRUSt功能預測與統(tǒng)計分析

基于古菌群落16S rRNA 基因代謝功能的預測通過PICRUSt v1.0.0 軟件實現(xiàn)[16]。參照Greengenes數(shù)據(jù)庫（13_5_release）將用于PICRUSt 預測的序列劃分為OTUs（97%相似性）。在比對過程中剔除與參照序列不符合的序列。根據(jù)OTU 對應基因組中16S rRNA 基因拷貝信息，將每個OTU 對應序列數(shù)除以其16S拷貝數(shù)進行標準化。最后，將標準化的數(shù)據(jù)乘以其對應的基因組中基因含量，從而實現(xiàn)代謝功能預測，產(chǎn)生KEEG（Kyoto Encyclopedia of Genes and Genomes） 同源 序 列 （KEEG orthologs，KOs）表。預測結(jié)果的準確度可以通過參數(shù)NSTI（Nearest Sequenced Taxon Index）值的大小來評估。NSTI 表征某樣品中所有微生物OTU 與其親緣關系最近的已測序基因組間系統(tǒng)進化距離的評價值。因此，該值越小表示預測結(jié)果越準確。

不同高程土壤理化性質(zhì)差異利用IBM SPSS Statistics 20.0 for Windows 的單因素方差分析（One-way ANOVA）進行比較。使用R 語言（version 3.0.1）的“vegan”包計算基于布雷-柯蒂斯距離（Bray-Curtis Distance）的不同高程土壤古菌群落的差異性（Analysis of Similarity，ANOSIM；置換檢驗等于999），并進行古菌群落組成變化與土壤性質(zhì)關系的Mantel檢驗。主坐標軸分析（Princi-pal Coordinates Analysis，PCoA）和冗余分析（Redundancy Analysis，RDA） 通過 CANOCO 5 軟件實現(xiàn)。圖形繪制使用ORIGIN 9.0軟件完成。

2 結(jié)果

2.1 土壤性質(zhì)沿高程的分布

消落帶土壤NH4+和NO3-濃度、pH 值以及含水率均表現(xiàn)出沿采樣高程的梯度分布規(guī)律（圖2）。其中，NH4+，pH 值和含水率表現(xiàn)出隨著采樣高程降低而增加的趨勢。但150 m高程樣點的NH4+濃度和含水率要明顯低于155 m 高程（P＜0.05），這種情況很可能是因為150 m 高程樣點位于河床附近，采集的土壤樣品中含有大量的砂石。此外，消落帶土壤中的NO3-濃度總體上隨著采樣高程的升高而增加。而土壤中OM，TC，TN 和TS 的變化都與采樣高程沒有關系。

圖2 消落帶4次采樣不同高程土壤性質(zhì)的均值Fig.2 The mean values of the soil properties at different elevations in four sampling surveys.

2.2 古菌群落組成變化及驅(qū)動因子

通過基于OTU 水平的ANOSIM 分析，比較不同高程樣點間古菌群落組成的差異（表1），可以發(fā)現(xiàn)各高程土壤古菌群落組成與相鄰高程沒有顯著性差異，但是與相隔較遠高程的土壤差異顯著。如175 m高程樣點的古菌群落組成與170 m沒有顯著差異（P＞0.05），而與165 m，155 m 和150 m 高程的土壤古菌群落都有顯著差異（P＜0.05）。170 m高程樣點與165 m 高程樣點的古菌群落無顯著差異 （P＞0.05），而與155 m 和150 m 高程呈現(xiàn)顯著差異（P＜0.05）。說明消落帶土壤古菌群落組成沿高程梯度呈現(xiàn)出明顯的空間分異特征。Mantel檢驗結(jié)果表明，消落帶古菌群落組成的變化主要與土壤NH4+，TN和含水率顯著正相關（P＜0.05，表2）。

表1 基于ANOSIM分析的不同高程間土壤古菌群落組成比較Tab.1 Comparison of soil archaeal community compositions between different sampling elevations based on ANOSIM analysis

表2 基于Mantel檢驗的土壤性質(zhì)與古菌群落組成變化的相關關系Tab.2 The correlations of soil properties with the variation of archaeal community composition based on Mantel test

2.3 代謝功能預測

本研究中每個樣品的NSTI值范圍為0.126～0.296，平均為0.238。采用PICRUSt 預測的功能基因組成的變化沿高程可以分為2 個主要的組，分別是170～175 m 和150～155 m 樣點組。組內(nèi)各自對應的兩個高程樣品的功能組成比較類似。然而，165 m高程樣點分布比較零散，此高程功能基因組成表現(xiàn)出一定的隨機性[圖3（a）]。此外，本研究將具有代謝功能的基因組從KOs 中篩選出來，對土壤中代謝功能組成進行PCoA 分析，發(fā)現(xiàn)與全功能組成分布類似，具有代謝功能的基因組成同樣沿高程可以分為170～175 m和150～155 m樣點2個組[圖3（b）]。

圖3 采用PICRUSt 預測的土壤樣品（a）KEGG orthologies（KOs）和（b）在KEGG功能層次結(jié)構(gòu)中被歸為“新陳代謝”的KOs的主坐標軸分析（PCoA）Fig.3 Principal coordinate analysis (PCoA) of soils using PICRUSt predicted (a) KEGG orthologies (KOs) and (b) KOs only grouped as “Metabolism” in KEGG functional hierarchies

為了解消落帶土壤中代謝功能變化的驅(qū)動因子，本研究采用RDA 方法分析代謝功能組成與土壤性質(zhì)的相關關系。RDA1 和RDA2 兩軸分別解釋了代謝功能變化的32.6%和3.6%（圖4）。土壤樣品預測的代謝功能變化與土壤的含水率、OM和TC顯著相關（圖4）。但是沒有任何一個測得的土壤性質(zhì)與某一高程樣點預測的新陳代謝功能的分布模式有明顯的相關關系。

圖4 預測的代謝功能群落與土壤性質(zhì)的冗余分析（RDA）Fig.4 Redundancy analysis (RDA) of predicted metabolism functional communites and soil properties

3 討論

三峽消落帶處于濕地生態(tài)系統(tǒng)次生演替初期階段，大量的本地物種消失、多樣性降低[11]。微生物群落作為消落帶生態(tài)系統(tǒng)的關鍵組成部分，處于適應變化環(huán)境的過程中，最終會形成一個更穩(wěn)定的群落來應對這種干擾[20-21]。因此微生物群落對于指示消落帶生態(tài)系統(tǒng)功能的恢復[22]以及生物地球化學循環(huán)[23]等方面具有重要作用。

古菌群落組成的PCoA 分析表明，消落帶土壤中古菌群落具有明顯的空間異質(zhì)性。這可能跟淹水導致不同高程土壤的氧環(huán)境存在較大差異有關。絕大多數(shù)的廣古菌門（如產(chǎn)甲烷菌）與土壤含水量顯著正相關[24]。此外，穩(wěn)定的氧化還原條件更有利于廣古菌門微生物的生長[25]。在三峽庫區(qū)消落帶，較低高程土壤的淹水時間較長，形成的相對穩(wěn)定的厭氧環(huán)境更適合廣古菌。而奇古菌門主要由好氧氨氧化古菌組成，短期的淹水導致較高高程的消落帶土壤形成較好的微氧生境，為好氧氨氧化古菌提供了良好的生境[26]。除此之外，好氧氨氧化古菌的分布還與氨氮濃度相關[27]。Mantel檢驗同樣證實消落帶古菌群落的分布與土壤氨氮濃度和含水率顯著正相關。

本研究中各樣品的NSTI均值為0.238，預測精度要低于已有研究，這可能是消落帶土壤環(huán)境包含很多未知的多樣性[16]。消落帶土壤古菌功能基因組成變化表現(xiàn)出與古菌群落組成變化類似的模式，即沿高程的空間分布。175 m 和170 m 高程土壤古菌功能基因組成類似，而155 m和150 m高程的功能基因組成類似，但較高高程和較低高程之間存在較明顯的差異。然而，不同高程功能基因組成的整體變化要弱于古菌群落組成的變化，古菌群落這種預測功能結(jié)構(gòu)可能與其功能性冗余有關[28]。總的來說，與古菌群落結(jié)構(gòu)組成的驅(qū)動因子一樣，含水率同樣是影響預測的代謝功能組成的主要環(huán)境因子。古菌中的產(chǎn)甲烷古菌和好氧氨氧化古菌是碳氮循環(huán)中的重要組成部分[8，10]，此外某些古菌還是硫循環(huán)的主要參與者[9]。其中，產(chǎn)甲烷古菌和氨氧化古菌的分布與土壤含水率密切相關，這在一定程度上反映了古菌代謝功能沿高程的空間分異，較高高程（175～170 m 高程組）主要行使氨氧化功能，而較低高程（155～150 m 高程組）的古菌主要參與碳循環(huán)。因此，長時間的淹水可能會導致三峽消落帶土壤氮循環(huán)和碳循環(huán)在空間上的分異。