韓 霜，余靜雅，韓 赟，張發(fā)起

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

高山繡線菊（Spiraea alpinaPall.）是青藏高原高寒灌叢主要建群種之一，隸屬于薔薇科（Rosaceae）繡線菊屬（Spiraea），主要生長在海拔2 000~4 000 米的向陽坡地或灌叢中［11］。該植物廣泛分布于整個橫斷山區(qū)，生長于年平均氣溫較低的高寒甸，具耐寒、耐旱、耐瘠薄、耐陰濕等特點［12-14］。高山繡線菊可用于治療咽喉腫痛、風熱癢癥，是一種廣泛應(yīng)用于民間的中草藥，其根、葉、果實可做獸藥［15］。目前，對高山繡線菊的研究主要集中在化學成分、生物活性、繁殖栽培及冰期演化歷史上［16-20］，學者們證實從該植物中分離純化的抗真菌化合物對植物病原真菌具有一定的抑制作用［17］；作為高海拔地區(qū)園林綠化的少數(shù)樹種之一［16］，不少學者對該植物的馴化和栽培技術(shù)進行探索，指出其育苗技術(shù)簡單，栽培后成活率高且適應(yīng)性強［21-22］。此外，一些學者依據(jù)高山繡性菊葉片膜透性及膜脂過氧化特征，驗證其抗寒性強度［23］。然而關(guān)于高山繡線菊分子生物學的研究較少，尤其是基因組及轉(zhuǎn)錄組方面的研究。

隨著高通量測序技術(shù)的發(fā)展，使得轉(zhuǎn)錄組測序成為挖掘功能基因、篩選分子標記、闡明代謝途徑的有效工具［24］。目前，已成功應(yīng)用到該技術(shù)的模式生物有水稻（Oryza sativa）［25］、玉米（Zea mays）［26］、擬南芥（Arabidopsis thaliana）［27］等。在現(xiàn)有的植物轉(zhuǎn)錄組分析研究報道中，對灌木沙棘（Hippophae rhamnoides）進行轉(zhuǎn)錄組測序并與擬南芥比較后，最終確定與沙棘脂肪酸生物合成相關(guān)的基因序列［28］；三江源地區(qū)灌木亞菊的轉(zhuǎn)錄組解析成功獲取了藥用活性成分相關(guān)的代謝通路，為其資源利用和保護、遺傳多樣性提供基礎(chǔ)數(shù)據(jù)［29］。此外，利用轉(zhuǎn)錄組測序進行植物低溫脅迫相關(guān)的研究不少，如模式生物水稻轉(zhuǎn)錄組信息分析明確了其苗期植物激素對低溫的應(yīng)答機制［30］；驟然低溫下草本植物川百合（Lilium davidii）的轉(zhuǎn)錄組分析成功篩選出與抗寒性、光合作用、代謝途徑等相關(guān)關(guān)鍵基因，為川百合的分子育種提供參考［31］；木本植物仁用杏花（Prunus armeniaca）的轉(zhuǎn)錄組分析發(fā)現(xiàn)了抗寒差異表達基因，為今后采取基因工程手段培育抗寒的仁用杏新品種提供基因資源［32］；這些成功的案例為今后高山繡線菊適應(yīng)脅迫的分子機制研究提供參考。

高山繡線菊資源豐富，適應(yīng)性強［14］。目前，對該植物的研究相對較少，尤其是非生物逆境脅迫方面的認識較淺薄。因此，我們對高山繡線菊進行高通量測序并獲得轉(zhuǎn)錄組數(shù)據(jù)，對基因進行功能注釋并分析低溫脅迫相關(guān)代謝通路，挖掘關(guān)鍵基因，對后續(xù)其低溫脅迫分子機制的研究具有一定指導意義。

1 材料與方法

1.1 材料

高山繡線菊新鮮幼葉于秋季（9 月份）采集于青海省黃南州同仁縣采集（地理坐標：35°13′58.30″N，101°51′05.50″E；海拔：3 532米；年平均氣溫：5.2℃），用超純水和75%酒精清洗后迅速置于液氮中，后轉(zhuǎn)移至-80℃的超低溫冰箱中保存?zhèn)溆?。憑證標本（標本號：Zhang2018047）存于中國科學院西北高原生物研究所青藏高原生物標本館（HNWP）。

1.2 方法

1.2.1 高山繡線菊RNA提取與建庫測序

利用Total RNA Extractor（Trizol）試劑法［33］提取高山繡線菊總RNA，并用瓊脂糖凝膠電泳和NanoDrop-TM2000c（Thermo Scientific，美 國）檢測RNA 質(zhì)量純度。高山繡線菊轉(zhuǎn)錄組測序文庫的構(gòu)建參照夏銘澤在多裂駱駝蓬中的建庫方法［34］。檢測合格的文庫用Illumina HiseqTM進行測序。

1.2.2 轉(zhuǎn)錄組拼接與Unigene功能注釋

利用FastQC v0.11.2（http：//www.bioinformatics.babraham.ac.uk/projects/fastqc/）對測序數(shù)據(jù)進行質(zhì)量評估和質(zhì)控，過濾原始數(shù)據(jù)（Raw reads）中含有帶接頭的、低質(zhì)量的序列。通過Trimmomatic v0.36［35］進行質(zhì)量剪切，得到Clean reads。使用Trinityv2.4.0［36］（參數(shù)min_kmer cov=2，其余為默認設(shè)置）對樣本有效數(shù)據(jù)進行混合拼接，Trinity 組裝后的轉(zhuǎn)錄本序列信息以FASTA 格式儲存，并對轉(zhuǎn)錄本Transcript 去冗余，取每條基因中最長的轉(zhuǎn)錄本作為Unigene，以此作為后續(xù)分析的參考序列。

不進則退，李高明在蒙自花了三萬塊，挨著一家生意紅火的理發(fā)店，街頭街尾地打起了“價格戰(zhàn)”。此時的李高明，是小店五六個師傅中技術(shù)最好的，凡事他都想親力親為。有時生意好，給客人做頭做到一兩點，客人看晚了回不去，竟也樂意在他的小店里睡一宿。生意不好時，他犯愁，整夜整夜睡不著。心情起起伏伏地過了三個月，實在太煎熬，他受不了了。

隨機從Clean 數(shù)據(jù)中抽取10 000 條序列，與多個數(shù)據(jù)庫進行比對，取evalue<=1e-10 并且相似度＞90%，coverage＞80%的比對結(jié)果作為后續(xù)分析基因功能注釋及分類的數(shù)據(jù)，所采用的數(shù)據(jù)庫有CDD（Conserved Domain Database）、KOG（EuKaryotic Ortholog Groups）、COG（Clusters of Orthologous Groups of Proteins）、NR（NCBI Non-redundant Protein Sequences）、NT（NCBI Nucleotide Sequences）、PFAM（Protein Family）、Swissprot（A Manually Annotated and Reviewed Protein Sequence Database）、TrEMBL等。通過與NR 數(shù)據(jù)庫的比對，可得到高山繡線菊轉(zhuǎn)錄本序列與相近物種的近似情況以及同源序列的功能信息。根據(jù)Unigenes 與Swissprot、TrEMBL 的注釋結(jié)果得到GO（Gene Ontology）功能注釋信息，統(tǒng)計分子功能、細胞成分、生物過程三大分類中注釋成 功 的Unigenes。利 用KAAS v2.1［37］得 到 轉(zhuǎn) 錄 本KEGG注釋信息。

1.2.3 SNP與SSR分析

利用BCFtools v1.5［38］（參數(shù)：質(zhì)量值大于20且覆蓋度大于8）將組裝好的Unigene 作為參考序列進行單核苷酸變異（Single nucleotide polymorphisms，SNP）分析，找出可能的單核苷酸變異位點并統(tǒng)計篩選出SNP 突變類型。利用MISA v1.0［39］進行微衛(wèi)星（Microsatellite）標記或簡單序列重復（Simple sequence repeat）標記分析，分別設(shè)置二核苷酸、三核苷酸、四核苷酸、五核苷酸以及六核苷酸重復單元的重復次數(shù)為至少6、5、5、5、5次。

2 結(jié)果與分析

2.1 轉(zhuǎn)錄本組裝

高山繡線菊轉(zhuǎn)錄組高通量測序共獲得51 340 802條Raw reads，經(jīng)過濾處理后，得到49 947 114 條Clean Reads，總長為7 212 398 740 bp，GC 含量為49.37%，根據(jù)Q20（堿基質(zhì)量在20 以上的序列，占97.71%）和Q30（堿基質(zhì)量在30 以上的序列，占93.68%）信息統(tǒng)計結(jié)果，可以說明測序所得序列的質(zhì)量滿足后續(xù)的轉(zhuǎn)錄組分析。經(jīng)轉(zhuǎn)錄本組裝，共獲得117 280 個Trancripts，53 892 個Unigenes（表1）。其中200~300 bp 的Transcript 和Unigene 數(shù)量最多，1 900~2 000 bp的最少（圖1）。

表1 高山繡線菊轉(zhuǎn)錄組拼接結(jié)果Tab.1 Splicing statistical of S.alpina transcriptome 單位：bp

圖1 高山繡線菊Transcript與Unigenes的長度分布圖Fig.1 Distribution of Trancript and Unigenes length for S.alpina

2.2 Unigenes的功能注釋

利用NCBI Blast+［40］將Unigenes 與CDD、KOG、COG、NR、NT、PFAM、Swissprot、TrEMBL、GO 和KEGG 等9 個數(shù)據(jù)庫進行比對（表2）。注釋到NR 數(shù)據(jù)庫的Unigene 最多，占39.49%，注釋到KEGG 數(shù)據(jù)庫的Unigene 最少，占6.45%。35 954 條Unigenes 注釋到至少一個數(shù)據(jù)庫中，占66.71%，2 310 條Unigenes 能注釋到所有數(shù)據(jù)庫中，占4.29%。此外，還有17 938 條Unigenes 并未注釋成功。注釋到NR 庫的Unigenes 有31 627 條，其中比對到薔薇科梅（Armeniaca mumeSieb.）物種的序列數(shù)最多（8 163 條），其次分別為桃（AmygdalusL.，6 720 條）、蘋果屬（MalusMill.，4 000 條）、梨屬（PyrusL.，3 544 條）、草莓屬（FragariaL.，858 條）等植物。此外，還有8 325 條序列（26.32%）比對到其他525 種植物，但每個物種所能比對上序列都較少，其原因可能與高山繡線菊近緣種中基因組數(shù)據(jù)庫匱乏有關(guān)。

GO 可以全面描述生物體中基因和基因產(chǎn)物的屬性［41］。Unigene 與Swissprot、TrEMBL 比對后，有29 111 條Unigenes 獲得223 924 條注釋信息。根據(jù)注釋結(jié)果對得到的基因進行GO 分類（圖2），注釋到分子功能（Molecular function）、所處的細胞位置（Cellular component）和參與的生物過程（Biological process）三個ontology 的term（GO 分類的單位），分別有19、26、22 個子類，合計67 個子類，其中注釋在細胞（Cell）、細胞部分（Cell part）、細胞過程（Cellular process）子 類 的Unigene 較 多，分 別 有21 176 條（72.7%）、21 131 條（72.6%）、18 058 條（62%），而在電子載體活動（Chemoattractant activity）、形態(tài)發(fā)生活性（Morphogen activity）、生物（Biological phase）等子類中注釋到的Unigene 相對最少，分別有4 條（0.014%）、1條（72.7%）、1條（0.0034%）。

Unigene 與KOG 數(shù)據(jù)庫比對后，16 166 條基因被注釋成功，按KOG 的group 可分為26 個類型（表3）。其中，翻譯后修飾、蛋白轉(zhuǎn)運、信號傳遞機制和只有一般功能預測分類下的基因較多，分別有1 962條（12.1%）、1 953 條（12.1%）和1 850 條（11.4%），而細胞外結(jié)構(gòu)、核結(jié)構(gòu)和細胞活性注釋的基因數(shù)量較少，分別有60 條（0.33%）、48 條（0.27%）和10 條（0.056%）。此外，還有779 條注釋成功的基因未知其功能。

根據(jù)KO 與Pathway 的關(guān)聯(lián)性對其進行KEGG代謝通路分類。根據(jù)與KEGG 數(shù)據(jù)庫的比對結(jié)果，成功注釋到3 475 條Unigene，占6.45%。3 475 條基因被分為了代謝（Metabolism）、遺傳信息處理（Genetic information processing）、細胞過程（Cellular processes）和環(huán)境信息處理（Environmental information processing）四大類23 個子類。在四個大類中，涉及基因最多的是代謝，共有2 367 條，占68.1%，其次為遺傳信息處理、細胞過程和環(huán)境信息處理，分別有1 450 條（41.7%）、606 條（17.4%）和539 條（15.5%）。在23 個子類中，與代謝相關(guān)的通路最多，有12 條，與細胞過程、環(huán)境信息處理和遺傳信息處理三大類涉及的通路較少，分別有4 條、3 條和4 條。代謝途徑中，注釋基因最多的通路是翻譯（Translation）、信號轉(zhuǎn)導（Signal transduction），碳水化合物代謝（Carbohydrate metabolism），分別為682 條（19.6%）、512條（14.7%）和411 條（11.8%）；而注釋基因最少的通路是細胞運動（Cell motility）、膜運輸（Membrane transport）、信號分子和互作作用（Signaling molecules and interaction），分別為40 條（0.81%）、26 條（0.52%）和1條（0.02%）。

基于KEGG 數(shù)據(jù)庫分析，共統(tǒng)計到213 代謝途徑，其涉及6 560 條Unigenes，按注釋基因數(shù)量從高到低排序，選擇前11個代謝通路進行分析并列于表4 中。注釋基因最多的代謝通路為核糖體（Ribosome），有379 條，占總數(shù)的5.78%，其次為碳代謝（Carbon metabolism）和氨基酸合成（Biosynthesis of amino acid），分 別 為164 條（2.5%）和155 條（2.36%）。

表4 高山繡線菊Unigene數(shù)量最多的11個代謝通路Tab.4 The top eleven metabolic pathways involved in the largest number of S. alpina Unigenes

2.3 高山繡線菊低溫脅迫代謝通路分析

在213 條代謝途徑中，對高山繡線菊低溫脅迫代謝通路進行統(tǒng)計及分析，可分為低溫脅迫生理代謝（Physiological metabolism of cold resistance）、冷調(diào)節(jié)信號通路（Cold regulation signal pathway）、光合作用（Photosynthesis）等主要途徑（表5），分別有14 條（碳代謝、氨基酸類的生物合成、淀粉和蔗糖代謝為主）、4 條（光合作用、光合作用生物中的碳固定為主）和3 條（植物激素信號轉(zhuǎn)導為主）。所有代謝通路中，碳代謝（Carbon metabolism）、氨基酸類的生物合成（Biosynthesis of amino acids）、植物激素信號轉(zhuǎn)導（Plant hormone signal transduction）、淀粉和蔗糖代謝（Starch and sucrose metabolism）代謝通路涉及的Unigenes最多，分別為164條、155條、137條、101條；而亞油酸代謝（Linoleic acid metabolism）、甜菜堿生物合成（Betalain biosynthesis）代謝通路涉及的最少，分別為19條和13條。

表5 高山繡線菊低溫脅迫代謝通路及基因統(tǒng)計Tab.5 Metabolic pathway and gene statistical table of active components for S.alpina

3 討論

本研究利用Illumina Hiseq 對高山繡線菊轉(zhuǎn)錄組進行高通量測序，共獲得49 947 114 條Clean Reads。經(jīng)組裝共獲得117 280個Transcripts和53 892個Unigenes，其平均長度為708.72 bp，N50 為1 340 bp。相比于其他植物轉(zhuǎn)錄組測序及組裝結(jié)果［42-44］，如杜鵑（Rhododendron simsii；平均長度為636 nt，N50為1 018 nt）、羅布麻（Apocynum venetum；平均長度為878 bp，N50 為1 663 bp）和西藏嵩草（Kobresia tibetica；平均長度890.1 bp，N50為1 342 bp），根據(jù)有效數(shù)據(jù)Q20 值、Q30 值，表明高山繡線菊轉(zhuǎn)錄組測序及組裝質(zhì)量較好，能夠滿足后續(xù)轉(zhuǎn)錄組信息分析的要求。

與9 個數(shù)據(jù)庫比對后，成功注釋的Unigenes 數(shù)量為53 892 條，而所得結(jié)果中仍然還有17 938 條Unigenes 未與已知基因匹配成功。這在其它植物的分析中也有發(fā)現(xiàn)，如珠子參（Panax japonicusvar.Major）［45］、虎杖根（Polygonum cuspidatum）［46］、金錢松（Pseudolarix amabilis）［47］等，這可能與基因庫中缺少該種相關(guān)的基因組信息有關(guān)；其次，缺乏轉(zhuǎn)錄組方面的研究作為參考，導致該種特有的某些基因和數(shù)據(jù)庫中序列的識別和比對十分困難。KOG 結(jié)果幾乎整合了高山繡線菊所有信息，為其基因功能研究提供了數(shù)據(jù)基礎(chǔ)，成功注釋到16 166條Unigenes，其中參與該物種翻譯后修飾、蛋白轉(zhuǎn)運和信號傳遞機制等生命活動通路的基因最多，可以推斷這些基因在該物種中表達較豐富，說明這三類生命活動在該物種的生長發(fā)育中具有重要的地位。NR 數(shù)據(jù)庫中比對到梅的序列數(shù)最多，表明二者具有較高的序列同源性，親緣關(guān)系較近，但僅有5個同科植物與該物種比對成功，其原因可能與薔薇科植物的轉(zhuǎn)錄組及基因組信息較匱乏或該種本身具有的特異基因較多有關(guān)。

在低溫脅迫下，植物通過改變生理生化、響應(yīng)抗寒分子機制等方式提高自身抗寒性。具體而言，植物細胞膜系統(tǒng)、光合作用、可溶性糖/蛋白含量、游離脯氨酸含量等均受到影響而發(fā)生變化，植物也會對低溫作出響應(yīng)并進行一系列的冷信號轉(zhuǎn)導途徑，如碳水化合物（蔗糖和淀粉）代謝，植物激素合成及轉(zhuǎn)導，次生代謝產(chǎn)物合成，信號轉(zhuǎn)導（Ca2+），光合系統(tǒng)，脂質(zhì)代謝等［10，48］。本研究依據(jù)KEGG 數(shù)據(jù)庫對低溫脅迫相關(guān)代謝通路的篩選結(jié)果與青海草地早熟禾一致［49］。將Unigene 映射到與低溫脅迫相關(guān)的代謝通路中，對其進行初步篩選，結(jié)果顯示低溫脅迫生理代謝涉及的通路占大多數(shù)，碳代謝、氨基酸類的生物合成等通路所涉及的Unigene 最多，進一步說明這些代謝通路與高山繡線菊的耐寒性特征密切相關(guān)。

此外，低溫脅迫下的植物是基于基因表達的迅速變化而進行的轉(zhuǎn)錄調(diào)控，從而合成相應(yīng)蛋白指導代謝物的生成，以此對脅迫做出響應(yīng)，如同科植物蘋果的轉(zhuǎn)錄組分析發(fā)現(xiàn)植物激素信號轉(zhuǎn)導、光合作用、糖酵解、淀粉和蔗糖等代謝通路在低溫處理后發(fā)生基因的差異表達，前者通路中差異基因上調(diào)的數(shù)目高于下調(diào)的數(shù)目，而后兩者通路中則呈現(xiàn)相反的情況［48，50］。與蘋果轉(zhuǎn)錄組相比，高山繡線菊映射到這些代謝通路的Unigene 數(shù)量較多，這可能與該植物在低溫環(huán)境下對這些代謝通路的依賴性較高有關(guān)，進一步說明了兩種植物受到的脅迫程度和響應(yīng)脅迫的基因數(shù)量具有差異，而這些代謝通路在高山繡線菊響應(yīng)低溫脅迫時發(fā)揮著重要作用，這為低溫脅迫相關(guān)關(guān)鍵差異表達基因的尋找提供了線索。

4 結(jié)論

本研究利用分子手段獲得了高山繡線菊的轉(zhuǎn)錄組數(shù)據(jù)，在無參考基因組的前提下，對這些數(shù)據(jù)進行分析，包括基因比對、功能注釋、代謝通路，為其后續(xù)的分子生物學研究提供了基礎(chǔ)數(shù)據(jù)，進一步豐富了薔薇科植物轉(zhuǎn)錄組數(shù)據(jù)庫，同時篩選出與低溫脅迫相關(guān)的代謝通路，挖掘關(guān)鍵基因，對后續(xù)其適應(yīng)低溫脅迫分子機制的研究具有一定指導意義。