海洋与湖沼  2023, Vol. 54 Issue (5): 1463-1475   PDF    
http://dx.doi.org/10.11693/hyhz20230300058
中国海洋湖沼学会主办。
0

文章信息

闫晗, 惠敏, 沙忠利, 程娇. 2023.
YAN Han, HUI Min, SHA Zhong-Li, CHENG Jiao. 2023.
深海化能极端环境中劳盆拟刺铠虾(Munidopsis lauensis)的全长转录组测序分析
SMRT SEQUENCING OF FULL-LENGTH TRANSCRIPTOME OF THE MUNIDOPSID SQUAT LOBSTER MUNIDOPSIS LAUENSIS ENDEMIC TO DEEP-SEA CHEMOSYNTHETIC ENVIRONMENTS
海洋与湖沼, 54(5): 1463-1475
Oceanologia et Limnologia Sinica, 54(5): 1463-1475.
http://dx.doi.org/10.11693/hyhz20230300058

文章历史

收稿日期:2023-03-09
收修改稿日期:2023-04-23
深海化能极端环境中劳盆拟刺铠虾(Munidopsis lauensis)的全长转录组测序分析
闫晗1,3, 惠敏1,2, 沙忠利1,2, 程娇1,2     
1. 中国科学院海洋研究所海洋生物分类与系统演化实验室 青岛市海洋生物多样性与保护重点实验室 山东青岛 266071;
2. 崂山实验室 山东青岛 266237;
3. 中国科学院大学 北京 100049
摘要:劳盆拟刺铠虾(Munidopsis lauensis)是一种十足目铠甲虾, 能够适应深海化能合成生态系统的极端环境, 但是其已知基因序列信息却很少。为此, 采用PacBio平台的SMRT测序技术对来自南海台西南冷泉劳盆拟刺铠虾的肝胰腺、鳃、肌肉和肠混合组织进行了三代全长转录组测序。通过聚类、校正、去冗余之后共获得28 811条高质量isoform, 平均长度为2 086 bp, N50长度为2 275 bp。使用Nr、SwissPort、KEGG、KOG数据库对转录本序列进行功能注释, 共有20 616 (71.56%)条isoform得到注释。进一步高级注释共获得21 848个蛋白编码序列, 共预测到537个转录因子、6 430个长链非编码RNA和10 060个SSRs位点。KEGG通路分析发现两条可能与劳盆拟刺铠虾适应深海化能极端生境相关的通路, 分别为过氧化物酶体通路和谷胱甘肽代谢通路。其中, 共有180条isoform被发现参与编码过氧化物酶体通路中的31个关键酶, 213条isoform参与谷胱甘肽代谢通路中19个关键酶的编码。基于全长转录组数据, 利用同源比对、功能结构域预测及系统发育树构建等方法对劳盆拟刺铠虾谷胱甘肽S-转移酶(GST)基因家族进行挖掘与鉴定, 共发现属于theta、delta、mu、kappa四种亚家族的20条GST候选序列。通过对劳盆拟刺铠虾全长转录组测序、功能注释和环境适应相关的重要基因及通路的分析, 不仅丰富了深海组学数据资源, 也为进一步揭示甲壳动物适应深海化能极端环境的分子机制奠定了基础。
关键词劳盆拟刺铠虾    深海化能生态系统    全长转录组    功能注释    极端环境适应    
SMRT SEQUENCING OF FULL-LENGTH TRANSCRIPTOME OF THE MUNIDOPSID SQUAT LOBSTER MUNIDOPSIS LAUENSIS ENDEMIC TO DEEP-SEA CHEMOSYNTHETIC ENVIRONMENTS
YAN Han1,3, HUI Min1,2, SHA Zhong-Li1,2, CHENG Jiao1,2     
1. Laboratory of Marine Organism Taxonomy and Phylogeny, Qingdao Key Laboratory of Marine Biodiversity and Conservation, Institute of Oceanology, Chinese Academy of Sciences, Qingdao 266071, China;
2. Laoshan Laboratory, Qingdao 266237, China;
3. University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: The munidopsid squat lobster Munidopsis lauensis (Baba and de Saint Laurent, 1992) is endemic to deep-sea chemosynthetic environments. However, no reference genome or transcriptome is available for M. lauensis, which impedes the studies of evolution and adaptation studies of the species. Therefore, we aimed to obtain the full-length (FL) transcriptome of M. lauensis using PacBio single molecule real-time sequencing (SMRT) technology. The total RNA extracted from hepatopancreas, gill, muscle, and intestinal tract of the species was mixed for SMRT sequencing. After isoform-level clustering, polishing, and redundancy removing, we identified 28 811 high-quality isoforms with mean length value of 2 086 bp and N50 value of 2 275 bp. A total of 20 616 (71.56%) isoforms had at least one significant hit in the Nr, SwissProt, KOG, or KEGG databases. We further predicted 21 848 coding sequences, 537 transcription factors, 6 430 long non-coding RNAs and 10 060 simple sequence repeats. Two pathways possibly associated with deep-sea extreme environment adaptation in M. lauensis were identified, including Peroxisome pathway with 31 key enzymes encoded by 180 isoforms and Glutathione metabolism pathway with 19 key enzymes encoded by 213 isoforms. A total of 20 GST-encoding isoforms, belonging to four GST classes (theta, delta, mu, and kappa), were identified from M. lauensis FL transcriptome by examining multiple-sequence alignments for conserved functional domains and by performing phylogenetic analyses. Overall, this study provided the first large-scale FL transcriptome in M. lauensis, allowing for functional, structural, and genomics studies of this endemic crustacean species found in deep-sea chemosynthetic ecosystems, and enriched the data resources of deep-sea genomics for further studies on molecular mechanisms of deep-sea extreme environment adaptation in crustaceans.
Key words: Munidopsis lauensis    deep-sea chemosynthetic ecosystem    full-length transcriptome    functional annotation    extreme environment adaptation    

自1977年美国深海潜水器阿尔文号在加拉帕戈斯裂谷附近发现深海热液活动和高密度的热液生物群落以来, 深海化能合成生态系统开始展现于世人面前(Corliss et al, 1979)。随后, 1983年美国研究人员在广阔的墨西哥湾3 200 m深处首次发现了海底冷泉系统(Anderson et al, 1983)。与依赖太阳光能的光合生态系统不同, 以热液和冷泉为代表的深海化能合成生态系统主要依靠化能自养微生物利用硫化氢(H2S)、元素硫(S)、甲烷(CH4)和氢气(H2)等还原性物质氧化产生的能量进行固碳、合成有机物, 从而维持生态系统的运转, 这种“黑暗食物链”的发现改变了人们对深海初级生产力主要来源的认知(Dubilier et al, 2008)。此外, 深海化能生态系统是认识特殊生命过程和研究生物适应性进化的天然实验场所。深海化能环境具有黑暗、高压、低氧、寡营养且富含高浓度硫化物、甲烷和重金属等特征, 大型生物却能在如此极端环境中生存并形成群落, 主要原因是它们经过长期的适应与历史进化过程, 形成了独特的生命特征和环境适应机制(程娇等, 2021)。因此, 对深海化能合成生态系统中生物的环境适应机制的研究, 不仅将有助于加深对深海极端环境中特殊生命过程的理解, 也将为深海生物资源的开发利用奠定理论基础。

近年来, 高通量测序技术的快速发展及大数据分析能力的提高, 使得在组学水平探究深海生物演化和环境适应机制成为可能。作为后基因时代应用最为广泛的技术手段, 转录组测序技术可以全面快速获得生物体在特定时空或生理条件下几乎所有转录本的序列信息和表达信息, 已经在深海生物适应性演化机制研究中逐步得到应用。例如, 在双壳纲贝类中, Zheng等(2017)通过与浅海贻贝的比较转录组分析发现两种深海偏顶蛤(Bathymodiolus platifronsB. manusensis)中与共生菌识别相关的免疫识别受体发生收缩, 物质转运、溶酶体活动以及细胞凋亡等相关基因受正选择或高表达, 可能有助于共生关系的建立和维持; 在多毛类环节动物中, Zhang等(2017)对分布于热液区深海多鳞虫(Branchipolynoe pettiboneaeLepidonotopodium sp.)与浅海近缘物种进行了比较转录组分析, 发现热液多鳞虫可能通过不同类型的血红蛋白的快速进化或高表达等策略实现对深海低氧环境的适应; 在十足目甲壳动物中, Cheng等(2019)对冲绳热液与南海冷泉共有优势物种柯氏潜铠虾(Shinkaia crosnieri)进行转录组测序和进化分析, 提出柯氏潜铠虾很可能是通过功能基因的适应性进化和改变基因的表达水平来适应热液与冷泉两种不同的深海化能环境。以上研究结果表明转录组学分析方法已成为探究深海极端环境中特殊生命过程的重要研究手段。然而, 由于二代测序读取长度短, 转录本易发生拼接错误, 无法准确获取完整转录本, 可变剪切和转录本多样性在很大程度上仍然是未知(Zhang et al, 2020), 这限制了深海生物环境适应机制及深海基因资源开发方面的深入研究。

单分子实时(single molecular real-time, SMRT)测序技术是近年逐渐成熟的新一代测序技术, 具有超长读长、无GC偏好性、无PCR扩增偏向性、数据覆盖度与准确度高等特点(Sharon et al, 2013)。基于PacBio平台的全长转录组测序技术无须拼接可直接获得转录本的全长序列, 可以弥补非模式生物参考基因信息不完整的劣势, 亦可实现有参考基因组物种可变剪切及融合基因等结构变异的深入研究。目前三代全长测序技术已成功地用于一些海洋生物全长转录组的获取, 例如三疣梭子蟹(Portunus trituberculatus, Zhang et al, 2022)、凡纳滨对虾(Litopenaeus vannamei, Zhang et al, 2019)、脊尾白虾(Exopalaemon carinicauda, Shi et al, 2020)等。但是全长转录组测序技术在动物对深海化能极端环境的适应性研究方面的应用相对较少, 仅见Wang等(2022)对冲绳热液的光足深海阿尔文虾(Shinkaicaris leurokolos)和长角阿尔文虾(Alvinocaris longirostris)的全长转录组测定和比较, 通过分析其适应不同热液微环境的关键基因的表达及其进化特征揭示出阿尔文虾科物种适应热液不同生态位的分子基础。

十足目铠甲虾类是深海多样性最高的甲壳动物类群之一。在深海化能生态环境中, 铠甲虾是底栖生态群落中的常见种或优势种。劳盆拟刺铠虾(Munidopsis lauensis)隶属于十足目(Decapoda)、异尾下目(Anomura)、铠甲虾总科(Galatheoidea)、拟刺铠虾科(Munidopsidae)、拟刺铠虾属(Munidopsis), 是深海化能生境铠甲虾类代表性物种, 在西太平洋劳海盆(Lau Basin)、马努斯海盆(Manus Basin)、北斐济海盆(North Fiji Basin)、南海东北部以及印度洋等海域的热液、冷泉区均有分布, 分布水深在1 120~2 000 m之间(Baba, 2005; Martin et al, 2005; Cubelio et al, 2007; Hwang et al, 2022)。目前对劳盆拟刺铠虾分子生物学方面的研究仅局限于线粒体基因组的测定(Sun et al, 2019), 以及微卫星分子标记的筛选(Boyle et al, 2013)。针对劳盆拟刺铠虾深海适应性的相关研究尚未开展, 且该物种基因序列信息相对匮乏, 限制了对其适应性分子机制的深入研究。本研究利用PacBio平台的SMRT测序技术开展了劳盆拟刺铠虾全长转录组测序, 对获得的转录本进行注释分析, 预测其编码序列、转录因子、长链非编码RNA等, 挖掘可能与其适应深海化能极端环境相关的关键基因和通路。研究旨在丰富劳盆拟刺铠虾的基因资源, 为进一步探究深海化能生态系统中大型底栖甲壳动物的进化和环境适应分子机制提供数据支持。

1 材料与方法 1.1 样品采集

用于本研究全长转录组分析的劳盆拟刺铠虾样品于2021年6月由“科学号”(中国科学院海洋研究所)科考船上搭载的遥控潜水器Quasar MkII从位于中国南海的台西南冷泉(119°17.01′E, 22°06.91′N, ~1 118 m)采集。样品收集到船上后, 通过形态学初步鉴定立即进行解剖。将肝胰腺、鳃、肌肉和肠组织分离并在液氮中速冻后保存在–80 ℃条件下, 以便用于后续RNA提取。保留部分肌肉组织储存于浓度为70%的乙醇中, 以备于DNA提取与COI序列的扩增, 从而进行分子鉴定。

1.2 RNA提取

参照使用说明书, 使用Trizol试剂盒(Invitrogen, Carlsbad, CA, USA)提取各组织的总RNA。使用Nanodrop2000微分光光度计(Thermo Fisher)测定总RNA的纯度和浓度, 使用Agilent 2100生物分析仪(Agilent Technologies, Palo Alto, CA, USA)检测总RNA的完整度以评估RNA的质量。不同组织检测合格后的总RNA等量混合, 用于下一步全长转录组cDNA文库的构建。

1.3 文库构建、测序与数据处理

首先用Oligo(dT)磁珠富集mRNA, 再用NEBNext® Single Cell/Low Input cDNA Synthesis & Amplification Module将mRNA反转录成cDNA。然后进行PCR周期优化, 用于确定下游大规模PCR反应的最佳扩增周期数, 之后进行大规模PCR富集合成cDNA。cDNA经过DNA损伤修复, 末端修复, 并与adapter连接。将SMRTbell模板与测序引物退火并与聚合酶结合形成完整的SMRTbell文库之后, 在PacBio Sequel Ⅱ平台(美国太平洋生物科学公司)上进行测序。

使用SMRT Link-v8.0.0 (Gordon et al, 2015)处理下机的PacBio序列数据。首先, 从subreads.bam文件中提取高质量的环形一致性序列(circular consensus sequence, CCS)。通过去除序列中包含的5'端引物、3'端引物和多聚体结构, 获得全长非嵌合体(full length non-concatemer, FLNC)序列。然后使用Minimap2 (Li, 2018)对同一转录本的FLNC序列进行聚类, 获得到一致性序列。再使用Quiver算法对序列进行校正, 获得高质量一致转录本。使用CD-HIT-v4.6.7 (Cluster Database at High Identity with Tolerance) (Fu et al, 2012)对获得的一致性序列去冗余, 得到劳盆拟刺铠虾的全长转录本。最后, 基于BUSCO (Benchmarking Universal Single-Copy Orthologs)的节肢动物单拷贝直系同源基因标准集数据库, 使用BUSCO-v3.0.2 (Simão et al, 2015)对劳盆拟刺铠虾的全长转录本集进行完整性评估。

1.4 转录本的功能注释

为了获得劳盆拟刺铠虾全长转录本的注释信息, 将isoform序列与NCBI非冗余蛋白质(Nr)数据库(http://www.ncbi.nlm.nih.gov)、Swiss-Prot蛋白质数据库(http://www.expasy.ch/sprot)、KEGG数据库(http://www.genome.jp/kegg)和KOG数据库(http://www.ncbi.nlm.nih.gov/COG)进行比对分析, 阈值设置为1E-5。根据Nr数据库的注释信息, 使用Blast2GO (Conesa et al, 2005)软件对全长转录本进行GO (Gene Ontology)注释与功能分类。

1.5 转录本的结构分析

使用上述数据库blast比对和ANGEL (Shimizu et al, 2006)软件预测两种方法获得全长转录本蛋白编码区(CDS)的核酸序列和氨基酸序列。通过hmmscan将Isoform的编码序列与动物TFdb数据库(http://www.bioguo.org/AnimalTFDB/)进行比对, 预测转录因子(transcription factor, TF)家族。使用CNCI-v2 (Sun et al, 2013)和CPC (Kong et al, 2007)软件预测编码能力, 并将两个软件判定为非编码序列的结果合集作为最终长链非编码RNA (LncRNA)的预测结果。利用Perl脚本MISA (http://pgrc.ipk-gatersleben.de/misa/misa.html)对劳盆拟刺铠虾的全长转录组进行简单重复序列检测。

1.6 候选基因序列分析

使用Genedoc软件(Nicholas et al, 1997)对劳盆拟刺铠虾全长转录组中注释到的谷胱甘肽S-转移酶(glutathione S-transferase, GST)候选序列与其他节肢动物不同亚家族的GST氨基酸序列进行多重序列比对。基于邻接法(neighbor-joining, NJ), 使用MEGA软件(Kumar et al, 2016)进行系统发育树的构建, 并设置1 000次重复。使用Genedoc软件对劳盆拟刺铠虾GST theta亚家族的蛋白序列与包括柯氏潜铠虾、家蚕(Bombyx mori)、美州东部熊蜂(Bombus impatiens)、二化螟(Chilo suppressalis)在内的四种节肢动物的GST-theta蛋白序列进行比对。利用SMART (http://smart.embl-heidelberg.de/)预测劳盆拟刺铠虾GST-theta蛋白的结构域。

2 结果与分析 2.1 全长转录本测序数据

经COI序列扩增、测序与NCBI数据库进行同源比对后, 从形态和分子层面鉴定所采样品为劳盆拟刺铠虾。采用PacBio Sequel Ⅱ平台对劳盆拟刺铠虾肝胰腺、鳃、肌肉和肠组织的RNA进行全长转录组测序, 共获得55 402 300条原始序列, 平均长度为1 506 bp, N50长度为1 587 bp。经零模波导孔(ZMW)循环过后, 通过环型一致性序列分析, 共得到1 210 829个CCS, 平均长度为1 793 bp, 循环数为43。对从CCS得到的FLNC序列进行聚类后获得36 573条一致性序列, 经过校正后的高质量一致转录本为36 037条。去冗余后, 获得28 811条高质量非冗余isoform, 平均长度为2 086 bp, N50长度为2 275 bp (表 1)。BUSCO结果显示劳盆拟刺铠虾全长转录组中包括节肢动物核心单拷贝同源基因数的60.51% (图 1a), 表明本研究三代全长转录组测序获得的非冗余转录本相对完整, 可继续进行下游的数据处理。

表 1 基于PacBio测序的劳盆拟刺铠虾全长转录组测序结果 Tab. 1 Summary for the full-length transcriptome of M. lauensis using PacBio sequencing
参数 劳盆拟刺铠虾 Munidopsis lauensis
总碱基数 83.47 G
Subread 数目 55 402 300
Polymerases read N50 长度 1 587 bp
CCS read 数目 1 210 829
一致性序列数目 36 573
校正后的高质量转录本数目 36 037
非冗余高质量转录本数目 28 811
Isoform N50 长度 2 275 bp

图 1 劳盆拟刺铠虾转录本的功能注释 Fig. 1 Functional annotation of M. lauensis full-length transcripts 注: a. BUSCO 评估结果; b. KOG 注释结果统计; c. GO 功能分类
2.2 功能注释 2.2.1 基本注释

为了获得劳盆拟刺铠虾更全面的遗传信息, 将获得的全长转录本在Nr、SwissProt、KOG和KEGG数据库中进行了比对分析, 共有20 616 (71.6%)条isoform得到注释。Nr数据库搜索比对显示, 共注释到348个物种, 其中凡纳滨对虾(13 363; 65.01%)、克氏原螯虾(Procambarus clarkii, 881; 4.29%)和中华绒螯蟹(Eriocheir sinensis, 458; 2.23%)是比对到同源序列数目前三的物种。

2.2.2 KOG功能分类

共有17 957条isoform注释到25个KOG分类中。其中, 占比前三位的功能分类首先是一般功能预测(general function prediction only), 具有2 990条(16.65%)注释信息, 其次是翻译后修饰、蛋白质周转、伴侣(posttranslational modification, protein turnover, chaperones)和信号转导机制(signal transduction mechanisms), 分别有2 095 (11.67%)和1 690 (9.4%)条基因注释信息(图 1b)。

2.2.3 GO功能注释

将测序获得的基因与GO数据库进行匹配, 共有17 339条isoform被划分到三个大类。其中基因划分最多的类别是分子功能(molecular function, 6 669个), 其次是生物过程(biological process, 5 955个)和细胞成分(cellular component, 4 715个)。分子功能类别中可以细分为18个亚类别, 其中占比最多的是结合活性(binding), 占比42.04%; 其次是催化活性(catalytic activity), 占比34.78%。生物过程类别中注释到26个亚类别, 占比最多的是细胞过程(cellular process), 占比12.40%; 其次是代谢过程(metabolic process), 占比11.30%。注释到细胞成分类别中的基因被划分为23个亚类别, 其中比例最大的为细胞(cell)和细胞组分(cell part), 占比分别为16.98%和16.95%; 其次为细胞器(organelle), 占比为14.14% (图 1c)。

2.2.4 KEGG功能注释

共有20 117条isoform被注释到144条已知的KEGG通路中, isoform数量分布前20条通路见表 2。其中isoform数量最多的通路是代谢途径(metabolic pathways), 包括3 210条isoform, 其次是溶酶体(lysosome)与RNA转运(RNA transport), 分别有551条和439条isoform。

表 2 前20个KEGG代谢途径 Tab. 2 Top 20 KEGG pathways
通路名称 通路编号
代谢途径 ko01100
溶酶体 ko04142
RNA转运 ko03013
碳代谢 ko01200
内质网中的蛋白质加工 ko04141
氧化磷酸化 ko00190
吞噬体 ko04145
剪接体 ko03040
内吞作用 ko04144
氨基酸的生物合成 ko01230
核糖体 ko03010
花生四烯酸代谢 ko00590
糖酵解/糖异生 ko00010
长寿信号通路-蠕虫 ko04212
自噬-动物 ko04140
谷胱甘肽代谢 ko00480
真核生物的核糖体发生 ko03008
脂肪酸代谢 ko01212
泛素调节的蛋白质降解 ko04120
过氧化物酶体 ko04146
2.3 编码序列预测与全长转录组结构分析

通过对劳盆拟刺铠虾全长转录本进行CDS预测, 共得到21 848个CDS序列。CDS序列长度分布在102~7 719 bp, 平均值为1 049 bp; 5'UTR的长度分布在3~4 831 bp, 平均值为266 bp; 3'UTR的长度分布在3~7 189 bp, 平均值为970 bp。

2.3.1 转录因子

转录因子是在转录过程中识别和结合特定核苷酸序列的蛋白质, 通过与其他相关蛋白的相互作用介导附近基因的表达, 发挥重要的调控作用(Babu, 2010)。在劳盆拟刺铠虾全长转录组中共发现了41个转录因子家族, 包含537个转录因子。其中最多的转录因子家族是zf-C2H2家族, 其次是转录因子ZBTB和CSD家族(图 2a)。

图 2 劳盆拟刺铠虾全长转录组结构分析 Fig. 2 Analysis of full-length transcriptome structure of M. lauensis 注: a. 转录因子预测结果; b. SSRs 分布图
2.3.2 长链非编码RNA

LncRNA被定义为长度超过200个核苷酸且未被翻译成蛋白质的RNA分子, 它们从根本上参与一些生物过程, 如转录、翻译、蛋白质定位、细胞结构完整性、重编码以及其他细胞活动(Ma et al, 2013)。基于CPC和CNCI两种预测方法, 在劳盆拟刺铠虾全长转录组中共发现了6 430个LncRNA。

2.3.3 简单重复序列

简单重复序列又称为微卫星DNA、短串联重复序列, 通常指以1~6个碱基为单位组成的短串联重复序列。本研究共检测到6 517条包含SSR位点的isoform, 其中1 765条包含超过1个SSR位点, 以复合形式存在的SSRs数目为2 044个。进一步根据短串联重复单元的类型对所检测到的10 060个SSRs进行分类(图 2b)。依据重复核苷酸数目, SSRs可分为5种类型, 核苷酸重复数目为2~6个, 其中劳盆拟刺铠虾全长转录组中数量最多的是三核苷酸重复序列(4 233, 42.08%), 其次是二核苷酸重复(4 199, 41.74%)。SSRs重复次数在4~99次之间, 主要集中在4~7次, 且SSRs数量随着重复次数的增加而降低。

2.4 候选基因及通路分析

对isoform数量分布最多的前20个KEGG通路进行分析, 发现其中有两个通路可能与劳盆拟刺铠虾适应深海化能极端生境相关, 分别为过氧化物酶体(Peroxisome, ko: 04146)和谷胱甘肽代谢(Glutathione metabolism, ko: 00480)。在过氧化物酶体通路中, 共有180条isoform编码该通路的31个关键酶, 其中脱氢酶/还原酶SDR家族的数量最多, 为27条; 其次是D-天冬氨酸氧化酶(21条)和过氧化氢酶(19条); 甲羟戊酸激酶、α-甲酰基辅酶A消旋酶、(S)-2-羟基-酸性氧化酶、ATP结合盒亚家族D(ALD)成员2、2-羟基酰基-CoA裂解酶1、谷胱甘肽硫转移酶kappa1等数量最少, 分别有1条(图 3a)。共有213条isoform参与谷胱甘肽代谢通路中19个关键酶的编码, 其中GST的数量最多, 为59条, 其次是前列腺素H2 D异构酶谷胱甘肽转移酶(25条)和谷胱甘肽过氧化物酶(18条), 最少的是γ谷氨酰环化转移酶、蛋白质二硫还原酶(谷胱甘肽)以及谷胱甘肽特异性γ-谷氨酰环转移酶1, 数量仅为1条(图 3b)。

图 3 代谢通路图 Fig. 3 Maps of the metabolic pathway 注: a. KEGG过氧化物酶体代谢通路(map04146)示意图; b. KEGG谷胱甘肽代谢通路(map00480)示意图。红色边框标注劳盆拟刺铠虾全长转录组中注释到的酶

GSTs属于多功能酶, 可被分为15个亚家族, 包括Alpha、Beta、Delta、Epsilon、Kappa、Lambda、Mu、Omega、Pi、Phi、Rho、Sigma、Theta、Tau和Zeta (Huang et al, 2017)。劳盆拟刺铠虾全长转录组中共有20条isoform被分别注释到theta、delta、mu、kappa四个GST亚家族中。我们利用邻接法对劳盆拟刺铠虾全长转录组中注释到的GST候选序列与其他节肢动物的GST氨基酸序列进行了系统发育树的构建(图 4a)。结果显示劳盆拟刺铠虾全长转录组中被注释为不同亚家族的GSTs分别聚类在相应亚家族的分支上。多重序列比对结果显示劳盆拟刺铠虾的GST-theta序列与柯氏潜铠虾的序列相似度最高, 第2~74位氨基酸为N端结构域, 第108~201位氨基酸为C端结构域(图 4b)。

图 4 劳盆拟刺铠虾GST序列分析结果 Fig. 4 Result of sequence analysis on the glutathione S-transferase (GST) from M. lauensis 注: a. 利用劳盆拟刺铠虾(菱形表示)和节肢动物GST氨基酸序列构建的系统发育树; b. 劳盆拟刺铠虾与其他四种节肢动物GST-theta氨基酸序列比对。N 端结构域用蓝线表示; C 端结构域用红线表示
3 讨论

近年来, 以Pacbio平台为代表的第三代测序技术因其强大的读长优势改变了二代测序序列读长短和碎片化的局限性, 极大地提高了获得完整基因信息的能力(Eid et al, 2009)。全长转录组测序目前已成为从超大、复杂基因组物种中获取海量基因序列信息的重要工具, 并逐步应用于水生生物生长发育、胁迫响应、物种适应性进化以及重要基因挖掘等研究(Mehjabin et al, 2019; 王菁等, 2022; Cheng et al, 2022; Wang et al, 2022)。本研究通过Pacbio Sequel Ⅱ平台首次对深海化能生境铠甲虾类的全长转录组进行测序和分析。与之前深海化能极端环境中的其他甲壳动物二代转录组结果相比(Hui et al, 2018; Cheng et al, 2019), 虽然本研究所获得的转录本数量(28 811条)相对较少(143 516条, Hui et al, 2018; 108 237条, Cheng et al, 2019), 但在转录本完整性和质量(N50, 2 275 bp)方面却远高于二代(1 068 bp, Hui et al, 2018; 1 125 bp, Cheng et al, 2019) (图 5), 可为后续劳盆拟刺铠虾基因资源的挖掘提供基础, 并为探究该物种在不同条件或不同组织的转录表达差异提供了参考序列。功能注释是了解基因功能、研究基因的表达调控机制与预测蛋白功能的前提。劳盆拟刺铠虾各数据库的转录本注释率为71.56%, 显著高于基于二代测序的长角阿尔文虾(40.27%, Hui et al, 2018)和柯氏潜铠虾(36.78%, Cheng et al, 2019)。然而仍有8 195个转录本未被注释, 推测可能是铠甲虾类或该物种的特有基因, 公共的数据库中缺乏其相关的基因信息等原因造成。通过Nr数据库的注释发现, 劳盆拟刺铠虾与凡纳滨对虾比对到的同源信息最多, 所占比例为65.01%。仅注释到两个拟刺铠虾科物种, 分别为劳盆拟刺铠虾(116条)和威氏拟刺铠虾(Munidopsis verrilli, 2条), 表明拟刺铠虾科的已有基因资源相对匮乏。GO注释和KEGG注释结果显示劳盆拟刺铠虾的转录组功能主要集中于细胞过程与代谢过程。这些基因注释和分类信息将有助于对劳盆拟刺铠虾基因功能的理解与深入研究。

图 5 劳盆拟刺铠虾PacBio三代测序与柯氏潜铠虾和长角阿尔文虾Illumina二代测序转录本长度分布比较 Fig. 5 Comparison of distribution in transcript length between the PacBio sequencing of M. lauensis and the Illumina sequencing of S. crosnieri and A. longirostris

转录因子通过特异性结合调控区域的DNA序列来识别下游靶基因, 从而在调控各种生物过程的激活或者抑制中发挥重要作用(Charoensawan et al, 2010)。本研究预测了劳盆拟刺铠虾的转录因子, 其中zf-C2H2和ZBTB是最丰富的转录因子家族。C2H2型是锌指蛋白最常见的类型之一, 在真核生物的生长发育以及应激反应过程中起重要的调控作用(Liu et al, 2015)。ZBTB蛋白是锌指蛋白家族的亚家族之一, 在DNA损伤反应和细胞发育以及先天淋巴细胞的发育和分化等过程中起重要作用(Ren et al, 2019)。本结果与马氏珠母贝(Pinctada fucata martensii, Zhang et al, 2020)、口虾蛄(Oratosquilla oratoria, Cheng et al, 2022)、大菱鲆(Scophthalmus maximus, Fu et al, 2022)的全长转录组解析得到的转录因子分析结果相似。近年来研究发现长链非编码RNA可以通过与DNA相互作用来调节转录、表观遗传修饰、翻译和翻译后修饰等生物过程(Bridges et al, 2021)。截至目前, 尚未有关于劳盆拟刺铠虾的长链非编码RNA的报道。本研究鉴定出劳盆拟刺铠虾的6 430个LncRNA, 将为LncRNA具体功能的进一步研究提供基础。SSRs广泛存在于真核生物体内, 因其分布广、种类多及重复单元的重复次数在个体间呈现高度特异性, 已成为检测物种遗传多样性的有效工具。通过劳盆拟刺铠虾的全长转录组分析共预测到了10 060个SSRs, 结合Boyle等(2013)对马努斯盆地的劳盆拟刺铠虾开发的10个SSRs, 可为劳盆拟刺铠虾在不同海盆的遗传连通性和种群遗传结构研究提供候选分子标记。

环境压力可导致细胞内活性氧(ROS)水平的增高(Lushchak, 2011), 而过氧化物酶体是在ROS的代谢中起主要作用的细胞器(Mannaerts et al, 1993)。据报道, 过氧化物酶体在各种化合物包括诱导抗氧化酶的污染物刺激下在生物体内发生增殖(Orbea et al, 2002)。本研究发现过氧化物酶体通路位于isoform数量分布最多的前20条KEGG通路中, 且有部分转录本注释到参与过氧化物酶体生物发生的蛋白质peroxin编码基因上, 因此我们推测劳盆拟刺铠虾极有可能在深海极端环境中受到外界因子刺激后发生了过氧化物酶体增殖现象, 以应对细胞内高ROS水平。GSTs广泛分布于生物体内, 主要通过协助对抗内源或者外来有害物质从而实现细胞解毒功能(杨海灵等, 2006)。真核生物中参与谷胱甘肽代谢的GSTs主要属于MAPEG超家族, 在内源性脂类新陈代谢以及细胞解毒方面发挥重要的生物学作用(Jakobsson et al, 1996, 1999)。郑佩华等(2020)通过对凡纳滨对虾的GSTs进行基因克隆, 发现MGST3在对虾抗逆境胁迫以及抵御病原菌感染的调控机制中发挥重要作用。Cheng等(2019)在对深海热液和冷泉的共有甲壳动物柯氏潜铠虾进行转录组比较分析时发现编码GSTs的基因在热液个体中显著上调表达, 推测其有助于柯氏潜铠虾适应理化条件更极端的热液环境。在劳盆拟刺铠虾中, 我们发现了59个编码GSTs的isoform, 参与谷胱甘肽代谢过程。此外还检测到了同样具有抗氧化作用的谷胱甘肽过氧化物酶以及参与谷胱甘肽合成的γ谷氨酰转肽酶等。系统发育树与多重序列比对结果显示劳盆拟刺铠虾与柯氏潜铠虾GSTs序列相似度高, 可能与两物种的分类地位较近以及生存环境相似有关。不同家族的GSTs的亚基在序列上高度进化, 但都由N末端和C末端两个结构域组成。与谷胱甘肽结合的G位点存在于N末端结构域中, 与外源性和内源性亲电子底物结合的H位点存在于C末端结构域中, 这两者构成了GSTs的催化反应活性中心(Armstrong, 1991; Ranson et al, 2005)。进一步分析发现劳盆拟刺铠虾的N末端比C末端拥有更高的序列相似性, 这种现象在GTS家族中也是普遍存在的(Sinning et al, 1993)。这些基因可能在劳盆拟刺铠虾适应高浓度硫以及重金属等的深海化能环境过程中发挥重要作用。

4 结论

鉴于目前对深海极端环境下特殊生命过程的认识不足, 与浅海近缘物种的比较分析是解析深海生物进化与环境适应机制的关键。本研究利用SMRT三代测序技术对劳盆拟刺铠虾进行了全长转录组测序与分析, 不仅丰富了劳盆拟刺铠虾的遗传信息, 同时为深入揭示深海化能生态系统大型甲壳动物的适应性进化机制提供了研究基础与数据支撑。后续可基于本研究进一步开展拟刺铠虾科物种与其浅海近缘生物的比较转录组分析, 从而加深十足目甲壳动物对深海化能极端环境适应机制的理解。

参考文献
王菁, 李桂英, 林楷琪, 等, 2022. 基于马氏珠母贝(Pinctada fucata)血细胞全长转录组数据补体样组分的鉴定与分析[J]. 海洋与湖沼, 53(1): 176-186.
杨海灵, 聂力嘉, 朱圣庚, 等, 2006. 谷胱甘肽硫转移酶结构与功能研究进展[J]. 成都大学学报(自然科学版), 25(1): 19-24.
郑佩华, 汪蕾, 张秀霞, 等, 2020. 凡纳滨对虾微粒体谷胱甘肽硫转移酶3基因克隆及其功能分析[J]. 南方农业学报, 51(10): 2311-2320.
程娇, 沙忠利, 孙邵娥, 等, 2021. 深海化能生态系统大型生物多样性分布格局及其起源演化研究进展[J]. 海洋与湖沼, 52(2): 508-521.
ANDERSON R K, SCALAN R S, PARKER P L, et al, 1983. Seep oil and gas in Gulf of Mexico slope sediment[J]. Science, 222(4624): 619-621. DOI:10.1126/science.222.4624.619
ARMSTRONG R N, 1991. Glutathione S-transferases: reaction mechanism, structure, and function[J]. Chemical Research in Toxicology, 4(2): 131-140. DOI:10.1021/tx00020a001
BABA K, 2005. Deep-sea chirostylid and galatheid Crustaceans (Decapoda: Anomura) from the Indo-Pacific, with a list of species[J]. Galathea Reports, 20: 1-317.
BABU M M, 2010. Structure, evolution and dynamics of transcriptional regulatory networks[J]. Biochemical Society Transactions, 38(5): 1155-1178. DOI:10.1042/BST0381155
BOYLE E A, THALER A D, JACOBSON A, et al, 2013. Characterization of 10 polymorphic microsatellite loci in Munidopsis lauensis, a squat-lobster from the southwestern Pacific[J]. Conservation Genetics Resources, 5(3): 647-649. DOI:10.1007/s12686-013-9872-1
BRIDGES M C, DAULAGALA A C, KOURTIDIS A, 2021. LNCcation: lncRNA localization and function[J]. Journal of Cell Biology, 220(2): e202009045. DOI:10.1083/jcb.202009045
CHAROENSAWAN V, WILSON D, TEICHMANN S A, 2010. Genomic repertoires of DNA-binding transcription factors across the tree of life[J]. Nucleic Acids Research, 38(21): 7364-7377. DOI:10.1093/nar/gkq617
CHENG J, HUI M, SHA Z L, 2019. Transcriptomic analysis reveals insights into deep-sea adaptations of the dominant species, Shinkaia crosnieri (Crustacea: Decapoda: Anomura), inhabiting both hydrothermal vents and cold seeps[J]. BMC Genomics, 20(1): 388. DOI:10.1186/s12864-019-5753-7
CHENG J, ZHANG L W, HUI M, et al, 2022. Insights into adaptive divergence of Japanese mantis shrimp Oratosquilla oratoria inferred from comparative analysis of full-length transcriptomes[J]. Frontiers in Marine Science, 9: 975686. DOI:10.3389/fmars.2022.975686
CONESA A, GÖTZ S, GARCÍA-GÓMEZ J M, et al, 2005. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research[J]. Bioinformatics, 21(18): 3674-3676. DOI:10.1093/bioinformatics/bti610
CORLISS J B, DYMOND J, GORDON L I, et al, 1979. Submarine thermal springs on the Galápagos Rift[J]. Science, 203(4385): 1073-1083. DOI:10.1126/science.203.4385.1073
CUBELIO S S, TSUCHIDA S, HENDRICKX M E, et al, 2007. A new species of vent associated Munidopsis (Crustacea: Decapoda: Anomura: Galatheidae) from the Western Pacific, with notes on its genetic identification[J]. Zootaxa, 1435(1): 25-36. DOI:10.11646/zootaxa.1435.1.3
DUBILIER N, BERGIN C, LOTT C, 2008. Symbiotic diversity in marine animals: the art of harnessing chemosynthesis[J]. Nature Reviews Microbiology, 6(10): 725-740. DOI:10.1038/nrmicro1992
EID J, FEHR A, GRAY J, et al, 2009. Real-time DNA sequencing from single polymerase molecules[J]. Science, 323(5910): 133-138. DOI:10.1126/science.1162986
FU L M, NIU B F, ZHU Z W, et al, 2012. CD-HIT: accelerated for clustering the next-generation sequencing data[J]. Bioinformatics, 28(23): 3150-3152. DOI:10.1093/bioinformatics/bts565
FU Q, ZHANG P, ZHAO S C, et al, 2022. A novel full-length transcriptome resource from multiple immune-related tissues in turbot (Scophthalmus maximus) using Pacbio SMART sequencing[J]. Fish & Shellfish Immunology, 129: 106-113.
GORDON S P, TSENG E, SALAMOV A, et al, 2015. Widespread polycistronic transcripts in fungi revealed by single-molecule mRNA sequencing[J]. PLoS One, 10(7): e0132628. DOI:10.1371/journal.pone.0132628
HUANG X L, FAN D S, LIU L, et al, 2017. Identification and characterization of glutathione S-Transferase genes in the antennae of codling moth (lepidoptera: tortricidae)[J]. Annals of the Entomological Society of America, 110(4): 409-416. DOI:10.1093/aesa/sax041
HUI M, CHENG J, SHA Z L, 2018. Adaptation to the deep-sea hydrothermal vents and cold seeps: insights from the transcriptomes of Alvinocaris longirostris in both environments[J]. Deep Sea Research Part Ⅰ: Oceanographic Research Papers, 135: 23-33. DOI:10.1016/j.dsr.2018.03.014
HWANG H S, CHO B, CHO J, et al, 2022. New record of hydrothermal vent squat lobster (Munidopsis lauensis) provides evidence of a dispersal corridor between the Pacific and Indian Oceans[J]. Journal of Marine Science and Engineering, 10(3): 400. DOI:10.3390/jmse10030400
JAKOBSSON P J, MANCINI J A, FORD-HUTCHINSON A W, 1996. Identification and characterization of a novel human microsomal glutathione S-transferase with leukotriene C4 synthase activity and significant sequence identity to 5-lipoxygenase-activating protein and leukotriene C4 synthase[J]. Journal of Biological Chemistry, 271(36): 22203-22210. DOI:10.1074/jbc.271.36.22203
JAKOBSSON P J, THORÉN S, MORGENSTERN R, et al, 1999. Identification of human prostaglandin E synthase: a microsomal, glutathione-dependent, inducible enzyme, constituting a potential novel drug target[J]. Proceedings of the National Academy of Science of the United States of America, 96(13): 7220-7225. DOI:10.1073/pnas.96.13.7220
KONG L, ZHANG Y, YE Z Q, et al, 2007. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine[J]. Nucleic Acids Research, 35(S2): W345-W349.
KUMAR S, STECHER G, TAMURA K, 2016. MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets[J]. Molecular Biology and Evolution, 33(7): 1870-1874. DOI:10.1093/molbev/msw054
LI H, 2018. Minimap2: pairwise alignment for nucleotide sequences[J]. Bioinformatics, 34(18): 3094-3100. DOI:10.1093/bioinformatics/bty191
LIU Q, WANG Z, XU X, et al, 2015. Genome-wide analysis of C2H2 zinc-finger family transcription factors and their responses to abiotic stresses in poplar (Populus trichocarpa)[J]. PloS One, 10(8): e0134753. DOI:10.1371/journal.pone.0134753
LUSHCHAK V I, 2011. Environmentally induced oxidative stress in aquatic animals[J]. Aquatic Toxicology, 101(1): 13-30. DOI:10.1016/j.aquatox.2010.10.006
MA L N, BAJIC V B, ZHANG Z, 2013. On the classification of long non-coding RNAs[J]. RNA Biology, 10(6): 924-933. DOI:10.4161/rna.24604
MANNAERTS G P, VAN VELDHOVEN P P, 1993. Metabolic pathways in mammalian peroxisomes[J]. Biochimie, 75(3/4): 147-158.
MARTIN J W, HANEY T A, 2005. Decapod crustaceans from hydrothermal vents and cold seeps: a review through 2005[J]. Zoological Journal of the Linnean Society, 145(4): 445-522. DOI:10.1111/j.1096-3642.2005.00178.x
MEHJABIN R, XIONG L, HUANG R, et al, 2019. Full-length transcriptome sequencing and the discovery of new transcripts in the unfertilized eggs of Zebrafish (Danio rerio). G3 Genes| Genomes[J]. Genetics, 9(6): 1831-1838.
NICHOLAS K B, NICHOLAS H B JR, DEERFIELD Ⅱ D W, 1997. GeneDoc: analysis and visualization of genetic variation[J]. EMBnet News, 4(2): 14.
ORBEA A, ORTIZ-ZARRAGOITIA M, SOLÉ M, et al, 2002. Antioxidant enzymes and peroxisome proliferation in relation to contaminant body burdens of PAHs and PCBs in bivalve molluscs, crabs and fish from the Urdaibai and Plentzia estuaries (Bay of Biscay)[J]. Aquatic Toxicology, 58(1/2): 75-98.
RANSON H, HEMINGWAY J, 2005. Mosquuito glutathione transferases[J]. Methods in Enzymology, 401: 226-241.
REN R, HARDIKAR S, HORTON J R, et al, 2019. Structural basis of specific DNA binding by the transcription factor ZBTB24[J]. Nucleic Acids Research, 47(16): 8388-8398. DOI:10.1093/nar/gkz557
SHARON D, TILGNER H, GRUBERT F, et al, 2013. A single-molecule long-read survey of the human transcriptome[J]. Nature Biotechnology, 31(11): 1009-1014. DOI:10.1038/nbt.2705
SHI K P, LI J T, LV J J, et al, 2020. Full-length transcriptome sequences of ridgetail white prawn Exopalaemon carinicauda provide insight into gene expression dynamics during thermal stress[J]. Science of the Total Environment, 747: 141238. DOI:10.1016/j.scitotenv.2020.141238
SHIMIZU K, ADACHI J, MURAOKA Y, 2006. Angle: a sequencing errors resistant program for predicting protein coding regions in unfinished cDNA[J]. Journal of Bioinformatics and Computational Biology, 4(3): 649-664. DOI:10.1142/S0219720006002260
SIMÃO F A, WATERHOUSE R M, IOANNIDIS P, et al, 2015. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs[J]. Bioinformatics, 31(19): 3210-3212. DOI:10.1093/bioinformatics/btv351
SINNING I, KLEYWEGT G J, COWAN S W, et al, 1993. Structure determination and refinement of human alpha class glutathione transferase A1-1, and a comparison with the Mu and Pi class enzymes[J]. Journal of Molecular Biology, 232(1): 192-212. DOI:10.1006/jmbi.1993.1376
SUN L, LUO H T, BU D C, et al, 2013. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts[J]. Nucleic Acids Research, 41(17): e166. DOI:10.1093/nar/gkt646
SUN S, SHA Z L, WANG Y R, 2019. The complete mitochondrial genomes of two vent squat lobsters, Munidopsis Lauensis and M. Verrilli: novel gene arrangements and phylogenetic implications[J]. Ecology and Evolution, 9(22): 12390-12407. DOI:10.1002/ece3.5542
WANG A Y, SHA Z L, HUI M, 2022. Full-length transcriptome comparison provides novel insights into the Molecular Basis of adaptation to different ecological niches of the Deep-Sea hydrothermal vent in Alvinocaridid Shrimps[J]. Diversity, 14(5): 371. DOI:10.3390/d14050371
ZHANG X J, LI G Y, JIANG H Y, et al, 2019. Full-length transcriptome analysis of Litopenaeus vannamei reveals transcript variants involved in the innate immune system[J]. Fish & Shellfish Immunology, 87: 346-359.
ZHANG Y, NI M Q, BAI Y H, et al, 2022. Full-Length transcriptome analysis provides new insights into the diversity of immune-related genes in Portunus trituberculatus[J]. Frontiers in Immunology, 13: 843347. DOI:10.3389/fimmu.2022.843347
ZHANG Y J, SUN J, CHEN C, et al, 2017. Adaptation and evolution of deep-sea scale worms (Annelida: Polynoidae): insights from transcriptome comparison with a shallow-water species[J]. Scientific Reports, 7: 46205. DOI:10.1038/srep46205
ZHANG H, XU H Z, LIU H R, et al, 2020. PacBio single molecule long-read sequencing provides insight into the complexity and diversity of the Pinctada fucata martensii transcriptome[J]. BMC Genomics, 21(1): 481. DOI:10.1186/s12864-020-06894-3
ZHENG P, WANG M X, LI C L, et al, 2017. Insights into deep-sea adaptations and host–symbiont interactions: a comparative transcriptome study on Bathymodiolus mussels and their coastal relatives[J]. Molecular Ecology, 26(19): 5133-5148. DOI:10.1111/mec.14160