大青轉錄組測序及生物信息學分析

劉敏，黃煒忠，何孟璐，梅瑜，王繼華

（1.廣東省農作物遺傳改良重點實驗室/廣東省農業科學院作物研究所，廣東廣州510640；2.廣東羅浮山國藥股份有限公司，廣東惠州516100）

大青（Clerodendron cyrtophyllumTurez）是馬鞭草科大青屬（Verbenaceae clerodendronL.）灌木或小喬木植物，又名大青、大青木、山大青和羊咪青等，多生于海拔1 700米以下的平原、山地林下或溪谷旁，分布于兩湖、兩廣及云貴高原等地，朝鮮、越南、馬來西亞也有分布[1-2]。大青全株皆可入藥，具有清熱解毒、祛風除濕的功效，主治乙腦、流腦、外感熱病熱盛煩渴、咽喉腫痛、口瘡、黃疸、熱毒痢、急性腸炎、癰疽腫毒、衄血、血淋、外傷出血等病癥?，F代藥理學研究表明，大青具有抑菌抗炎、抗病毒、利尿、鎮痛降壓的作用，是多種重要的中成藥及保健品的主要成分[3-4]。大青葉在我國南方也常作為食療藥材，其干燥葉也是涼茶的主要原料之一，在廣東多被制作成大青葉茶等保健產品[5]。

據《本草綱目》考證，民間藥用大青為馬鞭草科大青屬植物大青，但其種類繁多?，F根據地域區分就有廣東大青、江西大青、廣西大青、浙江大青等不同品種，這些不同品種之間的形態特征極其相似[6]。而現代中藥材的資源鑒定多以形態鑒定及氣味辨別為主，但不同大青的種類大多形態及氣味極其相近，極難辨別鑒定區分；且不同品種之間雖然化學成分相近，但部分藥效成分的含量卻差別極大[7]。隨著生物學技術的發展，基于高通量測序技術的轉錄組和生物信息學分析不受基因組信息的限制，在解析藥用活性成分合成途徑、挖掘相關功能基因、開發中藥分子標記等方面得到了廣泛的應用[8-12]。目前，許多藥用植物的轉錄組學研究已經開展，如金銀花[12]、地黃[13]、甘草[14]、茯苓[15]等，并鑒定到大量藥效活性物質的合成途徑相關基因，并開發了可靠的分子標記。目前，大青的研究多集中于化學成分鑒定及藥理活性研究等方面，而關于大青有效化學成分的合成代謝途徑關鍵基因的挖掘與鑒定尚未見報道。因此，本研究開展大青的轉錄組測序及生物信息學分析，以期為下一步解析其有效活性物質的合成代謝通路，挖掘其調控關鍵基因，開發分子標記提供基礎數據，現將研究結果報道如下。

1 材料與方法

1.1 實驗材料供試材料種植于廣東省農業科學院作物研究所特色作物與南藥資源圃。選取健康的大青植株的根、莖、葉，迅速用錫箔紙包裹并在液氮中冷凍，存儲于-80℃冰箱中備用。

1.2 大青總RNA的提取用TRIzol（上海生工）方法分別提取大青根、莖、葉總RNA，采用Qubit 2.0 RNA檢測試劑盒（Fluorometer Life Tech Invitrogen）分別對大青的根、莖、葉總RNA質量進行檢測和定量，選取光密度OD（260 nm）/OD（230 nm）值大于1.8的RNA樣品等量混合，建庫。

1.3 大青轉錄組測序與拼接組裝將檢測合格的RNA委托北京百邁客生物科技有限公司采用Illumina HiSeqTM2500的高通量測序平臺對大青的轉錄組進行測序。測序獲得原始數據后，按照標準分析流程（通過FastQC軟件進行質量評估，再使用Trimmomatic軟件去掉接頭和低質量序列）對數據進行處理，得到Clean Data用于de novo組裝，然后采用Trinily軟件將轉錄本拼接組裝獲得Unigene。

1.4 基因功能注釋采用Blastx程序將組裝的Unigene與保守域數據庫（Conserved Domain Database，CDD），非冗余蛋白序列（NCBI non-redundant protein sequences，NR），核酸序列數據庫（NCBI nonredundant nucleotide sequences，NT），蛋白結構域預測（Protein Families Database of Alignments and Hidden Markov Models，PFAM），真核生物蛋白質同源簇數據庫/蛋白質聚類（eu Karyotic ortholog groups/clusters of orthologous groups，KOG/COG），SwissProt，TrEMBL，KEGG直系同源數據庫（KEGGortholog database，KO），基因本體論（Gene Ontology，GO）等多個數據庫比對，獲得Unigene功能得到注釋信息、GO功能注釋信息。

1.5簡單重復序列（SSR）、單核苷酸多態性（SNP）檢測和引物設計采用微衛星識別工具（Microsatellite Identification Tool，MISA）軟 件 對大青轉錄組所有Unigene存在的SSR微點進行鑒定及分析，并使用Primer 3軟件（http://primer3.sourceforge.net/releases.php）設計SSR引物。

2 結果與分析

2.1 大青轉錄組測序與de novo組裝采用Illumina HiSeqTM2500高通量測序平臺進行大青轉錄組測序，獲得Raw Reads后，通過軟件去除接頭盒低質量序列后，共得到26 394 223條Clean Reads，長度為7 859 604 312 bp，GC含量為45.5%，Q20達98.17%，Q30達94.58%。該結果表明，大青轉錄組數據質量較高，能夠開展后續生物信息學分析。利用Trinity軟件將大青Clean Reads進行de novo組裝成轉錄本，共得到172 984條轉錄本，共計183 501 164 bp，平均長度1 060.8 bp，N50值為1 726 bp，其中序列長度大于500 bp的有105 429個，占總序列數的60.94%，序列長度大于1 000 bp的有64 345個，占總序列數的37.19%；去除冗余處理后得到100 191個Unigene，平均長度724.4 bp，N50值為1 055 bp，其中序列長度大于500 bp的Unigene數有42 108個，占總序列數的42.03%，序列長度大于1 000 bp的Unigene數有18 421個，占總序列數的18.39%。見表1、圖1。從組裝結果來看，轉錄本和Unigene的N50值均大于平均值，表明大青轉錄組的測序深度大，組裝效果較好。

圖1 大青Unigene序列長度分布Figure 1 Sequence length distribution of of unigenes of Clerodendron cyrtophyllum Turez

表1 大青轉錄組測序的結果Table 1 Results of transcriptome sequencing of Clerodendron cyrtophyllum Turez

2.2 大青轉錄組功能注釋及其分類使用Blast軟件將組裝的Unigene與九大生物信息學數據庫進行比對分析，得到大青轉錄組中的基因功能、代謝通路和轉錄因子等數據。通過比對分析，共有59 690（59.58%）個Unigene至少在一個數據庫中獲得信息注釋，長度介于300 bp到1 000 bp之間的Unigene有31 333個，長度大于1 000 bp的Unigene有16 325個。除COG數據庫外，在其他八大數據庫中能注釋的大青Unigene數目均在30%以上。其中：在NR和TrEMBL數據庫中注釋到的Unigene數目最多，分別達到58 952（58.84%）個和58 311（58.2%）個；在COG數據庫中注釋到的Unigene數目最少，僅為13 668（13.64%）個；49 278（49.18%）個Unigene在GO數據庫中得到注釋；38 260（38.19%）個Unigene在KEEG數據庫中得到注釋；32 453（32.39%）個Unigene在KOG數據庫中得到注釋；尚有40 501（40.42%）個Unigene尚未得到注釋。見表2。結果表明，大青轉錄組獲得的數據可以用于后續開展相應分析。

表2 大青Unigene注釋的統計結果Table 2 Statistical results of annotated unigenes of Clerodendron cyrtophyllum Turez

2.2.1 NR功能注釋大青的Unigene在NR數據庫中的比對結果表明，其與芝麻（Sesamum indicum）、玫瑰木屬（Rhodamnia argentea）、紫花風鈴木（Handroanthus impetiginosus）的基因組有一定的相似度。其中與芝麻的匹配度最高，共有11 140條Unigene與之相匹配，占整個NR數據庫注釋Unigene的18.9%，其次分別為玫瑰木屬和紫花風鈴木，分別占比為18.41%和10.25%，見圖2。從物種注釋數據可以看出，大青的Unigene序列全部可以與植物相匹配，且與唇形目（管狀花目）物種的序列相似度最高。但由于馬鞭草科大青屬植物的參考基因組數據比較少，因此并未從比對數據庫中發現同科屬的比對結果。

圖2 大青Unigene的NR注釋物種相似度分布Figure 2 Similarity distribution of NR annotated species of unigenes of Clerodendron cyrtophyllum Turez

2.2.2 KOG功能分類將組裝的大青Unigene與KOG數據庫比對分析，結果表明共有32 453條（32.39%）Unigene被注釋到25個KOG功能分類中。其中涉及一般功能預測的基因最多，共有6 359個（占19.59%），其后依次是翻譯后修飾，蛋白質轉運，伴侶基因（3 945個，12.16%），信號轉導機制基因（3 382個，10.42%），細胞內運輸，分泌和囊泡運輸（2 072個，6.38%），碳水化合物運輸與代謝（2 110個，6.5%）及翻譯，核糖體結構與生物起源（1 935個，5.96%），僅有11個基因注釋到細胞運動性，占比0.03%。此外，還有1 862個Unigene被注釋到未知功能。共有1 252個Unigene被注釋到次生代謝物生物合成、運輸和分解代謝，占總基因的比例為3.86%，這表明這些基因可能與大青中次生代謝物的生物合成、運輸和積累相關。見圖3。

圖3 大青Unigene的KOG功能分類Figure 3 KOG functional classification of unigenes of Clerodendron cyrtophyllum Turez

2.2.3 GO功能注釋將組裝的大青Unigene與GO數據庫比對分析，結果顯示共有49 278個（49.18%）Unigene被注釋到43個GO功能分類中。其中，生物過程注釋到的Unigene最多，共有76 354個，占76.2%，其中注釋較多的功能分別為細胞過程（26 146個，26.1%），代謝過程（22 970個，22.93%），生物調節（8 166個，8.15%），定位（5 359個，5.35%）及刺激響應（4 681個，4.67%）。注釋到分子功能中的Unigene數目為55 108個，占比為55%，其中注釋到較多功能的依次為結合（25 094個，25.05%），催化活性（22 737個，22.69%）及轉運活性（2 509個，2.5%）。注釋到細胞組分中的Unigene的數目為47 555個，占比為47.46%，其中注釋到較多的功能依次為細胞結構物質（27 133個，27.08%），細胞內物質（15 645個，15.62%）及蛋白質復合物（4 771個，4.76%）。見圖4。

圖4 大青Unigene的GO功能分類Figure 4 GO functional classification of unigenes of Clerodendron cyrtophyllum Turez

2.2.4 KEGG代謝通路分析大青中主要含有萜類、黃酮類、酚酸類及氨基酸等藥效活性成分，但大青缺少基因組信息作為參考，因此，對其基因功能注釋具有較高的應用價值。將組裝的大青Unigene與KEGG數據庫比對分析，結果顯示，共有38 260個Unigene得到注釋，并參與到細胞過程、環境信息過程、生物體系統、人類疾病、遺傳信息過程及代謝等六大類共136個代謝通路中。其中386個Unigene涉及苯丙烷生物合成（ko00940），165個Unigene涉及類黃酮生物合成（ko00941），37個Unigene參與異黃酮生物合成（ko00943），29個Unigene參與黃酮和黃酮醇生物合成（ko00944）；206個Unigene涉及萜類骨架生物合成（ko00900），130個Unigene涉及類胡蘿卜素生物合成（ko00906），80個Unigene涉及玉米素生物合成（ko00908），50個Unigene參 與 單 萜 生 物 合 成（ko00902），71個Unigene參 與 二 萜 生 物 合 成（ko00904），80個Unigene參與倍半萜和三萜生物合成（ko00909）。這些與次生代謝產物相關的Unigene的鑒定結果為進一步解析大青藥用成分物質的生物合成打下了基礎。見圖5。

圖5 大青Unigene的KEGG功能分類Figure 5 KEGG functional classification of unigenes of Clerodendron cyrtophyllum Turez

2.2.5 轉錄因子分析植物轉錄因子主要調節各種生物過程的基因表達模式，且普遍存在于植物中，如植物中常見的有bHLH、MYB和NAC等。通過研究這些TFs的作用機制是功能基因組學的重要環節。根據大青的對比結果顯示，共有3 543個Unigene被分為210個轉錄因子家族，其中C2H2轉錄因子類的Unigene數量最多，達到了115個，占比為3.24%，其次是bHLH、RLK-Pelle_DLSV、MYB-related、GRAS、C3H等。這些轉錄因子涉及到大青的生長發育、抗逆、次生代謝和合成等多種生物學途徑。大青轉錄因子的分析為進一步研究其基因調控與基因互作提供了一定的數據支持，也為提高其類黃酮、萜類酚酸類等有效成分的生物合成和抗逆性提供了理論依據。見圖6。

圖6 大青轉錄因子分類Figure 6 Classification of transcription factors of Clerodendron cyrtophyllum Turez

2.3 大青轉錄組中遺傳標記點位的分布SSR是檢測遺傳多樣性和構建遺傳圖譜的有效分子標記之一。利用MISA軟件從大青Unigene序列中鑒定潛在的SSR標記，并且進行統計。結果表明，在6 680條Unigene上共檢測到8 640個SSR位點。其中，400條Unigene中檢測到513個（5.8%）復雜重復類型的SSR位點。最豐富的重復類型是單堿基重復，共檢測到5 135個位點，占59.4%；其次為雙堿基重復（2 195個，25.4%），三堿基重復（745個，8.6%），混合堿基重復（500個，5.7%）和四堿基重復（42個，0.4%）；最少的為五堿基重復和六堿基重復，分別僅檢測到5個和5個位點，分別占0.06%和0.06%。根據8 640個SSR的位點和Unigene序列，使用Premier 3.0設計出相應的SSR擴增引物，為采用分子標記輔助育種技術的大青遺傳改良、群體結構的研究提供了基礎數據。見圖7。

圖7 大青Unigene的SSR位點分析Figure 7 SSR locus analysis of unigenes of Clerodendron cyrtophyllum Turez

3 討論

近年來，本草基因組學、轉錄組學技術的快速發展使得更多藥用植物的基因資源得以被保護及利用。通過轉錄組數據的分析研究，可挖掘重要植物功能基因和代謝通路，進一步構建中藥指紋圖譜，為天然藥物來源新途徑、種質資源鑒定、保存及良種選育提供分子基礎[16]；其次，通過對代謝通路關鍵基因的研究，可為中藥資源活性成分的生物合成與調控提供新的思路和方法；第三，通過轉錄水平的調節，提高藥用成分的產量與活性，尋找最佳中藥狀態，為重要的良種選育、規范種植、質量控制提供技術支撐[17]。

目前，對大青化學成分的分析和植株藥理作用的研究已有一定的經驗，但對大青的研究多限于整體水平實驗，還未見關于其基因組、轉錄組學的研究報道，這制約了大青次生代謝產物的生物合成途徑解析，不能充分挖掘其新的藥用價值以使其資源得以充分利用[18-19]。本研究采用Illumina HiSeq 2500的高通量測序平臺對大青的不同組織進行轉錄組測序并構建大青的Unigene庫。通過de novo組裝共獲得100 191條Unigene，N50為1 055 bp，長度大于1 000 bp的Unigene有18 421個，占總序列數的18.39%。N50的長度是評價轉錄組組裝質量的重要指標，大青轉錄組組裝的N50長度與茶樹（1 081 bp）的長度相近，但比中藥黃芩（797.64 bp）要長，說明大青的轉錄組組裝序列質量較高，能夠滿足后續數據分析的要求[20-21]。本研究通過與多個公共數據庫進行比對分析獲得Unigene的功能注釋，結果表明，大青組裝的Unigene共有59 690個（59.58%），Unigene在至少1個數據庫中得到注釋，其中以NR數據庫中比對到的Unigene最多，共有58 952個Unigene得到注釋，其結果顯示與管狀花目物種芝麻的相似度最高。這可能是由于馬鞭草科大青屬植物的種類少、研究基礎薄弱導致大青轉錄組Unigene并未比對到同科屬物種。另外，還有40 501個（40.42%）Unigene沒有在一個數據庫中獲得注釋，可能是由于組裝的Unigene序列太短而缺乏保守區域，或者數據庫中缺乏大青基因組、轉錄組信息導致未能匹配得到注釋。在KEGG數據庫中共注釋到38 260個Unigene參與136個代謝途徑，包括萜類、類胡蘿卜素、生物堿類、黃酮類等次級代謝的生物合成途徑。其中，407個Unigene參與萜類化合物的生物合成，231個Unigene參與黃酮類化合物的生物合成，386個Unigene參與苯丙烷素的生物合成。這些數據為進一步解析大青中活性成分的生物合成途徑提供了研究基礎。SSR標記操作簡單、重復性好。轉錄組數據中包含著大量SSR位點，隨著高通量測序技術的快速發展，利用這些位點開發的功能ESTSSR標記在群體遺傳學及標記輔助育種研究中顯得更加方便快捷。隨著大量藥用植物的轉錄組數據產生及功能基因得到注釋，也必然為藥用植物的功能EST-SSR標記開發奠定了基礎[22]。本研究利用MISA軟件共檢測到8 640個SSR位點，其中單堿基重復最為豐富（5 135個，59.4%），有很大一部分SSR位點的Unigene在各大數據庫中得到功能注釋，這也為大青后續功能EST-SSR分子標記開發提供了數據基礎。

本研究對大青轉錄組進行了初步的探究，彌補了大青基因組信息的不足，為解析大青次級代謝物質合成通路及分子生物學方面的研究奠定了基礎。大青SSR位點的發掘，可為其分子標記的開發、遺傳多樣性分析、種質資源鑒定與優選、分子標記育種等提供理論基礎，為利用分子手段鑒定、區分大青品種及評價其質量提供了依據。