徐 浩，夏銘澤，張發(fā)起

（1.中國科學院高原生物適應與進化重點實驗室，中國科學院西北高原生物研究所，青海西寧810001；2.中國科學院大學生命科學學院，北京100039）

灌木亞菊（Ajania fruticulose（Ledeb.）Poljak）隸屬菊科（Asteraceae）亞菊屬（Ajania），小亞灌木，主要分布于我國西北地區(qū)及俄羅斯中亞地區(qū)550～4 400 m的荒漠及荒漠草原［1］。灌木亞菊在青藏高原地區(qū)資源量豐富，全草入藥，有止咳止血、排毒消炎、鎮(zhèn)定、祛痰及抗蛔蟲等藥效，還可用于活血化瘀、咽喉腫痛等疾病［2］。該植物具縱深發(fā)達的根，且蒸騰作用較弱，使得灌木亞菊的抗旱能力較強，廣布于干旱和半干旱的荒漠地帶［3］。此外，因其適口性好，可作為牛、羊等牲畜的食料［4］。目前亞菊屬植物的研究主要集中在化學成分及其藥理學方面，如Li等［5］以灌木亞地上部分為原料提取出兩種新的愈創(chuàng)木內(nèi)酯和四種黃酮類化合物，對黃嘌呤氧化酶有較好的抑制作用。張占欣［6］通過亞菊屬植物細葉亞菊（A.tenuifolia）的全草提取到了以萜類、黃酮及香豆素為主的35個化合物。Liu等［7］對細裂亞菊（A.przewal?skii）精油研究發(fā)現(xiàn)其具有一定的抗菌活性。Shi等［8］在柳葉亞菊（A.salicifolia）中分離得到了47種次生代謝產(chǎn)物，以此來探討它們在亞菊屬中的化學分類學意義。此外，通過研究亞菊屬植物（A.shiwogi?kuvar.kinokuniense）內(nèi)抗氧化酶對銹病的抗性來探明其對銹病的抗性機制［9］。通過測定同科植物菊花（Chrysanthemum morifolium）及其炮制品中綠原酸、木犀草苷等藥用成分的含量，提供炮制品合適的質(zhì)量控制參考方法［10］。然而該屬的分子生物學研究涉及較少，僅有少數(shù)基于高通量測序獲取葉綠體基因組開展系統(tǒng)發(fā)育分析等研究工作［11］，而在轉(zhuǎn)錄組與基因組方面的研究接近空白。

近年來，基因組學、轉(zhuǎn)錄組學、蛋白組學、代謝組學等組學技術的發(fā)展突飛猛進，特別是以轉(zhuǎn)錄組等功能基因組學在植物分子生物學研究中的應用越來越廣泛［12-13］。對接種鏈格孢菌（Alternaria tenui?ssima）的菊花（C.morifolium）葉片進行轉(zhuǎn)錄組分析，觀察到病原體識別、細胞壁修飾等差異轉(zhuǎn)錄的基因在病原體侵入后被上調(diào)［14］。油料作物紅花（Cartha?mustinctorius）的轉(zhuǎn)錄組分析成功預測出黃酮類化合物和種子油脂生物合成的相關基因［15］。低溫誘導的灰白銀膠菊（Parthenium argentatum）表皮組織轉(zhuǎn)錄組研究發(fā)現(xiàn)天然橡膠的生物合成與夜間低溫相關［16］。

青藏高原地區(qū)的灌木亞菊資源量巨大，有著巨大的開發(fā)潛力，且其主要分布在荒漠區(qū)，也具有重要的生態(tài)價值。然而，目前對該植物的研究涉及還較少，尤其在灌木亞菊的藥理作用、功能成分及作用機理研究較少［4］。因此，本研究利用Illumina HiseqTM測序平臺對灌木亞菊葉片進行轉(zhuǎn)錄組測序、組裝和生物信息學分析，挖掘藥用活性相關基因及代謝通路，并對其分子標記進行預測分析，為灌木亞菊的資源利用和保護、遺傳多樣性等研究提供基礎數(shù)據(jù)。

1 材料與方法

1.1 試驗材料

灌木亞菊新鮮、健康的葉片于夏季（7月初，植物花期）采集于青海省海南州共和縣（地理坐標N36°10′53″，E100°59′14″；海拔2 835 m）。分別采集3株植物的葉片經(jīng)75%乙醇和超純水處理后迅速置于液氮，帶回實驗室后置于-80℃的超低溫冰箱下保存?zhèn)溆谩{證標本（憑證標本號：zhang2018023）保存于中國科學院西北高原生物研究所青藏高原生物標本館（HNWP）。

1.2 試驗方法

1.2.1 灌木亞菊RNA的提取與文庫構建

利用Total RNA Extractor（Trizol）提取試劑盒提取灌木亞菊葉片（3株植物的葉片混樣）中總RNA，并用瓊脂糖凝膠電泳與Nanodrop 2000檢測總RNA的濃度與純度。檢測合格的RNA參照生工生物工程（上海）股份有限公司的標準流程進行文庫構建。文庫構建完成后，采用Qubit 2.0、Agilent 2100進行文庫質(zhì)量檢測，并使用Q-PCR方法對文庫的有效濃度進行準確定量（文庫有效濃度＞2 nM），以保證文庫質(zhì)量。

1.2.2 序列的測序及拼接

基于Illumina HiseqTM高通量測序平臺進行測序，經(jīng)CASAVA堿基識別（Base calling）把數(shù)據(jù)文件轉(zhuǎn)換為原始測序序列（Raw reads）。利用FastQC v0.11.2［17］對原 始 數(shù)據(jù)進行質(zhì)量 評 估，Trimmomatic v0.36（http：//www.usadellab.org/cms/？page=trimmo?matic）對含有帶接頭的、低質(zhì)量的序列進行數(shù)據(jù)質(zhì)控，得到高質(zhì)量的序列（Clean reads）。最后使用Trinity v2.4.0［18］將Clean reads進行從頭組裝，得到的轉(zhuǎn)錄本去冗余后，將每個轉(zhuǎn)錄本聚類中長度最長的作為Unigenes用于后續(xù)分析。

1.2.3 基因的功能注釋

為獲得全面的基因功能信息，通過NCBIBlast+v2.60［19］將Unigenes與CDD、KOG、NR、NT、PFAM、Swissprot、TrEMBL、GO、KEGG等數(shù)據(jù)庫比對。根據(jù)Swissprot、TrEMBL的蛋白注釋結果得到GO功能注釋信息。利用KAASv2.1［20］將Unigenes與KEGG數(shù)據(jù)庫比對進行代謝通路分析。根據(jù)轉(zhuǎn)錄本與數(shù)據(jù)庫Blast的比對結果通過軟件TransDecoder v3.0.1進行編碼序列（Coding sequence，CDS）預測。

1.2.4 SSR與SNP的預測

基于組裝好的轉(zhuǎn)錄本Unigenes序列，利用工具MISA v1.0［21］來 識 別Unigenes中 的 簡 單 重 復 序 列（Simple sequence repeat，SSR），其中一至六重復單元SSR的至少重復次數(shù)依次為：10、6、5、5、5、5次。結果包含復合型與完美型SSR，BCFtools v1.5［22］（參數(shù)設置：質(zhì)量值＞20，覆蓋度＞8）抽提Unigenes中的單核苷酸多態(tài)性SNP（Single nucleotide polymor?phisms，SNP）進行統(tǒng)計分析。

2 結果與分析

2.1 轉(zhuǎn)錄組測序數(shù)據(jù)分析

經(jīng)瓊脂糖凝膠電泳檢測總RNA的產(chǎn)量與質(zhì)量膠圖如圖1所示（本研究樣品編號為23，其中Marker為DNA Marker，其只為檢測基因組污染，并不代表RNA條帶大小），檢測合格后的樣本再進行后續(xù)轉(zhuǎn)錄組分析。灌木亞菊轉(zhuǎn)錄組共得到Raw reads 47 084 098條（平均長度為150 bp），質(zhì)控后Clean reads 45 983 774條，平均長度為145 bp，平均GC含量46.13%，Q30（堿基質(zhì)量在30以上的序列）的序列占95.01%，表明測序所得序列質(zhì)量能夠滿足后續(xù)分析。

圖1 總RNA瓊脂糖凝膠電泳圖Fig.1 Agarose gel electrophoresis of total RNA

基于Trinity從頭測序組裝和拼接共獲得103 820個Transcripts，其中挑選最長的轉(zhuǎn)錄本共得到53 760條Unigenes，N50長度為1 155 bp（表1）。對Unigenes的長度進行統(tǒng)計（圖2），200～300 nt之間數(shù)量最多，有19 461條，占比36.20%；300 nt～2 000 nt之間的有31 158條，占比57.96%；大于2 000 nt的序列共有3 141條，占比5.84%。其中，共有38 215條Unigenes可作為CDS序列，長度在200～300 nt之間的CDS序列占比最多（25.20%），與Unigenes長度的分布趨勢大體一致（圖3）。

圖2 灌木亞菊Unigenes的長度分布圖Fig.2 Length distruption of A.fruticulose Unigenes

圖3 灌木亞菊CDS的長度分布圖Fig.3 Length distruption of A.fruticulose CDS

表1 灌木亞菊轉(zhuǎn)錄組數(shù)據(jù)拼接結果Table 1 Assembly result of transcriptome data from A.fruticulose （單位：bp）

2.2 灌木亞菊Unigenes的功能注釋

將從頭組裝的Unigenes共53 760條分別與CDD、KOG、NR、NT、PFAM、Swissprot、TrEMBL、GO、KEGG數(shù)據(jù)庫進行比對，發(fā)現(xiàn)至少注釋到一個數(shù)據(jù)庫的Unigenes數(shù)量為33 766（62.81%）條，有2 142條Unigenes在9個數(shù)據(jù)庫中都能匹配成功，但仍有19 994條未能與上述數(shù)據(jù)庫成功匹配（表2）。將Unigenes與NR數(shù)據(jù)庫比對，共獲得注釋信息32 505條。與灌木亞菊Unigenes的NR注釋匹配數(shù)最多的十個物種中，同科多年生草本植物洋薊（Cynara car?dunculusvar.scolymus）與灌木亞菊的Unigenes相關序列最多，共18 898條（占比58.14%），剩余物種匹配度不高且與灌木亞菊親緣關系較遠（圖4）。

圖4 灌木亞菊NR庫的物種分布Fig.4 Species distribution of A.fruticulosa in NRdatabase

表2 灌木亞菊轉(zhuǎn)錄組注釋結果Table 2 The resultsof transcriptome annotation of A.fruticu?lose

2.3 GO數(shù)據(jù)庫功能分類注釋

基于Swissprot和TrEMBL兩種數(shù)據(jù)庫的蛋白注釋結果，根據(jù)Uniprot所得注釋信息得到Unigenes的GO注釋，共有27 892條Unigenes被注釋，獲得171 904條注釋信息，分屬于細胞組分（Cellular com?ponent）、分子功能（Molecular function）以及生物學過程（Biological process）三大類（圖5）。在與藥用成分相關基因注釋中，涉及甾體（366條）與黃酮類化合物（204條）的居多，分別歸屬于生物學過程及分子功能類。苯丙素類主要涉及苯丙素的生物合成及代謝，有68條信息被注釋出，而中草藥中常見的萜類化合物也有48條信息被注釋。

圖5 灌木亞菊Unigenes的GO功能分類統(tǒng)計圖Fig.5 GOfunctional classification statisticsdiagramof A.fruticulose Unigene

2.4 KEGG數(shù)據(jù)庫的代謝通路分析

為了解灌木亞菊基因所參與的代謝途徑，將Unigenes序列與KEGG數(shù)據(jù)庫比對分析，共3 619條Unigenes獲得注釋，包含有5 892個注釋信息，共歸為細胞進程、環(huán)境信息處理、遺傳信息處理、新陳代謝、有機體系統(tǒng)五大類（圖6）。

圖6 灌木亞菊Unigenes的KEGG代謝通路分析Fig.6 Metabolic pathways analysis of A.fruticulose Unigenes in KEGGdatabase

為了解灌木亞菊化學成分的相關信息，篩選出了包含主要藥用活性成分的相關通路（表3）。通路涉及苯丙素類、黃酮類、醌類、生物堿類、萜類、有機酸及酚類這些常見藥用活性成分。在萜類化合物的生物合成途徑中，倍半萜合成通路較為顯著，其主要是由焦磷酸金合歡酯衍生而成，多以五元內(nèi)酯環(huán)的形式存在，是亞菊屬植物中特有的萜類成分。黃酮類化合物也是重要的藥用活性成分，且涉及該通路基因的表達比重上升，它通過抑制酪氨酸酶的活性而具有顯著的抗氧化能力。此外，反式肉桂酸等芳香族衍生物表達量的上升促進了苯丙素類生物合成的化學反應；通過接受呼吸鏈電子傳遞的脂溶性輔酶促進泛醌及其衍生物的合成。

表3 灌木亞菊藥用活性成分代謝通路分析Table 3 Metabolic pathway of active components of A.fruticulose

2.5 灌木亞菊轉(zhuǎn)錄組的SSR及SNP分析

灌木亞菊轉(zhuǎn)錄組共鑒定出6種SSR重復類型，不同類型簡單重復序列出現(xiàn)的數(shù)量差距顯著，其中單、二和三堿基重復的數(shù)量要遠高于其他類型。單核苷酸多態(tài)性分析共獲得42 013個SNP位點，4 470個InDels，轉(zhuǎn)換和顛換所發(fā)生的比率近似2：1（圖7）。

圖7 灌木亞菊Unigenes的SSR（A）和SNP突變類型（B）分布Fig.7 Distribution of SSR（A）and SNPmutation type（B）of A.fruticulose Unigenes

3 結論與討論

灌木亞菊廣泛分布于我國西北地區(qū)及俄羅斯中亞地區(qū)，被用于治療闌尾炎、肺結核、肺氣腫等疾?。?3］。灌木亞菊多作為提取倍半萜和酚類等藥用化學組分的主要原料［24-25］。為了解灌木亞菊藥用次級代謝產(chǎn)物的代謝途徑，本研究對灌木亞菊進行了轉(zhuǎn)錄組測序及分析，旨在挖掘藥用成分的相關功能基因，為亞菊屬植物的資源利用及保護提供基礎數(shù)據(jù)。植物基因組結構復雜和倍性變化等導致其基因組學研究進展緩慢。而轉(zhuǎn)錄組可以精確反映植物在特定時空狀態(tài)下的基因表達，能夠在不同的組織器官、發(fā)育階段和生理狀態(tài)下，發(fā)現(xiàn)植物基因的差異性表達［27-28］。

對灌木亞菊轉(zhuǎn)錄本進行測序、質(zhì)控、拼接共獲得53 760條Unigenes，N50長度為1 155 bp，平均長度為817.23 bp，與巴哈雀稗（Paspalum notatum）［29］、光皮樹（Swida wilsoniana）［30］、萬壽菊（Tagetes erec?ta）［31］、濱麥（Leymus mollis）［32］、碰碰香（Plectranthus tomentosus）［33］等物種的轉(zhuǎn)錄本在同一水平，表明本研究所得拼接序列Unigenes能夠滿足轉(zhuǎn)錄組生物信息學分析的要求。

約有1/3（19 994條）的灌木亞菊Unigenes在與CDD等主要數(shù)據(jù)庫比對中未獲得注釋。在云南松（Pinus yunnanensis）［34］、多 穗 柯（Lithocarpus polys?tachyrus）［35］、云錦杜鵑（Rhododendron fortune）［36］等植物的研究中也有類似結果。這可能由于其近緣屬植物的轉(zhuǎn)錄本測序研究相對較少，數(shù)據(jù)庫中沒有足夠相似的序列與之匹配，也可能是灌木亞菊本身的特異性較高所導致，這需要對其近緣類群開展進一步的研究。NR注釋中，多年生草本植物洋薊（Cy?nara cardunculusvar.scolymus）與灌木亞菊的相似序列最多，洋薊與灌木亞菊同為菊科植物，其親緣關系較近，且二者均分布于干旱的荒漠環(huán)境中［1］。與灌木亞菊相類似，洋薊的次生代謝產(chǎn)物主要有酚酸、類黃酮、倍半萜等，具備抗菌、抗氧化、抗癌等特性［37］。通過與GO數(shù)據(jù)庫進行功能分類注釋，獲得與灌木亞菊藥用成分相關的基因主要富集于甾體、黃酮類、苯丙素類和萜類化合物，涉及次級產(chǎn)物的生物合成和代謝調(diào)控過程，基因富集提供了這些化合物合成途徑的線索。

通過KEGG數(shù)據(jù)庫進行代謝通路分析得到與灌木亞菊藥用活性相關的代謝通路百余條，如萜類、黃酮類、苯丙素類等。這些代謝通路與灌木亞菊藥效，如抗炎［38］、抗氧化作用［39］、抗菌抗氧化［40-41］等密切相關［42］。蒼術（Atractylodeslancea）轉(zhuǎn)錄組分析發(fā)現(xiàn)萜類骨架和倍半萜主要經(jīng)甲羥戊酸途徑進行生物合成，并且在干旱脅迫下倍半萜的合成通路下調(diào)［43］。灌木亞菊的倍半萜合成通路的數(shù)目相比較少，可能是因為灌木亞菊對荒漠及荒漠草原干旱環(huán)境的適應從而導致倍半萜合成通路依賴性降低。紫菀屬植物（Aster spathulifolius）轉(zhuǎn)錄組分析發(fā)現(xiàn)該植物在抗脂肪生成、抗高血脂等方面的功效主要與其高表達的苯丙素途徑密切相關［44］，而本研究的分析結果與之一致，表明灌木亞菊植物在抗脂肪生成、抗高血脂等方面可能也有一定應用潛力。但是仍需要對灌木亞菊開展進一步的研究，尤其是次級產(chǎn)物的藥效方面，如抗菌消炎、止痛止血等，以確定灌木亞菊的藥物效能。

此外，在灌木亞菊轉(zhuǎn)錄本中鑒定的SSR要以單堿基、三堿基、二堿基重復序列較多，與懷菊（Chry?santhemum morifolium）［45］、藥用植物半楓荷（Semiliq?uidambar cathayensis）［46］、長梗秦艽（Gentiana walton?ii）和粗壯秦艽（Gentiana robusta）［47］的結果相近。分子標記能夠反映基因座的遺傳多態(tài)性，為繪制高密度的遺傳圖譜提供了可能，為今后灌木亞菊的資源保護及遺傳多樣性研究提供了一系列的有效工具。

本研究通過高通量測序獲得了三江源地區(qū)灌木亞菊轉(zhuǎn)錄組，預測得到53 760條Unigenes，通過生物信息學分析獲得了豐富的功能基因、代謝通路、分子標記等基礎數(shù)據(jù)，發(fā)現(xiàn)與藥用活性相關的基因主要富集在甾體、黃酮類、苯丙素類及萜類化合物。KEGG代謝通路分析中得到233條Unigenes與藥用活性相關的次生代謝通路，主要涉及具有鎮(zhèn)痛消炎的苯丙素類、萜類化合物等成分。此外還檢測出SSR位點4 586個及SNP多態(tài)性位點42 013個，為灌木亞菊的資源利用和保護、遺傳多樣性等研究提供了基礎數(shù)據(jù)。