枣疯病导致的花异常发育相关基因鉴定及表达分析

刘伟，霍辰思，王国平，王定仙，樊新萍

（山西农业大学果树研究所，山西 太谷 030801）

枣（Ziziphus jujubaMill.）属于鼠李科枣属，具有较强的抗逆性及适应性，其果实既可以鲜食亦可作为加工食品的原料，因其优良的营养价值而广受欢迎［1］。枣是我国特有的经济林及果树树种，也是我国第一个完成全基因组测序的木本干果树种［2］，具有较高的研究价值及较广的产业发展前景。

枣疯病作为我国枣树的严重病害之一，近些年给整个枣产业造成了巨大的损失。枣树染病后常表现出花变叶、丛枝和小叶症状，最终导致枣树绝产，甚至树体死亡。该病害在我国南北方各枣区均有发生，对生产影响极大，严重威胁枣产业的健康发展。枣疯病属于一种植原体病害，由于植原体具有韧皮部专性寄生的特性，目前还无法进行人工离体培养，因此使得该病害致病机理方面的研究十分困难。前人针对枣疯病的研究主要涉及不同品种的抗病性、病原物的提取与鉴定、染病植株的生理生化指标的改变等。郝新迪等人探索了不同枣树品种枣疯病的抗性差异以及病原物在染病植株体内迁移特性差异，初步确认了不同枣树品种对枣疯病植原体具有响应差异［3］。肖京等人以骏枣不同株系为试材，筛选得到了具有较高的枣疯病抗性的优良株系［4］。这些为抗病品种的选育奠定了坚实基础。针对植原体病原物的检测鉴定，目前已开发出常规 PCR 检测［5］、多重 PCR 检测［6］、荧光定量 PCR 检测［7，8］以及 LAMP 可视化检测［9］等多种检测方法。此外，研究表明，感病枣树叶片叶绿素含量显著降低，光合系统及保护酶系统均受到显著影响［10，11］。

近年来，测序技术的发展对植物病害研究领域产生了极大影响，促进了人们对病原物特点及致病分子机制的深入研究，提高了对病原物、植物与昆虫互作关系的更深层次的理解。一些枣疯病植原体相关基因序列得到进一步解析，比如可用于植原体亚组分类的延伸因子tuf基因、运输蛋白基因（secY）以及核糖体蛋白基因（rp）等［12~15］。枣疯病相关的部分致病因子被发现并得到初步验证。倪静发现了具有核定位的Zaofeng 8 致病因子并明确了其致病的表型［16］。屈晓逸基于高通量测序技术及生物信息学分析预测得到11 种枣疯病植原体的效应蛋白［17］。一些响应枣疯病植原体侵染的基因也被发现，比如丝氨酸/苏氨酸蛋白激酶家族 成 员ZjMAPKKKs［18］、谷 胱 甘 肽 S 移 换 酶 基 因ZjGSTU1［19］、AP2/ERF转录因子［20］、WRKY转录因子［21］、TCP转录因子［22］等。枣疯病致病机理也被进一步挖掘，Fan 等利用转录组分析枣疯病株，获得了植原体侵染枣树后叶片及花器官差异表达的基因［23］。Ye 等利用 iTRAQ 蛋白质组学与转录组测序分析枣疯病植原体侵染枣树后的多重调控水平［24］。但是枣疯病导致的枣花异常发育以及小叶、丛枝等异常性状的分子基础仍不清楚。

本研究基于转录组高通量测序，分析枣疯病染病植株异常发育花器官差异表达基因，进一步挖掘成花诱导与发育相关的关键基因并进行荧光定量PCR 表达验证。初步探讨枣花异常发育相关基因的表达调控规律，旨在为进一步挖掘植原体导致的器官异常发育的分子机制提供理论依据。

1 材料和方法

1.1 试验材料

本研究采集的植物样本均来源于山西农业大学果树研究所国家枣种质资源圃，供试枣品种为“伏脆蜜”。于2018 年6 月5 日在田间采集枣健康植株及枣疯病染病植株花器官样本，各取3 个重复，装入无菌管后，迅速投入液氮速冻，放入−80℃超低温冰箱保存备用。

1.2 试验方法

1.2.1 RNA 提取及质量检测

采用RNAprep Pure 多糖多酚植物总RNA 提取试剂盒（天根生化科技有限公司）提取各样本总RNA，利用DNaseⅠ去除基因组DNA。使用Agilent 2100 Bioanalyzer 生物分析仪检测RNA 的浓度以及完整性，使用Eppendorf BioPhotometer D30核酸蛋白测定仪检测RNA 纯度。

1.2.2 转录组测序与数据分析

转录组测序由上海派森诺生物科技有限公司完成。样品经过RNA 抽提、纯化、建库之后，采用第二代测序技术（Next-Generation Sequencing，NGS），基于 Illumina HiSeq 测序平台，对这些文库进行双末端（Paired-end，PE）测序。测序获得的原始数据经过滤得到高质量序列（Clean Data），将其与枣的参考基因组比对［2］。利用Geneontology数 据 库（http：//geneontology. org）［25，26］和 KEGG数据库（https：//www.kegg.jp/kegg/pathway.html）［27］进行基因功能注释及富集分析，用 R 软件作图。采用TBtools 软件绘制差异基因表达热图［28］。

1.2.3 实时荧光定量qRT-PCR 检测

利用与转录组测序相同的样本进行实时荧光定量PCR 检测，3 次生物学重复。将提取的RNA参照Promega 反转录试剂盒说明书制备成cDNA。根据筛选得到的差异表达基因序列设计并合成qRT-PCR 引物，以枣Actin-depolymerizing factor 1（ACT1）基因作为内参基因进行qRT-PCR 检测［29］，具体引物序列信息见表 1。qRT-PCR 反应液按照 GoTaq qPCR Master Mix（Promega）试剂盒说明书进行配制，反应液总体积为20 μL，反应体系中包含GoTaq qPCR Master Mix（2×）10 μL、上下游引物（10 μmol·L-1）各 0.4 μL、cDNA 模板 2 μL、Nuclease-Free Water 7.2 μL。利用罗氏Light Cycler 96 实时荧光定量PCR 仪进行qRTPCR 检测。qRT-PCR 反应程序：95 ℃预变性 2 min；95 ℃变性 15 s，60 ℃退火延伸 1 min，40 次循环。

表1 枣花器官生长发育相关基因qRT-PCR 引物Table 1 qRT-PCR primers of genes involved in floral organ development of jujube

1.3 统计分析

1.3.1 转录组差异表达基因分析

使用HTSeq 统计比对到每一个基因上Read Count 值，作为基因的原始表达量，采用RPKM 对表达量进行标准化。采用DESeq 对枣疯病染病植株与健康植株样品基因表达进行差异分析，筛选差异表达基因条件为：表达差异倍数|log2 Fold-Change|＞1 ，显著性P-value＜0.05 。

1.3.2 实时荧光定量PCR 数据统计分析

获取qRT-PCR 检测得到的Cq 值数据，采用2-ΔΔCt法计算基因的相对表达量及表达量的标准化［30］。利用GraphPad Prism 软件进行作图及表达差异显著性分析（T 检验，P＜0.01）。

2 结果与分析

2.1 测序数据组装与质量分析

为深入探索枣疯病在转录水平上对枣树花器官的影响，采集枣疯病染病植株及健康植株花器官样本，利用Illumina HiSeq 测序平台进行了转录组测序，共获得 37.21 Gb Clean Data，Q20 值平均为96.34%，Q30 值平均为91.75%。样品测序获得原始 reads 数为 39 804 662~43 115 716 条，通过对原始reads 进行数据过滤，将有效reads 比对到枣参考基因组，平均有78.51%有效reads 能够比对到枣基因组上。其中比对到多个位置的序列总数平均为5 520 171，占总序列的平均百分比为17.06%；只比对到一个位置的序列总数平均为26 836 114，占总序列的平均百分比为82.94%。具体测序数据信息如表2 所示。

表2 样品转录组测序数据统计及参考基因组比对情况Table 2 Samples of the transcriptome sequencing data statistics and comparison with the reference genome

对测序的6 个样品数据进行主成分分析（Principal Component Analysis，PCA），染病植株样本和健康植株样本被清晰的划分为2 类（图1A），同一样品类型的不同样品之间紧密排列，说明样品间重复性良好。PCA 分析结果表明枣花器官转录组测序结果整体评估优良，可用于后续分析。同时，进一步对样品间基因表达水平相关性进行分析。相关性分析是检验实验可靠性和样本选择是否合理的重要指标，也是进行差异表达分析的基础。由图1B 可知，生物学重复的样本之间相关系数均大于0.9，表现为极强的相关性，而不同分组的样本间的相关系数均低于0.8，表示不同分组样品之间的相关性较低。这进一步证明转录组测序结果的可靠性，可进行后续差异表达分析。

图1 转录组测序数据评估及差异基因分析Fig.1 The transcriptome sequencing data assessment and DEGs analysis

2.2 差异表达基因鉴定

采用DESeq 对染病植株和健康植株的花器官样本基因表达进行差异分析，结果如图1C 所示。枣疯病染病植株花器官与健康植株样本相比，有4 148 个基因存在表达差异，其中有2 305 个基因在染病植株中表达上调，占差异表达基因总量的55.57%。有1 843 个基因在染病植株中表达下调，占差异表达基因总量的44.43%。

2.3 差异表达基因的GO 分类及显著性富集分析

首先将比对到基因组的所有基因映射到Gene Ontology 数据库的各个条目，计算每个条目的差异基因数量，以整个基因组为背景，采用超几何分布计算差异基因显著富集的条目。其中14 472 个基因能够被GO 注释，其中差异表达基因被注释的数量为2 188 个，被注释的总Term 数为1 954 个。被注释的基因在GO 数据库中被分为细胞组件、分子功能及生物过程3 类，其中差异表达基因在细胞组分大类被注释到196 个GO Term，分子功能大类被注释到553 个GO Term，生物过程大类被注释到1 205 个GO Term。由此可见，枣疯病染病植株的许多生物反应过程发生了改变。

通过对差异表达基因的GO 条目进行富集分析，获得每个大类中富集程度较高的GO 条目。由图2 可知，细胞组件类别中，富集程度最高的条目（前10 条）依次是：类囊体部分（GO：0044436）、光合膜（GO：0034357）、类囊体（GO：0009579）、类囊体膜（GO：0042651）、光系统（GO：0009521）、叶绿体类囊体膜（GO：0009535）、质体类囊体膜（GO：0055035）、胞外区（GO：0005576）、叶绿体类囊体（GO：0009534）、质体类囊体（GO：0031976）。由此可见枣疯病染病植株的光合系统及生物膜系统受到较大影响。分子功能类别中，富集程度最高的条目（前10 条）依次是：氧化还原酶活性（GO：0016491）、叶绿素结合（GO：0016168）、四吡咯色素结合（GO：0046906）、催化活性（GO：0003824）、单加氧酶活性（GO：0004497）、木葡聚糖木葡糖基转移酶活性（GO：0016762）、二级主动跨膜转运蛋白活性（GO：0015291）、氧化还原酶活性（GO：0016679）、果胶裂解酶活性（GO：0030570）、碳-氧裂合酶活性（GO：0016837）。生物过程类别中，富集程度最高的条目（前10 条）依次是：光合作用（GO：0015979）、氧化还原过程（GO：0055114）、光合作用：光反应（GO：0019684）、光合作用：光捕获（GO：0009765）、外 部 封 装 结 构 组 织（GO：0045229）、甘氨酸分解过程（GO：0006546）、细胞壁组织发生（GO：0071554）、碳水化合物代谢过程（GO：0005975）、单 一 组 织 形 成 过 程（GO：0044699）、水分平衡（GO：0030104）。

2.4 差异表达基因的KEGG 富集分析

在 KEGG 数据库中，统计各个KEGG Pathway 不同层级上包含的差异表达基因数目，进而确定差异表达基因主要参与的代谢途径和信号通路。以整个基因组为背景，采用超几何分布计算枣疯病染病植株花器官差异表达基因显著富集的通路，结果显示1 845 个差异表达基因被注释到237 条代谢途径。其中，172 个差异表达基因富集到细胞进程类别下的20 条代谢途径；221 个差异表达基因富集到环境信息过程下的29 条代谢途径；126 个差异表达基因富集到遗传信息过程下的20条代谢途径；426 个差异表达基因富集到新陈代谢类别下的106 条代谢途径；305 个差异表达基因富集到有机体系统类别下的62 条代谢途径。从差异表达基因数量上看（图3A），枣疯病染病植株主要富集的代谢通路包括淀粉与蔗糖代谢、植物激素信号转导、苯丙烷生物合成、碳代谢、氨基酸生物合成、光合作用、戊糖和葡萄糖醛酸相互转化、糖酵解和糖异生、核糖体、氰基氨基酸代谢等；从显著性上来看（图3B），枣疯病染病植株在12 条代谢通路中具有显著差异，包括光合作用、光合作用-天线蛋白、淀粉与蔗糖代谢、苯丙烷生物合成、光合碳同化、氰基氨基酸代谢、甘氨酸、丝氨酸和苏氨酸代谢、5-羟色胺能突触、戊糖和葡萄糖醛酸相互转化、植物激素信号转导、乙醛酸和二羧酸的代谢、色氨酸代谢。由此可见，淀粉与蔗糖代谢途径、植物激素信号转导途径、苯丙烷生物合成途径、光合作用途径、戊糖和葡萄糖醛酸相互转化途径以及氰基氨基酸代谢途径中差异表达基因数量较多，且存在显著性差异，可能是参与枣疯病植原体响应的重要代谢途径。

图3 差异表达基因的KEGG 代谢途径富集图Fig.3 The KEGG pathways enrichment map of DEGs

2.5 花器官生长发育相关差异表达基因分析

为进一步挖掘枣疯病导致花器官异常发育的关键基因，根据已发表的与成花诱导与发育相关文献和枣基因组注释关键词对差异表达基因进行检索并获得基因编号。从中筛选获得27 个差异表达基因（表3），其中在染病植株中表达上调基因10个，表达下调基因17 个，包括MADS-box 转录因子、AP2 转录因子、SPL 转录因子等。MADS-box转录因子家族有多个基因的表达存在显著差异，其中上调表达的有4 个，即LOC107430150（AGL14）、LOC107421471（AGL12）、LOC10742 5139、LOC107427902；表达下调的 MADS-box 转录 因 子 有 9 个 ，即 LOC107416438（SEP1）、LOC107416339（AGL1）、LOC107421109、LOC 107417456（AGL15）、LOC107421424（AGL11）、LOC107431619（AGL18）、LOC107425386（AGL66）、LOC107425568 （FBP24-like） ，LOC112492135（AGL104）。2 个基因被注释为 AP2 转录因子，即LOC107413613 和 LOC107411289，均在枣疯病染病花器官中表达显著上调。4 个差异表达基因被注释为SPL 转录因子，其中LOC107432336 被注释为SPL9，在染病花器官中表达显著上调；LOC107418038、LOC107408885、LOC107408850均被注释为SPL13，在染病花器官中表达显著下调。图4 利用热图展示了筛选获得的差异表达基因在染病和健康植株花器官的表达情况，并根据差异基因的表达量对样本进行了聚类。

图4 花器官生长发育相关基因差异表达热图Fig.4 Heat map of differentially expressed genes involved in floral organ development

2.6 花器官生长发育相关差异表达基因qRTPCR 验证

为了验证转录组测序数据的准确性并进一步分析枣疯病导致的花器官异常生长的机制，从筛选得到的27 个差异表达基因中选择了8 个基因进行qRT-PCR 检测，分析健康样本和染病样本差异基因的表达变化情况。由图5 可见，8 个基因的FPKM 值与 qRT-PCR 验证结果的变化趋势基本一致，但是表达差异的变化幅度略有不同。与健康花器官相比，染病花器官LOC107421471（AGL12）、LOC107430150（AGL14）、LOC10741 3613（AP2）和 LOC107422311（bZIP27-like）基 因表达显著上调，分别提高170.45 倍、19.02 倍、15.46 倍 、5.41 倍 。 相 反 ，LOC107421424（AGL11）、LOC107425386（AGL66）、LOC10743 1619（AGL18）和 LOC107408885（AGL13）在染病花器官中的表达显著下调，分别下降49.77%、84.10%、75.25%和78.31%。由此可见，枣疯病造成许多花器官生长发育的基因的异常表达。

图5 花器官生长发育相关基因qRT-PCR 相对表达Fig.5 The relative expression by qRT-PCR of differentially expressed genes involved in floral organ development

3 讨论与结论

枣疯病作为影响枣树的严重病害，其病原物与植物体的互作机理一直不是很明确。高通量转录组测序是一种挖掘差异表达基因、探索生物学机理的重要手段。本研究通过对健康与染病的枣树花器官进行转录组测序分析，明确了枣疯病导致枣花异常发育的相关代谢通路，其中涉及光合作用、激素信号、氨基酸合成等众多方面，这与前人的研究报道基本一致［10，11］。本研究还针对枣疯病导致的花器官生长发育差异表达基因进行了挖掘分析，获得的相关基因多属于MADS-box 转录因子、AP2 转录因子、SPL 转录因子家族。MADS-box 基因家族是控制开花和花发育过渡网络的关键组成部分，该家族中有多个AGL基因参与枣疯病导致枣花异常发育的过程。AGL14在拟南芥中已被证明是复杂基因调控网络的重要组成部分，是茎尖分生组织转变的基础，AGL14过表达会导致花器官会有营养器官的特征，同时导致拟南芥早花［31］，枣疯病导致花器官变为营养器官，本研究中LOC107430150（AGL14）在枣树染病花器官中显著上调，这种枣花异常发育的特征与拟南芥的结果一致。对开花起抑制作用的AGL18表达量显著下调，这可能是其与其他AGL基因存在功能互补或存在其他调控途径［32］。AGL12已被证明正向调节根分生组织细胞分裂并与拟南芥的成花诱导相关［33］，本研究中AGL12在染病枣花器官中显著上调，其在枣疯病导致的花器官异常发育中的具体调控作用还有待进一步研究。AGL11作为经典的ABCDE 模型中的D 类基因，在心皮的发育中是必不可少的，本研究中AGL11在染病枣花器官中显著下调，因此很可能在枣花器官异常发育过程中起关键作用。AP2 转录因子参与成花诱导及花器官的分化与调控，是植物花发育调控的重要功能基因，本研究中被注释为AP2 的基因显著上调，这可能打破了枣花器官正常的发育过程。此外，miR156-SPL 调控模块也是控制花发育的重要途径，在植物营养生长和生殖生长阶段转变中发挥重要作用。在番茄中干扰SPL13的表达和miR156a的过表达均会导致果实产量降低［34］，本研究中有3 个下调表达的基因被注释为SPL13同源基因，这也可能是枣疯病染病植株花发育异常的一个关键因素。

e-lik 24 X1 X2BP X2 X1 rm X1 rm fo fo rm rm rm fo tein F fo ro Descriptio fo n A iso 14释13 9 iso A 1-like ox p 4A iso注12 1115CM 11 iso 18 iso 1366-b 1013因GL GL 2tein 2 tein GL GL GL GL tein GL DS GL tein e基tein A LA tein A TA ro LA ro 1-like olog S-like LA tein A 2FL esis M en ro MA ro tein A S -lik ES g-like p ro tein A og TA 7-like ES ro om tein A tein A g-like p e AL TA tein A ro 10 TL ro ro PE ro x p DS g-like p ro x p -lik ro g-like p GA r 2 S-bo in olog 33 PE in tein A S-bo 1 h IN TL EP x p MA IN x p YB T and T orph x p x p in ter-bind x p EIN：21 tein A OT x p in r T AD n facto JO EL ter-bind JO tein S ro r M AD S-bo om tein P S-bo AD R o r o S-bo f m ro S-bo S-bo PR defh S-bo AD scriptio AD tein tic p f F AD ro ter-bind tein B AD AD ter-bind n facto tein ro s-like M s-like M eo ro ro mo mo eo ox p tein ox p la/leafy h ro tic p ro mo tal p n facto HE tic p eo n facto ro OT LITY ro mo ro tran scriptio sa p om ro pm ment s-like M om ou ox p sa p ou Floral h scriptio ou -b en s-like M ou scriptio s-like M ou s-like M sa p uamo ox p s-like M ou ou UA-b W Q s-like M ou sa p DS am tran DS om am Ag IP因elop Ag Sq bZ meob MA Floral h DS Floricau Ho-b uamo am am am am uamo am am uamo IP基ev MA Sq Develo Tran Ag Floral h Ag Protein M Tran Ag Ag Sq Ag LO MA Ag Sq bZ达62070432030812120903 1109041613040202表an d e 120203092002150304异rg P 值E−E−E−E−E−E−E−E−E−E−E−E−E−E−E−E−E−E−E−E−E−E−E−E−E−E−E−差ral o P-valu 28544184930997503612 8171875817 2.1.4.73 9.6.1.20 3.4.9.70 3.3.4.66 4.7.2.64 4.4.84911340592309关5.4.7.3.1.6.4.6.3.2.相育发长olved in flo ge生数倍an 19984412284412926 4948463831151514100505020000000000官器EG s inv异ch ld .8.9.5 406.195.174.2.2.2.2.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.Fo 差花与f D 3 ble 3 L ist o表表 对量PK M相01573940306280712 69866668197076266636280000000000织9..7 44.5 169..9 299.3.3.9.9.7..0 489.9.2..9 191.2.6.7.0.3.0.0.0.0.0.Ta 组达63662841688875 84223967病FP：F 1 7染 表对相M 62494154661024476 501440107902336156226065织量PK .6 2.0..8 5..4 5.4.7.5..3.3 4..9.8 0..5.2.3 4.7.5.8..1.0.6 6.组达159768 14283833 2537 495575 13111964 1419163441康FH：F 1 73 61 02 1健e 50117113883936891402 3838453909569068241985866851355084号am 01231436335123120279 8064286311741929141688535516218854因43424241424243414342 4141414142414243424340424249494040基Gene N 07 C107 C10707 C10707 C107070707 0707070707070707070707070712120707 LO LO LO LO C1 LO C1 LO C1 LO C1 LO C1 LO C1 LO C1 LO C1 LO C1 LO C1 LO C1 LO C1 LO C1 LO C1 LO C1 LO C1 LO C1 LO C1 LO C1 LO C1 LO C1 LO C1 LO C1 LO式attern调调模上下达达达表pression p表表Ex

本试验获得的枣树花器官测序结果较为可靠准确，其中光合作用相关途径、淀粉与蔗糖代谢、苯丙烷生物合成等代谢途径是枣疯病导致花器官异常发育的重要代谢途径。此外，MADS-box 转录因子可能是导致染病枣树花器官异常发育的关键基因，具体调控模式有待进一步挖掘。本研究结果为揭示枣疯病病原物与植物之间互作机制提供了科学参考。