脊尾白虾响应环境胁迫的转录组学分析

王兴强， 卢雪旎， 沈 晔， 曹 梅， 韦寿永， 李庆国， 茆丹婷， 秦传新

(1.江苏海洋大学海洋科学与水产学院，江苏连云港 222005； 2.连云港旺岛旅游开发有限公司，江苏连云港 222000；3. 中国水产科学研究院南海水产研究所，广东广州 510300)

温度、盐度和溶解氧均是甲壳类生活环境中重要的环境因子。脊尾白虾()具有耐低温、低盐和低氧等特点，可以在温度2～35 ℃和盐度4～35的环境条件下生存，时刻关注这些环境因子的变化在白虾养殖中非常重要。低温胁迫可使凡纳滨对虾()血淋巴溶菌和抗菌能力显著降低，而低盐和低氧胁迫会降低其抗氧化能力和免疫酶活性。斑节对虾()体内酚氧化酶活性在低温胁迫下显著下降，而C型凝集素表达量在低盐胁迫下显著上升。低氧胁迫时日本沼虾()需要消耗大量糖原和磷酸精氨酸来维持能量需要，而高温时其幼体也会产生明显的应激反应。环境胁迫后的物种均可通过高通量转录组测序，经生物信息学分析来发掘其响应环境变化的调控基因。运用转录组测序技术，孙政等发现，脊尾白虾体内各种免疫相关基因表达上调以响应各种病原刺激。梅园通过尼罗罗非鱼()转录组测序，在温度胁迫产生的差异表达基因中发现一种共有结构域。本研究通过Illumina测序得到脊尾白虾低温、低盐和低氧胁迫下的转录组数据，统一组装，然后对不同胁迫条件下的差异表达基因进行注释和富集分析，以期阐明脊尾白虾在环境胁迫下的生理调控机制。

1 材料与方法

1.1 低盐试验

试验所用脊尾白虾购自连云港南极路水产品市场，在实验室暂养7 d以适应环境，暂养期间水体温度24～26 ℃、盐度31、pH值8.25以及水中溶解氧含量为7.6～8.1 mg/L。在低盐试验开始之际，随机挑选出120尾脊尾白虾[平均湿体质量(4.12±0.56) g]用于试验，设低盐胁迫(盐度0.2)、自然海水(盐度31)2个处理，每个处理3个重复，共6个水族箱(水体约30 L)，每个水族箱分别投放脊尾白虾20尾，然后对水体进行人工增氧。在进行低盐胁迫之后，水族箱中的白虾约在5.5 h之后虾体颜色变化大致可描述为从尾部逐渐变白，在6 h之后虾体整体呈现侧卧状，在这种状态下根据经验预测被低盐胁迫的脊尾白虾已处于昏厥状态。因此为了保证脊尾白虾的取样效果，对于低盐胁迫组的脊尾白虾，当其整体处于侧卧状态下5 min左右时开始进行取样操作，分别在低盐胁迫组与自然海水组的每个水族箱中随机选取脊尾白虾5尾，并使用液氮将整虾快速冷冻，转录组测序备用。

1.2 低温试验

脊尾白虾来源和暂养同“1.1”节，试验开始时随机挑选180尾脊尾白虾，平均湿体质量(1.27±0.15) g，设低温胁迫(0 ℃)、自然水温(25 ℃)2个温度梯度处理，2个处理各设置3个重复，总计水族箱6个(水体30 L左右)，其中每个水族箱分别放置脊尾白虾30尾，人工增氧，低温骤变处理4 h。对照组、低温胁迫组同时整虾取样，速冻于液氮中，转录组测序备用。

1.3 低氧试验

脊尾白虾来源和暂养同“1.1”节，试验开始时随机挑选60尾脊尾白虾，平均湿体质量(3.15±0.26) g，设低氧胁迫、自然溶解氧[溶解氧(7.59±0.66) mg/L]2个溶解氧梯度处理，2个处理各设置3个重复，总计水族箱6个(水体30 L左右)，每个水族箱分别放置脊尾白虾10尾。用2 mm厚液体石蜡封闭低氧胁迫组水族箱，胁迫约4 h后白虾处于昏厥状态，溶解氧降低到(1.13±0.26) mg/L，整虾取样，速冻于液氮中，对照组同样取样保存，转录组测序备用。

1.4 RNA提取与文库构建

用液氮将冻存的样品进行研磨，采取常规Trizol试剂法来提取总RNA，同时对它的纯度、浓度和完整性进行检测。样品达到要求后，富集mRNA并以其为模板合成双链cDNA链，最后进行PCR扩增，构建测序文库。

1.5 转录组测序与生物信息学分析

将样品数据送至上海元莘生物医药科技有限公司进行高通量测序，评估原始测序数据和文库构建质量，通过去除接头序列、剪切低质量读段、舍弃N率(测序中识别不出的碱基率)达到10%读数及长度小于25 bp片段的数据，最终得到高质量测序数据(Clean data)。通过Trinity对所有高质量数据进行从头组装，最终得到转录本和unigene。紧接着进行功能注释，运用BlastX分别与NR、STRING、Swiss-Prot和KEGG数据库进行比对，对基因的生物学功能进行进一步的探究。通过RSEM软件对基因的表达量进行评估，通过FPKM表示基因的表达水平，并运用edgeR软件进行基因的差异表达分析，主要富集到的功能和通路通过差异表达基因的GO富集分析和KEGG富集分析得到。

2 结果与分析

2.1 测序数据质量分析

经过滤后得到高质量测序数据，低盐、低温和溶解氧处理组和对照组均产出>35M条高质量序列，碱基数均>4 Gb，GC含量均>43%。测序质量的好坏是用质量值(Q，quality value)来评估的，Q30表示碱基测序出错的概率为0.001，值越高出错率越低。本次测序样品的Q30碱基百分比均大于92%(表1)。

表1 不同处理数据评估统计

使用Trinity对高质量数据进行序列组装，将Contigs合并分组生成Graph后得到转录本，最终组装得到转录本436 023条和unigene 305 682条，转录本和unigene的N50分别为814 bp和552 bp，其中序列长度为200～400 bp的unigene最多，占总数的67.79%，长度在1 kb以上的占8.11%(表2)。

表2 组装结果统计

2.2 unigene功能注释

unigene序列是通过BLAST在NR、Swiss-Prot、GO、COG、KEGG各大数据库中进行比对的，大部分unigene(56 898条)注释到NR数据库，31 695条注释到KEGG数据库， 28919 条注释到GO数据库。

COG功能分类结果见图1，主要注释到翻译、核糖体结构与生物生成(translation，ribosomal structure and biogenesis)，翻译后修饰、蛋白质周转与伴侣蛋白(posttranslational modification，protein turnover，chaperones)，仅通用功能预测(general function prediction only)和信号转导机制(signal transduction mechanisms)等功能中。

2.3 差异表达基因分析

差异表达基因中FDR<0.05并且log|FC|≥1的为差异显著。低温胁迫产生的差异表达基因中，具有显著性表达差异的有264条，包括215条上调基因和49条下调基因，其中159条得到注释。低盐胁迫产生差异表达显著的基因仅有1条，为显著下调的α-葡萄糖苷酶(α-glucosidase)，低氧胁迫没有产生差异表达显著的基因。

2.3.1 GO注释 低温胁迫产生的264条显著差异表达基因中，共有66条受到功能分类，GO注释结果见图2，其中，结构分子活性中富集到52个上调基因。

2.3.2 GO富集分析 使用Goatools软件对低温胁迫产生的差异表达基因进行GO功能显著性富集分析，为控制计算的假阳性率，对值进行校正，当校正后的<0.05时，GO功能显著富集。<0.001的有表皮结构组成(structural constituent of cuticle)、脂质转运活性(lipid transporter activity)和结构分子活性(structural molecule activity)等分子功能，其中表皮结构组成的显著性较高，所富集到的差异表达基因也较多，为51条。<0.01的有肽基二肽酶活性(peptidyl dipeptidase activity)和有机物运输(organic substance transport)等，<0.05的有单生物体运输(single organism transport)、单有机体定位(single organism localization)和底物特异性转运体活性(substrate specific transporter activity)等。同时对得到的功能节点标出编号和功能，以值进行显著性排序，最显著的10个节点信息见表3。

表3 GO富集(低温胁迫)

2.3.3 KEGG注释 使用KOBAS进行KEGG通路富集分析，经过校正的值以0.05为阈值，定义差异表达基因中显著富集的KEGG通路。低温胁迫产生的差异表达基因显著富集到肾素-血管紧张素系统(renin-angiotensin system，RAS)、美国锥虫病(american trypanosomiasis)、阿米巴原虫病(amoebiasis)、肾素分泌(renin secretion)、肥厚性心肌病(hypertrophic cardiomyopathy)、氨基糖和核苷酸糖代谢(amino and nucleotide sugar metabolism)、霍乱弧菌感染(vibrio cholerae)和补体与凝血级联(complement cascades with clotting)等通路中(表4)。

表4 KEGG富集部分结果(低温胁迫)

差异表达的基因在该通路的比例与所有注释基因在该通路的比例的比值越大，校正后的值越小，通路越具有参考价值。由肾素-血管紧张素系统信号通路(图3)可以看出，绿色边框有3个，表示3个下调基因。由表5展示信息发现，3个下调基因均为血管紧张素转换酶(angiotensin converting enzyme，简称ACE)基因，差异倍数对数值logFC分别为 -8.80、-6.65和-7.81。

表5 肾素-血管紧张素系统信号通路包含序列

3 讨论与结论

王余菊获取日本蟳()鳃组织转录组69.22 Gb高质量测序读数，Q30达93.46%以上，组装得到90 973条转录本和52 972条unigene，N50值为2 450、1 775 bp，注释率37.84%。董丽君等对凡纳滨对虾进行高通量测序，获得50 921条unigene。N50为1 589 bp，注释率37.28%。李喜莲等进行红螯螯虾()转录组测序，获得147 915 744条高质量读数，Q30达94.76%，组装得到67 369条unigene，N50长度为 1 376 bp，注释率30.83%。本次测序各样品均产出>35×10条序列数，碱基数均>4 Gb，GC含量均>43%，Q30碱基百分比均>92%，因而测序碱基识别出错率很低，测序质量较好。组装得到转录本436 023条和unigene 305 682 条，转录本和unigene的N50分别是814、552 bp。N50较短可能是本次测序数据量庞大，且200～400 bp长度范围占比较多所导致。56 898条unigene注释到NR数据库，31 695条注释到KEGG数据库，28 919条注释到GO数据库，注释率低的原因可能是在数据库中近源物种基因注释信息较少。

-葡萄糖苷酶又称为葡萄糖基转移酶(EC.3.2.1.20)，是糖苷水解酶家族的一员，水解碳水化合物中的-糖苷键生成葡萄糖。-葡萄糖苷酶具有双重作用，即水解和转糖苷，可以对-葡萄糖苷酶的活性进行抑制，从而降低生物体内葡萄糖的生成速率，同时可以对生物体内的脂肪水平进行调节。低盐胁迫仅产生1个显著下调的基因，为α-葡萄糖苷酶。低盐胁迫下α-葡萄糖苷酶表达量显著降低，可能导致脊尾白虾机体难以维持正常的能量需求。

肾素-血管紧张素系统(RAS)是重要的血压和水电解质平衡调节系统，而血管紧张素转移酶(ACE)是RAS中的一种重要的酶，又称为肽酰二肽水解酶，ACE分为ACE 1(原称ACE)和ACE 2，都是膜相关的锌肽家族成员，可以使肽链的C-末端二肽残基水解。血管紧张素Ⅱ受体1型(AT 1R)响应循环中的肽类激素血管紧张素Ⅱ(Ang Ⅱ)的系统性变化，ACE负责将Ang Ⅰ转换为AngⅡ，并通过调节血压和电解质稳态来使缓激肽失活。Ang Ⅱ在RAS中起到关键作用，能增强血管平滑肌张力，升高血压，还会造成组织损伤。AT 1R和ACE的特异性抑制剂对心血管系统有重要影响，血管紧张素转换酶抑制剂(ACEⅠ)具有使血压下降的作用。ACE 2水解AngⅡ生成血管紧张素 1-7[Ang(1-7)]，有舒张血管和抗炎的作用，抵消ACE、AngⅡ和AT 1R对脑神经元的损伤作用。低温胁迫产生的差异表达基因中，具有显著性表达差异的有264条，包括215条上调基因和49条下调基因，159条得到注释，其中66条受到GO功能注释，结构分子活性中富集到52个上调基因，推测与低温胁迫相关。GO富集分析发现，低温胁迫下最显著的功能节点是脂质转运活性和肽基二肽酶活性，而肽基二肽酶活性中包含的基因为ACE。同样，KEGG通路富集分析显示，差异表达基因显著富集到RAS中，且该通路富集显著性最高，具有参考价值。由信号通路图可以看出，该通路共包含3个表达下调的基因，均为ACE，推测低温通过抑制血管紧张素转换酶的表达，降低了血管紧张素Ⅱ的生成，从而降低血管平滑肌张力，抑制血压的升高和组织损伤；同时动员脂肪和蛋白质等能源物质，以响应低温胁迫。