连续施用化肥及秸秆还田对潮土酶活性、细菌群落和分子生态网络的影响

马 垒，李 燕，魏建林，周晓琳，李子双，李国生，吴小宾，刘兆辉，谭德水*

（1 山东省农业科学院农业资源与环境研究所 /农业农村部废弃物基质化利用重点实验室，山东济南 250100；2 德州市农业科学研究院，山东德州 253015）

潮土是黄淮海平原最主要的土壤类型，其砂粒含量高、土壤结构差，在常年麦玉两熟的高强度利用下，土壤肥力水平普遍低下[1]。为了满足粮食需求，大量化肥被施入土壤，为作物生长补充大量养分，迅速提升作物产量。然而不合理化肥施用导致土壤出现板结、酸化，还会引发农业面源污染等环境问题[2]。秸秆还田配施化肥是一种可持续农业生产措施，可有效缓解单施化肥引起的土壤质量下降问题[3]。越来越多的研究表明，秸秆与化肥配合施用可显著提升土壤有机质和养分含量，提高土壤酶活性，改良土壤结构，增加作物产量，同时降低秸秆焚烧引起的空气污染[4-7]。但化肥及秸秆还田对土壤微生物群落的影响仍有待明确。

土壤微生物是农田生态系统中最为活跃的组分，其在有机质形成与降解、养分循环转化、土壤团聚体结构形成中具有重要作用[8]。由于微生物对其生存环境变化极其敏感，可对细小的土壤环境变化迅速做出反应，因此研究微生物群落的变化可较为直观地反映土壤生态功能的变化[9]。关于化肥和秸秆还田对微生物群落影响的研究已有大量报道，大多数研究发现，长期施用化肥引起土壤酸化，导致微生物多样性降低[10-11]；而秸秆还田配施化肥可缓解酸化问题，从而提高土壤微生物多样性，并促进土壤中参与碳、氮、磷循环的特定微生物生长[12-14]。然而一些研究者发现长期秸秆还田后土壤微生物多样性并未发生明显改变[5,15]。我国安徽砂姜黑土长期定位试验结果表明，连续30年秸秆还田配施化肥后细菌多样性显著降低[3,16]。这种不一样的结论，一方面是由于试验地土壤养分状况和自然条件不同，另一方面可能与秸秆还田量、还田方式有关。此外，大多数研究的焦点在于微生物多样性或群落结构，忽视了对种间关系的研究。土壤中的微生物并非单独存在，不同物种通过频繁的相互作用，形成复杂的生态网络，共享土壤生态位和营养资源。微生物参与的生态功能不仅仅是不同个体的集体特征，还与它们之间的相互作用密切相关。研究物种间的网络互作关系对于理解微生物群落构建过程，了解物种间的生态位差异或环境异质性，预测土壤生态功能具有重要意义[17]。因此从网络共现模式角度解析长期化肥及秸秆还田对微生物群落影响的机制仍需探讨。

以位于山东德州的潮土长期施肥定位试验为研究平台，采用16S rRNA高通量测序技术，对连续10年不施肥、单施化肥、秸秆还田配施化肥后土壤细菌群落多样性、结构、组成展开研究；同时通过构建分子生态网络，研究不同施肥模式对细菌种间互作关系和网络共现模式的影响；探索细菌群落与土壤理化性质、酶活性之间的关系，以期为黄淮海潮土区土壤肥力提升和秸秆资源化利用提供理论依据。

1 材料与方法

1.1 长期试验地点概况

田间定位试验始于2010年10月，位于山东省德州市哨马营村 (116°17'44'' E，37°30'34'' N)。该地属暖温带半湿润季风气候，年平均气温12.9℃，年降水量547.5 mm，供试土壤为非石灰性潮土。试验地为冬小麦-夏玉米轮作，小麦品种为济麦22，玉米品种为郑单958。试验开始前土壤基本理化性质如下：pH 8.3、有机质20.3 mg/kg、全氮1.2 g/kg、全磷0.83 g/kg、全钾20.8 g/kg、有效磷37.4 mg/kg、速效钾282 mg/kg。

1.2 土壤样品采集与测定

2021年在该长期定位试验选择不施肥料对照(CK)、单施化肥(NPK)和秸秆还田配施化肥(NPKS)3个处理进行取样分析，每个处理4个重复，共计12个小区。NPK和NPKS处理小麦季施N 300 kg/hm2、P2O5120 kg/hm2、K2O 100 kg/hm2，玉米季施N 250 kg/hm2、P2O545 kg/hm2、K2O 45 kg/hm2，NPKS 处理小麦和玉米秸秆粉碎后全量还田，于下茬作物播种前旋耕翻入土壤。供试化肥为尿素(N 46%)、过磷酸钙 (P2O516%)和硫酸钾 (K2O 50%)。

土壤样品于2021年小麦灌浆期采集。截止样品采集时，共连续种植了11季小麦、10季玉米。按照S形采样法，使用直径5 cm的土钻采集0—20 cm土壤，每个小区采取7钻混合成1个样品重复。所有样品置于保温箱，采用冰袋保鲜运回实验室。土壤样品在去除大的石块和植物根系后过2 mm筛。所有样品分为3部分，一部分常温风干后用于土壤理化性质测定，一部分4℃保鲜用于土壤酶活性测定，另一部分-80℃保存用于土壤DNA提取。

土壤理化指标测定均采用常规方法[18]。pH采用水土比(体积质量比) 5∶1提取，pH计测定；有机质(SOM)采用重铬酸钾容量法测定；全氮(TN)采用浓硫酸消煮，全自动凯氏定氮仪测定；全磷(TP)采用浓硫酸-高氯酸消煮，钼锑抗比色法测定；碱解氮(AN)采用碱解扩散法测定；有效磷(AP)采用0.5 mol/L碳酸氢钠浸提—钼锑抗比色法测定；速效钾(AK)采用醋酸铵浸提—火焰光度计法测定。土壤酶活性采用北京索莱宝科技有限公司试剂盒测定：β-葡糖苷酶采用对硝基酚比色法测定，以单位时间单位质量土样催化对-硝基苯-β-D 吡喃葡萄糖苷生成的对-硝基苯酚表示；脲酶采用靛酚蓝比色法，以单位时间单位质量土样水解尿素产生的铵态氮量表示；碱性磷酸酶采用磷酸苯二钠比色法，以单位时间单位质量土样催化磷酸苯二钠水解生成的苯酚量表示；脱氢酶采用氯化三苯基四氮唑还原法，以单位时间单位质量土样还原氯化三苯基四氮唑(TTC)生成的三苯基甲臜量表示。土壤生化指标结果均以烘干土(105℃，24 h)为基准表达。

采用Fast DNA Spin Kit for Soil (MP Biomedicals，Santa Ana，CA，USA)试剂盒提取0.50 g鲜土DNA，具体步骤参见说明书。采用16S rRNA V4～V5区进行细菌高通量测序，使用515F/907R通用引物进行PCR扩增。每个样品前端引物均含有不同的7 bp Barcode用于区分不同样品。PCR扩增条件包括94℃5 min，90℃ 60 s，55℃ 60 s，72℃ 75 s，30 个循环，之后 72℃ 10 min。反应产物采用QIA quick PCR Purification kit (Qiagen)进行纯化。将不同样品的PCR扩增产物等摩尔混合后，采用Illumina公司MiSeq测序仪完成序列分析。

1.3 数据统计分析

高通量测序数据处理主要使用Usearch软件[19]完成：1)合并双端测试序列；2)切除引物并去除质量分数低于20的低质量序列；3)采用RDP数据库(http://rdp.cme.msu.edu/)去除嵌合体；4)使用UNOISE算法去噪生成扩增子序列变体(amplicon sequence variant，ASV)；5)与RDP数据库进行比对进行注释。经过质量控制后，总共得到673598条高质量序列(每个样品48929～68929条序列)，所有样品随机抽取48929条序列进行后续分析。测序原始数据已上传至NCBI的SRA数据库(SRR17968322～SRR 17968333)，项目编号为PRJNA805075。

细菌alpha多样性(包括Richness丰富度指数、Shannon多样性指数和Peilou’s evenness均匀度指数)、基于Bray-curtis距离的主坐标分析(PCoA)、冗余分析(RDA)和置换多元方差分析(Adonis)均使用R 4.0.5的“vegan”包进行计算。多元回归树分析(MRT)采用“mvpart”包进行计算。不同处理ASV水平指示种分析采用“indicspecies”包进行计算，属水平差异物种在Galaxy网站(http://huttenhower.sph.harvard.edu/galaxy/)采用线性判别分析(LEfSe)进行计算，其中线性判别分析(LDA)值大于2 的物种为显著差异物种。采用“WGCNA”包计算ASV间Spearman相关系数，为了减少网络复杂度，只选取|r|＞0.9，P＜0.001的关系用于网络构建；采用Fast Greedy法进行网络模块分析，使用Cytoscape 3.8.2软件进行网络图形绘制。单因素方差分析和Pearson相关性双尾检验采用SPSS 25.0软件完成，采用Fisher LSD法进行平均值多重比较。采用Microsoft Excel 2019 和Origin 2021b软件进行数据处理和绘图。

2 结果与分析

2.1 土壤理化性质

由表1可知，经过连续10年不同施肥处理后，土壤理化性质发生明显改变。与CK处理相比，NPK和NPKS处理的土壤pH显著降低了0.06个单位(P＜0.05)；NPK处理土壤有机质和碱解氮含量分别提升23.19%和24.22%，NPKS处理则分别提升34.82%和39.17%，较NPK处理分别提升9.20%和12.03%；与CK处理相比，NPK处理全氮、全磷、有效磷和速效钾分别升高了26.80%、27.06%、298.39%和26.32%，NPKS处理分别升高了30.93%、17.65%、310.43%和25.74%，但NPK和NPKS处理间差异未达显著水平(P＞0.05)。

表1 长期不同施肥处理下土壤理化性质Table 1 Physio-chemical properties of soil under different long-term fertilizations

2.2 土壤酶活性

通过测定β-葡糖苷酶(β-GC，降解纤维素)、脲酶(UE，水解尿素)、碱性磷酸酶(ALP，水解有机磷)和脱氢酶(DHA，表征微生物活性)等4种土壤胞外酶活性，表征长期施用化肥及秸秆还田对土壤微生物功能和活性的影响(图1)。β-GC活性高低顺序为NPKS＞NPK＞CK，NPKS处理β-GC活性较CK和NPK分别增加了78.31%和17.09%，NPK处理较CK处理升高47.91%。与CK处理相比，NPK和NPKS处理ALP活性分别显著增加了50.35%和46.53%；DHA在NPKS处理中最高，比CK和NPK处理分别提高50.91%和37.21%。UE活性3个处理间差异不显著(P＞0.05)。

图1 长期不同施肥处理下土壤酶活性Fig.1 Soil enzyme activities under different long-term fertilizations

2.3 土壤细菌Alpha 多样性

通过16S rRNA高通量测序总共得到673598条高质量序列，可划分为5993个ASV (每个样品4860～5088条ASV)，经过与RDP数据库进行比对，这些ASV归属于23个门，50个纲，69个目，145个科和323个属。将所有样品均一化为48929条序列后，不同施肥处理细菌Alpha多样性如图2a所示。细菌多样性(Diversity)和均匀度(Evenness)在CK处理中最高(11.11和0.90)，而在NPK处理(11.02和0.89)和NPKS处理(11.01和0.89)中均显著降低；丰富度(Richness)处理间差异未达显著水平。

相关性分析结果(图2b)表明，细菌多样性与土壤有机质 (r=0.76，P＜0.01)、全氮 (r=0.59，P＜0.05)、碱解氮 (r=0.59，P＜0.05)、有效磷 (r=0.81，P＜0.01)、速效钾 (r=0.60，P＜0.05)、β-葡糖苷酶(r=0.62，P＜0.05)、碱性磷酸酶 (r=0.60，P＜0.05)等生化指标呈显著负相关，与pH (r=0.74，P＜0.05)呈显著正相关；细菌均匀度与土壤有机质 (r=0.81，P＜0.01)、全氮 (r=0.59，P＜0.05)、碱解氮 (r=0.62，P＜0.05)、有效磷 (r=0.81，P＜0.01)、速效钾 (r=0.62，P＜0.05)、β-葡糖苷酶 (r=0.68，P＜0.05)、碱性磷酸酶(r=0.59，P＜0.05)等生化指标呈显著负相关，与pH(r=0.75，P＜0.05)呈显著正相关；细菌丰富度则与土壤生化指标相关性未达显著水平。以上结果说明长期施用化肥和秸秆还田后，土壤养分含量的升高反而导致细菌Alpha多样性显著降低。

图2 长期不同施肥处理下土壤细菌Alpha 多样性Fig.2 Soil bacteria alpha-diversity under different long-term fertilizations

2.4 土壤细菌群落结构和组成

长期施用化肥及秸秆还田下土壤细菌群落结构和组成发生明显改变。变形菌门(Proteobacteria)是所测潮土中的优势菌门(图3a)，其相对丰度在12个样品中介于40.40%～49.20%，其余丰度大于1%的物种包括酸杆菌门(Acidobacteria，相对丰度16.10%～21.30%)、放线菌门(Actinobacteria，相对丰度8.01%～11.80%)、浮霉菌门(Planctomycetes，相对丰度3.54%～5.20%)、拟杆菌门(Bacteroidetes，相对丰度3.07%～7.01%)、绿弯菌门(Chloroflexi，相对丰度1.96%～3.55%)、WPS-1 (相对丰度0.89%～1.41%)和厚壁菌门(Firmicutes，相对丰度0.67%～1.24%)。基于Braycurtis距离对细菌群落结构进行主坐标分析(PCoA)(图3b)，结果发现主坐标的前两轴总共解释了34.06%的群落差异(第一轴和第二轴分别解释了21.31%和12.75%的群落差异)：CK处理聚集在主坐标第一轴的左侧，NPK处理聚集在主坐标第一轴的右侧，而NPKS处理介于两者之间。置换多元方差分析(Adonis)结果表明，NPKS处理与CK (R2=0.27，P=0.03)和NPK处理(R2=0.20，P=0.03)细菌群落结构均具有显著差异。

图3 长期不同施肥处理下土壤的细菌群落结构Fig.3 Soil bacteria community structures under different long-term fertilizations

采用多元回归树分析(MRT)和冗余分析(RDA)检验细菌群落结构与环境因子间的关系(图3c, d)。结果发现MRT模型总共解释了35.31%的群落变异，基本将3个处理分开(除1个NPKS处理归属于第二片叶子)。细菌群落首先由土壤有机质(15.80 g/kg)分为不施肥处理(CK)和施用化肥处理(NPK和NPKS)两大分支，然后由土壤有机质 (18.18 g/kg)又分为单施化肥(NPK)和秸秆还田配施化肥(NPKS)两大分支。与MRT结果相似，RDA分析(图3d)结果同样发现土壤有机质是导致NPK和NPKS处理细菌群落变异的最主要因子。以上结果表明长期施用化肥及秸秆可明显改变细菌群落结构，而这主要是由于土壤有机质含量的升高造成的。

在ASV水平上对细菌群落进行指示种分析(图4a和图4b)，结果发现在CK、NPK和NPKS处理中分别有141、49和38个指示种，其中CK处理指示种中酸杆菌门、浮霉菌门相对丰度在3个处理中最高；NPK处理中变形菌门、放线菌门和拟杆菌门相对丰度在3个处理中最高，而NPKS处理中厚壁菌门相对丰度在3个处理中最高。属水平上的线性判别分析(LEfSe)发现(图4c)，CK处理中相对丰度显著升高的物种主要归属于酸杆菌门的Acidobacteria_Gp3纲，归属于变形菌门的生丝微菌科(Hyphomicrobiaceae)和伯克氏菌科(Burkholderiaceae)；NPK处理中相对丰度显著升高的物种主要归属于变形菌门的草酸杆菌科(Oxalobacteraceae)、红丹杆菌科(Rhodanobacteraceae)、黄单胞菌科(Xanthomonadaceae)；NPKS处理中显著升高的物种主要归属于厚壁菌门的芽胞杆菌科(Bacillaceae)。

图4 长期不同施肥处理下土壤细菌的群落组成Fig.4 Soil bacteria community composition under different long-term fertilizations

2.5 土壤细菌分子生态网络

通过构建细菌分子生态网络，研究不同施肥模式下土壤细菌物种间互作关系，结果如图5a所示。图中每一个节点代表一个ASV，各处理指示种采用不同颜色节点标识，两个节点间的连线表示存在显著相关关系。通过对网络进行模块(相同模块内物种高度相连，与其它模块物种关系较少)分析，发现整个网络由4个主要模块组成：CK处理中的指示种主要分布在模块1和模块2中；NPK处理中的指示种主要聚集在模块3；而NPKS处理中的指示种集中于模块4。对4个模块内的物种累计丰度进行统计(图5b)，结果发现模块1和模块2内物种累计丰度在CK处理中最高，较NPK处理分别升高136.61%和69.70%，较NPKS处理分别升高76.95%和103.74%；模块3内物种累计丰度在NPK和NPKS处理中分别是CK处理的3.57和2.57倍；模块4内物种累计丰度按照CK＜NPK＜NPKS的顺序升高，其中NPKS处理内物种累计丰度分别较NPK和CK处理提升52.69%和269.57%。

相关性分析结果(图5c)表明，模块1内物种累计丰度与土壤有机质、全氮、全磷、碱解氮、有效磷、速效钾、β-葡糖苷酶、碱性磷酸酶等生化指标显著负相关，与pH呈显著正相关；模块2内物种累计丰度与土壤有机质、全氮、碱解氮、有效磷、速效钾、β-葡糖苷酶、碱性磷酸酶和脱氢酶等生化指标显著负相关；模块3内物种累计丰度与土壤有机质、全氮、全磷、速效磷、速效钾和碱性磷酸酶等指标呈显著正相关，与pH呈显著负相关；模块4内物种累计丰度与土壤有机质、全氮、碱解氮、速效磷、速效钾、β-葡糖苷酶、碱性磷酸酶和脱氢酶等指标呈显著正相关，与pH呈显著负相关。以上结果说明不同施肥模式可通过改变土壤养分含量，调节细菌物种网络互作关系，促进特定微生物集群的生长，从而改变土壤碳、磷循环转化胞外酶活性。

图5 长期不同施肥处理下土壤细菌的分子生态网络Fig.5 Soil bacterial co-occurrence patterns under different long-term fertilizations

3 讨论

3.1 连续施用化肥及秸秆还田对土壤细菌群落的影响

微生物多样性是土壤质量和肥力评价的重要依据[3]。然而，在本研究中，长期施用化肥(NPK)和化肥配合秸秆处理(NPKS)虽然提高了土壤有机质含量和全量及速效氮磷钾含量，但细菌多样性分别降低0.90%和0.91%，均匀度均降低1.11%。根据微生物对碳底物利用能力和生存策略的不同，可分为富营养型物种(R-策略，Copiotrophs)和寡营养型物种(K-策略，oligotrophs)。富营养型物种适宜在养分含量和有效性高的土壤中生长，优先利用简单有机碳，表现出相对较快的生长速度；相反地，寡营养型物种则偏向于养分含量和有效性低的土壤，可分解利用顽固有机碳，具有较慢的生长速率[20-21]。因此施用化肥后土壤有机质和有效磷等养分的升高可能会淘汰环境中的K-策略物种，导致细菌多样性降低。Sun等[16]同样发现砂姜黑土中土壤肥力与细菌群落变化呈相反趋势，并认为微生物多样性的变化与土壤肥力的变化并不一致。我们的研究发现细菌群落多样性与多个土壤养分指标呈显著负相关关系，表明土壤中养分的富集是细菌多样性降低的重要因素。由于大多数微生物在土壤中占据特有的生态位并执行不同生态功能，因此微生物多样性的降低通常意味着功能的损失[22]。然而本研究发现细菌多样性与土壤β-葡糖苷酶和碱性磷酸酶等胞外酶呈显著负相关关系，说明即便是在生物多样性较低(相对于自然生态系统)的农田生态系统中，土壤微生物仍存在较强功能冗余现象。土壤生态系统中每种代谢功能，可由多个分类地位不同但可以共存的物种执行，因此物种多样性的降低并不会影响微生物群落执行特定的生态功能。

长期不同施肥处理导致细菌群落组成发生明显改变，CK处理土壤富集的群落大多归属于酸杆菌门(Acidobacteria_Gp3纲)、浮霉菌门、生丝微菌科和伯克氏菌科。其中酸杆菌门是典型的寡营养型物种，广泛分布在养分贫瘠和退化土壤中，其相对丰度与土壤有机质含量呈显著负相关[23]。浮霉菌门则被认为是一种温和的寡营养型物种(moderate oligotrophs)，对不同的环境具有广泛的适应性[24]，且在秸秆还田处理后丰度显著降低[13]。一个有趣的现象是归属于变形菌门(通常被认为是富营养型物种)的生丝微菌科和伯克氏菌科同样在CK处理中广泛富集，研究表明生丝微菌科和伯克氏菌科均具有较强的固氮功能，可吸收利用大气中N2，这说明固氮微生物可能在全氮和碱解氮含量较低的CK处理中具有较强的竞争力[25]。NPK处理中变形菌门(草酸杆菌科、红丹杆菌科和黄单胞菌科)、放线菌门和拟杆菌门相对丰度显著升高。其中草酸杆菌科是重要的几丁质溶解细菌[26]；红丹杆菌科和黄单胞菌科均属于黄单胞菌目，可以分解木质素作为碳源生长[27]；而放线菌门和拟杆菌门均具有较强的分解土壤中复杂大分子有机质的能力[28]。与前人研究结果[29]相同，长期秸秆还田处理中富集大量厚壁菌门(芽胞杆菌科)，芽胞杆菌科中的许多物种可以分泌纤维素酶，参与秸秆中纤维素的降解。由以上分析可以看出土壤有机质主导了潮土细菌群落结构和组成，这与位于河南郑州[12]和封丘[30]的潮土施肥定位试验结果相同。长期不施肥处理中富集了大量养分需求低的寡营养型细菌和固氮菌；单施化肥处理中由于C∶N较低，为了获得碳源，一些利用顽固碳组分的细菌类群大量富集；而常年秸秆还田配施化肥为土壤提供充足易利用碳源，激发了土壤中具有纤维素降解功能的芽胞杆菌科相关物种。

3.2 连续施用化肥及秸秆还田对土壤细菌分子生态网络的影响

Li等[12]通过网络分析发现，长期不同施肥后土壤中富集的物种间存在强烈的合作关系，且与土壤有机质含量呈显著相关关系，并认为有机无机肥配施可通过调节微生物种间关系，影响土壤有机碳累积。Fan等[31]发现长期施肥抑制了固氮微生物集群(模块)的生长，导致土壤固氮功能的降低，说明施肥对土壤生态功能的调控可能是通过影响特定微生物模块的生长。本研究中，我们将3种施肥方式下所有细菌物种进行网络构建，并对不同施肥模式下的指示物种在网络中的分布模式进行研究。发现施肥方式不同，网络中的指示种通过强烈的共生关系聚集在不同的模块中。模块1和模块2主要包含CK处理的优势菌群，且均与土壤养分指标、β-葡糖苷酶和碱性磷酸酶呈显著负相关，说明长期施用化肥及秸秆还田降低了这两个模块内的群落丰度，淘汰了那些占据低养分含量生态位，且碳、磷降解能力差的群落。模块3主要包含NPK施肥模式的指示种，与多个土壤养分指标正相关，但与β-葡糖苷酶相关性不高，说明该模块内的物种可能降解纤维素的能力不强。模块4内聚集了大量NPKS处理的指示种，且与土壤有机质、全氮、碱解氮、有效磷、速效钾、β-葡糖苷酶、碱性磷酸酶和脱氢酶等指标呈显著正相关，说明长期秸秆还田配施化肥刺激该模块内物种生长，这些微生物通过强烈的种间共生关系，形成碳、磷循环功能微生物集群。前人研究结果表明，秸秆还田后土壤微生物丰度的变化导致群落生态功能发生改变[12,32]，而本研究结果则进一步说明，长期不同施肥模式下土壤生态功能的改变，不仅与物种个体丰度变化有关，还与微生物种间共生关系密切相关。秸秆还田配施化肥可促进特定网络模块内物种的生长，形成强烈的共生关系，从而提高细菌功能和活性。

4 结论

化肥配合秸秆还田不仅提升了土壤肥力，还改变了细菌群落结构和组成，虽然细菌的多样性较不施肥有显著降低，但激活了土壤中芽胞杆菌等纤维素降解群落的活性，通过强烈共生关系形成碳、磷循环功能微生物集群，优化了细菌群落结构和组成，提升了土壤碳、磷的循环转化功能。