晁韶良，郇 環(huán)，周向同，劉 媚

1. 中國環(huán)境科學(xué)研究院，北京 100012

2. 生態(tài)環(huán)境部土壤與農(nóng)業(yè)農(nóng)村生態(tài)環(huán)境監(jiān)管技術(shù)中心，北京 100012

3. 江蘇大學(xué)環(huán)境與安全工程學(xué)院，江蘇 鎮(zhèn)江 212013

4. 中國地質(zhì)大學(xué)（武漢）環(huán)境學(xué)院，湖北 武漢 430078

20世紀(jì)90年代以來，我國農(nóng)田施用氮肥量不斷增加，土壤中殘留的硝態(tài)氮造成農(nóng)業(yè)面源氮素污染.土壤中硝態(tài)氮受灌溉和降水影響，易淋溶到深層包氣帶，進(jìn)而滲透到地下水危害地下水環(huán)境[1-4]. 包氣帶是復(fù)雜的非飽和帶，包含氣、液、固三相和多種微生物.包氣帶對污染物具有吸附作用，同時微生物的生命活動也會影響包氣帶中污染物的遷移[5]. 反硝化作用對減緩NO3

—向地下水遷移有重要影響，其強(qiáng)度可以指示包氣帶的NO3—防護(hù)效果[6-7]. 了解包氣帶反硝化微生物的垂直多樣性和豐度對于提高作物產(chǎn)量、應(yīng)對環(huán)境污染和調(diào)節(jié)生物地球化學(xué)氮循環(huán)至關(guān)重要.

反硝化作用是指反硝化菌將硝酸鹽(NO3—)依次還原為亞硝酸鹽(NO2—)、一氧化氮(NO)、一氧化二氮(N2O)和氮?dú)?N2)的過程. 亞硝酸鹽還原的過程是反硝化作用中重要的限速過程，nirK和nirS基因分別編碼催化該過程的酶，主要包括Cu型亞硝酸鹽還原酶(Cu-Nir)和細(xì)胞色素還原酶(cd1-Nir)，一氧化二氮還原酶(Nos)是催化反硝化作用的最后一步[8]. 通過研究nirS、nirK和nosZ等標(biāo)記基因[9-13]，可以反映出反硝化菌的生態(tài)行為和對環(huán)境的響應(yīng)，可以進(jìn)一步指示相關(guān)環(huán)境的反硝化強(qiáng)度及其影響因素.

包氣帶是一個具有多種環(huán)境梯度和高度異質(zhì)性的生境，微生物豐度、多樣性、群落組成和功能受到深度的控制，并在包氣帶形成垂直分層[14-16]. 目前，關(guān)于反硝化菌多樣性和群落結(jié)構(gòu)研究主要集中在表層包氣帶(＜2 m)，深層包氣帶的研究相對匱乏，且已有的針對深層包氣帶的研究主要集中在微生物豐度方面. 例如：Tang等[17]收集了森林0～10、10～20、20～40、40～60和60-80 cm的土壤樣品，發(fā)現(xiàn)氮循環(huán)基因nirS和nirK的豐度隨深度的增加而減少. Liu等[18]研究表明，nosZ相比nirS和nirK更容易出現(xiàn)在較深(40～80 cm)的土壤中. 相比表層土壤，深層土壤養(yǎng)分有限，溫度和水分相對恒定，微生物群落可能也明顯不同. He等[16]對40個土壤樣品進(jìn)行了分析，覆蓋了0～600 cm的深度范圍， 結(jié)果表明，包氣帶深度是微生物豐度分布的主要影響因素. 反硝化菌的多樣性和群落結(jié)構(gòu)與農(nóng)田演替、土壤氧化還原電位、土壤類型、氣候變化、pH、有機(jī)肥料組成等存在相關(guān)關(guān)系[9,11-13,19].Chen等[20]通過RDA分析發(fā)現(xiàn)TOC含量(P＜0.01)、土壤AP含量(P＜0.01)、Cu含量(P＜0.01)和Zn含量(P＜0.01)均與nirK反硝化群落組成呈顯著相關(guān). 瀏陽市紅壤水稻田間相關(guān)研究發(fā)現(xiàn)，土壤pH對nirK群落組成的解釋比例(32.7%，P＜0.05)和對nirS群落組成的解釋比例(31.0%，P＜0.05)最大[21]. Liang等[9]采集了熱帶、亞熱帶、暖溫帶和中溫帶氣候區(qū)域的40種水稻土，結(jié)構(gòu)方程模型結(jié)果表明，氣候因子和pH是影響反硝化酶活性差異的主要原因. Tao等[11]研究了不同施肥情況對農(nóng)田N2O排放和反硝化群落結(jié)構(gòu)的影響，結(jié)果表明，有機(jī)肥增加了反硝化基因拷貝數(shù).由此可知，土壤反硝化菌群落結(jié)構(gòu)受到不同環(huán)境因子的影響.

潮白河沖洪積扇區(qū)是北京第二大沖洪積扇，是北京市最大的地下水供水水源地[22]. 2003年北京密云、懷柔、順義地區(qū)地下水水源地的硝酸鹽含量持續(xù)上升[7,23-24]. 目前潮白河沖洪積扇區(qū)域的氮素研究多集中在硝態(tài)氮的時空變化[25-27]和反硝化強(qiáng)度垂向空間分布規(guī)律[28]方面，截至目前鮮見學(xué)者開展潮白河反硝化菌群落結(jié)構(gòu)及其與環(huán)境因子相關(guān)性研究.為此，本文選取潮白河區(qū)域2處典型包氣帶剖面(S6和S8)，采集深度10.0 m以內(nèi)的包氣帶樣品作為研究對象，首先采用高通量測序分析反硝化菌(nirS、nirK和nosZ)多樣性和群落結(jié)構(gòu)，再利用斯皮爾曼(Spearman)非參數(shù)分析方法研究影響反硝化細(xì)菌多樣性的環(huán)境因子，并使用冗余分析方法、隨機(jī)森林分析方法和曼特爾檢驗(yàn)(Mantel Test)相關(guān)性分析方法，研究影響反硝化菌群落結(jié)構(gòu)的環(huán)境因子. 本研究旨在深入了解包氣帶反硝化菌群落結(jié)構(gòu)，探究不同環(huán)境因子對反硝化菌群落結(jié)構(gòu)的影響程度，并為包氣帶對地下水污染防護(hù)能力的研究提供數(shù)據(jù)支持和理論依據(jù).

1 材料與方法

1.1 研究區(qū)概況

潮白河沖洪積扇區(qū)地下水資源豐富，是北京東北區(qū)域重要的水源地之一. 該地區(qū)水位埋深情況：密云地區(qū)(上游)平均埋深為15～30 m，懷柔和順義地區(qū)(中下游)為5～15 m，通州地區(qū)(下游)為5～10 m. 地下水總體從南到北流動，大氣降水是最主要的補(bǔ)給來源[21]. 潮白河洪積扇區(qū)上游(扇頂)地層主要由較厚的砂卵礫石組成，含水層比較單一. 中游(扇中)地層包含了礫砂、粗砂和黏土、亞黏土. 下游(扇緣)地層含有更細(xì)的顆粒細(xì)砂、粉砂與黏土、亞黏土，是多層含水層結(jié)構(gòu).

本研究選取潮白河沖洪積扇中下部典型包氣帶作為研究區(qū)域，并于2016年10月27日—2016年11月4日使用PowerProbe鉆機(jī)(9500-VTR，美國AMS公司)在S6和S8剖面(見圖1)采樣. S6剖面共11個樣品，依次記為S61、S62、S63、S64、S65、S67、S68、S69、S610、S611，分別對應(yīng)埋藏深度0.2、0.4、0.6、0.8、1.0、1.2、2.0、3.0、4.0、5.0、10.0 m. S8剖面缺少1.2 m和10.0 m深度的樣品，編號為S81～S89.采集的土壤樣品應(yīng)避免人員接觸污染，采集后密封，在4 ℃環(huán)境下保存，標(biāo)注好樣品的采樣時間、地點(diǎn)和采樣深度.

圖1 潮白河沖洪積扇包氣帶采樣點(diǎn)分布示意Fig.1 Distribution of sampling points in the aerated zone of the Chaobai River alluvial fan

1.2 土壤理化性質(zhì)和反硝化速率測定

根據(jù)《土壤環(huán)境質(zhì)量 農(nóng)用地土壤污染風(fēng)險(xiǎn)管控標(biāo)準(zhǔn)(試行)》(GB 15618—2018)[29]，測定土壤的理化性質(zhì)，包括含水率、pH、Eh，溶解性總固體(total dissolved solids，TDS)、總有機(jī)碳(TOC)、土壤有機(jī)質(zhì)(organic matter，OM)、HCO3—、NH4+、NO3—、NO2—、Ca、K、Mg、砂粒(sand)、粉砂粒(silt)和黏粒(clay)的含量. 反硝化速率測定方法參照文獻(xiàn)[28]：取風(fēng)干10 d后的20個樣品置于試管中，加入20 mg/L的KNO3和去離子水，條件為厭氧條件和25 ℃. 在第0、1、2、5、8、11、15、20、25、28、30、33、35天使用Smart Chem 200(法國Alliance公司)測定上清液的硝酸鹽濃度，通過一級動力學(xué)方程〔見式(1)〕擬合獲得兩個剖面的反硝化速率(denitrification).

式中：k1表示一級反硝化反應(yīng)強(qiáng)度，mg/(kg·d)；c0和ct分別表示初始時刻和t時刻液相中的硝酸鹽濃度，mg/L；t表示反應(yīng)時間，d.

1.3 DNA提取和基因測序

采用土壤power DNA提取試劑盒(MoBio Laboratories，Carlsbad，CA，USA)提取樣品總DNA，利用純化試劑盒(BioTeke Corporation，北京)對粗提總DNA進(jìn)行純化，采用分光光度計(jì)(NanoDrop Technologies，Wilmington，DE，USA)檢測DNA濃度.純化得到的DNA稀釋至10 ng/μL，以降低抑制聚合酶鏈?zhǔn)椒磻?yīng)的可能性，然后將DNA樣品保存于—20 ℃環(huán)境中. 將提取的樣品的DNA作為模板進(jìn)行PCR擴(kuò)增，功能性反硝化基因的PCR反應(yīng)引物如表1所示.總擴(kuò)增體積為20 μL，包括5 μL 5×Reaction Buffer、5 μL 5×High GC Buffer、0.5 μL dNTP(10 mmol/L)、1 μL模板DNA(10 ng/μL)、正向引物和反向引物各1 μL(10 μmol/L)、0.25 μL Q5 high-fidelity DNA polymerase和11.25 μL ddH2O[26]. 最后，根據(jù)測序要求，將純化的擴(kuò)增產(chǎn)物按等摩爾比混合，在中國上海派森諾Personalbio公司的Illumina MiSeq(2×300)平臺上進(jìn)行測序分析.

表1 功能基因PCR擴(kuò)增引物Table 1 primers for PCR amplification of functional genes

在FunGene平臺[33]對質(zhì)控后的序列[34]進(jìn)行翻譯，去除不能翻譯成亞硝酸還原酶和一氧化二氮還原酶的序列. 使用QIIME(v5.2.236，http://www.drive5.com/usearch)去除嵌合體序列和無效序列(＜150 bp的序列、8個以上重復(fù)堿基的序列以及包含模糊堿基的序列)的PCR擴(kuò)增子[35]，高質(zhì)量序列的數(shù)量如圖2所示. 高質(zhì)量序列被UCLUST按相似度97%聚類成Operational Taxonomic Unit(OTU)[36]，并且所有樣本中OTU含量＜0.001%的序列被丟棄. 使用BLAST的方法對比序列數(shù)據(jù)庫進(jìn)行比對注釋，得到每個OTU分類學(xué)信息. 對于不同類別的序列，分別采用各自特定的數(shù)據(jù)庫作為OTU分類地位鑒定的模板序列：細(xì)菌和古菌的16S rRNA基因—Greengenes數(shù)據(jù)庫(Release 13.8，http://greengenes.secondgenome.com)；真菌的18S rRNA基因—Silva數(shù)據(jù)庫(Realease119，http://www.arb-silva.de)；真菌ITS序列的基因—UNITE數(shù)據(jù)庫(Release 5.0，https://unite.ut.ee). 為了最小化樣本間的測序深度差異，在R語言中使用vegan包進(jìn)行抽平分析，后續(xù)多樣性分析基于抽平后的OTU表.

圖2 微生物樣品的高通量測序結(jié)果和α多樣性分析Fig.2 High-throughput sequencing results and α-diversity analysis of microbial samples

1.4 反硝化菌多樣性分析和群落結(jié)構(gòu)分析

基于MOTHUR(http://www.mothur.org)分析反硝化菌的α多樣性：Chao1豐富度指數(shù)；基于豐度的覆蓋度估計(jì)值(ACE)；Shannon-Wiener和Simpson指數(shù).基于Bray-Curtis距離，采用主坐標(biāo)分析(principal co-ordinates analysis，PCoA)反硝化菌屬水平上的群落組成結(jié)構(gòu). 使用R軟件進(jìn)行Adonis分析和Anosim相似度分析，并作999次置換檢驗(yàn)確定組間差異是否具有統(tǒng)計(jì)學(xué)顯著性. 使用R軟件，對豐度前20位的屬進(jìn)行聚類分析并繪制熱圖.

1.5 相關(guān)性分析

采用Spearman非參數(shù)分析方法研究微生物α多樣性與土壤理化性質(zhì)、環(huán)境因子的相關(guān)性，通過冗余分析(Redundancy analysis，RDA)方法、隨機(jī)森林分析方法和分析微生物群落結(jié)構(gòu)與土壤理化性質(zhì)、環(huán)境因子的相關(guān)性. 其中，Spearman非參數(shù)分析使用R軟件中Hmisc、vegan和reshape2包進(jìn)行；RDA分析、隨機(jī)森林和Mantel Test通過R軟件中vegan、dplyr等包進(jìn)行，使用anova函數(shù)評估統(tǒng)計(jì)學(xué)意義. 在R語言中繪圖使用ggplot2包進(jìn)行.

2 結(jié)果與討論

2.1 采樣點(diǎn)不同深度的理化性質(zhì)和反硝化速率

潮白河沖洪積扇S6和S8剖面不同深度的理化性質(zhì)和反硝化速率見表2. 結(jié)果顯示：包氣帶土樣環(huán)境為弱堿性，含水率隨深度的增加而逐漸降低，OM、TOC和TDS含量均呈現(xiàn)隨深度的增加而逐漸降低的趨勢. 從土壤質(zhì)地來看，兩個剖面的黏粒含量較低，主要由砂粒組成. S6和S8兩個剖面的包氣帶以氧化環(huán)境(Eh＞400 mV)為主. S6剖面硝酸鹽、氨氮和亞硝酸鹽含量總體高于S8剖面，S8剖面的HCO3—含量高于S6剖面. 在兩個包氣帶剖面深度不超過10.0 m的區(qū)域內(nèi)，反硝化速率隨著深度的變化呈現(xiàn)出先下降后升高再下降的趨勢，推測包氣帶剖面反硝化作用和硝化作用交替進(jìn)行，反硝化速率可能與氧化還原環(huán)境及反硝化微生物的分離富集有關(guān).

表2 S8和S6剖面不同深度的土壤理化性質(zhì)和反硝化速率[28]Table 2 Soil physical and chemical properties and denitrification intensity at different depths in S8 and S6 profiles[28]

2.2 反硝化菌多樣性和反硝化菌群落結(jié)構(gòu)

在Illumina Miseq測序平臺上，通過對低質(zhì)量read進(jìn)行篩選，得到包氣帶2個剖面(不同深度，共20個包氣帶土壤樣品)共669 762條nirS序列、631 917條nirK序列和843 648條nosZ序列，序列長度大部分在300～400 bp之間. 圖2展示了測序結(jié)果和α多樣性的不同指標(biāo)，包括序列數(shù)、屬水平上的物種數(shù)、ACE、Chao1、Shannon-Wiener和Simpson指數(shù). 含有nirS基因的微生物α多樣性的垂向規(guī)律如下： S8剖面在0.4 m處ACE指數(shù)(1 281.71)較高，隨后下降，在1.0 m處略有回升(1 020.64)，在3.0 m出現(xiàn)最小值(300.47)，4.0 m處達(dá)到峰值(1 606.09)；S6剖面總體的ACE指數(shù)隨深度變化而不斷波動，在0.8、1.2、3.0、5.0 m處均超過1 150，在0.6、1.0、2.0、4.0 m處均不足750. 含有nirK基因的微生物α多樣性的垂向規(guī)律如下：S8剖面的ACE指數(shù)在3 m處最低(1 482.05)，隨后增加，在4.0 m處趨于穩(wěn)定(2 268.25)；S6剖面的ACE指數(shù)在0～1 m內(nèi)出現(xiàn)波動變化，1.0～3.0 m內(nèi)出現(xiàn)明顯增加(1 905.14～3 576.27)，3.0～5.0 m內(nèi)開始降低，隨后5～10 m內(nèi)達(dá)到穩(wěn)定(2 415.59～2 447.46).含有nosZ基因的微生物α多樣性的垂向規(guī)律如下：S8剖面ACE指數(shù)隨包氣帶深度變化顯示出先減少后增加再減少的規(guī)律，0.2 m(1 068.37)和4.0 m(868.3)處出現(xiàn)峰值；S6剖面的ACE指數(shù)整體規(guī)律也為先減少后增加再減少，在0.8 m(1 045.43)和5.0 m(682.33)處出現(xiàn)峰值.

總體而言，S6剖面含有nirS和nirK基因的微生物類群屬的平均數(shù)量比S8剖面高出7.97%和17.41%，而在含有nosZ基因的微生物屬類別數(shù)量上，S6比S8剖面低7.23%. S8和S6剖面微生物α多樣性的中位數(shù)均表現(xiàn)出nirK＞nirS＞nosZ的規(guī)律. 在考察不同人工林間反硝化基因的多樣性也發(fā)現(xiàn)，nirK的25～67個OTU遠(yuǎn)高于nosZ的7～16個和2～14個的nirS[37].S8和S6剖面反硝化菌的α多樣性和包氣帶深度為非線性的關(guān)系，S6剖面多樣性變化相比S8剖面更加明顯. 含有nirK基因的微生物α多樣性指標(biāo)中， S8和S6剖面在包氣帶3.0 m埋深時ACE指數(shù)變化的趨勢不同，但在5.0 m處都趨近于同一個值2 500. 含有nosZ基因的微生物α多樣性指標(biāo)中，雖然S6剖面ACE指數(shù)出現(xiàn)波動，但整體規(guī)律與S8剖面相同，也表現(xiàn)出先減少后增加再減少的規(guī)律，這與反硝化速率的結(jié)果相吻合(兩個剖面反硝化速率在0～10 m深度呈現(xiàn)先下降后升高再下降的變化規(guī)律)，表明反硝化微生物群落的多樣性可以指示反硝化強(qiáng)度的大小，同時包氣帶的反硝化強(qiáng)度一定程度上可以表征包氣帶對地下水NO3—的防護(hù)能力[7]，研究反硝化菌多樣性對進(jìn)一步指示相關(guān)包氣帶的反硝化強(qiáng)度和對包氣帶氮素研究具有重要意義.

圖3展示了nirS、nirK和nosZ基因基于Bray-Curtis的PCoA以及RDA分析結(jié)果，置信區(qū)間為95%.基于Bray-Curtis距離的PCoA分析〔見圖3(a)(c)(e)〕結(jié)果顯示，nirS、nirK和nosZ的前兩個主坐標(biāo)合并分別解釋了OTU群落變化的34.54%、46.68%和56.25%.基于Adonis和Anosim分析，S6和S8剖面在統(tǒng)計(jì)學(xué)上具有分組的意義(P＜0.05)，而采樣點(diǎn)的垂直差異性不明顯(P＞0.05). 在解析崗南水庫微生物群落時空分布特征時，發(fā)現(xiàn)不同深度采樣點(diǎn)基于Adonis(R2=0.166 0，P＜0.001)和MRPP(P＜0.05)分析表現(xiàn)出顯著的垂直差異性[13]. 推測因本研究僅關(guān)注3種反硝化菌，所以群落結(jié)構(gòu)的垂直差異并不明顯，需要明確反硝化菌在不同包氣帶深度下的種類.

圖3 三種基因基于Bray-Curtis的細(xì)菌群落組成主坐標(biāo)分析(PCoA)以及生態(tài)冗余分析(RDA)Fig.3 Principal coordinate analysis (PCoA) and ecological redundancy analysis (RDA) of bacterial community composition based on Bray-Curtis

圖4為nirS、nirK和nosZ型反硝化菌在屬水平上的物種豐度，展示的物種不包含未定義(Incertae_sedis)和未注釋(No blast hit Other)的菌屬，且進(jìn)行了歸一化操作. S6剖面含有nirS基因豐度較高的菌屬包括Kocuria、Bos、Pseudogulbenkiania、Pseudomonas，這些菌屬在深度＞1.2 m的土壤中出現(xiàn)概率較大；而S8剖面在0.2、1.0、3.0 m處豐度較高的菌屬有Sphingomonas、Microbacterium、Caulobacter. S6剖面含有nirK基因的菌屬在3.0 m、5.0 m、10.0 m處的群落組成高度相似，主要來源于Sphingomonas、Halalkalicoccus、Geobacter和Frankia. 而S8剖面相較S6剖面菌屬的豐度減少，其主要來源于Aeromicrobium、Delftia和Rhodococcus. 相比于nirS和nirK，含有nosZ基因的菌屬種類較少，S8剖面的nosZ基因分別來源于Rhodanobacter、Geobacter、Thauera、Modestobacter；而S6剖面的nosZ基因分別來源于Achromobacter、Acidovorax、Roseiflexus.

圖4 結(jié)合聚類分析的屬水平群落組成熱圖Fig.4 Heat map of community composition at genus level combined with cluster analysis

2.3 反硝化菌多樣性、群落結(jié)構(gòu)和環(huán)境因子的相關(guān)性

根據(jù)方差膨脹因子(variance inflation factor，VIF)和前向選擇，在R語言中編輯代碼篩選出共線性的環(huán)境因子，包括Eh、含水率、TDS、砂粒和粉砂粒；保留非共線的環(huán)境因子，包括深度(VIF=3.299 204)、TOC(VIF=3.934 66)、OM(VIF=7.782 175)、HCO3—(VIF=3.551 133)、NH4+(VIF=6.421 39)、NO3—(VIF=2.653 721)、NO2—(VIF=1.238 354)、pH(VIF=6.605 247)、Ca(VIF=4.207 441)、K(VIF=5.213 417)、Mg(VIF=2.347 984)、黏粒(VIF=3.016 942)、反硝化速率(VIF=6.055 176)，后續(xù)的相關(guān)性分析基于篩選過后的環(huán)境因子[13].

Spearman分析結(jié)果表明環(huán)境因子與微生物多樣性之間存在顯著相關(guān)關(guān)系(見圖5). 在S8剖面，pH與nirK基因的ACE指數(shù)呈顯著正相關(guān)(R=0.80，P=0.009 2)，與屬水平上的物種數(shù)呈顯著正相關(guān)(R=0.74，P=0.022)；NH4+含量與nirS基因的ACE指數(shù)呈顯著正相關(guān)(R=0.73，P=0.024 5)，與Shannon-Wiener指數(shù)呈顯著正相關(guān)(R=0.82，P=0.007 2)，與屬水平的物種數(shù)呈顯著正相關(guān)(R=0.71，P=0.031 6)；NO2—含量與nirS基因的ACE指數(shù)呈顯著負(fù)相關(guān)(R=—0.69，P=0.041 2)；NO3—含量與nirK基因的Shannon-Wiener指數(shù)呈顯著負(fù)相關(guān)(R=—0.75，P=0.019 94)；包氣帶深度與nosZ基因的ACE指數(shù)呈顯著負(fù)相關(guān)(R=—0.65，P=0.058 0)，Shannon-Wiener指數(shù)呈顯著負(fù)相關(guān)(R=—0.72，P=0.029)，與nirK基因的Shannon-Wiener指數(shù)呈顯著正相關(guān)(R=0.75，P=0.019 9). 在S6剖面，pH與nosZ基因的ACE指數(shù)呈顯著正相關(guān)(R=0.83，P=0.001 4)，與Shannon-Wiener指數(shù)呈顯著正相關(guān)(R=0.81，P=0.094 6)，與屬水平上的物種數(shù)呈顯著正相關(guān)(R=0.61，P=0.045 3)；黏粒含量與nosZ基因的ACE指數(shù)呈顯著正相關(guān)(R=0.8，P=0.002 9)；包氣帶深度與nosZ基因的ACE指數(shù)呈顯著負(fù)相關(guān)(R=—0.71，P=0.014 5)，與Shannon-Wiener指數(shù)呈顯著負(fù)相關(guān)(R=—0.7，P=0.016 4).

圖5 微生物α多樣性與環(huán)境因子之間相關(guān)性Fig.5 correlation between microbial α diversity and environmental factors

RDA分析〔見圖3(b)(d)(f)〕可以在OTU等級上揭示土壤環(huán)境因子與nirS、nirK和nosZ型反硝化菌群落結(jié)構(gòu)的關(guān)系[9]，同時得到土壤環(huán)境因子對OTU變化的解釋比例. 本研究中，環(huán)境因子對nirS、nirK和nosZ基因OTU的變化解釋比例分別為67.27%、66.33%和80.55%；其中對nirS基因OTU解釋比例＞5%的環(huán)境因子包括反硝化速率(7.82%)、深度(7.30%)、TOC含量(7.14%)、HCO3—含量(6.77%)、Ca含量(6.50%)和OM含量(5.37%)，對nirK基因OTU解釋比例＞5%的環(huán)境因子包括TOC含量(14.83%)、pH(10.65%)、黏粒含量(5.71%)、深度(5.58%)，對nosZ基因OTU解釋比例＞5%的環(huán)境因子包括OM含量(19.19%)、TOC含量(13.04%)、反硝化速率(8.57%)、HCO3—含量(7.89%)、NO3—含量(7.86%). 隨機(jī)森林檢驗(yàn)(見圖6)結(jié)果表明，環(huán)境因子分別解釋了nirS、nirK和nosZ微生物多維尺度分析(multidimensional scaling，MDS)中NMDS1的51.20%、41.06%和35.43%. 結(jié)果表明，nirS和nosZ微生物群落的β多樣性變異主要由包氣帶深度驅(qū)動，nirK微生物群落的β多樣性變異主要由Eh來驅(qū)動.

圖6 隨機(jī)森林解析潮白河包氣帶反硝化基因β-多樣性的環(huán)境驅(qū)動因子Fig.6 Random forest analysis of environmental drivers of denitrification gene β-diversity in Chaobai River vadose zone

分析氣泡圖(見圖7)展示了Mantel Test分析的結(jié)果，即環(huán)境因子Euclidean矩陣和OTU基于Bray-Curtis矩陣的相關(guān)性. S8剖面上，0.2～5.0 m范圍內(nèi)深度(R=0.849)、OM含量(R=0.599)與nirK型反硝化菌群落結(jié)構(gòu)呈顯著正相關(guān)；0.2～1.0 m范圍內(nèi)，顯著正相關(guān)的環(huán)境因子增多，黏粒含量(R=0.734)、深度(R=0.915)、pH(R=0.66)與nosZ型反硝化菌群落結(jié)構(gòu)呈顯著正相關(guān)，K含量(R=0.754)、黏粒含量(R=0.663)分別與nirK、nirS型反硝化菌群落結(jié)構(gòu)呈顯著正相關(guān)；1.0～5.0 m范圍內(nèi)，Ca、Mg含量與nirK、nirS、nosZ型反硝化菌群落結(jié)構(gòu)均呈顯著正相關(guān).S6剖面上，0.2～10.0 m范圍內(nèi)黏粒含量(R=0.664)、pH(R=0.643)均與nosZ型反硝化菌群落結(jié)構(gòu)呈顯著正相關(guān)；0.8～2.0 m范圍內(nèi)，與3種反硝化菌群落結(jié)構(gòu)呈顯著正相關(guān)的環(huán)境因子數(shù)量最多，其中與nosZ型反硝化菌群落結(jié)構(gòu)相關(guān)性最高的環(huán)境因子為pH(R=0.929)，與nirS型反硝化菌群落結(jié)構(gòu)相關(guān)性最高的環(huán)境因子為深度(R=0.821)，與nirK型反硝化菌群落結(jié)構(gòu)相關(guān)性最高的環(huán)境因子為反硝化速率、NH4+含量和pH(R均為0.679)；2.0～5.0 m范圍內(nèi)反硝化速率(R=0.745)與nirK型反硝化菌群落結(jié)構(gòu)呈顯著正相關(guān)，Mg含量(R=0.779)與nirS型反硝化菌群落結(jié)構(gòu)呈顯著正相關(guān)，HCO3—含量(R=0.821)、OM含量(R=0.745)、pH(R=0.78)、TOC含量(R=0.805)與nosZ型反硝化菌群落結(jié)構(gòu)均呈顯著正相關(guān).

圖7 反硝化微生物群落結(jié)構(gòu)與環(huán)境因子間顯著相關(guān)的氣泡圖(P＜0.05)Fig.7 Bubble chart of significant correlation between denitrifying microbial community structure and environmental factors (P＜0.05)

不同剖面環(huán)境因子和微生物多樣性的相關(guān)性程度不同，表明反硝化基因存在水平空間差異性. 一般來說，氮循環(huán)微生物的數(shù)量隨著土壤剖面深度的增加而減少，氮循環(huán)微生物群落的多樣性和結(jié)構(gòu)在深層和淺層土壤也發(fā)生變化[17,38]. 在S6和S8剖面，nosZ型反硝化菌多樣性隨深度的增加而減少，Liu等[18]研究得出nosZ菌的豐度在較深的土壤中更高，推測深度變化會篩選掉一部分的nosZ型反硝化菌，并出現(xiàn)相關(guān)優(yōu)勢菌.nifH、AOA和nirS多樣性指數(shù)在0～60 cm顯著下降，60～80 cm不顯著[17]. Aamer等[39]發(fā)現(xiàn)土壤有機(jī)碳可以通過提高土壤pH的方式增加nosZ和nirK的豐度，進(jìn)而控制N2O的排放. 周夢娟等[40]發(fā)現(xiàn)通過碳源富集培養(yǎng)后的微生物多樣性有所下降，然而發(fā)揮反硝化作用的微生物豐度有所增加. 基于ANCOVA分析發(fā)現(xiàn)nirS基因豐度受土壤pH影響，細(xì)菌群落中nirK基因比例與土壤NO3—濃度有關(guān)，但是不同濕地類型土壤中的反硝化相關(guān)的微生物對相同的環(huán)境變量的反應(yīng)并不相似[41].

整體情況下，環(huán)境因子和微生物群落結(jié)構(gòu)呈現(xiàn)正相關(guān)的關(guān)系，且在不同深度范圍下，顯著相關(guān)的環(huán)境因子數(shù)量不同、類型不同，表明不同深度層次下微生物群落結(jié)構(gòu)差別很大. 冗余分析方法、隨機(jī)森林分析方法和曼特爾檢驗(yàn)分析三種方法研究微生物群落結(jié)構(gòu)與土壤理化性質(zhì)、環(huán)境因子的相關(guān)性結(jié)果相吻合，影響反硝化微生物群落結(jié)構(gòu)的主要環(huán)境因子包括包氣帶深度(Depth)、氧化還原電位(Eh)、土壤有機(jī)質(zhì)含量(OM)、土壤黏粒含量、pH等. 包氣帶深度對微生物群落結(jié)構(gòu)有直接控制作用：隨機(jī)森林分析得出，包氣帶深度驅(qū)動nirS和nosZ微生物群落的β多樣性變異最強(qiáng)；RDA分析得出，包氣帶深度對nirK微生物群落解釋比例為5.58%. Fierer等[42]研究表明土壤深度與微生物生物量的變化在兩個剖面中均呈顯著相關(guān)（P＜0.001），特定微生物群體的垂直分布主要可以歸因于土壤C的可利用性隨著土壤深度的增加而下降. 深度對反硝化菌造成影響的深層次原因需要進(jìn)一步發(fā)掘.同時酸化土壤已被證明會破壞土壤深度與細(xì)菌群落之間的關(guān)系[43-44]. 在本研究中，pH在不同的土壤深度對nosZ型反硝化菌群落結(jié)構(gòu)均呈顯著相關(guān)：0.2～1.0 m深度下，R=0.66；1.0～5.0 m深度下，R=0.643.與本研究結(jié)果(pH對nirS和nirK型反硝化菌群落結(jié)構(gòu)相不具有相關(guān)性)不同，Xiao等[21]發(fā)現(xiàn)，土壤pH解釋了nirS群落31.0%的變化(P＜0.05)，解釋了nirK群落34.0%的變化(P＜0.05)，推測該研究中nirS和nirK型反硝化菌群落結(jié)構(gòu)受到弱堿性土壤的影響很小，而nosZ型反硝化菌對pH變化更加敏感. 本研究使用Mantel Test分析探究基于Bray-Curtis矩陣的反硝化菌群落結(jié)構(gòu)和環(huán)境因子的關(guān)系，發(fā)現(xiàn)nosZ型反硝化菌群落和黏粒含量呈顯著正相關(guān)(P＜0.05)，與已有研究結(jié)果[35]相吻合.

本研究探究了不同包氣帶深度下反硝化菌的多樣性和群落結(jié)構(gòu)，并分析了它們與環(huán)境因素的相關(guān)性. 已有研究僅關(guān)注淺層(＜2.0 m)土壤，而本研究揭示了土壤深層(0.2～10.0 m)中反硝化作用的主控因素. 在對包氣帶進(jìn)行反硝化過程和氮循環(huán)研究時，應(yīng)考慮不同深度下環(huán)境因子的影響：在淺層包氣帶埋深下，多考慮土壤結(jié)構(gòu)、pH、亞硝酸鹽含量的影響，在深層包氣帶埋深下，多考慮土壤有機(jī)質(zhì)、金屬元素和硝酸鹽含量的影響. 本研究分析的包氣帶樣品和參數(shù)有限，反硝化微生物和土壤之間復(fù)雜的作用關(guān)系需要更進(jìn)一步深入研究.

3 結(jié)論

a) 總體來說，S6剖面的α多樣性整體高于S8剖面，同時不同類型的反硝化菌α多樣性表現(xiàn)出一定的規(guī)律：nirK＞nirS＞nosZ. S8和S6剖面的α多樣性和包氣帶深度為非線性的關(guān)系，S6剖面α多樣性變化相比S8剖面更加明顯. 反硝化微生物α多樣性存在垂直空間差異性和水平空間差異性.

b) 對微生物群落結(jié)構(gòu)的分析表明，nirS、nirK和nosZ基因主要來源于變形菌門和放線菌門中參與氮循環(huán)的屬，如固氮菌屬和根瘤菌屬. S8剖面相較于S6剖面，含有nirS、nirK基因的菌屬豐度較少，而且兩個剖面之間菌屬種類存在差異. 而含有nosZ基因的菌屬在兩個剖面的豐度均較低.

c) 不同反硝化菌α多樣性和環(huán)境因子的相關(guān)程度不同. 影響潮白河反硝化微生物群落結(jié)構(gòu)的主要環(huán)境因子包括包氣帶深度、氧化還原電位(Eh)、包氣帶有機(jī)質(zhì)含量(OM)、包氣帶黏粒含量、pH等. 不同深度范圍內(nèi)，影響反硝化菌群落結(jié)構(gòu)的環(huán)境因子數(shù)量和種類存在較大的差異. 因此，在進(jìn)行包氣帶的反硝化情況和反硝化菌群落功能結(jié)構(gòu)調(diào)查時，需要根據(jù)不同深度關(guān)注不同的環(huán)境因子. 對于淺層包氣帶，應(yīng)重點(diǎn)應(yīng)關(guān)注黏粒含量、深度、pH、碳酸鹽和亞硝酸鹽含量等環(huán)境因子；而對于深層包氣帶，則需要關(guān)注鈣、鎂、硝酸鹽和總有機(jī)碳含量等環(huán)境因子.