毛常麗，李 玲，張鳳良，李小琴，楊 湉，吳 裕

（云南省熱帶作物科學(xué)研究所，云南景洪666100）

低溫寒害一直是限制我國(guó)橡膠樹(shù)種植面積擴(kuò)大和越冬安全的主要因子，抗寒性狀也一直是我國(guó)橡膠樹(shù)良種培育的主要選擇指標(biāo)之一。根據(jù)《NY/T607-2002橡膠樹(shù)選育種技術(shù)規(guī)程》中的育種程序，選出一個(gè)優(yōu)良無(wú)性系一般需要25～30年。云南省熱帶作物科學(xué)研究所培育的三倍體無(wú)性系“云研77-2”和“云研77-4”具有抗寒、高產(chǎn)特性，在云南省累積推廣種植160 000 hm2。這兩個(gè)品種從1974雜交到2000年良種審定經(jīng)歷了26年（后經(jīng)細(xì)胞學(xué)鑒定為三倍體）。一般認(rèn)為，樹(shù)木三倍體具有良好的速生、抗逆和豐產(chǎn)性，實(shí)際上只能說(shuō)從三倍體中更有機(jī)會(huì)獲得性狀突出的個(gè)體。橡膠樹(shù)自然三倍體的發(fā)生率特別低，要想通過(guò)傳統(tǒng)的育種程序選出優(yōu)良三倍體就更不容易。

轉(zhuǎn)錄組學(xué)是從RNA水平研究基因表達(dá)的情況，是研究細(xì)胞表型和功能的重要手段。轉(zhuǎn)錄組測(cè)序技術(shù)（RNA-Seq）通過(guò)深度測(cè)序?qū)悠分械霓D(zhuǎn)錄因子進(jìn)行分析，能夠反映出轉(zhuǎn)錄水平上的基因表達(dá)情況，其已被廣泛應(yīng)用于植物種質(zhì)評(píng)價(jià)、性狀篩選等方面。篤斯越橘（Vaccinium uliginosum）、厚樸（Magnolia officinalis）、枇杷（Eriobotrya japonica）、櫻桃（Prunus avium）［1-4］等木本物種的轉(zhuǎn)錄組學(xué)抗逆性研究為橡膠樹(shù)抗逆性轉(zhuǎn)錄組學(xué)的研究提供了新的思路。

目前，橡膠樹(shù)轉(zhuǎn)錄組學(xué)主要基于第二代高通量測(cè)序平臺(tái)的RNA-Seq技術(shù)對(duì)病害、產(chǎn)膠、生殖發(fā)育等［5-7］方面進(jìn)行研究，且用于測(cè)序的材料為二倍體，而對(duì)三倍體材料進(jìn)行的轉(zhuǎn)錄組測(cè)序還未見(jiàn)報(bào)道?；谔镩g鑒定［8］、生理指標(biāo)測(cè)定［9-11］等結(jié)果已確定三倍體品種云研77-2和云研77-4抗寒性比其母本GT1強(qiáng)，但抗寒性增強(qiáng)是來(lái)自于母本的2n配子中抗寒效應(yīng)基因的加倍或是其他原因則還未見(jiàn)報(bào)道。如果能探清其抗寒機(jī)制，將有助于三倍體抗寒品種的培育。本研究以云研77-4為材料作全長(zhǎng)轉(zhuǎn)錄組測(cè)序及分析，旨在得到完整的轉(zhuǎn)錄組序列，為橡膠樹(shù)三倍體分子輔助選擇及抗寒三倍體育種奠定理論基礎(chǔ)。

1 材料和方法

1.1 實(shí)驗(yàn)材料

材料選擇為橡膠樹(shù)三倍體無(wú)性系“云研77-4”，其親本為GT1×PR107；苗木以GT1開(kāi)放授粉種子苗為砧木芽接培育；選擇長(zhǎng)到第3蓬葉且頂蓬葉穩(wěn)定，整株無(wú)病蟲(chóng)害，長(zhǎng)勢(shì)一致的苗木，準(zhǔn)備實(shí)施低溫脅迫處理。

1.2 方法

在4℃人工氣候箱下處理2 h、6 h和20 h，對(duì)照為不處理（0 h），每處理設(shè)3次重復(fù)。處理完后剪取第2個(gè)葉蓬中下部復(fù)葉的中間小葉，并迅速用液氮凍干后轉(zhuǎn)移至-80℃冰箱保存?zhèn)溆?。先?duì)各樣品分別提取RNA，再將不同處理的RNA混合，對(duì)混合樣進(jìn)行測(cè)序。樣品提取、文庫(kù)構(gòu)建及質(zhì)控、轉(zhuǎn)錄組測(cè)序及分析委托北京百邁客生物科技有限公司進(jìn)行。參考已報(bào)道的橡膠樹(shù)全基因組序列對(duì)本研究的轉(zhuǎn)錄組序列進(jìn)行注釋。

2 結(jié)果與分析

2.1 數(shù)據(jù)產(chǎn)出及統(tǒng)計(jì)

云研77-4冷處理后的樣品經(jīng)PacBio三代轉(zhuǎn)錄組測(cè)序，共獲得41.85 Gb clean data的數(shù)據(jù)，按照f(shuō)ull passes＞=3且序列準(zhǔn)確性大于0.9的篩選條件，從原始數(shù)據(jù)中提取到環(huán)形一致序列（Circular Consensus Sequencing，CCS）504 032條，平均每條序列的長(zhǎng)度為2 157 bp。對(duì)CCS序列進(jìn)行長(zhǎng)度分布分析，得到其長(zhǎng)度分布范圍主要為1 800～2 200 bp（圖1）。過(guò)濾掉CCS序列中不包含5’引物、3’引物及polyA尾的非全長(zhǎng)序列，共得到全長(zhǎng)非嵌合序 列（Full-length Non-Chimeric sequence,FLCN）413 878條，占總CCS序列的82.11%。

圖1 橡膠樹(shù)三倍體三代轉(zhuǎn)錄組全長(zhǎng)測(cè)序CCS長(zhǎng)度分布

2.2 FLCN的一致序列（consensus isoform）分析

使用SMRTLink軟件中的IsoSeq模塊對(duì)FLNC序列進(jìn)行相似性序列聚類(lèi)，得到一致序列157 638條，從其中篩選出準(zhǔn)確度大于99%的高質(zhì)量轉(zhuǎn)錄本154 484條，平均每條一致序列長(zhǎng)度為2 181 bp，其consensus isoform序列長(zhǎng)度分布如圖2。用GMAP（Genomic Mapping and Alignment Program）將校正后的一致序列與橡膠樹(shù)參考基因組進(jìn)行序列比對(duì)，用cDNA_Cupcake軟件將比對(duì)后的序列進(jìn)行去冗余處理，過(guò)濾identity小于0.9，coverage小于0.85的序列，共得到非冗余轉(zhuǎn)錄本序列96 437條。對(duì)非冗余的轉(zhuǎn)錄本進(jìn)行可變剪接分析，檢測(cè)到基因位點(diǎn)23 307個(gè)，其中新基因位點(diǎn)4 288個(gè)，新發(fā)現(xiàn)轉(zhuǎn)錄本70 055個(gè)。為了探究更多的轉(zhuǎn)錄組信息，同時(shí)對(duì)三倍體橡膠樹(shù)中的融合轉(zhuǎn)錄本進(jìn)行了分析，共得到融合轉(zhuǎn)錄本1 983個(gè)。

圖2 consensus isoform長(zhǎng)度分布

2.3 新轉(zhuǎn)錄本分析

2.3.1 SSR分析

從新發(fā)現(xiàn)的轉(zhuǎn)錄本中篩選500 bp以上的轉(zhuǎn)錄本，利用MISA軟件進(jìn)行SSR分析，結(jié)果見(jiàn)表1。對(duì)檢測(cè)到的SSR進(jìn)行類(lèi)型劃分，劃分為完美單堿基重復(fù)（p1）、完美雙堿基重復(fù)（p2）、完美三堿基重復(fù)（p3）、完美四堿基重復(fù)（p4）、完美五堿基重復(fù)（p5）、完美六堿基重復(fù)（p6）和混合SSR（c，即包含至少兩個(gè)完美SSR，且之間距離小于100 bp），并對(duì)其進(jìn)行密度統(tǒng)計(jì)，統(tǒng)計(jì)結(jié)果見(jiàn)圖3。

由表1、圖3可知，橡膠樹(shù)新轉(zhuǎn)錄本中存在6種類(lèi)型的SSR，每種重復(fù)類(lèi)型的數(shù)量差異較大，主要以單堿基重復(fù)為主，且在所有的堿基重復(fù)中，以完美單堿基重復(fù)類(lèi)型最多。

表1 橡膠樹(shù)新轉(zhuǎn)錄本SSR分析

圖3 橡膠樹(shù)三代轉(zhuǎn)錄組新轉(zhuǎn)錄本SSR類(lèi)型分布圖

2.3.2 新基因編碼區(qū)序列預(yù)測(cè)

對(duì)可變剪接分析中得到的新轉(zhuǎn)錄本使用TransDecoder軟件［12］預(yù) 測(cè) 其編 碼 區(qū) 序 列（Coding Sequence，CDS）及其對(duì)應(yīng)的氨基酸序列，共獲得開(kāi)放閱讀框（Open Reading Frame,ORF）66 182條，其中完整ORF57 516條。在獲得的ORF中，0～500 aa的有54 475條，占82.31%，500～1 000 aa的有10 537條，占15.91%，1 000～2 300 aa的有1 167條，占1.76%（圖4）。

圖4 橡膠樹(shù)轉(zhuǎn)錄組預(yù)測(cè)的CDS編碼蛋白長(zhǎng)度分布

2.3.3 長(zhǎng)鏈非編碼RNA（Longnon-coding RNA，lncRNA）預(yù)測(cè)及轉(zhuǎn)錄因子分析

通過(guò)對(duì)轉(zhuǎn)錄本中不編碼蛋白的lncRNA進(jìn)行編碼潛能篩選，從而判定該轉(zhuǎn)錄本是否為lncRNA。利用CPC［13］、CNCI［14］、pfam、CPAT［15］四種分析方法對(duì)新發(fā)現(xiàn)的轉(zhuǎn)錄本進(jìn)行l(wèi)ncRNA預(yù)測(cè)，分別預(yù)測(cè)到5 598、7 164、15 207、17 488個(gè)lncRNA，其中4種方法均預(yù)測(cè)到的lncRNA共3 575個(gè)（圖5）。同時(shí)，根據(jù)lncRNA在參考基因組注釋信息上的位置對(duì)lncRNA進(jìn)行分類(lèi)繪圖，結(jié)果見(jiàn)圖6。由圖6可知，本研究材料中的lncRNA主要以lincRNA為主，占53.19%，其次為sense_lncRNA，占27.30%。

圖5 四種篩選方法維恩圖

圖6 lncRNA位置分類(lèi)

LncRNA的存在為下一步研究橡膠樹(shù)表觀遺傳、轉(zhuǎn)錄水平及轉(zhuǎn)錄后水平調(diào)控基因的表達(dá)提供了新的思路。

轉(zhuǎn)錄因子是指能夠結(jié)合在某基因上游特異核苷酸序列上的蛋白質(zhì)，這些蛋白質(zhì)可以調(diào)控RNA聚合酶與DNA模板的結(jié)合，從而調(diào)控基因的轉(zhuǎn)錄。使用iTAK軟件［16］對(duì)三倍體橡膠樹(shù)三代測(cè)序轉(zhuǎn)錄因子進(jìn)行預(yù)測(cè)，共預(yù)測(cè)到轉(zhuǎn)錄因子12 064個(gè)，且預(yù)測(cè)到的轉(zhuǎn)錄因子包括了在其他物種中有報(bào)道的RLK-Pelle_DLSV、C3H、NAC、bHLH、MYBrelated及RLK-Pelle_LRR-VIII-1等 類(lèi) 型。C3H、bHLH、MYB-related等轉(zhuǎn)錄因子在植物抗逆、生長(zhǎng)發(fā)育及有效成分合成方面均有報(bào)道，這為今后開(kāi)展橡膠樹(shù)相關(guān)研究提供了重要信息。同時(shí)，為了更直觀地呈現(xiàn)全長(zhǎng)轉(zhuǎn)錄組測(cè)序所得到的注釋序列在全基因組上的分布密度及全長(zhǎng)轉(zhuǎn)錄組測(cè)序的完整性，我們將注釋到的轉(zhuǎn)錄本與全基因組進(jìn)行了Circos可視化繪圖（圖7）。從圖7可知，本研究材料全長(zhǎng)轉(zhuǎn)錄組序列中的基因密度、轉(zhuǎn)錄本密度和全基因組中的基因密度、轉(zhuǎn)錄本密度相當(dāng)，lncRNA在各條染色體上均有分布。

圖7 基因組水平Cicros可視化圖

2.4 新轉(zhuǎn)錄本功能注釋

為探清三倍體橡膠樹(shù)新轉(zhuǎn)錄本的功能，用BLAST軟件（version 2.2.26）［17］將得到的新轉(zhuǎn)錄本與NR［18］、Swissprot［19］、GO［20］、COG［21］、KOG［22］、Pfam［23］、KEGG［24］數(shù)據(jù)庫(kù)比對(duì)，進(jìn)行新轉(zhuǎn)錄本功能注釋?zhuān)沧⑨尩焦δ芑?6 684個(gè)。7個(gè)數(shù)據(jù)庫(kù)注釋到的轉(zhuǎn)錄本數(shù)量見(jiàn)表2。

從表2可以看出，7個(gè)功能數(shù)據(jù)庫(kù)中注釋到的轉(zhuǎn)錄本長(zhǎng)度大部分在1 000 bp以上，這也進(jìn)一步說(shuō)明全長(zhǎng)轉(zhuǎn)錄組測(cè)序能測(cè)得較長(zhǎng)的序列。綜合7個(gè)功能數(shù)據(jù)庫(kù)中注釋到的轉(zhuǎn)錄本，其功能分類(lèi)主要涉及細(xì)胞組分、分子功能及生物學(xué)過(guò)程三大方面，這為抗寒性相關(guān)基因的篩選奠定了基礎(chǔ)。

表2 7個(gè)數(shù)據(jù)庫(kù)注釋到的轉(zhuǎn)錄本數(shù)量

3 討論

轉(zhuǎn)錄組是研究動(dòng)植物生命過(guò)程的有效方法之一，隨著測(cè)序手段的不斷更新，近幾年發(fā)展起來(lái)的三代全長(zhǎng)轉(zhuǎn)錄組測(cè)序因其能讀取更長(zhǎng)的序列而被廣泛應(yīng)用于植物的生長(zhǎng)發(fā)育、非生物脅迫等領(lǐng)域。筆者所在單位云南省熱帶作物科學(xué)研究所在國(guó)際上首次獲得了達(dá)到染色體級(jí)別的高質(zhì)量巴西橡膠樹(shù)優(yōu)良品種GT1的參考基因組序列［25］這為本研究的轉(zhuǎn)錄組注釋、組裝提供了重要參考。

橡膠樹(shù)轉(zhuǎn)錄組測(cè)序已經(jīng)應(yīng)用于自根幼態(tài)無(wú)性系體胚發(fā)育［26］、橡膠樹(shù)叢枝病機(jī)理［5］及苗期生長(zhǎng)優(yōu)勢(shì)［27］等方面，并從中篩選了相關(guān)的調(diào)控基因，為橡膠樹(shù)分子輔助育種奠定了重要的理論基礎(chǔ)，但是現(xiàn)有的報(bào)道均基于Illumina HiSeq的二代轉(zhuǎn)錄組測(cè)序技術(shù)，且測(cè)序的樣品均為二倍體。本研究通過(guò)PacBio三代轉(zhuǎn)錄組測(cè)序技術(shù)對(duì)三倍體橡膠樹(shù)抗寒無(wú)性系進(jìn)行測(cè)序，共產(chǎn)生41.85 Gb clean data的數(shù)據(jù)，分析得到高質(zhì)量轉(zhuǎn)錄本154 484條，平均每條序列長(zhǎng)度為2 181 bp，數(shù)據(jù)得量、轉(zhuǎn)錄本總數(shù)以及單條轉(zhuǎn)錄本序列長(zhǎng)度均大于先前報(bào)道的橡膠樹(shù)二代轉(zhuǎn)錄組數(shù)據(jù)［28-29］。將得到的數(shù)據(jù)分析后得到新轉(zhuǎn)錄本70 055條。對(duì)94 084條轉(zhuǎn)錄本進(jìn)行分析，發(fā)現(xiàn)了潛在的80 432個(gè)SSR分子標(biāo)記，高于Wirulda等［30］用454高 通 量 測(cè) 序 技 術(shù) 獲 得 的24 000個(gè)SSRs，這些標(biāo)記可廣泛應(yīng)用于橡膠樹(shù)種質(zhì)資源評(píng)價(jià)、分子輔助育種、抗逆性鑒定等。利用CPC、CNCI、pfam、CPAT四種方法預(yù)測(cè)到共有的lncRNA 3 575個(gè)，同時(shí)對(duì)新轉(zhuǎn)錄本進(jìn)行了NR、Swissprot、GO、COG、KOG、Pfam、KEGG功能注釋?zhuān)沧⑨尩焦δ芑?6 684個(gè)，高于前人注釋到的37 432個(gè)功能基因［28］。轉(zhuǎn)錄因子可與基因5’端上游特定序列專(zhuān)一性結(jié)合，從而保證目的基因在特定的時(shí)間與空間進(jìn)行表達(dá)。本研究中共預(yù)測(cè)到轉(zhuǎn)錄因子12 064個(gè)，其中包含了與抗逆性相關(guān)的MYB、C3H轉(zhuǎn)錄因子，下一步將從中篩選出與三倍體橡膠樹(shù)抗寒性相關(guān)的轉(zhuǎn)錄因子，為研究三倍體抗寒性提供參考。從全長(zhǎng)轉(zhuǎn)錄組數(shù)據(jù)和全基因組數(shù)據(jù)Circos可視化圖可以看出，全長(zhǎng)轉(zhuǎn)錄組數(shù)據(jù)在全基因組測(cè)序得到的數(shù)據(jù)覆蓋率相當(dāng)，說(shuō)明全長(zhǎng)轉(zhuǎn)錄組序列可靠，可為下一步二代轉(zhuǎn)錄組序列的拼接、組裝及注釋提供參考。

4 結(jié)論

本研究基于三代測(cè)序技術(shù)建立了三倍體橡膠樹(shù)的轉(zhuǎn)錄組數(shù)據(jù)庫(kù)，并對(duì)序列進(jìn)行分析，為后續(xù)的橡膠樹(shù)分子輔助育種、二代轉(zhuǎn)錄組數(shù)據(jù)組裝，特別是三倍體橡膠樹(shù)的轉(zhuǎn)錄組分析提供重要的基礎(chǔ)，下一步將繼續(xù)從lncRNA、轉(zhuǎn)錄因子等方面分析三倍體橡膠樹(shù)抗寒性標(biāo)記及強(qiáng)抗寒性的分子機(jī)制，為三倍體選育種提供有用標(biāo)記及參考。