花桃流膠病桃膠真菌多樣性及基因差異表達分析

徐 曼，王燕，張文斗，吳超廣

（常德職業(yè)技術(shù)學(xué)院，湖南常德 415000）

桃屬薔薇科（Rosaceae）桃屬（Prunus Linn.）植物，原產(chǎn)于我國黃河上游海拔1200～2000m 的高原地帶，是我國人民普遍喜愛栽培的最古老果樹之一[1]。迄今為止，桃樹在我國已有4000 多年的栽培歷史?；ㄌ?，一般稱作“觀賞桃”，有時“桃花”即指此意[2]。但南方高溫高濕的氣候易使花桃感染流膠病，流膠病在一定程度上制約了桃花的發(fā)展和應(yīng)用[3]。流膠病又名瘤皮病、疣皮病，是在核果類果樹上普遍發(fā)生的一種病害，其分布十分廣泛，且以多雨地區(qū)發(fā)生更為嚴(yán)重。流膠病的發(fā)生會導(dǎo)致桃樹樹勢衰弱，產(chǎn)量銳減，壽命縮短，甚至樹體死亡，是桃種植業(yè)中的一大障礙[4]。本研究以二代高通量測序為基礎(chǔ)的18S rDNA 技術(shù)對編碼原核核糖體小亞基rRNA 的DNA 序列進行測序，分析桃膠中的真菌物種相對豐度，并采用轉(zhuǎn)錄組測序技術(shù)對真菌感染后的基因差異表達進行對比[5]，探究花桃桃膠真菌的物種注釋、分類學(xué)分析、組間差異顯著性分析，然后通過轉(zhuǎn)錄組測序技術(shù)，對花桃桃膠病感染位置的基因表達進行分析，為花桃抗桃膠病的分子基礎(chǔ)及品種創(chuàng)新應(yīng)用研究提供科學(xué)依據(jù)。

1 材料及方法

1.1 材料采集

采樣地為湖南省常德職業(yè)技術(shù)學(xué)院花桃種質(zhì)資源庫。該研究區(qū)保育了桃花種質(zhì)資源120 多種，且處于南方高溫高濕地區(qū)，為流膠病多發(fā)區(qū)。本研究以壽星（HSS）易感病品種和合歡二色（HHR）高抗病品種為研究對象，取樹干上的桃膠作為18S rRNA 測序樣品，感病位置樹皮作為mRNA 測序樣品，重復(fù)3 次[6]。

1.2 桃膠真菌的高通量分析

1.2.1 建庫測序。使用Illumina 公司TruSeq DNA PCR-Free Sample Preparation Kit 建庫試劑盒（Illumina,USA）進行文庫的構(gòu)建，構(gòu)建好的文庫經(jīng)過Qubit 定量和文庫檢測，合格后，使用NovaSeq 6000 PE250 進行上機測序[7]。

1.2.2 測序數(shù)據(jù)分析。使用ANCOM、ANOVA、Kruskal Wallis、LEfSe 和DESeq2 等方法鑒定分組和樣本間豐度有差異的真菌。

1.3 花桃感染位置mRNA 高通量分析

1.3.1 建庫測序。建庫使用llumina 的NEBNext UltraTMRNA Library Prep Kit 建庫試劑盒，純化后的雙鏈cDNA 經(jīng)過末端修復(fù)，最終獲得文庫。

1.3.2 測序數(shù)據(jù)分析。使用DESeq2 軟件進行2 個比較組合之間的差異表達分析（每個組2 個生物學(xué)重復(fù)）。DESeq2 提供了統(tǒng)計程序，用于使用基于負(fù)二項式分布的模型來確定數(shù)字基因表達數(shù)據(jù)中的差異表達。

2 結(jié)果與分析

2.1 OTU 物種組成分析

使用Qiime2 軟件對所有樣品的全部原始序列進行質(zhì)量控制，去噪，拼接，并且去嵌合體，形成OTU，選取OTU 的代表性序列，與18S 數(shù)據(jù)庫進行比對，獲得物種注釋信息，得到575 個OTU。通過繪制韋恩圖（Venn diagram），分析不同樣品組之間特有或共有的OTU，發(fā)現(xiàn)樣品HHR 的OTU 為382，HSS 的OTU 為278。如圖1 所示，共有的OTU 為85，HHR 特有的OTU 為297，HSS 特有的OTU 為193。因而可知，隨著品種的抗性差異，桃膠的真菌物種組成也存在著較大的差異。

2.2 真菌群落系統(tǒng)進化分析及差異分布

HHR 和HSS 2 個樣品進化樹的繪制由R 語言ggtree 包完成。選取關(guān)注的OTU 代表序列進行系統(tǒng)進化分析（每個屬選取一個豐度最高的OTU 作為代表OTU，之后再選取豐度最高的前50 個屬）繪制進化關(guān)系樹圖，同時結(jié)合OTU 在各個分組的絕對豐度進行熱圖可視化展示（圖2）。其中在門水平上可分為擔(dān)子菌門Basidiomycota 和子囊菌門Ascomycota 2 個分類，其擔(dān)子菌門包含了Colletotrichum、Nectria、arasarcopodium、Trichoderma、ljuhya、Cephalotheca、Microascus、Discosia、Cephalotrichum、Microdochium、Phomopsis、Cytospora、Penicillium、Aspergilus、Phaeosphaeria、Cuyularia、Phom-a、Sclerococcum、Saccharomyces、Candida 20 個屬，其中Colletotrichum、Nectria、arasarcopodium、Trichoderma、lju-hya、Cephalotheca、Microascus、Microdochium、Phomopsis、Phaeosphaeria、Cuyularia、Saccharomyces中的屬間豐度比較，HSS遠(yuǎn)高于HHR，而在Disco sia、Cephalotrichum、Cytospora、Penicillium、Aspergilus、Phom-a、Sclerococcum、Candida 的比較中，HSS 低于HHR；子囊菌門包含Banno a、Hygrocybe、Clavulina、Thanatephorus、Hebe loma、Psathyrella、Perenniporia、Russula、Serendipita、Sebacin a 10個屬，Banno a、Clavulina、Thanatephorus、Hebe loma、Serendipita、Sebacin a 中的屬間豐度比較，HSS 遠(yuǎn)高于HHR，而Hygrocybe、Psathyrella、Perenniporia、Russula 的比較中，HSS低于HHR。因此，花桃桃膠病的發(fā)生可能是多個真菌協(xié)同侵染的結(jié)果。

2.3 花桃基因的差異表達及注釋

通過TPM（Transcript Per Million）或FPKM（Fragment Per Kilo base Million reads）值來表示某個基因在某個樣本中的表達量，該值與mRNA 的表達量具有正相關(guān)性，而且是標(biāo)準(zhǔn)化處理后的值。如圖3 所示，對HSS 和HHR 2 個組進行比較，以火山圖的形式進行展示，對火山圖中的18502 個基因進行了比較，紅色標(biāo)記的基因表現(xiàn)出顯著性的差異（P＜0.05，log2FC＞=1），其中具有顯著差異的2129 個基因表現(xiàn)出HRR 比HSS高，724 個基因表現(xiàn)為HRR 比HSS 低。因此，花桃會因品種抗性間的差異表現(xiàn)出不同的基因表達模式。

然后，又從KEGG 物種基因組官網(wǎng)中收集到了桃樹的編碼蛋白的氨基酸序列作為參考系列，通過InterProScan 軟件對整個基因組上的所有蛋白序列做GO 詞條注釋，GO 詞條可分為生物學(xué)途徑（Biological Process）、細(xì)胞組件（Cellular Component）和分子功能（Molecule Function）3 個方面，選擇基因?qū)?yīng)的最長的蛋白序列，來對這個蛋白做GO 詞條注釋（圖5），作為這個基因的注釋結(jié)果。最后通過對上游顯著差異表達基因的GO 詞條注釋，各選出擁有顯著差異基因個數(shù)最多的前10 種GO 詞條，用R 語言作堆疊型條形統(tǒng)計圖（圖5），粉紅、紅色表示高表達，藍色表示低表達。其中生物學(xué)途徑中，蛋白質(zhì)磷酸化1（Protein Phosphorylation 1）、DNA 為模板的轉(zhuǎn)錄調(diào)節(jié)因子（Regulation Of Transcription,Dna-Templated）、防御感知生物學(xué)途徑（Defense Response Biological Process）、跨膜轉(zhuǎn)運因子（Transmembrane Transport）、碳水化合物代謝過程（Carbohydrate Metabolic Process）、信號轉(zhuǎn)導(dǎo)（Signal Transduction）、蛋白水解（Proteolysis）、脂質(zhì)代謝途徑（Lipid Metabolic Process）、蛋白泛素化（Protein Ubiquitination）、細(xì)胞葡聚糖代謝過程（Cellular Glucan Metabolic Process）等生物學(xué)途徑相關(guān)的基因表現(xiàn)出了明顯差異；而細(xì)胞元件中，則在核（Nucleus）、質(zhì)外體（Apoplast）細(xì)胞質(zhì)（Cytoplasm）、細(xì)胞壁（Cell Wall）、細(xì)胞間隙（Extracellular Region）、微管（Microtubule）、胞外（Exocyst）、核小體（Nucleosome）等相關(guān)基因表現(xiàn)出了顯著差異；最后在分子功能方面，蛋白結(jié)合（Protein Binding）、Atp 集合（Atp Binding）、蛋白激酶活動（Protein Kinase Activity）、Dna 結(jié)合（Dna Binding）、氧化還原酶活性（oxidoreductase activity）、DNA 結(jié)合轉(zhuǎn)錄因子活動（DNA-binding transcription factor activity）、ADP 結(jié)合（ADP binding）、催化活性（catalytic activity）、鋅離子結(jié)合（zinc ion binding）、血紅素結(jié)合（heme binding）等相關(guān)基因表現(xiàn)了明顯的差異。

2.4 差異表達基因的KEGG 代謝通路分析

KEGG 代謝通路數(shù)據(jù)庫中，有很多前人歸納的生物的代謝通路，一個代謝通路包含了一定數(shù)量的基因，為了上文中注釋的差異基因是否顯著地對應(yīng)著一些功能分類的基因的集合，實驗使用R 語言clusterProfiler包，通過富集分析（over-representation analysis），選取KEGG 通路富集分析結(jié)果中的最顯著的通路作圖，得到KEGG 通路dot-plot 氣泡圖（圖6），其中縱坐標(biāo)表示pathway 的名稱，橫坐標(biāo)GeneRatio 的含義，Count 表示差異基因中有多少個基因落到了某一個通路中，差異基因個數(shù)越多，則該通路圓點越大。如圖6 所示，具有顯著差異的KEGG 通路包括植物與病原菌互作、內(nèi)質(zhì)網(wǎng)蛋白加工、淀粉和蔗糖代謝、亞麻酸代謝、半乳糖代謝、光合作用，其中植物病原菌互作所含差異基因最多，這與本文研究的花桃抗流膠病樣本結(jié)果是一致的。

2.5 植物病原菌互作相關(guān)基因分析

根據(jù)差異基因的KEGG 代謝通路分析可知，HSS和HRR 差異表達基因KEGG 代謝通路注釋的差異基因最多，因而對該通路的基因進一步分析。由表1 可知，發(fā)現(xiàn)2 個樣本中被注釋為植物微生物互作通路（Plant-pathogen interaction）的基因共有59 個。如圖7所示，該代謝通路中的59 個基因分別參與到了真菌模式分子（Fungal Pamp）、細(xì)菌鞭毛蛋白（Bacterial Flagellin）、細(xì)菌EF-Tu 分子（Bacterial EF-Tu）、真菌效應(yīng)因子（fungal effector）、細(xì)菌分泌系統(tǒng)（Bac teria Secretion System）的識別，其中CDPK、Rboh、CaMCM、WRKY、Pti1/5、HSP90、PBS1、RPS2、KCS1/10 等注釋的基因有著顯著的升高，而CNGCs、RPM1、Rd19 等基因卻顯著下降，也從側(cè)面驗證了花桃桃膠病的發(fā)生與真菌有關(guān)。

表1 植物與病原菌互作通路基因富集

3 結(jié)論

通過對2 個不同花桃品種桃膠真菌18S 測序和感染位置轉(zhuǎn)錄組聯(lián)合分析得出：一是花桃桃膠病的發(fā)生與擔(dān)子菌門和子囊菌門不同種屬的真菌互作有關(guān)，且樣品的抗性差異導(dǎo)致了真菌豐度的差異；二是通過轉(zhuǎn)錄組測序，檢測59 個植物微生物互作通路的基因，即植物防御基因，這些基因參與到了真菌模式分子（Fungal Pamp）、細(xì)菌鞭毛蛋白（Bacterial Flagellin）、細(xì)菌EF-Tu 分子（bacterial EF-Tu）、真菌效應(yīng)因子（fungal effector）、細(xì)菌分泌系統(tǒng)（Bacteria Secretion System）的識別，其中CDPK、Rboh、CaMCM、WRKY、Pti1/5、HSP90、PBS1、RPS2、KCS1/10 等基因的表達有利于抵抗桃膠病的發(fā)生。