陈 洁 董 梦 崔 悦 房翠莲 陈旭堂 刘子明

(1. 丽水学院生态学院, 丽水 323000; 2. 丽水市生态经济研究所, 丽水 323000)

拟穴青蟹(Scylla paramamosain)属于甲壳纲、梭子蟹科、青蟹属, 主要分布于我国东南沿海海域,是青蟹属的优势种[1—3]。“螯足腕节外缘中部发达刺的有无”是鉴定青蟹属不同种的主要分类特征之一。对福建漳州拟穴青蟹进行形态学分析, 发现拟穴青蟹大多数个体(80%左右)一对螯足腕节外缘中部都不具有发达的刺, 少数个体具有发达的刺[4]。于是, 将“一对螯足腕节外缘中部都不具有发达的刺”和“至少一只螯足腕节外缘中部具有发达的刺”的拟穴青蟹分别定义为Sp1和Sp2群体。通过研究上述两个形态群体拟穴青蟹鳃、肝胰腺和肌肉线粒体呼吸速率的季节变化, 发现Sp1和Sp2群体线粒体呼吸代谢的季节变化存在差异, 其中冬季差异最为显著[5]。

至今, 有关温度对拟穴青蟹的生物学效应已有不少的研究报道。Kong等[6]的研究结果表明, 冬季青蟹鳃超氧化物歧化酶(SOD)、过氧化氢酶(CAT)和谷胱甘肽过氧化物酶(GPX)活性显著低于春、夏和秋季, 且低温驯化会降低拟穴青蟹肌肉SOD、CAT和GPX活性[7], 由此说明在低温条件下青蟹代谢水平降低, 活性氧生成减少进而引起抗氧化酶活性的下降。近年来, 热休克蛋白(HSPs)成为研究生物受温度胁迫的一大热点[8—10]。目前拟穴青蟹hsp90[11]、hsp70[12]、hsp60[13]和hsp10[14]等热休克蛋白基因已被克隆, 初步的研究结果表明拟穴青蟹hsp70、hsp60和hsp10表达对温度变化十分灵敏。

上述这些研究主要从个体水平考察温度对拟穴青蟹的影响, 在群体水平上探索温度胁迫对拟穴青蟹生物学效应的报道很少。本研究以Sp1和Sp2两种形态群体拟穴青蟹为研究对象, 通过转录组测序分析, 筛选低温胁迫下差异表达基因, 将初步揭示拟穴青蟹不同形态群体对低温胁迫的适应性差异。这项工作的开展将有利于拟穴青蟹耐寒优质品种的选育及其种质资源的保护。

1 材料与方法

1.1 拟穴青蟹样品采集与处理

拟穴青蟹采样于广东省汕头海域(纬度约23°22′N)。采集到的拟穴青蟹活体运回实验室, 按照本课题组建立的方法[4]筛选出Sp1和Sp2个体。Sp1和Sp2个体的壳宽平均值分别为 (12.56±1.37)和(12.28±1.85) cm, 雌雄比例为1﹕1。选取大小相似、肢体完整、体表无菌斑的个体暂养于200 L的培养箱中, 内含150 L沙滤海水, 盐度22, 温度20℃, 暂养期间不投喂饲料, 暂养3d。拟穴青蟹Sp1和Sp2群体, 每个群体48只, 设置低温胁迫组(8.0±0.5)℃和对照组(20.0±0.5)℃, 每组设置3个平行组(每组每个群体8只), 分别放于6个培养箱中。6h后取肝胰腺, 液氮速冻, 随后存放于–80℃备用。

1.2 转录组测序

利用Trizol(生工, 中国)提取拟穴青蟹肝胰腺总RNA。1%琼脂糖凝胶电泳和 NanoDrop 2000(Thermo,美国)检测总RNA质量和浓度。将检测合格的总RNA送杭州联川生物技术股份有限公司测序。采用mRNASeq sample preparation kit(Illumina, 美国)构建cDNA文库, 利用IlluminaHiseq4000高通量平台(LC Sciences, 美国), 通过Paired-end 300 bp双末端测序模式对12个样本进行测序。

1.3 序列组装和基因注释

将原始序列中低质量、带接头和N比例大于10%的序列去除后, 利用Trinity 2.4.0软件进行无参拼接, 将拼接得到的最长序列作为基因的参考序列(Unigene)[15]。使用BUSCO(v5.1.2)软件对转录组进行完整性评估。随后用DIAMOND软件将Unigene分别与NCBI NR(http://www.ncbi.nlm.nih.gov/)、eggNOG(http://eggnogdb.embl.de/)、KEGG(http://www.genome.jp/kegg/)、Pfam(http://pfam.sanger.ac.uk/)、GO(http://www.geneontology.org/)和Swiss-prot(http://www.expasy.ch/sprot/)数据库进行比对和注释。

1.4 基因差异表达分析

基因表达量采用TPM(Transcripts Per Million)方法计算[16]。差异表达基因采用edgeR[17]软件包筛选, 筛选阈值为log2(fold change)>1和FDR<0.05。差异表达基因的GO和KEGG富集分析分别利用GOseq R和KOBAS软件包进行(FDR<0.05)[18,19]。

1.5 荧光定量PCR

选取上调和下调基因各5个。用TRIzol抽提样本总RNA, 利用AMV反转录酶(TaKaRa, 日本)反转录获得cDNA第一链。使用BioRad CFX96 Touch实时荧光定量PCR系统(BioRad, 美国)进行qPCR, 反应体系为: SYBR Premix ExTaq(2×)缓冲液(TaKaRa)12.5 μL, cDNA模板0.5 μL, 引物(10 μmol/L)各1 μL,以ddH2O补足至25 μ L。反应程序为: 94℃预变性5min; 95℃ 30s, 60℃ 30s, 72℃ 30s, 40个循环。qPCR检测结果使用2–ΔΔCt方法分析[20]。以18S rRNA作为内参基因, 所用引物见表 1, 引物的扩增效率根据荧光定量PCR的MIQE指南, 使用 cDNA样品的标准曲线确定, 引物的特异性根据熔解曲线分析确定。每组每个群体4个样品, 每个样品技术重复4次。

表1 引物序列Tab. 1 Primers used in this study

2 结果

2.1 测序和组装结果

低温胁迫组和对照组2个群体各取3个肝胰腺经过转录组测序、原始数据过滤后共计获得66.00 Gb有效数据(表 2)。有效数据已上传至NCBI Sequence Read Archive (SRA)数据库(登录号: PRJNA 743858)。Trinity组装去冗余后获得81853个Unigene。Unigene序列长度在201—17151 bp, 平均长度为420 bp,N50为1460 bp。BUSCO分析显示,C值(Complete)为917占比90.5%, 其中S值(Complete and single-copy)544,D值(Complete and duplicated)373, 分别占比53.7%和36.8%, 说明拼接的转录组完整性较好。

表2 不同形态群体拟穴青蟹肝胰腺转录组数据Tab. 2 Hepatopancreatic transcriptome data of Scylla paramamosain with different morphological populations

2.2 转录组序列功能注释

本次转录组测序共计获得81853个Unigene, 将Unigene分别与NR、Pfam、Swissprot、eggNOG、GO和KEEG数据库比对, 其中比对到NR数据库的基因数最多为18278, 注释率为22.33%; 比对到swissprot库的基因数最少为12556, 注释率为15.34%(表 3)。GO分析结果显示, 共有13951个Unigene得到注释。总共匹配到50个GO类, 归为3大类, 分别是生物进程(25类, 50.00%)、细胞组分(15类, 30.00%)和分子功能(10类, 20.00%), 其中富集在细胞质和细胞核等细胞组分的基因数最多。有12848条Unigene注释到KEGG数据库, 分属于40个代谢通路, 归为6大类。其中注释序列数量较多的代谢途径分别为运输和分解代谢; 翻译、信号转导; 折叠、分选和降解和糖代谢等。

表3 DIAMOND注释结果Tab. 3 DIAMOND annotation results

2.3 差异表达基因筛选

采用edgeR分析发现Sp1群体在低温胁迫下共有15773个差异表达基因, 其中上调145个, 下调15628个(图 1A); Sp2群体在低温胁迫下共有323个差异表达基因, 其中上调114个, 下调209个(图 1B)。

图1 低温胁迫下不同形态群体拟穴青蟹肝胰腺基因表达差异火山图Fig. 1 The ‘volcano plot’ picture of differentially expressed genes of Scylla paramamosain with different morphological populations under cold challenge

依据FDR值越小, 差异越显著的原则, 从两种形态拟穴青蟹在低温胁迫下差异表达基因中筛选到差异表达最显著的基因各10个(表 4)。

2.4 差异基因GO和KEGG通路富集分析

对所有差异表达基因进行GO富集和KEGG通路分析, 以获得其潜在的功能注释。Sp1群体在低温胁迫下差异表达基因在细胞质基质、RNA结合、核糖体结构组分和翻译等富集量较高(图 2A);Sp2群体在低温胁迫下差异表达基因在甾体羟化酶活性、DNA整合和碳酸氢盐跨膜转运体活性等富集较多(图 2B)。

图2 差异表达基因的GO富集分布散点图Fig. 2 GO enrichment results of differentially expressed genes

对差异表达基因进行KEGG富集分析, 结果显示, Sp1群体在低温胁迫下差异表达基因在核糖体、剪接体和内质网蛋白质加工等信号通路富集较多(图 3A), 而Sp2群体在低温胁迫下差异表达基因在赖氨酸降解、甘氨酸、丝氨酸和苏氨酸的代谢、精氨酸和脯氨酸的代谢等代谢通路富集较多(图 3B)。

图3 差异表达基因的KEEG富集分布散点图Fig. 3 KEGG enrichment results of differentially expressed genes

2.5 差异表达基因RT-qPCR验证

为了验证转录组测序结果的可靠性, 从表 4中挑选了5个上调基因和5个下调基因。如图 4所示,转录组表达分析结果和荧光定量PCR结果基本一致, 表明本研究转录组测序分析结果是可靠的。

图4 RT-qPCR验证差异表达基因Fig. 4 Validation of differentially expressed genes by RT-qPCR

表4 低温胁迫下两种不同形态群体拟穴青蟹差异表达最显著的前10个基因Tab. 4 The top 10 differentially expressed genes of Scylla paramamosain with different morphological populations under cold challenge

3 讨论

肝胰腺作为甲壳动物重要的组成部分, 具有消化吸收、解毒、清除氧自由基和免疫防护等功能[21,22]。前期研究已明确温度胁迫可影响拟穴青蟹肝胰腺的膜脂肪酸组成、线粒体呼吸速率、线粒体酶活性、抗氧化酶活性和热休克蛋白表达等[4,23,24]。因此本研究以拟穴青蟹肝胰腺为材料, 通过转录组测序分析比较Sp1和Sp2群体在低温胁迫下转录水平的差异。测序共获得81853条Unigene, 其中22.33%为已知基因, 这在一定程度上丰富了拟穴青蟹转录组数据库。

变温动物在应对低温时, 通常通过温度补偿机制来提高代谢速率, 以维持正常生命活动。但当环境温度剧烈变化时, 一些变温动物则会停止代谢进入休眠状态, 以度过恶劣的环境[25―27]。水温从20℃骤降至8℃, 在Sp2群体拟穴青蟹肝胰腺中检测到

323个Unigene显著差异表达, 其中上调114个, 下调209个。这些Unigene在甾体羟化酶活性、DNA整合和碳酸氢盐跨膜转运体活性等GO富集较多, 并进一步富集在赖氨酸降解、甘氨酸、丝氨酸和苏氨酸的代谢、精氨酸和脯氨酸的代谢等代谢通路。水温从24℃缓慢降至18℃(每天降温2℃), 凡纳滨对虾(Litopenaeus vannamei)肝胰腺的转录组测序结果同样表明低温胁迫下差异表达基因富集在氨基酸的生物合成途径以及精氨酸和脯氨酸的代谢中[28]。凡纳滨对虾3个家系的差异表达Unigene分别为1973、1677和687个[28], 多于经历水温骤降12℃的拟穴青蟹Sp2群体, 由此推测Sp2群体的拟穴青蟹对水温快速下降的代谢响应相对迟缓。

对于Sp1群体, 差异表达基因筛选的结果显示低温胁迫下共有15773个Unigene显著差异表达。广温性草鱼(Ctenopharyngodon idellus)从27℃水温缓慢降至12℃和4℃, 脑转录组测序分别检测到1273和1987个差异表达Unigene[29]。对低温敏感的热带鱼虎皮鱼(Puntius tetrazona)从27℃水温骤降至13℃, 脑、鳃和肝脏转录组测序分别检测到10272、7771和4632个差异表达Unigene[30]。由此可看出,面对水温骤降, 拟穴青蟹Sp1群体的代谢响应远比Sp2群体敏感。Sp1群体肝胰腺下调表达的Unigene在细胞质基质、RNA结合、核糖体结构组分和翻译等GO富集量较高, 并进一步富集在核糖体、剪接体和内质网蛋白质加工等信号通路, 这说明Sp1群体拟穴青蟹蛋白质合成加工过程大幅度减弱。

面对不利环境, Sp1群体拟穴青蟹能快速地降低生理活动, 减缓低温对机体的损伤, 而Sp2群体拟穴青蟹在应对低温时, 仅调整了物质代谢, 这可能使得Sp2群体在应对低温时机体更容易受到伤害,甚至导致死亡现象的发生。在早期的研究中, 我们调研了福建漳州海区的Sp1和Sp2群体, 其中Sp1群体的数量是Sp2群体数量的3倍多, 这可能与两个群体对温度变化存在适应差异有关。在本研究的基础上, 我们将进一步挖掘转录组数据同时结合蛋白质组学分析寻找其他低温胁迫适应性指标来进一步揭示Sp1和Sp2群体对温度的适应差异。