周 潔，黃 婧，孫體如，張 忠，曹 瑤

(1.江蘇省林業(yè)科學(xué)研究院，江蘇 南京 211153; 2.沭陽縣生產(chǎn)力促進(jìn)中心，江蘇 沭陽 223600)

牡丹(PaeoniasuffruticosaAndrews.)為芍藥科芍藥屬落葉小灌木，是我國的十大傳統(tǒng)名花之一。我國牡丹資源豐富，目前國內(nèi)外共有牡丹 1 100余品種。中國栽培牡丹主要有中原品種群、西北品種群、西南品種群、江南品種群以及國外品種群[1]，在品種資源、優(yōu)異性狀和花型齊備等方面，我國居世界首位?；ㄐ褪悄档ぷ钪匾挠^賞價值之一，牡丹的花型有2大類(單花類、臺閣花類)、4個亞類、13個類型[2]。牡丹千層類花型花瓣花芽分化的特點為花瓣向心式層層增加，樓子類花型分化特點為雄蕊離心式分化和瓣化，臺閣類花型為同一花原基上分化出上下2朵重疊的花器官。

花器官特征基因決定花原基形成不同的花器官，從而形成形態(tài)多樣的花。目前已確定4條開花調(diào)控途徑，即光周期途徑、 春化作用、 自主途徑、赤霉素(GA3)途徑[3]，通過調(diào)控植物開花時間的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)等使莖端分生組織轉(zhuǎn)化為花序分生組織[6]。AP1/2屬A功能基因，決定萼片、花瓣特性, APETALA3(AP3)，TM6，PI屬B功能基因[7-8]，決定花瓣、雄蕊特性；AGAMOUS (AG)是唯一C功能基因，早期調(diào)控雌、雄蕊的發(fā)育，后期決定細(xì)胞分化；SHATTERPROOF(SHP)屬D功能基因，控制心皮發(fā)育；WUS為頂端分生組織特征基因，維持AG活性。

重瓣性是牡丹最要的觀賞價值之一，也是現(xiàn)代牡丹育種的重要目標(biāo)之一。重瓣牡丹品種已占據(jù)整個牡丹品種的80%以上。目前牡丹中與花器官發(fā)育的基因已有被分離和克隆，但是僅限于個別基因的克隆，關(guān)于牡丹花型調(diào)控的基因尚未系統(tǒng)研究。本研究利用Illumina HiSeq技術(shù)對牡丹2個不同花型的品種雌雄蕊正常的‘黑海撒金’和雌蕊瓣化的‘瓔珞寶珠’進(jìn)行測序，以期找到不同花型的差異表達(dá)基因，為花型的分子研究提供依據(jù)。研究牡丹花型調(diào)控的相關(guān)基因, 對深入了解牡丹花型發(fā)育機(jī)制具有重要理論意義。

1 材料與方法

1.1 試驗植物

本試驗選用牡丹品種‘黑海撒金’(tree peony ‘HEIHAISAJIN’)和‘瓔珞寶珠’(‘YINGLUOBAOZHU’)總花瓣為材料，采集于江蘇省林業(yè)科學(xué)研究院牡丹種質(zhì)資源圃，于4月待完全開花后采集2種牡丹的花瓣，每個品種采集3份，樣品迅速凍于液氮中，放于-80 ℃保存，備用。

1.2 轉(zhuǎn)錄組De novo測序

本研究利用無參轉(zhuǎn)錄組對‘黑海撒金’和‘瓔珞寶珠’的花瓣(各3份)進(jìn)行de novo測序，花瓣的總RNA采用TRIzol? Reagent(Invitrogen)提取后，利用Illumina Truseq TM RNA sample prep Kit方法構(gòu)建文庫，通過TBS380(Picogreen)定量后使用Illumina HiSeq進(jìn)行測序，對測序得到的原始數(shù)據(jù)進(jìn)行過濾，去除低質(zhì)量、接頭污染以及N比例超過10%的序列，使用的軟件為SeqPrep (https:∥github.com/jstjohn/SeqPrep), Sickle (https:∥github. com/najoshi/sickle)。將過濾后的干凈序列(clean reads)使用Trinity技術(shù)進(jìn)行de novo組裝，然后將組裝的轉(zhuǎn)錄本進(jìn)行聚類去冗余，得到unigene。

1.3 功能注釋與差異表達(dá)基因的篩選

對拼接得到的unigene序列進(jìn)行注釋，并與數(shù)據(jù)庫NR，Swiss-Prot，KEGG和COG比對(E value<10-5)。對每個轉(zhuǎn)錄本的表達(dá)量進(jìn)行樣本間的表達(dá)差異顯著性分析，找到相應(yīng)的差異表達(dá)基因，并對差異表達(dá)進(jìn)行可視化分析。差異表達(dá)基因的GO功能分析使用軟件Goatools (https:∥github.com/tanghaibao/GOatools) 進(jìn)行富集分析，使用方法為Fisher精確檢驗，p值(p_fdr)≤0.05。差異表達(dá)基因的KEGG分析使用KOBAS( http:∥kobas.cbi.pku.edu.cn/home.do ) 進(jìn)行KEGG PATHWAY富集分析，校正的p值(Corrected P-Value)以0.05為閾值。

2 結(jié)果與分析

2.1 轉(zhuǎn)錄組序列測序和拼接

通過對‘黑海撒金’和‘瓔珞寶珠’2個品種6個樣本進(jìn)行RNA-Seq轉(zhuǎn)錄組測序，獲得原始數(shù)據(jù)后，經(jīng)過過濾、去接頭、去冗余和低質(zhì)量序列后，獲得33 186 290到48 461 266條高質(zhì)量序列clean reads，Q20的比例為98.46%—98.6%，結(jié)果表明測序質(zhì)量較好(見表1)。通過trinity軟件(http:∥trinityrnaseq. sourceforge.net/)進(jìn)行de novo拼接，組裝后獲得49 191個unigenes，平均長度為753 bp，最長的為10 070 bp，N50為835 bp(見表2)。

表1 牡丹轉(zhuǎn)錄組測序組裝

表2 de novo序列拼接

2.2 unigenes功能注釋

如圖1所示，將拼接得到的unigenes進(jìn)行注釋，并與NR，GO，KEGG，SWSS數(shù)據(jù)庫進(jìn)行比對，其中比對到NR的unigenes為31 994，比對到GO數(shù)據(jù)庫的為10 451個unigenes，比對到COG數(shù)據(jù)庫的為24 605個unigenes，比對到KEGG的為13 279的unigenes。這4個數(shù)據(jù)庫中共有的基因為5 196個unigenes。

圖1 unigenes的數(shù)據(jù)庫比對和注釋

2.3 差異表達(dá)基因DEGs的篩選

采用RPKM的方法，將差異表達(dá)基因定義為FDR<0.001,且倍數(shù)差異在2倍以上(|Log2 |Ratio|>1)的基因。以散點圖表示1個基因在2個轉(zhuǎn)錄本中的表達(dá)量，橫坐標(biāo)為在‘黑海撒金’的表達(dá)量，縱坐標(biāo)為在‘瓔珞寶珠’中的表達(dá)量(見圖2-a)；火山圖(Volcano-plots)為基因或轉(zhuǎn)錄本在2個樣本間表達(dá)差異的倍數(shù)變化值，橫坐標(biāo)為倍數(shù)的對數(shù)，縱坐標(biāo)為P值的對數(shù)(見圖2-b)。通過‘黑海撒金’和‘瓔珞寶珠’2個轉(zhuǎn)錄組的表達(dá)量差異的比較，共得到2 735 個差異表達(dá)基因，其中1 082個為上調(diào)基因，1 653個差異表達(dá)基因下調(diào)(如圖2)。

注：a散點圖；b 火山圖。黑為‘黑海撒金’，瓔為‘瓔珞寶珠’; Log2(FC)為‘黑海撒金’和‘瓔珞寶珠’的表達(dá)量比值的對數(shù)值。圖2 差異表達(dá)基因可視化分析

2.4 差異表達(dá)基因的注釋

將差異表達(dá)基因(DEGs)進(jìn)行GO 功能注釋，主要分布在25個功能區(qū)，在細(xì)胞組分中顯著富集的為90 S preribosome， 在分子功能中最顯著富集的為taxane 10-beta-hydroxylase activity，在生物學(xué)過程中最為顯著的是response to wounding。其中ribosomal subunit中富集的基因最多，有14個unigenes(見圖3-a)。對DEG進(jìn)行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所示,對功能基因進(jìn)行挖掘，將得到的DEGs序列通過與NCBI進(jìn)行比對，發(fā)現(xiàn)與開花相關(guān)的功能基因有2個，為agamous-like MADS-box protein 65(AGL65)和floral homeotic protein APETALA 2(AP2)。AP2在‘黑海撒金’中幾乎無表達(dá)，在‘瓔珞寶珠’中表達(dá)量稍高。 AGL65在‘黑海撒金’中表達(dá)量稍高，在‘瓔珞寶珠’中幾乎無表達(dá)。同時還發(fā)現(xiàn)了與性別決定相關(guān)控制雌雄發(fā)育的基因EMBRYO SAC DEVELOPMENT ARREST 3(EDA3), 該基因在‘黑海撒金’中表達(dá)量較高，而在‘瓔珞寶珠’中幾乎無表達(dá)。說明在牡丹中雌蕊的瓣化由開花基因和性別決定基因共同調(diào)控。

表3 功能基因的篩選

3 結(jié)論與討論

目前從牡丹萼片、花瓣、雄蕊和心皮中已分離出很多關(guān)于花發(fā)育的相關(guān)基因。其中與ABCDE模型相關(guān)的基因有PsAP1，PsAP2，PsPI，PsMADS1，PsMADS9，PsAG和PsMADS5等[9]。MADS-box基因作為參與花器官發(fā)育的一類重要轉(zhuǎn)錄因子，在生長發(fā)育、信號傳導(dǎo)、開花調(diào)節(jié)和器官發(fā)育方面也發(fā)揮著重要的作用[10-11]。AP2屬于A類基因，主要負(fù)責(zé)萼片和花瓣的形成[12]。本研究中AP2基因在‘黑海撒金’中幾乎無表達(dá)，在‘瓔珞寶珠’中表達(dá)量稍高，說明AP2參與了雌蕊瓣化的過程，負(fù)責(zé)花瓣的形成，與ABCDE模型的理論一致。AGL6基因已在柳杉[13]、買麻藤[14]、春蘭[15]等植物中分離得到。水稻AGL6同源基因OsMADS6在雌蕊與分生組織中表達(dá)[16]。臘梅CpAGL6可能具有調(diào)控開花時間與花器官形成的雙重功能[17]。牡丹中PsAGL6在薔薇型牡丹中的表達(dá)量較高，在萼片、花瓣和雌蕊顯著表達(dá)[18]。而在本研究中AGL65在‘黑海撒金’中表達(dá)量稍高，在‘瓔珞寶珠’中幾乎無表達(dá)。說明AGL65參與了雌蕊的分化和發(fā)育，而在雌蕊瓣化的花型中無表達(dá)，且AGL65與AGL6功能存在差異，AGL65成員目前在功能方面還沒有報道，所以AGL65可能是新的調(diào)節(jié)雌蕊瓣化發(fā)育的基因成員，屬于C類基因。EMBRYO SAC DEVELOPMENT ARREST 3 (EDA3)參與胚囊的分化, 且僅在子房中表達(dá)，且EDA3能促進(jìn)胚珠的發(fā)育[19]。EDA3也僅在‘黑海撒金’中表達(dá)，表達(dá)量為7倍，所以這一結(jié)果也與報道的功能一致。所以可以推測EDA3為D類基因，調(diào)控胚珠的發(fā)育。綜上所述，牡丹花型雌蕊的瓣化不僅有ABCDE開花模型A類基因AP2促進(jìn)了花瓣的形成，C類基因AGL65調(diào)控雌蕊的發(fā)育，也可能由D類基因EDA3作用于胚珠，調(diào)控雌蕊的發(fā)育。本研究利用轉(zhuǎn)錄組測序技術(shù)研究了不同花型牡丹的差異表達(dá)基因，篩選出了調(diào)控牡丹花發(fā)育的新成員基因AGL65，AP2和EDA3，這些基因的挖掘為牡丹花器官的研究包括牡丹雌蕊發(fā)育和瓣化、胚珠發(fā)育的研究提供了很好的理論基礎(chǔ)。