基于SMRT 技術(shù)的水杉全長(zhǎng)轉(zhuǎn)錄組分析及基因功能注釋

劉小紅

（1.西華師范大學(xué)生命科學(xué)學(xué)院，四川 南充 637002；2.西南野生動(dòng)植物資源保護(hù)教育部重點(diǎn)實(shí)驗(yàn)室，四川 南充637002）

水杉（Metasequoia glyptostroboides） 屬孑遺植物，為杉科（Taxodiaceae）、水杉屬（Metasequoia） 的單一種, 在植物界有“活化石”之稱[1]。據(jù)史料記載，水杉起源于中生代的白堊紀(jì)，在新生代的第三紀(jì)時(shí)遍布北半球，但在第四紀(jì)的冰川期后全世界僅有極少量個(gè)體存活，直到20 世紀(jì)40 年代才在中國(guó)湖北省利川市的小河境內(nèi)被首次發(fā)現(xiàn)[2]。水杉生長(zhǎng)歷史悠久，對(duì)于研究古生物、古氣候、古地質(zhì)以及裸子植物的系統(tǒng)演化等都具有重要意義[3]。目前，水杉已被列入國(guó)家一級(jí)保護(hù)植物名錄。

水杉雖然具有重要的生物學(xué)作用，但關(guān)于它的研究主要集中在栽培引種[4-5]、生理生化[6-7]、系統(tǒng)發(fā)育[8]等方面，而在分子遺傳水平上的研究較少[9]。轉(zhuǎn)錄組測(cè)序可以獲得大量的基因信息[10-13]，能揭示潛在的代謝途徑和遺傳機(jī)制，可為其進(jìn)一步的分子生物學(xué)研究提供依據(jù)。

關(guān)于轉(zhuǎn)錄組測(cè)序，以Illumina/Solexa、Roche/454和ABI/SOLiD 為代表的二代測(cè)序平臺(tái)具有測(cè)序時(shí)間短、成本低、準(zhǔn)確性高、高通量等優(yōu)點(diǎn)，目前在差異基因表達(dá)方面被廣泛應(yīng)用[14-15]。但二代測(cè)序讀長(zhǎng)較短，在沒(méi)有參考基因組的情況下難以獲得基因的全長(zhǎng)序列信息。而近幾年發(fā)展起來(lái)的三代測(cè)序技術(shù)已經(jīng)成為更好的選擇，如PacBio 測(cè)序平臺(tái)的單分子測(cè)序SMRT（single-molecule real-time）技術(shù)已經(jīng)成為獲得全長(zhǎng)轉(zhuǎn)錄序列的首要選擇。三代測(cè)序技術(shù)的主要優(yōu)點(diǎn)是讀長(zhǎng)較長(zhǎng)（平均讀長(zhǎng)可達(dá)20 kb）, 對(duì)于逆轉(zhuǎn)錄生成的全長(zhǎng)轉(zhuǎn)錄本不需要對(duì)其進(jìn)行片段化處理，可以直接作單分子測(cè)序獲得全長(zhǎng)序列信息[16-17]。全長(zhǎng)轉(zhuǎn)錄組測(cè)序技術(shù)已被用于許多植物的全長(zhǎng)轉(zhuǎn)錄組分析[18-20]，但在水杉植物中還少見(jiàn)相關(guān)報(bào)道。基于此，本研究擬以來(lái)自湖北省利川市的原生水杉種為材料，采用SMRT 三代測(cè)序技術(shù)對(duì)全長(zhǎng)轉(zhuǎn)錄本進(jìn)行測(cè)序分析，以獲得相應(yīng)的遺傳信息，旨在為后續(xù)基因克隆及基因功能鑒定提供依據(jù)。

1 材料與方法

1.1 材料

選取來(lái)自湖北省利川市林業(yè)科學(xué)研究所（108°96′E，30°28′N(xiāo)） 的原生水杉種為材料，在植株長(zhǎng)至10 cm 左右時(shí)，分別取其根、莖、葉，液氮速凍保存。

1.2 方法

1.2.1 總RNA 的提取及其質(zhì)量檢測(cè)

水杉根、莖、葉的總RNA 的提取都用RNeasy Plus Mini Kit（Qiagen, Valencia, CA, USA） 完成。對(duì)提取得到的總RNA 先用Nanodrop 2000 和Agilent 2100 分析其濃度，測(cè)定OD260/280、OD260/230、25S/18S和RIN 值，再用1.0%的瓊脂糖凝膠電泳分析是否有DNA 污染和RNA 降解現(xiàn)象。

1.2.2 SMRT 測(cè)序

將檢測(cè)合格的水杉根、莖和葉的總RNA 等量混勻，經(jīng)mRNA 純化、反轉(zhuǎn)錄及SMRT 測(cè)序后得到全長(zhǎng)轉(zhuǎn)錄組原始數(shù)據(jù)（raw data）。該工作由北京諾禾致源科技股份有限公司的PacBio 測(cè)序平臺(tái)完成。

1.2.3 raw data 及subreads 統(tǒng)計(jì)

原始數(shù)據(jù)經(jīng)質(zhì)量控制處理后進(jìn)行subreads 的統(tǒng)計(jì)分析，包括有效插入片段subreads 數(shù)據(jù)量大小、有效插入片段subreads 條數(shù)、平均subreads 長(zhǎng)度及N50（將得到的subreads 按照長(zhǎng)度從大到小排序，依次累加subreads 的長(zhǎng)度，直至其長(zhǎng)度不小于總subreads 長(zhǎng)度的50%）的長(zhǎng)度。繪制subreads 長(zhǎng)度分布圖。

1.2.4 全長(zhǎng)轉(zhuǎn)錄本分析

選用SMRT Link 5.1 對(duì)全長(zhǎng)轉(zhuǎn)錄本進(jìn)行統(tǒng)計(jì)，按默認(rèn)參數(shù)設(shè)置，統(tǒng)計(jì)以下項(xiàng)目：環(huán)形一致性序列（CCS） 數(shù)量；全長(zhǎng)非嵌合reads 數(shù)量；全長(zhǎng)非嵌合reads 平均長(zhǎng)度；聚類分析之后得到的一致性序列的reads 數(shù)量。

1.2.5 轉(zhuǎn)錄本去冗余

利用CD-HIT[21]通過(guò)序列比對(duì)聚類法去除冗余，輸出1 個(gè)非冗余的序列文件，去冗余后的基因即為unigenes。對(duì)轉(zhuǎn)錄本去冗余前后長(zhǎng)度頻數(shù)分布情況和含有相同轉(zhuǎn)錄本拷貝數(shù)的基因數(shù)量進(jìn)行統(tǒng)計(jì)。

1.2.6 基因的功能注釋

為了得到全面的基因功能信息，對(duì)去冗余之后的基因進(jìn)行數(shù)據(jù)庫(kù)（包括NR、NT、Pfam、KOG、Swiss-Prot、KEGG 以及GO）功能注釋，分析這些基因在數(shù)據(jù)庫(kù)中的注釋情況，并利用NR 數(shù)據(jù)庫(kù)繪制物種分布圖和基因功能的GO 分類統(tǒng)計(jì)圖。

1.2.7 基因的結(jié)構(gòu)分析

利用Angel[22]進(jìn)行CDS（coding sequence） 預(yù)測(cè)分析。該軟件具有無(wú)錯(cuò)和容錯(cuò)2 種預(yù)測(cè)模式。本研究中，選用默認(rèn)的容錯(cuò)模式，對(duì)CDS 序列的長(zhǎng)度及其分布進(jìn)行預(yù)測(cè)，繪制序列長(zhǎng)度分布圖。運(yùn)用iTAK[23]對(duì)去冗余的unigenes 進(jìn)行植物轉(zhuǎn)錄因子預(yù)測(cè)，并將注釋到轉(zhuǎn)錄本數(shù)量最多的轉(zhuǎn)錄因子家族進(jìn)行柱形圖展示。

2 結(jié)果與分析

2.1 總RNA 的分離純化

水杉總RNA 分離純化后，用Nanodrop 2000和Agilent 2100 檢測(cè)的結(jié)果顯示：在最后加入溶劑都為30.00 μL 的條件下，葉的總RNA 濃度最高，達(dá)301.70 ng/μL，其后依次為莖的（268.80 ng/μL） 和根的（189.70 ng/μL），三者的總量分別為9.05、8.06、5.69 μg。根、莖和葉的總RNA 的OD260/280值分別為2.04、2.05、2.11，OD260/230值分別為1.25、1.10、1.47，25S/18S 值分別為1.50、1.50 和1.40，RIN 值分別為8.50、8.70、7.70。瓊脂糖凝膠電泳結(jié)果顯示，每個(gè)樣品都有2 條明亮的條帶，沒(méi)有拖尾現(xiàn)象，表明提取的RNA 質(zhì)量較好，沒(méi)有DNA 污染和RNA降解現(xiàn)象。上述檢測(cè)結(jié)果表明，提取的RNA 樣品能滿足建庫(kù)測(cè)序?qū)|(zhì)量的要求。

2.2 raw data 及subreads 結(jié)果分析

將提取的水杉根、莖、葉3 個(gè)部位的總RNA 樣品等量混勻，經(jīng)mRNA 純化、反轉(zhuǎn)錄及SMRT 測(cè)序，對(duì)原始數(shù)據(jù)進(jìn)行質(zhì)量控制處理后共得5 914 711 個(gè)subreads 和14.0 Gb subread base，平均subreads 長(zhǎng)度為2368 nt，N50 為2569 nt。subreads 的長(zhǎng)度分布如圖1 所示。在長(zhǎng)度為100～10 000 nt 的范圍內(nèi)，2300 nt附近的subreads 最多，而超過(guò)7000 nt 的則非常少。

圖1 水杉全長(zhǎng)轉(zhuǎn)錄組subreads 的長(zhǎng)度分布Fig.1 The length distribution of the s ubreads from the full-length transfriptome of Metasequoia glyptostroboides

2.3 全長(zhǎng)轉(zhuǎn)錄本數(shù)據(jù)統(tǒng)計(jì)分析

采用SMRT 技術(shù)測(cè)序后共獲得339 296 個(gè)CCS序列，其中用于后續(xù)試驗(yàn)分析的全長(zhǎng)非嵌合reads共計(jì)236 130 個(gè)，平均長(zhǎng)度為2598 nt，對(duì)其作一致性序列處理，最后得到一致性reads 97 626 個(gè)。

2.4 轉(zhuǎn)錄本的去冗余結(jié)果分析

對(duì)97 626個(gè)一致性reads進(jìn)一步作去冗余處理，去冗余前后的結(jié)果如表1 所示。小于500 bp 的轉(zhuǎn)錄本和基因數(shù)量較少，分別僅有739、363 個(gè)；在1～3 kb 的轉(zhuǎn)錄本和基因數(shù)量較多，約占總量的96.92%。去冗余后最終得到的基因數(shù)（unigenes 的數(shù)量） 為61 057 個(gè)，占總轉(zhuǎn)錄本的62.54%。

表1 水杉轉(zhuǎn)錄本去冗余前后的長(zhǎng)度分布Table 1 Length distribution of the transcripts before and af ter redundancy removal in Metasequoia glyptostroboides

對(duì)去冗余的61 057 個(gè)基因作轉(zhuǎn)錄本拷貝數(shù)分析，結(jié)果拷貝數(shù)最少的為1 個(gè)拷貝，最多的為10個(gè)拷貝，其中1 個(gè)拷貝的基因數(shù)達(dá)51 523 個(gè)，占總量的84.39%；含9 個(gè)拷貝的基因數(shù)最少，為134個(gè)，僅占總量的0.22%。

2.5 基因功能注釋結(jié)果分析

2.5.1 注釋基因數(shù)量

利用NR、Swiss-Prot、KEGG、KOG、GO、NT 和Pfam 7 個(gè)數(shù)據(jù)庫(kù)對(duì)水杉全長(zhǎng)unigenes 進(jìn)行功能注釋，共有54 099 個(gè)基因被成功注釋，被注釋的基因在7 個(gè)數(shù)據(jù)庫(kù)中的數(shù)量依次為52 714、42 752、51 073、32 559、24 045、31 455、24 045 個(gè)，在7個(gè)數(shù)據(jù)庫(kù)中均被注釋的有12 043 個(gè)。選取常用的NR、KEGG、KOG、GO 和NT 5 個(gè)數(shù)據(jù)庫(kù)的注釋情況繪制韋恩圖（圖2）。從圖2 可知，在NR 數(shù)據(jù)庫(kù)中注釋的基因數(shù)最多，其次為KEGG 數(shù)據(jù)庫(kù)，最少的為GO 數(shù)據(jù)庫(kù)，在5 個(gè)數(shù)據(jù)庫(kù)中均有注釋的基因數(shù)為12 209 個(gè)。

圖2 水杉基因功能注釋韋恩圖Fig.2 The Venn diagram of gene function annotation in Metasequoia glyptostroboides

2.5.2 NR 數(shù)據(jù)庫(kù)注釋

通過(guò)與NR 數(shù)據(jù)庫(kù)進(jìn)行比對(duì)注釋，共有52 714個(gè)基因被注釋到568 個(gè)物種中。匹配基因數(shù)最多的是北美云杉（Picea sitchensis），達(dá)12 519 個(gè)；其次為無(wú)油樟（Amborella trichopoda） ，5888 個(gè)；排名第三的為蓮（Nelumbo nucifera），3417 個(gè)；最少的為小油桐（Jatropha curcas），僅有442 個(gè)。

2.5.3 GO 分類

對(duì)上述unigenes 進(jìn)行GO 注釋，將注釋成功的基因按照GO 3 個(gè)大類的下一級(jí)進(jìn)行分類，結(jié)果顯示：與細(xì)胞組分（cellular component）有關(guān)的基因有21 966 個(gè)，歸屬于18 個(gè)下一級(jí)分類，平均每級(jí)分類為1220 個(gè)。與生物過(guò)程（biological process）有關(guān)的基因數(shù)為42 898 個(gè)，歸屬于25 個(gè)下一級(jí)分類，平均每級(jí)為1716 個(gè)。與分子功能（molecular function）相關(guān)的基因有28 015 個(gè)，歸屬于11 個(gè)下一級(jí)分類，平均每級(jí)為2547 個(gè)。

2.6 基因結(jié)構(gòu)分析

2.6.1 CDS 預(yù)測(cè)

CDS 預(yù)測(cè)結(jié)果如圖3 所示。共有61 259 個(gè)CDS序列，其長(zhǎng)度不一，最長(zhǎng)的達(dá)6 477 nt，最短的僅為144 nt，平均長(zhǎng)度約為679 nt?？偟内厔?shì)是轉(zhuǎn)錄本序列越長(zhǎng)，則轉(zhuǎn)錄本數(shù)量越少。

圖3 水杉基因的CDS 長(zhǎng)度分布情況Fig. 3 CDS length distribution of the genes in Metasequoia glyptostroboides

2.6.2 轉(zhuǎn)錄因子分析

在61 057 個(gè)unigenes 中，預(yù)測(cè)到共有2386 個(gè)表達(dá)產(chǎn)物為轉(zhuǎn)錄因子，這些轉(zhuǎn)錄因子可以分成29個(gè)家族，其含有的轉(zhuǎn)錄因子成員數(shù)最多的為C3H 家族，有234 個(gè)；其次為GRAS 家族，有194 個(gè)；最少的為HSF 家族，僅有28 個(gè)。

3 結(jié)論與討論

全長(zhǎng)轉(zhuǎn)錄組測(cè)序技術(shù)可有效獲取物種轉(zhuǎn)錄組信息且可信度高，目前該技術(shù)被廣泛應(yīng)用于許多物種[24-25]，有助于了解沒(méi)有參考基因組的非模式生物的遺傳信息[26]。本研究中，基于PacBio 平臺(tái)的SMRT 測(cè)序技術(shù)被用于分析水杉幼苗全株的全長(zhǎng)轉(zhuǎn)錄組測(cè)序，共獲得了5 914 711個(gè)subreads和14.0 Gb的數(shù)據(jù)量。經(jīng)過(guò)數(shù)據(jù)優(yōu)化和去冗余處理，最后獲得61 057 個(gè)基因的全長(zhǎng)cDNA 序列信息。本研究結(jié)果與在野草莓（Fragaria vesca）[27]、紫苜蓿（Medicago sativa）[28]和高山杜鵑（Rhododendron lapponicum）[29]等中的研究結(jié)果類似。本研究中，轉(zhuǎn)錄本的平均長(zhǎng)度達(dá)2598 nt，遠(yuǎn)高于二代測(cè)序的平均讀長(zhǎng)[30-31]。本研究結(jié)果還表明，不同的基因在轉(zhuǎn)錄組測(cè)序分析中轉(zhuǎn)錄本拷貝數(shù)可能不一樣，其中具有單拷貝轉(zhuǎn)錄本的最多，占全部轉(zhuǎn)錄本的84.39%，總的趨勢(shì)是拷貝數(shù)越多，對(duì)應(yīng)的基因數(shù)量越少。

水杉雖然不是模式植物，但可以基于已有參考植物的數(shù)據(jù)庫(kù)信息對(duì)本研究獲得的61 057 個(gè)全長(zhǎng)基因進(jìn)行功能注釋，結(jié)果共有54 099 個(gè)基因被注釋到7 個(gè)數(shù)據(jù)庫(kù)中，被注釋的基因占88.60%，其中有12 043個(gè)同時(shí)被7個(gè)數(shù)據(jù)庫(kù)注釋，占總量的19.72%。通過(guò)NR 數(shù)據(jù)庫(kù)的注釋，發(fā)現(xiàn)水杉和北美云杉之間的親緣關(guān)系較近。此外，通過(guò)GO 數(shù)據(jù)庫(kù)對(duì)基因功能進(jìn)行分類統(tǒng)計(jì)，并將這些基因分成了生物過(guò)程、細(xì)胞組分及分子功能3 大類。在基因結(jié)構(gòu)分析中，CDS 預(yù)測(cè)結(jié)果顯示，轉(zhuǎn)錄本序列越長(zhǎng)，則其在細(xì)胞中的拷貝數(shù)越少，這與前人[20,24]在其他物種作全長(zhǎng)轉(zhuǎn)錄組分析的結(jié)果一致。轉(zhuǎn)錄因子分析發(fā)現(xiàn)，獲得的部分基因的表達(dá)產(chǎn)物其實(shí)是與細(xì)胞中基因轉(zhuǎn)錄調(diào)控密切相關(guān)的轉(zhuǎn)錄因子，與細(xì)胞中很多基因能否轉(zhuǎn)錄及轉(zhuǎn)錄強(qiáng)度密切相關(guān)。該結(jié)果為分析水杉基因的功能提供了分子基礎(chǔ)。一小部分基因不能被注釋，可能是因?yàn)檫@些unigenes 是一些新基因，在數(shù)據(jù)庫(kù)中還沒(méi)找到與此相似的轉(zhuǎn)錄本序列，有待于對(duì)這些全長(zhǎng)轉(zhuǎn)錄本對(duì)應(yīng)的基因的結(jié)構(gòu)和功能作進(jìn)一步研究。