长期有机/无机肥配施对褐土细菌群落及抗生素抗性基因赋存的影响

刘志平，解文艳，杨振兴，胡雪纯，马晓楠，周怀平

(山西农业大学资源环境学院，省部共建有机旱作农业国家重点实验室(筹)，山西 太原 030001)

有机肥与无机肥配合施用是目前农业生产中所推崇的一种施肥方法，既可以改良培肥土壤又可以提高作物产量，这是因为有机肥可将大量的有益微生物带入土壤。作为土壤的重要组成部分，微生物是土壤有机质和养分循环的主要驱动力，在土壤养分转化、能量流动和土壤有机质形成中具有关键作用[1-4]。细菌常因对土壤环境变化反应迅速而被作为土壤质量的早期生物指标[5-6]。细菌多样性和群落组成受土壤pH、有机质、土壤水分、矿物质等因素影响[7-8]，而这些因素又受施肥管理措施的影响。因此，研究细菌群落及多样性对于指导合理施肥具有重大意义。

随着畜禽养殖的集约化，作为亚治疗剂量的食品添加剂，抗生素被广泛应用于畜禽养殖业以防治动物疾病、促进动物生长及提高饲料利用率。我国是世界上最大的抗生素生产国和消费国之一，统计表明，2013年中国的抗生素生产量为2.48×105t，使用量达1.62×105t，居世界第一，约占世界总量的一半，其中约7.8×104t为人类使用，占比48%，其余8.4×104t用于畜禽养殖[9]。大约30%～90%兽用抗生素通过粪便或尿液以母体化合物或代谢产物的形式排出[10-12]。未经充分腐熟的畜禽粪污还田又将抗生素带入土壤，虽然抗生素可以通过生物降解和光解等作用进行转化，但低浓度抗生素残留也会对微生物群落产生选择性压力，诱导产生抗生素抗性基因(Antibiotic Resistence Genes，ARGs)，一些ARGs甚至能够在选择性压力消失后持续存在，使得土壤成为丰富的ARGs储存库[13-15]。土壤中赋存的ARGs可以通过食物链的传递对人体健康造成潜在危害。与其他污染物相比，ARGs分布范围更广、传播速度更快、环境风险更高，是人类面临的最严峻的生态环境问题[16]，已被世界卫生组织列为21世纪威胁公共健康的最重大挑战之一[17]。随着检测技术的发展，越来越多的ARGs在污水处理厂、动物粪便、土壤、河流、饮用水、冰川等不同环境样品中被检测出来[18-22]。

研究土壤细菌及ARGs的方法众多，其中宏基因组技术作为一种新兴技术可直接从样品中提取总DNA构建宏基因组文库进行测序，进而获得样品中所包含的全部微生物及其遗传基因信息[23]，有助于全面了解不同环境中细菌群落及ARGs的多样性与丰度变化。在耕地中，检测长期施用有机肥对土壤中ARGs赋存特征的影响具有重大意义，可以从食物链源头评估施用有机肥对生态安全及食品安全的潜在风险，但目前此类研究尚较为缺乏。本研究以山西寿阳长达30年的有机/无机肥配施试验为基础，基于宏基因组测序技术分析耕层土壤ARGs及细菌多样性，旨在研究不同有机/无机肥配施对褐土中细菌群落及ARGs赋存特征的影响，为今后农业生产中有机肥的投入量提供理论依据。

1 材料与方法

1.1 研究区概况

本研究中长期定位试验位于山西省晋中市寿阳县宗艾镇宗艾村国家农业环境寿阳观测实验站(37°58′23″N，113°06′38″E)，海拔1 130 m，土壤类型为褐土，系统分类为简育干润雏形土[24]。该区域属中纬度暖温带半湿润偏旱区大陆性季风气候区，年均气温7.4℃，年均降水量500 mm。试验地初始pH为8.4，有机质23.8 g·kg-1；碱解氮117.69 mg·kg-1；有效磷4.84 mg·kg-1；有效钾100 mg·kg-1[25]。

该长期定位试验从1992年开始种植春玉米直至2021年，历时30 a。试验采用氮(N)、磷(P)、有机肥(M)三因素四水平正交设计，另设有对照和高量有机肥区，共18个处理，各处理3次重复，小区面积66.7 m2，随机排列。玉米品种为‘大丰30’，播种密度为66 000株·hm-2，玉米播种前进行旋地和施肥，方式为全部播前基施。本研究依据腐熟牛粪用量梯度，选择其中5个处理，即N0P0M0，N2P1M1，N4P2M2，N3P2M3和N0P0M6，研究长期施用牛粪对土壤ARGs及细菌多样性的影响，具体施肥量见表1。

表1 不同试验处理施肥量Table 1 Amount of fertilizer applied in different experimental treatments

1.2 样品采集

于2021年6月在各小区采用5点“S”取样法，用土钻采集0～20 cm的土壤样品，将同一小区的5个样混合为一个样品，去掉明显的石砾和残枝败叶等杂物后将土壤分成两部分，一部分风干用于土壤理化性质的检测，另一部分用无菌袋装好，并置于冰盒带回实验室，进行土壤DNA的提取及后续宏基因组测序。

1.3 测定方法

1.3.1 土壤理化性质测定 土壤理化性质的测定方法按照《土壤农业化学分析方法》[26]进行。其中，pH采用酸度计检测，水土比为2.5∶1；碱解氮(AN)采用比色法检测；有效钾(AK)采用1 mol·L-1乙酸铵火焰光度计检测；有效磷(AP)采用NaHCO3(pH 8.5、0.05 mol·L-1)浸提法测定；有机质(OM)采用重铬酸钾容量法检测；全氮(TN)采用凯氏定氮仪检测；全磷(TP)采用钼锑抗比色法检测；全钾(TK)采用火焰光度法检测。

1.3.2 土壤宏基因组测序 准确称取0.5 g土样，采用Fast DNA SPIN DNA提取试剂盒(MP Biomedicals,Santa Ana,CA,United States)，按照说明书的方法提取土壤DNA，并用0.7%的琼脂糖凝胶进行电泳，观察条带是否清晰、完整，而后用微量紫外分光光度计NanoDrop ND-1000 (Thermo Fisher Scientific,Waltham,MA,USA) 测定DNA的浓度。将合格DNA送至生工生物工程(上海)股份有限公司进行宏基因组测序。

1.4 数据统计

对测序原始数据通过FastQC进行质量评估，并通过Trimmomatic进行过滤处理[27]，得到有效数据。使用IDBA_UD对各样本有效序列进行拼接组装成长序列contigs，并对多个Kmer的组装结果进行综合评定，选择最佳Kmer组装结果。采用Prodigal对拼接结果进行ORF预测，选择长度≥100 bp的基因并将其翻译成氨基酸序列，对于各样本的基因预测结果，采用CD-HIT软件进行去冗余，以获得非冗余的基因集[28]。使用DIAMOND将基因集与ARDB(Antibotic Resistance Genes Database)数据库进行比对，获得基因的物种注释信息和功能注释信息，并根据基因集丰度得到功能丰度和物种丰度[29]。利用QIIME计算Chao1、ACE、Simpson和Shannon指数。利用MEGAN[30]进行细菌群落分类组成和丰度可视化。采用R软件绘制不同处理ARGs维恩图，采用主成分分析(Principal component analysis，PCA)法研究不同样品间ARGs的变化。采用置换多元方差分析(PERMANOVA)[31]评估各组间ARGs组成差异的显著性，通过CANOCO软件进行门水平细菌与ARGs的冗余分析(Redundancy analysis，RDA)。

采用SPSS 19.0软件进行单因素方差分析，检验土壤理化性质、细菌α-多样性指数等的正态分布和方差同质性。采用LSD检验分析数据间的显著性差异(P<0.05)，并进行细菌门及ARGs与土壤理化性质间的皮尔森相关性分析。

2 结果与分析

2.1 土壤理化性质分析

不同有机/无机肥配施对土壤理化性质的影响如表2所示，与对照处理N0P0M0相比，施肥显著降低了土壤pH值，且随总施肥量增加呈现显著降低趋势。与此相反，随施入养分量的增加，土壤OM、TN、TP、AN和AP的含量均呈增加趋势(P<0.05)。

2.2 土壤细菌多样性分析

2.2.1 细菌α-多样性α多样性指数反映了细菌物种的丰富度及均匀度，是鉴定土壤肥力的有效指标之一。由表3可见，不同有机/无机肥配施处理下，细菌Chao1指数在N2P1M1处理和N0P0M6处理中显著高于N4P2M2处理(P<0.05)，ACE指数在N2P1M1处理中显著高于N4P2M2处理(P<0.05)，Simpson指数和Shannon指数在4个处理中无显著差异。

表2 各处理土壤理化性质Table 2 Soil physicochemical properties in each treatment

表3 不同处理细菌α多样性Table 3 Bacterial α-diversity under different treatments

2.2.2 土壤门水平细菌丰度 由图1可见(见 172页)，不同处理土壤样品中门水平细菌种类相似，相对丰度前五位的分别为变形菌门(Proteobacteria，48.6%～59.1%)、酸杆菌门(Acidobacteria，7.19%～16.4%)、放线菌门(Actinobacteria，8.61%～15.1%)、芽单胞菌门(Gemmatimonadetes，5.19%～6.97%)和绿弯菌门(Chloroflexi，2.53%～4.21%)，平均相对含量占总量的82.6%～85.9%。变形菌门、放线菌门和绿弯菌门的相对丰度随着养分投入量的增加呈先升高后降低的趋势，N3P2M3处理变形菌门相对丰度最高且显著高于其他处理(P<0.05)；N4P2M2处理放线菌门和绿弯菌门相对丰度最高，且显著高于其他处理(P<0.05)。酸杆菌门的相对丰度随牛粪施用量的增加呈下降趋势，且N0P0M0处理显著高于N2P1M1处理，但在N4P2M2、N3P2M3和N0P0M6处理中无显著差异。芽单胞菌门的相对丰度随牛粪施用量增加呈升高趋势，且N0P0M6处理显著高于N0P0M0处理。

图1 不同处理门水平细菌相对丰度Fig.1 Relative abundance of bacteria at phylum level under different treatments

2.3 土壤抗生素抗性基因分析

2.3.1 抗生素抗性基因Venn图 由图2可见(见 172页)，5个处理土壤样品共有256个共同ARGs亚型，N0P0M0、N2P1M1、N4P2M2、N3P2M3和N0P0M6处理含有的独特ARGs亚型个数分别为0，0，1，5和3，说明高量施用牛粪可增加土壤独特ARGs亚型的数量，且随着养分投入量的增加呈先增加后降低趋势，并在N3P2M3处理达到峰值。

图2 不同处理抗生素抗性基因Venn图Fig.2 Venn diagram of antibiotic resistance for different treatments

2.3.2 抗生素抗性基因组成 15个土壤样品中检测到的ARGs亚型共有318个，共分为12个种类。其中亚型较多的种类有大环内酯-林可霉素-链阳性霉素B类(48个)、氨基糖苷类(39个)、四环素类(34个)、万古霉素类(34个)、多药类(31个)。图3展示了相对丰度前30位的ARGs亚型，平均相对丰度>3%的抗性基因亚型有13种，按照含量高低依次为MacB、BcrA、MexW、MexF、pbp1a、CarA、VanSD、CeoB、VanRE、VanRC、VanRA、tetPB和otrA，其中万古霉素类抗性基因居多，MacB基因和BcrA基因在N3P2M3处理中含量最高、在N0P0M0处理中含量最低。

图3 ARGs相对丰度Fig.3 Relative abundance of ARGs

2.3.3 主成分分析及置换多元方差分析 将注释到的ARGs进行降维处理，提取主成分进行分析。由图4可知，三个主成分共解释了62.71%(PC1：31.97%，PC2：15.86%，PC3：14.88%)的ARGs组成及变化。每个处理的三个重复聚拢在一起，说明样品平行性良好。置换多元方差分析进一步显示(图5)，不同施肥处理中土壤ARGs差异显著，检验结果可信度高(R=0.9926，P=0.001)。

图4 ARGs主成分分析Fig.4 Principal component analysis of ARGs

图5 置换多元方差分析Fig.5 Permanova

2.3.4 抗生素抗性基因与门水平细菌的冗余分析 为了评价ARGs与细菌群落结构之间的关系，选择相对丰度前十位的ARGs和门水平细菌进行冗余分析。由图6(见 172页)可知，细菌群落可以解释89.8%(RDA1：80.7%，RDA2：9.1%)的ARGs变化。变形菌门与MacB、CarA、BcrA、VanRC、VanRA等之间的夹角均为锐角，表明它们之间呈正相关关系；酸杆菌门(Acidobacteria)与tetPB、MexF、pbp1a等呈正相关关系，芽单胞菌门(Gemmatimonadetes)与MexW、CarA、MacB、BcrA等呈正相关关系，表明这些细菌门中的某些物种可能是相应ARGs的重要宿主菌。综上所述，ARGs的变化受细菌群落结构的影响。

2.3.5 门水平细菌及ARGs与土壤理化性质的相关性分析 为了进一步揭示细菌及ARGs与土壤理化性质间的关系，对排名前十位的细菌门和ARGs与8项土壤理化性质进行皮尔森相关性分析。如表4所示，除了放线菌门(Actinobacteria)、MexF和CeoB以外，其余细菌门及ARGs与某些理化性质均存在显著相关性。变形菌门(Proteobacteria)、拟杆菌门(Bacteroidetes)、MacB、BcrA、MexW、CarA、VanRA和VanRC除了与土壤pH呈负相关外，与其余因子均呈正相关关系；与此相反，酸杆菌门(Acidobacteria)、绿弯菌门(Chloroflexi)、Bacteria_noname、浮霉菌门(Planctomycetes)、厚壁菌门(Firmicutes)、MexF和tetPB除了与土壤pH呈正相关外，与其余因子均呈负相关关系。土壤OM、TN和AN与门水平细菌和ARGs的相关性呈现相似的规律；TP和AP与门水平细菌和ARGs的相关性呈现相同的规律。综上所述，土壤理化性质直接或间接地影响着细菌群落和ARGs的组成和变异。

3 讨 论

3.1 不同有机/无机肥配施对土壤细菌多样性的影响

不同有机/无机肥配施影响了细菌的α-多样性，细菌Chao1指数和ACE指数在N2P1M1处理中高于其他处理，分析认为该施肥量下土壤环境(碳源、水分、pH)为细菌生长和繁殖提供了良好的条件，更有利于保持细菌多样性。土壤养分的可利用性通常与贫营养微生物向富营养微生物的转化有关。例如，土壤中较高的养分利用率促进了富营养微生物的生长和繁殖，而在营养有限的土壤环境中，生长缓慢的贫营养微生物的数量会增加[32-33]。变形菌门、酸杆菌门、放线菌门、芽单胞菌门和绿弯菌门是本研究土壤样品中丰度最高的五大菌门，占细菌总数的82%以上，与之前不同系统报道的结果一致[34-36]。在这些菌门中，变形菌门被认为是一种广泛存在于土壤中的富营养微生物，在粪便处理中其丰度相对较高；酸杆菌门内的大多数类群被认为是生长缓慢的贫营养细菌[32]。本研究中，N3P2M3和N0P0M6处理中变形菌门的相对丰度显著高于N0P0M0处理(P<0.05)；与之相反，随着养分投入量增加，酸杆菌门丰度呈现下降趋势，贫营养微生物绿弯菌门的相对丰度在N0P0M6处理中降到最低。

图6 抗生素抗性基因与门水平细菌的冗余分析Fig.6 Redundancy analysis of ARGs and bacteria at phylum level

3.2 不同有机/无机肥配施对土壤抗生素抗性基因的影响

随着养殖业的集约化及兽用抗生素的广泛使用，畜禽粪便已经成为残留抗生素的主要储存库，其施用到耕地土壤后，会明显激发土壤微生物的抗性及ARGs的丰度[37]。本研究中，PCA和Permanova表明，不同有机/无机肥配施处理ARGs存在显著差异(P<0.001)，分析认为施肥通过改变土壤养分影响了土壤中细菌的活性，进而导致ARGs出现差异。

含量较多的ARGs有大环内酯类、多药类、万古霉素类、四环素类、氯霉素类，这些ARGs普遍存在于畜禽粪便、草鱼肠道、土壤及水体中[38-41]，相对应的这几类抗生素也是广谱性抗生素。其中万古霉素是一种临床上重要的糖肽类抗生素，对革兰氏阳性菌和一些多重耐药菌具有极好的杀灭效果，被称为人类健康的最后一道防线[42-43],但在本长期定位试验中，万古霉素类的抗性基因亚型最多，这不得不引起我们的警惕。

本研究中未施用牛粪的处理N0P0M0中仍然有大量的ARGs，这进一步说明ARGs的存在具有普遍性[37]。分析其原因认为土壤微生物本底中含有ARGs，也有可能是长年累月施肥中，ARGs通过水平转移到达空白土壤[44]。ARGs中的质粒、转座子、整合子等可移动因子可以通过接合、转化和转导等方式从一种菌株转移到其他菌种中，从而使后者携带相应的抗性，并通过微生物的大量繁殖将ARGs传递下去。因此当微生物菌群一旦建立对抗生素的抗性，即便抗生素的选择压力消失，这类菌群也会持续对此类抗生素产生抵抗作用[45]。

4 结 论

(1)不同有机/无机肥配施处理改变了细菌及ARGs的多样性，且不同处理间ARGs差异显著(P<0.001)。

表4 优势细菌门及抗生素抗性基因与土壤理化性质的皮尔森相关性分析Table 4 Pearson correlation analysis of dominant phylum-bacteria and ARGs with soil physicochemical properties

(2)不同有机/无机肥配施处理褐土中，5个主要的细菌门依次为变形菌门、酸杆菌门、放线菌门、芽单胞菌门和绿弯菌门。

(3)不同有机/无机肥配施处理褐土中，相对丰度较高的ARGs种类为大环内酯类、多药类、万古霉素类、四环素类、氯霉素类；其中相对丰度最高的ARGs为MacB、BcrA和MexW。

(4)土壤理化性质直接或间接地影响着细菌群落及ARGs的组成和变异，而ARGs的变化受细菌群落结构的影响。

ARGs已经成为威胁人类健康的潜在污染物，农业生产中施用畜禽粪污时应充分腐熟并与化肥配合施用，以保障生态环境安全。在未来研究中，探索有机肥与化肥的合理配比既可以保障粮食安全，又可以合理规避ARGs带来的潜在风险。