大青轉(zhuǎn)錄組測序及生物信息學(xué)分析

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

（1.廣東省農(nóng)作物遺傳改良重點(diǎn)實(shí)驗(yàn)室/廣東省農(nóng)業(yè)科學(xué)院作物研究所，廣東廣州510640；2.廣東羅浮山國藥股份有限公司，廣東惠州516100）

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

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

1 材料與方法

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

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

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

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

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

2 結(jié)果與分析

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

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

表1 大青轉(zhuǎn)錄組測序的結(jié)果Table 1 Results of transcriptome sequencing of Clerodendron cyrtophyllum Turez

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

表2 大青Unigene注釋的統(tǒng)計(jì)結(jié)果Table 2 Statistical results of annotated unigenes of Clerodendron cyrtophyllum Turez

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

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

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

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

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

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

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

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

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

圖6 大青轉(zhuǎn)錄因子分類Figure 6 Classification of transcription factors of Clerodendron cyrtophyllum Turez

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

圖7 大青Unigene的SSR位點(diǎn)分析Figure 7 SSR locus analysis of unigenes of Clerodendron cyrtophyllum Turez

3 討論

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

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

本研究對大青轉(zhuǎn)錄組進(jìn)行了初步的探究，彌補(bǔ)了大青基因組信息的不足，為解析大青次級代謝物質(zhì)合成通路及分子生物學(xué)方面的研究奠定了基礎(chǔ)。大青SSR位點(diǎn)的發(fā)掘，可為其分子標(biāo)記的開發(fā)、遺傳多樣性分析、種質(zhì)資源鑒定與優(yōu)選、分子標(biāo)記育種等提供理論基礎(chǔ)，為利用分子手段鑒定、區(qū)分大青品種及評價(jià)其質(zhì)量提供了依據(jù)。