包氣帶反硝化菌多樣性和群落結構垂向分布規律與主控環境因子研究：以北京市潮白河沖洪積扇為例

晁韶良，郇 環，周向同，劉 媚

1. 中國環境科學研究院，北京 100012

2. 生態環境部土壤與農業農村生態環境監管技術中心，北京 100012

3. 江蘇大學環境與安全工程學院，江蘇 鎮江 212013

4. 中國地質大學（武漢）環境學院，湖北 武漢 430078

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

—向地下水遷移有重要影響，其強度可以指示包氣帶的NO3—防護效果[6-7]. 了解包氣帶反硝化微生物的垂直多樣性和豐度對于提高作物產量、應對環境污染和調節生物地球化學氮循環至關重要.

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

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

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

1 材料與方法

1.1 研究區概況

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

本研究選取潮白河沖洪積扇中下部典型包氣帶作為研究區域，并于2016年10月27日—2016年11月4日使用PowerProbe鉆機(9500-VTR，美國AMS公司)在S6和S8剖面(見圖1)采樣. S6剖面共11個樣品，依次記為S61、S62、S63、S64、S65、S67、S68、S69、S610、S611，分別對應埋藏深度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.采集的土壤樣品應避免人員接觸污染，采集后密封，在4 ℃環境下保存，標注好樣品的采樣時間、地點和采樣深度.

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

1.2 土壤理化性質和反硝化速率測定

根據《土壤環境質量 農用地土壤污染風險管控標準(試行)》(GB 15618—2018)[29]，測定土壤的理化性質，包括含水率、pH、Eh，溶解性總固體(total dissolved solids，TDS)、總有機碳(TOC)、土壤有機質(organic matter，OM)、HCO3—、NH4+、NO3—、NO2—、Ca、K、Mg、砂粒(sand)、粉砂粒(silt)和黏粒(clay)的含量. 反硝化速率測定方法參照文獻[28]：取風干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公司)測定上清液的硝酸鹽濃度，通過一級動力學方程〔見式(1)〕擬合獲得兩個剖面的反硝化速率(denitrification).

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

1.3 DNA提取和基因測序

采用土壤power DNA提取試劑盒(MoBio Laboratories，Carlsbad，CA，USA)提取樣品總DNA，利用純化試劑盒(BioTeke Corporation，北京)對粗提總DNA進行純化，采用分光光度計(NanoDrop Technologies，Wilmington，DE，USA)檢測DNA濃度.純化得到的DNA稀釋至10 ng/μL，以降低抑制聚合酶鏈式反應的可能性，然后將DNA樣品保存于—20 ℃環境中. 將提取的樣品的DNA作為模板進行PCR擴增，功能性反硝化基因的PCR反應引物如表1所示.總擴增體積為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]. 最后，根據測序要求，將純化的擴增產物按等摩爾比混合，在中國上海派森諾Personalbio公司的Illumina MiSeq(2×300)平臺上進行測序分析.

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

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

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

1.4 反硝化菌多樣性分析和群落結構分析

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

1.5 相關性分析

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

2 結果與討論

2.1 采樣點不同深度的理化性質和反硝化速率

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

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

2.2 反硝化菌多樣性和反硝化菌群落結構

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

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

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

圖3 三種基因基于Bray-Curtis的細菌群落組成主坐標分析(PCoA)以及生態冗余分析(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)的菌屬，且進行了歸一化操作. S6剖面含有nirS基因豐度較高的菌屬包括Kocuria、Bos、Pseudogulbenkiania、Pseudomonas，這些菌屬在深度＞1.2 m的土壤中出現概率較大；而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 結合聚類分析的屬水平群落組成熱圖Fig.4 Heat map of community composition at genus level combined with cluster analysis

2.3 反硝化菌多樣性、群落結構和環境因子的相關性

根據方差膨脹因子(variance inflation factor，VIF)和前向選擇，在R語言中編輯代碼篩選出共線性的環境因子，包括Eh、含水率、TDS、砂粒和粉砂粒；保留非共線的環境因子，包括深度(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)，后續的相關性分析基于篩選過后的環境因子[13].

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

圖5 微生物α多樣性與環境因子之間相關性Fig.5 correlation between microbial α diversity and environmental factors

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

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

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

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

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

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

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

3 結論

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

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

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