牡丹花型基因的研究

周 洁，黄 婧，孙体如，张 忠，曹 瑶

(1.江苏省林业科学研究院，江苏 南京 211153; 2.沭阳县生产力促进中心，江苏 沭阳 223600)

牡丹(PaeoniasuffruticosaAndrews.)为芍药科芍药属落叶小灌木，是我国的十大传统名花之一。我国牡丹资源丰富，目前国内外共有牡丹 1 100余品种。中国栽培牡丹主要有中原品种群、西北品种群、西南品种群、江南品种群以及国外品种群[1]，在品种资源、优异性状和花型齐备等方面，我国居世界首位。花型是牡丹最重要的观赏价值之一，牡丹的花型有2大类(单花类、台阁花类)、4个亚类、13个类型[2]。牡丹千层类花型花瓣花芽分化的特点为花瓣向心式层层增加，楼子类花型分化特点为雄蕊离心式分化和瓣化，台阁类花型为同一花原基上分化出上下2朵重叠的花器官。

花器官特征基因决定花原基形成不同的花器官，从而形成形态多样的花。目前已确定4条开花调控途径，即光周期途径、 春化作用、 自主途径、赤霉素(GA3)途径[3]，通过调控植物开花时间的4条信号途径激活一套开花途径整合因子从而激活花分生组织特异基因启动开花。主要的开花途径整合因子有CO (CONSTANS)、FLOWERING LOCUS T (FT)、LEAFY (LFY) 和SUPPRESSOR OF OVEREXPRESSION OF CONSTANS 1 (SOC1)[4],SPL(squamosa promoter-binding protein-like)[5]，然后激活花分生组织特性基因如APETALA1/2 (AP1/2) 、FRUITFUL ( FUL )、SEEPALLATA4 (SEP4)等使茎端分生组织转化为花序分生组织[6]。AP1/2属A功能基因，决定萼片、花瓣特性, APETALA3(AP3)，TM6，PI属B功能基因[7-8]，决定花瓣、雄蕊特性；AGAMOUS (AG)是唯一C功能基因，早期调控雌、雄蕊的发育，后期决定细胞分化；SHATTERPROOF(SHP)属D功能基因，控制心皮发育；WUS为顶端分生组织特征基因，维持AG活性。

重瓣性是牡丹最要的观赏价值之一，也是现代牡丹育种的重要目标之一。重瓣牡丹品种已占据整个牡丹品种的80%以上。目前牡丹中与花器官发育的基因已有被分离和克隆，但是仅限于个别基因的克隆，关于牡丹花型调控的基因尚未系统研究。本研究利用Illumina HiSeq技术对牡丹2个不同花型的品种雌雄蕊正常的‘黑海撒金’和雌蕊瓣化的‘璎珞宝珠’进行测序，以期找到不同花型的差异表达基因，为花型的分子研究提供依据。研究牡丹花型调控的相关基因, 对深入了解牡丹花型发育机制具有重要理论意义。

1 材料与方法

1.1 试验植物

本试验选用牡丹品种‘黑海撒金’(tree peony ‘HEIHAISAJIN’)和‘璎珞宝珠’(‘YINGLUOBAOZHU’)总花瓣为材料，采集于江苏省林业科学研究院牡丹种质资源圃，于4月待完全开花后采集2种牡丹的花瓣，每个品种采集3份，样品迅速冻于液氮中，放于-80 ℃保存，备用。

1.2 转录组De novo测序

本研究利用无参转录组对‘黑海撒金’和‘璎珞宝珠’的花瓣(各3份)进行de novo测序，花瓣的总RNA采用TRIzol® Reagent(Invitrogen)提取后，利用Illumina Truseq TM RNA sample prep Kit方法构建文库，通过TBS380(Picogreen)定量后使用Illumina HiSeq进行测序，对测序得到的原始数据进行过滤，去除低质量、接头污染以及N比例超过10%的序列，使用的软件为SeqPrep (https:∥github.com/jstjohn/SeqPrep), Sickle (https:∥github. com/najoshi/sickle)。将过滤后的干净序列(clean reads)使用Trinity技术进行de novo组装，然后将组装的转录本进行聚类去冗余，得到unigene。

1.3 功能注释与差异表达基因的筛选

对拼接得到的unigene序列进行注释，并与数据库NR，Swiss-Prot，KEGG和COG比对(E value<10-5)。对每个转录本的表达量进行样本间的表达差异显著性分析，找到相应的差异表达基因，并对差异表达进行可视化分析。差异表达基因的GO功能分析使用软件Goatools (https:∥github.com/tanghaibao/GOatools) 进行富集分析，使用方法为Fisher精确检验，p值(p_fdr)≤0.05。差异表达基因的KEGG分析使用KOBAS( http:∥kobas.cbi.pku.edu.cn/home.do ) 进行KEGG PATHWAY富集分析，校正的p值(Corrected P-Value)以0.05为阈值。

2 结果与分析

2.1 转录组序列测序和拼接

通过对‘黑海撒金’和‘璎珞宝珠’2个品种6个样本进行RNA-Seq转录组测序，获得原始数据后，经过过滤、去接头、去冗余和低质量序列后，获得33 186 290到48 461 266条高质量序列clean reads，Q20的比例为98.46%—98.6%，结果表明测序质量较好(见表1)。通过trinity软件(http:∥trinityrnaseq. sourceforge.net/)进行de novo拼接，组装后获得49 191个unigenes，平均长度为753 bp，最长的为10 070 bp，N50为835 bp(见表2)。

表1 牡丹转录组测序组装

表2 de novo序列拼接

2.2 unigenes功能注释

如图1所示，将拼接得到的unigenes进行注释，并与NR，GO，KEGG，SWSS数据库进行比对，其中比对到NR的unigenes为31 994，比对到GO数据库的为10 451个unigenes，比对到COG数据库的为24 605个unigenes，比对到KEGG的为13 279的unigenes。这4个数据库中共有的基因为5 196个unigenes。

图1 unigenes的数据库比对和注释

2.3 差异表达基因DEGs的筛选

采用RPKM的方法，将差异表达基因定义为FDR<0.001,且倍数差异在2倍以上(|Log2 |Ratio|>1)的基因。以散点图表示1个基因在2个转录本中的表达量，横坐标为在‘黑海撒金’的表达量，纵坐标为在‘璎珞宝珠’中的表达量(见图2-a)；火山图(Volcano-plots)为基因或转录本在2个样本间表达差异的倍数变化值，横坐标为倍数的对数，纵坐标为P值的对数(见图2-b)。通过‘黑海撒金’和‘璎珞宝珠’2个转录组的表达量差异的比较，共得到2 735 个差异表达基因，其中1 082个为上调基因，1 653个差异表达基因下调(如图2)。

注：a散点图；b 火山图。黑为‘黑海撒金’，璎为‘璎珞宝珠’; Log2(FC)为‘黑海撒金’和‘璎珞宝珠’的表达量比值的对数值。图2 差异表达基因可视化分析

2.4 差异表达基因的注释

将差异表达基因(DEGs)进行GO 功能注释，主要分布在25个功能区，在细胞组分中显著富集的为90 S preribosome， 在分子功能中最显著富集的为taxane 10-beta-hydroxylase activity，在生物学过程中最为显著的是response to wounding。其中ribosomal subunit中富集的基因最多，有14个unigenes(见图3-a)。对DEG进行KEGG代谢通路的分析，显著富集的通路有13条(P<0.05), 前5个富集的通路为Ribosome，Chemical carcinogenesis，Drug metabolism-cytochrome P450，Oxytocin signaling pathway，Tyrosine metabolism, 其中Ribosome最为显著富集且基因最多，有36个(见图3-b)。

注：a GO功能富集图；b KEGG路径富集图图3 DEGs的GO功能和KEGG路径富集图

2.5 候选功能基因的挖掘

如表3所示,对功能基因进行挖掘，将得到的DEGs序列通过与NCBI进行比对，发现与开花相关的功能基因有2个，为agamous-like MADS-box protein 65(AGL65)和floral homeotic protein APETALA 2(AP2)。AP2在‘黑海撒金’中几乎无表达，在‘璎珞宝珠’中表达量稍高。 AGL65在‘黑海撒金’中表达量稍高，在‘璎珞宝珠’中几乎无表达。同时还发现了与性别决定相关控制雌雄发育的基因EMBRYO SAC DEVELOPMENT ARREST 3(EDA3), 该基因在‘黑海撒金’中表达量较高，而在‘璎珞宝珠’中几乎无表达。说明在牡丹中雌蕊的瓣化由开花基因和性别决定基因共同调控。

表3 功能基因的筛选

3 结论与讨论

目前从牡丹萼片、花瓣、雄蕊和心皮中已分离出很多关于花发育的相关基因。其中与ABCDE模型相关的基因有PsAP1，PsAP2，PsPI，PsMADS1，PsMADS9，PsAG和PsMADS5等[9]。MADS-box基因作为参与花器官发育的一类重要转录因子，在生长发育、信号传导、开花调节和器官发育方面也发挥着重要的作用[10-11]。AP2属于A类基因，主要负责萼片和花瓣的形成[12]。本研究中AP2基因在‘黑海撒金’中几乎无表达，在‘璎珞宝珠’中表达量稍高，说明AP2参与了雌蕊瓣化的过程，负责花瓣的形成，与ABCDE模型的理论一致。AGL6基因已在柳杉[13]、买麻藤[14]、春兰[15]等植物中分离得到。水稻AGL6同源基因OsMADS6在雌蕊与分生组织中表达[16]。腊梅CpAGL6可能具有调控开花时间与花器官形成的双重功能[17]。牡丹中PsAGL6在蔷薇型牡丹中的表达量较高，在萼片、花瓣和雌蕊显著表达[18]。而在本研究中AGL65在‘黑海撒金’中表达量稍高，在‘璎珞宝珠’中几乎无表达。说明AGL65参与了雌蕊的分化和发育，而在雌蕊瓣化的花型中无表达，且AGL65与AGL6功能存在差异，AGL65成员目前在功能方面还没有报道，所以AGL65可能是新的调节雌蕊瓣化发育的基因成员，属于C类基因。EMBRYO SAC DEVELOPMENT ARREST 3 (EDA3)参与胚囊的分化, 且仅在子房中表达，且EDA3能促进胚珠的发育[19]。EDA3也仅在‘黑海撒金’中表达，表达量为7倍，所以这一结果也与报道的功能一致。所以可以推测EDA3为D类基因，调控胚珠的发育。综上所述，牡丹花型雌蕊的瓣化不仅有ABCDE开花模型A类基因AP2促进了花瓣的形成，C类基因AGL65调控雌蕊的发育，也可能由D类基因EDA3作用于胚珠，调控雌蕊的发育。本研究利用转录组测序技术研究了不同花型牡丹的差异表达基因，筛选出了调控牡丹花发育的新成员基因AGL65，AP2和EDA3，这些基因的挖掘为牡丹花器官的研究包括牡丹雌蕊发育和瓣化、胚珠发育的研究提供了很好的理论基础。