锦绣龙虾全长转录组测序分析

周钱森，任宪云，史鲲鹏，于振兴，刘 萍，李 健

（1.水产科学国家级实验教学示范中心（上海海洋大学），上海 201306；2.中国水产科学研究院黄海水产研究所，农业农村部海洋渔业可持续发展重点实验室，青岛 266071；3.青岛海洋科学与技术试点国家实验室，海洋渔业科学与食物产出过程功能实验室，青岛 266071）

随着高通量测序技术的发展，转录组测序在基因表达水平分析和差异表达分析、新基因的挖掘、寻找单核苷酸多态性及应用、基因功能注释等方面都有所应用，在生物研究中发挥重要作用［1］。目前为止，大多数转录组测序研究利用以illumina平台为代表的二代测序技术，其测序长度远远低于真核生物RNA长度［2－4］。以PacBio平台为代表的第三代测序技术因其强大的读长优势，有效地解决了二代测序技术由于读长限制引起拼接不完整的问题［5］。PacBio平台单分子实时测序技术，又称SMRT（single molecule realtime）测序，SMRT cell可以在高GC含量区域完全跨过，确保序列的均匀覆盖度，并且DNA建库过程中，不需要进行PCR扩增，避免了PCR冗余和覆盖度不均一的问题，数据精准度可以达到99.9%［6－7］。PacBio SMRT测序借助其诸多优势，近些年在水产领域研究中得到广泛应用，张金勇等［8］对金乌贼（Sepiaesculenta）采用全长转录组测序获得高质量的生物学数据信息，筛选SSR位点并分析了其分布及组成特征。孔啸兰等［9］通过对竹筴鱼（Trachurusjaponicus）高通量测序，筛选获得27个具有多态性的微卫星标记，利用一个竹筴鱼群体对通过筛选的微卫星标记的种群遗传学特征进行评价。ZHANG等［10］通过三代与二代测序技术相结合挖掘施氏鲟（Acipenserschrenckii）性腺中早期配子发生的相关基因，更好地了解其生殖调控机制。本研究采用三代测序技术的PacBio SMRT平台对锦绣龙虾（Panulirusornatus）进行全长转录组测序，通过生物信息学方法对所测的序列进行拼接、分类、功能注释和代谢途径分析，获得锦绣龙虾丰富的序列信息，旨在为进一步挖掘锦绣龙虾相关功能基因、基因组学及开发分子标记等研究奠定基础。

锦绣龙虾隶属于节肢动物门（Arthropoda），甲壳纲（Crustacea），十足目（Decapoda），龙虾科（Palinuridae），龙虾属，在已知的19个种类的龙虾中锦绣龙虾是其中较为珍贵的优质水产品，其体型较大、营养丰富、脂肪含量低、味道鲜美［11－12］，主要分布在热带印度洋东部、东南亚、澳大利亚和西太平洋地区［13］，生长速度快于波纹龙虾（P.homarus）、中国龙虾（P.stimpsoni）、杂色龙 虾（P.versicolor）和 黄 斑 龙 虾（P.polyphagus）等［14］。作为宝贵的海洋渔业资源，锦绣龙虾因遭受高强度捕捞，资源量急剧减少。目前，锦绣龙虾育苗尚未取得成功，亟需开展人工繁殖相关研究以保护和增殖资源［15］。锦绣龙虾基因组测序亦未完成，其基因序列信息还比较匮乏，相关的分子遗传基础也很薄弱，因此本研究利用高通量测序技术对锦绣龙虾全长转录组测序，以期为后期相关功能基因的研究以及种质繁育等提供理论支持和数据支撑。

1 材料与方法

1.1 材料

实验所用锦绣龙虾购于海南省琼海市博鳌镇永贺生物科技有限责任公司，随机挑选3只体表完整无明显伤痕、活力好、规格整齐的龙虾进行取样，分别取鳃、肝胰腺、肌肉、胃、肠、脑和心脏组织，迅速放入液氮中保存，用于后续的RNA提取。

1.2 RNA提取与质控

采用Trizol法［16］分别提取锦绣龙虾鳃、肝胰腺、肌肉、胃、肠、脑和心脏组织总RNA，利用琼脂糖凝胶电泳分析RNA是否存在污染；通过Nanodrop测取OD260／280比值检测RNA浓度；Qubit精确定量RNA浓度；Agilent 2100判断RNA的完整性。选取符合标准的各组织RNA等量混合构建测序文库。

1.3 测序文库构建

将各组织混合后的RNA用Oligo（dT）的磁珠富集含有polyA的mRNA，然后通过SMARTer PCR cDNA Synthesis Kit反转录为cDNA，循环优化选取最合适的条件，通过PCR大规模扩增利用磁珠筛选的片段，获得一定数量的cDNA，全长cDNA进行末端修复、损伤修复、连接SMRT哑铃型接头，构建全长转录组文库。核酸外切酶消化，剔除cDNA左右端未连接的序列，最终通过DNA聚合酶结合引物，形成完整的文库。测序由北京诺禾致源科技股份有限公司完成。

1.4 全长转录组测序及序列分析

采用SMRTlink（PacBio）处理原始下机数据获得子序列，校正子序列获得环形一致性序列（circular consensus sequence，CCS），然后根据环形一致性序列是否含有3′端引物、5′端引物以及polyA尾将其分为全长序列与非全长序列；再对全长序列进行聚类，获得聚类一致序列［17］；最后对得到的全长序列进行改良，获得高质量的一致序列用于后续分析。

1.5 全长转录组序列的功能注释及结构分析

利用CD-HIT软件对获得的改良一致序列去冗余，将锦绣龙虾组装所得单基因簇与公共数据库进行比对，通过序列相似程度得出具有最高相似性的蛋白，从而对该单基因簇进行功能注释，利用NCBI非冗余蛋白序列数据库（nonredundant protein sequences，NR）、NCBI核酸序列数据库（nucleotide sequences，Nt）［18］、蛋白质真核同源数据库（eukaryotic orthologous groups，KOG）［19］、蛋白质序列数据库（SwissProt protein sequence database，SwissProt）［20］、蛋白质家族数据库（protein families database，Pfam）［21］、基因功能描述分类系统（gene ontology，GO）［22］及京都基因与基因组百科全书（Kyoto encyclopedia of genes and genomes，KEGG）［23］。

2 结果与分析

2.1 全长转录组测序与组装分析

2.1.1 测序结果与数据组装

采用PacBio Sequel测序平台对锦绣龙虾鳃、肝胰腺、肌肉、胃、肠脑和心脏等组织混样进行全长转录组测序，共获得39 828 440个子序列（71.1 Gb），平均子序列长度为1 786 bp，N50为2 492 bp。通过每个ZMW孔中的子序列得到CCS聚类之后得到的序列数为1 018 599个，序列平均长度为2 377 bp，N50为2 872 bp。同时含有3′引物和5′引物，以及3′引物前含有polyA尾的全长序列（full-length，FL）554 600个，全长非嵌合序列（full-length non-chimeric read，FLNC）543 552个，序列平均长度为2 055 bp，FLNC／CCS为53.36%。全长转录组得到改良后一致序列28 034个，序列平均长度为2 250 bp。

2.1.2 CDS预测

利用ANGEL软件对获得的基因片段进行蛋白编码区预测分析［24］，结果显示一共有12 660个基因片段可视为蛋白编码区，其序列长度的范围为0～5 000 bp，主要集中于250～1 250 bp（图1）。

图1 CDS长度分布图Fig.1 Statistics of sequence length of CDS

2.1.3 SSR预测

简单重复序列标记（simple sequence repeats，SSR），又称为短串联重复序列或微卫星标记。是一类由几个核苷酸为重复单位组成的长达几十个核苷酸的重复序列，长度较短，且广泛均匀分布于真核生物基因组中。通过检索Unigene序列，共发现14 277个SSR位点，分布在8 813个Unigene中。从表1可以发现，单核苷酸到三核苷酸重复最多，占总SSR位点数量的99.13%，其中单核苷酸重复以A／T重复基序最多；二核苷酸重复以AC／GT和AT／AT重复基序出现频率最高；三核苷酸重复则以ATT／AAT和ATC／GAT为主要类型；四核苷酸重复则以GAAA／TTTC较多。其他五核苷酸和六核苷酸重复基元类型较多，数量非常少。

表1 锦绣龙虾SSR不同重复基元分布情况Tab.1 Distribution of different repeat motifs in P.ornatus

2.1.4 lncRNA分析

长链非编码RNA（lncRNA）是一类转录本长度超过200 nt、不编码蛋白质的RNA分子。由于建库原理的限制，我们只能获得含有polyA尾的lncRNA。使用CNCI［25］、CPC2［26］、Pfam［27］以及PLEK［28］对CD-HIT去冗余得到的基因进行编码潜能预测，经数据库比对后预测其编码潜能，最终分析6 315个LncRNA序列，其中共有数目为1 342个（图2）。

图2 编码潜能预测维恩图Fig.2 Venn diagram of encoding potential prediction

2.2 Unigene的功能注释

2.2.1 注释结果统计

将锦绣龙虾转录组所获得的单基因簇序列在NR数据库中进行对比，共有11 086个Unigene获得注释，对比到424个物种，其中对比较多的前20个物种包括钩虾（Hyalellaazteca，3 907个）、湿木白蚁（Zootermopsisnevadensis，504个）、大型蚤（Daphniamagna，483个）、凡纳滨对虾（Litopenaeusvannamei，296个）、美洲鲎（Limulus polyphemus，266个）、斑节对虾（Penaeusmonodon，257个）、克氏原螯虾（Procambarusclarkii，196个）、中华绒螯蟹（Eriocheirsinensis，182个）、松叶蜂（Neodiprionlecontei，166个）、拟穴青蟹（Scylla paramamosain，160个）、淡水枝角水蚤（Daphnia pulex，139个）、日 本 囊 对 虾（Marsupenaeus japonicus，133个）、鸭嘴舌形贝（Lingulaanatina，121个）、中国明对虾（Fenneropenaeuschinensis，116个）、红螯螯虾（Cheraxquadricarinatus，115个）、美洲螯龙虾（Homarusamericanus，114个）、白氏文昌鱼（Branchiostomabelcheri，102个）、三疣梭子蟹（Portunustrituberculatus，85个）、罗氏沼虾（Macrobrachium rosenbergii，83个）、淡 水 螯 虾（Pacifastacusleniusculus，81个）。从NR数据库的注释结果来看，注释到钩虾的基因最多，推测锦绣龙虾与钩虾同源性较高。随着锦绣龙虾生物学研究的不断深入，GenBank的基因数据越来越丰富，后期将获得更多锦绣龙虾基因注释。

2.2.2 GO功能分类

GO是一套国际标准化的基因功能描述的分类系统，分为细胞组分（cellular component，CC）、分子功能（molecular function，MF）、生物学过程（biological process，BP）3大类［22］。将获得的锦绣龙虾Unigene与GO数据库进行匹配，分别有17 624个、10 398个和9 530个被划分到BP、CC及MF 3大类，涵盖上述功能类别的51个亚类（图3），其中BP 24个亚类，涉及细胞进程、代谢进程及单生物进程的较多，分别有3 768个、3 456个和2 918个；CC 18个亚类，涉及细胞、细胞部分及细胞器的较多，分别有1 894个、1 894个和1 428个；MF 9个亚类，涉及结合活性和催化活性的较多，分别有4 789个和3 398个。

图3 Unigenes的GO功能分类图Fig.3 GO functional categories of Unigenes

2.2.3 KOG分类统计

COG（cluster of orthologous groups of proteins）是根据细菌、藻类和真核生物完整基因组的编码蛋白系统进化关系分类构建而成的蛋白数据库，其中真核生物数据库称之为蛋白质真核同源数据库（eukaryotic orthologous groups，KOG）［19］。通过KOG数据库对比到锦绣龙虾单基因簇，共有9 433个获得注释信息，共分为26个功能组分（图4），其中获得注释较多的功能组分一般功能预测（1 413个，14.98%）、其次是信号转导机制（1 030个，10.92%）及翻译后修饰、蛋白转换、伴侣蛋白（917个，9.72%），其中未知蛋白仅有4个。

图4 Unigenes的KOG功能分类Fig.4 KOG function classification of Unigenes

2.2.4 KEGG功能注释

KEGG（Kyoto encyclopedia of genes and genomes）是系统分析基因产物和化合物在细胞中的代谢途径以及基因产物功能的数据库，已建立了一套完整KO注释的系统，可完成新测序物种的基因组或转录组的功能注释［29］。与KEGG数据库对比，最终锦绣龙虾单基因簇注释到5大类34小类（图5），其中基因数量较多的代谢通路有信号转导机制通路（1 015个），转运和分解代谢通路（726个），跨膜转运的代谢通路基因最少，仅有9个。这些获得注释的单基因簇可为后续的锦绣龙虾功能基因研究和应用提供数据基础。

2.2.5 转录因子分析

转录因子（transcription factor，TF）作为一类特殊的DNA结合蛋白，能与基因5′端上游的特定序列专一性结合，从而使目的基因能够在特定的时空进行表达，通过转录因子间以及与其他相关蛋白间的相互作用达到激活或抑制转录的效果，发挥着重要的调控作用。动物转录因子鉴定使用动物转录因子数据库［30］animalTFDB 2.0，预测到转录因子有415个，隶属于29个转录因子家族（图6），其中比较多的转录因子家族：zf-C2H2家族有109个、ZBTB家族有75个，E2F家族、NFYB家族、THR-like家族、HSF家族、IRF家族、MBD家族和STAT家族最少，均只有2个，这些转录因子家族成员的获取可为后期锦绣龙虾生长发育、代谢调节、免疫应答等相关研究奠定基础。

图6 锦绣龙虾转录因子分析Fig.6 Transcription factor analysis of P.ornatus

3 讨论

PacBio Sequel测序仪是美国太平洋生物技术公司（Pacific Biosciences）基于PacBio RSII平台最新推出的第三代测序平台，该系统以其测序通量高、单Gb数据成本低、周期短等特点广受青睐。PacBio第三代测序技术的最大特点是无需进行PCR扩增，可直接读取目标序列，是转录组从头测序的首选［31］。对于无基因组参考的水生生物来说，运用三代测序技术进行转录组分析以及功能基因的挖掘是较好的策略［32－34］。本研究对锦绣龙虾全长转录组测序及分析，共获得39 828 440个子序列，平均子序列长度为1 786 bp，N50为2 492 bp，结果显示转录组所获得的组装序列完整性较好。利用NR、Nt、KOG、SwissProt等7大公共数据库进行功能注释分类共获得序列28 034个，有11 086个获得注释，对比到424个物种，有10 463个注释到350个代谢通路，有9 433个获得KOG注释，分为26个功能组分。曾地刚等［35］对凡纳滨对虾肝胰腺采用二代转录组测序，结果 获得500 177个EST（expressed sequence tags），序列平均长度363 bp。拼接获得20 225个Unigene，序列平均长度507 bp。注释到NR、COG以及KEGG数据库的分别为13 676个、4 645个、4 104个。李喜莲等［36］以红螯螯虾肝脏、精巢以及卵巢为材料进行二代转录组测序，获得了6 736个单基因簇，注释到GO的为16 989个，注释到COG的为4 697个，注释到KEGG的为9 842个。对比研究结果显示，三代测序获得的序列质量、基因数以及注释到的基因信息均优于二代测序。

基于PacBio测序平台对构建甲壳动物转录组文库序列聚类和改良，ZHANG等［37］通过凡纳滨对虾全长转录组文库获得72 648条高质量序列，WANG等［38］通过中国明对虾进全长转录组文库获得10 795条高质量序列，POOTAKHAM等［39］通过斑节对虾全长转录组文库得到22 418条高质量序列。本研究通过锦绣龙虾全长转录组文库最终获得了13 297条高质量序列，相较于凡纳滨对虾和斑节对虾少，除了物种的差异外，也归因于转录组聚类过程中的参数设置。同时，测序过程中5′末端的降解也可能导致同一转录组序列进行不同的聚类划分。考虑到这一点，本研究将锦绣龙虾全长转录组序列聚类，获得高质量的一致序列用于后续分析。

本研究利用Nr、KEGG等数据库对锦绣龙虾转录组序列进行功能注释，结果显示有28 034个单基因簇获得注释信息，部分基因序列仍没有获取注释，可能由于锦绣龙虾是未测序物种，公共数据库中缺乏相关的功能注释信息，也可能是部分锦绣龙虾单基因簇序列较短造成。将获得的蛋白同源序列与NR数据库对比发现，与钩虾对比的同源信息最多，占35%，推测可能是由于钩虾与锦绣龙虾的进化史和生活史较为相似。将锦绣龙虾转录组单基因簇注释到GO数据库，可划分为3大类共计51个分支，与KOG数据库对比可分为26类，GO和KOG注释功能分类的结果显示锦绣龙虾单基因簇的功能涉及了各类生命活动。KEGG注释到5大类34小类，其中涉及代谢通路和次生物质的生物合成基因较多。从转录组共检测出14 277个SSR位点，这些SSR在提高物种遗传多样性潜能方面发挥着重要的作用，为下一阶段开发锦绣龙虾多态性SSR分子标记提供了基础数据。通过对锦绣龙虾转录组进行注释分析，初步阐明了锦绣龙虾基因的功能、参与的生物过程、所在的代谢途径或信号通路等，对于后期深入研究关键基因的功能提供了信息，为发掘锦绣龙虾功能基因、研究相关生理功能奠定了基础。