江 文，何芳練，黃誠梅，高美萍, 2*，董偉清，蔣慧萍，黃詩宇

(1. 廣西壯族自治區(qū)農(nóng)業(yè)科學(xué)院生物技術(shù)研究所，廣西 南寧 530007;2. 福建農(nóng)林大學(xué)園藝學(xué)院，福建 福州 350000)

【研究意義】荸薺(EleocharistuberoseSchulut)俗稱馬蹄，屬單子葉莎草科多年生淺水草本植物，是我國14種重點發(fā)展的特色蔬菜之一。廣西作為荸薺種植和加工主產(chǎn)區(qū)，種植面積(2×104hm2)和產(chǎn)量(7.5×105t以上)均占全國的一半以上，荸薺產(chǎn)業(yè)現(xiàn)已成為廣西農(nóng)業(yè)支柱產(chǎn)業(yè)之一。然而對荸薺分子遺傳、轉(zhuǎn)錄組及基因組信息的研究十分缺乏，因此，開展荸薺球莖發(fā)育過程轉(zhuǎn)錄組研究，對荸薺產(chǎn)業(yè)的發(fā)展具有重要意義?！厩叭搜芯窟M展】荸薺前期的研究主要集中在種植栽培技術(shù)[1-3]、活性成分提取與利用[4-5]、組培苗及愈傷組織誘導(dǎo)試驗[6-7]等。利用轉(zhuǎn)錄組研究植物相關(guān)基因功能和結(jié)構(gòu)的方法已被廣泛應(yīng)用，分別針對植物不同組織部位、不同發(fā)育時期、不同品種等進行轉(zhuǎn)錄組測序，分析挖掘功能基因。張歡等[8]對杜梨根莖葉特異表達基因進行了RNA-Seq分析，探討了梨屬植物生長發(fā)育及組織間功能差異的分子機制；李和平等[9]對黃秋葵果實進行轉(zhuǎn)錄組測序，分析果實代謝途徑信息，并發(fā)現(xiàn)830個SSR位點；對睡蓮[10]、獼猴桃[11]及馬鈴薯[12]等植物的類似研究也均有報道，在植物抗病、抗逆的過程也有諸多應(yīng)用[13-17]。【本研究切入點】針對荸薺轉(zhuǎn)錄組及基因組信息的研究十分缺乏。本研究利用RNA-Seq高通量測序技術(shù)，首次在荸薺優(yōu)良品種‘桂蹄3號’球莖發(fā)育的不同時期進行轉(zhuǎn)錄組測序，分析球莖發(fā)育過程相關(guān)基因表達信息及品質(zhì)相關(guān)代謝途徑關(guān)鍵基因的表達情況，以期彌補荸薺品質(zhì)形成相關(guān)分子機制的空白?！緮M解決的關(guān)鍵問題】對發(fā)育期荸薺球莖進行轉(zhuǎn)錄組測序，采用生物信息學(xué)方法對得到的Uigene進行分類和功能注釋，篩選荸薺球莖發(fā)育過程相關(guān)基因，為進一步開展荸薺品質(zhì)形成分子機制等功能研究提供理論依據(jù)。

1 材料與方法

1.1 供試材料

試驗材料為荸薺優(yōu)良品種‘桂蹄3號’，從2017年10月20日開始，選擇生長健壯無病蟲害的球莖，每隔10 d取1次球莖，為球莖膨大初期(T01)，球莖膨大中期(T02)和膨大后期(T03)，清水沖洗干凈，液氮速凍后放于-80 ℃超低溫保存?zhèn)溆谩?/p>

1.2 總RNA提取和檢測

總RNA 的提取選用Trizol Reagent 方法，提取的RNA用無rnas DNaseI處理(TaKaRa生物科技有限公司)去除殘留的DNA。采用Nanodrop檢測RNA樣品的純度、Qubit 2.0和Aglient 2100檢測濃RNA濃度和完整性，以備合格的樣品進行轉(zhuǎn)錄組測序。

1.3 構(gòu)建cDNA文庫和illumina測序

樣品檢測合格后，構(gòu)建cDNA文庫，用帶有Oligo(dT)的磁珠富集mRNA；加入Fragmentation Buffer使得mRNA隨機被打斷為短片段；以打斷后的mRNA為模板，用六堿基隨機引物合成cDNA第一鏈，然后加入緩沖液、RNase 、HdNTPs和DNA polymerase I，合成第二條cDNA鏈。純化cDNA后，緩沖液洗脫，再經(jīng)過末端修復(fù)、加堿基A，連接測序接頭，再用AMPure XP beads選擇片段，瓊脂電泳回收目的片段，最后通過PCR擴增，完成cDNA文庫制備，采用HiSeq 2500進行測序，測序讀長為PE150[18]。

1.4 序列拼接和組裝

通過去除重復(fù)測序，未知N含量和低質(zhì)量讀序，獲得干凈讀序數(shù)據(jù)。轉(zhuǎn)錄組使用Trinity軟件對干凈序列進行從頭組裝獲得轉(zhuǎn)錄本，使用Tgicl 軟件進行聚類去冗余，獲得非冗余的unigenes。

1.5 組裝后的unigenes的功能注釋

利用BLAST系列軟件與網(wǎng)站公布的蛋白數(shù)據(jù)庫進行比對分析，針對Nt(核酸數(shù)據(jù)庫)、Nr(非冗余蛋白數(shù)據(jù)庫)、Swiss-Prot(經(jīng)注釋的蛋白質(zhì)序列數(shù)據(jù)庫)和TrEMBL(核酸序列數(shù)據(jù)庫)(E-value≤10-5)基于序列相似性檢索蛋白質(zhì)功能進行注釋。利用KEGG(Kyoto Encyclopedia of Genes and Genomes)、GO、(Gene Ontology)COG(Clusters of Orthologous Groups)分別進行比對預(yù)測分析分子功能及相關(guān)的代謝通路[19]。

1.6 差異基因表達(DGEs)分析

FPKM方法消除了不同基因長度和序列水平對基因表達計算的影響，因此直接采用FPKM值比較不同樣本間的基因表達差異，獲得樣品之間的差異表達基因集。采用校正后的P值，即FDR作為差異表達基因篩選的標(biāo)準，將FDR≤0.01和差異表達倍數(shù)絕對值log2比率≥1作為2個樣本之間基因表達的不同水平標(biāo)準。

2 結(jié)果與分析

2.1 荸薺球莖轉(zhuǎn)錄組測序及組裝

荸薺轉(zhuǎn)錄組測序獲得39.96 Gb干凈讀序，Q30堿基百分比均不小于95.90 %；將Clean Data 進行組裝，獲得荸薺轉(zhuǎn)錄本223 182條，Unigene 90 542條，轉(zhuǎn)錄本和Unigene 的N 50分別為2032和1119(表1～2)。Unigene長度分布情況為：300～500 bp 序列占51.3 %，500～1000 bp序列占26.89 %，1000～2000 bp序列占13.99 %，大于2000 bp序列占7.82 %(表2)。荸薺測序數(shù)據(jù)與組裝結(jié)果比對率達75 %以上(表3)，組裝完整性較高，效果較好，可用于進一步生物信息學(xué)分析。

2.2 基因注釋和功能分類

通過BLAST搜索比對，共有50 583個Unigene成功注釋(表4)，由于缺少參考基因組信息，部分序列未能注釋。與NCBI的Nr 數(shù)據(jù)庫進行比對獲得的基因注釋信息較全面，占全部注釋序列的92.10 %。與Swiss-prot 和Pfam數(shù)據(jù)庫成功比對的序列分別有30 811和36 554條。

表1 ‘桂蹄3號’荸薺樣品測序數(shù)據(jù)評估統(tǒng)計

表2 ‘桂蹄3號’荸薺組裝結(jié)果統(tǒng)計

表3 ‘桂蹄3號’荸薺測序數(shù)據(jù)與組裝結(jié)果的比對統(tǒng)計

2.3 荸薺球莖差異表達基因分析

2.3.1 荸薺球莖發(fā)育過程中差異表達基因比較 荸薺球莖發(fā)育過程中不同時期間基因表達差異見表5。結(jié)果顯示，膨大初期T01與膨大中期T02相比，有4027條Unigene差異表達，其中2157個基因上調(diào)表達，1870個基因下調(diào)表達；膨大中期T02與膨大后期T03相比，共檢測到3326個差異表達基因，其中1657個上調(diào)和1704個基因下調(diào)；T01與T03相比，檢測到2931條差異表達基因，其中1675個基因上調(diào)，1256個基因下調(diào)。由上述結(jié)果可知，球莖膨大初期的差異表達基因所占比例高于其他比較，說明DEGs在早期發(fā)育階段起著至關(guān)重要的作用。

表4 ‘桂蹄3號’荸薺Unigene注釋統(tǒng)計

表5 ‘桂蹄3號’荸薺差異表達基因數(shù)目統(tǒng)計

2.3.2 荸薺球莖發(fā)育過程中差異表達基因GO功能富集 ‘桂蹄3號’荸薺球莖發(fā)育過程中的3個發(fā)育階段T01、T02和T03經(jīng)兩兩比對，獲得的GO注釋種類和數(shù)量大體相同。由圖1可知，生物過程大多集中在代謝過程(metabolic process)、細胞活動過程(cellular process)和單生物代謝過程(single-organism process)。在細胞組分(cellular component)中，大多基因參與了細胞(cell)、細胞組成(cell part)和細胞器(organelle)的形成；在分子功能中，差異基因大多數(shù)富集到催化活性(catalytic activity)、結(jié)合(binding)和運轉(zhuǎn)活性(transporter activity)。

2.4 荸薺球莖發(fā)育過程中差異表達基因COG功能分類

對荸薺球莖進行COG數(shù)據(jù)庫比對分析結(jié)果(圖2)表明，有17 743條Unigene比對成功，根據(jù)功能分類可分為25類。一般功能預(yù)測(General function prediction only)，包含的Unigene最多，共4341個，占17.99 %；其次是轉(zhuǎn)錄、核糖體結(jié)構(gòu)和生物合成功能(Translation, ribosomal structure and biogenesis)，Unigene 2443個，占10.13 %；復(fù)制、重組和修復(fù)功能(Replication, recombination and repair)注釋到Unigene 1881個，占7.8 %； 碳水化合物運輸和代謝功能(Carbohydrate transport and metabolism)含Unigene1490，占6.18 %。最少的是核結(jié)構(gòu)(Nuclear structure)，注釋到9個Unigene，占0.4 %?？梢?，在球莖發(fā)育過程中，除復(fù)制、翻譯及轉(zhuǎn)錄等最基本的生命活動外，碳水化合物代謝占重要地位，表明碳水化合物代謝在荸薺球莖形成與發(fā)育中起到重要作用。

2.5 荸薺球莖發(fā)育過程中差異表達基因通路分析

對差異表達基因的代謝通路注釋分析有利于進一步解讀基因功能。結(jié)果表明，T01與T02 2個階段的差異表達基因共參與116種代謝，T02與T03、T01與T03階段分別涉及110和114種代謝途徑。差異表達基因參與的主要代謝途徑分類圖如圖3所示，其中淀粉-蔗糖代謝在3個文庫比對中占主要作用，其它主要代謝途徑是苯丙氨酸代謝(Phenypropanoid biosynthesis)、植物激素信號轉(zhuǎn)導(dǎo)(Plant hormone signal transduction)、氨基酸合成(Biosynthesis of amino acids)。

圖1 ‘桂蹄3號’荸薺球莖發(fā)育過程中差異表達基因GO功能富集Fig.1 Gene Ontology functional enrichment of differentially expressed genes in expansion of corms in water chestnuts ‘Guiti 3’

圖2 ‘桂蹄3號’荸薺球莖發(fā)育過程中差異表達基因COG注釋分類統(tǒng)計Fig.2 Clusters of Orthologous Groupsannotation classification of differentially expressed genes in expansion of corms in water chestnut ‘Guiti 3’

圖3 ‘桂蹄3號’荸薺球莖發(fā)育過程中差異表達基因KEGG分類Fig.3 Kyoto Encyclopedia of Genes and Genomes classification of differentially expressed genes in expansion of corms in water chestnut ‘Guiti 3’

3 討 論

RNA-Seq技術(shù)具有高通量、成本低、完成速度快等優(yōu)點，目前該技術(shù)已廣泛應(yīng)用于特異基因的挖掘及其功能預(yù)測。國內(nèi)外對荸薺研究主要集中在繁殖方式、栽培技術(shù)、病蟲害防治等生理栽培等研究，對荸薺的分子生物學(xué)方面特別是功能基因方面研究較少。本研究利用轉(zhuǎn)錄組測序分析荸薺球莖發(fā)育，為荸薺新基因挖掘及代謝途徑研究奠定了良好基礎(chǔ)。

淀粉是荸薺球莖中碳水化合物的重要貯藏物質(zhì)，荸薺的生長發(fā)育過程中淀粉不斷累積。研究表明，淀粉的生物合成和降解與蔗糖有著密切關(guān)系。本研究中注釋到的Unigene較多的是生物學(xué)過程中的代謝過程，大部分基因與荸薺的球莖發(fā)育代謝相關(guān)。KEGG功能注釋獲得20 667條Unigene，涉及116條代謝途徑，注釋到最多的代謝途徑是淀粉-糖代謝，共注釋到155條；其次是植物激素信號轉(zhuǎn)到途徑，有100條；苯丙氨酸代謝途徑共注釋到81條。可見，荸薺球莖發(fā)育過程中激素代謝通路也被激活。利用轉(zhuǎn)錄組測序技術(shù)研究淀粉-糖代謝在其他植物中也有相關(guān)報道。劉玉林等[19]應(yīng)用Illumina Solexa Hiseq 2000高通量測序技術(shù)對遼東櫟進行轉(zhuǎn)錄組測序，發(fā)掘出67條參與淀粉合成的Unigene 以及15 901個SSR 位點。Wang等[20]對馬鈴薯研究發(fā)現(xiàn)，AGPase基因很大程度上調(diào)節(jié)淀粉-糖代謝，該基因的沉默，引起淀粉產(chǎn)物的減少和可溶性糖的積累。程立寶等[21]利用轉(zhuǎn)錄組測序技術(shù)在蓮藕根狀莖膨大階段進行淀粉積累研究，發(fā)現(xiàn)了10個貯藏蛋白合成相關(guān)基因和5個淀粉合成相關(guān)基因(Lrgbss、Lrsbe1、Lrsbe2、LrsbeⅡ和LrsbeⅢ)，篩選出對蓮藕根狀莖的膨大起到重要作用的基因Lrplp8和Lrgbss。李雪艷等[22]采用二代高通量測序技術(shù)對百合小鱗莖形成與發(fā)育過程進行轉(zhuǎn)錄組測序，發(fā)現(xiàn)淀粉合成相關(guān)酶在小鱗莖中的表達豐度要高于母鱗片中的表達豐度，且以支鏈淀粉合成相關(guān)酶為主，而淀粉分解相關(guān)酶的基因表達量在母鱗片中較高；并在此基礎(chǔ)上挖掘百合小鱗莖形成與發(fā)育過程中差異表達的淀粉-蔗糖關(guān)鍵酶SuSy和INV，探討其在小鱗莖發(fā)育不同階段的表達模式。

后期可對淀粉-糖生物合成相關(guān)基因進行深入研究，分析出其關(guān)鍵基因。深入挖掘影響荸薺淀粉-糖代謝途徑的相關(guān)因素，可以為闡明荸薺淀粉生物合成機理提供理論依據(jù)，為荸薺的分子育種奠定基礎(chǔ)。

4 結(jié) 論

本研究首次利用高通量轉(zhuǎn)錄組測序技術(shù)建立了荸薺優(yōu)良品種‘桂蹄3號’球莖的轉(zhuǎn)錄組數(shù)據(jù)庫，組裝共得到223 182條Transcript和90 542條Unigene，平均長度為809 bp，N50為1119。共有50 583條Unigene在7個數(shù)據(jù)庫成功注釋。經(jīng)KEGG代謝通路比對，有20 667個基因獲得功能注釋，共有116條代謝途徑，其中注釋到最多的代謝途徑是淀粉-糖代謝，共注釋到155條。