桑葚病原菌高通量鉴定及种属分类方法及其应用与流程

本发明涉及植物病原菌鉴定及种属分类技术领域，更具体地，涉及一种桑葚病原菌高通量鉴定及种属分类的方法。

背景技术：

近年来，随着蚕桑多元化发展和桑产业的兴起，桑不仅是以桑叶作为蚕桑产业中特种经济动物——蚕的饲料，更因其桑树本身的药理及营养价值，使得桑的多用途资源利用变得更加丰富。桑树全身都是宝，桑叶、桑椹和桑枝是目前种桑的主要目的收获物，桑椹作为桑科桑属植物的果实，又称桑果、桑葚或桑枣，被誉为“民间圣果”，自古以来就是百姓常采用的一种利尿，保健，消暑的鲜果，可制作果干、果汁，也可来泡酒。桑椹性味甘寒，具有补肝益肾、生津润燥、乌发明目等功效。因此桑果是人们植桑时，因桑葚具有较高的营养及食用价值，故桑果是除收获桑叶外具有较大经济价值的桑产物(陈冬梅，2013)。

对桑果威胁最大的病害为桑菌核病。桑菌核病是果桑一种主要病害，俗称桑白果病，属真菌类病害；桑树开花时病菌开始侵入，桑果快成熟时病状显现，颜色呈白色，失去经济价值，大多数的果桑品种均较易发病，发病率可高达90％以上(蒯元璋等，2012)。

桑椹菌核病(popcorn disease)分桑椹肥大性菌核病、桑椹缩小性菌核病和桑椹小粒性菌核病3种，俗称白果病。桑椹菌核病病菌以菌核在土壤中越冬，越冬后春季温暖、多雨、土壤潮湿利于菌核萌发，产生子囊盘多，病害重。在通风透光差，低洼多湿的桑园更适合病原菌的萌发与生长；花果多，树龄老的桑园，随着果桑养成年份的增加，菌核病病原(菌核或子囊盘)的田间残留量增加，也加大了菌核病爆发的可能性(娄利峰等，2016)。

桑椹肥大性菌核病的植株花被厚肿，颜色呈现灰白，果实比较膨大，果实的中心有个黑色菌核，得病的桑椹破碎后会闻到一股臭味；感染缩小性菌核病的果实缩小，明显颜色呈现灰白，质地变得比较坚硬，表面会有暗褐色的细斑长出，内部存在黑硬菌核；小粒性菌核病的症状是果实体积膨大，果实内部长有小粒形菌核，果实呈现出病态的灰黑色，易脱落。

桑椹肥大性菌核病病菌Ciboria shiraiana.称白杯盘菌(桑实杯盘菌)，属子囊菌门真菌。分生孢子梗丛生，基部粗顶端细小，上生分生孢子。分生孢子单胞，卵形，无色。菌核萌发产生1-5个子囊盘，盘内生子囊，侧丝细长，内有8个子囊孢子，子囊孢子椭圆形，无色单胞，具隔膜1-2个。

桑椹缩小性菌核病病原菌为白井地杖菌Mitrula shiraiana，其分生孢子梗细丝状，具分枝，端生，卵至椭圆形、无色。单胞型分生孢子菌核会以两种方式长出子实体。子实体有灰褐色、扁平长柄，长有茸毛；子实体头部长椭圆形，具数条浅褐色纵向排列的稻纹，子囊生在头部外侧子实层里，长棍棒状，先端圆，基部细，内生8个子囊孢子，子囊孢子单胞、无色、椭圆形。

桑椹小粒性菌核病病原菌为肉阜状杯盘菌Ciboria canrunculoides，分生孢子接近于球形，较圆，子囊盘杯状，具长柄，子囊的形状近似于长棍棒，内部含有8个子囊孢子，子囊孢子单胞，无色，肾形，有半球形小体附着。侧丝有分枝，有隔或无隔。

桑菌核病致病菌除上述的桑实杯盘菌、白井地杖菌、肉阜状杯盘菌外，最近陈倩茜(2015)报道从桑菌核病的样本中分离确定了一种新的病原菌，鉴定为拟茎点霉属(Phomopsis)的病原真菌Phomopsis sp.SC1104，该病原菌可能与桑椹菌核病流行爆发有关。桑菌核病致病菌的寄主多样，病原可危害32科160多种植物，而桑菌核病还可与油菜等十字花科植物菌核病病原菌交叉感染(范锦等，2014；吕蕊花等，2015)，这给桑菌核病的防治带来了困难。

植物病原菌分类一般以寄主的名称作为真菌种名的命名原则，同科同属且存在于同一寄主上的真菌无法在命名上和形态上区分开来，而若同种病原菌在不同的寄主间出现交叉感染情况，该病原菌的命名变得复杂而难以鉴别或区分，类似桑葚菌核病有多种病原菌感染的情况，对于其主要感染病原菌的鉴别，则需要分子生物学的技术加以辅助鉴别。

真菌的系统分类方法中，更多的是以遗传物质的信息来佐证传统形态学的分类，并纠正一些主观分类的错误。现有的研究方法中，常用核糖体DNA(ribosome DNA,rDNA)序列进行测序比对。rDNA中编码18S、5.8S、28S的rDNA序列较为保守，可用于真菌目、科、属等分类阶元的系统发育研究；内转录间隔区ITS位于核糖体rDNA 18S、5.8S及28S之间，由于不需要加入成熟核糖体，因此在进化过程中能够承受更多的变异，其进化速率为18S rDNA的10倍，属于中度保守的区域，利用它可研究种及种以下的分类阶元。

全基因组测序作为已经成熟且推广开来的一种测序，多用于基因的功能注释或进一步的比较基因组学研究。全基因组的高通量测序多用于已知菌株的功能基因的查找鉴定，用于真菌的种属鉴定，其准确性毋庸置疑。基于分子水平判定物种的技术已较为成熟，而测序的快速与成本低廉也使得这项技术变得更大众化。针对现有的桑葚菌核病病原菌的鉴定及种属分类研究中，受传统研究方法的限制，缺乏高通量技术、结果往往欠准确性，对桑葚菌核病病原菌进行定量检测方面几乎空白，因此本发明的方法具有前瞻性。

技术实现要素：

本发明要解决的技术问题是针对现有对桑葚病原菌鉴定不够准确的技术不足，提供一种准确性更高同时种属分类明确的鉴定及物种分类方法，并且检测对象是基于桑葚。

本发明要解决的另一技术问题是提供所述方法在对桑葚病原菌进行定量分析方面的应用。

本发明还一要解决的技术问题是提供所述方法在分析微生物物种丰度方面的应用。

本发明还一要解决的技术问题是提供所述方法在计算物种间遗传距离方面的应用。

本发明的目的通过以下技术方案予以实现：

提供一种桑葚病原菌高通量鉴定及种属分类的方法，利用高通量测序手段及生物信息学分析方法，建立一种桑葚病原菌的快速、准确、高效的鉴定及种属分类方法。

具体地，所述鉴定及种属分类方法包括以下步骤：

S1.收集桑葚病果；

S2.提取桑葚病果的总DNA；

S3.Illumina DNA文库构建；

S4.Illumina高通量测序；

S5.去除测序数据中桑树基因组序列；

S6.组装微生物基因组序列；

S7.组装完整核糖体序列；

S8.筛选标记微生物核糖体DNA序列；

S9.对比分析核糖体DNA序列进行种属分类；

其中，步骤S2所述提取桑葚病果的总DNA的方法为：剪取桑病果中的患病区域，加入液氮充分研磨，用真菌DNA提取试剂盒按照其操作说明提取病叶总DNA，所述总DNA保存于-20℃冰箱；

步骤S3所述Illumina DNA文库构建的方法为：依照Illumina文库构建流程，将所述步骤S2中所述总DNA构建为片段大小400-600bp范围内的的双末端高通量测序文库；

步骤S4所述Illumina高通量测序的方法为：用Illumina Hiseq 2500测序仪对所述步骤S3中所述DNA文库进行高通量测序。

步骤S5所述去除测序数据中桑树基因组序列的方法为：利用比对软件对步骤S4中所述高通量测序数据进行数据比对分析。选择比对算法，将所述测序数据与桑树参考基因组进行比对，并将比对上参考基因组的所述测序数据判定为桑树基因组序列。使用编写的计算机程序从所述测序数据中去除桑树基因组序列。

S5所述去除桑树的基因组DNA序列选用的参考基因组序列为：

桑属川桑(Morus notabilis)全基因组序列(GCA_000414095.2)和叶绿体基因组序列(NC_027110.1)。

优选地，所述比对软件为bwa(0.7.12-r1039)软件；

优选地，所述比对算法为mem比对算法；

优选地，所述将测序数据与桑树参考基因组进行比对选用的是双末端比对方法和bwa(0.7.12-r1039)软件的默认参数；

优选地，所述编写的计算机程序是用python计算机语言编写。

步骤S6所述组装微生物基因组序列的方法为：使用组装软件对步骤S5所述已去除桑树基因组序列的测序数据进行组装。

优选地，所述组装软件为MetaVelvet(v1.2.01)。

步骤S7所述组装完整核糖体DNA序列的方法为：采用比对软件，对所述组装序列进行比对，根据比对结果从测序数据中获取双末端测序片段，使用组装软件对序列进行组装和延伸，经多个循环操作，直至获得完整的核糖体DNA序列；

优选地，所述比对软件为bwa(0.7.12-r1039)软件；

优选地，所述比对方法采用无错配(0mismatch)和无空缺(0gap)比对；

优选地，所述组装软件为MetaVelvet(v1.2.01)软件。

步骤S8所述筛选标记微生物核糖体DNA序列的方法为：使用序列比对分析软件，设定比对期望值，将步骤S6所述微生物基因组序列与数据库进行比对，根据比对结果对序列标签进行注释及物种分析；

优选地，所述比对期望值<1e-20；

优选地，所述序列比对分析软件为blastn(2.2.31+)软件；

优选地，所述数据库为NCBI数据库中的nt数据库。

步骤S9所述对比分析核糖体DNA序列进行种属分类的方法为：使用软件1对所述微生物核糖体DNA序列进行统计分析，选取系统树构建模型，将所述模型代入软件2中，构建系统进化树。

优选地，所述软件1为jModelTest2(https://github.com/ddarriba/jmodeltest2)在线软件；

优选地，所述系统树构建模型为取AIC(Akaike Information Criterion)值最小的替代模型GTR+G；

优选地，所述软件2为RaxML(8.1.5)；

优选地，所述构建系统进化树的方法为使用最大似然法构建的系统进化树，进化树迭代次数为1000次。

本发明同时提供所述桑葚病原菌高通量鉴定及种属分类方法在桑树病原菌定量分析方面的应用。

本发明还提供所述桑葚病原菌高通量鉴定及种属分类方法在分析微生物物种丰度方面的应用。

所述应用的方法为：使用分析软件，计算核糖体DNA片段的平均测序深度，以此作为该物种的丰度值。

优选地，所述分析软件为bwa(0.7.12-r1039)与samtools(v1.2)。

本发明还提供所述桑葚病原菌高通量鉴定及种属分类方法在计算物种间遗传距离方面的应用。

附图说明

图1桑葚菌核病病果所检测到的真菌微生物分类树。

图2验证核糖体组装结果的两对引物与两对通用引物ITS1/ITS4、ITS4/ITS5电泳图。

图3核糖体DNA序列组成与碱基GC比例分布图。

图4菌核病病原体Ciboria carunculoides 18S进化树。

图5菌核病病原体Ciboria carunculoides ITS进化树。

图6菌核病病原体Ciboria carunculoides28S进化树。

图7菌核病病原体rRNA全基因序列进化树。

具体实施方式

下面结合具体实施例进一步说明本发明方法。下述实施例和附图仅用于示例性说明，不能理解为对本发明的限制。除非特别说明，下述实施例中使用的试剂原料为常规市购或商业途径获得的生试剂原料，除非特别说明，下述实施例中使用的方法和设备为本领域常规使用的方法和设备

本发明以对桑葚菌核病病原菌的鉴定及种属分类为例，根据本发明思想，本领域技术人员可以参照现有技术，对其他病症的病原菌进行鉴定及种属分类。

实施例1

在发病的桑园寻找到具有典型桑葚菌核病病征的桑果，收集，剪取桑病果中的患病区域，剪下的材料使用液氮充分研磨，总DNA的提取使用康为世纪公司的真菌DNA提取试剂盒，具体按照其操作说明进行，提取后的总DNA保存于-20℃。依照Illumina文库构建流程，将总DNA构建为片段大小450bp的双末端高通量测序文库，使用Illumina Hiseq2500测序仪对构建好的DNA文库进行高通量测序，共测得9.03M对测序片段，测序读长为双末端125bp，总测序数据量2.26Gb。

由于DNA提取过程中包含桑树的桑果材料，为减少测序数据中桑树基因组数据对微生物序列组装的影响，在进行微生物序列组装前首先去除桑树的基因组DNA序列。选取桑树Morus notabilis全基因组序列(GCA_000414095.2)和叶绿体基因组序列(NC_027110.1)作为参考基因组序列，采用bwa(0.7.12-r1039)比对软件进行数据比对分析。比对选择mem比对算法，使用双末端比对方法和软件的默认参数，将测序数据与桑树参考基因组进行比对，并将比对上参考基因组的测序片段判定为桑树基因组序列。使用python编写的计算机程序从fastq测序数据中去除桑树测序数据，然后再进入微生物序列组装。微生物序列的组装使用MetaVelvet(v1.2.01)组装软件进行。

序列标签注释使用blastn(2.2.31+)序列比对分析软件，组装的序列标签序列与NCBI的nt数据库进行比对，blastn比对设定期望值<1e-20，根据比对结果对序列标签进行注释。核糖体DNA序列是细菌和真菌鉴定的重要的最常用的分子标记，因此物种分类鉴定和定量以核糖体DNA为主要分子标记。根据序列标签注释结果，选择核糖体DNA序列作为微生物鉴定与定量分析依据。使用bwa(0.7.12-r1039)+samtools(v1.2)分析软件，计算测序数据中核糖体DNA片段的平均测序深度，并以此作为该物种的丰度值。

结果显示，共注释228条rRNA序列标签，数据包含4条细菌序列标签，82条真菌序列标签，分别对应3种细菌，3种真菌。通过对序列标签注释结果的查询，共发现桑病果内存在3个属的真菌微生物，相对丰度最高的为Ciboria属的真菌，见附图1，其分子标签测序深度为1115，相对丰度为98.50，所占比例最高，相关文献报道该属绝大部分为植物病原，大部分导致植物出现果实或种子僵化膨大等病征，与本发明研究的病征相符合。

真菌的核糖体DNA由18S区段，ITS1区段，5.8S区段，ITS2区段和28S区段组成，序列总长度约为5.5Kb。MetaVelvet(v1.2.01)初始组装的序列标签为断裂的核糖体标签，为得到完整的核糖体DNA序列，分析采用序列捕获和从头组装策略，以组装完整的核糖体DNA。选择包含目标病原菌ITS序列的核糖体DNA序列为参考序列，采用bwa(0.7.12-r1039)软件进行0mismatch和0gap比对，根据比对结果从测序数据中获取双末端测序片段，进一步采用MetaVelvet(v1.2.01)组装软件对序列进行组装和延伸，经多个循环操作，获取完整的核糖体DNA序列。

根据组装完整的核糖体DNA序列，采用真菌通用引物ITS1和ITS4、ITS4和ITS5(White，1990)，和基于本研究组装获得的完整核糖体DNA序列设计了PCR验证引物见表1，然后进行PCR扩增和Sanger测序方法，验证组装的正确性。通用引物序列见表2所示，PCR反应体系见表3所示。

表1.PCR验证引物序列表

表2真菌ITS区域通用引物序列

表3 PCR反应体系(20μL)

真菌通用引物ITS1和ITS4引物组PCR程序：预变性94℃5min；94℃30s，55℃30s，72℃1min，30个循环；72℃10min。

真菌通用引物ITS4和ITS5引物组PCR程序：预变性94℃5min；94℃30s，55℃30s，72℃1min，30个循环；72℃10min。

验证核糖体组装结果的P1-F/R引物组PCR反应体系(20μL)PCR程序：94℃5min；94℃1min，56℃1min，72℃4min，30个循环；72℃10min。

验证核糖体组装结果的P2-F/R引物组PCR反应体系(20μL)PCR程序：94℃5min；94℃1min，56℃1min，72℃2min30s，30个循环；72℃10min。

取3μLPCR扩增产物用1.2％的琼脂糖凝胶(EB染色)电泳检测。通过琼脂糖凝胶电泳回收大小对应的PCR产物片段，见附图1所示，图1中，M：Takara DL5000Marker；1.真菌通用引物ITS1和ITS4；2.真菌通用引物ITS4和ITS5；3.验证核糖体组装结果引物P1-F和P1-R；4.验证核糖体组装结果引物P2-F和P2-R。

PCR验证结果见附图2，从图中可以看出，本次样品较纯净，泳道1和2的结果说明用真菌通用引物ITS1和ITS4、ITS4和ITS5扩增的模板DNA是属于真菌的DNA；泳道3和4的结果说明验证核糖体组装结果的P1引物组和P2引物组均扩增到相应的目的条带，进一步测序后的结果与组装的结果高度一致。

结果显示，根据高通量测序数据组装与实验验证结果，组装得到完整核糖体DNA序列，见附图3，长度6202bp(GC比例52.98％)。序列包含18S区域，ITS1区域，5.8S区域，ITS2区域，28S区域。其中18S区域长度1782bp(GC比例50.05％)，ITS1区域长度101bp(GC比例43.56％)，5.8S区域长度151bp(GC比例46.36％)，ITS2区域长161bp(GC比例56.52％)，28S区域长4007bp(GC比例54.63％)。18S区域和28S区域长度较大，占序列总长度的96.37％；而在GC比例上，ITS2区域(56.52％)的GC比例则明显高于其他区域。

构建了对应Ciboria carunculoides的rRNA分子系统树模型，菌核病病原体Ciboria carunculoides 18S进化树、菌核病病原体Ciboria carunculoides ITS进化树、菌核病病原体Ciboria carunculoides 28S进化树、菌核病病原体Ciboria carunculoides的rRNA全基因序列进化树分别见附图4～附图7所示。

使用blastn(2.2.31+)比对软件将核糖体DNA与nt数据库进行比对，对核糖体DNA序列进行注释，共注释出28S区域，ITS1区域，5.8S区域，ITS2区域和28S区域。使用python语言编写的计算机程序对序列的GC比例进行分析，计算各区域的平均GC比例。同时程序以100bp为计算窗口，以10bp为步进，计算核糖体DNA的GC比例特征。

SEQUENCE LISTING

<110> 华南农业大学

<120> 桑葚病原菌高通量鉴定及种属分类方法及其应用

<130>

<160> 7

<170> PatentIn version 3.3

<210> 1

<211> 20

<212> DNA

<213> 引物P1-F

<400> 1

taaataccga tgcggggctc 20

<210> 2

<211> 20

<212> DNA

<213> 引物P1-R

<400> 2

cactaggtcg ggacctctca 20

<210> 3

<211> 20

<212> DNA

<213> 引物P2-F

<400> 3

cgcttagcgt gctacccata 20

<210> 4

<211> 20

<212> DNA

<213> 引物P2-R

<400> 4

acaggaccgt gtaacgcaat 20

<210> 5

<211> 19

<212> DNA

<213> 通用引物ITS1

<400> 5

tccgtaggtg aacctgcgg 19

<210> 6

<211> 20

<212> DNA

<213> 通用引物ITS4

<400> 6

tcctccgctt attgatatgc 20

<210> 7

<211> 22

<212> DNA

<213> 通用引物ITS5

<400> 7

ggaagtaaaa gtcgtaacaa gg 22