基于高通量測序的野生毛葡萄轉錄組數(shù)據(jù)分析

方輝++蔣勝理++曲俊杰++周思泓

摘要：利用RNA-seq技術對所構建的野生毛葡萄（Vitis quinquangularis Rehd）葉片的轉錄組進行測定，對原始reads進行過濾和組裝，得到了35 238條質量較高的unigene，平均長度為1 081 nt，N50為1 735 nt?；贜CBI蛋白質數(shù)據(jù)庫（Nr）、蛋白質序列數(shù)據(jù)庫（Swiss-Prot）、京都基因與基因組百科全書（KEGG）數(shù)據(jù)庫和直系同源基因簇（COG）進行相似性比對，共注釋了26 751條unigene，另有8 487條unigene未被注釋。物種同源性顯示與葡萄的同源性最高為74.48%。利用COG數(shù)據(jù)庫將unigene分成25類，通過GO分類和KEGG富集性分析，將unigene分別歸類于44個GO類別和122個代謝途徑。此外，在35 238條unigene中共搜索到4 428個SSR位點，二核苷酸的SSR數(shù)目最多（1 906條），其次為三核苷酸（1 762條）。這些信息為毛葡萄功能基因、相關候選基因的發(fā)掘以及分子標記輔助育種提供了重要依據(jù)。

關鍵詞：毛葡萄（Vitis quinquangularis Rehd）；轉錄組；高通量測序

中圖分類號：S663.1 文獻標識碼：A 文章編號：0439-8114（2016）18-4827-05

DOI：10.14088/j.cnki.issn0439-8114.2016.18.050

野生毛葡萄（Vitis quinquangularis Rehd）是葡萄科（Vitaceae）葡萄屬（Vitis）真葡萄亞屬（Euvitis）東亞種群，盛產(chǎn)于中國，是中國葡萄屬東亞種群中分布最廣的一個野生種[1]。廣西壯族自治區(qū)作為中國野生毛葡萄之鄉(xiāng)和中國野生毛葡萄酒原產(chǎn)地，在桂中、桂北、桂西的喀斯特地貌地區(qū)野生毛葡萄較集中。野生毛葡萄含多種維生素、各種氨基酸、礦物質等營養(yǎng)物質，還含有豐富的花色素雙糖苷、超氧化歧化酶（SOD）及白藜蘆醇（RES）物質，具有軟化血管、抗衰老、防癌的作用[2]。多數(shù)野生毛葡萄表現(xiàn)抗逆性強、耐病性好、耐南方高溫多雨潮濕天氣、耐旱耐瘠等優(yōu)良特性，是南方石漠化治理、生態(tài)重建的優(yōu)良樹種，具有重要的經(jīng)濟價值[3，4]。目前對毛葡萄的研究主要集中在栽培技術和雜交育種方面，缺乏對其分子生物學方面的研究，對其基因信息的研究極少。

轉錄組是在特定發(fā)育階段和一定生理條件下，細胞內所有轉錄出來的RNA的集合。通過高通量的轉錄組分析，可以獲得機體在生命過程中基因的表達模式[5]。RNA-seq是2008年建立起來的基于深度測序的轉錄組分析技術，它能夠在單核苷酸水平上對任何物種進行整體轉錄活動的檢測，是當前在全基因組水平上研究基因表達模式的主要技術[6]。本研究擬采用RNA-seq技術對所構建的毛葡萄葉片的轉錄組進行測定，并在基因組水平上進行轉錄組分析，旨在獲得毛葡萄更多的轉錄本和更為全面的轉錄組信息，為重要性狀相關基因的克隆及功能分析、鑒定次級代謝物生物合成相關基因奠定基礎。

1 材料與方法

1.1 材料

植物材料為毛葡萄五峒-1，生長于廣西壯族自治區(qū)河池市都安瑤族自治縣高嶺鎮(zhèn)五峒村，樹齡數(shù)十年、生長旺盛，葉片生長狀況良好，果實產(chǎn)量高，經(jīng)調查與試驗發(fā)現(xiàn)其對多種病害高抗或免疫。采集后液氮速凍，并轉移至-80 ℃水箱保存?zhèn)溆谩?/p>

1.2 方法

1.2.1 野生毛葡萄測序 采取野生毛葡萄的嫩葉并提取總RNA，送至北京諾禾致源生物信息科技有限公司進行文庫構建與高通量測序，測序平臺為Illumina HiSeq TM 2000，每個樣品測序量約為4.75 G。

1.2.2 測序數(shù)據(jù)的組裝 對測序后得到的原始數(shù)據(jù)total raw reads進行質量分析，去除重復、含接頭、測序質量低的reads，獲得clean reads。使用短reads組裝軟件Trinity做轉錄組從頭組裝[7]。首先將具有一定長度重疊的reads連成更長的片段，通過reads重疊關系得到組裝片段contig。然后將reads比對回contig，通過paired-end reads確定來自同一轉錄本的不同contig以及這些contig間的距離，Trinity軟件將這些contig連在一起，最后得到兩端不能再延長的序列，即為unigene。

1.2.3 功能注釋 利用Blastx將unigene序列與Nr（Non-redundant Protein Sequence Database in GenBank）、Swiss-Prot（Swiss-Prot Protein Sequence Database）、KEGG（Kyoto Encyclopedia of Genes and Genomes）和COG（Cluster of Orthologous Groups of proteins）數(shù)據(jù)庫進行比對，獲取與unigene具有最高序列相似性的蛋白，從而得到該unigene的蛋白功能注釋信息。根據(jù)Nr注釋信息，使用Blast2GO[8]軟件進行GO注釋，得到每個unigene的GO信息后，用WEGO[9]軟件進行GO功能分類統(tǒng)計。

按照Nr、Swiss-Prot、KEGG、COG的優(yōu)先級順序，將unigene序列與以上蛋白庫進行Blastx比對，取比對結果中等級最高的蛋白確定該unigene的編碼區(qū)序列，然后根據(jù)標準密碼子表將編碼區(qū)序列翻譯成氨基酸序列，從而得到該unigene編碼區(qū)的核酸序列（序列方向5′→3′）和氨基酸序列。最后，將與以上蛋白庫都比對不上的unigene再用軟件ESTScan預測其編碼區(qū)，得到其編碼區(qū)的核酸序列（序列方向5′→3′）和氨基酸序列[10]。

1.2.4 SSR位點的篩選 利用MIcroSAtellite（MISA） 軟件在所有unigene中搜索SSR位點，設置二核苷酸重復次數(shù)至少為6次，三核苷酸、四核苷酸、五核苷酸和六核苷酸重復次數(shù)至少為4次。

2 結果與分析

2.1 轉錄組測序和組裝

通過Illumina Hiseq2000平臺測序，總計產(chǎn)出56 374 280條reads，去除低質量和含有接頭的reads后，得到53 697 490條clean reads，共計5 369 749 000個核苷酸（Nucleotides，nt），GC含量45.52%、Q20為95.63%。利用Trinity軟件對這些reads進行組裝得到35 238條unigene，平均長度為1 081 nt，N50為 1 735 nt。unigene的長度分布（圖1）顯示，長度大于1 000 nt的unigene 13 651條，占全部unigene的38.74%。這些說明本研究中轉錄組文庫的測序和組裝結果都較好，能夠進行后續(xù)生物信息學分析。

2.2 功能注釋

為了預測unigene功能，分別將unigene與Nr、Swiss-Prot、COG、KEGG等生物學數(shù)據(jù)庫進行比對。通過Blast比對（表1），共有26 751條unigene獲得了基因注釋，占所有unigene的75.92%；有8 487條unigene（24.08%）未被注釋，這些未被注釋的基因可能認為是新基因。

利用blastx將組裝出來的unigene序列與Nr數(shù)據(jù)庫進行比對后，取每個unigene在Nr庫中比對結果最好（E值最低）的那一條序列為對應同源序列（如有并列，取第一條）確定同源序列所屬物種，統(tǒng)計比對到各個物種的同源序列數(shù)量，共找到26 719個unigene與其他近緣生物的已知基因具有不同程度的同源性，其中有19 901條（74.48%）與葡萄（Vitis vinifere）同源；1 371條（5.13%）與可可（Theobroma cacao）序列同源；865條（3.24%）與樹棉（Gossypium arboreum）序列同源。35 238條unigene中還有8 519條（24.18%）unigene沒有比對到近源物種 （圖2）。

將所有unigene基因與COG和GO數(shù)據(jù)庫進行比對，并進行功能注釋與分類。26 816條unigene被注釋到25個COG類別中（圖3），其中一般功能基因（Generalfunction prediction only）是最大類別，包含5 280條unigene，占被注釋unigene總數(shù)的19.69%； 其次是信號傳導機制（Signal transduction mechanisms），包含3 293條（12.28%）unigene，其余依次為翻譯后修飾，蛋白質周轉以及分子伴侶類（Posttranslational modification，protein turnover，chaperones）（2 980條，11.11%），RNA加工與修改（RNA processing and modification）（1 591條，5.93%），而細胞運動（Cell motility）是最小的類別，僅包含10條unigene。

在已經(jīng)得到的Nr注釋信息基礎上，采用Blast2GO獲得毛葡萄unigene的GO分類信息，共有14 692條unigene得到GO注釋（圖4）。在GO分類體系中，生物學過程（biological process）、細胞組分（cellular component）和分子功能（Molecular function）。這3個大的類別又被劃分為詳細的44個小類別，其中“代謝過程（Metabolic process）”（8 819， 60.03%）、“細胞過程（Cellular process）”（8 352， 56.84%）和“結合（Binding）”（8 027，54.64%），3個類群占了主要部分，而“翻譯調節(jié)活性（Translation regulator activity）”（2，0.01%）、“核酸結合轉錄因素活性（Nucleic acid binding transcription factor activity）”（3，0.02%）和“內核（nucleoid）”（4，0.04%）僅有非常少的基因歸入，這一分類結果顯示了毛葡萄基因表達譜的總體情況。

2.3 代謝通路分析

對野生毛葡萄葉片的轉錄組進行KEGG分析發(fā)現(xiàn)，有6 871條unigene注釋到KEGG數(shù)據(jù)庫中，分布于122條已知的通路中，包括淀粉和蔗糖代謝（Starch and sucrose metabolism）（165條，ko00500）、植物病原物相互作用（Plant-pathogen interaction）（145條，ko04626）和光合生物碳固定（Carbon fixation in photosynthetic organisms）（78條，ko00710）。注釋序列數(shù)目較多的5個通路分別是代謝途徑（Metabolic pathways）（1 743條，ko01100）、次生代謝物生物合成（Biosynthesis of secondary metabolites）（812條，ko01110）、核糖體（Ribosome）（278條，ko03010）、植物激素信號轉導（Plant hormone signal transduction）（220條，ko04075）和真核細胞核糖體合成（Ribosome biogenesis in eukaryotes）（209，ko03008）（表2）。

2.4 SSR分析

利用MISA軟件在毛葡萄的35 238條unigene中共搜索到4 428個SSR位點，占unigene總序列數(shù)的12.57%，平均每8.6 kb出現(xiàn)1個SSR，其中包含有兩個及兩個以上SSR的unigene共有576條。二核苷酸和三核苷酸重復類型的SSR數(shù)目最多，其中二核苷酸1 906條、三核苷酸1 762條；而五核苷酸和六核苷酸重復類型在毛葡萄中轉錄組序列中含量較少，僅有165條和201條（表3）。除此之外，不同核苷酸的重復次數(shù)也有很大的變化。在檢測到的SSR中出現(xiàn)頻率最高的5類基序為AG/CT 1 425個、AAG/CTT 552個、AT/TA 345個、ATC/ATG 247個、AGC/CTG 235個，上述SSR特征分析有助于開展毛葡萄及其同屬物種的基因組差異分析通用性標記開發(fā)和遺傳圖譜構建的研究。

3 討論

基于高通量測序技術的轉錄組學研究是一種高效、可靠的發(fā)掘功能基因的手段，Illumina高通量測序的數(shù)據(jù)量大、速度快、成本低、效率高[11]。本研究應用Illumina高通量測序技術對毛葡萄轉錄組進行測序，對毛葡萄轉錄組測序獲得了56 374 280個reads序列，去除低質量的和含有接頭的reads以后，得到53 697 490條clean reads，5 369 749 000個核苷酸，GC含量45.52%、Q20為95.63%。利用Trinity軟件對這些reads進行組裝得到35 238條unigene，長度從224 nt到15 540 nt，平均長度為1 081 nt，N50 為1 735 nt，N50值越大反映組裝得到的長片段越多組裝效果就越好。測序數(shù)據(jù)產(chǎn)量和數(shù)據(jù)組裝質量是評價轉錄組測序完成情況的重要指標。以上研究結果表明，此次序列組裝的質量和長度可以滿足轉錄組分析的基本要求。

結合生物信息學分析方法對毛葡萄unigene與Nr、Swiss-Prot、COG、KEGG數(shù)據(jù)庫進行比對，進行序列相似性和功能注釋分析。26 719個unigene與其他近緣生物的已知基因具有不同程度的同源性，有19 901條（74.48%）與葡萄同源，與葡萄的親緣關系最近。KEGG分析中，19 274條unigene被注釋到128個代謝通路。根據(jù)KEGG數(shù)據(jù)庫對上述unigene進行代謝途徑分析，涉及112個具體的代謝途徑分支，參與到毛葡萄體內的碳水化合物代謝、脂類代謝、次生物質代謝等過程中，為進一步大量挖掘毛葡萄生長發(fā)育過程中的重要表達基因，開展毛葡萄的基因克隆及功能驗證等研究提供了基礎數(shù)據(jù)。轉錄組數(shù)據(jù)分析結果表明，毛葡萄葉片為植物生長提供了物質、能量、信號轉導以及抗病條件。

SSR分子標記具有操作簡便、重復性好、多態(tài)性豐富、遺傳信息量大、共顯性遺傳等優(yōu)點，已在遺傳多樣性分析、遺傳圖譜構建、功能基因發(fā)掘、分子標記輔助育種等研究中得到了廣泛應用。采取實驗室方法開發(fā)SSR標記費時、耗力、成本高、試驗復雜，基于轉錄組數(shù)據(jù)信息進行SSR分子標記開發(fā)是一種經(jīng)濟有效的方法。目前毛葡萄可利用的分子標記數(shù)量非常有限，轉錄組產(chǎn)生的海量數(shù)據(jù)為SSR分子標記的開發(fā)提供了豐富和極有價值的可利用資源。本研究在35 238條unigene中共搜索到4 428個SSR位點，占unigene總序列數(shù)的12.57%。二核苷酸和三核苷酸重復類型的SSR數(shù)目最多，其中二核苷酸1 906條（43.04%）、三核苷酸1 762條（39.79%），前人的研究結果顯示有些物種二核苷酸占優(yōu)勢，另一些物種三核苷酸占優(yōu)勢，這種現(xiàn)象的出現(xiàn)與MISA搜索時相關參數(shù)的設置有很大的關系[12]，后續(xù)可對這些SSR進行引物設計和擴增檢測，篩選出擴增穩(wěn)定、條帶清晰、多態(tài)性好的引物為進一步開發(fā)新的SSR標記奠定了基礎。

本研究采用Illumina Hiseq2000高通量測序技術建立了毛葡萄轉錄組數(shù)據(jù)庫，獲得了大量的轉錄本信息，并對表達基因進行了序列組裝、功能注釋、代謝途徑等分析，為今后更深入研究毛葡萄功能基因組、基因克隆及抗逆機理研究提供了極大的方便。

參考文獻：

[1] 林 玲，張 瑛，盧 江，等.中國葡萄屬野生資源毛葡萄的研究與利用進展[J].中外葡萄與葡萄酒，2015（2）：66-69.

[2] 涂 佳，鄧學良，周文化，等.野生毛葡萄常規(guī)營養(yǎng)成分和原花青素含量分析[J].食品與機械，2010（2）：82-85.

[3] 鄒 瑜，吳代東，張 瑛，等.毛葡萄霜霉病防治試驗[J].廣西農(nóng)業(yè)科學，2002（6）：315-317.

[4] 鄒 瑜. 野生毛葡萄種質資源及利用研究進展[J].南方農(nóng)業(yè)學報，2008，39（5）：664-667.

[5] WANG Z，GERSTEIN M，SNYDER M. RNA-Seq：A revolutionary tool for transcriptomics[J].Nature Reviews Genetics，2009， 10（1）：57-63.

[6] COSTA V，ANGELINI C，DE FEIS I，et al. Uncovering the complexity of transcriptomes with RNA-Seq[J]. Journal of Biomedicine and Biotechnology，2010：853916.

[7] HAAS B J，PAPANICOLAOU A，YASSOUR M，et al. De novo transcript sequencere construction from RNA-seq using the Trinity platform for reference generation and analysis[J].Nature Protocols，2013，8（8）：1494-1512.

[8] 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（18）：3674-3676.

[9] YE J，F(xiàn)ANG L，ZHENG H，et al. WEGO：A web tool for plotting GO annotations[J].Nucleic Acids Research，2006，34（1）：293-297.

[10] BIAO J，DASEN X，WENRUI L，et al. De Novo assembly and characterization of the transcriptome， and development of SSR markers in wax gourd（Benicasa hispida）[J].Plos One，2013， 8（8）：653-661.

[11] SUN C， LI Y，WU Q，et al. De novo sequencing and analysis of the American ginseng root transcriptome using a GS FLX Titanium platform to discover putative genes involved in ginsenoside biosynthesis[J]. BMC Genomics，2010，11：262.

[12] WEI W L，QI X Q，WANG L H，et al. Characterization of the sesame（Sesamum indicum L.） global transcriptome using Illumina paired-end sequencing and development of EST-SSR markers[J].BMC Genomics，2011，12：451.