菊花根状茎发育的转录组分析

徐俊 叶雨晴 牛雅静 黄河 张蒙蒙

（1. 北京林业大学园林学院 花卉种质资源创新与分子育种北京市重点实验室 国家花卉工程技术研究中心 城乡生态环境北京实验室，北京 100083；2. 北京市植物园管理处 北京市花卉园艺工程技术研究中心，北京 100093）

菊花（Chrysanthemum × morifolium）是中国传统十大名花之一，有着3 000多年的栽培历史，具有极高的观赏和经济价值。在规模化生产中，菊花主要通过由根状茎发育而来的脚芽进行扦插繁殖。根状茎（rhizome）是指植物在地下水平生长的变态茎，是许多多年生植物的营养器官之一［1-2］。作为一种繁殖策略和营养器官，根状茎能够帮助植物抵御非生物胁迫，例如，根状茎型多年生禾本科植物草地早熟禾（Poa pratensis）和结缕草（Zoysia japonica），比非根状茎型禾本科植物多年生黑麦草（Lolium perenne）表现出更好的耐旱性［3］。同时，根状茎能增加植物的宿根性和繁殖能力，在普通水稻（Oryza sativa）中引入长雄野生稻（O.longistaminata）的根状茎性状能培育出一次种植多年采收的多年生水稻［4］。此外，荷花（Nelumbo nucifera）、生姜（Zingiber officinale）等农作物的根状茎还是其主要的采收器官。有研究表明具有更多根状茎的菊花品种抗寒性更强，根状茎数量是抗寒菊花品种快速筛选的重要指标［5］。并且菊花根状茎还有助于抵御干旱胁迫，帮助菊花在干旱后更好的恢复［6］，但相对于在花色、开花期等方面的表型分析和分子调控机制研究［7-9］，目前对于菊花根状茎形成关键基因的克隆及其功能还少有报道。

目前控制高等植物根状茎形成的关键基因尚报道较少，随着转录组测序技术的日趋成熟，该技术被广泛应用到根状茎发育研究当中，在长雄野生稻、拟高粱（Sorghum propinquum）、芦苇（Phragmites australis）、早竹（Phyllostachys praecox）等植物中，通过转录组测序发现了一些可能与根状茎发育相关的基因，包括激素合成和信号传导、光周期响应以及调控分生组织形成相关的基因。植物激素是植物体内产生的在极低浓度下能够产生生理效应的信号分子，参与植物生长发育各个阶段的调控。多种植物的根状茎转录组报道中筛选出了与植物激素合成和信号传导相关的基因，如在荷花根状茎生长发育过程中，赤霉素信号转导基因Gibberellic Acid Insensitive（GAI）、脱落酸受体基因Pyrabactin Resistant‑Like（PYL）和生长素响应因子Auxin Response Factor（ARF）在伸长期和膨大期的表达差异显著［10］，赤霉素合成基因Gibberellin 20‑oxidase（GA20ox）的表达量在根状茎发育过程中降低［11］。此外，对强根状茎草地早熟禾QH和弱根状茎草地早熟禾SN的根状茎芽、根状茎节和根状茎节间3个部位分别进行了转录组测序，发现脱落酸受体基因PYL和脱落酸信号传导基因PP2C在QH和SN的节中差异表达［1］。以上研究表明，根状茎的生长发育可能受到赤霉素、脱落酸等多种激素的调控。由于根状茎起源于茎基部的腋生分生组织［2］，因此分生组织形成相关基因也可能参与控制根状茎的形成。在拟高粱转录组中找到了在根状茎中特异表达的MYB36，该基因在拟南芥（Arabidopsis thaliana）中控制腋生分生组织的形成［12］；同时，在早竹根状茎转录组中也筛选到了REVOLUTA（REV）、CLAVATA1（CLV1）等和分生组织形成有关的基因［13］。另外，有报道显示，光周期诱导途径关键基因也可能参与根状茎的形成和发育，FT基因的表达在温代莲和热带莲的根状茎发育过程中下降［10］。

综上，目前的研究对高等植物根状茎的发育过程有了一定了解，但根状茎形成和生长发育的分子基础尚不清楚，诱导根状茎形成的关键基因仍未找到。为了探究菊花根状茎形成的分子机制，筛选可能参与菊花根状茎发育的基因。本研究以能稳定形成根状茎的菊花品种‘2017XS’为转录组测序材料，选取根状茎尖、根状茎中部、根状茎下部、叶片、茎段、根系、茎尖、舌状花8个组织进行转录组测序。通过对测序结果的生物信息学分析，筛选有可能控制根状茎发育的关键基因，并利用RT‑qPCR对其进行表达验证，从而为解析菊花根状茎形成机理，培育具有丰富根状茎、抗逆性强的菊花新品种奠定理论基础。

1 材料与方法

1.1 材料

前期通过对大量的菊花种质资源进行筛查，发现‘2017XS’和‘2005042’两个株系能够稳定形成根状茎，并且数量较多（图1）。本研究以‘2017XS’为转录组测序材料，为了获取全面的转录组数据并且获得根状茎中差异表达基因，除了根状茎的3个部位根状茎尖（RH）、根状茎中部（RT）、根状茎下部（RB）还选取了茎尖（SA）、叶片（L）、茎段（S）、根（R）、舌状花（F）5个部位进行测序（图1‑C）。所有的样品取自无性繁殖的植株，每个部位取3次重复。样品用液氮速冻，并储存在‑80℃条件下。

图1 ‘2017XS’和‘2005042’根状茎示意图及转录组取材示意图Fig. 1 Schematic diagram of ‘2017XS’ and ‘2005042’rhizomes and schematic diagram of 8 transcriptome sampling

以‘2005042’为验证品种，同样取其根状茎尖、根状茎中部、根状茎下部、茎尖、叶片、茎段、根和舌状花8个部位，验证转录组测序中获得差异表达基因的表达模式。

1.2 方法

1.2.1 RNA 提取与文库构建 使用植物RNA快速提取试剂盒（北京华越洋生物科技有限公司）分别提取总RNA。用带有Oligo（dT）的磁珠富集mRNA，加入破碎缓冲液将mRNA 随机打断，并以mRNA为模板，用六碱基随机引物合成第1条cDNA链，然后加入缓冲液、dNTPs、RNase H和DNA polymerase I合成第2条cDNA链，利用AMPure XP beads纯化cDNA。纯化的双链cDNA再进行末端修复，加A尾并连接测序接头，然后进行片段大小的选择，最后进行PCR扩增。使用Agilent 2100生物分析仪和ABI StepOnePlus实时PCR系统对文库进行质检。

1.2.2 转录组测序与组装 利用Illumina Hi‑seqTM2000测序平台进行测序。过滤raw data中的低质量数据以获得clean data，随后使用Trinity软件［14］进行从头组装，得到Unigene序列。获得的Unigene序列通过BLAST［15］（E‑value ＜ 10-5）与NR、Swiss‑Prot、GO、COG、KOG、eggNOG、KEGG数据库进行比对，使用KOBAS2.0［16］得到Unigene在KEGG中的KEGG Orthology结果，预测完Unigene的氨基酸序列之后使用HMMER［17］（E‑value ＜ 10-10）软件与Pfam数据库比对，获得Unigene的注释信息。利用clusterProfiler［18］软件对差异表达基因进行GO富集和KEGG通路富集分析，并选取Padj ＜ 0.05的显著富集的GO类别和KEGG代谢通路。

采用Bowtie［19］将测序得到的Reads与Unigene库进行比对，根据比对结果，结合RSEM［20］进行表达量水平估计。利用FPKM（Fragments Per Kilobase of transcript per Million mapped reads）值表示对应Unigene的表达丰度。采用DESeq2进行样品组间的差异表达分析，以FDR（False Discovery Rate）＜ 0.01且差异倍数FC（Fold Change）≥2为筛选标准。差异表达基因的聚类分析使用Mfuzz软件进行［21］。使用Tbtools软件制作基因表达量热图［22］。Venn图使用在线平台绘制（https://www.bioladder.cn/web/#/chart/17）。

1.2.3 加权基因共表达网络分析 首先通过FPKM将低表达的和低变异度的基因去除。数据过滤之后对样本进行层次聚类分析。利用pickSoftThreshold函数计算最佳软阈值，选取无尺度网络拟合指数R2＞ 0.8时power值最小的数为最佳power值。利用WGCNA包的blockwiseModules函数构建共表达矩阵，模块相似度阈值设置为0.25，合并相似度为0.8的模块，模块内最小基因数设置为30，其他参数按照默认设置［23］。

1.2.4 根状茎高表达基因的筛选 使用Excel按照根状茎尖的FPKM表达量高于其他组织3倍的标准筛选根状茎尖高表达基因，按照根状茎3个部位（根状茎尖、根状茎中部、根状茎下部）的FPKM表达量均高于其他部位3倍的标准筛选根状茎高表达基因［24］。

1.2.5 RT‑qPCR验证 根据转录组筛选结果，选取6个在根状茎中特异表达的基因，在‘2017XS’和另一个稳定形成根状茎菊花株系‘2005042’的8个部位（根状茎尖、根状茎中部、根状茎下部、茎尖、叶片、茎段、根和舌状花）进行基因表达的荧光定量PCR（RT‑qPCR）分析，使用 primer premier 5.0设计引物（表1）。用德国耶拿实时荧光定量 PCR 仪qTOWER 2.2（Analytik Jena，German）进行反应，每个样品重复3次，反应体系20 μL。选用SAND作为内参基因，用 2-ΔΔCt算法计算相对表达量［25］。

表1 RT-qPCR引物序列Table 1 Primer sequences for RT-qPCR

2 结果

2.1 测序和从头组装

为了筛选得到根状茎尖以及整个根状茎中高表达的基因，对来自叶片、茎、根、舌状花、茎尖、根状茎尖、根状茎中部和根状茎下部的24个 cDNA文库使用 Illumina HiSeq 2000 高通量测序平台进行测序。在过滤掉接头序列、不明确和低质量的reads之后，所有样品总的clean reads大约159.51 GB。使用Trinity程序，将所有clean reads从头组装成196 467个转录本，平均长度为905 bp，N50长度为1 295 bp。这些转录本进一步组装成100 235个Uni‑gene。Unigene平均长度为832 bp，N50长度为1 280 bp。在这些Unigene中，28 114个（28.05%）大于1 000 bp，46 596个（46.48%）小于500 bp（图2‑A）。

2.2 根状茎差异基因的功能注释和分类

为了对组装的Unigene进行注释，对7个公共蛋白质/核苷酸数据库进行了E值为10-5的BLAST搜索，最终获得64 956（64.80%）个有注释信息的Unigene。其中3 1951条（49.19%）在KOG中得到注释，40 653条（62.47%）在Pfam中得到注释，34 388条（52.86%）在Swissprot中得到注释，42 639条（65.64%）在eggNOG中得到注释，62 005条（95.46%）在Nr数据库中得到注释。

为深入了解根状茎发育基因的功能分类和代谢途径，将根状茎和其他组织间进行比较获得差异基因（FDR＜ 0.01，FC≥ 2），共获得3 274个差异基因。将这些基因根据序列同源性，分为25个KOG类别（图2‑B）。一般功能预测基因（273，18.03%）、翻译后修饰，蛋白质转运和蛋白伴侣（159，10.5%）、信号传导机制（136，8.98%）和碳水化合物转运和代谢（136，8.98%）类别包含基因较多。相比之下，染色体结构和动态（9，0.59%）和细胞外结构（7，0.46%）的包含基因较少。在3 274个差异基因中，共有1 190个Unigene在KEGG数据库中匹配，并匹配到50条通路，覆盖五大KEGG类别（图2‑C）。最具代表性的3种途径是植物-病原体互作（91, 8.95%）、植物激素信号传导（55, 5.41%）和碳代谢（53, 5.21%），其次是植物MAPK信号通路（50, 4.92%）和苯丙酸生物合成（48, 4.72%）途径，而氨基酰基tRNA生物合成（11, 1.08%）和油菜素生物合成（11, 1.08%）途径包含基因最少。

2.3 加权基因共表达网络分析

为了鉴定根状茎起始和发育过程中高表达的基因，在去除低表达和低变异度的基因后构建相关性聚类树，采用动态切割法将产生的聚类树进行切割，把表达模式相似的基因合并在同一分支上，每个分支代表1个共表达模块，根据模块相似度对表达模式相似的模块合并（合并相似度为0.8的模块）后进行模块划分，最终获得17个共表达模块（图3‑A）。Meroyalblue模块中共有460个基因，该模块与根状茎正相关；Meplum1模块中共有632个基因，该模块与根状茎尖正相关（图3‑B），基于这两个模块我们进行了根状茎尖和根状茎高表达基因的筛选。利用FPKM值衡量基因的表达水平，按照方法1.2.4中的标准筛选高表达基因。基于此，在MEroyalblue模块中共筛选到20个在根状茎中高表达的基因，在MEplum1模块中共筛选到16个在根状茎尖中高表达的基因，共计36个基因，根状茎尖中高表达的16个基因主要包括和抗逆相关的基因（LEA1、TIP1‑1、TIP2‑1）、脱落酸代谢相关的基因（ABA8ox）、红光受体基因（PHYB）以及1个bZIP类转录因子（bZIP34）（表2）。根状茎中高表达的基因同样包含大量抗逆相关的基因（LEA1、PIP1‑1、TIP2‑1和DREB1D），同时还包含了光周期转录因子（TFL1）、紫外光受体基因（UVR8）、生物节律控制基因（LHY）基因以及锌指蛋白（zinc finger）（表3）。

表3 ‘2017XS’根状茎中高表达基因汇总Table 3 Summary of highly expressed genes in the rhizomes of‘ 2017XS’

图3 加权基因共表达网络分析Fig. 3 WGCNA analysis

2.4 基于韦恩图筛选差异基因

为了鉴定根状茎起始和发育过程中高表达的基因，进一步着重于根状茎3个部位（根状茎尖、根状茎中部、根状茎下部）中DEGs的鉴定。在根状茎尖和其他非根状茎部位间的比较中有693个基因都差异表达（图4‑A），618个基因在根状茎中部和其他非根状茎部位都差异表达（图4‑B），有763个基因在根状茎下部和其他非根状茎部位都差异表达（图4‑C）。在这2 074个基因中筛选根状茎尖中高表达基因和根状茎高表达基因。按照方法1.2.4的标准筛选到16个在根状茎尖中高表达的基因，其中除了11个WGCNA筛选到的基因，还包含了抗逆相关转录因子（DREB2F）、锌指蛋白（CZF2）、ABA降解基因（ABA8ox）、ABA信号转导基因（PP2C）和糖基转移酶基因（UGT83A1）等（表2）。根状茎高表达的基因在WGCNA筛选到20个差异基因的基础上新筛选到5个基因分别是锌指蛋白（CZF4）、NADH脱氢酶基因（NDA1）、乙烯响应因子基因（ERF C3）、β-糖苷酶基因（BG 46）和胚胎晚期发育蛋白基因（LEA Dc3）（表3）。

图4 根状茎3个部位与其他组织的差异表达基因Venn图Fig. 4 Venn diagram of DEGs in three parts of rhizome and other tissues

2.5 基因共表达趋势分析

为了筛选根状茎发育的关键基因，采用K‑means聚类对Venn图中涉及的差异基因进行分类（图5）。共鉴定出25个具有不同表达模式的聚类，其中有8个聚类在单个组织中高表达。聚类3和16在舌状花中高表达，聚类5、13和18在根中高表达，聚类22在叶片中高表达，聚类14在茎尖中高表达，聚类1在根中高表达，没有发现在根状茎尖中高表达的聚类模块。但发现聚类17中的473个基因在根状茎3个样品中高于其他组织。利用FPKM值衡量基因的表达水平，按照根状茎3个部位中的表达量均高于其他组织3倍的标准筛选根状茎高表达基因。在WGCNA和韦恩图筛选的基础上新筛选到10个在根状茎中高表达的基因，主要包括光周期响应基因（COL12）、紫外光受体基因（UVR8）、非生物胁迫基因（LEA）和糖基转移酶基因（UGT）（表3）。

图5 根状茎稳定型株系‘2017XS’基因表达模式的Cluster分析Fig. 5 Cluster analysis of gene expression patterns of rhizome stable line ‘2017XS’

综上，基于以上3种方式共筛选到了20个在根状茎尖中高表达的基因，36个在根状茎中高表达的基因。这些基因涉及了激素代谢和信号传导、光周期、非生物胁迫等重要的生物过程。

2.6 RT‑qPCR验证差异基因的表达情况

为了验证 RNA‑seq 结果准确性，选择6个筛选出来的高表达基因进行RT‑qPCR 验证，分析 RT‑qPCR 与 RNA‑seq 结果之间的相关性，每个基因进行 3 次生物学重复。转录组中，bZIP34和ABA8ox在根状茎尖中高表达，其余4个在根状茎中高表达。结果显示，6个基因的表达量变化与转录组测序结果的变化趋势一致，说明转录组的测序数据和RT‑qPCR 的数据之间具有良好的一致性（图6‑A），表明本研究 RNA‑seq 数据具有较高的可靠性。在‘2005042’中，6个基因的表达在根状茎尖中的表达均显著高于其他部位，结果和转录组数据以及‘2017XS’的表达验证趋势较为一致（图6‑B）。

3 讨论

作为一种在园林中大量应用的地被花卉，菊花具有开花繁密、株型整齐等特点。在园林应用中，根状茎特性赋予了菊花极强的宿根性和无性繁殖能力。但目前关于菊花根状茎的研究较少，菊花根状茎形成的分子机制仍不清楚。本研究对包括菊花根状茎在内的8个组织进行了转录组测序分析，以期筛选可能参与根状茎发育的基因。

3.1 菊花根状茎及其在抗逆中的作用

根状茎是植物在地下水平生长的地下茎，源于腋生分生组织［2］。根状茎在顶端分生组织在发育的过程中既能形成一个新的节间维持原来的状态，也能使根状茎尖向上弯曲形成母株的一个营养克隆［2］。同时根状茎的节间能形成不定根从而拓展植物的根系。根状茎作为一种重要的地下贮藏器官有助于植物在恶劣环境下生存，在胁迫条件下根状茎能够在地下生存，一旦条件合适根状茎就能发育形成地上茎。有大量文献报道根状茎和植物非生物胁迫有关。在草地早熟禾中，根状茎中贮藏的非结构性碳水化合物有助于其干旱后恢复［26］。在菊花中根状茎同样在非生物胁迫中起着重要作用，Anderson等［5］发现根状茎越多的菊花冬季存活率越高，并且可以通过根状茎表型进行抗寒菊花的快速筛选。Zhang等［6］认为菊花根状茎有助于菊花在干旱后的恢复［27］，并且发现干旱胁迫能降低CmRH56的表达从而较少根状茎的数量。

本研究中在根状茎尖和根状茎中找到了大量和抗逆相关的基因包括TIP、PIP、DREB和LEA。目前大部分关于DREB报道都和植物抗逆有关，在小麦（Triticum aestivum）中过表达大豆（Glycine max）的GmDREB1增加了小麦的抗旱能力，并且在干旱条件下的转基因小麦的产量高于野生型［28］；过表达番茄（Solanum lycopersicum）SlDREB3改善了番茄在冷胁迫下的生理状态，使植株在4℃环境下有更好的生长表现［29］。在辣椒（Capsicum annuum）中抑制LEA同源基因CaDIL1表达，降低了辣椒在干旱胁迫下的表现［30］。AQPs在植物的抗逆胁迫方面有大量的报道，同时也有文献报道其参与种子的萌发。目前未见这些基因在根状茎发育中起作用的报道，但有研究表明干旱能够诱导菊花的根状茎产生［27］，我们推测抗逆类基因有可能在干旱胁迫和根状茎发育间起到桥梁作用。

3.2 脱落酸影响菊花根状茎发育

越来越多的证据表明，激素影响植物发育过程的各个方面，脱落酸和赤霉素对植物的地下茎的休眠萌发和生长发育起到了重要作用，袁娅娟［31］对草地早熟禾不同生育时期的根状茎扩展指标和激素含量进行了分析，认为ABA含量越低，越有利于根状茎扩展；Li等［32］发现屏边空竹（Cephalostachyum pingbianense）的根状茎芽在休眠到萌发期间ABA含量下降，并且与ABA合成和信号转导途径相关的调控基因，如NCED、PYL、SnRK2和ABF，在根状茎芽萌发期间也显著下调。此外，赤霉素也会影响根状茎发育，外施赤霉素可以促进高羊茅和草地早熟禾根状茎伸长［33］；在菊花中，外施赤霉素能增加其根状茎数量［6］。

在本研究中找到了2个ABA代谢基因（ABA8ox）、5个调控ABA稳态基因（UGT和BG）以及1个ABA信号传导基因（PP2C），在根状茎尖或整个根状茎中高表达。但未找到和赤霉素相关的基因。ABA浓度和活性的调控是一个复杂的过程，ABA可以通过复杂的生化过程由2‑C‑甲基‑D‑赤藻糖醇‑4‑磷酸（Methylerythritol‑4‑phosphate pathway，MEP）起始合成ABA，整个过程中ABA合成的关键限速酶是NCED［34］。ABA的分解代谢主要通过8′羟基化实现，主要的酶是ABA8ox［35］。除外UDP‑糖基转移酶基因UGTs和β-糖苷酶基因BGs在调控ABA活性方面发挥作用，UGTs可以将葡萄糖Glu结合到ABA上形成ABA葡萄糖酯（ABA‑GE），使ABA失活，BGs 基因可以将Glu去除从而恢复 ABA活性［31］。在葡萄中过表达ABA8ox同源基因VvA8H‑CYP707A4能降低ABA含量，并促进葡萄侧枝生长［35］。在番茄中过表达柿树（Diospyros kaki）的UGT3能降低番茄ABA含量，使番茄出现植株变矮、种子休眠程度降低等表型［36］。柿树DkBG1可以水解ABA‑GE，从而释放游离ABA，在番茄中过表达DkBG1可以促进果实提前成熟［37］。本研究找到了在根状茎尖或者根状茎中高表达的2个ABA8ox、3个UGT和2个BG。因此ABA很有可能参与到菊花根状茎形成和发育过程中。

3.3 光周期及光照影响菊花根状茎发育

光周期对高等植物地下根状茎的发育具有调控作用，在荷花中，短日照促进根状茎增粗，而长日照促进根状茎伸长［38］。在高羊茅中较长的光周期有利于根状茎的形成［39］。在竹类植物中光周期核心转录因子FT与其根状茎类型有密切关系［40］。在本研究中我们筛选到了两个光周期核心转录因子COL12和TFL1，一个红光受体基因PHYB和两个紫外光受体基因UVR8在根状茎中特异高表达。

目前有大量文献报道了光周期对根状茎发育具有重要影响，其中光周期对于马铃薯（Solanum tuberosum）块茎和草莓（Fragaria vesca）匍匐茎形成的分子机制研究得较为清楚，在马铃薯中共有3个CO‑like基因（StCOL1/StCOL2/StCOL3），其中StCOL1在长日照下优先积累。StCOL1在RNAi系中的下调加速马铃薯块茎的形成。StCOL1可以直接激活长日照下的FT同源基因StSP5G。激活后，StSP5G通过负调节FT同源基因StSP6A转录来抑制非诱导条件下的块茎形成［41］。

叶片感受光信号是通过光受体完成的，如光敏色素和隐花色素等［42］。在本研究中还筛选得到了两个紫外光受体基因UVR8在根状茎中高表达，一个红光受体基因PHYB在根状茎尖中高表达。这些光受体基因可能通过感受光周期进而影响植物地下茎发育。在马铃薯基因组中鉴定的5个光敏色素基因中，StPHYB和StPHYF在响应长日照抑制块茎形成中发挥最显著的作用。抑制StPHYF表达发现马铃薯能在长日照下形成块茎并且StPHYF和StPHYB能够互作［43］。

4 结论

本研究利用 Illumina Hi‑seqTM2000 高通量测序平台对菊花种质‘2017XS’的8个组织进行了转录组测序，共获得100 235个Unigene，其中64 956个Unigene 得到注释。通过WGCNA、Venn图筛选和K‑means一共筛选得到20个在根状茎尖中高表达的基因和36个在根状茎中高表达的基因，主要包括ABA分解代谢、光周期相关的基因以及部分抗非生物胁迫类基因。