长期的矿区开采改变了原有地貌形态和土壤组成结构, 而且形成排土场占用土地, 对原生态环境造成了毁灭性的破坏^{[1, 2]}.对矿区废弃地进行复垦, 可有效改善土壤基质的结构和肥力, 还能提高土壤中微生物的活性^[3].微生物作为土壤中最活跃的成分, 通过合成土壤酶来参与土壤中的生物化学过程, 在有机质的转化、养分矿化和同质化等方面发挥着重要作用, 影响着土壤生态系统的能量流动和物质循环, 土壤酶也由于其对环境因子的变化敏感而成为土壤环境质量的重要指标^[4~7].功能基因通过编码特定功能蛋白, 参与介导微生物的各项生理生化过程和土壤中的元素循环, 因此具有指示土壤中生命活动的功能, 具有标志性的功能基因也常被用作分子标记物来表征土壤中某一类群微生物的群落结构、丰度和多样性^{[8, 9]}.

有研究表明, 不同植被类型由于凋落物和根际分泌物的性质不同, 对土壤理化和微生物特征均有显著影响, 表现在微生物群落结构、土壤理化性质和土壤酶活差异等方面^[10].目前常用高通量测序表征土壤中微生物多样性和群落差异^[11], 而对不同植被类型对土壤微生物功能基因影响方面研究还较少.本研究在山西省大同市晋华宫矿山公园, 不同植被样地采集土壤样品, 对不同植被类型的土壤微生物碳氮磷硫功能基因丰度进行测定, 分析植被类型对酶活性和碳氮磷硫功能基因丰度的影响, 以及土壤酶活和基因丰度之间的相关性, 通过揭示不同复垦植被类型对土壤质量的影响, 以期为矿区废弃地恢复工程植被选择提供科学依据.

研究区位于大同煤田东北边缘的十里河狭长谷地的晋华宫国家矿山公园, 距离大同市12.5 km, 整体呈由南向北的上升趋势.位于40°10′N, 113°13′E; 平均海拔1 056 m.地处中温带大陆性半干旱季风气候, 全年温差较大, 年均降雨量约为400 mm.本研究包括5个植被均是2005年公园建设时栽种, 栽种前覆盖客土1 m左右, 样地之间距离均在6 m以上, 树苗间距为3 m×3 m, 分别是侧柏(Platycladus orientalis, 盖度65%)、桧柏(Sabina chinensis, 盖度65%)、樟子松(Pinus sylvestris, 盖度80%)、云杉(Picea asperata, 盖度70%)和油松(Pinus tabuliformis, 盖度70%).每个植被复垦类型下, 设立3个重复样地, 每个样地为10 m×10 m.

2019年7月在侧柏、桧柏、油松、樟子松和云杉复垦区的样地分别采集土样, 用土钻采用多点混合土样样品处理法, 去除地表植被和覆盖物, 每个样地分别取6个样点挖穴(0~10 cm)并混匀.一部分带回实验室自然风干, 研磨过2 mm的筛以备用; 另一部分4℃低温保存用于土壤碳氮磷硫功能基因的测定.

土壤pH采用1∶5土和水溶液测定; 土壤全碳、全氮和全硫采用vario MACRO cube元素分析仪测定.土壤微生物碳氮磷硫功能基因丰度采用QMEC(quantitative microbial element cycling)基因芯片测定(广东美格基因科技有限公司), 具体方法是：按照MagaBio土壤DNA提取试剂盒方法, 取0.25 g土壤样品, 提取土壤DNA, 用Qubit 4.0 (Thermo Fisher Scientific, Waltham, USA)仪器对提取的土壤DNA总量及纯度进行检测, 将DNA浓度稀释20 ng·μL^-1, 分别将样品板和引物板试剂添加至高通量qPCR芯片SmartChip MyDesign Chip(Takara Biomedical Technology, Clontech)微孔中, 在SmartChip Real-Time PCR System(WaferGen Biosystems USA)中进行qPCR反应及荧光信号检测.土壤酶活性均采用酶活试剂盒(生工生物工程股份有限公司), 酶活性测定的方法分别是, 土壤碱性磷酸酶活性采用磷酸苯二钠比色法; 土壤脲酶活性采用靛酚蓝比色法; 土壤脱氢酶活性采用三苯基甲臜比色法; 土壤碱性蛋白酶活性采用钨蓝比色法.

采用单因素方差分析(One-way ANOVA)分析不同植被恢复类型对土壤理化性质、土壤酶活性和土壤微生物功能基因丰度的影响; 采用非度量多维尺度分析(NMDS)反映植被类型对微生物群落的影响; 土壤细菌群落功能多样性指数分析采用香农-威纳指数(Shannon-Wiener index)和辛普森指数(Simpson index)^[12].土壤酶活性和土壤微生物功能基因丰度间的相关性进行Pearson相关性分析.数据通过Microsoft Excel 2013和SPSS 24.0进行处理和分析, 作图采用Microsoft Excel 2018和R 3.5.1等完成.

对理化指标和微生物特征进行主成分分析法(PCA), 根据特征向量值识别最重要的指标, 提取特征值大于1的主成分(PC)进行初始解释^[13], 土壤综合肥力指数(integrated fertility index, IFI)具体计算公式是：

式中, W_i为第i个主成分的特征向量值, S_i为第i个变量的指标贡献率, i=1, 2, 3, …, n.

从表 1可以看出, 侧柏容重最高, 而樟子松容重最小, 且与其它3种样地差异显著; 复垦地土壤pH均呈弱碱性, 且不同恢复方式会对土壤pH产生显著的影响.侧柏林土壤pH最高, 云杉最低.桧柏具有最高的总氮量, 其次是樟子松, 桧柏与樟子松差异不显著, 均显著高于其余3种样地, 油松、侧柏和云杉林具有相近的总氮量; 桧柏具有显著高于其它样地的总碳量, 其次是油松和樟子松, 樟子松和油松差异不显著, 侧柏总碳量最低; 桧柏具有最高总硫量, 其次是樟子松, 侧柏和油松总硫量差异不显著, 云杉林总硫量最低.

表 1 (Table 1)

表 1 不同植被恢复类型下的土壤理化性质1) Table 1 Soil physiochemical characteristics under different regeneration scenarios 植被 容重/g·mL-1 pH ω(总氮)/g·kg-1 ω(总碳)/g·kg-1 ω(总硫)/g·kg-1 侧柏 0.47±0.02a 8.17±0.02a 0.50±0.18b 12.39±2.73d 0.64±0.04c 油松 0.41±0.03b 8.12±0.02bc 0.80±0.02b 23.65±0.63b 0.86±0.09c 云杉 0.42±0.02b 8.08±0.01c 0.80±0.02b 17.53±0.83c 0.33±0.03d 樟子松 0.36±0.01c 8.09±0.01c 1.18±0.02a 25.86±0.44b 1.13±0.05b 桧柏 0.43±0.03b 8.13±0.02ab 1.38±0.22a 29.91±2.94a 1.82±0.15a 1)小写字母表示不同植被恢复地之间的差异(P < 0.05), n=3

表 1 不同植被恢复类型下的土壤理化性质1) Table 1 Soil physiochemical characteristics under different regeneration scenarios

不同植被恢复类型对土壤中碱性蛋白酶、脲酶、脱氢酶和碱性磷酸酶活性影响见图 1.从中可知, 不同土壤酶对植被类型的响应并不一致.

图 1

Fig. 1

不同小写字母表示不同植被恢复地之间的差异(P < 0.05), n=3 图 1 不同植被恢复类型下土壤酶活性 Fig. 1 Soil enzyme activities in different regeneration scenarios

樟子松具有最高的碱性蛋白酶活性, 油松林次之, 桧柏林与侧柏林活性差异较小, 云杉林的碱性蛋白酶活性最低, 各样地之间的碱性蛋白酶活性均具有显著差异[图 1(a)]; 云杉林具有显著高于其他样地的土壤脲酶活性, 其次是桧柏, 云杉林的土壤脲酶活性是桧柏林的2.5倍左右, 侧柏林具有最低的土壤脲酶活性, 桧柏和樟子松、樟子松和油松之间差异不显著[图 1(b)]; 云杉林具有最高的土壤脱氢酶活性, 油松林次之, 云杉林土壤脱氢酶活性显著高于侧柏、樟子松和桧柏等3种样地, 云杉林和油松林、油松林和侧柏林之间差异不显著[图 1(c)]; 樟子松与油松具有相近的碱性磷酸酶活性, 其中樟子松略高于油松, 二者均显著高于其余3种样地, 侧柏、云杉和桧柏林地彼此间碱性磷酸酶活性差异不显著[图 1(d)].

不同植被恢复地土壤微生物功能基因绝对丰度如图 2所示.从中可知, 16S rRNA、6种类型的功能基因丰度和全部功能基因丰度之和均表现为云杉样地基因丰度最高, 樟子松林次之, 其他3种恢复地较低; 除甲烷代谢功能基因外, 桧柏、油松和侧柏这3种样地均无显著差异.与碳代谢有关的功能基因之和远超其余基因, 其中碳固定基因丰度最高, 碳降解基因最低, 甲烷代谢基因最低.

图 2

Fig. 2

不同小写字母表示不同植被恢复地之间的差异(P < 0.05), n=3 图 2 不同植被恢复地土壤微生物功能基因绝对丰度 Fig. 2 Absolute abundance of different soil microbial functional genes under different vegetation restoration

在碳固定相关功能基因中, 丰度最高的是acsA(脱氢酶)和mct(转移酶)基因, 且云杉林这两种基因丰度均远超其他样地[图 2(c)].在丰度最高的碳降解功能基因中, 云杉林和樟子松林具有显著高于其他样地的xylA(木糖异构酶)丰度, 且云杉林abfA(阿拉伯呋喃糖苷酶)基因丰度远高于其他样地; 除樟子松林外, 其他样地均没有检出pox(酚氧化酶)基因[图 2(d)].在甲烷代谢功能基因中, pmoA(甲烷/氨单加氧酶)基因丰度在云杉林和樟子松林最高, 其余3种林地差异不显著[图 2(e)].氮循环相关基因中最主要的是ureC(脲酶)基因, ureC基因丰度在各样地之间的相关关系与总功能基因丰度一致[图 2(f)].磷循环功能基因中丰度最高的是phnK(磷酸转运系统蛋白), 尽管油松的总功能基因丰度高于桧柏, 油松的phnK却略大于桧柏; cphy(半胱氨酸植酸酶)基因仅在油松林中检出[图 2(g)].在与硫循环过程有关的功能基因中, 丰度最高的3种基因soxY(硫氧化蛋白)、yedZ(亚硫酸盐氧化酶)和apsA(腺苷-5′-磷酸还原酶)均以云杉林中的丰度最高, 且云杉林略高于樟子松, 但显著高于其他样地; 除云杉林和樟子松林外, 其余样地未检测出dsrB(亚硫酸盐还原酶)基因[图 2(h)].

为了评价样地中的基因多样性水平, 图 3展示的是各样地的多样性指数.从中可知, 云杉林多样性指数要低于其他4种林地, 除云杉林外的其他4种样地差异不显著[图 3(a)和3(b)].图 3(c)是利用功能基因丰度绘制的非度量多维尺度分析(NMDS)图, 排序分析的平均胁强系数为0.069 7, 表明结果可信.各样地间的功能基因丰度差异显著(P < 0.01), 可见各恢复地土壤微生物群落差异显著.由图 3(d)可知, 云杉和樟子松在聚类分析上距离较近, 油松和桧柏较近.

图 3

Fig. 3

(c)5种样地功能基因非度量多维尺度分析, (d)5种样地功能基因丰度热图; 不同小写字母表示不同植被恢复地之间的差异(P < 0.05), n=3 图 3 不同植被恢复下土壤功能微生物多样性指数 Fig. 3 Soil function microbial index under different vegetation restoration

表 2是对土壤酶活性和土壤微生物碳氮磷硫功能基因丰度间相关性分析, 从中可知, 4种酶类彼此之间并无相关性, 且除土壤脲酶外, 其余3种酶类与微生物碳氮磷硫功能基因丰度均不相关.土壤脲酶与六大类功能基因均有极显著相关性(P < 0.01).

表 2 (Table 2)

表 2 土壤酶活性与微生物功能基因的相关检验1) Table 2 Correlation coefficients between soil enzyme activities and microbial functional genes 蛋白酶 脲酶 脱氢酶 磷酸酶 碳固定 0.295 0.974** 0.691 0.692 碳降解 0.284 0.984** 0.647 0.669 氮循环 0.297 0.964** 0.720 0.703 磷循环 0.300 0.988** 0.636 0.675 硫循环 0.270 0.992** 0.494 0.602 甲烷代谢 0.290 0.995** 0.567 0.641 1)*表示相关水平为0.05(双尾), **表示相关水平为0.01(双尾)

表 2 土壤酶活性与微生物功能基因的相关检验1) Table 2 Correlation coefficients between soil enzyme activities and microbial functional genes

各样地综合肥力指数见表 3.利用主成分分析, 抽取特征向量值大于1的主成分3个, 其特征向量值分别为10.432、5.016和2.164, 这3个主成分的累计贡献率达到97.844%, 第一主成分的贡献率达57.953%, 其中16s rRNA、碳固定和氮循环基因丰度的荷载系数大.第二主成分的贡献率为27.868%, 该主成分在容重和蛋白酶上的荷载最大.第三主成分贡献率为12.022%, 主要受磷酸酶支配.通过各个样地主成分得分与贡献率加权得到综合肥力指数, 从高到低排序依次是; 樟子松>云杉>油松>桧柏>侧柏.

表 3 (Table 3)

表 3 各样地综合肥力指数(IFI) Table 3 Integrated fertility index of each field 植被类型 PC1 PC2 PC3 IFI 侧柏 1.58 0.17 0.10 1.85 油松 1.63 0.23 0.12 1.98 云杉 1.78 0.13 0.10 2.01 樟子松 1.72 0.27 0.12 2.11 桧柏 1.64 0.18 0.10 1.91

表 3 各样地综合肥力指数(IFI) Table 3 Integrated fertility index of each field

土壤微生物分解是重要的养分循环过程之一, 其中微生物功能基因编码参与物质循环的各种酶, 能够反映土壤养分循环活动的强度.土壤微生物量越多, 其基因丰度越高^[14].土壤酶由微生物或植物根系分泌产生, 其活性的高低是土壤生物活性强弱的综合体现^[15].在本试验中涉及的16S rRNA基因和6大类功能基因中, 云杉林均具有最高的绝对丰度, 这表明云杉林的土壤环境可能适宜土壤微生物的生长.同时, 云杉林也具有最高的脱氢酶活性, 脱氢酶是一种具有代表性的胞内酶, 参与土壤微生物的呼吸活动, 可酶促脱氢反应, 在土壤化学反应中作为氢元素的供体, 其在土壤中的活性可在一定程度上表征土壤微生物的多少及相应的活性^{[16, 17]}.

土壤脲酶和蛋白酶是与土壤氮素循环相关的酶^{[18, 19]}.脲酶是存在于多种原核和真核生物中的金属结合酶, 主要参与土壤中氮素的循环, 能促进土壤中含氮有机物如尿素分子酰胺肽键的水解, 生成可被植物根系吸收利用的铵和CO₂, 因此脲酶活性常被用来表征土壤中N素的情况^[20~22].在本试验中, 云杉林具有最高的脲酶活性.有研究认为, 土壤脲酶与土壤肥力有密切的相关关系^[23], 本研究中侧柏样地脲酶活性偏低可能是由于侧柏林具有较低的综合肥力导致的. ureC是编码尿素酶功能亚基的最大的基因, 同时具有若干个高度保守的区域, 因此常被用来作为脲酶活性分析的靶基因^{[24, 25]}.本试验中云杉林具有最高的ureC基因丰度和脲酶活性, 且显著高于其他样地.土壤脲酶与6大类功能基因均具有极强的相关性(表 2), 可见土壤脲酶能更敏感地指示土壤微生物功能基因丰度的变化, 是土壤微生物活性的良好指标.

土壤中的磷元素往往会和金属离子形成沉淀或被土壤吸附固定, 从而形成难溶性磷, 碱性磷酸酶促进土壤中的难溶性磷转化为可溶性磷酸盐, 提高土壤磷的生物有效性, phoD基因参与碱性磷酸酶的编码, 对土壤有效磷的释放有重要作用^[26~28].本研究中油松和樟子松林碱性磷酸酶活性显著高于其他3种样地[图 1(d)], 这表明油松和樟子松的生长需更多磷元素供应.尽管云杉林具有最高的phoD功能基因丰度, 其碱性磷酸酶活性却显著低于油松和樟子松, 这可能是由于phoD基因只有在磷限制条件下才会激活^{[29, 30]}, 而云杉林并不具备激活phoD基因的条件.

云杉林土壤中总硫含量显著低于其他4种样地, 由于各样地的初始条件一致, 这表明云杉的复垦有效减少了土壤中的硫元素.本研究中硫循环有关的功能基因中, 与硫氧化相关的yedZ和soxY基因丰度要远大于与硫还原相关的dsrA和dsrB^{[31, 32]}[图 2(h)], 说明该地区土壤的硫氧化过程要强于硫还原过程.硫氧化菌是一类以硫化物或单质硫为能源的化能自养菌, 具有氧化产酸能力, 可以对无机硫和有机硫能够进行氧化脱除^[33].云杉林可能具有较高的硫氧化菌丰度, 因此对土壤硫的脱除能力较强, 所以土壤总硫较低.

β多样性指数反映群落间微生物物种组成上的差异, 是微生物对环境异质性的响应, 可以利用非度量多维尺度分析(NMDS)来反映^{[34, 35]}.物种组成越相似, 反映在图中就越接近.本研究以5种样地为分类单元, 以各样点的微生物功能基因丰度为依据, 进行排序.结果表明, 各样地间的功能基因丰度差异显著, 较好地体现出了微生物群落的植被差异性[图 3(c)].对功能基因绝对丰度聚类分析得到热图[图 3(d)], 由图可知樟子松和云杉、油松和桧柏在聚类分析上距离很近, 这与NMDS分析的结果是一致的.尽管云杉林具有最高的功能基因丰度, α多样性指数却显著低于其他4种恢复方式的林地[图 3(a)和3(b)].这可能是由于α多样性指数是在综合考虑种数和均匀度的基础上, 反映生物多样性水平导致的^[36].由于云杉林各功能基因间丰度差异较大, mct、ureC和acsA等3种基因丰度之和占全部功能基因的52%以上, 才表现出了绝对丰度高而多样性指数低的现象.

(1) 植被恢复类型对土壤理化、酶活性、功能基因丰度和土壤微细菌丰度均有显著影响, 这些指标最差的均出现在侧柏样地, 樟子松的理化特征最好, 而云杉的细菌和功能基因丰度均显著高于其它样地.

(2) 脲酶与六大类功能基因丰度均具有极显著相关性, 可见脲酶和功能基因丰度能较好指示土壤恢复情况.

(3) 构建了以土壤细菌、碳固定和氮代谢基因为主的土壤综合肥力指标模型, 从土壤综合肥力指标来看, 樟子松和云杉林适宜于该地区复垦树种, 而侧柏并不适应于该地区.