基于转录组测序的枇杷腋芽成花相关基因的筛选

章加应,肖海涛,2*,张学英**,安海山,徐芳杰,周博强

(1上海市农业科学院林木果树研究所,上海 201403;2上海应用技术大学生态技术与工程学院,上海 201418)

成花是植物生长发育过程中一个重要的生长阶段,是对周围环境的适应性表现,标志着植物从营养生长向生殖生长转变[1]。植物成花的本质为茎顶端分生组织由营养生长阶段分化产生的叶片转变为生殖发育阶段形成花、果实和种子的过程[2],这一过程主要分为3个阶段,即成花诱导、花芽分化和花器官形成阶段[3]。植物成花过程是一个高度复杂的生理过程,在外部和内部各种因素形成的成花调控网络下完成。通过外部环境和内部信号对拟南芥开花过程调控的分子遗传机理研究鉴定了多种调控路径[4-5],发现光周期途径、春化途径、温度途径和赤霉素途径等传递光、温度和激素等信号调控成花;而自主途径和年龄途径在很大程度上以内源信号为主。这些途径通过主要的整合基因,实现对植物成花和开花时间的精准调控[6-7]。FLOWERING LOCUS T(FT)和SUPPERSSOR OF OVEREXPRESSION OF CONSTANS1(SOC1)作为最重要的开花整合基因,在苹果[8-9]、梨[10]、桃[11]和枇杷[12-13]等果树开花过程中发挥着至关重要的作用。枇杷(Eriobotrya japonicaLindl.)属蔷薇科(Rosaceae)枇杷属植物,是原产于我国的亚热带常绿树种,枇杷果实富含大量营养物质,备受消费者青睐,目前已成为我国南方地区重要的经济果树之一[14]。枇杷以顶芽成花为主,秋冬开花,初夏果实成熟,关于枇杷成花的相关机理研究集中在枇杷顶芽成花方面。在枇杷顶芽花芽诱导阶段相对较高水平的玉米素(ZT)和脱落酸(ABA)利于花芽分化;形态分化期,高水平的ABA利于花芽的形成[15]。张玲等[16]对台湾枇杷‘恒春’变型品种顶芽花期调控进行研究,成功克隆到FT、CO、GI、SOC和PIF4等开花相关基因,并在拟南芥中验证了EdFT、EdFD1、EdFD2、EdCO、EdGI和EdSOC1等基因促进花期提前的生物学功能。EjFT、EjAP1和EjSOC1基因是枇杷开花过程中的重要成花整合因子,在顶芽花芽萌动和花芽诱导过程中发挥着重要作用[13-14,17]。但有关枇杷腋芽的成花生理和分子遗传机制的相关研究鲜有报道。本课题组在2010年定植于金山果树试验站的枇杷实生后代中发现一棵具有腋芽成花特点的枇杷单株(暂定名:春花),顶花芽秋冬开花,腋芽3—5月份开花结果,但有关春花枇杷腋芽成花过程中的转录组信息和基因表达谱尚未不清。随着基因组测序和分子生物技术的不断进步,高通量测序近年来发展迅速,RNA测序(RNA-seq)可实现转录本和差异表达基因的同步分析,其结果可以在动态范围准确鉴定基因表达情况,能够快速发现差异表达基因,并对转录组进行更加敏感和准确的分析,从而更好地解析分子机制[18]。近年来,该技术也被应用到枇杷生物学性状中差异候选基因的挖掘,成功揭示了枇杷生长发育[19]、果实成熟[20]和抗逆性[21]的分子机理,通过RNA-seq筛选与枇杷腋芽成花相关的关键基因,有助于解析枇杷腋芽成花的分子调控机制。

本研究对枇杷品系春花(具有腋芽特性)和品种‘火炬’(不具有腋芽特性)腋芽进行转录组测序和分析,从转录组数据中筛选与成花相关的差异表达基因,以期为枇杷腋芽成花调控提供重要的候选基因,并为阐明枇杷腋芽成花分子遗传机制奠定理论基础。

1 材料与方法

1.1 试验材料

供试枇杷品系春花(上海市农业科学院通过实生选育而成的大果红肉枇杷优系,具有腋芽成花特性)和‘火炬’(上海市农业科学院通过杂交选育的大果红肉枇杷品种,不具有腋芽成花特性)(图1),均为5年生果树,生长状态基本一致,定植于上海市农业科学院金山果树试验站。于2020年于2月中旬—4月底每10 d对腋芽进行取样(表1),每个时期将10个独立的腋芽收集一起作为一个重复,3次生物学重复。采集后将样品放于冰盒中带回实验室,液氮速冻后,-80℃保存,备用。

图1 春花和‘火炬’腋芽发育进程Fig.1 Development process of axillary buds of Chunhua and‘Huoju’

表1 春花和‘火炬’腋芽取样时间Table 1 Sampling time of axillary buds of Chunhua and‘Huoju’

1.2 试验方法

1.2.1 RNA提取和RNA-seq

总RNA提取按照mirVana miRNA提取分离试剂盒(Ambion)说明书上的步骤进行操作提取。使用Agilent 2100生物分析仪(Agilent Technologies,Santa Clara,CA,USA)对RNA完整性进行评估,RNA Integrity Number(RIN)≥7.0的样本进行后续分析。通过使用TurSeq Strand mRNA LTSample PrepKIT(Illumina,San Diego,CA,USA)构建文库。由北京百迈客生物科技有限公司采用Illumina-Hiseq测序平台进行转录组测序,生成125 bp∕150 bp的配对端读数。

1.2.2 转录组数据的组装

为保证后续组装质量,提高生物信息分析结果的正确性,应用SeqPrep和Sickle软件去除原始测序数据中的测序接头、引物序列和低质量值的读段,处理后得到高质量测序数据。后期对原始数据过滤得到的RNA-seq测序数据进行从头组装。使用Trinity和Inchworm软件搜索方法建立K-mer graph(K=25)中全部可能路径[22];采用Chrysalis对上述得到可能路径进行分组后,对Chrysalis生成路径进行修剪和整合,保存主要路径。同时利用Trinity(版本:trinityrnaseq_r20131110)将clean reads组装成表达序列标签簇(contigs),重新组装成转录本[23]。根据序列的相似性和长度选择最长的转录本作为unigene进行后续分析。

1.2.3 基因功能注释

通过对unigenes与NCBI non-redundancy(NR)、SwissProt和真核完整基因组(KOG)数据库的同源群进行比对(使用阈值e值为10-5的Blastx),对unigenes的功能进行注释。对unigenes点击最高的蛋白进行功能注释。在SwissProt注释的基础上,利用SwissProt与GO术语间的映射关系进行基因本体(GO)分类。unigenes被映射到京都基因和基因组百科全书(KEGG)数据库,以注释其潜在的代谢途径。

1.2.4 腋芽成花相关差异表达基因的鉴定

将测序得到的Reads与Unigene库进行比对,根据比对结果结合RSEM进行表达量水平估计。利用RPKM(reads per kilobase per millionreads)表示对应Unigene的表达丰度。在差异表达分析中采用校正后的P值,即FDR(false discovery rate)＜0.05,且差异倍数FC(fold change)≥2作为差异表达基因(DEGs)的筛选标准,通过聚类分析探讨转录产物的表达模式。基于超几何分布,分别用R法对DEGs进行GO富集和KEGG通路富集分析。

1.2.5 qRT-PCR基因表达量分析

提取枇杷腋芽每个时期的总RNA。为了进行定量逆转录PCR(qRT-PCR)分析,使用SYBR Prime Script RNA RT-PCR试剂盒(日本TaKaRa)进行第一链cDNA合成。qRT-PCR使用Light Cycler 480(Roche,USA)进行分析。采用SYBR-Green PCR试剂盒(日本TaKaRa)进行qRT-PCR,反应液为20μL,其中2×SYBR-PreMix EX Taq为10μL,cDNA为50 ng,每个引物为0.25μmol∕L。qRT-PCR程序设置如下:94°C变性5 min;94变性10 s,60°C变性30 s,72℃变性30 s;72℃变性3 min。每个基因有3个生物重复序列和3个技术重复序列。采用2-ΔΔCt算法计算各基因的相对表达量。以β-ACTIN为内参基因。利用SnapGene软件设计引物(表2)。

表2 qRT-PCR引物设计Table 2 Primer design of qRT-PCR

1.3 数据分析

运用Excel 2013和SPSS 17.0软件对试验数据进行整理和分析。采用Graphpad prism 8.0软件(Graphpad Software,USA)进行数据统计和图表制作。利用SPSS 17.0软件中单因素方差分析(ANOVA)检测各数据间的显著水平(P＜0.05)。

2 结果与分析

2.1 枇杷腋芽成花时期确定

枇杷是以顶芽成花为主,本课题组前期对枇杷顶芽的花芽分化期进行RNA-seq测序,筛选到调控顶芽成花的大量候选基因,包括EjFT、EjFY、EjFLK和EjCAL1-like等(数据待发表),本研究于2018年12月28日—2019年5月13日对春花枇杷的腋芽进行连续取样并分析了这些基因在枇杷品种春花腋芽中的时空表达模式(图2)。结果显示:EjFT基因在3月13日—4月19日出现表达峰值,EjFY基因在4月19日出现表达峰值,EjFLK和EjCAL1-like基因在2月11日和4月11日出现2次表达高峰。由此可见,2月中旬至4月底可能为枇杷腋芽成花关键期。

图2 开花基因在春花枇杷腋芽中的时空表达水平Fig.2 The spatial-temporal expression level of flowering genes in axillary bud of Chunhua loquat

2.2 枇杷腋芽转录组测序

2.2.1 转录组测序数据分析

为进一步筛选调控春花枇杷腋芽成花的候选基因,于2月中旬—4月底每10 d对春花和‘火炬’的腋芽进行取样并进行RNA-seq测序。结果显示(表3):所测样品的读长长度为19 832 347—24 142 070 bp,碱基数范围为5 911 148 778—7 209 277 186 bp,其中所有样品的Q30值都在93%以上,同时GC含量为47.26%—48.09%,表明测序数据质量较高,可以进行下一步的基因组组装和差异表达基因分析。

表3 测序所得Clean Data统计Table 3 Clean data statistics from sequencing

2.2.2 转录组数据的组装

利用Bowtie2(v2.1.0)软件将转录组测序数据Clean Data进行基因组有参比对(参考基因组见http:∕∕creativecommons.org∕licenses∕by∕4.0),各样品的总比对数介于33 680 111—44 131 623 bp(比对率:76.88%—91.50%),比对到参考基因组多个位置的比对数在1 100 109—1 419 663 bp,占所有比对数的2.41%—3.00%。将比对到基因组唯一位置的比对数用于后续的基因表达分析,其数量为32 557 477—42 853 791 bp,占所有比对数的74.39%—88.85%(表4)。

表4 基因比对结果Table 4 Blast results of genes

2.2.3 差异表达基因的筛选

对春花和‘火炬’枇杷腋芽在不同时间点转录组的表达量进行计算,共鉴定到66 197个差异表达基因(表5),其中2月27日—4月8日样品的上调表达基因比较丰富,暗示这段时间可能是枇杷腋花芽分化的关键期(表1、图3a)。维恩图显示,不同处理同时上调表达的差异基因有160个(图3b),同时下调表达的差异基因有198个(图3c)。

图3 枇杷腋芽转录组测序筛选到的差异表达基因Fig.3 Differential expression genes from transcriptome sequencing screening in loquat axillary buds

表5 枇杷腋芽转录组筛选的差异表达基因统计Table 5 Statistics of differential expression genes from transcriptome screening in loquat axillary buds

2.2.4 枇杷腋芽差异表达基因GO富集分析

对春花和‘火炬’不同时间腋芽中的差异表达基因进一步进行GO富集分析,确定差异表达基因主要行使的生物学功能。枇杷腋芽中的差异表达基因在GO注释中共分为3个功能组(图4),包括20个生物过程、16个细胞组分和15个分子功能。生物过程中,注释分类到与生物代谢、细胞转化、组织进化和生物调节等相关过程的差异表达基因较多;在细胞组分中,差异表达基因主要分布在细胞器、细胞膜和质体等部位;在分子功能方面,差异表达基因主要集中于催化活性、物质转运和信号转导等功能。

图4 枇杷腋芽中差异表达基因GO功能注释Fig.4 GO functional annotation of differential expression genes in loquat axillary buds

2.2.5 枇杷腋芽差异表达基因KEGG通路显著性富集分析

腋芽差异表达基因被注释到KEGG数据库中主要的20个代谢通路上(图5)。植物激素信号传导(Plant hormone signal transduction)、淀粉和糖代谢(Starch and sucrose metabolism)、苯基丙烷代谢(Phenylpropanoid biosynthesis)、病原体互作反应(Plant-pathogen interaction)和谷胱甘肽代谢(Glutathione metabolism)为最重要的5种显著特异性富集KEGG通路,表明这5种通路在枇杷腋芽开花过程中发挥着重要作用。激素作为植物生长发育过程中的重要调节物质,其代谢影响着植物的发育进程,激素信号传导通路可作为枇杷腋芽成花过程中的一个重点研究方向。同时,KEGG代谢通路中表达量上调的基因数量显著高于表达量下调的基因数量。

图5 枇杷腋芽中差异表达基因KEGG通路显著性富集分析Fig.5 Analysis of KEGG pathways significance enrichment of differential expression genes in loquat axillary buds

2.2.6 腋芽中差异基因的时空表达模式

经GO富集、KEGG富集和NR功能注释后,进一步从枇杷品系春花和品种‘火炬’腋芽转录组筛选到32个成花基因、1个ABA受体基因PYL4、1个ABA脱氢酶基因ABH4和大量功能未知的新基因(表6)。基于RPKM值,分析差异表达基因在春花和‘火炬’腋芽中的时空表达模式。基因表达模式分析显示(图6和图7):FT2基因在‘火炬’和春花腋芽中的时空表达模式存在显著差异。‘火炬’腋芽中开花基因FT2在2月27日和4月8日出现两次表达峰值,其他取样点表达较低,但在春花腋芽中,开花基因FT2在2月27日开始上调表达,3月9日达到表达峰值,3月底到4月初持续上调表达(图6),表明3月9日可能是春花腋花芽分化的关键期,FT2基因在春花腋芽中起重要作用;与‘火炬’相比,GIGANTEA-like和AP1基因在春花腋芽中显著上调表达,模式植物拟南芥中大量研究已证实AP1基因是花序分生组织向花分生组织转化过程中的关键基因及花器官决定的重要基因,而AP1基因在2月27日—3月19日一直在春花腋芽中保持较高水平(图6),表明2月27日—3月19日可能是春花腋芽转变为腋花芽的关键时期;开花促进因子FPF1-like基因在‘火炬’腋芽中表达水平较高,尤其是在2月27日和3月19日显著上调表达,但其在春花腋芽中始终保持较低表达水平,无明显波动(图6);2月份时早花基因ELF3-like的表达水平在‘火炬’和春花腋芽中无差异,但在3月9—30日ELF3-like基因在春花腋芽中显著下调表达(图6),暗示FPF1-like和ELF3-like基因的下调表达可能有利于春花腋芽成花;FY-like、FPA和LF等开花基因以及开花抑制基因FLC在‘火炬’和春花腋芽中的表达水平基本一致,无明显差异。ABA受体基因PYL4在‘火炬’腋芽中表达水平较低,在春花腋芽中显著上调表达(图6);ABA脱氢酶基因ABH4在‘火炬’腋芽中表达水平较高,在春花腋芽中显著下调表达(图6),表明PYL4基因的上调表达和ABH4基因的下调表达可能在春花腋芽成花中起重要调控作用。

表6 枇杷腋芽转录组筛选到的34个候选基因Table 6 Thirty-four candidate genes from transcriptome screening of loquat axillary buds

图6 开花基因和ABA信号基因在春花和‘火炬’腋芽中的表达水平Fig.6 Expression levels of flowering-related genes and ABA signaling genes in axillary buds of Chunhua and‘Huoju’

此外,从‘火炬’和春花腋芽中筛选到大量差异表达的功能未知基因,其中EVM0019086、EVM0037327、EVM0024514、EVM0036212、EVM0018511和Eriobotrya_japonica_newGene_7354等基因在春花腋芽显著上调表达,这些基因在‘火炬’腋芽中不表达或表达水平极低(图7);EVM0004793、EVM0025921和EVM0037211基因在春花腋芽中几乎不表达,但在‘火炬’腋芽中表达量极高,表明这些功能未知基因的上调或下调表达可能在春花腋芽成花过程中起重要调控作用。

图7 功能未知基因在春花和‘火炬’腋芽中的表达水平Fig.7 Expression levels of function-unknown genes in axillary buds of Chunhua and‘Huoju’

3 讨论

成花是植物生命周期中从营养生长阶段向生殖生长阶段转变的一个实质性标志,其过程主要包括花芽诱导、成花启动和花发育。同时,植物成花是一个高度复杂的生理进程,在外部和内部各种因素相互调节的网络下完成。内部因素的调节影响着植物成花的进程和成功与否,关键候选基因的筛选以及功能验证、分子遗传机理解析已成为成花研究的一个热点领域[24-26]。枇杷属于顶芽成花植物,有关枇杷成花机理的研究主要集中于顶芽成花机制解析,但关于枇杷腋芽的成花生理和分子遗传机制的相关研究鲜有报道,RNA测序(RNA-seq)作为转录组测序和差异表达基因筛选的重要手段,已被广泛应用于关键候选基因的筛选以及分子遗传机制的解析[17]。本研究以枇杷品系春花(腋芽可成花)和品种‘火炬’(腋芽不成花)为试验材料,共鉴定到66 197个差异表达基因,其中上调表达基因主要集中在2月27日—4月8日,维恩图表明在不同取样时间点同时上调表达的基因有160个,同时下调表达的基因有198个。陈晨等[27]以两个不同时期的杏花花芽为样本进行转录组测序,共筛选出6 850个显著差异表达基因,其中在花芽萌动期上调表达基因有2 784个,下调表达基因有4 066个,花芽萌动期特异表达基因有392个,盛花期特异表达基因有346个。李奥等[28]对不同分化时期的蝴蝶兰成花差异基因进行筛选与分析,所检测的16 181个基因中,2个花芽分化时期的蝴蝶兰共获得6 088个差异表达基因,其中上调基因有3 319个,下调基因有2 769个。本研究对6个取样时间点的枇杷腋芽进行转录组测序,筛选得到不同时期同时上调表达基因160个,下调基因198个,基因数量远低于上述研究结果。不同物种在不同处理下进行转录组测序筛选到的差异表达基因存在差异。

转录组分析是筛选成花途径中差异表达基因的有效途径,并在毛竹[29]、苹果[30]和菊花[31-33]等植物中广泛应用,并鉴定了不同开花阶段或与开花能力相关的关键候选基因。进一步对差异表达基因进行GO功能分类以及GO富集分析和KEGG通路显著性富集分析,深入解析差异基因的功能和参与的代谢途径。芒果不同花芽分化时期的转录组筛选到的差异表达基因主要富集于次生物质的生物合成和积累、糖代谢和光合作用等通路上[34]。植物开花过程中植物激素、糖类代谢以及次生物质的合成等途径对成花进程具有重要作用。本研究对枇杷不同成花时期腋芽中的差异基因进行GO和KEGG分析,代谢通路主要富集于植物激素信号传导、淀粉和糖代谢、苯基丙烷代谢、病原体互作反应和谷胱甘肽代谢等途径,研究结果与前期报道相一致,对该富集途径中关键候选基因进行功能验证和分子机制解析,可深入阐明植物成花调控网络。

大量研究对模拟植物拟南芥成花过程调控的分子遗传机理进行探讨,并且鉴定了多种调控路径,包括光周期途径、春化途径、温度途径和赤霉素途径等传递光、温度和激素等信号调控成花,同时这些途径通过整合一些重要的基因,实现对植物成花时间的精确调控。FT蛋白是长期以来寻找的成花诱导物质,成花素的主要成分,在茎顶端分生组织中,与一个basic leucine zipper(bzip)转录因子FD相互作用激活另一个开花整合基因SOC1和一个花分生组织特异基因APETALA1(AP1)表达启动花的发育[35-36]。FY、FPA、FLD基因为自主途径中的重要调控基因,可通过抑制MADS转录因子基因FLOWERING LOCUS C(FLC)的表达来促进开花。植物激素作为植物体内重要的调节物质,对植物的成花具有重要影响[37-38]。在拟南芥中,ABA作为植物的抑制因子,通过外源ABA的喷施,促进ABSCISIC ACID-INSENSITIVE MUTANT5(ABI5)基因的过表达和FLC基因表达量的上调,延迟植物的开花时间[39]。目前,通过对花期不同的枇杷材料花芽分化时期的叶片进行转录组测序,成功克隆到相关开花基因(EdFT、EdCO、EdSOC1、EdPIF4、EdFD1、EdFD2和EdSVP)[16],但有关枇杷腋芽成花的机理研究较少。本研究在枇杷腋芽转录组中筛选到与枇杷腋芽成花相关的差异表达基因(FT、FY、ELF4、AP2、FLK),分析在春花(腋芽可成花)和‘火炬’(腋芽不成花)中的时空表达模式,差异表达基因在春花中表达量显著高于‘火炬’,基因表达模式与开花时间呈现相一致的趋势,表明筛选的差异基因在枇杷腋芽成花过程中发挥重要作用。同时筛选到的ABA信号途径基因(PYL4和ABH4),表达水平在春花和‘火炬’中差异显著,PYL4基因的上调表达和ABH4基因的下调表达可能在春花腋芽成花中起重要调控作用。

4 结论

本研究从春花(腋芽可成花)和‘火炬’(腋芽不成花)腋芽转录组数据中共鉴定到66 197个差异表达基因,维恩图表明在不同取样时间点同时上调表达的基因有160个,同时下调表达的基因有198个;GO功能注释分子和KEGG通路显著性富集分析筛选到植物激素信号传导、淀粉和糖代谢、苯基丙烷代谢、病原体互作反应和谷胱甘肽代谢等重要代谢途径;FT、GIGANTEA-like、AP1、FPF1-like、ELF3-like基因和PYL4、ABH4 ABA信号途径相关差异基因时空表达模式与春花腋芽开花进程呈现相一致的趋势;同时,差异基因在枇杷品系春花中的表达量显著高于‘火炬’,表明差异表达基因在枇杷品系春花腋芽成花过程中发挥至关重要的作用,调控春花腋芽的成花进程。本研究有助于更好地了解枇杷腋芽成花过程中的基因转录水平,并可为阐明枇杷腋芽成花的分子遗传机制奠定理论基础。