荸薺優良品種‘桂蹄3號’球莖發育過程轉錄組測序分析

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

(1. 廣西壯族自治區農業科學院生物技術研究所，廣西 南寧 530007;2. 福建農林大學園藝學院，福建 福州 350000)

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

1 材料與方法

1.1 供試材料

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

1.2 總RNA提取和檢測

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

1.3 構建cDNA文庫和illumina測序

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

1.4 序列拼接和組裝

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

1.5 組裝后的unigenes的功能注釋

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

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

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

2 結果與分析

2.1 荸薺球莖轉錄組測序及組裝

荸薺轉錄組測序獲得39.96 Gb干凈讀序，Q30堿基百分比均不小于95.90 %；將Clean Data 進行組裝，獲得荸薺轉錄本223 182條，Unigene 90 542條，轉錄本和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)。荸薺測序數據與組裝結果比對率達75 %以上(表3)，組裝完整性較高，效果較好，可用于進一步生物信息學分析。

2.2 基因注釋和功能分類

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

表1 ‘桂蹄3號’荸薺樣品測序數據評估統計

表2 ‘桂蹄3號’荸薺組裝結果統計

表3 ‘桂蹄3號’荸薺測序數據與組裝結果的比對統計

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

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

表4 ‘桂蹄3號’荸薺Unigene注釋統計

表5 ‘桂蹄3號’荸薺差異表達基因數目統計

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

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

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

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

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

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

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

圖3 ‘桂蹄3號’荸薺球莖發育過程中差異表達基因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技術具有高通量、成本低、完成速度快等優點，目前該技術已廣泛應用于特異基因的挖掘及其功能預測。國內外對荸薺研究主要集中在繁殖方式、栽培技術、病蟲害防治等生理栽培等研究，對荸薺的分子生物學方面特別是功能基因方面研究較少。本研究利用轉錄組測序分析荸薺球莖發育，為荸薺新基因挖掘及代謝途徑研究奠定了良好基礎。

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

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

4 結 論

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