葉興狀 劉 丹 羅佳佳 范輝華 張國防 劉 寶* 陳世品

(1.福建農(nóng)林大學(xué)林學(xué)院，福州 350002； 2.福建省林業(yè)科學(xué)研究院，福州 350012)

轉(zhuǎn)錄組(Transcriptome)能夠在缺乏基因資源的條件下，運(yùn)用轉(zhuǎn)錄組測序技術(shù)獲得物種的生長和代謝規(guī)律，并揭示其生物學(xué)特性與基因內(nèi)在關(guān)聯(lián)，同時(shí)可獲得物種絕大多數(shù)的轉(zhuǎn)錄產(chǎn)物[1]。轉(zhuǎn)錄組測序不僅能夠發(fā)掘得到新的代謝通路及功能基因還是開展EST-SSR分子標(biāo)記的前提，而EST-SSR分子標(biāo)記是目前最重要的群體遺傳多樣性研究方法之一。雖然第三代高通量測序已問世，但第二代高通量測序相較第三代測序具有成本低、操作簡便、運(yùn)行時(shí)間短等特點(diǎn)，是一種性價(jià)比高的基因序列研究手段。近年來已有多種藥用植物完成轉(zhuǎn)錄組測序，如鄭少華[2]等對淫羊藿(Epimediumbrevicornu)進(jìn)行轉(zhuǎn)錄組分析，找到類黃酮素代謝通路；Roberto A Barrero等[3]對狼毒(Stellerachamaejasme)進(jìn)行轉(zhuǎn)錄組分析，提出了前列腺素合成上游途徑的新見解，為批量產(chǎn)生這種化合物奠定理論基礎(chǔ)，將可用于HIV研究或治療患者；李響等[4]對冬蟲夏草(Chinese Caterpillar Fungus)進(jìn)行轉(zhuǎn)錄組分析，開發(fā)了一種用于合成主要藥用化合物(如蟲草素)的模型；Xu Guo等[5]對鐵皮石斛(Dendrobiumofficinale)莖進(jìn)行轉(zhuǎn)錄組分析，揭示了代表參與生物堿骨架合成的25個(gè)基因及69個(gè)獨(dú)特序列。

半楓荷(Semiliquidambarcathayensis)屬金縷梅科(Hamamelidaceae)半楓荷屬(Semiliquidambar)是我國特有的瀕危國家二級保護(hù)植物[6]，半楓荷在《中國生物多樣性紅色名錄—高等植物卷》中被評為VU級[7]，具有楓香屬(Liquidambar)和蕈樹屬(Altingia)兩屬的綜合性狀特征，異型葉是其重要形態(tài)特征，葉子兼具有楓香(Liquidambarformosana)和阿丁楓(Altingiachinensis)特征，半楓荷嫩葉部分紫紅，極為美觀，是南方為數(shù)不多的變色樹種，具有良好的園林綠化應(yīng)用前景。此外其材質(zhì)優(yōu)良、樹干通直、旋刨性良好、材質(zhì)優(yōu)良，可作旋刨制品及家具，不僅是珍貴用材樹種[8～9]，而且藥用價(jià)值極高，其根、枝、葉，都可入藥，有祛風(fēng)、除濕、活血、止痛、消腫等功效，可治療腰肌勞損、跌打淤積、風(fēng)濕性關(guān)節(jié)炎，且外傷常用其葉止血[10～12]。目前半楓荷的研究處于起步階段，主要涉及種群生態(tài)特征和藥理研究等領(lǐng)域，鮮有涉及分子生物學(xué)領(lǐng)域，胡君等[13]在四川省發(fā)現(xiàn)了半楓荷新分布，并調(diào)查半楓荷種群生境特征；廖娜等[14]發(fā)現(xiàn)半楓荷中的多酚具有抗菌和顯著抗氧化能力；周光雄等[15]首次從金縷半楓荷中分離得齊墩果酸、3-羰基齊墩果酸等9種有抗炎藥效的化合物，為半楓荷藥理研究提供基礎(chǔ)數(shù)據(jù)；此外孫靜[16]等還發(fā)現(xiàn)半楓荷根的活性成分齊墩果酸對病毒性肝炎的抗原具有抑制活性的作用。由于半楓荷用途廣、經(jīng)濟(jì)價(jià)值高，再加上天然更新困難，對生境要求苛刻，同時(shí)因根可入藥、樹根常被藥農(nóng)挖斷，人為破壞嚴(yán)重，目前僅在貴州、湖南、福建、廣西、江西、廣東等省份零星分布[11～12,17]，已瀕臨滅絕，亟需保護(hù)。

鑒于此，本研究選擇福建省沙縣、順昌野生半楓荷做為研究材料，借助高通量測序平臺Illumina HiSeq 2500首次對半楓荷葉進(jìn)行轉(zhuǎn)錄組測序，建立半楓荷基因文庫，然后利用生物信息學(xué)的方法對半楓荷進(jìn)行功能注釋，代謝途徑、CDS預(yù)測及SNP檢測等分析，以期初步了解半楓荷基因表達(dá)、分布情況，如參與藥效合成的酶的基因分布與表達(dá)情況，為半楓荷遺傳多樣性研究提供基礎(chǔ)數(shù)據(jù)，同時(shí)也為今后半楓荷功能基因的研究奠定基礎(chǔ)。

1 材料與方法

1.1 實(shí)驗(yàn)材料

本實(shí)驗(yàn)選取半楓荷葉為實(shí)驗(yàn)材料。于2018年3月份到福建省順昌縣、沙縣各選取1棵半楓荷，采集一年生葉子放入預(yù)先準(zhǔn)備的冰盒中帶回實(shí)驗(yàn)室，將葉片擦干去除污漬，后用錫箔紙將葉片包好做好標(biāo)記，并用液氮處理10～20分鐘，再用超低溫零下80℃冰箱保存，用于半楓荷引物開發(fā)及遺傳多樣性分析實(shí)驗(yàn)。

1.2 方法

1.2.1 半楓荷總RNA的提取

將采集的2棵半楓荷葉樣品混合在一起，然后用生工生物(上海)有限公司的EZ-10DNA away RNA Mini-prep Kit試劑盒提取總RNA，使用杭州奧盛儀器有限公司的Nano-200微量分光光度計(jì)測定RNA的濃度和質(zhì)量，如果OD260NM/OD280NM的范圍為1.8～2.2，28S/18S>2，RNA完整值(RIN)>8.5，OD260NM/OD230NM>1.8，則認(rèn)為提取的RNA符合建庫要求。

1.2.2 半楓荷測序與拼接組裝

采用Illumina Hi Seq2500的高通量測序技術(shù)對半楓荷的轉(zhuǎn)錄組進(jìn)行測序，測序得到的原始數(shù)據(jù)，里面含有帶接頭的、低質(zhì)量的序列。為了保證信息分析質(zhì)量，對原始數(shù)據(jù)過濾，得到Clean數(shù)據(jù)，參數(shù)如下：

①去除帶N堿基的序列；

②去除reads中的接頭序列；

③從reads 3′到5′方向開始去除低質(zhì)堿基(Q值<20)；

④使用滑窗法去除reads尾部質(zhì)量值在20以下的堿基(窗口大小為5 bp)；

⑤去除reads長度小于35nt的reads本身及其配對reads。

⑥從reads 5′到3′方向開始去除低質(zhì)堿基(Q值<20)；

將reads回帖到全部轉(zhuǎn)錄本，并歸一化后得到reads在轉(zhuǎn)錄本上的平均分布，發(fā)現(xiàn)有少許偏移，但未發(fā)生明顯偏移，符合建庫要求。然后將Clean數(shù)據(jù)利用Trinity[18]軟件(參數(shù)min_kmer_cov 2，其余默認(rèn))進(jìn)行拼接組裝。建庫過程中的PCR擴(kuò)增會(huì)產(chǎn)生冗余序列(duplicate reads)，指的是堿基排列完全相同的序列，如果建庫擴(kuò)增時(shí)不均勻，可能會(huì)導(dǎo)致部分基因冗余序列過高，造成冗余序列分布的異常，使用RSeQC軟件去除冗余序列，得到Unigenes。

1.2.3 基因功能注釋

運(yùn)用NCBI Blast+將轉(zhuǎn)錄本與CDD、KOG、NR、NT、GO、Swissprot、PFAM、KEGG、TrEMBL等9個(gè)數(shù)據(jù)庫進(jìn)行比對，得到其功能注釋信息。根據(jù)轉(zhuǎn)錄本與Swissprot、TrEMBL的注釋結(jié)果得到GO功能注釋信息。利用KAAS[19]與基因和蛋白質(zhì)數(shù)據(jù)庫進(jìn)行比對，得到轉(zhuǎn)錄本KEGG注釋信息。利用以上數(shù)據(jù)庫的分析結(jié)果，篩選Unigenes最優(yōu)比對片段作為該Unigenes的CDS,使用getorf[20]檢測Unigenes的開放閱讀框(Open Reading Frame，ORF)，采用hmmsearch[21]將ORF比對到轉(zhuǎn)錄因子蛋白結(jié)構(gòu)域(數(shù)據(jù)來源于PlantfTFDB)，然后根據(jù)PlantTFDB描述的轉(zhuǎn)錄因子家族特征對Unigenes進(jìn)行能力鑒定。最后使用BCFtools[22]根據(jù)Mapping結(jié)果找出可能的單核苷酸多態(tài)性(Single Nucleotide Polymorphsims，SNP)位點(diǎn)，進(jìn)行SNP分析。

2 結(jié)果與分析

2.1 半楓荷轉(zhuǎn)錄本組裝與CDS預(yù)測分析

通過使用軟件將clean數(shù)據(jù)Denovo組裝成轉(zhuǎn)錄本共得到127 249個(gè)Transcripts。對半楓荷轉(zhuǎn)錄本長度進(jìn)行分析(圖1A,表1)，經(jīng)組裝、去冗長等處理得到77 629個(gè)Unigenes，200～500 bp范圍內(nèi)的Unigenes最多，為36 630條，其次是300～400 bp、400～500 bp和大于2 000 bp。最短長度為201 bp，最長10 975 bp，N50長度1 105 bp，總長度48 512 368 bp，Unigene平均長度為624.93 bp，其中長度分布在200～500 bp的Unigenes占比67.88%，分布在500～1 000 bp的有11 296個(gè)，占比14.55%，長度大于1 000 bp的有13 635個(gè)，占比17.57%。同時(shí)，為了更好的發(fā)掘半楓荷功能基因和便于遺傳指紋圖譜的構(gòu)建，本研究對半楓荷編碼序列(Coding Sequence，CDS)進(jìn)行預(yù)測，獲取57 671個(gè)Unigenes作為CDS。CDS序列1 000 nt以上的有8 050個(gè)，僅占13.96%，100～500 nt的序列有40 608個(gè)，占70.41%，其中200～300 nt的序列有20 615個(gè)，占總CDS序列的37.75%，500～1 000 nt的序列僅占15.63%。將圖1B與圖1A對比，序列分布更加集中于100～600 nt。

圖1 半楓荷Unigene的長度和CDS長度的分布圖Fig.1 Distribution of unigene lengths and CDS length for S.cathayensis

項(xiàng)目Item數(shù)量No.≥500bp≥1000bpN50N90最長值Max len最小值Min Len總長度Total len平均長度Average lenTranscript127249572803214413052851079520197669341767.55Unigene77629249311363511052451079520148512368624.93

2.2 功能注釋及COG分類

將denovo組裝得到的77 629個(gè)基因序列分別與CDD、KOG、NR、NT、GO、Swissprot、PFAM、KEGG、TrEMBL這9個(gè)數(shù)據(jù)庫進(jìn)行比對，成功注釋到9個(gè)數(shù)據(jù)庫中(表2)；在9個(gè)數(shù)據(jù)庫中至少1個(gè)數(shù)據(jù)庫注釋成功的Unigenes數(shù)為45 293，占總Unigenes數(shù)的58.35%；共有2 692個(gè)Unigenes與9個(gè)數(shù)據(jù)庫都能匹配成功，占總Unigenes數(shù)的3.47%(表2)；而只能注釋到NR、KEGG、COG、Swissprot蛋白庫的Unigenes分別只有6 971、0、120、3 362個(gè)(圖2)；此外，還有32 336(41.65%)個(gè)Unigenes序列未能獲得功能注釋信息，這有可能是測序技術(shù)不同引起的。

表2 Unigene注釋結(jié)果統(tǒng)計(jì)表

圖2 半楓荷Unigene的功能注釋維恩圖Fig.2 Vean diagram functional annotation distribution of unigenes of transcriptome for S.cathayensis

2.3 KOG功能注釋

將組裝得到的Unigenes比對到KOG蛋白質(zhì)庫中，半楓荷KOG功能注釋獲信息結(jié)果顯示細(xì)分為25個(gè)子類、共有22 802個(gè)Unigenes被注釋(占Unigenes總數(shù)的29.37%)，其各個(gè)子類的基因表達(dá)豐度差異較大，共包括25 253個(gè)功能注釋信息，基本涵蓋了半楓荷大部分的生命活動(dòng)。其中，信號傳遞機(jī)制，翻譯后修飾、蛋白質(zhì)折疊和分子伴侶，一般功能基因最多，分別有2 858(11.32%)、2 699(10.69%)和2 552(10.11%)。而核結(jié)構(gòu)和細(xì)胞活性較少，分別僅有55(0.22%)和6(0.02%)個(gè)。另有971(占3.85%)個(gè)功能注釋信息未明確其生物學(xué)功能(圖3)。

圖3 半楓荷轉(zhuǎn)錄組Unigenes的KOG功能分布圖 A.RNA加工和修飾；B.染色體結(jié)構(gòu)和動(dòng)力學(xué)；C.能源產(chǎn)生與轉(zhuǎn)化；D.細(xì)胞周期調(diào)控，細(xì)胞分裂，染色體分離；E.氨基酸轉(zhuǎn)運(yùn)和代謝；F.核酸轉(zhuǎn)運(yùn)和代謝；G.碳水化合物轉(zhuǎn)運(yùn)和代謝；H.輔酶轉(zhuǎn)運(yùn)和代謝；I.脂類轉(zhuǎn)運(yùn)和代謝；J.翻譯，核糖體結(jié)構(gòu)和生物發(fā)生；K.轉(zhuǎn)錄；L.復(fù)制，重組和修飾；M.細(xì)胞壁/細(xì)胞膜生物發(fā)生；N.細(xì)胞活性；O.翻譯后修飾,蛋白翻轉(zhuǎn),伴侶；P.無機(jī)離子轉(zhuǎn)運(yùn)和代謝；Q.次生代謝物合成,轉(zhuǎn)運(yùn)和代謝；R.只有一般功能預(yù)測；S.未知功能；T.信號傳遞機(jī)制；U.細(xì)胞間運(yùn)輸,分泌物和囊泡運(yùn)動(dòng)；V.防御機(jī)制；W.細(xì)胞外結(jié)構(gòu)；Y.核結(jié)構(gòu)；Z.細(xì)胞骨架Fig.3 KOG functional annotation distribution of unigenes of transcriptome for S.cathayensis A. RNA processing and modification; B. Chromatin structure and dynamics; C. Energy production and conversion; D. Cell cycle control,cell division,chromosome partitioning; E. Amino acid transport and metabolism; F. Nucleotide transport and metabolism; G. Carbohydrate transport and metabolism; H. Coenzyme transport and metabolism; I. Lipid transport and metabolism; J. Translation,ribosomal structure and biogenesis; K. Transcription; L. Replication,recombination and repair; M. Cell wall/membrane/envelope biogenesis; N. Cell motility; O. Posttranslational modification,protein turnover,chaperones; P. Inorganic ion transport and metabolism; Q. Secondary metabolites biosynthesis,transport and catabolism; R. General function prediction only; S. Function unknown; T. Signal transduction mechanisms; U:Intracellular trafficking,secretion,and vesicular transport; V. Defense mechanisms; W. Extracellular structures; Y. Nuclear structure; Z. Cytoskeleton

圖4 半楓荷轉(zhuǎn)錄組Unigenes的NR注釋物種分布圖Fig.4 NR annotated species distribution of Unigenes of transcriptome for S.cathayensis

2.4 NR功能注釋

將獲得的半楓荷Unigenes比對到NR庫中(E≤1e-5)，NR功能注釋獲信息結(jié)果顯示有39 076個(gè)Unigenes被注釋(占50.34%)，且與其它物種已知基因序列具有不同程度的匹配性。如注釋物種分布圖所示(圖4)：注釋為葡萄(Vitisvinifera)相關(guān)基因的序列最多，達(dá)10 111條，占NR庫中被注釋Unigenes的25.89%，其次，序列同源性大于2%的近緣物種還有大麥亞種(Hordeumvulgaresubsp. Vulgare)、荷花(Nelumbonucifera)、可可(Theobromacacao)、棗(Ziziphusjujuba)、麻風(fēng)樹(Jatrophacurcas)、蓖麻(Ricinuscommunis)、甜橙(Citrussinensis)、毛果楊(Populustrichocarpa)和梅(Prunusmume)，其余46.07%的注釋Unigenes分布于其他648個(gè)物種中。由于缺乏半楓荷基因組信息和技術(shù)局限性，仍有29(0.074%)個(gè)Unigenes在NR庫中未能獲得注釋。

2.5 GO功能注釋及分類

根據(jù)轉(zhuǎn)錄本與Swissprot、TrEMBL的注釋結(jié)果得到GO功能注釋信息。綜合描述半楓荷葉片中相關(guān)基因的生物學(xué)特征，從宏觀上了解其表達(dá)基因的功能分布情況，對獲得相應(yīng)GO條目的Unigenes進(jìn)行統(tǒng)計(jì)分析(圖5)，共有39 293個(gè)Unigenes獲得317 475個(gè)GO注釋，平均每條8.08個(gè)。把獲得GO注釋的Unigenes劃分為細(xì)胞組分、分子功能和生物學(xué)過程三個(gè)類型，分別各獲得133 668(42.10%)、53 236(16.77%)和130 571(41.13%)個(gè)注釋，可細(xì)分為65個(gè)子類，其中，細(xì)胞、細(xì)胞組分、結(jié)合、細(xì)胞器、催化活性、細(xì)胞過程、代謝過程和單一有機(jī)體過程等獲得注釋較多。

2.6 KEGG代謝通路分析

將半楓荷的Unigenes序列映射到KEGG代謝庫(E≤1e-5)，根據(jù)注釋信息對其可能參與或涉及的代謝途徑進(jìn)行統(tǒng)計(jì)分析。結(jié)果表明，4 029條Unigenes獲得注釋(5.19%)；其參與的代謝通路可歸為5大類別、33個(gè)子類。由圖6可知，5種代謝通路大類中，代謝相關(guān)的通路獲得2 786個(gè)注釋，占41.99%，其次是遺傳信息處理和生物系統(tǒng)分別獲得1 562(23.55%)、1 036(15.62%)，再次是細(xì)胞過程與環(huán)境信息處理分別獲得650(9.80%)、600(9.04%)。進(jìn)一步細(xì)分為33子類代謝途徑，其中，翻譯獲得注釋最多，為722(17.92%)，其次為信號傳導(dǎo)、折疊排序和退化、碳水化合物代謝。共有8 923條次Unigenes歸入286條代謝途徑，按數(shù)量排列，靠前的有15個(gè)代謝通路(表3)。核糖體代謝途徑注釋到的Unigenes數(shù)量最多，有389(4.36%)條Unigenes，其次為碳代謝(221，2.48%)、氨基酸的生物合成(193，2.16%)和內(nèi)質(zhì)網(wǎng)蛋白處理(171，1.92%)等途徑。

在半楓荷的轉(zhuǎn)錄組數(shù)據(jù)中，有64條Unigenes映射到苯丙素生物合成(ko00940)；有65條Unigenes映射到萜類物質(zhì)生物合成(ko00904，ko00902，ko00909，ko00900)；有27條Unigenes映射到黃酮類物質(zhì)生物合成(ko00941，ko00942，ko00944，ko00965)；有2條Unigenes映射到硫甙生物合成(ko00966)，可合成糖苷類物質(zhì)；有18條Unigenes映射到了合成生物堿類(ko00232、ko00950)通路；有1條unigene映射到了黃酮和黃酮醇的生物合成(ko00944)通路，有1條Unigene映射到了異黃酮的生物合成(ko00943)通路(表4)。半楓荷中多條次生代謝途徑及相關(guān)Unigenes的發(fā)現(xiàn)，說明其次生代謝生物合成途徑的復(fù)雜性，通過對半楓荷化學(xué)成分的分析研究，為后續(xù)從半楓荷中分離新化合物提供了線索，也為闡述半楓荷藥效物質(zhì)基礎(chǔ)提供了理論依據(jù)。

圖5 半楓荷轉(zhuǎn)錄組Unigenes的GO功能分類統(tǒng)計(jì)圖 1.抗氧化活性；2.結(jié)合；3.催化活性；4.電子載體活動(dòng)；5.metallochaperone activity；6.金屬伴活動(dòng)分子功能調(diào)節(jié)器；7.分子轉(zhuǎn)導(dǎo)活性；8.核酸結(jié)合轉(zhuǎn)錄因子活性；9.營養(yǎng)庫活性類；10.蛋白質(zhì)標(biāo)記；11.信號轉(zhuǎn)導(dǎo)活性；12.結(jié)構(gòu)分子活性；13.轉(zhuǎn)錄因子活動(dòng)與蛋白質(zhì)結(jié)合；14.傳譯調(diào)治活性；15.轉(zhuǎn)運(yùn)蛋白活性；16.受體活性；17.酶調(diào)節(jié)活性；18.受體調(diào)節(jié)活性；19.運(yùn)轉(zhuǎn)狀態(tài)；20.生物粘附；21.生物調(diào)節(jié)；22.細(xì)胞聚集；23.細(xì)胞殺傷；24.細(xì)胞成分組織或生物合成；25.細(xì)胞過程；26.解毒；27.發(fā)育過程；28.生長；29.免疫系統(tǒng)過程；30.定位；31.運(yùn)轉(zhuǎn)；32.代謝過程；33.多組織進(jìn)程；34.多細(xì)胞生物過程；35.生物過程負(fù)調(diào)控；36.生物過程的積極調(diào)節(jié)；37.生物過程的調(diào)節(jié)；38.繁殖；39.reproductive process；40.生殖過程對刺激的反應(yīng)；41.節(jié)律過程；42.信號傳導(dǎo)；43.單一有機(jī)體過程；44.建立斑點(diǎn)定位；45.細(xì)胞；46.細(xì)胞連接；47.細(xì)胞部分；48.細(xì)胞外基質(zhì)；49.胞外基質(zhì)成分；50.胞外區(qū)域；51.胞外區(qū)域部分；52.大分子復(fù)合體；53.膜；54.膜部分；55.膜封閉腔；56.擬核；57.細(xì)胞器；58.細(xì)胞器部分；59.其他有機(jī)體；60.其他有機(jī)體部分；61.超分子纖維；62.共質(zhì)體；63.突觸；64.突觸部分；65.病毒；66.病毒部分Fig.5 GO functional classification of unigenes of transcriptome for S.cathayensis 1.antioxidant activity； 2.binding； 3.catalytic activity； 4.electron carrier activity； 5.metallochaperone activity； 6.molecular function regulator； 7.molecular transducer activity； 8.nucleic acid binding transcription factor activity； 9.nutrient reservoir activity； 10.protein tag； 11.signal transducer activity； 12.structural molecule activity； 13.transcription factor activity,protein binding； 14.translation regulator activity； 15.transporter activity； 16.receptor activity； 17.enzyme regulator activity； 18.receptor regulator activity； 19.behavior； 20.biological adhesion； 21.biological regulation； 22.cell aggregation； 23.cell killing； 24.cellular component organization or biogenesis； 25.cellular process； 26.detoxification； 27.developmental process； 28.growth； 29.immune system process； 30.localization； 31.locomotion； 32.metabolic process； 33.multi-organism process； 34.multicellular organismal process； 35.negative regulation of biological process； 36.positive regulation of biological process； 37.regulation of biological process； 38.reproduction； 39.reproductive process； 40.response to stimulus； 41.rhythmic process； 42.signaling； 43.single-organism process； 44.establishment of localization； 45.cell； 46.cell junction； 47.cell part； 48.extracellular matrix； 49.extracellular matrix component； 50.extracellular region； 51.extracellular region part； 52.macromolecular complex； 53.membrane； 54.membrane part； 55.membrane-enclosed lumen； 56.nucleoid； 57.organelle； 58.organelle part； 59.other organism； 60.other organism part； 61.supramolecular fiber； 62.symplast； 63.synapse； 64.synapse part； 65.virion； 66.virion part

圖6 半楓荷轉(zhuǎn)錄組Unigenes的KEGG功能注釋分布統(tǒng)計(jì)圖 1.細(xì)胞生長和死亡；2.細(xì)胞運(yùn)動(dòng)；3.細(xì)胞群體；4.運(yùn)輸和分解代謝；5.膜運(yùn)輸；6.信號轉(zhuǎn)導(dǎo)；7.信號分子和相互作用；8.折疊排序與退化；9.復(fù)制和修復(fù)；10.轉(zhuǎn)錄；11.翻譯；12.氨基酸代謝；13.其他次級代謝產(chǎn)物的生物合成；14.碳水化合物代謝；15.能量代謝；16.多糖合成與代謝；17.類脂(化合)物代謝作用；18.輔助因子和維生素的代謝；19.其他氨基酸的代謝；20.萜類化合物和聚酮化合物的代謝；21.核苷酸代謝；22.概觀；23.異生素生物降解和新陳代謝； 24.老化；25.循環(huán)系統(tǒng)；26.進(jìn)化；27.消化系統(tǒng)；28.內(nèi)分泌系統(tǒng)；29.環(huán)境適應(yīng)；30.排泄系統(tǒng)；31.免疫系統(tǒng)；32.神經(jīng)系統(tǒng)；33.感官系統(tǒng)Fig.6 KEGG functional annotation distribution of unigenes of transcriptome for S.cathayensis 1.Cell growth and death； 2.Cell motility； 3.Cellular community； 4.Transport and catabolism； 5.Membrane transport； 6.Signal transduction； 7.Signaling molecules and interaction； 8.Folding sorting and degradation； 9.Replication and repair； 10.Transcription； 11.Translation； 12.Amino acid metabolism； 13.Biosynthesis of other secondary metabolites； 14.Carbohydrate metabolism； 15.Energy metabolism； 16.Glycan biosynthesis and metabolism； 17.Lipid metabolism； 18.Metabolism of cofactors and vitamins； 19.Metabolism of other amino acids； 20.Metabolism of terpenoids and polyketides； 21.Nucleotide metabolism； 22.Overview； 23.Xenobiotics biodegradation and metabolism； 24.Aging； 25.Circulatory system； 26.Development； 27.Digestive system； 28.Endocrine system； 29.Environmental adaptation； 30.Excretory system； 31.Immune system； 32.Nervous system； 33.Sensory system

Table3TopsixteenmetabolicpathwaysinvolvedinS.cathayensisunigenes

代謝通路PathwayUnigenes數(shù)量Number of unigenesUnigene比例Percent of unigenes(%)代謝通路編號Pathway ID核糖體Ribosome3894.36ko03010碳代謝Carbon metabolism2212.48ko01200氨基酸合成Biosynthesis of amino acids1932.16ko01230內(nèi)質(zhì)網(wǎng)蛋白加工Protein processing in endoplasmic reticulum1711.92ko04141剪接Spliceosome1581.77ko03040RNA轉(zhuǎn)運(yùn) RNA transport1541.73ko03013氧化磷酸化Oxidative phosphorylation1501.68ko00190嘌呤代謝Purine metabolism1361.52ko00230植物激素信號轉(zhuǎn)導(dǎo)Plant hormone signal transduction1281.43ko04075胞吞作用Endocytosis1101.23ko04144泛素介導(dǎo)的蛋白水解Ubiquitin mediated proteolysis1081.21ko04120植物病原體相互作用Plant-pathogen interaction1081.21ko04626淀粉和蔗糖代謝Starch and sucrose metabolism1001.12ko00500糖酵解/葡糖糖生成Glycolysis/Gluconeogenesis991.11ko00010mRNA監(jiān)測通路mRNA surveillance pathway961.08ko03015嘧啶代謝Pyrimidine metabolism941.05ko00240

2.7 TF編碼能力及SNP檢測

本研究共預(yù)測出88個(gè)基因家族共1 547個(gè)編碼轉(zhuǎn)錄因子(Transcription Factor，TF)的Unigenes(圖7)，另有57個(gè)轉(zhuǎn)錄因子未知具體基因家族；其中鋅指蛋白C2H2類TF家族占比最多，為151個(gè)(9.76%)，其次為AP2/ERF-ERF、MYB-related和bHLH家族TF，分別占75(4.85%)、69(4.46%)和68(4.40%)個(gè)TF，而HRT、SOH1、STAT、ULT所占比例最少，均為1(0.065%)個(gè)。同時(shí)預(yù)測出同源異型盒(Homeobox，HB)家族TF、包括HB-HD-ZIP、HB-PHD、HB-KNOX、HB-WOX、HB-BELL、HB-other；可能參與調(diào)控藥效成分合成有關(guān)的TF有AP2/ERF-ERF、WRKY、bHLH、bZIP、MYB等；還發(fā)現(xiàn)可能參與調(diào)節(jié)種子發(fā)育、萌發(fā)有關(guān)的TF有bZIP、MYB、MYB-related、C2C2-Dof、NAC、NF-YB、NF-YA、NF-YC、MADS-M-type、MADS-MIKC；與逆境抗性相關(guān)TF主要有：bZIP、AP2/ERF-ERF、AP2/ERFWRKY、HB類、MYB、MYB-related和HSF。這些TF功能的發(fā)現(xiàn)表明半楓荷轉(zhuǎn)錄調(diào)控的復(fù)雜性、交叉性和多樣性，也為后續(xù)研究半楓荷瀕危原因、藥效藥理研究、遺傳多樣性、指紋圖譜、連鎖圖譜構(gòu)建奠定基礎(chǔ)。

表4 半楓荷藥用活性成分相關(guān)的次生代謝產(chǎn)物生物合成途徑

圖7 半楓荷轉(zhuǎn)錄組Unigene主要的轉(zhuǎn)錄因子家族分類圖Fig.7 The main transcription factor family classification of unigenes of transcriptome for S.cathayensis

圖8 半楓荷轉(zhuǎn)錄組Unigene的SNP變異類型分布圖Fig.8 SNP variant type distribution of unigenes of transcriptome for S.cathayensis

此外，本研究檢測出125 917個(gè)單核苷酸多態(tài)性(Single Nucleotide Polymorphsims，SNP)位點(diǎn)(圖8)，分析發(fā)現(xiàn)半楓荷Unigenes序列上SNP分布不均勻，轉(zhuǎn)換突變類型(藍(lán)色)的SNP數(shù)量比顛換突變(紅色)的2倍略多。其中，轉(zhuǎn)換突變類型(A→G、C→T、G→A、T→C)有84 445個(gè)，占67.06%，顛換突變類型(A→C、A→T、C→A、C→G、G→C、G→T、T→A、T→G)有41 472個(gè)，占32.94%。半楓荷SNP分析，結(jié)合Unigenes的表達(dá)量及表型，可在mRNA層面的基因型差異區(qū)別于其它物種，通過構(gòu)建半楓荷SNP圖譜，對了解其親屬關(guān)系和物種起源進(jìn)化關(guān)系具有重要的生物學(xué)意義。

3 討論

本研究采用RNA-Seq技術(shù)首次對半楓荷葉片進(jìn)行了高通量轉(zhuǎn)錄組測序分析，拼接組裝后獲得77 629條Unigene，N50長度為1 105，平均長度為624.93 bp，均比川芎(Ligusticumchuanxiong)[23]、鷹嘴豆(Cicerarietinum)[24]、葫蘆巴(Trigonellafoenum-graecum)[25]更長，說明本研究所得序列中長片段較多，拼接效果較好。半楓荷SSR出現(xiàn)頻率為13.74%、比同是中藥材的野三七(Panaxstipuleanatus)[26](16.82%)、半夏(Pinelliaternata)[27](16.24%)更低，但比黨參(Codonopsispilosula)[28](12.22%)、蓮霧[29](Syzygiumsamarangense)(9.27%)、川芎(Ligusticumchuanxiong)[23](8.11%)更高，其中有13 635條Unigenes的長度在1 kb以上，進(jìn)一步說明半楓荷SSR位點(diǎn)分布較密集，可開發(fā)位點(diǎn)多。

所有的Unigene與CDD、KOG、NR、NT、GO、Swissprot、PFAM、KEGG、TrEMBL等9個(gè)數(shù)據(jù)庫進(jìn)行比對，共有45 293條(58.35%)成功獲得功能注釋，但仍有32 337條(41.65%)序列定位不清楚，得到功能注釋基因相對較少，這與香樟(Cinnamomumcamphora)[30]，紫背天葵(Begoniafimbristipula)[31]，柳樹(Salixbabylonica)[32]均出現(xiàn)大量未獲得注釋Unigenes的結(jié)果相似，可能是由于這些Unigenes本身為為非編碼序列，或是技術(shù)上的限制及已測序的半楓荷近緣種較少有關(guān)[33]。成功注釋到NR數(shù)據(jù)庫中的Unigenes中，注釋為葡萄同源種的有25.89%，遠(yuǎn)高于其他648個(gè)被注釋物種，出現(xiàn)這種情況，可能是由于葡萄與半楓荷的進(jìn)化史和生活史較為接近的緣故，亦可能因?yàn)槠咸丫哂袇⒖蓟蚪M。而與半楓荷能比對上僅有6個(gè)同科植物，共有43條Unigenes，其中楓香(Liquidambarformosana)、北美楓香(Liquidambarstyraciflua)、缺萼楓香(Liquidambaracalycina)、馬蹄荷(Exbucklandiapopulnea)、北美金縷梅(Hamamelisvirginiana)、紅花檵木(Loropetalumchinense)，分別有31、8、1、1、1和1，之所以會(huì)出現(xiàn)這種情況，可能是由于金縷梅科植物轉(zhuǎn)錄組、基因組信息嚴(yán)重缺乏、也可能是因?yàn)楸敬螠y序得到的Unigenes中半楓荷特有的新功能基因較多導(dǎo)致未能獲得相應(yīng)匹配。

對拼接組裝獲得的45 293條Unigenes進(jìn)行了代謝途徑分析和功能分類，并根據(jù)KEGG數(shù)據(jù)庫分析其可能參與或涉及的代謝通路，共獲得4 029(5.19%)個(gè)Unigenes注釋，注釋類別可細(xì)分為5個(gè)代謝通路大類、33類代謝途徑，其中定位到代謝相關(guān)通路的基因數(shù)最多，約占總注釋量的42%，進(jìn)一步證實(shí)半楓荷具有較強(qiáng)的代謝活動(dòng)能力。值得一提的是，本研究共發(fā)現(xiàn)286個(gè)具體KEGG代謝途徑分支參與到半楓荷碳水化合物代謝、生命活動(dòng)代謝、次生物質(zhì)代謝和遺傳信息處理過程中，其中，232個(gè)Unigenes參與到合成藥效次生代謝產(chǎn)物中、包括生物堿類、萜類、苯丙素類、黃酮類、醌類和酚類等，基于以上半楓荷藥用代謝產(chǎn)物合成途徑，初步構(gòu)建了半楓荷活性成分基因調(diào)控網(wǎng)絡(luò)。我們知道評價(jià)中藥品質(zhì)的優(yōu)劣，不止是看某一成分含量的高低，而是成分與成分間相互影響的整體特征，因此本研究的結(jié)果及后續(xù)將要開展的EST-SSR分子標(biāo)記為半楓荷藥效質(zhì)量評價(jià)標(biāo)準(zhǔn)提供基礎(chǔ)數(shù)據(jù)。另外，本研究用半楓荷葉片做轉(zhuǎn)錄組分析，而半楓荷根藥效比葉更大，可以預(yù)測半楓荷根相關(guān)藥效次級代謝產(chǎn)物代謝強(qiáng)度更大，但需進(jìn)一步驗(yàn)證。本研究根據(jù)組裝結(jié)果預(yù)測出88個(gè)基因家族共1 547個(gè)編碼轉(zhuǎn)錄因子的Unigenes，這些家族在半楓荷生長發(fā)育過程中發(fā)揮著不同的作用，其中C2H2鋅指蛋白基因家族數(shù)量最多，可能參與激素的信號傳導(dǎo)、生長發(fā)育、抗病信號途徑、細(xì)胞分化、脅迫應(yīng)答等相關(guān)活動(dòng)[34]。此外，由于條件限制，本研究僅僅是發(fā)現(xiàn)了轉(zhuǎn)錄因子與藥效成分合成有關(guān)，但具體調(diào)控的靶基因和作用機(jī)制還不清楚，需要進(jìn)一步研究。本研究通過轉(zhuǎn)錄組測序分析，共檢測出125 917個(gè)SNP多態(tài)位點(diǎn)，包括84 445個(gè)轉(zhuǎn)換突變類型，41 472個(gè)顛換突變類型，即轉(zhuǎn)換突變類型約占SNP的三分之二，與杜瑋南[35]等得出的結(jié)論一致。通過SNP研究將為今后開展半楓荷遺傳圖譜構(gòu)建、物種起源進(jìn)化、遺傳多樣性分析、品種鑒定以及分子標(biāo)記輔助選擇育種等研究提供重要參考。

通過生物信息學(xué)方法對半楓荷進(jìn)行CDS預(yù)測、功能注釋、代謝相關(guān)通路、SNP檢測和TF編碼能力預(yù)測，綜合KOG、GO和KEGG三大功能數(shù)據(jù)庫的注釋結(jié)果，表明半楓荷轉(zhuǎn)錄、復(fù)制和翻譯等基因表達(dá)豐富度較高，具有較強(qiáng)的代謝活動(dòng)和遺傳信息處理能力，同時(shí)找出半楓荷藥效物質(zhì)可能相關(guān)的次生代謝產(chǎn)物生物合成途徑及轉(zhuǎn)錄因子。所得轉(zhuǎn)錄組信息可作為金縷梅科植物基因組序列的重要組成部分，進(jìn)一步豐富了該科植物的基因信息庫，為進(jìn)一步研究藥用植物半楓荷的代謝途徑，以及該物種的分子遺傳育種和基因工程提供基礎(chǔ)數(shù)據(jù)，也將為保護(hù)和利用這一瀕危藥用植物資源奠定基礎(chǔ)。