青藏高原地区3种蒿属植物转录组比较分析

余静雅 夏铭泽 徐 浩 张发起

（1. 中国科学院西北高原生物研究所高原生物适应与进化重点实验室，西宁 810001；2. 中国科学院大学生命科学学院，北京 100049）

蒿属（）隶属于菊科（Asteraceae），本属物种数目在不同的文献记录中出入较大，一般在350～500 种，多为草本和半灌木，主要分布在北半球温暖干旱的地区。米蒿（-Krasch.）属于蒿亚属（Subgen.）艾蒿组（Sect.），为半灌木状草本，主要生长于内蒙古、甘肃、青海和西藏等地区。冷蒿（Willd.）属于蒿亚属莳萝蒿组（Sect.），为多年生草本，有时略成半灌木状，广泛分布于干旱和半干旱地区；臭蒿（Ostenf.et Pauls.）属于蒿亚属艾蒿组，为一年生草本，具有浓烈臭味，从低海拔到高海拔地区均有分布。蒿属植物的挥发油具有杀菌抑菌，防治蚊虫、消灭害虫的功效。此外，蒿属植物常常作为荒漠化治理的先锋物种，米蒿具有较强的阻沙能力，在自然恢复沙地的草本植物中占有较大的生态位宽度，生态适应能力较强；冷蒿对干旱和强风胁迫都有良好的抗性，能够在退化草场形成优势种群，阻止草原沙漠化进一步发展。但是，冷蒿能够向周围环境释放化感物质，其提取液对多种牧草均有抑制生长的作用，被认为可能是草原退化原因之一。

蒿属中多种植物作为传统中草药，具有清热凉血，退虚热，解暑，驱寒除湿，活血通经等功效，主要用于治疗疟疾、发热感冒、自身免疫疾病等。臭蒿作为常用藏药，其化学成分具有良好的抗菌消炎活性，可用于治疗“赤巴”病、急性黄疸性肝炎、胆囊炎等。冷蒿作为蒙药“阿给”，全草入药，主要用于治疗各种出血和关节肿胀。蒿属植物的化学成分主要为萜类、黄酮类、香豆素类、咖啡酰奎宁酸、甾醇类和炔类化合物。由于蒿属植物的生物化学和药理功能的多样性，以及独特的生物活性成分合成方式，使其在药物开发领域受到越来越多的关注。主产青蒿素的黄花蒿（L.）已经发表了全基因组序列，蒿属植物的转录组和叶绿体基因组也已公布。目前对于蒿属的大多数研究集中在蒿属植物的药理学和生态功能方面，但是大部分相关基因还处于未知状态，研究进展缓慢。

本研究以3种蒿属植物（米蒿、冷蒿和臭蒿）的成熟叶片为材料，采用高通量测序技术获取转录组数据并组装，并通过分析基因结构、基因表达水平和差异表达基因，筛选高表达的差异基因和活性成分相关基因，以期为蒿属植物的物种鉴定和资源利用提供一定的分子依据。

1 材料与方法

1.1 实验材料

米蒿采集于青海省海南州共和县铁盖乡（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建库测序及结果评估

利用Illumina Hiseq2500 高通量测序平台进行测序，得到的原始图像文件由CASAVA 碱基识别（Base Calling）后分析转化为原始测序数据（Raw reads）。该数据文件包含测序得到的碱基序列（reads）信息及其对应的碱基测序质量信息。使用FastQC v0.11.2（http：//www.bioinformatics.babra‐ham.ac.uk/projects/fastqc/）对测序的原始数据质量进行可视化评估。原始测序数据含有带接头的、低质量的序列，利用Trimmomatic v0.36对原始数据进行质量剪切得到过滤后的序列（Clean reads）。利用HISAT2 v2.1.0将Clean reads 比对到已发表的黄花蒿基因组上，比对结果使用StringTie v1.3.3b进行组装。

1.2.3 表达水平和基因结构分析

使用StringTie 1.3.3b和已知的基因模型评估基因表达量。使用BCFtools v1.5根据比对结果进行单核苷酸变异（SNP）分析，过滤条件为：质量值大于20 且覆盖度大于8，寻找可能的SNP 位点。通过SnpEff v2.36统计变异位点在基因组结构上的分布情况。使用ASprofile v1.04根据每个样本预测出来的基因模型对可变剪切（Alterna‐tive Splicing，AS）进行分类。

1.2.4 表达差异基因分析

使用DESeq2进行基因表达差异分析，筛选差异表达基因（Differentially Expressed Genes，DEGs），参数为：值<0.05且差异倍数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 结果与分析

2.1 转录组测序结果评估

对米蒿（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 测序结果统计与质量评估Table 1 Statistics and quality assessment of sequencing data

2.2 SNP分析

在ADL、AF和AH转录组中分别找到了378 757、327 045和298 104个可能的SNP位点，5 259、2 558和2 822 个InDel。转录突变谱系显示在3 个样品的转录组中的SNP分布不均匀（见图1）。3个样品中，转换突变SNP的数量分别为247 811、212 161和192 855 个（按数量从多到少分别为C→T、G→A、T→C、A→G），颠换突变SNP的数量分别为133 102、117 465 和105 563 个（按数量由多到少分别为T→A、A→T、C→A、G→T、A→C、T→G、G→C、C→G）。并且在3 个样品中转换突变类型的SNP 数量约为颠换突变类型的2倍。

图1 转录组SNP突变类型分布图Fig.1 Mutation pedigree of transcriptome SNP

2.3 可变剪切分析

可变剪切事件一般分为五类：外显子替换、内含子保留、可变外显子、可变转录起始位点和可变转录终止位点。通过ASprofile 软件，将可变剪接事件细分为12 类，以ADL 为例发生概率从大到小分别为：第一个外显子可变剪切（alternative tran‐scription start site，TSS）、最后一个外显子可变剪切（alternative transcription termination site，TTS）、单外显子跳跃（模糊边界）（approximate exon skip‐ping，XSKIP）、单内含子保留（retention of single in‐tron，IR）、外显子替换（alternative exon ends，AE）、外显子替换（模糊边界）（approximate alternative ex‐on ends，XAE）、外显子跳跃（exon skipping，SKIP）、多外显子跳跃（cassette exons，MSKIP）、单内含子保留（模糊边界）（approximate retention of single in‐tron，XIR）、多外显子跳跃（模糊边界）（approximate cassette exons，XMSKIP）、多内含子保留（retention of multiple introns，MIR）和多内含子保留（模糊边界）（approximate retention of multiple introns，XMIR）（见表2）。可变剪接类型的数量模式在3个样品中基本一致，可变剪切事件数量为ADL>AF>AH，TSS和TTS类型的可变剪切数量明显高于其他类型的可变剪切。

表2 可变剪切事件统计表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 的上调表达基因数量分别为4 822、6 107和3 755，下调表达的基因数量分别为3 369、3 133 和5 079。除了AF vs AH，在其余两个对比组中上调表达基因均大于下调表达基因（见图2B）。经过筛选，一共鉴定出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 转录组测序差异表达基因A.样品间差异基因表达量聚类热图；B.样品间上下调表达基因数目分布情况；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 分类结果（见图3）表明，在ADL vs AH、ADL vs AF 和AF vs AH 这3组对比中，差异基因富集存在共同的分类特征，都分别注释到了GO 分类生物过程（biological process）的21 个条目、细胞组分（cellular component）的12个条目和分子功能（molecular function）的14 个条目。在生物过程分类中，3组对比中差异基因富集最多的两个节点分别为代谢过程（metabolic process）和细胞过程（cellular process），其中ADL vs AF 和ADL vs AH对比组在这两个代谢过程中主要富集上调基因。在细胞组分分类中，3个对比组中差异基因富集最多的两个节点分别为细胞（cell）和细胞部分（cell part）。在AF vs AH 对比组中，胞外部分（extracel‐lular region）和细胞外基质（extracellular matrix）都主要富集下调基因，这点与其他2组有所不同。在分子功能分类中，ADL vs AH 和AF vs AH 对比中差异基因富集最多的两个节点分别为酶调节活性（enzyme regulator activity）和金属伴侣活性（elec‐tron transfer activity）。ADL vs AF 对比组在金属伴侣活性节点中主要富集上调基因，而在AF vs AH对比组在该节点主要富集下调基因。在ADL vs AF 对比组中，差异基因富集最多的两个节点分别为催化活性（catalytic activity）和结合（binding）。在结合这一节点处，ADL vs AF 对比组主要富集上调基因，而ADL vs AH 和AF vs AH 对比组在该节点主要富集下调基因。

图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数据库。在ADL vs AF 对比组中，DEGs 富集到198 条通路，显著富集程度前三的代谢通路依次为MAPK 信号通路（MAPK signaling pathway）、类黄酮生物合成（Flavonoid biosynthesis）和硫代谢（Sulfur metabo‐lism），得到注释的DEGs最多的代谢通路是植物激素信号转导（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）和牛磺酸和亚牛磺酸代谢（Taurine and hypotaurine metab‐olism），得到注释的DEGs 最多的代谢通路是氨基酸的生物合成（Biosynthesis of amino acids），共注释到78 个DEGs，和植物激素信号转导（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 高表达差异基因

目前关于蒿属植物各基因的分子作用和机制未有系统性的研究，存在未知的功能基因，但对于未知的高表达差异基因也应该给予关注。以log2（FC）值作为筛选标准，选取各对比组内log2（FC）值最高的10 个基因进项功能分析（见表3）。ADL vs AF 对比组中存在4 个未知的功能基因，6 个已知的功能基因主要与应激和氧化还原过程相关。ADL vs AH 对比组中存在6 个未知的功能基因，4个已知的功能基因主要与氧化还原过程和应激反应相关。AF vs AH 对比组中存在3 个未知的功能基因，7个已知的功能基因主要与氧化还原过程和转运过程相关。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 共表达活性成分相关基因

根据KEGG pathway 富集分析结果，共筛选到25 个在3 个对比组中共表达的活性成分相关基因，根据其功能主要分为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种蒿属植物均为西北干旱、半干旱地区常见物种，在改善生态环境方面作用良多。并且作为传统中药材，具有良好的临床疗效和广泛的应用前景。这3 种蒿属植物在西北地区资源量较大，通过转录组比较分析，了解其序列信息和基因表达的差异情况，进一步明确差异基因的代谢途径和生物功能，以挖掘关键的功能基因，为植物资源利用提供分子基础。

通过对3 种蒿属植物成熟叶片进行转录组测序分析后，发现不同蒿属植物的转录本差异较大。通过差异表达基因分析，ADL vs AF 对比组拥有更多的差异基因，差异表达更明显，含有最多的上调差异表达基因，说明米蒿和冷蒿差异较大。米蒿和臭蒿同属艾蒿组，而冷蒿属于莳萝蒿组，同时差异基因表达量聚类热图也显示出米蒿和臭蒿关系较近，而与冷蒿关系较远。在3 种蒿属植物ADL、AF和AH转录组中分别找到了378 757、327 045和298 104 个可能的SNP 位点，转换突变类型的SNP数量约为颠换突变类型的2倍，这种情况在其他植物转录组的研究中均有发现。SNP 位点可用于物种鉴定，开发分子标记，本研究发现的SNP 可以为蒿属植物的物种鉴定、系统发育关系及物种亲缘关系分析提供分子依据。

可变剪切在生物中普遍存在，是调控基因互作的重要分子机制，在植物生长发育和非生物胁迫响应中起到重要的作用：调节开花时间、保证生物钟正常运转、响应温度胁迫和渗透胁迫。可变剪切能够增强植物对干旱和温度胁迫的适应能力，激发特定的分子功能。本研究中，测序数据的大小不同导致3 种蒿属植物可变剪切数量的不同，但是3个样品的可变剪切事件数量的分布呈现出一定的相似性，说明在3种蒿属植物中进行着较为一致的生理反应。本研究发现的可变剪接事件可为后续继续探寻蒿属植物相关生物反应的分子机制提供一定的生物信息学基础。

3 个对照组的高表达的差异基因与防御反应和氧化还原反应紧密相关，防御反应相关基因的高表达可能与野外采摘样品损伤叶片导致应激表达有关。氧化还原反应中高表达的C-4 甾醇甲基氧化酶活性在之前的研究中被证明与脱水胁迫密切相关，其相关基因可以作为培育耐旱植物的候选。同时，HSP20 家族蛋白在高表达已被证明参与高温、盐和干旱胁迫的反应，相关基因的高表达对应着3种蒿属植物生活的炎热干旱的环境，可为蒿属植物耐热机制研究提供相关的分子依据。

GO 富集分析结果显示3 个对照组在代谢过程、细胞过程和细胞这3个条目共同富集大量差异基因，而在胞外部分、结合、催化活性这3个通路明显不同。这说明3种蒿属植物具有一定的差异，冷蒿体内可能进行着更多的挥发性物质的相关反应。KEGG 富集分析结果显示，3 个对比组共同富集最多的通路为类黄酮代谢和甘油酯生物合成。类黄酮具有良好的抗氧化作用，可以提高植物对环境胁迫的抗性，研究证明干旱和炎热会导致蒿属植物中类黄酮的累积。甘油酯主要包括甘油三酯、甘油二酯、甘油一酯、糖脂、磷脂和甜菜碱脂，参与细胞膜脂的形成，是植物体内重要的初生代谢产物，对植物生长发育和抵抗环境胁迫具有重要作用。研究表明，在盐胁迫下，随着盐度的增加，双对栅藻（）中甘油三酯的含量呈现先升后降的趋势；在氮磷胁迫下，微藻中糖脂的累积量增加。本研究中的米蒿、冷蒿和臭蒿生长于青海省海西州的干旱盐碱地，这可能是导致类黄酮代谢和甘油酯生物合成这两个通路的共同富集的原因，推测相关基因与这蒿属植物对环境胁迫的适应有关。

通过GO和KEGG富集分析，共筛选出25个相关的共表达差异基因。相比于其他共表达差异基因，CTI12_AA555350 基因表达量非常高，该基因属于GERD 家族，其被证明是响应昆虫取食的关键基因，其产物大根香叶烯具有良好的驱虫作用。ALDH 能够使有毒的醛类氧化为无毒的羧酸，同时该基因能够响应植物非生物胁迫，在黄花蒿的转录组研究中，ALDH基因参与青蒿素生物合成途径。青蒿素是从黄花蒿中的一种内过氧化物倍半萜内酯（Endoperoxide sesquiterpene lac‐tone），在蒿属的多种植物中均有存在，相关基因及通路已有详细报道，然而其他蒿属植物的活性成分研究较少。本研究筛选的与活性成分相关的差异基因可为后续蒿属植物活性成分的研究提供参考信息。

本研究利用Illumina 测序技术对米蒿、冷蒿和臭蒿的成熟叶片进行了转录组测序，使用生物信息学的方法对测序结果进行分析。通过SNP、可变剪切、差异表达基因和功能富集分析，对3 种蒿属植物的基因结构和代谢通路有了一定的了解，并筛选了活性成分相关的差异基因。研究结果为今后研究蒿属植物生长发育、抗逆和活性物质和植物资源的开发与利用提供一定的科学依据。