青藏高原地區(qū)3種蒿屬植物轉(zhuǎn)錄組比較分析

余靜雅 夏銘澤 徐 浩 張發(fā)起

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

蒿屬（）隸屬于菊科（Asteraceae），本屬物種數(shù)目在不同的文獻記錄中出入較大，一般在350～500 種，多為草本和半灌木，主要分布在北半球溫暖干旱的地區(qū)。米蒿（-Krasch.）屬于蒿亞屬（Subgen.）艾蒿組（Sect.），為半灌木狀草本，主要生長于內(nèi)蒙古、甘肅、青海和西藏等地區(qū)。冷蒿（Willd.）屬于蒿亞屬蒔蘿蒿組（Sect.），為多年生草本，有時略成半灌木狀，廣泛分布于干旱和半干旱地區(qū)；臭蒿（Ostenf.et Pauls.）屬于蒿亞屬艾蒿組，為一年生草本，具有濃烈臭味，從低海拔到高海拔地區(qū)均有分布。蒿屬植物的揮發(fā)油具有殺菌抑菌，防治蚊蟲、消滅害蟲的功效。此外，蒿屬植物常常作為荒漠化治理的先鋒物種，米蒿具有較強的阻沙能力，在自然恢復沙地的草本植物中占有較大的生態(tài)位寬度，生態(tài)適應能力較強；冷蒿對干旱和強風脅迫都有良好的抗性，能夠在退化草場形成優(yōu)勢種群，阻止草原沙漠化進一步發(fā)展。但是，冷蒿能夠向周圍環(huán)境釋放化感物質(zhì)，其提取液對多種牧草均有抑制生長的作用，被認為可能是草原退化原因之一。

蒿屬中多種植物作為傳統(tǒng)中草藥，具有清熱涼血，退虛熱，解暑，驅(qū)寒除濕，活血通經(jīng)等功效，主要用于治療瘧疾、發(fā)熱感冒、自身免疫疾病等。臭蒿作為常用藏藥，其化學成分具有良好的抗菌消炎活性，可用于治療“赤巴”病、急性黃疸性肝炎、膽囊炎等。冷蒿作為蒙藥“阿給”，全草入藥，主要用于治療各種出血和關節(jié)腫脹。蒿屬植物的化學成分主要為萜類、黃酮類、香豆素類、咖啡酰奎寧酸、甾醇類和炔類化合物。由于蒿屬植物的生物化學和藥理功能的多樣性，以及獨特的生物活性成分合成方式，使其在藥物開發(fā)領域受到越來越多的關注。主產(chǎn)青蒿素的黃花蒿（L.）已經(jīng)發(fā)表了全基因組序列，蒿屬植物的轉(zhuǎn)錄組和葉綠體基因組也已公布。目前對于蒿屬的大多數(shù)研究集中在蒿屬植物的藥理學和生態(tài)功能方面，但是大部分相關基因還處于未知狀態(tài)，研究進展緩慢。

本研究以3種蒿屬植物（米蒿、冷蒿和臭蒿）的成熟葉片為材料，采用高通量測序技術獲取轉(zhuǎn)錄組數(shù)據(jù)并組裝，并通過分析基因結(jié)構、基因表達水平和差異表達基因，篩選高表達的差異基因和活性成分相關基因，以期為蒿屬植物的物種鑒定和資源利用提供一定的分子依據(jù)。

1 材料與方法

1.1 實驗材料

米蒿采集于青海省海南州共和縣鐵蓋鄉(xiāng)（36°06′56.61″N，100°41′41.74″E，海拔2 594 m），冷蒿和臭蒿采集于青海省海南州共和縣龍羊峽（36°10′53.90″N，100°59′41.35″E，海拔2 837 m），在野外選擇生長狀況良好的植株（每個物種選擇一株）進行葉片采集，將采集到的3 種蒿屬植物花期的葉片放入液氮速凍后保存于-80 ℃。憑證標本米蒿（zhang2018014）、冷蒿（zhang2018019）和臭蒿（zhang2018028）存放于中國科學院西北高原生物研究所青藏高原生物標本館（HNWP）。

1.2 實驗方法

1.2.1 RNA提取

將樣品于液氮中研磨，利用TRIzol Reagent（Invitrogen，美國）提取總RNA，用瓊脂糖凝膠電泳檢驗RNA的降解程度以及是否有污染，使用Nano‐drop 2000 初步檢測RNA 的純度（A/A應在1.9～2.0），最后使用Agilent 2100 測定樣品的RIN 值以確定RNA完整性。

1.2.2 RNAseq建庫測序及結(jié)果評估

利用Illumina Hiseq2500 高通量測序平臺進行測序，得到的原始圖像文件由CASAVA 堿基識別（Base Calling）后分析轉(zhuǎn)化為原始測序數(shù)據(jù)（Raw reads）。該數(shù)據(jù)文件包含測序得到的堿基序列（reads）信息及其對應的堿基測序質(zhì)量信息。使用FastQC v0.11.2（http：//www.bioinformatics.babra‐ham.ac.uk/projects/fastqc/）對測序的原始數(shù)據(jù)質(zhì)量進行可視化評估。原始測序數(shù)據(jù)含有帶接頭的、低質(zhì)量的序列，利用Trimmomatic v0.36對原始數(shù)據(jù)進行質(zhì)量剪切得到過濾后的序列（Clean reads）。利用HISAT2 v2.1.0將Clean reads 比對到已發(fā)表的黃花蒿基因組上，比對結(jié)果使用StringTie v1.3.3b進行組裝。

1.2.3 表達水平和基因結(jié)構分析

使用StringTie 1.3.3b和已知的基因模型評估基因表達量。使用BCFtools v1.5根據(jù)比對結(jié)果進行單核苷酸變異（SNP）分析，過濾條件為：質(zhì)量值大于20 且覆蓋度大于8，尋找可能的SNP 位點。通過SnpEff v2.36統(tǒng)計變異位點在基因組結(jié)構上的分布情況。使用ASprofile v1.04根據(jù)每個樣本預測出來的基因模型對可變剪切（Alterna‐tive Splicing，AS）進行分類。

1.2.4 表達差異基因分析

使用DESeq2進行基因表達差異分析，篩選差異表達基因（Differentially Expressed Genes，DEGs），參數(shù)為：值<0.05且差異倍數(shù)Fold Change>2。

1.2.5 差異基因富集分析

利 用topGO v2.24.0（http：//bioconductor.org/biocLite.R）對差異基因進行GO（Gene Ontology）富集 分 析，使 用clusterProfiler v3.0.5進 行KEGG（Kyoto Encyclopedia of Genes and Genomes）path‐way富集分析。

2 結(jié)果與分析

2.1 轉(zhuǎn)錄組測序結(jié)果評估

對米蒿（ADL）、冷蒿（AF）和臭蒿（AH）的植物葉片進行Illumina高通量測序，將得到的raw data過濾后分別獲得54 268 322，46 434 864和43 971 646條Clean reads。每條reads平均長度分別為143.37，144.70和143.68 bp，Q30 分別為94.63%、95.08%和94.82%。通過參考序列比對分析，3 種蒿屬植物分別有61.27%、58.76%和70.16%的Clean reads 比對到參考基因組上（見表1）。GC 含量分別為46.70%、46.30%和45.89%。

表1 測序結(jié)果統(tǒng)計與質(zhì)量評估Table 1 Statistics and quality assessment of sequencing data

2.2 SNP分析

在ADL、AF和AH轉(zhuǎn)錄組中分別找到了378 757、327 045和298 104個可能的SNP位點，5 259、2 558和2 822 個InDel。轉(zhuǎn)錄突變譜系顯示在3 個樣品的轉(zhuǎn)錄組中的SNP分布不均勻（見圖1）。3個樣品中，轉(zhuǎn)換突變SNP的數(shù)量分別為247 811、212 161和192 855 個（按數(shù)量從多到少分別為C→T、G→A、T→C、A→G），顛換突變SNP的數(shù)量分別為133 102、117 465 和105 563 個（按數(shù)量由多到少分別為T→A、A→T、C→A、G→T、A→C、T→G、G→C、C→G）。并且在3 個樣品中轉(zhuǎn)換突變類型的SNP 數(shù)量約為顛換突變類型的2倍。

圖1 轉(zhuǎn)錄組SNP突變類型分布圖Fig.1 Mutation pedigree of transcriptome SNP

2.3 可變剪切分析

可變剪切事件一般分為五類：外顯子替換、內(nèi)含子保留、可變外顯子、可變轉(zhuǎn)錄起始位點和可變轉(zhuǎn)錄終止位點。通過ASprofile 軟件，將可變剪接事件細分為12 類，以ADL 為例發(fā)生概率從大到小分別為：第一個外顯子可變剪切（alternative tran‐scription start site，TSS）、最后一個外顯子可變剪切（alternative transcription termination site，TTS）、單外顯子跳躍（模糊邊界）（approximate exon skip‐ping，XSKIP）、單內(nèi)含子保留（retention of single in‐tron，IR）、外顯子替換（alternative exon ends，AE）、外顯子替換（模糊邊界）（approximate alternative ex‐on ends，XAE）、外顯子跳躍（exon skipping，SKIP）、多外顯子跳躍（cassette exons，MSKIP）、單內(nèi)含子保留（模糊邊界）（approximate retention of single in‐tron，XIR）、多外顯子跳躍（模糊邊界）（approximate cassette exons，XMSKIP）、多內(nèi)含子保留（retention of multiple introns，MIR）和多內(nèi)含子保留（模糊邊界）（approximate retention of multiple introns，XMIR）（見表2）?？勺兗艚宇愋偷臄?shù)量模式在3個樣品中基本一致，可變剪切事件數(shù)量為ADL>AF>AH，TSS和TTS類型的可變剪切數(shù)量明顯高于其他類型的可變剪切。

表2 可變剪切事件統(tǒng)計表Table 2 Statistics of differentially happened AS events

2.4 差異基因表達分析

對3 種蒿屬植物的基因表達量相關性分析表明，AH 和ADL 最為相近，聚為一類（見圖2A）。以ADL、AF 和AH作相互對比，ADL vs AH、ADL vs AF和AF vs AH 的上調(diào)表達基因數(shù)量分別為4 822、6 107和3 755，下調(diào)表達的基因數(shù)量分別為3 369、3 133 和5 079。除了AF vs AH，在其余兩個對比組中上調(diào)表達基因均大于下調(diào)表達基因（見圖2B）。經(jīng)過篩選，一共鑒定出13 755個DEGs，其中1 593 個（11.58%）DEGs 為三者共有，3 135 個（22.94%）DEGs 為ADL vs AF 和ADL vs AH 共有，3 501 個（25.45%）DEGs 為ADL vs AF 和AF vs AH共有，2 688 個（19.54%）DEGs 為ADL vs AH 和AF vs AH 共 有，ADL vs AF 特 有 的DEGs 為1 011 個（7.35%），ADL vs AH 特有的DEGs 最少，為775 個（5.63%），AF vs AH 特有的DEGs 最多，為1 052 個（7.65%）（見圖2C）。

圖2 轉(zhuǎn)錄組測序差異表達基因A.樣品間差異基因表達量聚類熱圖；B.樣品間上下調(diào)表達基因數(shù)目分布情況；C.差異表達基因的韋恩圖Fig.2 Differentially expressed genes obtained of RNA-sequencingA.The cluster heat map of differentially expressed genes among samples；B.Number of up-and down-regulated genes between samples；C.Venn dia‐gram of differentially expressed genes among ADL vs AF，ADL vs AH，AF vs AH sample set

2.5 GO功能富集分析

為了分析DEGs 的功能，對其進行GO 分類和KEGG 富集分析。GO 分類結(jié)果（見圖3）表明，在ADL vs AH、ADL vs AF 和AF vs AH 這3組對比中，差異基因富集存在共同的分類特征，都分別注釋到了GO 分類生物過程（biological process）的21 個條目、細胞組分（cellular component）的12個條目和分子功能（molecular function）的14 個條目。在生物過程分類中，3組對比中差異基因富集最多的兩個節(jié)點分別為代謝過程（metabolic process）和細胞過程（cellular process），其中ADL vs AF 和ADL vs AH對比組在這兩個代謝過程中主要富集上調(diào)基因。在細胞組分分類中，3個對比組中差異基因富集最多的兩個節(jié)點分別為細胞（cell）和細胞部分（cell part）。在AF vs AH 對比組中，胞外部分（extracel‐lular region）和細胞外基質(zhì)（extracellular matrix）都主要富集下調(diào)基因，這點與其他2組有所不同。在分子功能分類中，ADL vs AH 和AF vs AH 對比中差異基因富集最多的兩個節(jié)點分別為酶調(diào)節(jié)活性（enzyme regulator activity）和金屬伴侶活性（elec‐tron transfer activity）。ADL vs AF 對比組在金屬伴侶活性節(jié)點中主要富集上調(diào)基因，而在AF vs AH對比組在該節(jié)點主要富集下調(diào)基因。在ADL vs AF 對比組中，差異基因富集最多的兩個節(jié)點分別為催化活性（catalytic activity）和結(jié)合（binding）。在結(jié)合這一節(jié)點處，ADL vs AF 對比組主要富集上調(diào)基因，而ADL vs AH 和AF vs AH 對比組在該節(jié)點主要富集下調(diào)基因。

圖3 差異表達基因的GO富集分析Fig.3 Function classification of the DEGs using GO analysis.

2.6 KEGG功能富集分析

分別有1 003（ADL vs AH）、1 084（ADL vs AF）和1 195（AF vs AH）個DEGs注釋到KEGG數(shù)據(jù)庫。在ADL vs AF 對比組中，DEGs 富集到198 條通路，顯著富集程度前三的代謝通路依次為MAPK 信號通路（MAPK signaling pathway）、類黃酮生物合成（Flavonoid biosynthesis）和硫代謝（Sulfur metabo‐lism），得到注釋的DEGs最多的代謝通路是植物激素信號轉(zhuǎn)導（Plant hormone signal transduction），共注釋到75 個DEGs，和氨基酸的生物合成（Biosyn‐thesis of amino acids），共注釋到62 個DEGs（見圖4A）。在ADL vs AH 對比組中，DEGs富集到198條通路，顯著富集程度前三的代謝通路依次為到脂肪酸代謝（Fatty acid metabolism）、β-丙氨酸代謝（beta-Alanine metabolism）和脂肪酸降解（Fatty ac‐id degradation），得到注釋的DEGs最多的代謝通路是二氧化碳代謝（Carbon metabolism），共注釋到65個DEGs，和氨基酸的生物合成（Biosynthesis of amino acids），共注釋到52 個DEGs（見圖4B）。在AF vs AH 對 比 組 中，DEGs 富 集 到197 條 代 謝 通路，顯著富集程度前三的代謝通路依次為氨基酸的生物合成（Biosynthesis of amino acids）、α-亞麻酸鈣代謝（alpha-Linolenic acid metabolism）和牛磺酸和亞?；撬岽x（Taurine and hypotaurine metab‐olism），得到注釋的DEGs 最多的代謝通路是氨基酸的生物合成（Biosynthesis of amino acids），共注釋到78 個DEGs，和植物激素信號轉(zhuǎn)導（Plant hor‐mone signal transduction），共注釋到74 個DEGs（見圖4C）。

圖4 差異表達基因KEGG富集分析A.ADL vs AF的KEGG富集分析；B.ADL vs AH的KEGG富集分析；C.AF vs AH的KEGG富集分析Fig.4 KEGG enrichment analysis of differentially expressed gene（sDEGs）A.KEGG enrichment analysis of the DEGs in ADL vs AH；B.KEGG enrichment analysis of the DEGs in ADL vs AF；C.KEGG enrichment analysis of the DEGs in AF vs AH

2.7 高表達差異基因

目前關于蒿屬植物各基因的分子作用和機制未有系統(tǒng)性的研究，存在未知的功能基因，但對于未知的高表達差異基因也應該給予關注。以log2（FC）值作為篩選標準，選取各對比組內(nèi)log2（FC）值最高的10 個基因進項功能分析（見表3）。ADL vs AF 對比組中存在4 個未知的功能基因，6 個已知的功能基因主要與應激和氧化還原過程相關。ADL vs AH 對比組中存在6 個未知的功能基因，4個已知的功能基因主要與氧化還原過程和應激反應相關。AF vs AH 對比組中存在3 個未知的功能基因，7個已知的功能基因主要與氧化還原過程和轉(zhuǎn)運過程相關。CTI12_AA551550、CTI12_AA6194 70、CTI12_AA626460 和CTI12_AA097740 在ADL vs AF 和ADL vs AH 對比組中表達量相同。CTI12_AA106760 和CTI12_AA274180 在ADL vs AH對比組中表達量比在AF vs AH對比組稍高。

表3 高表達基因功能分析Table 3 Functional analysis of high expression genes

2.8 共表達活性成分相關基因

根據(jù)KEGG pathway 富集分析結(jié)果，共篩選到25 個在3 個對比組中共表達的活性成分相關基因，根據(jù)其功能主要分為3類，分別是生物堿、萜烯類和類黃酮（見表4）。生物堿相關表達基因為ALDH、E2.1.1.104、OGDH、sucA、GOT2 和TAT。萜烯類相關表達基因為DXS、VTE3、APG1、wrbA、E2.1.1.95、GERD、DHDDS、RER2、SRT1、TAT、E1.1.1.208 和CYP82G1。類黃酮相關表達基因為E2.1.1.104、CYP98A、C3′H、FLS和CHS。

表4 活性成分相關基因分類Table 4 Gene classification of active ingredients

3 討論

米蒿、冷蒿和臭蒿這3種蒿屬植物均為西北干旱、半干旱地區(qū)常見物種，在改善生態(tài)環(huán)境方面作用良多。并且作為傳統(tǒng)中藥材，具有良好的臨床療效和廣泛的應用前景。這3 種蒿屬植物在西北地區(qū)資源量較大，通過轉(zhuǎn)錄組比較分析，了解其序列信息和基因表達的差異情況，進一步明確差異基因的代謝途徑和生物功能，以挖掘關鍵的功能基因，為植物資源利用提供分子基礎。

通過對3 種蒿屬植物成熟葉片進行轉(zhuǎn)錄組測序分析后，發(fā)現(xiàn)不同蒿屬植物的轉(zhuǎn)錄本差異較大。通過差異表達基因分析，ADL vs AF 對比組擁有更多的差異基因，差異表達更明顯，含有最多的上調(diào)差異表達基因，說明米蒿和冷蒿差異較大。米蒿和臭蒿同屬艾蒿組，而冷蒿屬于蒔蘿蒿組，同時差異基因表達量聚類熱圖也顯示出米蒿和臭蒿關系較近，而與冷蒿關系較遠。在3 種蒿屬植物ADL、AF和AH轉(zhuǎn)錄組中分別找到了378 757、327 045和298 104 個可能的SNP 位點，轉(zhuǎn)換突變類型的SNP數(shù)量約為顛換突變類型的2倍，這種情況在其他植物轉(zhuǎn)錄組的研究中均有發(fā)現(xiàn)。SNP 位點可用于物種鑒定，開發(fā)分子標記，本研究發(fā)現(xiàn)的SNP 可以為蒿屬植物的物種鑒定、系統(tǒng)發(fā)育關系及物種親緣關系分析提供分子依據(jù)。

可變剪切在生物中普遍存在，是調(diào)控基因互作的重要分子機制，在植物生長發(fā)育和非生物脅迫響應中起到重要的作用：調(diào)節(jié)開花時間、保證生物鐘正常運轉(zhuǎn)、響應溫度脅迫和滲透脅迫?？勺兗羟心軌蛟鰪娭参飳Ω珊岛蜏囟让{迫的適應能力，激發(fā)特定的分子功能。本研究中，測序數(shù)據(jù)的大小不同導致3 種蒿屬植物可變剪切數(shù)量的不同，但是3個樣品的可變剪切事件數(shù)量的分布呈現(xiàn)出一定的相似性，說明在3種蒿屬植物中進行著較為一致的生理反應。本研究發(fā)現(xiàn)的可變剪接事件可為后續(xù)繼續(xù)探尋蒿屬植物相關生物反應的分子機制提供一定的生物信息學基礎。

3 個對照組的高表達的差異基因與防御反應和氧化還原反應緊密相關，防御反應相關基因的高表達可能與野外采摘樣品損傷葉片導致應激表達有關。氧化還原反應中高表達的C-4 甾醇甲基氧化酶活性在之前的研究中被證明與脫水脅迫密切相關，其相關基因可以作為培育耐旱植物的候選。同時，HSP20 家族蛋白在高表達已被證明參與高溫、鹽和干旱脅迫的反應，相關基因的高表達對應著3種蒿屬植物生活的炎熱干旱的環(huán)境，可為蒿屬植物耐熱機制研究提供相關的分子依據(jù)。

GO 富集分析結(jié)果顯示3 個對照組在代謝過程、細胞過程和細胞這3個條目共同富集大量差異基因，而在胞外部分、結(jié)合、催化活性這3個通路明顯不同。這說明3種蒿屬植物具有一定的差異，冷蒿體內(nèi)可能進行著更多的揮發(fā)性物質(zhì)的相關反應。KEGG 富集分析結(jié)果顯示，3 個對比組共同富集最多的通路為類黃酮代謝和甘油酯生物合成。類黃酮具有良好的抗氧化作用，可以提高植物對環(huán)境脅迫的抗性，研究證明干旱和炎熱會導致蒿屬植物中類黃酮的累積。甘油酯主要包括甘油三酯、甘油二酯、甘油一酯、糖脂、磷脂和甜菜堿脂，參與細胞膜脂的形成，是植物體內(nèi)重要的初生代謝產(chǎn)物，對植物生長發(fā)育和抵抗環(huán)境脅迫具有重要作用。研究表明，在鹽脅迫下，隨著鹽度的增加，雙對柵藻（）中甘油三酯的含量呈現(xiàn)先升后降的趨勢；在氮磷脅迫下，微藻中糖脂的累積量增加。本研究中的米蒿、冷蒿和臭蒿生長于青海省海西州的干旱鹽堿地，這可能是導致類黃酮代謝和甘油酯生物合成這兩個通路的共同富集的原因，推測相關基因與這蒿屬植物對環(huán)境脅迫的適應有關。

通過GO和KEGG富集分析，共篩選出25個相關的共表達差異基因。相比于其他共表達差異基因，CTI12_AA555350 基因表達量非常高，該基因?qū)儆贕ERD 家族，其被證明是響應昆蟲取食的關鍵基因，其產(chǎn)物大根香葉烯具有良好的驅(qū)蟲作用。ALDH 能夠使有毒的醛類氧化為無毒的羧酸，同時該基因能夠響應植物非生物脅迫，在黃花蒿的轉(zhuǎn)錄組研究中，ALDH基因參與青蒿素生物合成途徑。青蒿素是從黃花蒿中的一種內(nèi)過氧化物倍半萜內(nèi)酯（Endoperoxide sesquiterpene lac‐tone），在蒿屬的多種植物中均有存在，相關基因及通路已有詳細報道，然而其他蒿屬植物的活性成分研究較少。本研究篩選的與活性成分相關的差異基因可為后續(xù)蒿屬植物活性成分的研究提供參考信息。

本研究利用Illumina 測序技術對米蒿、冷蒿和臭蒿的成熟葉片進行了轉(zhuǎn)錄組測序，使用生物信息學的方法對測序結(jié)果進行分析。通過SNP、可變剪切、差異表達基因和功能富集分析，對3 種蒿屬植物的基因結(jié)構和代謝通路有了一定的了解，并篩選了活性成分相關的差異基因。研究結(jié)果為今后研究蒿屬植物生長發(fā)育、抗逆和活性物質(zhì)和植物資源的開發(fā)與利用提供一定的科學依據(jù)。