普通油茶花發(fā)育過程數(shù)字表達譜測序及分析

胡玉玲，姚小華，任華東，王開良，張山

(1.銅仁學院農(nóng)林工程與規(guī)劃學院，貴州 銅仁 554300;2.湖南應(yīng)用技術(shù)學院, 湖南 常德 415000;3.中國林業(yè)科學研究院亞熱帶林業(yè)研究所, 浙江 杭州 311400)

高等植物開花是植物生命活動的重要現(xiàn)象，這不僅關(guān)系著該物種的延續(xù)，而且與農(nóng)業(yè)生產(chǎn)密切相關(guān)[1]。植物成花的整個過程中較關(guān)鍵的是成花誘導、花芽分化和花器官形成，這些過程直接決定了植物生育周期長短，與植物的產(chǎn)量及品質(zhì)關(guān)系緊密[2]。油茶(CamelliaoleiferaAbel.)是典型多花而不高產(chǎn)的物種，一方面由于油茶大部分品種自交可育性低，另一方面，油茶生產(chǎn)周期接近1年，茶果經(jīng)歷了4個物候差異非常明顯的季節(jié)，特別是秋冬季低溫潮濕和夏季的高溫少雨。因此培育出自交可育性高的油茶品種是當務(wù)之急，由于油茶是多年生木本植物，育種周期較長，可見培養(yǎng)優(yōu)良新品種非一蹴而就的事情。有研究表明，油茶花期遇到寒冷潮濕天氣會嚴重影響油茶產(chǎn)量，同時冬天溫度過低也會影響油茶幼果過冬從而影響油茶產(chǎn)量[3]。

基因工程是較有效的植物遺傳改良方式，尤其對多年生木本植物。針對上述影響油茶產(chǎn)量限定因子實施基因工程，需要明確并解決如下幾個問題：首先控制油茶成花轉(zhuǎn)變的關(guān)鍵基因，以有效實現(xiàn)人工成花轉(zhuǎn)變，縮短育種周期；其次控制花器官形成基因，如與雄蕊發(fā)育有關(guān)基因，以人工獲得雄性不育系；再次影響花開放進程相關(guān)基因，以改變油茶花期或者通過栽培措施改變花期以避開不利天氣，其中揭示油茶花發(fā)育分子機制最為重要。所有生物的整個生命活動就是不同編碼蛋白的基因進行選擇性開啟和關(guān)閉結(jié)果，為準確獲得需要遺傳改良成花相關(guān)基因，必須了解基因表達過程。細胞的特定發(fā)育階段，一般只有15%的基因進行了表達，這種在生物個體發(fā)育的不同階段、不同組織器官中，按照一定時空順序，基因的有序表達過程，就是基因的差別表達(Differential expressing)[4-5]。采用Solexa技術(shù)對植物樣本中數(shù)以百萬計的cDNA標簽進行序列測定，能夠?qū)φ麄€轉(zhuǎn)錄組和表達譜進行數(shù)字化的分析，數(shù)字的列信號可以準確、特異地反映對應(yīng)基因的真實表達情況。這種技術(shù)甚至可以精確地檢測低至一兩個拷貝的稀有轉(zhuǎn)錄本，并精確定量高達十萬個拷貝的轉(zhuǎn)錄本的表達量變化。由于序列無需事先設(shè)計，數(shù)字表達譜數(shù)據(jù)具有極佳的實時性，可以充分利用當前爆發(fā)式增長的信息資源，并與未來研究相銜接，數(shù)字表達譜可以檢測到許多未曾注釋的基因和基因組部位，為新基因的發(fā)現(xiàn)提供了良好的線索，然后通過生物信息學的分析及樣品間基因表達的比對可以初步判定基因表達過程[6-8]。

因此，本研究選擇普通油茶花芽分化前、萼片形成期、花瓣形成期、雌雄蕊形成期、子房及花藥形成期和雌雄蕊成熟期的花芽作為試驗材料，對該6個時期樣品進行表達譜測序并進行分析，以期揭示普通油茶花發(fā)育過程的分子機制及成花關(guān)鍵基因，為培育自交可育性高的油茶品種提供強有力的理論依據(jù)。

1 試驗材料與方法

1.1 試驗材料

普通油茶長林4號(Camelliaoleifera‘ChangLin 4’)優(yōu)良無性系。

1.2 采樣方法

根據(jù)相關(guān)文獻報道[9-11]，分別在油茶花芽分化前(5月8日)，萼片形成期(5月28日)、花瓣形成期(6月15日)、雌雄蕊形成期(7月3日)、子房和花藥形成期(7月25日)和雌雄蕊成熟期(8月25日)采集花芽，樣品采集后都在液氮帶回保存，放入-80℃的超低溫冰箱中保存?zhèn)溆谩?/p>

1.3 RNA提取

RNA提取采用艾德萊公司EASYspin Plus植物RNA快速提取試劑盒(目錄號為RN38)使用手冊進行，采用微量分光光度計Nanodrop 2000及瓊脂糖凝膠電泳檢測RNA質(zhì)量，保證RNA樣品濃度400ng/μL，28S︰18S大于1.8。

1.4 上機測序前準備和方法

利用富含ploy(T)低吸附磁珠對Total RNA中的mRNA進行富集純化處理，在高溫條件下，利用2價陽離子打斷mRNA以選取合適大小的目的片段，利用反轉(zhuǎn)錄酶和隨機引物將打斷后mRNA片段反轉(zhuǎn)錄形成cDNA第一條鏈，然后加入緩沖液、dNTPs，在RNase H和DNA polymerase I的作用下合成第二條cDNA鏈。對反轉(zhuǎn)錄合成的雙鏈cDNA進行末端修復、3’末端加A、連接接頭，通過瓊脂糖凝膠電泳，篩選一定范圍大小的片段。PCR擴增，上機測序。

1.5 生物信息學分析

用The BLAST-Like Alignment Tool對參考基因序列比對，基因表達定量分析，表達量用RPKM(Reads Per kb per Million reads)表示，RPKM值的計算公式,RPKM=(map the gene reads)/〔map all gene reads(million)×the gene length(kb)〕。

用IDEG6進行差異表達分析，用Cluster 3.0[12-13]進行基因聚類分析，從參考基因注釋信息里面提取出差異表達基因的信息，同時對差異表達基因做GO和KEGG Pathway的顯著性富集分析。通過GO功能顯著性富集分析能確定差異表達基因行使的主要生物學功能；通過Pathway顯著性富集能確定差異表達基因參與的主要生化代謝途徑和信號轉(zhuǎn)導途徑。整個過程見圖1。

圖1 數(shù)字表達譜生物信息分析流程圖

2 結(jié)果與分析

2.1 表達譜測序

2.1.1 測序數(shù)據(jù)統(tǒng)計

采用Solexa HiSeq 2000對每個樣品進行單端測序，統(tǒng)計測序數(shù)據(jù)產(chǎn)出、數(shù)據(jù)質(zhì)量以及cycleQ20%等情況。從表1可以看出6個發(fā)育時期的測序數(shù)據(jù)數(shù)量都超過8M，完全滿足后續(xù)的分析，E2樣品的GC含量最低(45.19%)，E5樣品的GC含量最高(54%)，測序的cycleQ20%都達到100%，保證測序結(jié)果準確性。

表1 測序數(shù)據(jù)

注：E1為 5月9日，E2為7月3日，E3為8月25日，E4為5月29日，E5為6月15日，E6為7月25日，下同。

2.1.2 各樣品質(zhì)量評估及相關(guān)基因表達特點

測序數(shù)據(jù)質(zhì)量評估分別有cDNA片段隨機性基因和測序飽和度分析。評價mRNA打斷隨機程度一般是以reads在參考基因上的分布情況來確定的。由于不同參考基因有不同的長度，必須先標準化reads在參考基因上的位置(reads在基因上的位置與基因長度的比值)，然后統(tǒng)計基因的不同位置比對上的reads數(shù)。如果打斷隨機性好，reads在基因各部位應(yīng)分布得比較均勻。從圖2 cDNA打斷可以看出，測序的reads數(shù)沒有3’或5’嚴重的偏向性，就可以保證結(jié)果的準確。

E1E4E5E2E6E3

圖2cDNA片段隨機性檢驗

Fig.2 Random testing of cDNA fragments

樣品中測序飽和度分析是用來衡量一個樣品的測序量多少的標準，隨著測序量(reads數(shù)量)的增多，檢測到的基因數(shù)也隨之上升，當測序量達到某個值時，其檢測到的基因數(shù)增長速度趨于平緩，說明檢測到的基因數(shù)趨于飽和。從圖3可以看出，各個樣品中檢測到的基因數(shù)基本趨于飽和，說明6個樣品的測序量已經(jīng)覆蓋了絕大部分表達的基因，基本可以挖掘全部的差異基因。

E1E4E5E2E6E3

圖3測序飽和度分析

Fig.3 Saturation analysis of sequencing

基因表達定量分析是根據(jù)各樣品比對到參考基因序列上的reads數(shù)，來計算各樣品的基因表達豐度信息，利用RPKM值來反映樣品的基因表達豐度(圖4)?；虮磉_豐度分布分析是根據(jù)各樣品中基因的RPKM值對各基因的表達豐度進行估算，從整體水平了解樣品中基因的表達模式。一般來說，不均一性、冗余性是細胞中mRNA的顯著特征，少量種類的基因表達豐度極高，而大部分種類的基因表達水平很低甚至極低。從6個樣品的表達基因中，RPKM值>100的高表達基因數(shù)約占總基因數(shù)量的1.5%，而此區(qū)間內(nèi)基因的總表達豐度約占全部基因表達豐度的75%；從基因表達量分布圖(圖5)可以看出，6個樣品中基因表達的RPKM值在0-5區(qū)間內(nèi)的低表達基因數(shù)約占總基因數(shù)量的 70%，而此區(qū)間內(nèi)基因的總表達豐度約占全部基因表達豐度的5%；其表達模式符合 mRNA不均一性、冗余性的特征。從圖5中的E5和E6可以看出RPKM值小于100的非常少，絕大部分基因表達量都大于100。

E1E4E5E2E6E3

圖4 基因表達分布圖

圖5基因表達量分布圖

Fig.5 Distribution of gene expression level

2.1.3 比對效率統(tǒng)計

由于參考基因序列是轉(zhuǎn)錄組部分構(gòu)建的Unigenes庫序列，表達譜測序的樣品通過統(tǒng)計比對效率，可以檢查數(shù)據(jù)是否正常，這也是后續(xù)基因表達及差異分析的結(jié)果準確性的保證。通過使用Blat比對軟件，將各樣品測序得到的reads分別與參考基因序列進行比對，樣品的Mapped percentage(表2)均超過60%，因此可以判定比對效率較高，測序結(jié)果可以進行后續(xù)分析。

表2 樣品與參考基因的比對效率

2.2 差異表達基因分析

差異表達基因分析見表3。

表3 兩兩樣品間比較差異表達基因數(shù)量統(tǒng)計

從表3可以看出，5月9日(E1)和8月25日(E3)的差異表達基因明顯較少，可見這兩個時期油茶都在進行最基本并相近的生命活動，各個時期和7月3日(E2)與6月15日(E5)差異基因最多，其中上調(diào)基因也最多，達到19 408個，其次是5月29日(E4)與6月15日(E5)，E4對E5的上調(diào)基因有16 644，下調(diào)僅107個，可見這兩個時期生命代謝活動最為活躍。

2.3 差異表達基因生物信息學分析

根據(jù)測序Reads與參考基因序列的比對，得到對應(yīng)基因在樣品中的表達豐度，基因在不同樣品間的表達豐度的不同，尋找差異表達的基因。結(jié)合差異基因的功能注釋，進行pathway(KEGG)與GO的富集分析，以及對差異表達基因進行模式聚類的分析。富集分析采用fisher精確檢驗，通過Bonferroni校正法進行校正，得到差異基因顯著富集的pathway和GO功能類，以便進行進一步的分析研究。

2.3.1 花發(fā)育6個時期樣品差異基因聚類分析

根據(jù)基因表達豐度，使用cluster聚類軟件對所有樣品間的差異表達基因做層次聚類分析，將具有相同或相似表達行為的基因進行聚類。從圖6各列聚類結(jié)果可以看出，在所有的差異基因中絕大部分基因表達行為相似，只有少部分基因扮演著關(guān)鍵的生物功能。從不同樣品聚類來看，E1、E3(5月9日和8月25日)和E5、E6(6月15日和7月25 日)基因表達水平較接近，E2、E5(6月15日和7月3日)基因表達水平差異最為明顯。

圖6 差異基因聚類圖

注：圖中列代表不同的樣品；行代表不同的基因；相應(yīng)的顏色代表基因在樣品中表達豐度的高低。黃色表示基因在樣品中高表達，藍色表示基因在樣品中低表達

Fig.6 Differences in gene clustering

2.3.2 花發(fā)育6個時期樣品差異表達基因生物功能分析

從圖7-A可以看出，執(zhí)行細胞組成功能中細胞外基質(zhì)和細胞核與差異基因關(guān)系緊密；分子功能中核酸和蛋白綁定轉(zhuǎn)錄因子活性與差異基因關(guān)系緊密；生物學過程中的抗氧化性，細胞增殖，病毒繁殖和各種生命運動與差異基因關(guān)系緊密。

從圖7-B可以看出，現(xiàn)較多基因的生命功能有翻譯，核糖體結(jié)構(gòu)和生物合成為245條，轉(zhuǎn)錄因子為372條，復制，重組和修復為320條，細胞信號轉(zhuǎn)導機制為285條，翻譯后修飾，蛋白質(zhì)周轉(zhuǎn)，伴侶蛋白為283條，能源生產(chǎn)和轉(zhuǎn)換為114條，碳水化合物運輸和代謝為136條，氨基酸的轉(zhuǎn)運和代謝為122條。

AB

圖7差異基因生物功能分析

注：A為GO的注釋圖，B為Cog分類圖

Fig.7 Biological attributes analysis on differential gene

2.4 花器官形成相關(guān)unigenes表達

花器官形成是植物成花的重要過程，ABC模式是目前花器官形成最成熟的理論，從表4可以看出，A類基因Unigene數(shù)量最多，其中Unigene67783在不同的發(fā)育階段差異較明顯，表達水平也較高；B類基因有9條Unigene，其中Unigene56059 表達水平較高，各發(fā)育階段差異也較明顯；屬于C類基因的Unigene數(shù)量有9條，總體表達都不高，差異不明顯。

3 結(jié)論與討論

3.1 數(shù)字表達譜測序

數(shù)字基因表達譜測序方法是直接測定每個基因特異性的表達標簽序列，以計數(shù)表達標簽序列數(shù)目來確定該基因的表達量，從而大大提高了定量分析的準確性[14-15]。采用Solexa方法測序，對于那樣表達差異不大的基因，也能夠檢測其表達的差異性，同時也能夠發(fā)現(xiàn)新的轉(zhuǎn)錄本或者基因組表達調(diào)控區(qū)域等。另外該技術(shù)不用事先設(shè)計探針，采用直接測序的方式，因此即使不了解物種基因信息，也可以直接對該物種進行測序，包括未知基因在內(nèi)的全基因組表達譜分析。從本文來看，從花發(fā)育的6個時期樣品測序reads與轉(zhuǎn)錄組部分測序構(gòu)建的Unigene庫的序列對比[16]，對比效率都超過60%，通過生物信息學的分析更有效地了解花發(fā)育的基因表達。本次試驗采樣嚴格把握花發(fā)育的不同階段，在高通量轉(zhuǎn)錄組測序數(shù)據(jù)基礎(chǔ)上，進行表達譜測序，本文表達譜測序數(shù)據(jù)量最少者超過8M，最多數(shù)據(jù)量近30M。因此，表達譜數(shù)據(jù)基本覆蓋了所有表達標簽，對準確后續(xù)研究及生物信息學分析奠定了基礎(chǔ)，為有效及準確分析基因表達提供了保障，通過樣品間兩兩比較獲得了15組的比較數(shù)據(jù)，同時得到差異表達基因26861個，對差異基因生物信息學的分析為有效獲得成花相關(guān)基因建立重要基礎(chǔ)。

表4 花器官形成相關(guān)Unigenes

3.2 生物信息學分析

普通油茶的成花決定機理是進行花期調(diào)控和縮短花期研究的重要依據(jù)，對加速育種周期及解決花期不遇有重要意義。對花器官形成的研究同樣也非常重要，一方面通過了解相關(guān)因素及基因?qū)ㄆ鞴儆绊憴C制，可以為研究雄性不育及人工調(diào)控、花色及形態(tài)及親和性研究提供重要依據(jù)。AP1和AP2屬于MADS-box基因家族，已經(jīng)證明不僅和植物花器官形成有關(guān)，同時也與植物花芽分化關(guān)系密切[17-18]。從本文分析可以看出，AP1在花芽分化前期表達量較高，AP2在花芽分化期及花器官形成后期表達較突出，結(jié)合之前研究的FT基因表達結(jié)果，F(xiàn)T基因促進AP1和AP2基因表達，而不經(jīng)LFY基因。通過本文分析PI基因與雄蕊發(fā)育關(guān)系有關(guān)十分密切，這與花器官發(fā)育的ABC模型一致[19-21]，通過研究發(fā)現(xiàn)雄蕊形成期7月3日與花瓣形成期6月15日差異基因最多，上調(diào)基因達到19 408個，可見6月份到7月份是花器官形成的關(guān)鍵時期，通過對該時期基因調(diào)控可以有效控制花器官發(fā)育，如獲得單性花植株，觀賞性好的油茶。

參考文獻：

[1] 陳曉亞,湯章城.植物生理與分子生物學[M].北京:高等教育出版社,2007.

[2] Paul K B,Ruth M B,Joshua S M,etal.Multiple pathways in the decision to flower:enabling,promoting,and resetting[J].The Plant Cell Online,2004,16(suppl 1):18-31.

[3] 曾燕如,黎章矩,戴文圣.油茶開花習性的觀察研究[J].浙江農(nóng)林大學學報,2009,26(6):802-809.

[4] Martin C,Yuan T,Ghia E,etal.Green fluorescent protein as a marker for gene expression[J].Science,1994,263(5148):802-805.

[5] 寧順斌,王玲.幾種差別表達基因顯示技術(shù)及其在植物方面的應(yīng)用[J].生命科學,1999,11(3):140-143.

[6] Morrissy A S,Ryan D M,Allen D,etal.Next-generation tag sequencing for cancer gene expression profiling[J].Genome Research,2009,19(10):1825-1835.

[7] Elaine R M.Next-generation DNA sequencing methods[J].Annu.Rev.Genomics Hum.Genet.,2008(9):387-402.

[8] Rebecca C,Olivia A,Pamela A H.The gext generation:using new sequencing technologies to analyse gene regulation[J].Respirology,2011,16(2):210-222.

[9] 袁德義,鄒鋒,譚曉風,等.油茶花芽分化及雌雄配子體發(fā)育的研究[J].中南林業(yè)科技大學學報,2011,31(3):65-70.

[10] 王湘南,蔣麗娟,陳永忠,等.油茶花芽分化的形態(tài)解剖學特征觀測[J].中南林業(yè)科技大學學報,2011,31(8):22-27.

[11] 羅塔.油茶花芽分化過程形態(tài)及內(nèi)含物質(zhì)變化的研究[D].廣州:華南農(nóng)業(yè)大學,2012.

[12] Michiel J L,Seiya I,John N,etal.Open source clustering Software[J].Bioinformatics,2004,20(9):1453-1454.

[13] Michael B E,Paul T S,Patrick O B,etal.Cluster analysis and display of genome-wide expression patterns[J].Proceedings of the National Academy of Sciences,1998,95(25):14863-14868.

[14] 寧順斌,王玲,宋運淳.幾種差別表達基因顯示技術(shù)及其在植物方面的應(yīng)用[J].生命科學,1999,11(3):140-144.

[15] Kim S Y,Randal B,Jerome G,etal.Temporal aspects of DNA and RNA synthesis during human immunodeficiency virus infection:evidence for differential gene expression.[J].Journal of Virology,1989,63(9):3708-3713.

[16] 胡玉玲,姚小華,任華東,等.油茶花發(fā)育轉(zhuǎn)錄組測序及相關(guān)基因表達分析[J].林業(yè)科學,2014,50(9):36-43.

[17] Soraya P,Gary S D,Elvira B,etal.B and C Floral organ identity functions require sepallata mads-box genes[J].Nature,2000,405(6783):200-203.

[18] Detlef W,Elliot M M.The abcs of floral homeotic genes[J].Cell,1994,78(2):203-209.

[19] Heike W,Erica M,Marco T,etal.On reconciling the interactions between apetala2,Mir172 and agamous with the Abc model of flower development[J].Development,2010,137(21):3633-3642.

[20] Enrico C.Goethe and the Abc model of flower development[J].Comptes rendus De L’académie Des Sciences-series Iii-sciences De La Vie,2001,324(6):523-530.

[21] Elena R A,Eugenio A,Rafael B,etal.From Abc genes to regulatory networks,epigenetic landscapes and flower morphogenesis:making biological sense of theoretical approaches[C].Seminars in Cell & Developmental Biology,Elsevier,2010:108-117.