基于轉錄組測序的枇杷腋芽成花相關基因的篩選

章加應,肖海濤,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信號途徑相關差異基因時空表達模式與春花腋芽開花進程呈現相一致的趨勢;同時,差異基因在枇杷品系春花中的表達量顯著高于‘火炬’,表明差異表達基因在枇杷品系春花腋芽成花過程中發揮至關重要的作用,調控春花腋芽的成花進程。本研究有助于更好地了解枇杷腋芽成花過程中的基因轉錄水平,并可為闡明枇杷腋芽成花的分子遺傳機制奠定理論基礎。