花桃流胶病桃胶真菌多样性及基因差异表达分析

徐 曼，王燕，张文斗，吴超广

（常德职业技术学院，湖南常德 415000）

桃属蔷薇科（Rosaceae）桃属（Prunus Linn.）植物，原产于我国黄河上游海拔1200～2000m 的高原地带，是我国人民普遍喜爱栽培的最古老果树之一[1]。迄今为止，桃树在我国已有4000 多年的栽培历史。花桃，一般称作“观赏桃”，有时“桃花”即指此意[2]。但南方高温高湿的气候易使花桃感染流胶病，流胶病在一定程度上制约了桃花的发展和应用[3]。流胶病又名瘤皮病、疣皮病，是在核果类果树上普遍发生的一种病害，其分布十分广泛，且以多雨地区发生更为严重。流胶病的发生会导致桃树树势衰弱，产量锐减，寿命缩短，甚至树体死亡，是桃种植业中的一大障碍[4]。本研究以二代高通量测序为基础的18S rDNA 技术对编码原核核糖体小亚基rRNA 的DNA 序列进行测序，分析桃胶中的真菌物种相对丰度，并采用转录组测序技术对真菌感染后的基因差异表达进行对比[5]，探究花桃桃胶真菌的物种注释、分类学分析、组间差异显著性分析，然后通过转录组测序技术，对花桃桃胶病感染位置的基因表达进行分析，为花桃抗桃胶病的分子基础及品种创新应用研究提供科学依据。

1 材料及方法

1.1 材料采集

采样地为湖南省常德职业技术学院花桃种质资源库。该研究区保育了桃花种质资源120 多种，且处于南方高温高湿地区，为流胶病多发区。本研究以寿星（HSS）易感病品种和合欢二色（HHR）高抗病品种为研究对象，取树干上的桃胶作为18S rRNA 测序样品，感病位置树皮作为mRNA 测序样品，重复3 次[6]。

1.2 桃胶真菌的高通量分析

1.2.1 建库测序。使用Illumina 公司TruSeq DNA PCR-Free Sample Preparation Kit 建库试剂盒（Illumina,USA）进行文库的构建，构建好的文库经过Qubit 定量和文库检测，合格后，使用NovaSeq 6000 PE250 进行上机测序[7]。

1.2.2 测序数据分析。使用ANCOM、ANOVA、Kruskal Wallis、LEfSe 和DESeq2 等方法鉴定分组和样本间丰度有差异的真菌。

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

1.3.1 建库测序。建库使用llumina 的NEBNext UltraTMRNA Library Prep Kit 建库试剂盒，纯化后的双链cDNA 经过末端修复，最终获得文库。

1.3.2 测序数据分析。使用DESeq2 软件进行2 个比较组合之间的差异表达分析（每个组2 个生物学重复）。DESeq2 提供了统计程序，用于使用基于负二项式分布的模型来确定数字基因表达数据中的差异表达。

2 结果与分析

2.1 OTU 物种组成分析

使用Qiime2 软件对所有样品的全部原始序列进行质量控制，去噪，拼接，并且去嵌合体，形成OTU，选取OTU 的代表性序列，与18S 数据库进行比对，获得物种注释信息，得到575 个OTU。通过绘制韦恩图（Venn diagram），分析不同样品组之间特有或共有的OTU，发现样品HHR 的OTU 为382，HSS 的OTU 为278。如图1 所示，共有的OTU 为85，HHR 特有的OTU 为297，HSS 特有的OTU 为193。因而可知，随着品种的抗性差异，桃胶的真菌物种组成也存在着较大的差异。

2.2 真菌群落系统进化分析及差异分布

HHR 和HSS 2 个样品进化树的绘制由R 语言ggtree 包完成。选取关注的OTU 代表序列进行系统进化分析（每个属选取一个丰度最高的OTU 作为代表OTU，之后再选取丰度最高的前50 个属）绘制进化关系树图，同时结合OTU 在各个分组的绝对丰度进行热图可视化展示（图2）。其中在门水平上可分为担子菌门Basidiomycota 和子囊菌门Ascomycota 2 个分类，其担子菌门包含了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远高于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 远高于HHR，而Hygrocybe、Psathyrella、Perenniporia、Russula 的比较中，HSS低于HHR。因此，花桃桃胶病的发生可能是多个真菌协同侵染的结果。

2.3 花桃基因的差异表达及注释

通过TPM（Transcript Per Million）或FPKM（Fragment Per Kilo base Million reads）值来表示某个基因在某个样本中的表达量，该值与mRNA 的表达量具有正相关性，而且是标准化处理后的值。如图3 所示，对HSS 和HHR 2 个组进行比较，以火山图的形式进行展示，对火山图中的18502 个基因进行了比较，红色标记的基因表现出显著性的差异（P＜0.05，log2FC＞=1），其中具有显著差异的2129 个基因表现出HRR 比HSS高，724 个基因表现为HRR 比HSS 低。因此，花桃会因品种抗性间的差异表现出不同的基因表达模式。

然后，又从KEGG 物种基因组官网中收集到了桃树的编码蛋白的氨基酸序列作为参考系列，通过InterProScan 软件对整个基因组上的所有蛋白序列做GO 词条注释，GO 词条可分为生物学途径（Biological Process）、细胞组件（Cellular Component）和分子功能（Molecule Function）3 个方面，选择基因对应的最长的蛋白序列，来对这个蛋白做GO 词条注释（图5），作为这个基因的注释结果。最后通过对上游显著差异表达基因的GO 词条注释，各选出拥有显著差异基因个数最多的前10 种GO 词条，用R 语言作堆叠型条形统计图（图5），粉红、红色表示高表达，蓝色表示低表达。其中生物学途径中，蛋白质磷酸化1（Protein Phosphorylation 1）、DNA 为模板的转录调节因子（Regulation Of Transcription,Dna-Templated）、防御感知生物学途径（Defense Response Biological Process）、跨膜转运因子（Transmembrane Transport）、碳水化合物代谢过程（Carbohydrate Metabolic Process）、信号转导（Signal Transduction）、蛋白水解（Proteolysis）、脂质代谢途径（Lipid Metabolic Process）、蛋白泛素化（Protein Ubiquitination）、细胞葡聚糖代谢过程（Cellular Glucan Metabolic Process）等生物学途径相关的基因表现出了明显差异；而细胞元件中，则在核（Nucleus）、质外体（Apoplast）细胞质（Cytoplasm）、细胞壁（Cell Wall）、细胞间隙（Extracellular Region）、微管（Microtubule）、胞外（Exocyst）、核小体（Nucleosome）等相关基因表现出了显著差异；最后在分子功能方面，蛋白结合（Protein Binding）、Atp 集合（Atp Binding）、蛋白激酶活动（Protein Kinase Activity）、Dna 结合（Dna Binding）、氧化还原酶活性（oxidoreductase activity）、DNA 结合转录因子活动（DNA-binding transcription factor activity）、ADP 结合（ADP binding）、催化活性（catalytic activity）、锌离子结合（zinc ion binding）、血红素结合（heme binding）等相关基因表现了明显的差异。

2.4 差异表达基因的KEGG 代谢通路分析

KEGG 代谢通路数据库中，有很多前人归纳的生物的代谢通路，一个代谢通路包含了一定数量的基因，为了上文中注释的差异基因是否显著地对应着一些功能分类的基因的集合，实验使用R 语言clusterProfiler包，通过富集分析（over-representation analysis），选取KEGG 通路富集分析结果中的最显著的通路作图，得到KEGG 通路dot-plot 气泡图（图6），其中纵坐标表示pathway 的名称，横坐标GeneRatio 的含义，Count 表示差异基因中有多少个基因落到了某一个通路中，差异基因个数越多，则该通路圆点越大。如图6 所示，具有显著差异的KEGG 通路包括植物与病原菌互作、内质网蛋白加工、淀粉和蔗糖代谢、亚麻酸代谢、半乳糖代谢、光合作用，其中植物病原菌互作所含差异基因最多，这与本文研究的花桃抗流胶病样本结果是一致的。

2.5 植物病原菌互作相关基因分析

根据差异基因的KEGG 代谢通路分析可知，HSS和HRR 差异表达基因KEGG 代谢通路注释的差异基因最多，因而对该通路的基因进一步分析。由表1 可知，发现2 个样本中被注释为植物微生物互作通路（Plant-pathogen interaction）的基因共有59 个。如图7所示，该代谢通路中的59 个基因分别参与到了真菌模式分子（Fungal Pamp）、细菌鞭毛蛋白（Bacterial Flagellin）、细菌EF-Tu 分子（Bacterial EF-Tu）、真菌效应因子（fungal effector）、细菌分泌系统（Bac teria Secretion System）的识别，其中CDPK、Rboh、CaMCM、WRKY、Pti1/5、HSP90、PBS1、RPS2、KCS1/10 等注释的基因有着显著的升高，而CNGCs、RPM1、Rd19 等基因却显著下降，也从侧面验证了花桃桃胶病的发生与真菌有关。

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

3 结论

通过对2 个不同花桃品种桃胶真菌18S 测序和感染位置转录组联合分析得出：一是花桃桃胶病的发生与担子菌门和子囊菌门不同种属的真菌互作有关，且样品的抗性差异导致了真菌丰度的差异；二是通过转录组测序，检测59 个植物微生物互作通路的基因，即植物防御基因，这些基因参与到了真菌模式分子（Fungal Pamp）、细菌鞭毛蛋白（Bacterial Flagellin）、细菌EF-Tu 分子（bacterial EF-Tu）、真菌效应因子（fungal effector）、细菌分泌系统（Bacteria Secretion System）的识别，其中CDPK、Rboh、CaMCM、WRKY、Pti1/5、HSP90、PBS1、RPS2、KCS1/10 等基因的表达有利于抵抗桃胶病的发生。