芒果果實轉(zhuǎn)錄組數(shù)據(jù)組裝及基因功能注釋

武紅霞　許文天　羅純　姚全勝　王松標　馬小衛(wèi)　詹儒林

摘 要 以紅芒6號為試材，對生長發(fā)育各時期果皮與果肉樣品提取RNA后等量混合進行轉(zhuǎn)錄組測序，共獲得了68 419 722個reads，包含6 157 744 980個核苷酸序列信息，平均讀長90 bp，序列信息全部登陸NCBI數(shù)據(jù)庫（SRP035450）。對reads進行拼接，共獲得124 002個Contig片段，進一步組裝獲得54 207個Unigene片段，平均長度為838 bp，閾值小于10-5的序列中共有42 515（78.43%）個unigenes被注釋到公共蛋白數(shù)據(jù)庫。其中，35 198個被注釋到生物學過程、分子功能和細胞組分GO類別的44個功能組，14 619個Unigene匹配到25類COG功能組，23 741個Unigene被富集到128個KEGG代謝通路，包括代謝途徑、次生代謝的生物合成、植物-病原物互作、植物激素信號轉(zhuǎn)導、苯丙氨酸生物合成、淀粉和蔗糖代謝、類黃酮類化合物合成等。本研究將為進一步的芒果功能基因克隆、基因表達分析、分子標記開發(fā)等研究奠定基礎(chǔ)。

關(guān)鍵詞 芒果；轉(zhuǎn)錄組；功能注釋

中圖分類號 S682.2 文獻標識碼 A

芒果是著名的熱帶水果，素有“熱帶果王”之美譽。中國是主要的芒果生產(chǎn)國之一，2014年栽培面積14萬hm2，產(chǎn)量130萬t，主要分布在海南、廣西、云南、四川、廣東和福建等省（區(qū)）。前期研究主要集中于栽培生理研究，如果實成熟過程[1-4]、果實香氣揮發(fā)物質(zhì)[5-6]、抗氧化活性[7-8]、采后處理和果實品質(zhì)研究[9-10]，盡管芒果是中國熱區(qū)農(nóng)業(yè)產(chǎn)業(yè)的經(jīng)濟支柱，其基因組較小，但基因組信息相當缺乏。近年來，高通量轉(zhuǎn)錄組測序技術(shù)的不斷發(fā)展和完善為開展不同生物功能基因組學研究提供了全新的思路和方法，轉(zhuǎn)錄組測序是基于下一代高通量測序技術(shù)建立的一種高效、快捷分子生物學研究手段，使科技工作者能在組學水平上研究基因組序列未知的非模式生物，從整體水平了解植物在特定階段的基因功能及基因結(jié)構(gòu)，更加便利地提示特定生物學過程的分子機制[11-13]。近年來，已廣泛應用于楊梅[14]、蘋果[15]、梨[16]、葡萄[17]、柑橘[18]等果樹，芒果研究方面，Dautt-Castro等[19]用轉(zhuǎn)錄組測序技術(shù)研究了與 Kent芒成熟過程相關(guān)的基因，在青熟和成熟果實中，2 306個基因差異表達。Luria 等[20]用轉(zhuǎn)錄組測序揭示了熱水處理對芒果影響的分子機制，基于基因表達譜測序技術(shù)鑒定了與熱處理相關(guān)基因，Azim等[21]報道了芒果葉綠體基因組，這些研究豐富了芒果生物信息數(shù)據(jù)。有關(guān)芒果果實生長發(fā)育過程的轉(zhuǎn)錄組信息鮮有報道。本研究利用高通量測序技術(shù)平臺illumina Hiseq 2 000對芒果發(fā)育成熟過程的果實進行轉(zhuǎn)錄組測序、拼接組裝，再用生物信息學的方法對得到的Unigene進行注釋和功能分類，以期為芒果功能基因的發(fā)掘利用、特異miRNA的鑒定及功能分析等奠定基礎(chǔ)。

1 材料與方法

1.1 材料

試驗于2012年3～7月在廣東省雷州市覃斗鎮(zhèn)芒果基地進行，試驗材料為15年生的紅芒6號芒果樹，試驗地土壤屬磚紅壤，株行距4 m×4 m，樹體生長良好。分別于幼果期（花后50 d）、膨大期（花后80 d）、采收期（青熟）與成熟期（常溫放置后成熟）采集芒果果實，每個時期采10～15個果，采集后迅速帶回實驗室，將果皮與果肉分開并用液氮速凍后保存于-80 ℃冰箱中備用。

1.2 方法

1.2.1 芒果果皮與果肉RNA提取 RNA提取參照Shan等[22]建立的方法。

1.2.2 芒果果實轉(zhuǎn)錄組測序 基于轉(zhuǎn)錄組學研究的優(yōu)勢，以紅芒6號芒果為材料，取不同發(fā)育階段的果皮與果肉樣本，分別提取總RNA后等量混合用于測序，轉(zhuǎn)錄組測序委托深圳華大基因公司完成。提取的總RNA使用DNaseI消化DNA后，用帶有Oligo（dT）的磁珠富集真核生物mRNA；加入打斷試劑在Thermomixer中適溫將mRNA打斷成短片段，以打斷后的mRNA為模板合成一鏈cDNA，然后配制二鏈合成反應體系合成二鏈cDNA，并使用試劑盒純化回收、粘性末端修復、cDNA的3′末端加上堿基“A”并連接接頭，然后進行片段大小選擇，最后進行PCR擴增；構(gòu)建好的文庫用Agilent 2 100 Bioanalyzer和ABI StepOnePlus Real-Time PCR System質(zhì)檢合格后，使用Illumina HiSeqTM 2000進行測序。

1.2.3 數(shù)據(jù)組裝及基因功能注釋 將測序獲得的原始數(shù)據(jù)去除接頭序列、或N含量過高或低質(zhì)量的序列，得到凈讀長，用Trinity軟件按Min contig length 100，Group pairs distance 250，path reinforcement distance 85，Min kmer cov 2參數(shù)進行數(shù)據(jù)組裝。將具有重疊區(qū)域的reads連成更長的Contig，采用Overlap的方法進一步拼接成Unigene。

通過Blastx程序?qū)nigene序列比對到蛋白數(shù)據(jù)庫NR、Swiss-Prot、KEGG和COG（E值<0.00 001），并通過Blastn程序?qū)nigene與核酸數(shù)據(jù)庫NT（E值<0.00 001）進行比對，獲得Unigene最高序列相似性的蛋白，便得到該Unigene的蛋白功能注釋信息。將Unigene與COG數(shù)據(jù)庫進行比對，獲得Unigene可能的功能注釋及功能分類。根據(jù)nr注釋信息，利用Blast2GO軟件得到Unigene的GO注釋信息，并對所有Unigene用WEGO軟件進行GO功能分類統(tǒng)計，根據(jù)KEGG注釋信息進一步分析得到Unigene的Pathway途徑注釋。

預測編碼蛋白框：按NR、Swiss-Prot、KEGG和COG順序?qū)nigene序列與以上數(shù)據(jù)庫Blastx比對（E值<10-5），取Blast比對結(jié)果中排列最高的蛋白確定該Unigene的編碼區(qū)序列，將編碼區(qū)序列翻譯成氨基酸序列，得到該基因編碼區(qū)的核苷酸序列（序列方向5′→3′）和氨基酸序列。若Unigene比對不上以上數(shù)據(jù)庫，用軟件ESTScan[23]預測編碼區(qū)，得到其編碼區(qū)的核苷酸序列和氨基酸序列。

2 結(jié)果與分析

2.1 轉(zhuǎn)錄組組裝結(jié)果

2.1.1 果實轉(zhuǎn)錄組組裝概要 轉(zhuǎn)錄組測序共獲得82 817 950讀段，經(jīng)過去除雜質(zhì)和冗余處理，得到68 419 722條凈讀段，6.1Gb（6 157 744 980 bp）的有效數(shù)據(jù)，13×的測序深度，平均讀長90 bp。其中讀長大于20個堿基比例為95.27%，GC%值為45.52%，可以看出此次轉(zhuǎn)錄組測序結(jié)果較好，可為后續(xù)的數(shù)據(jù)組裝提供很好的原始數(shù)據(jù)。轉(zhuǎn)錄組測序組裝得到124 002個Contig，平均長度為338 bp，N50值為623 bp。進一步拼接得到54 207個Unigene，平均長度838 bp，N50值為1 328。按序列相似度對Unigene進行基因家族聚類，相似度大于70%的命為Clusters，共26 413個，Singleton序列有27 794個，其中E值小于10-5的序列共43 751個（表1）。54 207個Unigene的注釋信息見http：//www.sciencedirect.com/science/article/pii/S1874391914001523。

2.1.2 Contig和Unigene序列長度分布 組裝的124 002個Contig中，核苷酸長度在100～200 bp的序列數(shù)有74 842，比例達到了60.36%；200～300 bp有17 520，占14.13%；300～400 bp有8 714，占7.03%；400～500 bp 的有4 538條，占3.66%；而≥ 500 nt的共有18 388條，比例為14.83%。由此可見，Contig片斷長度以100～200 bp為主（表2）。在Contig數(shù)據(jù)的基礎(chǔ)上，采用over-lap進一步拼接，共獲得了54 207個unigene片段，其中核苷酸長度在100～500 bp的序列數(shù)有25 509，比例高達47.06%，500～1 000 bp有12 727，比例為23.48%；1 000～1 500 bp有7 312，占13.49%；1 500～2 000 bp有 4 334條，占8.01%；≥2 000 bp有4 315，占7.96%（表3）。

2.1.3 Unigene的序列同源性分析 經(jīng)Blast比對無冗余核酸數(shù)據(jù)庫NR后，對Unigene進行分析，從比對的E值來看（圖1-A），發(fā)現(xiàn)29.68%的序列具有一定的同源性（10-30

序列的相似性比較結(jié)果表明，相似性高于80%的序列高達27.8%，相似性在60%～80%的序列高達45.7%，而26.5%的序列相似性在18%～60%間（圖1-B）。對獲得的 Unigene進行Blastn分析（圖1-C），大部分Unigene序列和葡萄、蓖麻、毛果楊、大豆、苜蓿、擬南芥的同源序列比對上。統(tǒng)計顯示29.9%的Unigene可比對到葡萄；27.6%的Unigene與蓖麻同源；依次為毛果楊（22.9%），大豆（6.3%），擬南芥（1.1%）等，這說明芒果與葡萄的進化關(guān)系較近。

2.2 Unigene功能注釋

通過Blastx程序在不同的數(shù)據(jù)庫進行搜索比對（表4），42 515條unigene比對到NR數(shù)據(jù)庫，26 380條Unigene比對到Swissprot數(shù)據(jù)庫，23 741條Unigene用于KEGG途徑分析，14 619條進行COG分析，35 198條被用于GO分類。

2.2.1 Unigene的GO功能分類分析 12 923條Unigene與GO數(shù)據(jù)庫中的基因具有相似性，且較多的單條Unigene與多種基因相對應，建立了302 479條對應關(guān)系，從而得到盡可能多的注釋和分類。芒果果實轉(zhuǎn)錄組中的Unigene根據(jù)GO功能可分為生物學過程、細胞組分和分子功能3大類58分支（表5）。生物學過程類別中，細胞代謝過程和代謝過程涉及的基因最多，分別有23 131和22 003 條，細胞組分類別中，組成細胞和細胞部分涉及的基因最多，均為28 286條，分子功能類別中，結(jié)合功能涉及的基因最多，有17 458條，其次是具有催化活性功能的基因，達到了16 842條，其他種類基因的表達豐度不盡相同。

2.2.2 Unigene的COG功能分類 將Unigene和COG數(shù)據(jù)庫進行比對，27 869條Unigene與數(shù)據(jù)庫中的基因具有相似性，且較多的單條Unigene能夠與多種基因相對應，建立了65 536條對應關(guān)系（表6）。Unigene根據(jù)功能大致可分為25類，Unigene的COG功能種類比較全面，涉及了大多數(shù)的生命活動，其中整體功能類的基因數(shù)量最多，有4 774 條，其次是轉(zhuǎn)錄功能（2 596），復制、重組和修復功能（2 268），翻譯后加工、蛋白質(zhì)折疊的伴侶蛋白（2 098），碳水化合物運輸和代謝（1 627），核結(jié)構(gòu)相關(guān)基因類的數(shù)量最少，只有4條，而未知功能的序列多達1 361條。COG分類結(jié)果表明這些Unigene涉及許多的生物學功能，這些未知基因可能參與芒果果實生長發(fā)育和成熟的調(diào)控。

2.2.3 Unigene的KEGG代謝途徑分類 芒果果實轉(zhuǎn)錄組數(shù)據(jù)有23 741個Unigene能注釋到128類代謝途徑中（表7），包括代謝途徑、次生代謝的生物合成通路、植物-病原互作、植物激素信號轉(zhuǎn)導、DNA剪切、RNA運輸、蛋白質(zhì)在內(nèi)質(zhì)網(wǎng)的加工、淀粉和蔗糖代謝、苯丙氨酸的生物合成、萜類化合物合成，脂類代謝，RNA降解、黃酮和黃酮醇的生物合成等，其中代謝途徑和次生代謝的生物合成所注釋的Unigene數(shù)最多，分別為4 969和2 450。

2.3 編碼蛋白框（CDS）的核苷酸和氨基酸序列分析

[5] Pandit， S S， Chidley， H G， Kulkarni， R S， et al. Cultivar relationships in mango based on fruit volatile profiles[J]. Food Chem， 2009， 114（1）： 363-372.

[6] Pandit S S， Kulkarni R S， Chidley H G， et al. Changes in volatile composition during fruit development and ripening of ‘Alphonso mango[J]. J Sci Food Agr， 2009， 89（12）： 2 071-2 081.

[7] Ma X W， Wu H X， Liu L Q， et al. Polyphenolic compounds and antioxidant properties in mango fruits[J]. SciHortic， 2011， 129（1）： 102-107.

[8] MantheyJ A， Perkins-Veazie P. Influences of harvest date and location on the levels of beta-carotene， ascorbic acid， total phenols， the in vitro antioxidant capacity， and phenolic profiles of five commercial varieties of mango（Mangifera indica L.）[J]. J Agric Food Chem， 2009， 57（22）： 10 825-10 830.

[9] Ornelas-Paz J D， Yahia E M， Gardea A A. Changes in external and internal color during postharvest ripening of ‘Manila and ‘Ataulfo mango fruit and relationship with carotenoid content determined by liquid chromatography-APcI（+）-time-of-flight mass spectrometry[J]. Postharvest Biol Tec， 2008， 50（2-3）： 145-152.

[10] Jiang Y， Joyce D C. Effects of 1-methylcyclopropene alone and in combination with polyethylene bags on the postharvest life of mango fruit[J]. Ann Appl Biol， 2000， 137（3）： 321-327.

[11] Schuster S C. Next-generation sequencing transforms today's biology[J]. Nat Methods， 2008， 5（1）： 16-18.

[12] Shendure J， Ji H. Next-generation DNA sequencing[J]. Nat Biotechnol， 2008， 26（10）： 1 135-1 145.

[13] Wang Z， Gerstein M， Snyder M. RNA-Seq： a revolutionary tool for transcriptomics[J]. Nat Rev Genet， 2009， 10（1）： 57-63.

[14] Feng C， Chen M， Xu C J， et al. Transcriptomic analysis of Chinese bayberry （Myrica rubra） fruit development and ripening using RNA-Seq[J]. BMC Genomics， 2012， 13： 19.

[15] Vimolmangkang S， Zheng D， Han Y， et al. Transcriptome analysis of the exocarp of apple fruit identifies light-induced genes involved in red color pigmentation[J]. Gene， 2014， 534（1）： 78-87.

[16] Liu G Q， Li W S， Zheng P H， et al. Transcriptomic analysis of ‘Suli pear （Pyrus pyrifolia white pear group） buds during the dormancy by RNA-Seq[J]. BMC Genomics， 2012， 13（1）： 1-18.

[17] Sweetman C， Wong D C， Ford C M， et al. Transcriptome analysis at four developmental stages of grape berry （Vitis vinifera cv. Shiraz） provides insights into regulated and coordinated gene expression[J]. BMC Genomics， 2012， 13（1）： 105.

[18] Yu K Q， Xu Q， Da X L， et al. Transcriptome changes during fruit development and ripening of sweet orange（Citrus sinensis）[J]. BMC Genomics， 2012， 13（1）： 10.

[19] Dautt-Castro M， Ochoa-Leyva A， Contreras-Vergara C A， et al. Mango （Mangifera indica L.） cv. Kent fruit mesocarp de novo transcriptome assembly identifies gene families important for ripening[J]. Front plant sci， 2015， 6（10）： 62.

[20] Luria N， Sela N， Yaari M， et al. De-novo assembly of mango fruit peel transcriptome reveals mechanisms of mango response to hot water treatment[J]. BMC Genomics， 2014， 15（1）： 957.

[21] Azim M K， Khan I A， Zhang Y. Characterization of mango （Mangifera indica L.） transcriptome and chloroplast genome[J]. Plant mol biol， 2014， 85（1-2）： 193-208.

[22] Shan L L， Li X， Wang P， et al. Characterization of cDNAs associated with lignification and their expression profiles in loquat fruit with different lignin accumulation[J]. Planta， 2008， 227（6）： 1 243-1 254.

[23] Iseli C， Jongeneel C V， Bucher P. ESTScan： a program for detecting， evaluating， and reconstructing potential coding regions in EST sequences.[C]// Proc Int Conf Intell Syst Mol Biol， 1999： 138-148.

[24] Li X Y， Sun H Y， Pei J B， et al. De novo sequencing and comparative analysis of the blueberry transcriptome to discover putative genes related to antioxidants[J]. Gene， 2012， 511（1）： 54-61.

[25] Ong W D， Voo L Y， Kumar V S. De novo assembly， characterization and functional annotation of pineapple fruit transcriptome through massively parallel sequencing[J]. PLOS One， 2012， 7（10）： e46937.

[26] Pandit S S， Kulkarni R S， Giri A P， et al. Expression profiling of various genes during the fruit development and ripening of mango[J]. Plant Physiol Biochem， 2010， 48（6）： 426-433.

[27] Sherman A， Rubinstein M， Eshed R， et al. Mango （Mangifera indica L.） germplasm diversity based on single nucleotide polymorphisms derived from the transcriptome[J]. BMC Plant Biol， 2015， 15： 277.

[28] 羅 純， 武紅霞， 姚全勝， 等. 芒果轉(zhuǎn)錄組中SSR位點信息分析與引物篩選[J]. 熱帶作物學報， 2015， 36（7）： 1 261-1 266.

[29] 王 明， 應東山， 王琴飛， 等. 基于454 GSFLX高通量測序的芒果微衛(wèi)星標記特征分析[J]. 南方農(nóng)業(yè)學報， 2015（10）： 1 745-1 749.

[30] AltermannE， Klaenhammer T R. PathwayVoyager： pathway mapping using the Kyoto Encyclopedia of Genes and Genomes（KEGG） database[J]. BMC Genomics， 2005， 6： 60.

[31] Blencowe B J， Ahmad S， Lee L J. Current-generation high-throughput sequencing： deepening insights into mammalian transcriptomes[J]. Genes Dev， 2009， 23（12）： 1 379-1 386.