穿心莲全长转录组测序及特性分析

孙铭阳，徐世强，顾 艳，梅 瑜，周 芳，李静宇，王继华

（广东省农业科学院作物研究所/广东省农作物遗传改良重点实验室，广州510640）

0 引言

植物次生代谢产物是生物活性物质的主要来源，它们具有治疗多种疾病的临床应用价值[1]。穿心莲(Andrographis Paniculata)是一类重要的“寒性”中草药，其地上部分入药后可降温消炎的记载历史悠久[2]。穿心莲次生代谢提取物中含有20多种结构类似的二萜类化合物和10多种黄酮类化合物。其中，含量最多的二萜类物质为穿心莲内酯，药理活性丰富[3]。目前，新型冠状病毒2019(2019 novel coronavirus，2019-nCoV)肺炎疫情已成为全球性的卫生危机。因穿心莲内酯及其荧光衍生物可通过与主要蛋白酶功能位点形成共价键的方式抑制2019-nCoV的活性，现已被建议用于新冠肺炎的治疗[4]。

穿心莲内酯合成相关的分子机理研究已取得部分进展。穿心莲内酯可通过下调磷脂酰肌醇3激酶/蛋白激酶b信号通路基因来抑制肿瘤坏死因子-α诱发的炎症[5]。13C核磁共振常用同位素富集光谱显示，脱氧木酮糖途径是二萜类前体被修饰为穿心莲内酯所经历的主要途径[6]。随高通量测序技术的兴起，二代Illumina HiSeq短序列测序(Second Generation Sequencing,SGS)已成为分析基因表达水平和调控途径的常用方法。Tong等[7]利用SGS数据揭示了激素对穿心莲种子萌发过程的代谢及信号转导进程的影响。Cherukupalli等[8]将SGS输出数据进行二次组装后，生成了大小为88 Mb的穿心莲叶片转录组数据库，共83800条序列片段，为其次生代谢途径研究提供了初始平台。然而，SGS技术的读长有限；组装数据时需拼接；无法精准评估可变剪切产生亚型的转录丰度及较长转录本的全长信息。

如今，以美国太平洋生物技术公司的PacBio平台为主的第三代超长单分子实时(Single Molecule Real Time,SMRT)测序技术已经成熟。平均15 kb的超长读长使数据无需拼接，可直接得到高质量的转录本全长序列和AS亚型信息[9]。为进一步明晰穿心莲次生代谢途径基因的响应情况，本研究利用PacBio平台的三代SMRT测序技术对生长60天的穿心莲根、茎和倒三叶进行全长转录组测序和分析。结果表明，穿心莲体内功能基因主要参与次生代谢进程；预测功能为转录因子(Transcription Factors,TF)的基因中，basic helixloop- helix(bHLH)、 myeloblastosis(MYB) 及WRKYGQK(WRKY)等响应植物次生代谢的主要TF家族成员含量排名前10；合成二萜类前体的甲羟戊酸途径(MVA)和2-C-甲基-D-赤藓糖醇-4-磷酸途径(MEP)基因发生的主要AS事件为内含子保留；在搜索到的简单重复序列(Simple Sequence Repeat,SSR)中，AT/AT为优势重复单元。以上结果为穿心莲次生代谢基因调控网络研究、AS亚型的功能注释以及穿心莲遗传多样性的评价提供分子基础。

1 材料与方法

1.1 实验材料

选用福建漳州生产用穿心莲种，种于草炭土(HAWITA)中，保持土壤湿润，于光照培养箱（光照day/night：16 h/8 h；温度：28℃）中培养60天。将幼苗从土中移出，清洗根部后轻柔擦干，剪取根、茎和倒三叶并迅速混合置于液氮速冻，每个部位取3次独立的生物学重复并混合，待RNA提取。穿心莲生长于广东省农业科学院作物研究所农作物遗传改良重点实验室，于2020年4—6月进行。

1.2 全长转录组文库构建

利用塔克拉生物科技公司的RNAiso Plus试剂提取样品的总RNA；使用Nanodrop检测样品浓度；利用Agilent 2100检测样品片段完整度。质量检测合格的RNA样品用于穿心莲全长转录组文库构建，建库步骤参考Yang等[10]。

1.3 全长转录组测序数据筛选

建库完成后，上机（Pacbio三代测序仪）测序。利用SMRT Link v6.0软件对输出的原始序列进行分析。将原始数据中的环形一致性序列(Circular consensus sequence，CCS)提取后，通过CCS对序列中是否存在测序引物序列或者是否存在嵌合序列进行筛选、去冗余和校正[11]。最终得到高质量全长转录本序列，用于后续功能注释等分析。

1.4 全长转录组功能注释

将转录本序列与National Center for Biotechnology Information(NCBI)网站中的blastx作比对，得到相似性最高的蛋白序列，用于进行蛋白功能注释。利用Eukaryotic ortholgous groups(KOG)进行蛋白的直系同源分类；利用swiss prot protein database(SwissProt)数据库进行去冗余后的蛋白功能注释；利用kyoto encyclopedia of genes and genomes(KEGG)数据库进行基因功能分类和代谢通路聚类；与nonredundant protein database(Nr)数据库进行比对，取比对结果中E值最低（同源性最高）的序列所属物种进行各物种的同源序列数量统计[12]。

1.5 基因结构分析

转录本在以上4个数据库的蛋白比对结果中选取E值小于e-5，且rank最高的多肽序列对应的碱基区域定为编码区(Coding Sequence,CDS)。利用ANGEL软件预测与以上4个数据库比对不上的转录本序列的CDS区和蛋白序列[13]。对没有注释到以上4个数据库的转录本序列进行lncRNA分析，利用cnci软件和CPC软件进行编码能力的预测，取2个软件均预测为“非编码”的序列作为最终的lncRNA结果[14-15]。利用Signal P网站进行蛋白信号肽预测。利用MISA软件搜索SSR，配置参数信息为：1、2 nts至少重复6次；3 nts至少重复5次；4 nts至少重复4次；5 nts至少重复4次和6 nts至少重复4次才会被认为是SSR。2、如果2个SSR的距离短于100 bp，就会被合并作为一个SSR。对于AS事件，本研究利用Cogent软件组装出经聚类后转录本的CDS序列，并以组装CDS序列作为参考，利用SUPPA软件进行可变剪切分析[16]。

2 结果与分析

2.1 穿心莲全长转录组测序的数据筛选

浓度合格的样品总RNA经双链cDNA合成及加接头，形成哑铃型单分子双链cDNA文库后上机测序。通过提取下机序列中的CCS，共得到全长序列487309条，其中全长非嵌合(Full-Length Non-Chimeric,FLNC)序列354407条。FLNC序列在去掉5’引物、3’引物和polyA序列后即为物种原始全长转录本序列。由图1可知，穿心莲幼苗的FLNC序列长度主要集中在200～ 3200 nts之间。因基因组时刻处于转录mRNA的动态中，所以测序得到的单分子FLNC序列存在大量冗余。将FLNC序列进行聚类、矫正和去冗余后，最终得到76812条高质量转录本信息，用于后续功能注释和结构分析。

图1 FLNC序列长度分布图

2.2 穿心莲全长转录组的功能注释

利用四大功能注释数据库KOG、Nr、Swiss Prot及KEGG对76812条穿心莲全转录组信息进行注释，各数据库比对情况见图2。共有63855条序列得到了注释，其中Nr数据库注释到的信息最多，共63822条序列；Swissport数据库注释到54430条序列；KOG数据库注释到43656条序列；KEGG数据库注释到30412条序列。由4个数据库共同注释到的序列有25010条。

图2 四大数据库注释穿心莲全转录组序列韦恩图

Nr注释结果显示，与穿心莲最相近的物种中的前三位是唇形目的芝麻(Sesamum indicum)（比对到34684条序列）、管状花目的猴面花(Erythranthe guttata（)比对到7101条序列）以及管状花目的旋蒴苣苔(Dorcoceras hygrometricum（)比对到1790条序列）。

KOG注释结果发现，穿心莲蛋白除参与维持正常生命活动所必须的途径（RNA进程、转录进程和常规功能）外，主要富集的功能途径为：翻译后修饰及蛋白质转换、信号转导及胞内运输、分泌和囊泡运输。众所周知，植物次生代谢产物是由公共前体经特异性蛋白修饰而成，如抗生素、激素和色素等。这些产物形成后可贮存在液泡和细胞壁中，部分也可运输到体外环境。在贮存和运输的过程中，需要细胞内的信号转导和囊泡运输等功能的帮助。因此，该结果体现了穿心莲细胞生命活动主要为次生代谢的特点。KEGG聚类进一步证明KOG数据库的注释结果。如表1。穿心莲功能基因聚类最多的途径为代谢途径、次生代谢产物合成途径和抗生素合成途径。

表1 穿心莲全转录组信息的KEGG注释

2.3 穿心莲全长转录组的基因结构分析

经过CDS和蛋白预测，穿心莲全长转录组信息得到进一步注释。本研究共识别出64401条mRNA序列和12411条长链非编码RNA(lncRNA)序列。其中，4008条序列的预测翻译产物含信号肽结构。

AS事件可在转录后水平丰富基因的转录本种类。利用生物信息学软件将全长转录组数据中的序列进行聚类，可比较同一聚类下的任意两条转录本发生AS事件的位置。经统计，共有11273个穿心莲基因发生了AS事件，体现出其转录本的多样性。其中，产生2个AS亚型的基因最多，超过总数的1/4（2845个，25.24%）。数量次之的为产生3个和4个AS亚型的基因，分别有1935个(17.16%)和1389个(12.32%)。本研究中出现次数最多的AS事件为内含子保留事件，其中启动子可变型内含子保留事件发生频率最高。其次为5’端可变事件。

穿心莲内酯的二萜类前体骨架(E,E,E)-香叶基香叶基二磷酸酯[(E,E,E)-geranylgeranyl pyrophosphate，GGPP]的合成原料为质体内的MEP途径和胞质中的MVA途径产生的异戊烯基二磷酸(IPP)和二甲基烯丙基二磷酸(DMAPP)[17]。本研究对MEP和MVA途径基因进行AS事件分析。由图3可知，两途径的基因均发生了AS事件。GGPP合成前期的基因主要以内含子保留和5’端可变为主。生成IPP和DMAPP后，香叶基香叶基焦磷酸合成酶(GGPPS)基因和法呢基二磷酸酯合成酶(FPPS)基因发生了内含子保留、3’端可变和外显子跳跃AS事件。由表2可知，在发生内含子保留的 GGPPS 基因中，Isoform0050046、Isoform0042454、Isoform0041790、Isoform0053892、Isoform0054050 和Isoform0054809为启动子可变式内含子保留AS事件产生的转录本亚型。

表2 GGPP合成相关基因AS事件统计表

图3 穿心莲MEP和MVA途径的AS事件示意图

有2803条序列的预测蛋白为TF家族成员，含量最多的前10类如图4所示。有研究表明，bHLH、basic leucine zipper(bZIP)、MYB及WRKY家族成员可直接调节次生代谢途径基因的表达水平[18]。

图4 含量排名前十位的TF家族

SSR位点广泛存在于基因组中，具有重复性好、共显性高、稳定性强且多态性丰富等特点，可作为分子标记来辅助育种[19]。本研究在穿心莲全长转录组中共找到45550处SSR位点，分布在27236条序列上。包含2个或以上SSR位点的序列有10569条。含2 nts重复单元的位点有20235个；含3 nts重复单元的位点有14403个；含4 nts重复单元的位点有6638个；含5 nts重复单元的位点有2411个；含6 nts重复单元的位点有1863个。单元重复次数主要集中在4～ 7次。包含不同碱基数重复单元的复合式SSR位点有7076个。由图5可知，包含2 nts的AT/AT是穿心莲转录组SSR中占比最高(27.7%)的重复单元，其次为AG/CT重复单元(12.8%)。

图5 SSR位点的数量占比

3 讨论与结论

穿心莲内酯作为穿心莲的主要二萜类次级代谢物质，具抗炎、抗感染及抗过敏等多重功效[20]。从转录水平分析穿心莲次生代谢途径的分子机制可帮助解析其药用价值。基于Pacbio平台的SMRT转录组测序无需组装即可捕获单个RNA分子的全长序列。相较于SGS技术，SMRT测序能够更准确地反应被测序物种转录组的实时信息，如发现更多的AS事件和剪切位点、发现新基因以及丰富基因组注释等[10]。本研究利用三代SMRT测序技术对穿心莲不同组织进行全长转录组测序分析，共得到高质量序列76812条。本次测序N50值为1414 bp，大于1000 bp，说明测序结果组装完整性较高[11]。功能注释结果显示，穿心莲蛋白主要参与翻译后修饰、信号转导及运输过程，体现次生代谢产物的前体在形成后的修饰和贮存过程。

TF是基因表达的开关，可通过与启动子结合来加速或缓和次生代谢进程[21]。在中药鼠尾草(Salvia Sclarea)中，SsWRKY18、SsWRKY40及SsMYC（bHLH家族成员）基因促进二萜类物质合成，过表达株系具细菌和真菌抗性[22]。番茄(Solanum lycopersicum)bZIP型TF可调控花青素积累[23]。中药黄芩(Scutellaria Baicalensis)的SbMYB8基因可促进黄酮类物质合成[24]。同时，腺毛体具有合成和分泌次生代谢产物的能力。烟草(Nicotiana Benthamiana)Cys2-His2(C2H2)锌指TF可正调控腺毛体的发育起始进程[25]。本研究对穿心莲全转录本信息进行TF家族聚类。结果表明，以上参与次生代谢途径的重要TF家族成员基因在穿心莲转录组中的含量占比位居前10。该结果为后续穿心莲次生代谢调控途径的TF挖掘提供了参考依据。

含内含子的前体mRNA可通过高度受控的AS过程从单个基因上产生多种转录本，从而增加胞内或胞间转录组及蛋白组的复杂性[26]。AS事件主要分为：外显子跳跃型、3’端可变型、5’端可变型及内含子保留型四大类。其中包括外显子互斥型（不同亚型跳跃保留不同的外显子）、可变启动子型（5’端内含子保留导致首个外显子不同的多种亚型）及可变终止子型（3’端内含子保留导致最后一个外显子不同的多种亚型）[27]。内含子保留型是植物中较为常见的AS事件，拟南芥(Arabidopsis Thaliana)和水稻(Oryza sativa)的内含子保留事件发生概率高达64%和55%[28]。Xu等[29]利用三代SMRT测序技术从丹参的根中鉴定出参与迷迭香酸和丹参酮生物合成的酶编码基因中的AS事件。本研究显示，穿心莲内酯二萜类前体GGPP的合成途径基因发生的主要AS类型为内含子保留，与前人研究一致。其中，1个GGPPS基因出现了可变启动子式内含子保留AS事件，产生6个首位外显子不同的AS亚型。该结果体现穿心莲全长转录组的蛋白编码多样性。与本研究类似，类胰岛素生长因子1基因的转录由两个启动子控制。启动子的交替使用配合其他位置的AS事件，使该基因产生若干第一外显子不同的mRNA亚型[27]。

SSR分子标记技术的核心序列为1～ 6个核苷酸的串联重复序列，具物种特异性，可用于种质资源鉴定[30]。本研究在穿心莲全长转录组中共找到45550处SSR位点；分布于27236条序列上；包含2个或以上SSR位点的序列有10569条；单元重复次数主要为4～ 7次；包含2 nts的AT/AT和AG/CT为优势重复单元。该结果有助于穿心莲的SSR分子标记开发。

本研究利用SMRT超长单分子测序技术来解析穿心莲幼苗的全长转录组信息，发现其体内功能基因主要参与次生代谢物质的合成与贮存相关的生物途径。活跃的次生代谢使其主要药效成分穿心莲内酯及其衍生物得到充分的积累。功能基因启动子通过结合直接促进次生代谢的转录因子以及转录后形成多种可变剪切亚型来维持上述次生代谢活性。此外，搜索到的穿心莲SSR位点可帮助学者对相关基因的遗传特性进行评估。综上，本研究为穿心莲内酯的合成机制研究及遗传多样性评价提供了分子基础。