寒溫帶落葉松天然林不同海拔土壤線蟲群落結構的季節動態

張巖，楊立學

(東北林業大學 林學院，哈爾濱 150040)

0 引言

以氣溫和降水格局變化為代表的氣候變化能夠對土壤環境產生深刻的影響。氣候變化引起土壤微生境變化，包括初級生產力、土壤微生物量、食物來源和土壤水分的變化，是影響土壤線蟲的主要因素[8]。隨后，線蟲群落也會影響土壤生態系統的相關生態功能和食物網[9]。植物寄生線蟲可促進光吸收的碳向根系分配，導致根系分泌物增加和根際微生物活性增加[10]。食細菌線蟲和食真菌線蟲與土壤中營養的礦化相關[11]。Long等[12]對線蟲多樣性指數的研究結果表明，氣候變暖可能會影響土壤線蟲敏感屬的豐度。然而，氣溫升高后，線蟲群落功能的變化不完全等同于土壤線蟲屬的變化。季節變化主要是通過溫度和濕度的變化、食物的供應和種內競爭等對土壤線蟲群落造成影響[13]。Verschoor等[14]認為，土壤線蟲的多樣性和種群結構呈現年周期變化，與濕度和溫度的季節性變化相關。土壤線蟲的豐度大致呈現出春季、冬季、秋季和夏季依次遞減的季節分布規律[15]，夏季降雨量大，土壤含水量和土壤溫度高，導致線蟲數量減少至最低，而春、冬兩季降雨少且溫度適中，適宜土壤線蟲生存和大量繁殖，故使其豐度達到較高的水平。現階段對于不同時間動態的土壤線蟲研究較為廣泛，但是對于山地森林系統土壤線蟲群落的研究為數不多，因此了解不同海拔梯度土壤線蟲群落季節動態變化的時空規律尤為重要。

土壤線蟲的分離和鑒定方法多種多樣[16]，傳統的土壤線蟲形態學分離鑒定法，雖然能比較準確地反映出線蟲群落狀況，但費時費力[17]，且分類比較粗糙，往往只到科屬水平，還容易造成樣本量的缺失。近年來，高通量測序、功能基因芯片等高通量宏基因組學技術以及生物信息分析技術的迅速發展，人們可以通過基因水平對微生物群落組成和多樣性進行比較研究[18]。目前，高通量測序技術也在線蟲研究中得到一定應用，如：Kerfahi等[19]對熱帶地區和北極地區的土壤線蟲群落多樣性的研究均通過DNA測序的方法將土壤線蟲活體作為研究的生物指標；Du等[20]對高緯度階地與森林交錯帶額爾古納濕地的土壤線蟲群落組成進行了研究。Geisen等[17]研究發現，采用通用引物3NDF(5’-GGCAAGTCTGGTGCCAG-3’)和1132modR(5’-TCCGTCAATTYCTTTAAGT-3’)針對18S rDNA可變性最強的V4區，作為真核生物的條形碼區域，可提高分類學水平的分辨率，并提高線蟲分類群組的覆蓋率。結合Illumina MiSeq技術提供的測序深度，得到的結果比以前從形態學和其他針對線蟲的HTS研究中鑒定出線蟲的多樣性要高得多[21-22]。同時，該引物在海拔方面應用較少。

北方地區的森林由于對氣候變化高度敏感，并且面積巨大，對全球氣候系統影響極大，存在極強的交互作用。大興安嶺林區地處北半球高緯度區域，是響應全球氣候變化最敏感的地區之一[23]。興安落葉松(Larixgmelinii)作為我國北方森林主要的建群種和優勢種，對改善東北乃至全國的森林生態環境、促進林業可持續發展具有重要作用[24-25]。為了能夠預測土壤線蟲群落對不同氣候條件的響應機制，了解不同海拔土壤線蟲多樣性是至關重要的[26]。本研究以大興安嶺北部地區寒溫帶不同海拔落葉松天然林為研究對象，利用Illumina MiSeq測序的方法，選擇3NDF和1132 modR作為引物，對寒溫帶不同季節海拔上土壤線蟲群落結構特征及其多樣性進行比較，假設：①土壤線蟲群落沿海拔梯度的分布模式均為隨海拔的升高而降低；②土壤pH和水分是影響不同季節線蟲沿海拔分布的主要土壤因子。本研究的目的是為了了解寒溫帶線蟲群落的海拔分布格局的季節動態，探究在大興安嶺北部地區不同季節變化影響土壤線蟲群落分布的主要土壤因子。

1 研究地概況

本研究位于內蒙古大興安嶺北部阿龍山林業局先鋒林場的奧克里堆山(51°50′N，122°01′E)，屬寒溫帶大陸性氣候，具有明顯的季風特征。夏季受太平洋東南季風的影響，溫暖濕潤；冬季受西伯利亞冷空氣影響，寒冷干燥。晚霜為6月上旬，無霜期為80 d，年平均氣溫-5.1 ℃，極端最高溫度為35.4 ℃，極端最低溫度為-49.6 ℃，年降水量為437.4 mm，降水集中在7月和8月。

奧克里堆山海拔為1 523 m，是大興安嶺主脈北段最高的山峰，植被仍保持完整原始狀態。主要植被為東西伯利亞針葉林向南延伸部分，植被垂直分布較明顯且占絕對優勢，其大致可分為3個植物垂直帶：亞高山矮曲林(海拔為1 300 m以上)；針葉林或針闊混交林帶(海拔為750～1 300 m)；谷底沼澤、灌叢和草甸帶(海拔為750 m左右)，見表1。興安落葉松的分布范圍最廣且面積最大，從海拔700 m到海拔1 350 m均能良好生長。

表1 奧克里堆山不同海拔植物垂直帶變化

2 研究方法

2.1 樣地設置和土樣采集

2018年5月底，以沿海拔上不同物種組成為依據，在奧克里堆山陽坡落葉松天然林中選取6個海拔作為樣點，分別是750、830、950、1 100、1 300、1 420 m。每個海拔分別設置3個20 m×30 m的標準地樣方，進行土壤樣品采集，采用隨機點取樣法在每個大樣方內選取8個點，去除地表凋落物并挖取土壤剖面，使用土鉆采集剖面上0～10 cm層的土壤樣品，混合均勻。將每個海拔、每塊樣地的部分土壤樣品裝入便攜式保溫箱中帶回實驗室分別于4 ℃保存(用于速效養分等指標測定)、-80 ℃條件下保存(用于土壤線蟲高通量測序)；其余的樣品在實驗室進行風干處理，用于土壤理化性質含量的測定。

2.2 土壤理化性質測定

2.3 DNA抽提和PCR擴增

根據E.Z.N.A.?soil試劑盒(Omega Bio-tek, Norcross, GA, U.S.)說明書進行總DNA抽提，DNA濃度和純度利用NanoDrop2000進行檢測，利用1%瓊脂糖凝膠電泳檢測DNA提取質量；用3NDF(5’-GGCAAGTCTGGTGCCAG-3’)和1132 modR(5’-TCCGTCAATTYCTTTAAGT-3’)引物對18S rDNA的最可變部分V4可變區進行PCR擴增，擴增程序為：95 ℃預變性3 min，27個循環(95 ℃變性30 s，55 ℃退火30 s，72 ℃延伸30 s)，最后72 ℃延伸10 min(PCR儀：ABI GeneAmp?9700型)。擴增體系為20 μL，4 μL 5×Fast Pfu緩沖液，2 μL 2.5 m mol/L dNTPs，0.8 μL引物(5 μmol/L)，0.4 μL Fast Pfu聚合酶；10 ng DNA模板。

2.4 Illumina MiSeq測序

使用2%瓊脂糖凝膠回收PCR產物，利用AxyPrep DNA Gel Extraction Kit(Axygen Biosciences, Union City, CA, USA)進行純化，Tris-HCl洗脫，2%瓊脂糖電泳檢測。利用QuantiFluorTM-ST(Promega, USA)進行檢測定量。根據Illumina MiSeq平臺(Illumina, San Diego, USA)標準操作規程將純化后的擴增片段構建PE2×300的文庫。構建文庫步驟: ① 連接“Y”字形接頭；② 使用磁珠篩選去除接頭自連片段；③利用PCR擴增進行文庫模板的富集；④ 氫氧化鈉變性，產生單鏈DNA片段。利用Illumina公司的MiSeq PE300平臺進行測序(上海美吉生物醫藥科技有限公司)。原始數據上傳至NCBI數據庫中。

2.5 數據處理

原始測序序列使用Trimmomatic軟件質控，使用FLASH軟件進行拼接。

(1)設置50 bp的窗口，如果窗口內的平均質量值低于20，從窗口開始截去后端堿基，去除質控后長度低于50 bp的序列。

(2)根據重疊堿基overlap將2端序列進行拼接，拼接時overlap之間的最大錯配率為0.2，長度需大于10 bp。

(3)根據序列首尾2端的barcode和引物將序列拆分至每個樣本，barcode需精確匹配，引物允許2個堿基的錯配，去除存在模糊堿基的序列。

使用的UPARSE軟件(version 7.1 http://drive5.com/uparse/)，根據97%的相似度對序列進行OTU聚類；使用UCHIME軟件剔除嵌合體。利用RDP classifier(http://rdp.cme.msu.edu/)對每條序列進行物種分類注釋，比對UNITE數據庫(UNITE 8.0)，設置比對閾值為70%。

2.6 數據分析

土壤線蟲按照攝食習慣在屬水平進行營養類群的分類[27]。利用QIIME對Illumina MiSeq測序數據的α多樣性指數進行分析[28]。采用R(version 3.6.1)中的“vegan”包[29]，基于Bray-Curtis距離進行β多樣性主坐標分析(PCoA)，分析線蟲群落的相似性。通過Bray-Curtis距離的相似性分析和排列多元方差分析(PERMANOVA)，研究了不同海拔土壤線蟲群落特征的差異。采用Microsoft Excel 2016和SPSS(vers.23.0; SPSS, Chicago, IL, USA)進行數據處理。利用CANOCO 4.5軟件對群落結構與環境因素的關系進行冗余分析(RDA)。所有數據為平均值±標準誤(SE)。

3 結果與分析

3.1 線蟲測序特征

在Illumina MiSeq平臺對土壤線蟲的3NDF和1 132 modR片段測序，經過優化和抽平，54個土壤樣品共生成738 915條高質量有效序列。每個樣本包含的序列范圍為7 845～21 276條(平均13 684條)。18S rRNA基因區域的平均讀取長度為274 bp，大于Good覆蓋率的99%。在97%相似度下OTU水平上，根據可觀測OTU數量做稀釋曲線，所有海拔和不同土層土壤樣品稀釋曲線趨于平坦，樣品OTU覆蓋度基本飽和，這表明測序深度能夠基本反映土壤線蟲群落結構。如圖1所示。

圖1 土壤線蟲群落的操作分類單位數(OTU)的Rarefaction曲線

3.2 線蟲群落組成

在97%序列相似性水平上，共獲得53個屬和86個OTU分布在3個綱和11個目中。在屬水平方面如下：

5月的優勢屬分別為真單宮屬(Eumonhystera)、巨宮屬(Tylolaimophorus)、棱咽屬(Prismatolaimus)、真茅線屬(Eudorylaimus)和g__unclassified_f__Plectidae，總體相對豐度分別為18.92%、10.44%、7.55%、6.29%和6.06%，其中，真單宮屬(Eumonhystera)在海拔830、950、1 100 m均有較高的相對豐度，巨宮屬(Tylolaimophorus)在海拔750 m和1 300 m有較高相對豐度。

7月的優勢屬分別為棱咽屬(Prismatolaimus)、真單宮屬(Eumonhystera)、真茅線屬(Eudorylaimus)、巨宮屬(Tylolaimophorus)和g__unclassified_f__Plectidae，總體相對豐度分別為12.68%、11.34%、9.61%、9.24%和6.50%，其中，棱咽屬(Prismatolaimus)在830 m和1 420 m最大，且超過這2個海拔相對豐度的四分之一，巨宮屬(Tylolaimophorus)在750、1 100、1 300 m最大。

9月的優勢屬分別為巨宮屬(Tylolaimophorus)、真單宮屬(Eumonhystera)、Bunonema、棱咽屬(Prismatolaimus)和Teratocephalus，總體相對豐度分別為19.26%、9.76%、8.14%、7.34%和5.51%，其中，巨宮屬(Tylolaimophorus)在830、950、1 300 m海拔處僅有較高分布，Bunonema在海拔1 100 m和1 420 m處有最大值。如圖2所示。

圖2 沿海拔梯度不同季節土壤線蟲屬的相對豐度

真單宮屬(Eumonhystera)、巨宮屬(Tylolaimophorus)、棱咽屬(Prismatolaimus)在各個季節的海拔梯度上均有較多分布，因此為該地區的優勢土壤線蟲屬。3個月共測定出53個土壤線蟲屬，其中，食細菌線蟲15屬，食真菌線蟲7屬，植食性線蟲21屬，雜食-捕食性線蟲8屬，未分類線蟲2屬。見表2。

表2 總季節動態土壤線蟲群落組成分類及營養類群c-p值

3.3 土壤線蟲多樣性和群落結構

海拔對土壤線蟲的香農(Shannon)指數存在顯著影響(P<0.05)。所有季節中，Sobs指數(樣本總量的指數)、香農(Shannon)指數、Chao指數(群落豐富度指數之一)和Faith’s PD指數(譜系多樣性的指數)的變化范圍分別為6.33±1.67～112.33±1.45、1.34±0.43～2.35±0.11、8.08±2.89～19.64±6.96和1.99±0.47～3.11±0.28。其中，Sobs指數在3個月中均在最低海拔750 m處最高。9月時，在950 m海拔土壤的Sobs指數、Shannon指數、Chao指數和Faith’s PD指數均為最低，見表3。總體上，隨海拔的升高，土壤線蟲Sobs指數和Faith’s PD指數呈降逐漸低的趨勢，如圖3所示。不同季節間各項α多樣性指數均無顯著影響，但是海拔和季節對Sobs指數和Shannon指數的綜合影響顯著(P<0.05)。

注：ns, 不顯著。

表3 不同季節土壤線蟲群落α多樣性指數

對各季節不同海拔土壤線蟲測序數據進行基于Bray-Curtis(通用的官方距離矩陣算法)距離的主坐標分析，所有季節中，2個主成分軸共解釋了總變異的40%以上。其中，5月2個主成分坐標軸共解釋了總變異的58.26%，其中,PC1(圖中表達方式為PCoA axis1)軸解釋了36.46%，PC2(圖中表達方式為PCoA axis2)軸解釋了21.8%；7月共解釋了總變異的45.76%，其中,PC1軸解釋了27.15%，PC2軸解釋了18.61%；9月共解釋了總變異的51.04%，其中,PC1軸解釋了31.69%，PC2軸解釋了19.35%。不同海拔的樣本都被分為3組，即低海拔(750 m)、中海拔(830、950、1 100、1 300 m)和高海拔(1 420 m)。相似性分析(ANOISM)和置換多元方差分析(PERMANOVA)表明，不同海拔的土壤線蟲群落結構存在顯著差異，如圖4所示。所有樣本的PERMANOVA結果表明，海拔對土壤線蟲群落結構的影響強于季節變化(表4)。

圖4 5月、7月和9月的土壤線蟲群落主坐標分析

表4 基于Bray-Curtis距離的海拔和季節對土壤線蟲群落的非參數多元方差分析(PERMANOVA)

3.4 土壤因子與線蟲群落結構的關系

通過冗余分析評價土壤因子和線蟲群落間的關系。RDA結果表明，所有季節的總解釋量均能達到線蟲群落結構變異的50%以上。其中，5月總解釋量為56.7%，土壤容重是主要的影響因子，土壤容重、全碳、土壤溫度、pH、銨態氮的解釋量分別為9.7%、7.7%、7.2%、6.9%、6.6%；7月總解釋量為54.6%，pH發揮顯著影響作用，pH、土壤溫度、可溶性有機氮、土壤容重、銨態氮的解釋量分別為11.0%、7.9%、6.8%、6.7%、6.2%；9月的總解釋量為64.2%，土壤可溶性有機碳影響顯著，土壤可溶性有機碳、pH、硝態氮、全碳、銨態氮的解釋量分別為13.8%、9.9%、7.9%、7.3%、7.1%，如圖5所示。不同季節間主要影響線蟲群落變化的因子存在差異，但整體上銨態氮、pH、可溶性有機碳和可溶性有機氮在一定程度上與線蟲群落保持較高的相關性。在全部月份中，大多數線蟲屬分布在圖5的中心位置，如：Calcaridorylaimus、Filenchus、Odontolaimus、Sphaeronema等，圖中使用的是前3位字母簡寫。

圖5 基于土壤線蟲屬水平和土壤因子的冗余分析

4 結論與討論

(1)本研究發現，寒溫帶大興安嶺北部奧克里堆山落葉松天然林中土壤線蟲的豐度和多樣性均與海拔呈極顯著負相關，在奧克里堆山沿海拔變化呈單調遞減的模式。該結果與Yeates等[4]的研究一致，類似的結果在捷克共和國的貝斯科迪山脈[30]和斯洛伐克的維奧拉克山脈[31]的落葉林中也得到證實；5月、7月、9月均符合此規律，該結果證實了筆者的第一個假設，即所有季節均表現為土壤線蟲的豐度和多樣性與海拔呈極顯著負相關。相比于季節變化，土壤線蟲群落結構對海拔更加敏感，沿海拔變化的非生物因子主導著土壤線蟲群落的變化，土壤溫度、土壤硝態氮、pH、可溶性有機碳和氮等在各季節中發揮重要作用。在7月份，土壤pH與大部分土壤線蟲群落的變異密切相關。已有研究表明，pH為4～12時，土壤線蟲都能成活[32]，而當pH比較低時，更易于吸引根結線蟲集聚，同樣，土壤pH過高會抑制線蟲生長發育[33]。土壤容重、pH、可溶性有機碳在5月、7月、9月分別發揮著主導作用，而土壤容重和可溶性養分與土壤含水量也與土壤線蟲群落結構密切相關，證實了筆者的第二個猜測。土壤線蟲個體對溫度和降水的響應強烈，土壤線蟲群落組成和功能與土壤食物網的聯系也很重要[34]。在本研究中，土壤溫度沿海拔升高而升高，雖然土壤含水量沒有明顯的下降，但是溫度升高導致土壤線蟲群落豐度和多樣性下降的結果，與之前的研究相一致[35]。但是溫度影響土壤線蟲群落結構變化的機制還需要進一步分析。Thakur等[9]發現，氣候變暖條件下，土壤線蟲群落結構的變化是由植物群落組成和生產力的變化引起的。然而，不同的營養類群對溫度變化的反應不同，食細菌線蟲和植食性線蟲對溫度感應尤其顯著。Zhang等[36]的研究表明，溫度的升高使線蟲群落中食細菌線蟲和植食性線蟲的豐度顯著降低，但食細菌線蟲仍然占主導地位，該結果與本文研究結果一致。升溫對食真菌線蟲豐度的影響較小，表明其對溫度升高具有一定抗性。

(2)本研究發現不同季節間的土壤線蟲的α和β多樣性不存在明顯差異，此結果與呂文娜[37]對黃河三角洲的研究結果不同。可能原因是本研究以食細菌線蟲為主要優勢營養類群，而黃河三角洲地區以植食性線蟲為主，季節變化導致的植被變化影響較大。土壤線蟲群落隨時間的變化不能僅用環境變化來解釋，還需考慮生物相互作用。土壤線蟲的豐度極大程度上也依賴于食物的來源，因此土壤細菌和真菌等豐度也會影響食細菌線蟲和食真菌線蟲的豐度[38-40]。土壤線蟲是在土壤微食物網中極其重要的一環[41]，土壤微食物網具有調節作用，能夠通過縮短食物鏈和選擇更廣泛的食物網的方法，使土壤線蟲群落的變化不會異常劇烈。通過PERMANOVA的結果，進一步驗證了在寒溫帶山地生態系統中，落葉松天然林海拔對土壤線蟲群落結構的影響要強于季節變化。

(3)本研究共測定出53個土壤線蟲屬，其中，食細菌線蟲15屬，食真菌線蟲7屬，植食性線蟲21屬，雜食-捕食性線蟲8屬，未分類線蟲2屬。在5月、7月和9月中，食細菌線蟲總豐度最高，因此食細菌線蟲為優勢營養類群。研究表明，在寒溫帶落葉松天然林中，不同時間動態的土壤線蟲的豐度和多樣性均表現出沿海拔升高逐漸降低的顯著負相關模式，且相較于季節的變化，海拔對土壤線蟲群落結構和多樣性表現出更為顯著的影響。主要受到土壤硝態氮、pH、可溶性有機碳和氮等的影響，其中pH是主要影響因子。