基于高通量测序的昭通核心烟区植烟土壤真菌多样性分析

倪霞，龚林，陈华，宋明健，周玉忠，陈敏，付业明，龚浩瀚，盘文政，张孟升，罗启鹏，李建云，鲍宏明，李德文，张海鹍

（1云南省烟草公司昭通市公司，云南昭通 657000；2昭通市植保植检站，云南昭通 657000；3云南云叶化肥股份有限公司，昆明 650217）

0 引言

烟草作为云南的支柱产业之一，多年来由于化肥的大量施用及不合理的耕作方式，使得植烟土壤板结、酸化等现象日突出，制约了烟叶品质的提升[1]。要提升烟叶品质，首先就得改良土壤结构，增强土壤肥力，增加土壤微生物多样性，以促进农业可持续发展[2-4]。真菌作为土壤微生物的重要组成部分，其能分解有机物质，促进植株生长[5]。因此，研究植烟土壤真菌多样性对提升烟叶品质具有重要意义。研究表明，长期马铃薯连作地块土壤中的真菌数显著高于其他种植模式，真菌病害易多发、高发[6]。张艺洁等[7]利用高通量测序技术，测试施肥对连作植烟土壤微生物多样性的影响，结果表明，有机质是影响土壤真菌群落多样性的主要因素。陈丹梅等[8]利用454高通量测序技术，探究轮作对植烟土壤真菌群落的影响，发现长期轮作可极大的改变土壤中真菌的组成及类型。杨灵朋等[9]采用Biolog-ECO 微平板技术测定烟草根腐病发生的土壤微生物功能多样性，无病土壤的土壤微生物的丰富度和均匀度高于发病土壤，且无病土壤中的微生物对碳源的利用能力更强。近年来，微生物多样性的研究方法层出不穷出不穷，包括平板培养法、BIOLOG生态平板法、PCR-DGGE法以及高通量测序技术等[10-15]，其中高通量测序较其它方法能快速、客观的反应微生物的数量及多样性[16]。本研究采用高通量测序技术，定量、定向的分析研究昭通植烟土壤真菌组成、指示物种、多样性，并对真菌群落的功能进行预测。以期全面深入了解昭通核心烟区土壤真菌的多样性，针对其表征特点进行合理施肥耕作规划，提升烟叶品质。

1 材料与方法

1.1 供试土壤样品

供试土壤样品采集于昭通市昭阳区、鲁甸、巧家、镇雄、彝良、威信、大关和永善8 个县的388 个土壤样品。具体情况如表1所示。

表1 供试土壤样品

1.2 试验设计

不同植烟土壤按新烟区（A组）、连作区（B组）、轮作区（E组）分类进行组间比较。

1.3 试验方法

1.3.1 DNA 提取微生物DNA 的提取使用HiPure Soil DNA提取试剂盒(Magen,Guangzhou,China)按照操作指南进行。

1.3.2 真菌ITS rDNA 测序提取基因组DNA，用带有barcode（测序接头）的特异引物扩增rDNA 的保守区，引物序列信息如表2。扩增体系：50 μL 混合物，包含10 μL 5 × Q5@ Reaction Buffer,10 μL 5 × Q5@ High GCEnhancer，1.5 μL 2.5 mM dNTPs，1.5 μL上下游引物(10 μM),0.2 μL Q5@ High-Fidelity DNAPolymerase,50 ng 模板DNA。扩增条件：95℃5 min，然后95℃1 min，60℃1 min，72℃1 min 进行30 个循环，最后72℃7 min。从2%琼脂糖凝胶中搜集扩增子，使用凝胶提取试剂盒(Axygen Biosciences,Union City,CA,USA)进行纯化，并使用ABI StepOnePlus 实时PCR 系统(Life Technologies,Foster City,USA)进行定量，将纯化后的扩增子在Illumina 平台上进行双端测序(PE250)。

表2 引物序列信息

1.4 生物信息分析

1.4.1 质控统计用FASTP[17]对Illumina 平台的原始数据进行过滤，获得的clean reads 用于组装分析。使用FLASH[18]将clean reads 按最小重叠为10 bp，错配率最高2%的阈值合并为tag。依据Bokulich 等[19]的过滤条件，对低质量tag 进行过滤，得到高质量的clean tag。使用UPARSE[20]流程将clean tag 按≥97%相似度聚类为OTUs（operational taxonomic units，操作分类单元）。

1.4.2 物种组成分析使用Krona[21]显示每个物种分类的丰度统计数据。物种丰度堆叠图使用R 语言ggplot2包[22]展示。

1.4.3 指示物种分析使用R 语言VennDiagram 包[23]进行Venn 分析和R 语言UpSetR 包[24]分析组间共有特有物种/OTU。两组间物种丰度差异分析使用R 语言Vegan 包[25]进行Welch’st 检验和Wilcoxon 秩和检验计算。使用LEfSe 软件[26]对各分组的Biomarker 物种进行筛选。

1.4.4 Alpha 多样性分析 Chao1、ACE、Shannon、Simpson、Good’s coverage多样性指数在QIIME[27]中计算。使用R 语言ggplot2 包[22]绘制多样性指数稀释曲线。多组间的Alpha 指数比较使用R 语言Vegan 包[25]（版本2.5.3）进行Tukey’s HSD 检验和Kruskal-Wallis H检验。

1.4.5 功能预测使用FUNGuild[28]推断真菌的功能组(guild)。使用PICRUSt 进行样本真菌MetaCyc 数据库(MetaCycMetabolic Pathway Database)pathway 功能预测。

2 结果与分析

2.1 土壤真菌群落组成分析

结合OTUs 分类关系，核心烟区土壤真菌群落涵盖30个门、79个纲、202个目、443个科和976个属。

由图1和表3可以看出，真菌群落门水平相对丰度在0.05%以上的门分别为子囊菌门(73.95%～81.36%)、被孢霉门(9.48%～16.05%)、担子菌门(6.68%～7.72%)、壶菌门(0.48%～0.86%)、毛霉门(0.23%～0.45%)和球囊菌门(0.06%～0.11%)。子囊菌门为昭通核心烟区土壤真菌群落的绝对优势门，相对丰度占比73.95%～81.26%。三类土壤真菌群落相对丰度前3 的门均为：子囊菌门、被孢霉门和担子菌门。

图1 不同真菌类群在门水平的相对丰度

表3 不同真菌类群在门水平的相对丰度物种相对丰度表

由表4 和图2 可知，属水平相对丰度前10 的分别为被孢霉属(9.48%～16.05%)、青霉菌属(7.98%～9.21%)、拟青霉菌属(1.55%～6.72%)、曲霉菌属(1.80%～4.00%)、假裸囊菌属(1.33%～4.11%)、镶刀菌属(1.75%～3.65%)、腐质霉属(1.89%～2.56%)、壳二胞菌属(1.29%～2.55%)、篮状菌属(0.81%～1.89%)和柱孢属(0.51%～2.15%)。腐质霉属和壳二胞菌属新烟区的土壤中相对丰度较高，青霉菌属、拟青霉菌属、假裸囊菌属和篮状菌属在连作区相对丰度较高，被孢霉属、曲霉菌属、镶刀菌属和柱孢属在轮作区相对丰度较高。被孢霉属和青霉菌属为昭通核心烟区土壤真菌群落的优势属，相对丰度占比分别为9.48%～16.05%和7.98%～9.21%。三类土壤真菌群落相对丰度前2的属均为被孢霉属和青霉菌属，新烟区土壤真菌群落相对丰度第3的为腐质霉属，连作区的为拟青霉菌属，轮作区的为镶刀菌属。

图2 不同真菌类群在属水平的相对丰度

表4 不同真菌类群在属水平的相对丰度物种相对丰度表

2.2 真菌群落指示物种分析

2.2.1 不同类型土壤真菌群落中优势物种差异性分析由图3可知，在门水平新烟区土壤真菌群落被孢霉门、担子菌门、球囊菌门的相对丰度显著高于连作区，毛霉门的相对丰度显著高于轮作区；子囊菌门连作区相对丰度显著高于新烟区，子囊菌门和毛霉门的连作区相对丰度显著高于轮作区；壶菌门的轮作区相对丰度显著高于新烟区，被孢霉门和壶菌门轮作区相对丰度显著高于连作区。由图4 可知，在属水平壳二胞菌属的相对丰度新烟区显著高于连作区和轮作区，拟青霉菌属、假裸囊菌属和篮状菌连作区相对丰度显著高于新烟区和轮作区，轮作区镶刀菌属的相对丰度显著高于新烟区、被孢霉属的相对丰度显著高于连作区。

图3 真菌群落门水平相对丰度Wilcoxon秩和检验

图4 真菌群落属水平相对丰度Wilcoxon秩和检验

2.2.2 不同类型土壤真菌群落中特有物种分析由表5、图5可知，新烟区土壤真菌群落特有物种主要涵盖：1个门、4个纲、5个目、24个科和67个属；连作区特有物种涵盖：2个目、13个科和30个属；轮作区特有物种涵盖1个门、3个纲、5个目、13个科和32个属。整体来看，新烟区土壤真菌群落特有物种多于连作区和轮作区。

图5 物种VENN图

表5 昭通核心烟区不同类型土壤真菌群落特有物种列表

2.2.3 不同类型土壤真菌群落中差异指示物种分析由图6可以看出，新烟区、连作区和轮作区土壤样本群落组成大部分重叠，利用LEfSe软件对新烟区、连作区和轮作区土壤进行线性判别分析(LDA)，未探寻到引起不同类型植烟土壤产生显著性差异的指示物种，进一步表明新烟区、连作区和轮作土壤的真菌群落结构组成差异不显著。

图6 昭通核心烟区不同类型土壤真菌群落组间比较分析

2.3 真菌群落多样性特征分析

从表6中可知，新烟区土壤真菌群落的Sobs指数、Chao指数和ACE指数最高，轮作区次之，连作区最低；新烟区香浓指数最高，连作区次之，轮作区最低；连作区辛普森指数最高，新烟区次之，轮作区最低。由图7可以看出，新烟区的Sobs指数、Chao指数和ACE指数与连作区和轮作区均有极显著差异，香农指数与轮作区有极显著差异、连作区有显著差异，辛普森指数与连作区和轮作区差异均不显著；连作区的Sobs指数和香浓指数与轮作区有显著差异，Chao 指数、ACE 指数和辛普森指数与轮作区有极显著差异。由此可见新烟区土壤真菌群落的生物多样性最高，物种丰富度显著高于连作区和轮作区；轮作区土壤真菌群落的物种丰富度显著高于连作区，物种均匀度显著低于连作区，真菌群落生物多样性略高于连作区。

图7 昭通核心烟区土壤真菌群落的生物多样性指数统计检验结果

表6 昭通核心烟区土壤真菌群落的生物多样性指数

由图8可知A、B、E 3个分组的稀释曲线均趋于平缓，说明包括的稀有物种在内的土壤真菌均得到较好分析，比较真实地反映了该研究区的真菌群落组成。

图8 昭通核心烟区土壤真菌的稀释曲线

2.4 真菌群落的FUNGuild功能预测

采用FUNGuild 对昭通核心烟区土壤真菌群落的营养型和功能群进行预测的结果见表7 和图9。从表和图中可以看出，昭通核心烟区土壤真菌群落可被区分为6个不同的营养型，分别为病理营养型、病理—腐生营养型、病理—腐生—共生营养型、病理—共生营养型、腐生营养型、腐生—共生营养型和共生营养型，其中病理营养型、腐生营养型和腐生—共生营养型的真菌群落为主要组成，占比分别为7.50%～14.17%、10.28%～15.02%和10.55%～17.61%。不同营养型的真菌群落，在不同类型的土壤中丰度值有差异。如图10所示，病理营养型和腐生营养型的真菌群落以连作区的最高，病理—腐生营养型、病理—共生营养型和腐生—共生营养型以新烟区和轮作区的较高，病理—腐生—共生营养型以轮作区的最高，共生营养型以新烟区的最高。整体来看，新烟区真菌群落占比前1至前3的营养型分别为腐生—共生营养型、腐生营养型和病理营养型，连作区分别为腐生营养型、病理营养型和腐生—共生营养型，轮作区分别为腐生—共生营养型、病理—腐生—共生营养型和腐生营养型，说明轮作有利于降低病理营养型的真菌群落，能减少真菌性病害发生的几率。

图9 昭通核心烟区真菌群落功能分布图

图10 真菌群落功能丰度热图

表7 昭通核心烟区真菌群落功能分布总览表

3 结论与讨论

3.1 昭通核心烟区土壤微真菌生物群落组成

昭通核心烟区土壤真菌群落涵盖30 个门、79 个纲、202 个目、443 个科和976 个属，其中子囊菌门、被孢霉门、担子菌门、壶菌门、毛霉门和球囊菌门真菌群落门水平相对丰度均在0.05%以上。昭通核心烟区土壤真菌群落的绝对优势门为子囊菌门，相对丰度占比73.95%～81.26%。三类土壤真菌群落相对丰度前3 的门均为：子囊菌门、被孢霉门和担子菌门。属水平相对丰度前10 的分别为被孢霉属、青霉菌属、拟青霉菌属、曲霉菌属、假裸囊菌属、镶刀菌属、腐质霉属、壳二胞菌属、篮状菌属和柱孢属。被孢霉属和青霉菌属为昭通核心烟区土壤真菌群落的优势属，相对丰度占比分别为9.48%～16.05%和7.98%～9.21%。三类土壤真菌群落相对丰度前2的属均为被孢霉属和青霉菌属，新烟区土壤真菌群落相对丰度第3的为腐质霉属，连作区的为拟青霉菌属，轮作区的为镶刀菌属。

3.2 昭通核心烟区土壤真菌群落指示物种

土壤真菌群落壳二胞菌属的相对丰度新烟区显著高于连作区和轮作区，拟青霉菌属、假裸囊菌属和篮状菌属连作区相对丰度显著高于新烟区和轮作区，轮作区镶刀菌属的相对丰度显著高于新烟区、被孢霉属的相对丰度显著高于连作区。新烟区土壤真菌群落特有物种多于连作区和轮作区，新烟区、连作区和轮作区土壤样本真菌群落组成大部分重叠，不同类型土壤的真菌群落结构组成差异不显著。

3.3 昭通核心烟区土壤真菌群落生物多样性特征

新烟区土壤真菌群落的生物多样性最高，物种丰富度显著高于连作区和轮作区；土壤真菌群落的物种丰富度轮作区显著高于连作区，物种均匀度轮作区显著低于连作区，真菌群落生物多样性轮作区略高于连作区。

3.4 昭通核心烟区土壤真菌群落的功能

昭通核心烟区土壤真菌群落可被区分为6个不同的营养型，分别为病理营养型、病理—腐生营养型、病理—腐生—共生营养型、病理—共生营养型、腐生营养型、腐生—共生营养型和共生营养型，其中病理营养型、腐生营养型和腐生—共生营养型的真菌群落为主要组成，占比分别为7.50%～14.17%、10.28%～15.02%和10.55%～17.61%。不同营养型的真菌群落，在不同类型的土壤中丰度值有差异。整体来看，新烟区真菌群落占比前1至前3的营养型分别为腐生—共生营养型、腐生营养型和病理营养型，连作区分别为腐生营养型、病理营养型和腐生—共生营养型，轮作区分别为腐生—共生营养型、病理—腐生—共生营养型和腐生营养型，说明轮作有利于降低病理营养型的真菌群落，能减少真菌性病害发生的几率。