基于高通量測序的女貞果實轉錄組研究△

徐德宏，羅月芳，江靈敏，孟蕾，譚朝陽

(湖南中醫(yī)藥大學 藥學院 生物工程實驗室，湖南 長沙 410208)

女貞LigustrumlucidumAit.為木犀科女貞屬植物，其干燥成熟果實(即女貞子)是一種常用中藥，具有滋補肝腎、明目烏發(fā)的功能，用于肝腎陰虛、眩暈耳鳴、腰膝酸軟、須發(fā)早白、目暗不明、內(nèi)熱消渴、骨蒸潮熱[1]。經(jīng)化學成分分析和活性檢測發(fā)現(xiàn)，女貞果實中含有多種具有藥理活性的次級代謝物[2-3]，如萜類、黃酮類等，因此近年來女貞果實在醫(yī)藥及保健方面?zhèn)涫苋藗兊年P注。目前，女貞果實的研究主要集中在次級代謝物的成分分析、提取和活性檢測，而關于其功能基因發(fā)現(xiàn)和注釋、物質(zhì)合成相關途徑分析、分子標記開發(fā)等轉錄組學方面的研究卻未有報道，故現(xiàn)階段有必要采用適當?shù)姆椒▽σ陨蟽?nèi)容展開深入系統(tǒng)的研究。

近年來，各種組學研究相繼興起，其中第二代高通量測序技術特別是RNA測序技術已廣泛應用于生物體轉錄組基因表達分析，采用該技術能全面快速地獲取研究對象在某一狀態(tài)下基因組轉錄信息，從中挖掘重要的功能基因，揭示不同生物學性狀的分子機制[4-5]。Illumina/Solexa是目前轉錄組測序中最常使用的一種RNA高通量測序技術，與傳統(tǒng)的EST sanger測序法相比，具有高通量、低成本和高靈敏度等無法比擬的優(yōu)勢，同時可對低豐度的基因表達進行檢測[6]。此外，隨著生物信息學的不斷發(fā)展，特別是序列組裝軟件的開發(fā)，主要蛋白質(zhì)數(shù)據(jù)庫的更新和轉錄組分析方法的建立[7-9]，為非模式生物進行轉錄組高通量測序、從頭拼接組裝和unigene的功能注釋提供了極大的方便，因而Illumina/Solexa第二代高通量測序技術已廣泛應用到了非模式生物的轉錄組研究中。

本研究將利用Illumina/Solexa Hiseq2000測序平臺，對女貞果實進行轉錄組測序，將測序得到的大量數(shù)據(jù)進行拼接與組裝，并應用生物信息學方法對所得序列進行功能注釋、功能分類、代謝途徑分析和簡單重復序列(simple sequence repeat，SSR)位點查找等研究，從而在轉錄組水平系統(tǒng)研究女貞果實，并為進一步開展女貞果實的分子標記開發(fā)和重要次級代謝產(chǎn)物的生物合成奠定基礎。

1 材料與方法

1.1 女貞果實轉錄組文庫構建、測序和組裝

女貞果實于2015年8月采摘于湖南中醫(yī)藥大學湘杏學院校區(qū)內(nèi)的女貞樹上，用Trizol試劑分別提取3棵女貞樹成熟果實的總RNA，每個樣本取等量混合組成RNA池，然后用磁珠富集含poly(A)的mRNA并將其打斷成短序列，以這些短序列為模板用六堿基隨機引物依次合成第1條cDNA鏈和第2條cDNA鏈。cDNA鏈經(jīng)試劑盒純化并加EB緩沖液洗脫之后做末端修復、加poly(A)并連接測序接頭，然后用瓊脂糖凝膠電泳選擇200～500 bp片段進行PCR擴增，最后對建好的測序文庫用Illumina/Solexa Hiseq2000測序平臺進行測序。鑒于測序平臺數(shù)據(jù)錯誤率對結果的影響，對原始reads數(shù)據(jù)進行質(zhì)量預處理得到clean reads，即先使用cutadapt軟件去除3′端測序接頭[10]；再利用滑動窗口法去除低質(zhì)量片段：質(zhì)量閾值20(錯誤率=1%)，窗口大小10 bp，長度閾值35 bp；最后切除reads中含N部分序列：長度閾值35 bp。參照Grabherr 等提出方法[11]，將所得clean reads通過Trinity軟件(版本號trinityrnaseq_r2013-02-25)進行從頭拼接組裝，即通過clean reads之間的重疊信息組裝得到重疊群(Contigs)，然后局部組裝得到轉錄本(Transcripts)，最后用 Easycluster軟件對Transcripts進行同源聚類以得到單基因簇(Unigene)[12-13]。

1.2 Unigene功能注釋、GO分類和代謝路徑分析

通過Blastx在e值<0.000 01的條件下，將Unigene序列比對到蛋白數(shù)據(jù)庫Nr(Non-redundant protein database，非冗余蛋白數(shù)據(jù)庫)、Swiss-Prot(SwissProt protein database，蛋白質(zhì)序列數(shù)據(jù)庫)，并通過Blastn在e值<0.000 01的條件下，將Unigene比對到核酸數(shù)據(jù)庫Nt(Nucleotide collection，核酸數(shù)據(jù)庫)，得到與給定Unigene具有最高序列相似性且功能已知的蛋白，從而實現(xiàn)Unigene的功能注釋。根據(jù)Nr數(shù)據(jù)庫的功能注釋結果，使用Blast2GO軟件得到unigene的GO條目[14]，然后用WEGO軟件對所有的unigene進行GO功能分類統(tǒng)計[15]。KEGG 通路注釋是在evalue<0.000 01的條件下，利用blastx把unigene比對到KEGG數(shù)據(jù)庫，然后根據(jù)比對結果查詢與次級代謝產(chǎn)物合成通路相關的unigene。

1.3 SSR位點搜索和分析

利用MISA軟件在所有unigene 中搜索SSR位點，參數(shù)設置如下：二核苷酸至少重復次數(shù)為6，三核苷酸、四核苷酸、五核苷酸和六核苷酸至少重復次數(shù)均為4[16]。

2 結果與分析

2.1 測序結果評估和讀序列(reads)拼接組裝

采用Illumina Hiseq2000測序平臺，對女貞果實mRNA逆轉錄構建的cDNA文庫進行測序，共獲得32 856 614條clean reads，且本次測序總的cDNA堿基數(shù)約為4 928 492 100 bp(4.93 Gb)，故每條reads的平均長度為150 bp。Q20及GC含量百分比依次為97.94%和43.97%(見表1)。由以上結果表明，轉錄組測序數(shù)據(jù)量和質(zhì)量都較高，為后續(xù)的組裝提供了較好的原始數(shù)據(jù)。

利用Trinity軟件按paired-end方法對reads進行拼接組裝，獲得229 286個重疊群(contig)，平均長為322 bp，總長為73 933 206 bp，N50為475 bp(見表1)，其中小于200 bp 的重疊群數(shù)量占59.32%，200～3000 bp以及3000 bp以上各占40.16%和0.52%，具體長度分布見圖1。所得重疊群再次組裝共獲得163 092條unigene，平均長為665.68 bp，總長為108 567 136 bp(108.57 Mb)，N50為1345 bp(見表1)，其中100～500、500～2000、2000 bp以上的各占64.01%、28.67%、7.32%，具體長度分布見圖2。

表1 女貞果實轉錄組測序和組裝結果

圖1 Contig長度分布圖

圖2 Unigene長度分布圖

2.2 序列比對和功能注釋

將所獲得的unigene與公共數(shù)據(jù)庫Nr和Swiss-Prot進行Blastx比對，并利用Blastn比對Nt核酸數(shù)據(jù)庫，結果顯示共有83 903條unigene獲得了基因注釋，占unigene總數(shù)的51.45%。

2.3 功能分類和代謝途徑分析

2.3.1 GO和COG功能分類 為進一步揭示女貞果實中unigene的功能分類，組裝產(chǎn)生的163 092條unigene在GO(gene ontology)和COG(cluster of orthologous groups)數(shù)據(jù)庫中分別進行比對分析。在GO分類中，共有16 876 條unigene被注釋，其中12 203條unigene歸入生物學過程類別，3015條unigene歸入分子功能類別以及1658條unigene歸入細胞組分類別(見圖3)。在這3個類別中，生物學過程類別包括22個功能組，分子功能類別包括15個功能組，細胞組分類別則包括16個功能組(見圖3)。代謝過程是生物學過程類別中最大功能組(9468條)；細胞和細胞成分是細胞組分類別中最大功能組(都為5690條)；催化活性是分子功能類別中最大功能組(9519條)。通過以上GO分析，使我們在轉錄水平對女貞果實功能基因的分類情況有了一個系統(tǒng)的了解。在COG分類中可將unigene分為25個類別(見圖4)，其中unigene條數(shù)最多的10個類別依次是一般功能預測類(6875條)；復制、重組與修復類(3969條)；轉錄類(3733條)；信號轉導機制類(2926條)；蛋白質(zhì)翻譯后修飾、折疊和分子伴侶類(2632條)；碳水化合物運輸與代謝類(2198條)；氨基酸運輸與代謝類(1667條)；功能未知類(1371條)；細胞周期調(diào)控與分裂、染色體重排類(1201條)和無機離子運輸與代謝類(1197條)，而細胞外結構則是最小的類，僅包含4條unigene。通過以上COG分析，我們可以看到在女貞果實中，unigene涉及的功能類別主要是基因表達、蛋白質(zhì)合成和物質(zhì)代謝相關方面，這就正好符合果實作為一種營養(yǎng)器官，利用以上功能活動為自身積累營養(yǎng)物質(zhì)這一基本功能。

2.3.2 KEGG代謝途徑分析 KEGG(Kyoto encyclopedia of genes and genomes)可對細胞內(nèi)unigene參與的代謝途徑進行分析，本研究將獲得的unigene比對到此數(shù)據(jù)庫中，結果顯示女貞果實的unigene主要參與了21類代謝途徑(見圖5)，每條途徑中又包括下一級具體的代謝過程，其中涉及主要的次級代謝產(chǎn)物黃酮類合成途徑(Ko00941)含有258條unigene，萜類合成途徑(Ko00900，Ko00902，Ko00904，Ko00909)含有929條unigene。以上這些unigene首次在女貞代謝途徑中被發(fā)現(xiàn)，這為日后進一步驗證它們的功能、并利用它們在體外生物合成相關的藥用活性物質(zhì)奠定了基礎。

注：1.生物黏附；2.生物調(diào)節(jié)；3.細胞成分組織或生物發(fā)生；4.細胞進程；5.發(fā)育進程；6.定位的建立；7.生長；8.免疫系統(tǒng)進程；9.定位；10.運動；11.代謝進程；12.多有機體進程；13.多細胞有機體進程；14.生物進程的負調(diào)控；15.生物進程的正調(diào)控；16.生物進程調(diào)控；17.復制；18.生殖進程；19.刺激應答；20.節(jié)律進程；21.信號；22.單有機體進程；23.細胞；24.細胞連接；25.細胞成分；26.細胞外基質(zhì)；27.細胞外基質(zhì)成分；28.胞外區(qū)域；29.大分子復合物；30.膜；31.膜成分；32.膜關閉內(nèi)腔；33.核區(qū)；34.細胞器；35.細胞器成分；36.共質(zhì)體；37.病毒體；38.病毒體成分；39.抗氧化活性；40.結合；41.催化活性；42.電子載體活性；43.酶調(diào)節(jié)活性；44.鳥苷酸交換因子活性；45.金屬伴侶活性；46.分子轉導活性；47.核酸結合的轉錄因子活性；48.營養(yǎng)庫活性；49.蛋白結合的轉錄因子活性；50.蛋白尾；51.受體活性；52.結構分子活性；53.轉運活性。圖3 女貞果實unigene的GO分類

注：1.翻譯，核糖體結構與生物發(fā)生；2.轉錄；3.信號傳導機制；4.次級產(chǎn)物合成，運輸及代謝；5.加工與修飾；6.復制、重組與修復；7.蛋白質(zhì)翻譯后修飾與轉運，分子伴侶；8.核苷酸運輸與代謝；9.核結構；10.脂類轉運與代謝；11.細胞內(nèi)轉運，分泌和小泡運輸；12.無機離子轉運和代謝；13.一般功能預測；14.功能未知；15.胞外結構；16.能量生成和轉換；17.防御機制；18.細胞骨架；19.輔酶運輸和代謝；20.染色質(zhì)結構與變化；21.細胞壁/膜生物發(fā)生；22.細胞運動；23.細胞周期控制，細胞分裂，染色體區(qū)分；24.碳水化合物運輸和代謝；25.氨基酸運輸與代謝。圖4 女貞果實unigene的COG分類

注：1.運輸與代謝；2.膜轉運；3.信號轉導；4.折疊，分類和降解；5.復制和修復；6.轉錄；7.翻譯；8.藥物抵抗；9.內(nèi)分泌代謝??；10.氨基酸代謝；11.其他次生物質(zhì)代謝；12.碳水化合物代謝；13.能量代謝；14.總覽圖；15.糖生物合成和代謝；16.脂質(zhì)代謝；17.輔助因子和維生素代謝；18.其他的氨基酸代謝；19.萜類和聚酮類化合物的代謝；20.核苷酸代謝；21.環(huán)境適應。圖5 女貞果實unigene的KEGG分類

2.4 SSR分析

利用SSR分析軟件，從14 270條unigene中共搜索到15 925個SSR位點。搜索到的SSR位點類型豐富，單核苷酸至六核苷酸類型均有(見表2)。其中，二核苷酸重復所占比例最高，達到了43.95%；比例最低的是四核苷酸重復，僅為2.37%。在檢測到的SSR中，出現(xiàn)頻率最高的10類基序為：AG/GT(3856個)、AT/AT(183個)、AC/GT(125個)、AAT/ATT(940個)、AAG/CTT(736個)、ATC/ATG(501個)、ACC/GGT(491個)、AGC/CTG(316個)、AGG/CCT(312個)、CCG/CGG(133個)。

表2 女貞果實unigene中的SSR位點數(shù)量與分布

3 討論

中藥材是人類醫(yī)藥衛(wèi)生領域的一個瑰寶，如何開發(fā)利用它使其不斷的造福當代世人，一直是現(xiàn)今許多中醫(yī)藥科研工作者不斷追求的目標。隨著生物醫(yī)藥科研技術的不斷發(fā)展，許多新技術新方法應用于中藥材研究，使得以上目標逐步被實現(xiàn)，其中基于高通量測序技術的轉錄組學研究便是一個典型例子[17-18]。目前，轉錄組學在中藥材研究領域主要集中于藥用植物功能基因組的研究，即通過轉錄組學發(fā)現(xiàn)藥用植物次生代謝產(chǎn)物生物合成關鍵酶基因、闡明次生代謝途徑及其調(diào)控機制等，這將為利用功能基因和代謝途徑體外生物合成中藥材中藥用活性成分，甚至是為創(chuàng)造性的研發(fā)新型藥物奠定堅實的基礎[19-20]。本研究正是基于以上描述，采用Illumina/Solexa高通量測序技術對女貞果實這一中藥材進行了轉錄組測序，所得結果顯示共有4.93 Gb有效數(shù)據(jù)產(chǎn)生，可組裝拼接得到163 092個unigene，其中有79 189條unigene無法獲得功能注釋，究其原因：一是因為女貞缺乏基因組、EST和蛋白質(zhì)序列信息，不能進行基因組條件下的比對分析，故導致較多的unigene無法匹配到合適的同源序列上；二是由于不能注釋的unigene多為小于1000 bp的小片段，因而它們很難與公共數(shù)據(jù)庫中已知的序列獲得良好的比對；三是小片段的unigene有些可能是短的非編碼序列或者新基因的某些不完整片段，從而導致它們無法被已有數(shù)據(jù)庫所注釋。另外，由測序數(shù)據(jù)組裝拼接成unigene的過程中，N50值(指從組裝最長的unigene依次向下求長度的總加和，當累加長度達到組裝長度的一半時對應的unigene長度)是一個判斷組裝程度好壞的重要指標，其值越大說明組裝得到的長片段越多，組裝效果越好。本次研究unigene的N50值為1345 bp，大于1000 bp，說明本次序列組裝的質(zhì)量和長度可以滿足轉錄組分析的基本要求，這為研究中進行以序列為基礎的相關轉錄組分析提供了可靠保障。

轉錄組研究的一個重要目標是發(fā)現(xiàn)功能基因，并闡明功能基因編碼產(chǎn)物所處的代謝途徑。要實現(xiàn)這一目標，研究者需利用GO數(shù)據(jù)庫、COG數(shù)據(jù)庫和KEGG數(shù)據(jù)庫對組裝得到的unigene進行分析。本研究中筆者通過GO數(shù)據(jù)庫的分析，首次系統(tǒng)地闡明了女貞果實功能基因的分類情況，并利用COG數(shù)據(jù)庫對女貞果實unigene進行了功能注釋，使其從分子水平找尋自己的直系同源體，從而預測它們的生物學功能，最后經(jīng)KEGG數(shù)據(jù)庫把一些預測的功能基因定位在了涉及主要次級代謝產(chǎn)物黃酮類和萜類合成的途徑中。在此基礎上，為進一步驗證這些預測基因的真正功能，筆者將在后續(xù)的研究中對它們進行體外克隆和異源表達實驗，并在體外和體內(nèi)兩個方面對表達產(chǎn)物進行功能活性檢測，從而用實驗結果檢驗轉錄組分析結果的正確與否。

SSR是真核生物基因組非編碼區(qū)中存在的一些簡單重復序列，在物種進化過程中呈現(xiàn)較為保守的特性，已廣泛應用于中藥材物種、種質(zhì)資源以及農(nóng)家種鑒定[21]。目前，女貞缺乏可使用的分子標記，此次果實轉錄組研究所獲大量數(shù)據(jù)，可為開發(fā)女貞SSR標記提供豐富的篩選資源。在今后的研究中，可利用已獲得的15 925個SSR位點，通過引物設計軟件設計合理的SSR引物，然后利用SSR引物進行擴增檢測，從中篩選出擴增穩(wěn)定、條帶清晰、多態(tài)性好的引物，為進一步開發(fā)女貞的SSR標記奠定基礎[22]。

本研究利用Illumina/Solexa Hiseq2000測序平臺首次對女貞果實進行了轉錄組測序研究，獲得了大量的轉錄組數(shù)據(jù)。對這些數(shù)據(jù)從拼接組裝、功能注釋、代謝途徑和SSR位點查詢4個方面進行分析研究，使我們在轉錄水平對女貞果實有了比較詳實的認識，同時也為今后開發(fā)利用女貞果實的藥用次生代謝產(chǎn)物和分子標記物，以及為女貞基因組的測序組裝提供了具有參考價值的基礎數(shù)據(jù)。

[1] 國家藥典委員會.中華人民共和國藥典：一部[S].北京：中國醫(yī)藥科技出版社，2015：45-46.

[2] 尹輝.中藥女貞子化學成分研究綜述[J].九江學院學報(自然科學版)，2015，30(1)：74-75.

[3] 劉亭亭，王萌.女貞子化學成分與藥理作用研究進展[J].中國實驗方劑學雜志，2014，20(14)：228-234.

[4] Nagalakshmi U，Waern K，Snyder M.RNA-Seq：A Method for Comprehensive Transcriptome Analysis[J].Curr Protoc Mol Biol，2010(suppl 89):Unit 4.11.1-13.

[5] Hoeijmakers W A，Bártfai R，Stunnenberg H G.Transcriptome analysis using RNA-Seq[J].Methods Mol Biol，2013，923：221-239.

[6] 田英芳，張曉政，周錦龍.轉錄組學研究進展及應用[J].中學生物教學，2013(12)：29-31.

[7] Martin J A，Wang Z.Next-generation transcriptome assembly[J].Nature Reviews Genetics，2011，12(10)：671-682.

[8] Chen C，Huang H，Wu C H.Protein Bioinformatics Databases and Resources[J].Methods Mol Biol，2011，694:3-39.

[9] Collins L J，Biggs P J，Voelckel C，et al.An approach to transcriptome analysis of non-model organisms using short-read sequences[J].Genome Informatics，2008，21：3-14.

[10] Martin M.Cutadapt removes adapter sequences from high-throughput sequencing reads[J].Embnet Journal，2011，17(1)：10-12.

[11] Grabherr M G，Haas B J，Yassour M，et al.Full-length transcriptome assembly from RNA-Seq data without a reference genome[J].Nat Biotechnol，2011，29：644-652.

[12] Picardi E，Mignone F，Pesole G.EasyCluster：a fast and efficient gene-oriented clustering tool for large-scale transcriptome data[J].BMC Bioinformatics，2009，10(6)：1-12.

[13] Bevilacqua V，Pietroleonardo N，Giannino E，et al.EasyCluster2：an improved tool for clustering and assembling long transcriptome reads[J].BMC Bioinformatics，2014，15(15)：1-10.

[14] Conesa A，G?tz S，García-Gómez J M，et al.Blast2GO：A universal tool for annotation，visualization and analysis in functional genomics research[J].Bioinformatics，2005，21：3674-3676.

[15] Ye J，F(xiàn)ang L，Zheng H，et al.WEGO：A web tool for plotting GO annotations[J].Nucleic Acids Res，2006，34：293-297.

[16] Wang H，Jiang J，Chen S，et al.Next-Generation Sequencing of the Chrysanthemum nankingense(Asteraceae)Transcriptome Permits Large-Scale Unigene Assembly and SSR Marker Discovery[J].PloS One，2013，8(4)：e62293.

[17] 張曉萌，李健春，王瓊，等.轉錄組測序技術在中醫(yī)藥領域的應用[J].中國現(xiàn)代中藥，2016，18(8)：1084-1087.

[18] 陳惠，辛麗麗，龔婕寧.基于全轉錄組測序技術的轉錄組學在中醫(yī)藥領域的應用前景分析[J].環(huán)球中醫(yī)藥，2013，6(10)：759-763.

[19] 王堯龍，黃璐琦，袁媛，等.藥用植物轉錄組研究進展[J].中國中藥雜志，2015，40(11)：2055-2061.

[20] 吳瓊，孫超，陳士林，等.轉錄組學在藥用植物研究中的應用[J].世界科學技術—中醫(yī)藥現(xiàn)代化，2010，12(3)：457-462.

[21] 周駿輝，袁媛，黃璐琦.SSR標記在中藥材分子身份證體系構建中的應用[J].中國現(xiàn)代中藥，2016，18(10)：1233-1236.

[22] 劉峰，王運生，田雪亮，等.辣椒轉錄組SSR挖掘及其多態(tài)性分析[J].園藝學報，2012，39(1)：168-174.