高寒灌叢建群種高山繡線菊的低溫適應(yīng)轉(zhuǎn)錄組信息分析

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

（1.中國(guó)科學(xué)院西北高原生物研究所高原適應(yīng)與進(jìn)化重點(diǎn)實(shí)驗(yàn)室，青海西寧810001；2.中國(guó)科學(xué)院大學(xué)生命科學(xué)學(xué)院，北京100039）

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

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

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

1 材料與方法

1.1 材料

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

1.2 方法

1.2.1 高山繡線菊RNA提取與建庫(kù)測(cè)序

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

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

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

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

隨機(jī)從Clean 數(shù)據(jù)中抽取10 000 條序列，與多個(gè)數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，取evalue<=1e-10 并且相似度＞90%，coverage＞80%的比對(duì)結(jié)果作為后續(xù)分析基因功能注釋及分類的數(shù)據(jù)，所采用的數(shù)據(jù)庫(kù)有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等。通過(guò)與NR 數(shù)據(jù)庫(kù)的比對(duì)，可得到高山繡線菊轉(zhuǎn)錄本序列與相近物種的近似情況以及同源序列的功能信息。根據(jù)Unigenes 與Swissprot、TrEMBL 的注釋結(jié)果得到GO（Gene Ontology）功能注釋信息，統(tǒng)計(jì)分子功能、細(xì)胞成分、生物過(guò)程三大分類中注釋成 功 的Unigenes。利 用KAAS v2.1［37］得 到 轉(zhuǎn) 錄 本KEGG注釋信息。

1.2.3 SNP與SSR分析

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

2 結(jié)果與分析

2.1 轉(zhuǎn)錄本組裝

高山繡線菊轉(zhuǎn)錄組高通量測(cè)序共獲得51 340 802條Raw reads，經(jīng)過(guò)濾處理后，得到49 947 114 條Clean Reads，總長(zhǎng)為7 212 398 740 bp，GC 含量為49.37%，根據(jù)Q20（堿基質(zhì)量在20 以上的序列，占97.71%）和Q30（堿基質(zhì)量在30 以上的序列，占93.68%）信息統(tǒng)計(jì)結(jié)果，可以說(shuō)明測(cè)序所得序列的質(zhì)量滿足后續(xù)的轉(zhuǎn)錄組分析。經(jīng)轉(zhuǎn)錄本組裝，共獲得117 280 個(gè)Trancripts，53 892 個(gè)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的長(zhǎng)度分布圖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 個(gè)數(shù)據(jù)庫(kù)進(jìn)行比對(duì)（表2）。注釋到NR 數(shù)據(jù)庫(kù)的Unigene 最多，占39.49%，注釋到KEGG 數(shù)據(jù)庫(kù)的Unigene 最少，占6.45%。35 954 條Unigenes 注釋到至少一個(gè)數(shù)據(jù)庫(kù)中，占66.71%，2 310 條Unigenes 能注釋到所有數(shù)據(jù)庫(kù)中，占4.29%。此外，還有17 938 條Unigenes 并未注釋成功。注釋到NR 庫(kù)的Unigenes 有31 627 條，其中比對(duì)到薔薇科梅（Armeniaca mumeSieb.）物種的序列數(shù)最多（8 163 條），其次分別為桃（AmygdalusL.，6 720 條）、蘋果屬（MalusMill.，4 000 條）、梨屬（PyrusL.，3 544 條）、草莓屬（FragariaL.，858 條）等植物。此外，還有8 325 條序列（26.32%）比對(duì)到其他525 種植物，但每個(gè)物種所能比對(duì)上序列都較少，其原因可能與高山繡線菊近緣種中基因組數(shù)據(jù)庫(kù)匱乏有關(guān)。

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

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

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

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

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

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

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

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

3 討論

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

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

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

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

4 結(jié)論

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