低溫脅迫下枇杷幼果轉(zhuǎn)錄組的De novo組裝和功能注釋

楊 偉，龔榮高，，石佳佳，賴 靜，廖明安，梁國魯

(1 四川農(nóng)業(yè)大學(xué) 園藝學(xué)院，四川 雅安625014；2 成都龍泉驛區(qū)農(nóng)村發(fā)展局，四川 成都610100；3 西南大學(xué) 園藝園林學(xué)院，重慶400700)

枇杷是調(diào)節(jié)市場淡季的重要亞熱帶水果，其以花和幼果過冬，開花坐果期正值一年中溫度最低的時(shí)期，因而枇杷的抗寒性已成為影響其產(chǎn)量和品質(zhì)的重要指標(biāo)之一。目前，有關(guān)枇杷抗寒性方面的研究對象相對單一，大多數(shù)集中在葉片凍害后的結(jié)構(gòu)、內(nèi)含物及光響應(yīng)等方面，而針對果實(shí)的研究相對較少[1]，尚未見有關(guān)枇杷幼果轉(zhuǎn)錄組測序等方面的研究報(bào)道。由于轉(zhuǎn)錄組表現(xiàn)的是生物或組織細(xì)胞在特定狀態(tài)下的所有RNA的總和，它包含了時(shí)間和空間的限定，是連接基因組遺傳信息與生物功能蛋白質(zhì)組的必然紐帶，因而，轉(zhuǎn)錄組的測序研究成為探究基因結(jié)構(gòu)及功能的基礎(chǔ)。運(yùn)用新一代的Illumina測序技術(shù)，可以得到基因組轉(zhuǎn)錄區(qū)域信息、鑒定轉(zhuǎn)錄發(fā)生位點(diǎn)、進(jìn)行可變剪切等，其精確的計(jì)數(shù)方法更可對基因進(jìn)行精確的定量分析，已廣泛應(yīng)用于生物學(xué)研究領(lǐng)域[2-3]。本研究通過對低溫脅迫下枇杷幼果轉(zhuǎn)錄組的分析，探索枇杷幼果低溫凍害的機(jī)理，并挖掘枇杷抗寒性相關(guān)基因，以期為枇杷引種、生態(tài)區(qū)劃和優(yōu)質(zhì)高產(chǎn)栽培提供相關(guān)的理論和實(shí)踐依據(jù)。

1 材料與方法

1.1 材 料

在四川省成都市龍泉驛區(qū)園藝場，選取10～15年生樹勢健壯、常規(guī)栽培管理、生長相對一致的5株“大五星”枇杷(EriobotryajaponicaLindley cv.Dawuxing)作為代表性植株,砧木為“解放鐘”(E.japonicaLindley cv.Jiefangzhong)，盛花后開始掛牌，并在易受凍害的花后75 d(75 days after full bloom，75 DAFB；正常發(fā)育果的直徑約1 cm)后，取樹冠中部外圍帶無病蟲害幼果的果枝20支，其中10支置于3 ℃的人工氣候箱內(nèi)低溫脅迫處理12 h，另外10支置于12 ℃的人工氣候箱內(nèi)作為對照。

1.2 方 法

1.2.1 RNA 的提取及質(zhì)量控制 枇杷幼果總RNA的提取參照Garg等[4]的方法。以幼果中心柱為對稱軸縱剖，在可食部分稱取0.3 g樣品置于液氮中，研磨成粉末狀。RNA樣本的質(zhì)量和濃度采用Garg等[4]的方法評判，用紫外分光光度計(jì)測定OD230、OD260和OD280，要求OD260/OD280和OD260/OD230的比值均≥1.8。運(yùn)用Agilent 2100 生物分析儀檢測RNA的完整性，RNA完整性值在8.6～10.0 的用于進(jìn)一步試驗(yàn)。即將樣品用RNase Free Water 在冰上溶解后，離心并充分混勻，取1 μL 樣品于70 ℃變性2 min后，利用Agilent 2100生物分析儀進(jìn)行檢測。

1.2.2 文庫制備與檢測 經(jīng)過DNase 消化、mRNA分離、mRNA 打斷，進(jìn)行反轉(zhuǎn)錄一鏈的合成，隨后進(jìn)行反轉(zhuǎn)錄二鏈的合成和末端修復(fù)。在末端修復(fù)產(chǎn)物cDNA 的3′末端加上A 堿基，讓接頭與A 堿基相連接。連接產(chǎn)物經(jīng)膠回收后進(jìn)行PCR 反應(yīng)及產(chǎn)物回收。使用Agilent 2100 生物分析儀和ABI實(shí)時(shí)熒光定量PCR系統(tǒng)對文庫質(zhì)量進(jìn)行檢測。

1.2.3 測序分析及Denovo組裝 由深圳華大基因研究中心完成轉(zhuǎn)錄組測序及Denovo組裝。應(yīng)用Illumina Hiseq 2000測序平臺，轉(zhuǎn)錄組文庫測序?yàn)镻E 90，reads為180 bp。Denovo組裝方法如下：使用短reads組裝軟件SOAPdenovo進(jìn)行轉(zhuǎn)錄組從頭組裝，得到Unigene。將Unigene序列與蛋白數(shù)據(jù)庫nr、Swiss-Prot、KEGG和COG做blastx比對(E-value<0.000 01)，確定Unigene的序列方向。如果不同庫之間的比對結(jié)果有矛盾，則按nr、Swiss-Prot、KEGG和COG的優(yōu)先級確定Unigene的序列方向；與以上4個(gè)庫皆比不上的Unigene，用軟件ESTScan預(yù)測其編碼區(qū)并確定序列的方向，給出其從5′到3′方向的序列；對于無法確定序列方向的Unigene，給出組裝軟件得到的序列。

1.2.4 基因功能注釋 由深圳華大基因研究中心完成基因測序及其功能注釋。功能注釋信息給出Unigene的蛋白功能注釋、COG功能注釋、Pathway注釋，以及GO(GO；http://www.geneontology.org)分類。

1.2.5 低溫脅迫下基因的表達(dá) 由深圳華大基因研究中心完成差異基因的表達(dá)分析。方法如下：Unigene表達(dá)量的計(jì)算采用FPKM法，其表達(dá)量用來比較低溫脅迫處理后的基因表達(dá)差異。隨后進(jìn)行差異Unigene的GO和Pathway分析。

2 結(jié)果與分析

2.1 枇杷幼果cDNA序列分析與組裝

運(yùn)用Illumina測序平臺對枇杷幼果cDNA樣品進(jìn)行測序，以弄清幼果的轉(zhuǎn)錄組及基因狀況。由表1可知，枇杷幼果中總體產(chǎn)生30 911 414條原始序列數(shù)據(jù)(reads)，經(jīng)過嚴(yán)格的質(zhì)量評估和數(shù)據(jù)過濾，得到26 024 270條clean reads。質(zhì)量要求Q(質(zhì)量值)>20的堿基百分比超過95%且G+C含量為35%～65%；數(shù)據(jù)過濾方法為去除測序接頭序列、重復(fù)冗余序列(含N比例大于10%)和質(zhì)量差的序列(Q≤10 的堿基數(shù)占整條read 的50%以上)。過濾后轉(zhuǎn)錄物Q>20的堿基百分比為97.47%，對其作進(jìn)一步分析。過濾后堿基G和C的數(shù)量占總堿基數(shù)的比例為20.3%～70.71%，平均為48.94%。

表 1 枇杷幼果轉(zhuǎn)錄組測序數(shù)據(jù)質(zhì)量統(tǒng)計(jì)結(jié)果

從表2可知，對reads進(jìn)行拼接，共組裝 116 723 條Contigs，Contig的平均長度為312 nt(核苷酸)，N50的長度則為487 bp(堿基數(shù))。組裝的Contigs長度為200～7 653 nt(圖1A)。

表 2 枇杷幼果基因序列組裝統(tǒng)計(jì)結(jié)果

圖 1 枇杷幼果基因序列組裝的Contigs和Unigenes的長度分布

由表2還可見，對Contigs進(jìn)一步拼接后，共組裝獲得轉(zhuǎn)錄物Unigenes序列64 814條，Unigenes的平均長度為541 nt，其N50的長度則是797 bp。Unigenes的長度為200～8 895 nt，總長度為35 036 073 nt(圖1B)。此外，長度大于1 000 nt的Contigs和Unigenes序列分別占到23.32% 和 13.79%。

2.2 枇杷幼果基因序列功能注釋信息

為了整體上了解轉(zhuǎn)錄物(Unigenes)的序列信息，將所有Unigene與nr、Swiss-Prot、KEGG、COG等蛋白或基因數(shù)據(jù)庫進(jìn)行blastx比對,以得到基因的方向、閱讀框及基因的功能注釋，并給出對應(yīng)的與數(shù)據(jù)庫中基因同源程度的衡量標(biāo)準(zhǔn)(E-value值)。利用nr、Swiss-Prot等蛋白數(shù)據(jù)庫對64 814條轉(zhuǎn)錄物數(shù)據(jù)進(jìn)行blastx比對的結(jié)果表明，有42 830條轉(zhuǎn)錄物序列被NCBI nr數(shù)據(jù)庫功能注釋。從圖2A的E-value分布圖可知，有42.5%的轉(zhuǎn)錄物序列具有極強(qiáng)的同源比對信息，而且有57.5%的轉(zhuǎn)錄物序列相似性E-value值在1.0E-5～1.0E-60。由圖2B的Identity分布圖可見，與具體蛋白匹配相似性超過80%的轉(zhuǎn)錄物占29.7%，而與具體蛋白匹配相似性范圍在19%～80%的轉(zhuǎn)錄物占70.3%，相似性低于40%的轉(zhuǎn)錄物占5.9%。由圖2C物種分布圖可見，有32.5%的轉(zhuǎn)錄物序列與葡萄序列高度匹配，與蓖麻匹配的占18.4%，與毛果楊匹配的占16.4%，與大豆匹配的占11.3%，與蘋果匹配的占4.1%，與蒺藜苜蓿匹配的占3.2%，與鼠耳芥匹配的占1.1%。

圖 2 枇杷幼果組裝轉(zhuǎn)錄物與nr數(shù)據(jù)庫比對后的分類分布圖

根據(jù)nr注釋、GO分析對枇杷幼果基因功能進(jìn)行可能的分類。由圖3可知，在GO功能分類體系中，有34 200條轉(zhuǎn)錄物具有具體功能定義，并獲得了273 645條轉(zhuǎn)錄物的功能注釋，其中有143 799條轉(zhuǎn)錄物注釋歸類為生物學(xué)過程，有90 020條歸為細(xì)胞組分，有39 826條歸為分子功能。在轉(zhuǎn)錄物的生物學(xué)過程功能類型中，細(xì)胞過程和代謝過程占絕大部分，其次是刺激反應(yīng)和生物調(diào)節(jié)過程。而在轉(zhuǎn)錄物的細(xì)胞組分功能類型中，細(xì)胞和細(xì)胞分離所占比例最高，其次是細(xì)胞器和生物膜。蛋白結(jié)合和催化活性是轉(zhuǎn)錄物分子功能類型中的主體部分。

為了進(jìn)一步預(yù)測和分類轉(zhuǎn)錄物，運(yùn)用COG庫比對注釋了枇杷幼果的14 538條轉(zhuǎn)錄物，證明其有蛋白功能。COG功能注釋了至少25個(gè)功能類別的具體蛋白功能，它包括了細(xì)胞結(jié)構(gòu)、生物化學(xué)代謝、分子過程、信號傳遞等生理生化過程(圖4)。其中，在COG轉(zhuǎn)錄物功能注釋分類中，一般功能基因所占比例最大(4 622，15%)，所占比例5%以上的有轉(zhuǎn)錄(3 181，10.32%)，復(fù)制、重組與修復(fù)(2 422，7.86%)，蛋白質(zhì)翻譯后修飾與轉(zhuǎn)運(yùn)、分子伴侶(2 360，7.66%)，信號傳導(dǎo)機(jī)制(2 060，6.69%)，翻譯、核糖體結(jié)構(gòu)與生物合成(1 969，6.42%)，碳水化合物運(yùn)輸與代謝(1 889，6.13%)和功能未知(1 855，6.02%)，而防御機(jī)制(321，0.96%)、染色體結(jié)構(gòu)與變化(314，0.94%)、核苷酸運(yùn)輸與代謝(274，0.82%)、RNA加工與修飾(189，0.57%)、胞外結(jié)構(gòu)(16，0.05%)和核酸結(jié)構(gòu)(3，0.01%)所占比例則均在1%以下。

圖 3 根據(jù)GO分析功能注釋枇杷幼果轉(zhuǎn)錄物

2.3 低溫脅迫下枇杷幼果的差異表達(dá)基因

為探索低溫脅迫對枇杷幼果轉(zhuǎn)錄物轉(zhuǎn)錄水平的影響，對低溫脅迫處理的2個(gè)樣本的基因差異表達(dá)，按照Audic等[5]根據(jù)轉(zhuǎn)錄物測序的基因表達(dá)差異方法進(jìn)行篩選，結(jié)果總共檢測到57 792條轉(zhuǎn)錄物與對照相比在表達(dá)上發(fā)生明顯變化，這些轉(zhuǎn)錄物定位到4 017個(gè)基因中，有1 619個(gè)基因表達(dá)上調(diào)，有2 398個(gè)基因表達(dá)下調(diào)。

2.4 低溫脅迫下枇杷幼果差異表達(dá)基因的功能注釋

通過對枇杷幼果低溫脅迫處理后的轉(zhuǎn)錄物的分析和定位，找到差異表達(dá)的基因。隨后通過對這些差異基因進(jìn)行GO分類和KEGG通徑分析，弄清低溫脅迫差異表達(dá)基因的功能，特別是在細(xì)胞的生物代謝途徑方面的基因功能?？傮w上看，枇杷幼果低溫脅迫誘導(dǎo)產(chǎn)生了2 007個(gè)差異表達(dá)的轉(zhuǎn)錄物。這些差異表達(dá)轉(zhuǎn)錄物可以被歸于113條生化代謝或信號轉(zhuǎn)導(dǎo)途徑中，但最主要的途徑是參與生化代謝和信號轉(zhuǎn)導(dǎo)。這113條途徑絕大部分是代謝途徑、生物合成次生代謝物、細(xì)胞的內(nèi)吞作用、甘油代謝、植物激素信號傳導(dǎo)、醚酯代謝、淀粉和糖代謝、植物病原菌相互作用等途徑。本試驗(yàn)結(jié)果中，值得注意的是枇杷幼果低溫脅迫產(chǎn)生油菜素內(nèi)酯生物合成和磷脂酰肌醇信號系統(tǒng)的幾個(gè)基因差異表達(dá)上調(diào)，其中，共有6個(gè)參與調(diào)控油菜素內(nèi)酯生物合成的基因差異表達(dá)上調(diào)(圖5、表3)。而PLC(磷脂酶 C)基因則在磷脂酰肌醇信號系統(tǒng)中表達(dá)上調(diào)(圖6、表3)。

圖 4 根據(jù)COG分析功能注釋枇杷幼果轉(zhuǎn)錄物

3 討 論

轉(zhuǎn)錄組是生物或細(xì)胞特定狀態(tài)下轉(zhuǎn)錄的所有RNA，含mRNA及非編碼RNA。轉(zhuǎn)錄組的測序逐漸成為研究基因結(jié)構(gòu)及功能的基礎(chǔ)。根據(jù)轉(zhuǎn)錄測序信息可以弄清許多生物過程，比如基因表達(dá)情況[6]、試驗(yàn)處理或病害后的基因表達(dá)譜、基因調(diào)控、組織中生物標(biāo)識物、基因挖掘、基因含量、分子系統(tǒng)學(xué)中的保守同源基因分離等[7-9]。因此，通過轉(zhuǎn)錄組的測序技術(shù)可以為研究對象提供許多候選功能基因。高通量技術(shù)是目前轉(zhuǎn)錄組測序技術(shù)中一種快捷而又可靠的方法，許多研究通過Denovo組裝進(jìn)行無參考基因組的非模式生物轉(zhuǎn)錄組分析[4]。然而，尚無有關(guān)枇杷果實(shí)轉(zhuǎn)錄組信息方面的研究報(bào)道。本研究中，運(yùn)用Illumina高能量技術(shù)，對無參考基因組的枇杷果實(shí)通過Denovo測序儀進(jìn)行了組裝轉(zhuǎn)錄分析，鑒定和分析了大量表達(dá)的基因，得到了2千多萬條序列數(shù)據(jù)，用多種互補(bǔ)方法組裝了42 830個(gè)轉(zhuǎn)錄物。功能注釋為枇杷研究提供了有價(jià)值的特定細(xì)胞過程、分子功能和生物代謝途徑。本試驗(yàn)結(jié)果表明，高能量測序技術(shù)是研究無參考基因組生物的有效平臺。這些枇杷果實(shí)轉(zhuǎn)錄組數(shù)據(jù)為今后開發(fā)枇杷分子和細(xì)胞標(biāo)記提供了理想的序列資源。

圖 5 通過KEGG 推定低溫脅迫下枇杷幼果油菜素內(nèi)酯的生物合成途徑

圖 6 通過KEGG 推定低溫脅迫下枇杷幼果磷脂酰肌醇的信號系統(tǒng)途徑

表 3 低溫脅迫下枇杷幼果油菜素內(nèi)脂生物合成和磷脂酰肌醇信號系統(tǒng)途徑中的差異基因

低溫脅迫影響植物的生長發(fā)育，但植物在低溫脅迫下會(huì)產(chǎn)生一個(gè)調(diào)節(jié)系統(tǒng)，有千百個(gè)基因參與到脅迫下的信號傳遞和下游的脅迫反應(yīng)過程中[10-12]。本研究通過枇杷低溫脅迫轉(zhuǎn)錄組分析，試圖鑒定可能參與到低溫脅迫的一組基因，經(jīng)鑒定發(fā)現(xiàn)有4 017個(gè)基因有可能參與到枇杷低溫脅迫中，其中有1 619個(gè)基因表達(dá)上調(diào)，而2 398個(gè)基因表達(dá)下調(diào)。這些基因促進(jìn)了人們對枇杷幼果耐低溫脅迫分子機(jī)制方面的認(rèn)識。

通過分子功能分析可全面探討枇杷幼果對低溫脅迫的反應(yīng)。本研究結(jié)果表明，有6個(gè)差異基因表達(dá)上調(diào)，其參與調(diào)控油菜素內(nèi)酯的生物合成。油菜素內(nèi)酯(BRs)是植物甾醇類激素，在植物的正常生長發(fā)育中是不可缺少的。目前許多植物生理學(xué)家已將其列為植物的第6大類激素。而且，油菜素內(nèi)酯還能增加植物的環(huán)境抗逆性。有研究發(fā)現(xiàn)，油菜素內(nèi)酯能提高擬南芥幼苗忍耐低溫的能力[13]。本試驗(yàn)結(jié)果表明，油菜素內(nèi)酯參與了枇杷幼果低溫脅迫調(diào)節(jié)系統(tǒng)。

本試驗(yàn)也發(fā)現(xiàn)，PLC(磷脂酶 C)的基因差異表達(dá)上調(diào)參與調(diào)控了磷脂酰肌醇信號系統(tǒng)，它是磷脂酰肌醇信號轉(zhuǎn)換系統(tǒng)非常重要的組成部分。磷脂酰肌醇信號轉(zhuǎn)換產(chǎn)生磷脂酶 C和Ca2+，而它們使胞外信號轉(zhuǎn)換為胞內(nèi)信號，并產(chǎn)生各種非生物脅迫響應(yīng)。肌醇磷脂被磷脂酶 C水解，在受體刺激的同時(shí)，打開鈣庫釋放Ca2+到細(xì)胞質(zhì)中與鈣調(diào)蛋白結(jié)合，隨后參與一系列的反應(yīng)。研究發(fā)現(xiàn)，AtPLC1(類磷脂酶 C突變體)很可能參與環(huán)境脅迫信號轉(zhuǎn)換途徑，提高AtPLC1水平可增強(qiáng)信號，這樣有助于提高植物對環(huán)境脅迫的適應(yīng)能力[14]。當(dāng)PI-PLC的協(xié)同作用受抑制后，RD29A和COR47(冷應(yīng)答信號通路的基因)2個(gè)基因表達(dá)下調(diào)，滲透脅迫誘導(dǎo)基因也出現(xiàn)下調(diào)[15]。因而，由本試驗(yàn)結(jié)果可以推測磷脂酶 C很可能參與調(diào)節(jié)枇杷幼果低溫抗性反應(yīng)。

本研究首次報(bào)道了運(yùn)用高能量測序技術(shù)并通過Denovo組裝，進(jìn)行枇杷幼果轉(zhuǎn)錄組分析和低溫脅迫下基因差異表達(dá)分析。本研究是在沒有參考基因的基礎(chǔ)上完成轉(zhuǎn)錄本的組裝，得到了高質(zhì)量的枇杷轉(zhuǎn)錄組序列，并進(jìn)一步通過對枇杷轉(zhuǎn)錄組的注釋，深入了解了枇杷幼果基因及其功能。這些發(fā)現(xiàn)將更有利于探索枇杷基因組資源，也能促進(jìn)今后枇杷基因組和育種方面的研究。此外，本研究運(yùn)用高能量測序技術(shù)分析低溫脅迫下枇杷幼果基因的表達(dá)情況，表明低溫脅迫下功能基因的研究方法是可行的，這也必然會(huì)促進(jìn)對低溫脅迫下枇杷幼果分子機(jī)理的深入了解。

[參考文獻(xiàn)]

[1] 鄭國華,張賀英,鐘秀容.低溫脅迫下枇杷葉片細(xì)胞超微結(jié)構(gòu)及膜透性和保護(hù)酶活性的變化 [J].中國生態(tài)農(nóng)業(yè)學(xué)報(bào)，2009,17(4):739-745.

Zheng G H,Zhang H Y,Zhong X R.Changes in cell ultra-structure,membrane permeability and protective enzyme activity inEriobotryajaponicaLindl.leaves under cold stress [J].Chinese Journal of Eco-Agricultre,2009,17(4):739-745.(in Chinese)

[2] Krost C,Petersen R,Schmidt E R.The transcriptomes of columnar and standard type apple trees (Malus×domestica):A comparative study [J].Gene,2012,498:223-230.

[3] Varshney R K,Hoisington D A,Tyagi A K.Advances in cereal genomics and applications in crop breeding [J].Trends Biotechnol,2006,24:490-499.

[4] Garg R,Patel R K,Tyagi A K,et al.Denevoassembly of chic-kpea transcriptome using short reads for gene discovery and marker identification [J].DNA Res,2011,18:53-63.

[5] Audic S,Claverie J M.The significance of digital gene expression profiles [J].Genome Res,1997,7:986-995.

[6] Torres T T,Metta M,Ottenwalder B,et al.Gene expression pr-ofiling by massively parallel sequencing [J].Genome Res,2008,18:172-177.

[7] Hegedus Z,Zakrzewska A,Agoston V C,et al.Deep sequencing of the zebrafish transcriptome response to mycobacterium infection [J].Mol Immunol,2009,46:2918-2930.

[8] Hughes J,Longhorn S J,Papadopoulou A,et al.Dense taxono-mic EST sampling and its applications for molecular systematics of theColeoptera(beetles) [J].Mol Biol Evol,2006,23:268-278.

[9] He R.Next-generation sequencing-based transcriptomic and pr-oteomic analysis of the common reed,Phrafmites australis (Poaceae),reveals genes involved in invasiveness and rhizome specificity [J].Amer J Bot,2012,99:232-247.

[10] Kreps J A,Wu Y,Chang H S,et al.Transcriptome changes forArabidopsisin response to salt,osmotic,and cold stress [J].Plant Physiol,2002,130:2129-2141.

[11] Shinozaki K,Yamaguchi-Shinozaki K,Seki M.Regulatory network of gene expression in the drought and cold stress responses [J].Cur Opin Plant Biol,2003,6:410-417.

[12] Xiong L,Schumaker K S,Zhu J K.Cell signaling during cold,drought,and salt stress [J].Plant Cell,2002,14:165-183.

[13] Divi U K,Krishna P.Overexpression of the brassinosteroid biosynthetic gene AtDWF4 inArabidopsisseeds overcomes abscisic acid-induced inhibition of germination and increases cold tolerance in transgenic seedlings [J].Plant Growth Reg,2010,29:385-393.

[14] Hirayama T,Ohto C,Mizoguchi T,et al.A gene encoding a phosphatidylinositol-specific phospholipase C is induced by dehydration and salt stress inArabidopsisthaliana[J].PNAS,1995,92:3903-3907.

[15] Takahashi S.Hyperosmotic stress induces a rapid and transient increase in inositol 1,4,5-trisphosphate independent of abscisic acid inArabidopsiscell culture [J].Plant Cell Physio,2001,42:214-222.