高寒灌丛建群种高山绣线菊的低温适应转录组信息分析

韩 霜，余静雅，韩 赟，张发起

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

高山绣线菊（Spiraea alpinaPall.）是青藏高原高寒灌丛主要建群种之一，隶属于蔷薇科（Rosaceae）绣线菊属（Spiraea），主要生长在海拔2 000~4 000 米的向阳坡地或灌丛中［11］。该植物广泛分布于整个横断山区，生长于年平均气温较低的高寒甸，具耐寒、耐旱、耐瘠薄、耐阴湿等特点［12-14］。高山绣线菊可用于治疗咽喉肿痛、风热痒症，是一种广泛应用于民间的中草药，其根、叶、果实可做兽药［15］。目前，对高山绣线菊的研究主要集中在化学成分、生物活性、繁殖栽培及冰期演化历史上［16-20］，学者们证实从该植物中分离纯化的抗真菌化合物对植物病原真菌具有一定的抑制作用［17］；作为高海拔地区园林绿化的少数树种之一［16］，不少学者对该植物的驯化和栽培技术进行探索，指出其育苗技术简单，栽培后成活率高且适应性强［21-22］。此外，一些学者依据高山绣性菊叶片膜透性及膜脂过氧化特征，验证其抗寒性强度［23］。然而关于高山绣线菊分子生物学的研究较少，尤其是基因组及转录组方面的研究。

随着高通量测序技术的发展，使得转录组测序成为挖掘功能基因、筛选分子标记、阐明代谢途径的有效工具［24］。目前，已成功应用到该技术的模式生物有水稻（Oryza sativa）［25］、玉米（Zea mays）［26］、拟南芥（Arabidopsis thaliana）［27］等。在现有的植物转录组分析研究报道中，对灌木沙棘（Hippophae rhamnoides）进行转录组测序并与拟南芥比较后，最终确定与沙棘脂肪酸生物合成相关的基因序列［28］；三江源地区灌木亚菊的转录组解析成功获取了药用活性成分相关的代谢通路，为其资源利用和保护、遗传多样性提供基础数据［29］。此外，利用转录组测序进行植物低温胁迫相关的研究不少，如模式生物水稻转录组信息分析明确了其苗期植物激素对低温的应答机制［30］；骤然低温下草本植物川百合（Lilium davidii）的转录组分析成功筛选出与抗寒性、光合作用、代谢途径等相关关键基因，为川百合的分子育种提供参考［31］；木本植物仁用杏花（Prunus armeniaca）的转录组分析发现了抗寒差异表达基因，为今后采取基因工程手段培育抗寒的仁用杏新品种提供基因资源［32］；这些成功的案例为今后高山绣线菊适应胁迫的分子机制研究提供参考。

高山绣线菊资源丰富，适应性强［14］。目前，对该植物的研究相对较少，尤其是非生物逆境胁迫方面的认识较浅薄。因此，我们对高山绣线菊进行高通量测序并获得转录组数据，对基因进行功能注释并分析低温胁迫相关代谢通路，挖掘关键基因，对后续其低温胁迫分子机制的研究具有一定指导意义。

1 材料与方法

1.1 材料

高山绣线菊新鲜幼叶于秋季（9 月份）采集于青海省黄南州同仁县采集（地理坐标：35°13′58.30″N，101°51′05.50″E；海拔：3 532米；年平均气温：5.2℃），用超纯水和75%酒精清洗后迅速置于液氮中，后转移至-80℃的超低温冰箱中保存备用。凭证标本（标本号：Zhang2018047）存于中国科学院西北高原生物研究所青藏高原生物标本馆（HNWP）。

1.2 方法

1.2.1 高山绣线菊RNA提取与建库测序

利用Total RNA Extractor（Trizol）试剂法［33］提取高山绣线菊总RNA，并用琼脂糖凝胶电泳和NanoDrop-TM2000c（Thermo Scientific，美 国）检测RNA 质量纯度。高山绣线菊转录组测序文库的构建参照夏铭泽在多裂骆驼蓬中的建库方法［34］。检测合格的文库用Illumina HiseqTM进行测序。

1.2.2 转录组拼接与Unigene功能注释

利用FastQC v0.11.2（http：//www.bioinformatics.babraham.ac.uk/projects/fastqc/）对测序数据进行质量评估和质控，过滤原始数据（Raw reads）中含有带接头的、低质量的序列。通过Trimmomatic v0.36［35］进行质量剪切，得到Clean reads。使用Trinityv2.4.0［36］（参数min_kmer cov=2，其余为默认设置）对样本有效数据进行混合拼接，Trinity 组装后的转录本序列信息以FASTA 格式储存，并对转录本Transcript 去冗余，取每条基因中最长的转录本作为Unigene，以此作为后续分析的参考序列。

不进则退，李高明在蒙自花了三万块，挨着一家生意红火的理发店，街头街尾地打起了“价格战”。此时的李高明，是小店五六个师傅中技术最好的，凡事他都想亲力亲为。有时生意好，给客人做头做到一两点，客人看晚了回不去，竟也乐意在他的小店里睡一宿。生意不好时，他犯愁，整夜整夜睡不着。心情起起伏伏地过了三个月，实在太煎熬，他受不了了。

随机从Clean 数据中抽取10 000 条序列，与多个数据库进行比对，取evalue<=1e-10 并且相似度＞90%，coverage＞80%的比对结果作为后续分析基因功能注释及分类的数据，所采用的数据库有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 数据库的比对，可得到高山绣线菊转录本序列与相近物种的近似情况以及同源序列的功能信息。根据Unigenes 与Swissprot、TrEMBL 的注释结果得到GO（Gene Ontology）功能注释信息，统计分子功能、细胞成分、生物过程三大分类中注释成 功 的Unigenes。利 用KAAS v2.1［37］得 到 转 录 本KEGG注释信息。

1.2.3 SNP与SSR分析

利用BCFtools v1.5［38］（参数：质量值大于20且覆盖度大于8）将组装好的Unigene 作为参考序列进行单核苷酸变异（Single nucleotide polymorphisms，SNP）分析，找出可能的单核苷酸变异位点并统计筛选出SNP 突变类型。利用MISA v1.0［39］进行微卫星（Microsatellite）标记或简单序列重复（Simple sequence repeat）标记分析，分别设置二核苷酸、三核苷酸、四核苷酸、五核苷酸以及六核苷酸重复单元的重复次数为至少6、5、5、5、5次。

2 结果与分析

2.1 转录本组装

高山绣线菊转录组高通量测序共获得51 340 802条Raw reads，经过滤处理后，得到49 947 114 条Clean Reads，总长为7 212 398 740 bp，GC 含量为49.37%，根据Q20（碱基质量在20 以上的序列，占97.71%）和Q30（碱基质量在30 以上的序列，占93.68%）信息统计结果，可以说明测序所得序列的质量满足后续的转录组分析。经转录本组装，共获得117 280 个Trancripts，53 892 个Unigenes（表1）。其中200~300 bp 的Transcript 和Unigene 数量最多，1 900~2 000 bp的最少（图1）。

表1 高山绣线菊转录组拼接结果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 个数据库进行比对（表2）。注释到NR 数据库的Unigene 最多，占39.49%，注释到KEGG 数据库的Unigene 最少，占6.45%。35 954 条Unigenes 注释到至少一个数据库中，占66.71%，2 310 条Unigenes 能注释到所有数据库中，占4.29%。此外，还有17 938 条Unigenes 并未注释成功。注释到NR 库的Unigenes 有31 627 条，其中比对到蔷薇科梅（Armeniaca mumeSieb.）物种的序列数最多（8 163 条），其次分别为桃（AmygdalusL.，6 720 条）、苹果属（MalusMill.，4 000 条）、梨属（PyrusL.，3 544 条）、草莓属（FragariaL.，858 条）等植物。此外，还有8 325 条序列（26.32%）比对到其他525 种植物，但每个物种所能比对上序列都较少，其原因可能与高山绣线菊近缘种中基因组数据库匮乏有关。

GO 可以全面描述生物体中基因和基因产物的属性［41］。Unigene 与Swissprot、TrEMBL 比对后，有29 111 条Unigenes 获得223 924 条注释信息。根据注释结果对得到的基因进行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）、形态发生活性（Morphogen activity）、生物（Biological phase）等子类中注释到的Unigene 相对最少，分别有4 条（0.014%）、1条（72.7%）、1条（0.0034%）。

Unigene 与KOG 数据库比对后，16 166 条基因被注释成功，按KOG 的group 可分为26 个类型（表3）。其中，翻译后修饰、蛋白转运、信号传递机制和只有一般功能预测分类下的基因较多，分别有1 962条（12.1%）、1 953 条（12.1%）和1 850 条（11.4%），而细胞外结构、核结构和细胞活性注释的基因数量较少，分别有60 条（0.33%）、48 条（0.27%）和10 条（0.056%）。此外，还有779 条注释成功的基因未知其功能。

根据KO 与Pathway 的关联性对其进行KEGG代谢通路分类。根据与KEGG 数据库的比对结果，成功注释到3 475 条Unigene，占6.45%。3 475 条基因被分为了代谢（Metabolism）、遗传信息处理（Genetic information processing）、细胞过程（Cellular processes）和环境信息处理（Environmental information processing）四大类23 个子类。在四个大类中，涉及基因最多的是代谢，共有2 367 条，占68.1%，其次为遗传信息处理、细胞过程和环境信息处理，分别有1 450 条（41.7%）、606 条（17.4%）和539 条（15.5%）。在23 个子类中，与代谢相关的通路最多，有12 条，与细胞过程、环境信息处理和遗传信息处理三大类涉及的通路较少，分别有4 条、3 条和4 条。代谢途径中，注释基因最多的通路是翻译（Translation）、信号转导（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 数据库分析，共统计到213 代谢途径，其涉及6 560 条Unigenes，按注释基因数量从高到低排序，选择前11个代谢通路进行分析并列于表4 中。注释基因最多的代谢通路为核糖体（Ribosome），有379 条，占总数的5.78%，其次为碳代谢（Carbon metabolism）和氨基酸合成（Biosynthesis of amino acid），分 别 为164 条（2.5%）和155 条（2.36%）。

表4 高山绣线菊Unigene数量最多的11个代谢通路Tab.4 The top eleven metabolic pathways involved in the largest number of S. alpina Unigenes

2.3 高山绣线菊低温胁迫代谢通路分析

在213 条代谢途径中，对高山绣线菊低温胁迫代谢通路进行统计及分析，可分为低温胁迫生理代谢（Physiological metabolism of cold resistance）、冷调节信号通路（Cold regulation signal pathway）、光合作用（Photosynthesis）等主要途径（表5），分别有14 条（碳代谢、氨基酸类的生物合成、淀粉和蔗糖代谢为主）、4 条（光合作用、光合作用生物中的碳固定为主）和3 条（植物激素信号转导为主）。所有代谢通路中，碳代谢（Carbon metabolism）、氨基酸类的生物合成（Biosynthesis of amino acids）、植物激素信号转导（Plant hormone signal transduction）、淀粉和蔗糖代谢（Starch and sucrose metabolism）代谢通路涉及的Unigenes最多，分别为164条、155条、137条、101条；而亚油酸代谢（Linoleic acid metabolism）、甜菜碱生物合成（Betalain biosynthesis）代谢通路涉及的最少，分别为19条和13条。

表5 高山绣线菊低温胁迫代谢通路及基因统计Tab.5 Metabolic pathway and gene statistical table of active components for S.alpina

3 讨论

本研究利用Illumina Hiseq 对高山绣线菊转录组进行高通量测序，共获得49 947 114 条Clean Reads。经组装共获得117 280个Transcripts和53 892个Unigenes，其平均长度为708.72 bp，N50 为1 340 bp。相比于其他植物转录组测序及组装结果［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），根据有效数据Q20 值、Q30 值，表明高山绣线菊转录组测序及组装质量较好，能够满足后续转录组信息分析的要求。

与9 个数据库比对后，成功注释的Unigenes 数量为53 892 条，而所得结果中仍然还有17 938 条Unigenes 未与已知基因匹配成功。这在其它植物的分析中也有发现，如珠子参（Panax japonicusvar.Major）［45］、虎杖根（Polygonum cuspidatum）［46］、金钱松（Pseudolarix amabilis）［47］等，这可能与基因库中缺少该种相关的基因组信息有关；其次，缺乏转录组方面的研究作为参考，导致该种特有的某些基因和数据库中序列的识别和比对十分困难。KOG 结果几乎整合了高山绣线菊所有信息，为其基因功能研究提供了数据基础，成功注释到16 166条Unigenes，其中参与该物种翻译后修饰、蛋白转运和信号传递机制等生命活动通路的基因最多，可以推断这些基因在该物种中表达较丰富，说明这三类生命活动在该物种的生长发育中具有重要的地位。NR 数据库中比对到梅的序列数最多，表明二者具有较高的序列同源性，亲缘关系较近，但仅有5个同科植物与该物种比对成功，其原因可能与蔷薇科植物的转录组及基因组信息较匮乏或该种本身具有的特异基因较多有关。

在低温胁迫下，植物通过改变生理生化、响应抗寒分子机制等方式提高自身抗寒性。具体而言，植物细胞膜系统、光合作用、可溶性糖/蛋白含量、游离脯氨酸含量等均受到影响而发生变化，植物也会对低温作出响应并进行一系列的冷信号转导途径，如碳水化合物（蔗糖和淀粉）代谢，植物激素合成及转导，次生代谢产物合成，信号转导（Ca2+），光合系统，脂质代谢等［10，48］。本研究依据KEGG 数据库对低温胁迫相关代谢通路的筛选结果与青海草地早熟禾一致［49］。将Unigene 映射到与低温胁迫相关的代谢通路中，对其进行初步筛选，结果显示低温胁迫生理代谢涉及的通路占大多数，碳代谢、氨基酸类的生物合成等通路所涉及的Unigene 最多，进一步说明这些代谢通路与高山绣线菊的耐寒性特征密切相关。

此外，低温胁迫下的植物是基于基因表达的迅速变化而进行的转录调控，从而合成相应蛋白指导代谢物的生成，以此对胁迫做出响应，如同科植物苹果的转录组分析发现植物激素信号转导、光合作用、糖酵解、淀粉和蔗糖等代谢通路在低温处理后发生基因的差异表达，前者通路中差异基因上调的数目高于下调的数目，而后两者通路中则呈现相反的情况［48，50］。与苹果转录组相比，高山绣线菊映射到这些代谢通路的Unigene 数量较多，这可能与该植物在低温环境下对这些代谢通路的依赖性较高有关，进一步说明了两种植物受到的胁迫程度和响应胁迫的基因数量具有差异，而这些代谢通路在高山绣线菊响应低温胁迫时发挥着重要作用，这为低温胁迫相关关键差异表达基因的寻找提供了线索。

4 结论

本研究利用分子手段获得了高山绣线菊的转录组数据，在无参考基因组的前提下，对这些数据进行分析，包括基因比对、功能注释、代谢通路，为其后续的分子生物学研究提供了基础数据，进一步丰富了蔷薇科植物转录组数据库，同时筛选出与低温胁迫相关的代谢通路，挖掘关键基因，对后续其适应低温胁迫分子机制的研究具有一定指导意义。