余静雅，夏铭泽，张发起

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

甘青蒿（Artemisia tanguticaPamp.）属于菊科（Asteraceae）蒿属（ArtemisiaL.）蒿亚属（SubgenusArtiemisia）腺毛蒿组（SectionViscidipubes），广泛分布于青藏高原东南边缘中高海拔地区，为青藏高原特有种之一［1］，民间入药具有多种功效［2］。蒿属中多个物种被发现具有驱寒除湿、活血通经的功效，可用于治疗疟疾和感冒发烧等疾病［3］。

蒿属植物的挥发油具有抗菌、抗病毒和杀虫的作用［4-6］，其次生代谢产物与各种药理特性密切相关［7］。作为挥发油、树脂和蜡的重要成分，许多萜类次生代谢产物同时也在植物生理过程中起着重要作用［8］。青蒿素作为一种内过氧化物倍半萜内酯（Endoperoxide sesquiterpene lactone），在蒿属多个物种中均有存在［9］。甘青蒿整株密被腺毛或粘质柔毛，叶片中含有大量的挥发油。然而，目前对于蒿属植物的研究主要集中于生态功能、化学和药理学研究［10-11］，其生物分子研究远远落后［4］，仅有部分蒿属植物的叶绿体和转录组被测序［12-13］，以及黄花蒿（Artemisia annuaL.）的基因组被测序［14］。作为青藏高原广布的蒿属植物之一，由于地理分布和缺乏基因组序列，甘青蒿相关研究受到限制。对于缺乏生物分子信息的非模式物种，高通量RNA测序和生物信息学是挖掘特定生物功能基因的有力工具［15］。

二代测序技术具有通量高和准确性较高的优势，被广泛应用，但是由于其读长较短（＜ 300 bp），基于第二代高通量测序平台往往不能准确得到和组装出完整转录本。相对于一代和二代测序技术，第三代测序技术具有超长读长（平均读长可达15 kb）和直接读取目标序列的优势［16］，可以弥补二代测序技术的不足，提高基因注释的准确性［17］。目前，全长转录组测序技术已广泛应用于动植物研究中［18-20］，如利用PacBio测序技术揭示了艾蒿（Artemisia argyiL.）不同部位中萜类化合物生物合成相关基因的活性差异，并筛选出12种编码重要酶的候选基因［12］；利用牛津纳米孔技术（Oxford nanopore technology，ONT）获得茶树［Camellia sinensis（L.） O. Ktze.］不同部位的全长转录组，并检测到茶叶风味合成相关的途径［21］；对花生（Arachis glabrataBenth.）的全长转录组分析揭示了其对干旱、寒冷和盐等非生物胁迫的基因表达谱［19］；地黄［Rehmannia glutinosa（Gaert.） Libosch. ex Fisch. et Mey.］的全长转录组检测到了毛蕊花糖苷生物合成途径的关键基因［22］。目前尚无关于甘青蒿的类似研究，因此，本研究以甘青蒿叶片为材料，利用三代测序平台ONT对全长转录组进行测序及组装，鉴定和注释新转录本，寻找可能与蒿属萜类化合物生物合成相关的候选基因，并对甘青蒿全长转录组的转录因子（Transcription factor，TF）进行鉴定和分析，为甘青蒿生物资源利用提供一定的科学依据。

1 材料与方法

1.1 材料

用于测序的甘青蒿鲜嫩叶片采集于青海省海西蒙古族藏族自治州乌兰县（地理坐标：N 36°56′37"，E 98°27′20"；海拔2 939 m），植株叶片采集后迅速放入液氮中保存，用于后续RNA提取。凭证标本（chen2019104）存放于中国科学院西北高原生物研究所青藏高原生物标本馆（HNWP）。

1.2 方法

1.2.1 RNA提取与建库测序

采用Trizol法［23］提取甘青蒿叶片总RNA，琼脂糖凝胶电泳用于检查提取的RNA是否存在污染，并通过Nanodrop检测RNA浓度，使用NanoPhotometer spectrophotometer和安捷伦2 100生物分析仪（Agilent Technologies，Palo Alto，CA）评估RNA质量和完整性。检测合格的RNA样品使用Oligo DT作为引物对目标mRNA进行反转，通过低循环PCR扩增全长cDNA。NEBNext End Repair/dA-Tailing Module被用于末端修读及加A，ONT SQKLSK109试剂盒及NEBNext Quick Ligationg Module用于测序接头的连接。测序平台为PromethION（Oxford Nanopore Technologies公司，英国）。

1.2.2 全长转录组序列分析

对ONT原始测序数据进行质量控制，使用NanoFilt v2.8.0［24］过滤质量小于7和长度小于500 bp的序列（参数：-q 7 -l 50），过滤得到的有效数据使用SeqKit v0.12.0［25］默认参数进行统计。已发表的黄花蒿全基因组作为参考，使用minimap2 v2.17-r941［26］将过滤后的全长序列与参考基因比对，参数设置为：-ax splice -uf -k 14，samtools v1.11［27］用于结果统计，参数为flagstat。使用Pinfish（https：//tracker.debian.org/pkg/pinfish，版本：0.1.0；参数：default）软件对全长序列快速构建非冗余转录本集，得到一致序列。将得到的一致序列与参考基因组进行比对，使用StringTie v2.1.4［28］对比对结果去冗余，参数设置为：-conservative -L -R，合并仅5’端外显子有差异的比对，得到非冗余转录本序列。Gffcompare v0.12.1［29］用于新转录本及新基因的预测，补充现有注释。TransDecoder v5.5.0（https：//github.com/TransDecoder/TransDecoder）对新鉴定的转录本进行蛋白质编码（Coding Sequence，CDS）序列预测，参数为：-m 50，-single_best_only。

为获得转录本全面的功能信息，基于序列相似性和motif相似性，将新转录本注释到公共数据库，包括非冗余蛋白数据库（Non-Redundant Protein atabase，NR）［30］、蛋白质真核同源数据库（Eukaryotic Orthologous Groups，KOG）［28］、非冗余蛋白质序列数据库（Universal Protein，uniprot）［31］、蛋白质原核同源数据库（Cluster of Orthologous Groups of Proteins，COG）［32］、基因本体论数据库（Gene Ontology，GO）［33］、蛋白质家族域数据库（Protein Families Database，Pfam）［34］和东京基因与基金组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）［35］。

1.2.3 全长转录组转录因子分析

利用iTAK［36］结合数据库PlantTFDB v5.0［37］实现转录因子的预测。

2 结果与分析

2.1 全长转录组测序与组装

甘青蒿纳米孔测序产出的raw data达到15.66 Gb。NanoFilt对原始数据过滤后获得10.68 Gb的clean data，共得到14 730 754条序列，N50为999 bp，平均长度为779 bp。将过滤后的全长序列与参考基因组比对，共比对上9 557 431条序列，比对率达到81.6%。进一步过滤冗余全长转录组序列，得到18 634条非冗余全长转录本，N50为1 185 bp，最大长度为5 211 bp，平均长度为842 bp。

将非冗余转录本与基因组一致转录本进行比较，发现4 358条新转录本，主要为未知的新转录本（图1A）。对新鉴定的转录本进行CDS预测，获得1 482条CDS序列，N50为816 bp，最大长度为3 492 bp，平均长度为695 bp，长度主要分布在300 ～ 1 000 bp（图1B）。

图1 新鉴定转录本Fig.1 Newly identified transcript

2.2 新鉴定转录本功能注释

将新鉴定的转录本分别注释到7个数据库中，4 358条转录本中共有2 184条（50.11%）被成功注释（表1）。其中共有2 172条转录本比对到NR数据库的122个物种，比对率较高的前3个物种都属于菊科，分别是黄花蒿（979条，48.51%）、除虫菊［Tanacetum cinerariifolium（Treviranus） Schultz Bipontinus］（649条，32.17%）和向日葵（Helianthus annuusL.）（116条，5.71%）。

表1 新鉴定转录本注释结果Tab. 1 The result of newly identified transcript annotation

选取NR、GO、KEGG、KOG和uniprot这5个常用的数据库注释情况绘制韦恩图（图2），915条转录本在5个数据库中均有注释。

图2 新鉴定转录本功能注释韦恩图Fig.2 Venn diagram of function annotated newly identified transcript

2.2.1 GO注释

GO数据库中生物学过程、细胞组分和分子功能3个分支注释到的新鉴定转录本的数量分别为1 064、1 968和1 322条，共注释到955个功能组，选取每个分支下前20个注释最多的功能组绘图（图3）。在生物学过程中注释数量最多的是翻译（82条）、碳水化合物代谢过程（31条）和甲基化（21条）。在细胞组分中注释数量最多的是膜的组成（423条）、细胞核（134条）和细胞质（100条）。在分子功能中注释数量最多的是ATP结合（121条）、金属离子结合（114条）和核糖体结构组成部分（92条）。

图3 新鉴定转录本GO注释结果Fig. 3 GO annotation classification of newly identified transcript

2.2.2 KEGG注释

KEGG数据库注释结果显示，新鉴定转录本共注释到5条主通路（细胞过程、环境信息过程、遗传信息过程、代谢和组织系统）和18条子通路（表2），其中注释数量最多的子通路为全局和概览图（673条）、碳水化合物代谢（161条）和翻译（140条）。共有70条转录本被注释到可能与萜类化合物代谢相关子通路中，其中29条转录本被注释到萜类和多酮类化合物代谢子通路，41条转录本被注释到辅助因子和维生素代谢子通路。

表2 KEGG注释结果Tab. 2 KEGG annotation classification of A.tangutica

对注释到可能与萜类化合物代谢相关子通路的转录本进行进一步筛选，去除不完整的转录本及注释到不相关子通路上的基因（诸如卟啉和叶绿素代谢、烟酸和烟酰胺代谢、核黄素代谢等）。新鉴定转录本中共有9条完整转录本涉及萜类生物合成途径（表3），可能是与甘青蒿叶片萜类化合物生物合成相关的候选基因，经过序列比对之后，其中涉及类单萜生物合成的转录本可能编码ctcad3和CTI12_AA295060基因。涉及萜类化合物骨架生物合成的转录本可能编码CTI12_AA344030、MCS和hdr基因。涉及泛素酮和其他萜类醌的生物合成的转录本可能编码CTI12_AA436250、CTI12_AA496620和CTI12_AA199400基因。

表3 萜类化合物生物合成相关基因Tab. 3 Genes related to terpenoid biosynthesis

2.3 转录因子分析

iTAK用于甘青蒿全长转录组序列的转录因子预测，结果获得29 154个转录因子，分别属于58个不同的转录因子家族，选择数量最多的前20个转录因子家族进行可视化（图4），前3位分别是bHLH家族（2 936个，10.07%）、NAC（2 325个，7.97%）和MYB_related（1 841个，6.31%）。这些转录因子信息为甘青蒿次生代谢产物的生物合成和抗逆性研究提供了依据。

图4 甘青蒿全长转录组转录因子家族分析Fig. 4 Analysis of transcription factor family in full-length transcriptome of A. tangutica

3 讨论

本研究利用三代测序平台ONT结合生物信息学分析，获得了甘青蒿全长转录组数据。通过聚类、去冗余和校正后最终得到18 634条非冗余全长转录本，N50为1 185 bp。N50大于1 000 bp且GC含量稳定，说明序列组装完整性较好，能够满足后续生物信息学分析的要求。

本研究中，共有4 358条新转录本被鉴定，其中2 184条（50.11%）得到7个数据库的注释。相比于艾蒿全长转录组高达93.97%的注释率［12］，甘青蒿注释率低原因可能为该物种与黄花蒿及艾蒿亲缘关系较远。此外，黄花蒿及艾蒿多分布于低海拔区域，而甘青蒿主要分布于青藏高原，由于生境不同，甘青蒿可能拥有独特的生理过程。这一情况也反映出公共数据库中菊科植物遗传信息的匮乏。因此，本研究获得的甘青蒿全长转录组能够丰富菊科植物遗传数据库；为探究青藏高原植物独特适应性提供基础数据；为挖掘相关功能基因提供参考依据。

萜类化合物能够抑制多种细菌和蚜虫，具有显著的药理活性，如免疫调节、抗肿瘤、降血脂等［38］，广泛应用于制药及食品工业。蒿属植物富含萜类化合物，对其萜类化合物进行制取，可以提高蒿属植物的经济价值［12］。通过将转录本注释到KEGG数据库，共筛选出9个可能参与萜类生物合成的完整候选基因，其中6个为新鉴定出的候选基因。2-C-甲基-D-赤藓糖醇4-磷酸（2-C-methyl-D-erythritol 4-phosphate，MEP）途径为生物合成萜类化合物的主要途径之一。MCS基因和hdr基因参与MEP途径，同时也是MEP途径中的两个关键酶［39］。其中，MCS基因则参与MEP途径的第五步，该基因编码2-C-甲基-D-赤藓糖醇2，4-环二磷酸酯合成酶（2-C-methyl-D-erythritol 2，4-cyclodiphosphate synthase），负责催化2-C-甲基-D-赤藓糖醇2，4-环二磷酸酯（2-C-methyl-D-erythritol 2，4-cyclodiphosphate）的合成，是二萜类化合物生物合成的关键基因之一［40］。hdr基因参与MEP途径的最后一步，该基因编码4-羟基-3-甲基丁-2-烯基二磷酸还原酶（4-hydroxy-3-methylbut-2-enyl diphosphate reductase），负责将4-羟基-3-甲基丁-2-烯基焦磷酸酯（4-Hydroxy-3-methylbut-2-enyl pyrophosphate）催化成焦磷酸异戊酯（Isopentenyl pyrophosphate）和二甲基烯丙基焦磷酸酯（Dimethyl allyl pyrophosphate）。研究表明该基因的过量表达不仅可以提高青蒿素的产量，也增加了黄花蒿中倍半萜类化合物和单萜类化合物的含量［41］。同时，该基因在艾蒿全长转录组分析中也被鉴定为萜类化合物生物合成的关键基因［12］。CTI12_AA199400基因编码细胞色素P450s，主要参与萜类化合物生物合成过程中多种修饰酶的合成，在植物萜类生物合成中起着重要的作用［42］。然而，其他候选基因对萜类化合物生物合成影响的研究则少见报道。

在此前的研究中，黄花蒿基因组鉴定出2 717个转录因子［14］，艾蒿全长转录组鉴定出5 604个转录因子［12］，本研究共鉴定出29 154个转录因子，丰富了蒿属植物的转录因子信息。已发现AP2、bHLH、MYB、NAC、WRKY和bZIP转录因子家族通过激活或抑制萜类生物合成途径中关键基因的表达来调节萜类代谢［43-44］。丹参（Salvia miltiorrhizaBunge）MYB转录因子能够促进相关基因的表达，增加丹参酮（Tanshinone）的产量［45］。黄花蒿bHLH转录因子通过激活关键酶的启动子来调节青蒿素的生物合成［46］，并参与植物的萜类化合物的生物合成［47］；WRKY［48］、HD-ZIP［49］和MYB［50］与青蒿素的合成密切相关。本研究中，同样鉴定到了这些转录因子家族的表达，由此推测甘青蒿体内可能存在青蒿素的生物合成，具体过程还需要进一步研究。

青藏高原不同寻常的极端环境通常会诱导植物的各种抗逆反应［51］。本研究涉及萜类生物合成的候选基因中，ctcad3编码乙醇脱氢酶，其过量表达能够提高植物对盐、干旱、寒冷和病原体感染的抗逆性［52］。CTI12_AA295060基因编码的葡萄糖/核糖醇脱氢酶作为一种醇-多糖脱氢酶（Alcohol-polyolsugar dehydrogenase），可能具有碳水化合物代谢和抗旱能力［53］，在干旱和盐胁迫下会大量表达，以增强植物的耐盐性［54］。此外，有研究发现，在非生物胁迫下，植物中NAC、bZIP、MYB和WRKY家族的转录因子会大量表达［55-56］。本研究相应基因及转录因子的表达，可能反映了甘青蒿对青藏高原极端环境的适应性。

本研究利用ONT测序技术对甘青蒿全长转录组进行测序和组装，鉴定新基因及新转录本，并对其进行功能注释，筛选出的9个与萜类化合物生物合成相关候选基因，有助于为甘青蒿中新基因的发现和功能基因分析提供可靠数据。分析甘青蒿全长转录组的转录因子，能够为甘青蒿功能基因挖掘和生物资源利用提供参考。