溪黃草轉錄組測序及生物信息學分析

王繼華，蔡時可，梅瑜，李漢章，楊少海

（1.廣東省農作遺傳改良重點實驗室，廣東省農業科學院作物研究所，廣東廣州 510640；2.廣東建鵬南藥種養有限公司，廣東廣州 510070）

溪黃草Rabdosia serra（Maxim.）Hara 為唇形科香茶菜屬多年生草本植物，喜陰涼濕潤環境，俗稱溪溝草、黃汁草等，多產于華南地區的廣東、廣西及臺灣[1-2]。全草均可入藥，味苦、性寒，具有清熱利濕、退黃、涼血散瘀的功效，用于治療濕熱瀉痢、跌打瘀腫、急性黃疸型肝炎、急性膽囊炎、口腔炎、腸炎等病證[3-5]。藥理學研究表明，溪黃草具有保肝、護肝、抗炎、抗癌作用，是多種中成藥和保健品的主要成分[3，6-7]。溪黃草是廣東省連州市的特產，連州溪黃草在2007 年獲得國家地理標志產品（國家質量監督檢驗檢疫總局，2007年第220號公告）。

溪黃草來源復雜，其基原植物為唇形科溪黃草、線紋香茶菜、以及線紋香茶菜的變種細花線紋香茶菜和狹基線紋香茶菜[8]。資源鑒定主要根據根、莖、葉、花、果實等形態鑒定及味道等進行區分[9-11]。溪黃草中主要含有萜類、黃酮類、多糖和酚酸類等化學成分，各種有效成分還具有協同增效作用[12]。不同的基源植物之間化學成分相近，但部分藥效成分的含量相差甚遠[13]。溪黃草基因組數據仍不完整，不能對有效成分的代謝途徑進行深入的挖掘。目前，高通量測序平臺的RNA-seq技術在解析藥用成分代謝通路、挖掘關鍵基因、開發分子遺傳標記等方面得到廣泛應用[14-16]。因此，本研究開展溪黃草的轉錄組測序及生物信息學分析，以期為解析其藥用物質合成代謝通路、挖掘關鍵調控基因以及開發分子標記奠定基礎。現將研究結果報道如下。

1 材料與方法

1.1實驗材料實驗材料為廣東省農業科學院作物研究所南藥資源圃栽培的溪黃草。于2018 年12 月份取樣，采集健壯植株的葉片，迅速用錫箔紙包裹并浸入液氮處理，隨后置于超-80 ℃冰箱保存。

1.2溪黃草RNA的提取采用生工生物工程（上海）股份有限公司的總RNA 提取試劑盒（B511311-0025）提取溪黃草總RNA，然后通過1%電泳凝膠檢測提取RNA的完整性。應用Invitrogen Qubit?2.0熒光計及試劑盒（Fluorometer Life Tech Invitrogen，Q32886）對總RNA進行定量。

1.3轉錄組測序與拼接組裝委托生工生物工程（上海）股份有限公司采用Illumina HiSeq2500 的高通量測序平臺進行轉錄組測序。測序得到的原始數據，通過FastQC軟件進行質量評估和Trimmomatic進行質量剪切，過濾掉接頭、低質量的序列（reads ＜ 35nt）、帶N 堿基的序列、低質堿基（Q值＜20）得到高質量的clean data[17]。應用Trinily 軟件對clean data 進行de novo拼接組裝，再采用RSeQC軟件去除轉錄本中的冗余序列，得到非冗余通用基因（universal gene，unigene）[18]。

1.4基因功能注釋采用基于局部比對算法的搜索工具（Basic Local Alignment Search Tool，BLAST）將組裝的unigene 與保守域數據庫（Conserved Domain Database，CDD）、真核生物蛋白質同源簇數據庫/蛋白相鄰的聚類（eu Karyotic ortholog groups/clusters of orthologous groups，KOG/COG）、非冗余（Non-redundant，NR）、核酸序列數據庫（Nucleotide Sequence Database，NT）、蛋白結構域預測（Protein Families Database of Alignments and Hidden Markov Models，PFAM）、Swissprot、TrEMBL 等多個數據庫比對得到功能注釋信息。使用Transdecoder進行編碼序列（Coding Sequence，CDS）預測。根據unigene 與Swissprot、TrEMBL 的注釋結果得到基因本體論（Gene Ontology，GO）功能注釋信息，利用京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）自動注釋服務器（KEGG Automatic Annotation Server，KAAS）得到KEGG注釋信息。

1.5基因結構分析使用微衛星識別工具（Microsatellite Identification Tool，MISA）軟件鑒定溪黃草unigene上存在的SSR 位點，并利用Primer 3軟件（http：//primer3.sourceforge.net/releases.php）設計SSR引物[19]。

2 結果與分析

2.1轉錄組測序與de novo組裝溪黃草cDNA 文庫的構建由生工生物工程（上海）股份有限公司完成。應用Illumina Hiseq 2500測序平臺測序，共獲得61 944 850條raw reads，總堿基數為9 291 727 500 bp。使用Trimmomatic 對原始測序數據進行處理，去掉含有帶接頭、低質量的序列，共得到60 234 786 clean reads，總堿基數目為8 704 764 735 bp，GC含量為51.21%，Q30 bases ratio達到95.64%，表明文庫構件質量良好，測序得到的數據準確可靠。使用Trinity 將clean reads 進行de novo組裝成轉錄本，共得到86 204 條轉錄本，平均長度為1 251.93 bp，N50為1 894 bp，序列長度大于500 bp的有59 403條，占總序列數目的68.91%。見表1。對Trinity拼裝得到的轉錄本去冗余，共獲得37 418 條unigene，平均長度為1 054.1 bp，N50為1 840 bp，其中20 426條序列長度大于500 bp，占總序列數目的54.59%，序列長度在1 000 bp以上的有14 226條，占總數的38.02%見。見表1、圖1。

表1 溪黃草轉錄組測序的結果Table 1 Summary of assembled transcripts and unigenes of the Rabdosia serra（Maxim.）Hara transcriptome

圖1 溪黃草unigene序列長度分布Figure 1 Length distribution of assembled unigenes of Rabdosia serra（Maxim.）Hara

2.2 Unigene功能注釋將溪黃草組裝后的unigene序 列 與 CDD、 KOG、 COG、 NR、 NT、 PFAM、Swissprot、TrEMBL 等多個數據庫進行比對，共有23 978 條（64.08%）unigene 在至少一個數據庫中獲得功能注釋，2 429 條（6.49%）unigene 在所有數據庫種均能獲得注釋，見表2。其中NR 數據庫中注釋到的unigene 數目（23 623 條）最多，占總unigene的63.13%，其次為TrEMBL 數據庫，注釋到的unigene 為 23 393 條，占 62.52%。KEGG 數據庫中注釋到的unigene 數目（3 674 條，9.82%）最少。尚有13 440條unigene沒有得到有效注釋，占35.92%。見表2。通過與NR 庫的比對，獲得溪黃草unigene序列與近緣種屬的近似情況并獲得同源序列的功能信息，共有23 616條unigene獲得注釋。匹配較多的物種主要有Sesamum indicum、Erythranthe guttata、Dorcoceras hygrometricum和Salvia miltiorrhiza，分別占58.88%、20.39%、3.16%和1.97%，見圖2。表明溪黃草與唇形目（管狀花目）Sesamum indicum的序列相似度最高。

表2 Unigene 的功能注釋Table 2 Functional annotation of assembled unigenes

2.3 KOG功能注釋將組裝得到的溪黃草unigene與KOG數據庫比對，共有11 473條unigene獲得注釋，按其功能共分為26 類。其中：信號轉導機制注釋到的unigene 數目最多，共有1 974 條，占15.35%；一般功能預測、翻譯后修飾，蛋白質轉運，分子伴侶伴侶、轉錄和碳水化合物的轉運和代謝注釋到的unigene 占比分別為11.42%、10.26%、5.50%和5.21%；僅有34條和4條unigene注釋對應到核結構和細胞運動性，分別占比0.26%和0.03%。見圖3。此外，676 條unigene 注釋到未知功能。由此可見，溪黃草unigene 涉及信號轉導功能最多，可為今后溪黃草代謝物質的調控研究提供寶貴資源。

圖2 NR數據庫的同源物種分類Figure 2 Species distribution of the top BLAST hits against the NR database for the assembled unigenes

2.4 GO功能注釋GO 數據庫是全面描述生物體中基因及其產物屬性的分類系統，主要分為生物過 程（biological process）、 細 胞 組 分（cellular component）及分子功能（molecular funtion）三大類，見圖4。根據GO 數據庫注釋結果，共有20 222 條溪黃草unigene注釋成功。這些unigene總共被劃分為65 個功能分類，生物過程注釋到的unigene 最多，共有42 872 個，占35.95%，其中注釋較多的功能分別為細胞過程（11 757個，7.88%）、代謝過程（9 935個，6.66%）、刺激響應（5 455個，3.66%）和生物調節（4 952 個，3.32%）。注釋到分子功能的unigene 數目為38 425個，占32.22%，其中注釋較多的功能分別為結合（11 625 個，7.79%）、催化活性（9 633，6.45%）和轉運活性（1 276個，0.86%）。37 943 個（31.82%）unigene 注釋到細胞組分，其中注釋較多的功能分別為細胞（14 297 個，9.58%）、細胞組分（14 265 個，9.56%）、細胞器（10 568 個，7.08%）、膜（8 331 個，5.58%）和膜組分（6 481 個，4.34%）。

圖3 溪黃草unigene的KOG功能分類Figure 3 KOG function classification of Rabdosia serra（Maxim.）Hara unigenes

圖4 溪黃草unigene的GO功能分類Figure 4 GO classification of Rabdosia serra（Maxim.）Hara unigenes

2.5 KEGG功能注釋根據KEGG 數據庫注釋結果，主要分為五大類功能，包括代謝（2 529 個，41.67%）、遺傳信息過程（1 298 個，21.39%）、細胞過程（638 個，10.51%）、環境信息過程（623 個，10.27%）和有機系統（981 個，16.16%）。根據unigene 參與的代謝過程，進一步劃分為32個功能分類，涉及289 個代謝通路，具體見圖5。在代謝中，與氨基酸代謝和脂質代謝相關的unigene 分別為305 個（5.03%）和275 個（4.53%）。在有機系統中，涉及環境適應的unigene 有162 個，占2.67%。在遺傳信息過程中注釋到unigene 最多的代謝過程為轉錄（212 個，3.49%），而在細胞過程中，與轉運和分解代謝的unigene 最多（310個，5.11%）。信號轉導在植物發育過程以及對外界刺激的響應過程起重要作用，本研究共檢測到590 條unigene 涉及信號轉導，可為下一步鑒定溪黃草的信號因子提供研究基礎。溪黃草主要含萜類、黃酮類、酚類、氨基酸等化學成分，其中二萜類化合物非常豐富。本研究鑒定到19個unigene涉及黃酮類物質的生物合成，12個unigene 涉及倍半萜類化合物和三萜類化合物生物合成，40個unigene 與萜類物質骨架生物合成相關，13個unigene 參與雙萜的生物合成，51個unigene 涉及苯丙素的生物合成。這些與次生代謝相關unigene 的鑒定結果為進一步解析溪黃草藥用成分物質的生物合成提供了可能。

圖5 溪黃草unigene的KEGG功能分類Figure 5 KEGG classification of Rabdosia serra（Maxim.）Hara unigenes

2.6 SSR分析見圖6。采用MISA 對組裝的unigene 進行SSR 檢測，并對SSR 的類型和密度進行統計。結果表明，在7 809條unigene中共鑒定到9 489 個 SSR 位點。其中，921 條 unigene 中檢測到944 個（9.95%）復雜重復類型的SSR 位點。最豐富的重復類型是雙堿基重復，共檢測到4 208 個位點，占44.35%；其次為單堿基重復（2 524 個，26.60%）、三堿基重復（1 685 個，17.76%）、六堿基重復（70 個，0.74%）和四堿基重復（42 個，0.44%）；最少的為五堿基重復，僅檢測到16 個位點，占0.17%。見圖6-A。在SSR 分布的密度上，雙堿基重復最高，達到142.2個/Mbp，五堿基重復最低，僅為0.48 個/Mbp。見圖6-B。基于9 489 個SSR 位點，使用Primer 3.0 設計引物，為進一步開發溪黃草的遺傳標記和近緣種屬的遺傳圖譜提供了基礎數據。

圖6 溪黃草unigene的SSR位點分析Figure 6 The analysis of SSR sites of Rabdosia serra unigenes

3 討論

近年來，隨著高通量測序技術的發展，憑借其產出數據量大、成本低、不需要參考基因組等優勢，在新基因發掘、功能基因鑒定和分子標記開發上的應用越來越廣泛。本研究基于Illumina HiSeq2500的高通量測序平臺對溪黃草葉片進行轉錄組測序，共獲得60 234 786 條clean reads，總堿基數目為8 704 764 735 bp，GC 含量為51.21%，Q30 bases ratio 達到95.64%，表明文庫構件質量良好。N50 是評價組裝序列完整性的重要指標。通過de nove組裝，本研究共獲得37 418 條unigene，平均長度為1 054.1 bp，N50為1 840 bp，與已經構建云南松（1 818 bp）、香榧（1 702 bp）的N50 長度相接近，但較中藥黃芩（797.64 bp）的N50 長，表明溪黃草轉錄組序列組裝質量較高[20-22]。

溪黃草組裝后的unigene 序列與多個公共數據庫比對，共有 23 978 條（64.08%）unigene 在至少1 個數據庫中獲得功能注釋，2 429 條（6.49%）unigene在所有數據庫中均能獲得注釋，但仍有約35.92%的unigene 沒有獲得注釋信息，可能是由于部分序列組裝長度過短，缺少保守的核心序列以及溪黃草基因組信息匱乏暫時無法獲得準確的功能注釋。在KOG 數據庫和KEGG 數據庫中分別鑒定到1 974條和590條unigene涉及信號轉導，可為下一步解析溪黃草響應外界刺激以及代謝物質的調控提供靶標基因。同時，還鑒定到19個unigene涉及黃酮類物質的合成，51個unigene 涉及苯丙素的生物合成，以及大量與萜類物質合成相關的unigene，這為明晰溪黃草黃酮類和萜類物質的合成途徑和代謝網絡提供了數據支撐。

SSR標記操作簡單、重復性好，與傳統方法相比，高通量能夠挖掘出大量的SSR 位點。甘草、鐵皮石斛、枸杞等中藥材均利用轉錄組測序開發了SSR 標記用于分子輔助育種和遺傳圖譜的構建[23-25]。本研究應用MISA 軟件共檢測到9 489 個SSR 位點，最豐富的重復類型是雙堿基重復（4 208 個位點），占44.35%，達到142.2 個/Mbp，該研究結果與中藥黃芩類似[20]。本研究對溪黃草轉錄組進行了初步的探究，彌補了溪黃草基因組信息的不足，為解析次級代謝物質合成通路及分子生物學方面的研究打下了基礎。溪黃草SSR 位點的發掘，可為溪黃草分子標記的開發、遺傳多樣性分析、遺傳圖譜構建等奠定理論基礎，為利用分子手段鑒定和區分溪黃草及其基源植物提供了依據。