基于PacBio平台的七叶一枝花全长转录组测序

宋发军，杨瑞霜，吕昕芮，刘祖懿，王红莹，孟艳艳

（中南民族大学 生命科学学院& 生物技术国家民委重点实验室& 武陵山区特色资源植物种质保护与利用湖北省重点实验室，武汉 430074）

重楼是我国名贵中药材，重楼皂苷是其主要活性成分［1］. 因其药用价值和经济价值高，且重楼药材市场的供需矛盾突出，当前重楼属多个物种资源已濒临枯竭. 同时，由于重楼的药用部位生长周期长，种子深度休眠，加之重楼皂苷的生物合成途径仍不甚清楚，完全依赖于植物中提取，因此加剧了重楼属植物的资源危机［2］.为此，研究者们从分子角度对重楼属植物展开了相关研究，以期为该属植物资源的保护和可持续利用提供科学参考［3］. 然而现阶段重楼属植物中关于重楼皂苷合成、种子休眠与萌发等的基因信息仍较少.这严重阻碍了人们对上述问题的解析. 故而，快速获得重楼属植物的大量基因信息并研究其功能特征，从而解析重楼皂苷的合成途径、种子萌发等机制，为后续人工调控重楼皂苷的生产与合成奠定基础，显得至关重要.

近年来，随着高通量测序技术的快速发展，转录组测序技术在获取大量基因信息、新基因挖掘、基因家族鉴定、代谢途径及系统进化关系的分析等方面发挥着重要作用［4］.前人也利用转录组测序技术分别以重楼属植物的药用组织、不同萌发状态的种子等为材料，展开了相关研究［5-8］.这些研究均为二代测序技术，且主要以滇重楼为对象，而对另一个重楼正品药源植物——七叶一枝花的研究较少.

在诸多测序技术中，第三代PacBio 技术在转录组测序方面有着更好的应用前景，它在植物优良品种选育、分析次生代谢产物生物合成途径等研究中发挥着重要作用，尤其适用于基因信息资源匮乏的药用植物领域［9-11］. 这是因为作为最新一代的测序技术，其最大读长可达40 kb，且结果不受序列特异性影响，也不会产生由PCR 过程中碱基突变等原因造成的测序错误［12］. 另外，利用PacBio 测序技术还可以检测到多种DNA 修饰如5mC、5hmC、m6A、硫代磷酸等［12-14］. 上述诸多优点，使得PacBio 技术在对无基因组信息物种的研究上有着极大的优势，该技术已经被应用于多个药用植物的相关研究［10，15］.

本研究以中药重楼正品药源之一的七叶一枝花为研究材料，进行了PacBio全长转录组的测序工作，并对所得的转录本进行分析了和注释.本研究将为重楼属植物的分子研究提供数据支撑和基因资源，同时为重楼皂苷合成相关基因的筛选及其种子休眠机制等方面的研究奠定基础.

1 材料与方法

1.1 实验材料

2016 年4 月在湖北省巴东县采集的4 年生七叶一枝花的5 个组织（根状茎、茎、叶、花、果荚）及2016 年秋季采收的种子.每种样品采集3份，其中新鲜根状茎洗净后切成小块立即进行生物反应灭活处理（液氮冷冻），其他部位直接进行灭活处理，于-80 ℃冻存.

1.2 RNA提取、反转录及测序

分别提取根状茎、茎、叶、花、果荚和种子6个部位的RNA［RNeasy Plus Mini Kit（#74134），Qiagen公司］，等量混合后，使用Clontech SMARTer PCR cDNA 试剂盒合成第一条链cDNA，再进行PCR 优化，并利用BluePippinTM将测序数据按照长度分为1～2 kb、2～3 kb、3～6 kb、5～10 kb 四个库，进行扩大PCR 反应（若转录本大于3～6 kb 则需要利用BluePippinTM进行大小筛选）后即可进行转录组测序.转录组测序工作委托深圳华大基因公司完成.

1.3 Reads Of Insert序列及全长转录本的获取

Sequel 测序 平 台共 包 含16 个SMRT cell，每 个SMRT cell包含100万个ZMW（zero-mode waveguides），每个ZMW 为一个测序单元.同一ZMW 中的所有subreads 来自同一个转录本，其中碱基出错率随机，通过subreads间比对以提高碱基质量，并利用SMRT软件获取Reads Of Insert.使用pbtranscript.py脚本识别全长转录本序列，检测这些序列是否包含5′端引物、3′末端引物和polyA 尾巴，区分全长序列和非全长序列.

1.4 高质量一致性序列（isoform）的获取和校正

同一孔中所有subreads 先进行矫正得到环形一致性序列，然后孔与孔之间矫正得到高质量一致性序列，即isoform.利用SMRT 进行isoform 水平的聚类，使用非全长序列对聚类的isoform 进行Quiver 质量校正.使用cd-hit分别对高质量（QV＞99%）和低质量（QV＜99%）isoform 进行去冗余分析，计算覆盖度并对其序列长度及分布进行统计.

1.5 转录本序列的功能注释、CDS预测及SSR检测

利 用NR、NT、GO、COG、KEGG、SwissProt、InterPro七个数据库对得到的转录本序列进行功能注释，并根据聚类到相同isoform 的序列数量，计算isoform 覆盖度.随后选取与数据库匹配度较高的片段，进行Blast注释以预测CDS序列，而未能进行功能注释的片段则利用EST Scan来预测CDS序列.然后设计特定引物，通过PCR技术检测七叶一枝花简单序列重复（Simple sequence repeat，SSR）的物种特异性.

2 结果与分析

2.1 测序结果及数据质量控制

基于Pacific Biosience RS Ⅱ平台，按照测序片段长度构建4 个库进行测序.如表1 所示，整合4 个库的数据，共获得的1219115492 个碱基分布于357362 个reads 中，QC（Quality control）值在86%～93%之间.由于PacBio 测序错误分布随机，一般情况下高质量的文库读段更长，本研究构建的文库中有较多3 kb 以上的片段，且短读段的QC 值较高，由此可以看出此次转录组测序结果较好，可为后续的数据组装提供很好的原始数据.

表1 测序数据概览Tab.1 Sequencing data overview

2.2 高质量isoform的获取及分析

整合4个库共获得58763个isoforms，去除冗余后进行质量分析，筛选出QV＞99%的高质量isoforms共52537个（表2），所得序列的相关信息均已提交在NCBI Bioproject 数据库中（https：//www.ncbi.nlm.nih.gov/bioproject/，BioSample accession：SAMN09762366）.

表2 七叶一枝花中isoform的数量Tab.2 Number of isoform in Paris polyphylla Smith var. chinensis

对52537 个高质量isoforms 进行分析（表3），平均长度为2607 bp，这远大于一代测序（600 bp）和二代测序（200 bp）长度.其N50 达到2998 bp，表明本次测序结果中序列拼接质量较好.对高质量isoforms的序列长度进行统计，可见这些isoforms 主要聚集于1800 bp和3000 bp（图1）.

图1 七叶一枝花中高质量isoform长度分布直方图Fig.1 Histogram for high quality isoform length in Paris polyphylla Smith var. chinensis

表3 七叶一枝花中高质量isoform数Tab.3 Number of high quality isoform in Paris polyphylla Smith var. chinensis

综上所述，本研究获得的高质量isoform 普遍为较长的序列，说明文库质量较好，而短序列质量相较长序列高，进一步从侧面反映出本次构建isoform文库整体质量优良，可为后续基因筛选及相关分析提供更好的参考.

2.3 转录本的功能注释及覆盖度分析

使 用NR、NT、GO、COG、KEGG、SwissProt、InterPro 共7 个数据库，对所得转录本进行功能注释，在52537 个高质量isoforms 中总计有40725 个被成功注释，涵盖了七叶一枝花所有高质量isoforms的77.52%.由图2 可知，在5 个数据库中共同检索到的isoforms 有11901 个，在NR 中单独检索到3878 个isoforms，在KEGG 中单独检索到52 个isoforms，在COG中单独检索到的isoforms有5个，在SwissProt中单独检索到的为38 个，在InterPro 中单独检索到isoforms139个.

图2 七叶一枝花中isoform注释结果的韦恩图Fig.2 Venn diagrams of isoform annotation results in Paris polyphylla Smith var. chinensis

其中，COG 功能注释的结果显示（表4），参与复制、重组和修复的isoforms 最多，达5094 个；属于一般功能预测的有3815 个.此外，还有少部分isoforms分布于防御机制、胞外机构等途径.与此同时，还发现了653 个分布于次级代谢产物的生物合成、转运与代谢等途径的isoforms，这些isoforms 极有可能与重楼皂苷合成途径相关.

表4 七叶一枝花中isoform的COG注释分布Tab.4 COG annotation distribution of isoform in Paris polyphylla Smith var. chinensis

GO 注释中的isoforms 共分为生物学过程（蓝色）、细胞组成（绿色）和分子功能（红色）3 大类（图3），49个小类.在生物学过程这一大分类分布于代谢过程的isoforms达2724个，居首位；在细胞组成中分布于细胞和细胞部分的isoforms 最多，均为2367 个；而在分子功能这一大分类中，大多数isoforms被归为催化活性一类，达2606个.

图3 七叶一枝花中isoform的GO注释分布图Fig.3 GO annotation distribution map of isoform in Paris polyphylla Smith var. chinensis

KEGG 注释中，52537 个isoforms 分布于碳水化合物的代谢、翻译、折叠、分类和降解、脂质代谢等分类中.为进一步了解基因的生物学功能，对注释到的isoforms 进行代谢路径分析，结果见表5，共有31934 个isoforms 分属于134 个代谢通路.其中，前3个代谢通路分别为代谢通路、次生代谢物的生物合成、剪接体，其注释到的isoforms 数分别为7412、4631、2235 个.此外，ABC 转运体、RNA 转运等代谢途径中，也有大量的isoforms富集.

表5 七叶一枝花KEGG分析中部分代谢通路Tab. 5 Some metabolic pathways in KEGG analysis of Paris polyphylla Smith var. chinensis

2.4 CDS预测

为进一步了解所得isoform 的基因信息，利用Blast 和EST Scan 注释预测其CDS 序列，结果如表6所示.本研究最终获得39343条CDS序列，平均长度759 bp.对预测的CDS 序列长度进行整合，结果如图4. 这些CDS 序列长度分布在200 nt 到3000 nt 之间，其中长度≥500 nt 及≥1000 nt 的序列占总序列数的比例分别为70.45%和30.82%，这表明本研究的测序结果较为准确可靠，并且所获得的长片段isoform 可直接用于基因的全长克隆及功能分析等研究.

图4 七叶一枝花CDS序列长度分布Fig.4 CDS sequence length distribution in Paris polyphylla Smith var. chinensis

表6 七叶一枝花中CDS预测结果Tab.6 CDS forecast results in Paris polyphylla Smith var. chinensis

2.5 SSR检测

转录本聚集后，继续对所得序列进行SSR 分析（图5）.SSR 检测中双核苷酸重复的含量最高，达到9925 个，占全部SSR 位点的46.30%；其次是单核苷酸重复，为8360 个.此外三核苷酸重复为2498 个，占总SSR 位点的37.54%.该项结果为后续七叶一枝花及重楼属植物分子标记的开发提供了依据.

图5 七叶一枝花SSR检测结果Fig.5 SSR test results in Paris polyphylla Smith var. chinensis

3 讨论

通过转录组测序技术可获得大量的转录本信息，这对于七叶一枝花等没有基因组信息的物种而言，是获取其基因序列和功能信息的重要手段. 因此，近年来人们对重楼属植物的转录组研究逐渐增多.

但是对比相关研究发现，大多数重楼属植物的测序研究主要以云南重楼为对象，对七叶一枝花的研究较少，并且由于测序目的不同，选取的组织和采用的测序技术也有所差异.前人研究中选取的大多为根状茎或种子等单一组织，重楼属植物中果荚和花等组织的测序研究尚未见报道.本研究利用七叶一枝花的根状茎、茎、叶、花、果荚及种子6个组织进行混合测序，所得的基因信息系统和全面.并且，已有的研究均采用二代测序技术，虽然获得的unigenes 数量较多，但大多为200～300 bp 左右的短片段，难以拼接和有效使用.如LING 等人为研究重楼种子的休眠机制，对云南重楼种皮和种子进行Illumina 测序，共鉴定得到146671 个平均长度为923 bp 的unigenes［16］.LIU 等对云南重楼4 年生根和8 年生根进行NGS 测序，总共获得87577 个平均长度为614 bp 的unigenes［17］. 相比二代测序所得的短片段序列信息，本研究使用PacBio 测序获取的isoform 大多数为1800～3000 bp 的片段，平均长度为2607 bp.因此，本研究利用PacBio 测序技术所获得的长片段的isoforms，可直接用于基因全长预测与克隆及功能鉴定等实验，有效减少了扩增目标基因全长等工作量，并可为二代测序数据的拼接提供参考模板，为后续研究提供了极大的便利.

转录本注释中，本研究所得的isoforms 有77.52%在NR，NT等数据库中得到注释，在LI等［18］的滇重楼基因组测序分析结果中，转录本长度多为500～900 nt，有94.36%的序列得到注释；两者都在NR 数据库中得到了最大占比的注释. 另外，KEGG分析中，本研究所得的isoforms更多的注释到代谢通路、次生代谢物的生物完成、剪接体等通路，而LI 的研究中，更多的isoforms被注释到代谢通路、RNA 转运、内质网中蛋白质加工等通路［18］.这可能是因为，两个研究工作中所选用的测序方法不同，因此获得转录本的数量和类型都产生了差异，也可能跟实验材料不同有一定的关系. 此外，与LIU 等的研究结果［19-20］对比，虽然本研究获取并注释到的isoform 总量 略 少（40725/65535），但 注 释到NR、SwissProt、KEGG 等数据库的isoforms 数量更多，且KEGG 和GO 的注释功能分布中各类别占比基本一致.综上，本研究所得的isoforms 数量和注释的质量都较为可靠，为后期相关研究中的基因挖掘、筛选和功能验证提供了有效资源.

植物种子休眠的解除伴随着胚乳中大量贮藏蛋白质的降解，同时与DNA 复制、糖酵解、信号转导、能量代谢等途径有关［19-20］. 此外，重楼皂苷的合成途径主要分为异戊烯焦磷酸的合成、甾体碳骨架的形成、甾体皂苷元的形成与修饰以及甾体皂苷元糖基化重楼皂苷等4部分［21-22］.本研究不仅获得了大量与糖类代谢、信号转导、能量转换等与种子萌发有关的转录本，同时也富集到了大量参与次生代谢物的生物合成及降解相关的isoform，对上述isoforms进行更加详细的筛选和分析，将为深入研究重楼皂苷合成相关基因以及重楼种子萌发等问题奠定良好的基础.

由图5 可知，SSR 检测结果中双核苷酸重复占据优势，其次为单核苷酸重复，而三苷酸重复较少，这与张成才、李俊仁等人的研究结果有所差异［8］，这可能是由于物种差异性和实验样本的不同而造成的. 这些差异的SSR 位点，将进一步为重楼属植物的遗传多样性分析、目标基因标定、遗传图谱构建和分子标记辅助育种等提供基础. 另外，由于PacBio 测序技术获得的转录本长度较长，更接近真实的转录本分布.本研究所预测到的CDS 长度分布能较为真实地显示七叶一枝花中CDS 的长度分布情况，这为后续重楼属植物的基因研究提供了重要参考.

综上，本研究利用PacBio 测序技术获取的七叶一枝花全长转录组测序的数据质量较为可靠和全面，该结果进一步丰富了七叶一枝花乃至整个重楼属植物的基因资源，为重楼属植物的二代测序数据拼接提供了良好的参考模板，同时也为解析重楼皂苷合成途径、种子萌发机制等问题提供了良好的依据.