王 丹，龔榮高

(四川農(nóng)業(yè)大學(xué) 園藝學(xué)院，四川 成都 611130)

枇杷(EriobotryajaponicaLindl.)是亞熱帶果樹中比較重要的一種，喜溫暖氣候和肥水濕潤(rùn)的土壤，在我國(guó)有著豐富的種植資源和優(yōu)良品種[1]。枇杷對(duì)干旱有一定的耐受性，但在生長(zhǎng)過程中易發(fā)生花期、幼果期和果實(shí)成熟期旱害，導(dǎo)致枇杷的坐果率和產(chǎn)量嚴(yán)重降低[2-3]。而干旱是不可避免的，它是影響植物生長(zhǎng)發(fā)育最重要的非生物因素之一[4]。關(guān)于枇杷干旱的研究在生理生化方面較多[5]，羅華建等[2]用盆栽的方法證明：干旱抑制枇杷的生長(zhǎng)和光合作用；也有研究表明，輕度的干旱脅迫有助于誘導(dǎo)枇杷花芽分化，使果實(shí)成熟期提前，對(duì)枇杷的生長(zhǎng)有促進(jìn)作用[6-7]。但由于枇杷基因組背景信息缺乏，尚未見有關(guān)枇杷葉片轉(zhuǎn)錄組測(cè)序等方面的研究報(bào)道，其抗逆機(jī)理的深入研究及抗旱分子標(biāo)記開發(fā)等發(fā)展相對(duì)滯后。

轉(zhuǎn)錄組測(cè)序(RNA sequencing)利用第2代高通量測(cè)序技術(shù)進(jìn)行cDNA測(cè)序，能全面快速地獲得某一狀態(tài)下的全部轉(zhuǎn)錄本信息，可以用于新轉(zhuǎn)錄本預(yù)測(cè)[8]。對(duì)于缺乏基因組信息的物種而言，采用轉(zhuǎn)錄組測(cè)序技術(shù)可獲得大量的轉(zhuǎn)錄本信息，從中發(fā)掘重要的功能基因，是揭示植物優(yōu)良性狀的重要研究手段[9]。轉(zhuǎn)錄組技術(shù)還可以研究不同器官、不同環(huán)境脅迫下基因表達(dá)的差異，過去10年里，上千個(gè)基因和數(shù)十條代謝和信號(hào)通路被證明與植物的耐旱性有關(guān)，擬南芥[10]、水稻[11]和短柄草[12]中的許多基因在轉(zhuǎn)錄水平都被證實(shí)參與干旱的調(diào)控。而參與植物光合作用、ABA合成、信號(hào)轉(zhuǎn)導(dǎo)、滲透物質(zhì)的生物合成和轉(zhuǎn)錄因子(WRKY、bZIP家族等)調(diào)控等的基因在干旱脅迫下表達(dá)異常，表明植物對(duì)水分刺激應(yīng)答的過程受眾多分子和細(xì)胞通路的調(diào)控[13]。但在其他植物中，類似的研究非常有限，植物對(duì)干旱耐性的分子機(jī)制在很大程度上還是未知的。

本研究利用Illumina HiSeq高通量測(cè)序技術(shù)對(duì)干旱前后的枇杷葉片轉(zhuǎn)錄組進(jìn)行研究，結(jié)合生物信息學(xué)方法對(duì)所獲得的差異表達(dá)基因(DEGs，Differently expressed genes)進(jìn)行功能的注釋、分類和代謝途徑分析，在轉(zhuǎn)錄水平上研究干旱脅迫下枇杷關(guān)鍵基因的表達(dá)，挖掘與干旱相關(guān)的基因，為枇杷抗旱相關(guān)基因的克隆以及功能分析和抗旱育種等研究提供理論基礎(chǔ)。

1 材料和方法

1.1 試驗(yàn)材料

試驗(yàn)材料為四川龍泉驛區(qū)長(zhǎng)勢(shì)良好的3年生大五星枇杷嫁接苗，定植于口徑40 cm、高30 cm的花盆中，每盆種植1株。所用土壤為菜園土和沙壤土適量混合，并加入10%腐熟的有機(jī)肥和0.1%的復(fù)合肥，后期處理過程中不再施肥。待馴化后選取高度、長(zhǎng)勢(shì)一致的材料進(jìn)行干旱脅迫控制性試驗(yàn)。試驗(yàn)于2015年6-9月在四川農(nóng)業(yè)大學(xué)農(nóng)場(chǎng)進(jìn)行。試驗(yàn)開始時(shí)將材料移至大棚，每盆充分灌溉，使土壤含水量一致，停止灌溉后采用自然耗水進(jìn)行干旱脅迫處理。

1.2 試驗(yàn)方法

1.2.1 試驗(yàn)處理 以土壤含水量(占田間持水量的百分?jǐn)?shù))設(shè)置正常供水(對(duì)照，田間持水量的61%～65%)、輕度干旱脅迫(田間持水量的50%～54%)、中度干旱脅迫(田間持水量的39%～43%)和重度干旱脅迫(田間持水量的28%～32%)4個(gè)水分處理。將輕度、中度和重度3個(gè)脅迫的葉片取樣后混合和對(duì)照葉片液氮速凍后于-70 ℃貯存，用于RNA的提取。

1.2.2 總RNA的提取與測(cè)序 采用TRIzol法提取總RNA，用Agilent 2100 Technologies檢測(cè)樣品RNA的完整性和濃度，用帶有Oligo(dT)的磁珠富集mRNA，再加入緩沖液將其隨機(jī)打斷成片段，以mRNA片段為模板用隨機(jī)引物合成cDNA第一鏈，再加入緩沖液、dNTPs、RNase H和DNA Polymerase Ⅰ合成cDNA第二鏈，將洗脫純化后的雙鏈cDNA進(jìn)行末端修復(fù)、加堿基A、加測(cè)序接頭處理，然后經(jīng)瓊脂糖凝膠電泳回收目的大小片段最后再進(jìn)行PCR擴(kuò)增，完成整個(gè)文庫(kù)制備后用Illumina HiSeq進(jìn)行測(cè)序，測(cè)序策略為PE100。

1.2.3 測(cè)序數(shù)據(jù)的組裝和功能注釋 采用Trinity(20140717版)[14]將除去接頭序列、兩端低質(zhì)量序列和低度復(fù)雜序列的高質(zhì)量序列數(shù)據(jù)(Clean reads)組裝出全長(zhǎng)轉(zhuǎn)錄本。對(duì)組裝出來的轉(zhuǎn)錄本利用TransDecoder(20140717版)鑒定編碼區(qū)域，預(yù)測(cè)開放閱讀框ORF(Opening reading fragment)。針對(duì)預(yù)測(cè)出來的ORF和Contig序列用Trinotate(20131110版)對(duì)其進(jìn)行功能注釋，包括與同源性搜索(NCBI-Blast)、蛋白質(zhì)結(jié)構(gòu)域鑒定(HMMER/PFAM)、蛋白質(zhì)信號(hào)預(yù)測(cè)(SingalP/TmHMM)、Uniprot(Universal Protein：Swiss-Prot、TrEMBL和PIR-PSD)、Blast2GO軟件進(jìn)行的GO(Gene Ontology)、COG(Cluster of Orthologous Groups)和KEGG(Kyoto Encyclopedia of Genes and Genomes)進(jìn)行序列比對(duì)，獲得基因的功能注釋信息。

1.2.4 差異基因篩選及富集分析 基因表達(dá)量的計(jì)算使用標(biāo)準(zhǔn)化轉(zhuǎn)錄本數(shù)(RPKM)法，選取|log2Ratio|≥1和q<0.05的基因作為DEGs。應(yīng)用超幾何檢驗(yàn)，以q<0.05 為閾值，篩選出在DEGs中顯著富集的條目和通路，分析DEGs行使的主要生物學(xué)功能和代謝通路。

1.2.5 實(shí)時(shí)熒光定量qRT-PCR 利用TRIzol試劑提取RNA后，采用Aidlab公司反轉(zhuǎn)錄試劑盒(TUREscript 1st Stand cDNA SYNTHESIS Kit)進(jìn)行反轉(zhuǎn)錄操作，選取測(cè)序結(jié)果中表達(dá)數(shù)量較多且在各條目和通路中顯著富集的8個(gè)差異基因，依據(jù)基因的核酸序列設(shè)計(jì)熒光定量引物，采用20 μL反應(yīng)體系，對(duì)熒光定量PCR的引物濃度、退火溫度進(jìn)行優(yōu)化，得到最佳的熒光定量PCR，每個(gè)樣品3個(gè)重復(fù)，引物見表1。采用2-ΔΔCT法計(jì)算各基因相對(duì)表達(dá)量。

表1 熒光定量PCR引物

2 結(jié)果與分析

2.1 測(cè)序結(jié)果和注釋數(shù)據(jù)

枇杷葉片通過新一代高通量測(cè)序技術(shù)進(jìn)行轉(zhuǎn)錄組測(cè)序，共獲得114 021 518個(gè)原始reads，過濾掉接頭污染、低質(zhì)量、含N比例大于50%的序列后，得到90 517 644個(gè)干凈reads。堿基Q30為92.5%，表明測(cè)序的數(shù)據(jù)量和質(zhì)量都比較高。利用Trinity將干凈reads進(jìn)行從頭組裝形成轉(zhuǎn)錄本，共發(fā)現(xiàn)88 530條轉(zhuǎn)錄本，平均長(zhǎng)度740.64 bp，N50長(zhǎng)度為1 189 bp，其中最長(zhǎng)的轉(zhuǎn)錄本長(zhǎng)度為10 189 bp，最短的長(zhǎng)度為201 bp。結(jié)果表明，長(zhǎng)度在200～400 bp的轉(zhuǎn)錄本數(shù)目最多，為41 141條，占所有轉(zhuǎn)錄本的46.47%(表2)。

組裝出的轉(zhuǎn)錄本共預(yù)測(cè)到41 748條ORF，占所有轉(zhuǎn)錄本總數(shù)的47.16%，ORF平均長(zhǎng)度為891.54 bp，N50長(zhǎng)度為1 137 bp，長(zhǎng)度400～600 bp所占的比例最高。將Contigs和ORF序列注釋到Uniprot、GO和KEGG數(shù)據(jù)庫(kù)中，對(duì)注釋到每個(gè)庫(kù)中的數(shù)目進(jìn)行統(tǒng)計(jì)，BlastX數(shù)據(jù)庫(kù)注釋的最多，為22 039條，占所有轉(zhuǎn)錄本的24.89%；BlastP數(shù)據(jù)庫(kù)注釋的數(shù)量為16 352條；KEGG數(shù)據(jù)庫(kù)注釋的數(shù)量為5 484條，占6.19%；注釋到GO數(shù)據(jù)庫(kù)的有20 912條。

表2 枇杷葉片轉(zhuǎn)錄組Transcript和ORF長(zhǎng)度統(tǒng)計(jì)

2.2 COG功能分類

COG分類結(jié)果表明，13 285條轉(zhuǎn)錄本具有功能注釋信息且被注釋到24個(gè)功能類別。涵蓋基因最多的是一般功能預(yù)測(cè)(18.14%)，其次是翻譯后修飾、蛋白質(zhì)折疊和分子伴侶(10.65%)，翻譯、核糖體結(jié)構(gòu)和生物合成(10.64%)，復(fù)制、重組和修復(fù)(9.73%)，碳水化合物運(yùn)輸與代謝(7.19%)，信號(hào)轉(zhuǎn)導(dǎo)機(jī)制(6.86%)和能量產(chǎn)生與轉(zhuǎn)化(6.57%)。所有類別中，細(xì)胞運(yùn)動(dòng)(0.11%)和核酸結(jié)構(gòu)(0.01%)的數(shù)量最少(圖1)。

A.加工與修飾RNA ；B.染色體結(jié)構(gòu)與變化；C.能量產(chǎn)生與轉(zhuǎn)化；D.細(xì)胞周期調(diào)控與分裂、染色體重排；E.氨基酸運(yùn)輸與代謝；F.核苷酸運(yùn)輸與代謝；G.碳水化合物運(yùn)輸與代謝；H.輔酶運(yùn)輸與代謝；I.脂質(zhì)運(yùn)輸與代謝；J.翻譯、核糖體結(jié)構(gòu)與生物合成；K.轉(zhuǎn)錄；L.復(fù)制、重組與修復(fù)；M.細(xì)胞壁、膜生物合成Cell；N.細(xì)胞運(yùn)動(dòng)；O.翻譯后修飾、蛋白質(zhì)折疊和分子伴侶；P.無機(jī)離子運(yùn)輸與代謝；Q.次生代謝物合成、運(yùn)輸和代謝；R.一般功能預(yù)測(cè)；S.未知功能；T.信號(hào)轉(zhuǎn)導(dǎo)機(jī)制；U.細(xì)胞內(nèi)轉(zhuǎn)運(yùn)、分泌和囊泡運(yùn)輸；V.防御機(jī)制；Y.核酸結(jié)構(gòu)；Z.細(xì)胞骨架。

2.3 差異基因的GO分類和富集分析

通過計(jì)算基因的RPKM值，將兩樣品間RPKM的log2倍數(shù)絕對(duì)值大于1和q<0.05的基因作為差異基因，共篩選出25 197個(gè)差異表達(dá)的基因，其中上調(diào)基因13 482個(gè)，下調(diào)基因11 715個(gè)。54個(gè)GO分類注釋結(jié)果(圖2)表明：生物學(xué)過程分類中，注釋數(shù)量最多的是細(xì)胞進(jìn)程(上調(diào)基因?yàn)?.23%，下調(diào)基因?yàn)?.84%)和代謝途徑(上調(diào)基因7.40%，下調(diào)基因8.0%)。在細(xì)胞組分中，細(xì)胞和細(xì)胞部分擁有的基因數(shù)最多(上調(diào)基因均為10.08%，下調(diào)基因均為11.94%)。分子功能中結(jié)合作用(上調(diào)基因分別為7.55%，下調(diào)基因?yàn)?.51%)的數(shù)量最多。

通過顯著富集計(jì)算，生物學(xué)過程中有345條GO條目顯著富集，細(xì)胞組分有87條，分子功能中有173條。其中分子功能中富集差異基因最多的前3條GO條目分別是催化活性、轉(zhuǎn)移酶活性和水解酶活性(表3)。

表3 分子功能中DEGs富集最多的前10 類GO注釋

2.4 代謝通路分析

KEGG分析結(jié)果表明，1 618個(gè)DEGs在291個(gè)標(biāo)準(zhǔn)代謝通路中得到注釋。對(duì)KEGG中每個(gè)代謝通路用超幾何檢驗(yàn)進(jìn)行富集分析，找出DEGs顯著性富集的30條代謝通路(表4)，主要有：代謝途徑、次生代謝物的二次代謝、不同生境微生物代謝作用和核糖體。這些途徑為研究枇杷干旱脅迫的具體過程、功能和途徑提供了珍貴的資源。其中，與激素代謝相關(guān)的3條代謝途徑：雙萜類、油菜素類固醇和類胡蘿卜素生物合成途徑中的差異基因呈現(xiàn)出較為一致的表達(dá)。雙萜類生物合成途徑的13個(gè)基因中有6個(gè)貝殼杉烯酸氧化酶和4個(gè)貝殼杉烯氧化酶基因均下調(diào)表達(dá)。類胡蘿卜素生物合成途徑共14個(gè)DEGs，其中有2個(gè)9-順式環(huán)氧類胡蘿卜素雙加氧酶基因和10個(gè)上調(diào)表達(dá)的脫落酸(ABA)-8′羥化酶基因。油菜素內(nèi)固醇生物合成中DEGs共14個(gè)，12個(gè)均表達(dá)下調(diào)，其中有9個(gè)基因?yàn)镃YP734A1。

圖2 干旱脅迫下枇杷葉片DEGs的GO分類

表4 DEGs富集的30條KEGG通路

2.5 差異基因及其實(shí)時(shí)熒光定量驗(yàn)證

通過功能分析，除了部分未知功能的基因外，可能與干旱脅迫響應(yīng)相關(guān)的基因主要有信號(hào)轉(zhuǎn)導(dǎo)相關(guān)基因：促分裂原活化蛋白激酶(MAPK)、富含亮氨酸重復(fù)序列類受體蛋白激酶(LRR-RLKs)、G-type-S受體蛋白、ERF、MYB、WRKY轉(zhuǎn)錄因子和熱激蛋白等；抗氧化酶基因：過氧化氫酶、谷胱甘肽過氧化物酶、過氧化物酶；滲透調(diào)節(jié)物質(zhì)合成相關(guān)基因：海藻糖-6-磷酸合酶、乙醛脫氫酶、谷氨酸脫氫酶；激素合成相關(guān)基因：細(xì)胞色素P450家族、貝殼杉烯酸氧化酶、貝殼杉烯氧化酶、黃烷酮醇-4-還原酶、三角狀五肽重復(fù)(PPR)蛋白等。

選擇了差異表達(dá)數(shù)量豐富且在GO分類和KEGG代謝通路中都顯著富集的8個(gè)差異基因(表1)進(jìn)行實(shí)時(shí)熒光定量分析，驗(yàn)證這8個(gè)基因在葉片中的表達(dá)情況(圖3)，其中的6個(gè)基因c1047(過氧化物酶-上調(diào))、c65100(蛋白激酶byr2-上調(diào))、c56685(絲氨酸蘇氨酸蛋白激酶-上調(diào))、c68887(脫落酸8′-羥化酶-上調(diào))、c22815(細(xì)胞色素P450 734A1-下調(diào))和c68870(吲哚-3-乙酰乙酸合酶-上調(diào))的表達(dá)結(jié)果與轉(zhuǎn)錄組測(cè)序一致，這也證明了測(cè)序方法的正確性。

圖3 八個(gè)差異表達(dá)基因的定量PCR表達(dá)

3 討論

抗旱性一個(gè)復(fù)雜的多基因特征[15]，高通量測(cè)序能夠?yàn)橹参锘虮磉_(dá)的全面分析提供合理的數(shù)據(jù)資源，尤其適合沒有參考基因組信息的物種[16]。本研究采用高通量測(cè)序共獲得枇杷葉片轉(zhuǎn)錄本88 530條，預(yù)測(cè)出41 748條ORF序列，將測(cè)序結(jié)果與Uniprot、GO和KEGG等數(shù)據(jù)庫(kù)進(jìn)行比對(duì)。通過篩選，得到25 197個(gè)DEGs(P<0.01)，然而有大量的DEGs未得到注釋且未被報(bào)道過，這很大程度上豐富了枇杷的基因信息資源庫(kù)[17-18]。注釋結(jié)果表明，大部分差異表達(dá)的基因參與了如催化酶、轉(zhuǎn)移酶和水解酶活性的調(diào)節(jié)，大分子物質(zhì)的代謝，滲透調(diào)節(jié)物質(zhì)合成，信號(hào)轉(zhuǎn)導(dǎo)和激素水平的調(diào)節(jié)。

受體蛋白激酶(RLKs)是一組膜蛋白，盡管其功能和詳細(xì)的分子機(jī)制還不清楚，但RLK基因轉(zhuǎn)錄變化表明他們可能在干旱中起到重要的作用[19]。有研究表明，蛋白激酶信號(hào)基因在干旱時(shí)被誘導(dǎo)[20]。水稻受體蛋白激酶4在脅迫中被發(fā)現(xiàn)上調(diào)表達(dá)[11]。本研究中的絲氨酸/蘇氨酸蛋白激酶、LRR-RLKs等大部分基因均不同程度的下調(diào)，MAPK則表現(xiàn)為上調(diào)，說明這幾類蛋白激酶對(duì)干旱脅迫應(yīng)答調(diào)控起不同的作用。蔡國(guó)華[21]的研究結(jié)果表明，在干旱條件下，過表達(dá)玉米MAPKK在擬南芥中通過ROS和依賴ABA途徑增強(qiáng)了對(duì)干旱脅迫的抗性，另有研究表明，一些RLKs，如擬南芥RPK1在干旱脅迫時(shí)與ABA和BR(油菜素內(nèi)固醇)信號(hào)通路有聯(lián)系，而不同受體蛋白類型在信號(hào)通路中的作用還有待進(jìn)一步探討[22]。

參與脅迫信號(hào)傳遞網(wǎng)絡(luò)的轉(zhuǎn)錄因子，尤其是ERF、WRKY、MYB等家族是信號(hào)轉(zhuǎn)導(dǎo)通路上的關(guān)鍵調(diào)控因子，對(duì)干旱脅迫響應(yīng)和水分刺激應(yīng)答發(fā)揮了重要調(diào)節(jié)作用[23]。ERF 類轉(zhuǎn)錄因子是AP2/EREBP 轉(zhuǎn)錄因子家族的一個(gè)亞家族，通過參與乙烯、ABA、茉莉酸等多種信號(hào)轉(zhuǎn)導(dǎo)途徑，在非生物脅迫應(yīng)答調(diào)控中發(fā)揮重要的作用，一些ERFs在轉(zhuǎn)基因植物中能加強(qiáng)對(duì)干旱的耐受性，過表達(dá)的ERF類轉(zhuǎn)錄因子WXP1、WXP2可以增加葉表皮蠟質(zhì)的積累，顯著減少葉片脫水、提高植株對(duì)干旱脅迫的耐受性[24]。也有研究表明，在耐旱玉米基因型中的轉(zhuǎn)錄因子MYB33和MYB101能增強(qiáng)玉米的耐旱性[25]，在柑橘胚珠的研究中與信號(hào)傳遞相關(guān)的基因WRKY被驗(yàn)證在脅迫發(fā)生時(shí)表達(dá)下調(diào)[26]。本試驗(yàn)中ERF、WRKY、MYB是干旱脅迫處理下數(shù)量豐富、表達(dá)量變化最顯著的幾類轉(zhuǎn)錄因子，下調(diào)顯著的主要是ERF類轉(zhuǎn)錄因子，上調(diào)顯著的是MYB相關(guān)家族編碼基因，說明這些轉(zhuǎn)錄因子對(duì)干旱脅迫的應(yīng)答顯著，這與Chen等[27]的研究結(jié)果一致，原因可能是通過上調(diào)基因來實(shí)現(xiàn)葉片對(duì)水分刺激的應(yīng)答調(diào)節(jié)，通過下調(diào)基因數(shù)量的平衡來實(shí)現(xiàn)葉片對(duì)干旱脅迫的響應(yīng)[28]。

ABA的生物合成和分解代謝率決定ABA的積累以及對(duì)脅迫的反應(yīng)強(qiáng)度，盡管有各種各樣的途徑合成ABA，但類胡蘿卜素生物合成途徑在被子植物中是唯一已經(jīng)確定的途徑[29]。本研究中類胡蘿卜素生物合成通路中的10個(gè)ABA 8′-羥化酶上調(diào)表達(dá)，而ABA 8′-羥化酶屬于細(xì)胞色素P450的CYP707A家族，它是合成植物激素ABA的前體[30]且CYP707A8′-羥基化基因在擬南芥的研究中參與抗旱[31]、種子休眠和萌發(fā)[32]。ABA 8′-羥化酶可能參與合成類胡蘿卜素，并積累額外的ABA以響應(yīng)脅迫。ABA與赤霉素(GA)有拮抗作用，在干旱脅迫下，激素的比值會(huì)發(fā)生不同程度的改變，有研究表明，GA3/ABA值能體現(xiàn)玉米萌發(fā)期適應(yīng)脅迫的能力[33]。在本研究中，雙萜類生物合成途徑中的GA合成上游的2種關(guān)鍵酶均下調(diào)表達(dá)：貝殼杉烯酸氧化酶和貝殼杉烯氧化酶，這可能與ABA 8′-羥化酶基因的上調(diào)表達(dá)相互平衡，共同響應(yīng)干旱，但二者的具體變化情況還不明確。油菜素內(nèi)酯(BRs)是一種植物甾醇類激素，有研究發(fā)現(xiàn)，油菜素內(nèi)酯能提高擬南芥幼苗忍耐低溫的能力，還能增加植物的環(huán)境抗逆性[34]。而本研究中，油菜素內(nèi)固醇生物合成途徑中下調(diào)基因中有9個(gè)都是細(xì)胞色素P450家族的CYP734A1基因，而擬南芥中的CYP734A1能夠使體內(nèi)BR失活[35]，說明此時(shí)植物可能通過降低CYP734A1基因的表達(dá)量來調(diào)節(jié)BR的活性以應(yīng)對(duì)干旱脅迫，關(guān)于CYP734A1基因和其他未知功能基因與植物響應(yīng)干旱之間的具體關(guān)系都還有待更深入的研究。