本發(fā)明涉及植物病原菌鑒定及種屬分類技術領域，更具體地，涉及一種桑葚病原菌高通量鑒定及種屬分類的方法。

背景技術：

近年來，隨著蠶桑多元化發(fā)展和桑產(chǎn)業(yè)的興起，桑不僅是以桑葉作為蠶桑產(chǎn)業(yè)中特種經(jīng)濟動物——蠶的飼料，更因其桑樹本身的藥理及營養(yǎng)價值，使得桑的多用途資源利用變得更加豐富。桑樹全身都是寶，桑葉、桑椹和桑枝是目前種桑的主要目的收獲物，桑椹作為桑科桑屬植物的果實，又稱桑果、桑葚或桑棗，被譽為“民間圣果”，自古以來就是百姓常采用的一種利尿，保健，消暑的鮮果，可制作果干、果汁，也可來泡酒。桑椹性味甘寒，具有補肝益腎、生津潤燥、烏發(fā)明目等功效。因此桑果是人們植桑時，因桑葚具有較高的營養(yǎng)及食用價值，故桑果是除收獲桑葉外具有較大經(jīng)濟價值的桑產(chǎn)物(陳冬梅，2013)。

對桑果威脅最大的病害為桑菌核病。桑菌核病是果桑一種主要病害，俗稱桑白果病，屬真菌類病害；桑樹開花時病菌開始侵入，桑果快成熟時病狀顯現(xiàn)，顏色呈白色，失去經(jīng)濟價值，大多數(shù)的果桑品種均較易發(fā)病，發(fā)病率可高達90％以上(蒯元璋等，2012)。

桑椹菌核病(popcorn disease)分桑椹肥大性菌核病、桑椹縮小性菌核病和桑椹小粒性菌核病3種，俗稱白果病。桑椹菌核病病菌以菌核在土壤中越冬，越冬后春季溫暖、多雨、土壤潮濕利于菌核萌發(fā)，產(chǎn)生子囊盤多，病害重。在通風透光差，低洼多濕的桑園更適合病原菌的萌發(fā)與生長；花果多，樹齡老的桑園，隨著果桑養(yǎng)成年份的增加，菌核病病原(菌核或子囊盤)的田間殘留量增加，也加大了菌核病爆發(fā)的可能性(婁利峰等，2016)。

桑椹肥大性菌核病的植株花被厚腫，顏色呈現(xiàn)灰白，果實比較膨大，果實的中心有個黑色菌核，得病的桑椹破碎后會聞到一股臭味；感染縮小性菌核病的果實縮小，明顯顏色呈現(xiàn)灰白，質(zhì)地變得比較堅硬，表面會有暗褐色的細斑長出，內(nèi)部存在黑硬菌核；小粒性菌核病的癥狀是果實體積膨大，果實內(nèi)部長有小粒形菌核，果實呈現(xiàn)出病態(tài)的灰黑色，易脫落。

桑椹肥大性菌核病病菌Ciboria shiraiana.稱白杯盤菌(桑實杯盤菌)，屬子囊菌門真菌。分生孢子梗叢生，基部粗頂端細小，上生分生孢子。分生孢子單胞，卵形，無色。菌核萌發(fā)產(chǎn)生1-5個子囊盤，盤內(nèi)生子囊，側(cè)絲細長，內(nèi)有8個子囊孢子，子囊孢子橢圓形，無色單胞，具隔膜1-2個。

桑椹縮小性菌核病病原菌為白井地杖菌Mitrula shiraiana，其分生孢子梗細絲狀，具分枝，端生，卵至橢圓形、無色。單胞型分生孢子菌核會以兩種方式長出子實體。子實體有灰褐色、扁平長柄，長有茸毛；子實體頭部長橢圓形，具數(shù)條淺褐色縱向排列的稻紋，子囊生在頭部外側(cè)子實層里，長棍棒狀，先端圓，基部細，內(nèi)生8個子囊孢子，子囊孢子單胞、無色、橢圓形。

桑椹小粒性菌核病病原菌為肉阜狀杯盤菌Ciboria canrunculoides，分生孢子接近于球形，較圓，子囊盤杯狀，具長柄，子囊的形狀近似于長棍棒，內(nèi)部含有8個子囊孢子，子囊孢子單胞，無色，腎形，有半球形小體附著。側(cè)絲有分枝，有隔或無隔。

桑菌核病致病菌除上述的桑實杯盤菌、白井地杖菌、肉阜狀杯盤菌外，最近陳倩茜(2015)報道從桑菌核病的樣本中分離確定了一種新的病原菌，鑒定為擬莖點霉屬(Phomopsis)的病原真菌Phomopsis sp.SC1104，該病原菌可能與桑椹菌核病流行爆發(fā)有關。桑菌核病致病菌的寄主多樣，病原可危害32科160多種植物，而桑菌核病還可與油菜等十字花科植物菌核病病原菌交叉感染(范錦等，2014；呂蕊花等，2015)，這給桑菌核病的防治帶來了困難。

植物病原菌分類一般以寄主的名稱作為真菌種名的命名原則，同科同屬且存在于同一寄主上的真菌無法在命名上和形態(tài)上區(qū)分開來，而若同種病原菌在不同的寄主間出現(xiàn)交叉感染情況，該病原菌的命名變得復雜而難以鑒別或區(qū)分，類似桑葚菌核病有多種病原菌感染的情況，對于其主要感染病原菌的鑒別，則需要分子生物學的技術加以輔助鑒別。

真菌的系統(tǒng)分類方法中，更多的是以遺傳物質(zhì)的信息來佐證傳統(tǒng)形態(tài)學的分類，并糾正一些主觀分類的錯誤?，F(xiàn)有的研究方法中，常用核糖體DNA(ribosome DNA,rDNA)序列進行測序比對。rDNA中編碼18S、5.8S、28S的rDNA序列較為保守，可用于真菌目、科、屬等分類階元的系統(tǒng)發(fā)育研究；內(nèi)轉(zhuǎn)錄間隔區(qū)ITS位于核糖體rDNA 18S、5.8S及28S之間，由于不需要加入成熟核糖體，因此在進化過程中能夠承受更多的變異，其進化速率為18S rDNA的10倍，屬于中度保守的區(qū)域，利用它可研究種及種以下的分類階元。

全基因組測序作為已經(jīng)成熟且推廣開來的一種測序，多用于基因的功能注釋或進一步的比較基因組學研究。全基因組的高通量測序多用于已知菌株的功能基因的查找鑒定，用于真菌的種屬鑒定，其準確性毋庸置疑。基于分子水平判定物種的技術已較為成熟，而測序的快速與成本低廉也使得這項技術變得更大眾化。針對現(xiàn)有的桑葚菌核病病原菌的鑒定及種屬分類研究中，受傳統(tǒng)研究方法的限制，缺乏高通量技術、結(jié)果往往欠準確性，對桑葚菌核病病原菌進行定量檢測方面幾乎空白，因此本發(fā)明的方法具有前瞻性。

技術實現(xiàn)要素：

本發(fā)明要解決的技術問題是針對現(xiàn)有對桑葚病原菌鑒定不夠準確的技術不足，提供一種準確性更高同時種屬分類明確的鑒定及物種分類方法，并且檢測對象是基于桑葚。

本發(fā)明要解決的另一技術問題是提供所述方法在對桑葚病原菌進行定量分析方面的應用。

本發(fā)明還一要解決的技術問題是提供所述方法在分析微生物物種豐度方面的應用。

本發(fā)明還一要解決的技術問題是提供所述方法在計算物種間遺傳距離方面的應用。

本發(fā)明的目的通過以下技術方案予以實現(xiàn)：

提供一種桑葚病原菌高通量鑒定及種屬分類的方法，利用高通量測序手段及生物信息學分析方法，建立一種桑葚病原菌的快速、準確、高效的鑒定及種屬分類方法。

具體地，所述鑒定及種屬分類方法包括以下步驟：

S1.收集桑葚病果；

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

S3.Illumina DNA文庫構(gòu)建；

S4.Illumina高通量測序；

S5.去除測序數(shù)據(jù)中桑樹基因組序列；

S6.組裝微生物基因組序列；

S7.組裝完整核糖體序列；

S8.篩選標記微生物核糖體DNA序列；

S9.對比分析核糖體DNA序列進行種屬分類；

其中，步驟S2所述提取桑葚病果的總DNA的方法為：剪取桑病果中的患病區(qū)域，加入液氮充分研磨，用真菌DNA提取試劑盒按照其操作說明提取病葉總DNA，所述總DNA保存于-20℃冰箱；

步驟S3所述Illumina DNA文庫構(gòu)建的方法為：依照Illumina文庫構(gòu)建流程，將所述步驟S2中所述總DNA構(gòu)建為片段大小400-600bp范圍內(nèi)的的雙末端高通量測序文庫；

步驟S4所述Illumina高通量測序的方法為：用Illumina Hiseq 2500測序儀對所述步驟S3中所述DNA文庫進行高通量測序。

步驟S5所述去除測序數(shù)據(jù)中桑樹基因組序列的方法為：利用比對軟件對步驟S4中所述高通量測序數(shù)據(jù)進行數(shù)據(jù)比對分析。選擇比對算法，將所述測序數(shù)據(jù)與桑樹參考基因組進行比對，并將比對上參考基因組的所述測序數(shù)據(jù)判定為桑樹基因組序列。使用編寫的計算機程序從所述測序數(shù)據(jù)中去除桑樹基因組序列。

S5所述去除桑樹的基因組DNA序列選用的參考基因組序列為：

桑屬川桑(Morus notabilis)全基因組序列(GCA_000414095.2)和葉綠體基因組序列(NC_027110.1)。

優(yōu)選地，所述比對軟件為bwa(0.7.12-r1039)軟件；

優(yōu)選地，所述比對算法為mem比對算法；

優(yōu)選地，所述將測序數(shù)據(jù)與桑樹參考基因組進行比對選用的是雙末端比對方法和bwa(0.7.12-r1039)軟件的默認參數(shù)；

優(yōu)選地，所述編寫的計算機程序是用python計算機語言編寫。

步驟S6所述組裝微生物基因組序列的方法為：使用組裝軟件對步驟S5所述已去除桑樹基因組序列的測序數(shù)據(jù)進行組裝。

優(yōu)選地，所述組裝軟件為MetaVelvet(v1.2.01)。

步驟S7所述組裝完整核糖體DNA序列的方法為：采用比對軟件，對所述組裝序列進行比對，根據(jù)比對結(jié)果從測序數(shù)據(jù)中獲取雙末端測序片段，使用組裝軟件對序列進行組裝和延伸，經(jīng)多個循環(huán)操作，直至獲得完整的核糖體DNA序列；

優(yōu)選地，所述比對軟件為bwa(0.7.12-r1039)軟件；

優(yōu)選地，所述比對方法采用無錯配(0mismatch)和無空缺(0gap)比對；

優(yōu)選地，所述組裝軟件為MetaVelvet(v1.2.01)軟件。

步驟S8所述篩選標記微生物核糖體DNA序列的方法為：使用序列比對分析軟件，設定比對期望值，將步驟S6所述微生物基因組序列與數(shù)據(jù)庫進行比對，根據(jù)比對結(jié)果對序列標簽進行注釋及物種分析；

優(yōu)選地，所述比對期望值<1e-20；

優(yōu)選地，所述序列比對分析軟件為blastn(2.2.31+)軟件；

優(yōu)選地，所述數(shù)據(jù)庫為NCBI數(shù)據(jù)庫中的nt數(shù)據(jù)庫。

步驟S9所述對比分析核糖體DNA序列進行種屬分類的方法為：使用軟件1對所述微生物核糖體DNA序列進行統(tǒng)計分析，選取系統(tǒng)樹構(gòu)建模型，將所述模型代入軟件2中，構(gòu)建系統(tǒng)進化樹。

優(yōu)選地，所述軟件1為jModelTest2(https://github.com/ddarriba/jmodeltest2)在線軟件；

優(yōu)選地，所述系統(tǒng)樹構(gòu)建模型為取AIC(Akaike Information Criterion)值最小的替代模型GTR+G；

優(yōu)選地，所述軟件2為RaxML(8.1.5)；

優(yōu)選地，所述構(gòu)建系統(tǒng)進化樹的方法為使用最大似然法構(gòu)建的系統(tǒng)進化樹，進化樹迭代次數(shù)為1000次。

本發(fā)明同時提供所述桑葚病原菌高通量鑒定及種屬分類方法在桑樹病原菌定量分析方面的應用。

本發(fā)明還提供所述桑葚病原菌高通量鑒定及種屬分類方法在分析微生物物種豐度方面的應用。

所述應用的方法為：使用分析軟件，計算核糖體DNA片段的平均測序深度，以此作為該物種的豐度值。

優(yōu)選地，所述分析軟件為bwa(0.7.12-r1039)與samtools(v1.2)。

本發(fā)明還提供所述桑葚病原菌高通量鑒定及種屬分類方法在計算物種間遺傳距離方面的應用。

附圖說明

圖1桑葚菌核病病果所檢測到的真菌微生物分類樹。

圖2驗證核糖體組裝結(jié)果的兩對引物與兩對通用引物ITS1/ITS4、ITS4/ITS5電泳圖。

圖3核糖體DNA序列組成與堿基GC比例分布圖。

圖4菌核病病原體Ciboria carunculoides 18S進化樹。

圖5菌核病病原體Ciboria carunculoides ITS進化樹。

圖6菌核病病原體Ciboria carunculoides28S進化樹。

圖7菌核病病原體rRNA全基因序列進化樹。

具體實施方式

下面結(jié)合具體實施例進一步說明本發(fā)明方法。下述實施例和附圖僅用于示例性說明，不能理解為對本發(fā)明的限制。除非特別說明，下述實施例中使用的試劑原料為常規(guī)市購或商業(yè)途徑獲得的生試劑原料，除非特別說明，下述實施例中使用的方法和設備為本領域常規(guī)使用的方法和設備

本發(fā)明以對桑葚菌核病病原菌的鑒定及種屬分類為例，根據(jù)本發(fā)明思想，本領域技術人員可以參照現(xiàn)有技術，對其他病癥的病原菌進行鑒定及種屬分類。

實施例1

在發(fā)病的桑園尋找到具有典型桑葚菌核病病征的桑果，收集，剪取桑病果中的患病區(qū)域，剪下的材料使用液氮充分研磨，總DNA的提取使用康為世紀公司的真菌DNA提取試劑盒，具體按照其操作說明進行，提取后的總DNA保存于-20℃。依照Illumina文庫構(gòu)建流程，將總DNA構(gòu)建為片段大小450bp的雙末端高通量測序文庫，使用Illumina Hiseq2500測序儀對構(gòu)建好的DNA文庫進行高通量測序，共測得9.03M對測序片段，測序讀長為雙末端125bp，總測序數(shù)據(jù)量2.26Gb。

由于DNA提取過程中包含桑樹的桑果材料，為減少測序數(shù)據(jù)中桑樹基因組數(shù)據(jù)對微生物序列組裝的影響，在進行微生物序列組裝前首先去除桑樹的基因組DNA序列。選取桑樹Morus notabilis全基因組序列(GCA_000414095.2)和葉綠體基因組序列(NC_027110.1)作為參考基因組序列，采用bwa(0.7.12-r1039)比對軟件進行數(shù)據(jù)比對分析。比對選擇mem比對算法，使用雙末端比對方法和軟件的默認參數(shù)，將測序數(shù)據(jù)與桑樹參考基因組進行比對，并將比對上參考基因組的測序片段判定為桑樹基因組序列。使用python編寫的計算機程序從fastq測序數(shù)據(jù)中去除桑樹測序數(shù)據(jù)，然后再進入微生物序列組裝。微生物序列的組裝使用MetaVelvet(v1.2.01)組裝軟件進行。

序列標簽注釋使用blastn(2.2.31+)序列比對分析軟件，組裝的序列標簽序列與NCBI的nt數(shù)據(jù)庫進行比對，blastn比對設定期望值<1e-20，根據(jù)比對結(jié)果對序列標簽進行注釋。核糖體DNA序列是細菌和真菌鑒定的重要的最常用的分子標記，因此物種分類鑒定和定量以核糖體DNA為主要分子標記。根據(jù)序列標簽注釋結(jié)果，選擇核糖體DNA序列作為微生物鑒定與定量分析依據(jù)。使用bwa(0.7.12-r1039)+samtools(v1.2)分析軟件，計算測序數(shù)據(jù)中核糖體DNA片段的平均測序深度，并以此作為該物種的豐度值。

結(jié)果顯示，共注釋228條rRNA序列標簽，數(shù)據(jù)包含4條細菌序列標簽，82條真菌序列標簽，分別對應3種細菌，3種真菌。通過對序列標簽注釋結(jié)果的查詢，共發(fā)現(xiàn)桑病果內(nèi)存在3個屬的真菌微生物，相對豐度最高的為Ciboria屬的真菌，見附圖1，其分子標簽測序深度為1115，相對豐度為98.50，所占比例最高，相關文獻報道該屬絕大部分為植物病原，大部分導致植物出現(xiàn)果實或種子僵化膨大等病征，與本發(fā)明研究的病征相符合。

真菌的核糖體DNA由18S區(qū)段，ITS1區(qū)段，5.8S區(qū)段，ITS2區(qū)段和28S區(qū)段組成，序列總長度約為5.5Kb。MetaVelvet(v1.2.01)初始組裝的序列標簽為斷裂的核糖體標簽，為得到完整的核糖體DNA序列，分析采用序列捕獲和從頭組裝策略，以組裝完整的核糖體DNA。選擇包含目標病原菌ITS序列的核糖體DNA序列為參考序列，采用bwa(0.7.12-r1039)軟件進行0mismatch和0gap比對，根據(jù)比對結(jié)果從測序數(shù)據(jù)中獲取雙末端測序片段，進一步采用MetaVelvet(v1.2.01)組裝軟件對序列進行組裝和延伸，經(jīng)多個循環(huán)操作，獲取完整的核糖體DNA序列。

根據(jù)組裝完整的核糖體DNA序列，采用真菌通用引物ITS1和ITS4、ITS4和ITS5(White，1990)，和基于本研究組裝獲得的完整核糖體DNA序列設計了PCR驗證引物見表1，然后進行PCR擴增和Sanger測序方法，驗證組裝的正確性。通用引物序列見表2所示，PCR反應體系見表3所示。

表1.PCR驗證引物序列表

表2真菌ITS區(qū)域通用引物序列

表3 PCR反應體系(20μL)

真菌通用引物ITS1和ITS4引物組PCR程序：預變性94℃5min；94℃30s，55℃30s，72℃1min，30個循環(huán)；72℃10min。

真菌通用引物ITS4和ITS5引物組PCR程序：預變性94℃5min；94℃30s，55℃30s，72℃1min，30個循環(huán)；72℃10min。

驗證核糖體組裝結(jié)果的P1-F/R引物組PCR反應體系(20μL)PCR程序：94℃5min；94℃1min，56℃1min，72℃4min，30個循環(huán)；72℃10min。

驗證核糖體組裝結(jié)果的P2-F/R引物組PCR反應體系(20μL)PCR程序：94℃5min；94℃1min，56℃1min，72℃2min30s，30個循環(huán)；72℃10min。

取3μLPCR擴增產(chǎn)物用1.2％的瓊脂糖凝膠(EB染色)電泳檢測。通過瓊脂糖凝膠電泳回收大小對應的PCR產(chǎn)物片段，見附圖1所示，圖1中，M：Takara DL5000Marker；1.真菌通用引物ITS1和ITS4；2.真菌通用引物ITS4和ITS5；3.驗證核糖體組裝結(jié)果引物P1-F和P1-R；4.驗證核糖體組裝結(jié)果引物P2-F和P2-R。

PCR驗證結(jié)果見附圖2，從圖中可以看出，本次樣品較純凈，泳道1和2的結(jié)果說明用真菌通用引物ITS1和ITS4、ITS4和ITS5擴增的模板DNA是屬于真菌的DNA；泳道3和4的結(jié)果說明驗證核糖體組裝結(jié)果的P1引物組和P2引物組均擴增到相應的目的條帶，進一步測序后的結(jié)果與組裝的結(jié)果高度一致。

結(jié)果顯示，根據(jù)高通量測序數(shù)據(jù)組裝與實驗驗證結(jié)果，組裝得到完整核糖體DNA序列，見附圖3，長度6202bp(GC比例52.98％)。序列包含18S區(qū)域，ITS1區(qū)域，5.8S區(qū)域，ITS2區(qū)域，28S區(qū)域。其中18S區(qū)域長度1782bp(GC比例50.05％)，ITS1區(qū)域長度101bp(GC比例43.56％)，5.8S區(qū)域長度151bp(GC比例46.36％)，ITS2區(qū)域長161bp(GC比例56.52％)，28S區(qū)域長4007bp(GC比例54.63％)。18S區(qū)域和28S區(qū)域長度較大，占序列總長度的96.37％；而在GC比例上，ITS2區(qū)域(56.52％)的GC比例則明顯高于其他區(qū)域。

構(gòu)建了對應Ciboria carunculoides的rRNA分子系統(tǒng)樹模型，菌核病病原體Ciboria carunculoides 18S進化樹、菌核病病原體Ciboria carunculoides ITS進化樹、菌核病病原體Ciboria carunculoides 28S進化樹、菌核病病原體Ciboria carunculoides的rRNA全基因序列進化樹分別見附圖4～附圖7所示。

使用blastn(2.2.31+)比對軟件將核糖體DNA與nt數(shù)據(jù)庫進行比對，對核糖體DNA序列進行注釋，共注釋出28S區(qū)域，ITS1區(qū)域，5.8S區(qū)域，ITS2區(qū)域和28S區(qū)域。使用python語言編寫的計算機程序?qū)π蛄械腉C比例進行分析，計算各區(qū)域的平均GC比例。同時程序以100bp為計算窗口，以10bp為步進，計算核糖體DNA的GC比例特征。

SEQUENCE LISTING

<110> 華南農(nóng)業(yè)大學

<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