急性高温低氧胁迫下刺参体腔液细胞miRNA表达谱的构建和靶标基因的功能分析❋

尹英超，刘富祥，张玉钦，陈慕雁❋❋

(1.中国海洋大学水产学院，山东 青岛 266003；2.海水养殖教育部重点实验室(中国海洋大学)，山东 青岛 266003；3.烟台中集蓝海洋科技有限公司，山东 烟台 264000)

刺参(Apostichopusjaponicus)属棘皮动物门(Echinodermata)海参纲(Holothuroidea)辛那参目(Synallactida)刺参科(Stichopodidae)，广泛分布于中国、日本、韩国和俄罗斯等东北亚沿海[1]。刺参作为底栖生物以底泥、碎屑为食，在维持生态系统平衡中发挥着重要的作用[2]。此外，刺参具有很高的营养价值和药用价值，长期以来在亚洲市场上作为一种食品和民间药物备受欢迎，是中国重要的海水养殖经济品种之一，目前已经成为北方水产养殖中的支柱性产业[1,3]。

全球气候的变化导致海水温度升高和缺氧现象已成为一个世界性环境问题，严重影响到海洋生物的生存。“厄尔尼诺”是导致太平洋海水表面温度异常升高的重要因素之一，可以引发严重而广泛的天气变化，如高温和暴雨[4]。2013年夏末以来，中国山东、河北两省因“厄尔尼诺”现象引发的连续高温暴雨天气，导致海水温度高达30 ℃，溶解氧水平降至2 mg/L以下；由于刺参长期营底栖生活，活动缓慢，受环境因素影响非常显著，极端高温、低氧胁迫可导致刺参大量死亡，严重制约刺参养殖产业的可持续发展。为了给刺参应对逆境关键指标筛选和抗逆品种选育提供理论支持，刺参在急性高温、低氧胁迫下内在分子调控机制的科学研究亟待开展。

棘皮动物通常缺乏特异性免疫，体腔液是棘皮动物进行机体防御的重要免疫组织。刺参体腔液中悬浮着各种体腔液细胞，同时也含有凝集素、水解酶、酚氧化酶、超氧化物歧化酶类白介素和类补体样因子等体液免疫因子，两者共同作用负责刺参的天然免疫防御[5-6]。已有大量研究表明棘皮动物在响应外界温度、盐度、溶解氧和pH等环境胁迫因子时，体腔液及其细胞的相关免疫抗病能力会发生变化[7-10]。

miRNA是一类长度为21～24 nt的内源性非编码RNA，miRNA能够通过种子序列与特定的靶基因位点相结合，引起mRNA降解或抑制mRNA翻译，从而对靶基因的表达起到调节作用。目前，关于miRNA介导高温、低氧胁迫的调节机制在水生生物中进行了广泛研究。Zheng等在日本对虾(Marsupenaeusjaponicus)体内共发掘了15个miRNA参与高温胁迫应答，靶标基因的预测和功能分析表明这些miRNA参与了信号传导、转录、抗胁迫应激等细胞过程[11]。Zhou等研究发现虹鳟鱼(Oncorhynchusmykiss)在致死温度22.5 ℃胁迫48 h下，头、肾组织中共发掘出12个上调和5个下调的miRNA，6个差异表达的miRNA共预测到22个靶基因，注释结果表明在热应激下miRNA主要参与了虹鳟鱼的代谢和免疫反应[12]。Hadj-Moussa等发现miRNA调控的代谢抑制是保障洪堡鱿鱼(Dosidicusgigas)能够在缺氧水域生存的重要机制，在缺氧诱导下鱿鱼的大脑、外套膜肌肉和心脏中表达量变化的miRNA主要参与了关键细胞功能维持和细胞保护的过程[13]。Wang等在南美白对虾(Litopenaeusvannamei)应对缺氧的研究中，鉴定出许多与缺氧耐受相关的miRNA，其差异表达的miRNA可能通过调节应激反应基因(如HIF1a和p53)的表达来调节宿主对缺氧的反应，并影响与缺氧适应相关的关键细胞事件[14]。Sun等在黒鲈(Micropterussalmoides)缺氧胁迫的研究中筛选出的13个差异表达的miRNAs的靶基因在血管内皮生长因子(VEGF)信号通路、丝裂原活化蛋白激酶(MAPK)信号通路和磷脂酰肌醇信号系统中显著富集[15]。关于刺参miRNA研究起步较晚，最早于2012年Li等从刺参的体腔液细胞中鉴定到miRNAs[16]，之后的研究表明miRNA广泛参与了刺参补体系统、脂多糖以及Toll样受体的调控[17-22]。近年来，关于miRNA在刺参响应环境胁迫中的调控作用研究也相继开展。Li等构建了刺参在高温胁迫下的miRNA表达谱，发现大量差异表达的miRNA参与免疫过程相关的调节[23]；Huo等在刺参缺氧胁迫研究中发现miRNA可能在氧化还原、运输以及耐缺氧胁迫中发挥作用[24]。Tian等研究发现miR-10可能通过靶标TBC1D5基因参与低盐胁迫下刺参的能量调节，并在维持细胞稳态方面起到了关键作用[25]。Huo等以刺参呼吸树为研究对象，构建了在高温、低氧胁迫下刺参的miRNA表达谱，并在高温(26 ℃)、低氧(2 mg/L)和高温低氧复合胁迫条件下，分别鉴定出21、26和22个差异表达miRNA，GO和KEGG通路分析显示差异miRNA的潜在靶基因参与了生物合成、代谢、免疫、细胞生长和死亡、翻译及信号转导等过程[26]。

本研究以刺参高温、低氧胁迫下的体腔液细胞为研究对象(由于高温低氧协同组的刺参全部死亡，无法进行后续的研究)，采用Illumina测序技术，结合生物信息学分析，构建了刺参在高温、低氧胁迫下的miRNA表达谱，使用miRanda软件，以构建的三代全长转录组为参考文库，预测了差异表达miRNA的靶标基因，并对候选靶标基因进行GO和KEGG富集分析，发掘了刺参在高温、低氧胁迫下的差异miRNA参与调控的关键通路，探讨了miRNA在刺参体腔液细胞响应高温、低氧胁迫中的潜在调控机制，为筛选刺参抗逆关键指标提供理论支持。

1 材料与方法

1.1 材料

1.1.1 实验动物暂养 实验所用的成熟刺参于2019年4月采集于威海刺参养殖场，低温运到实验室。挑选120头健康无病、发育良好的刺参(体质量(120±20)g)分别放入12个大小、形状、材质一致的蓝色避光的养殖水箱中，每个水箱中放置10头刺参，为了消除环境因素的差异，在海边收集大小不一的碎石经消毒清洗后，放入水箱中供刺参附着，箱内水温(T)为(15±1)℃，盐度为30±1，pH为8.0±0.5，充气保持溶解氧(DO)为8 mg/L，在此环境中驯化2周。期间每日上午定时清除水箱底部的粪便和剩余的饲料，饲料与底泥按1∶3比例混合投喂。

1.1.2 实验处理 对照组的条件设置为T=15 ℃、DO=8 mg/L，高温处理组T=30 ℃，DO=8 mg/L，低氧处理组T=15 ℃，DO=2 mg/L，高温低氧处理组T=30 ℃，DO=2 mg/L。高温设置方法为使用功率1 000 W的加热棒加热来达到实验所需的水温，并通过控温器来保持水温，水温升高速率为2 ℃/h。低氧设置方法为向水体中充入氮气，控制氮气充入速率，使溶解氧下降速率保持在0.8 mg/L。达到实验所需溶氧条件后，使用塑料薄膜覆盖水箱的水面，在低氧胁迫期间使用便携式溶氧仪实时监测溶解氧水平。对照组和实验组分别设置3个平行，对照组和实验组在达到所需温度溶氧条件后48 h进行取样。

1.1.3 样品收集 每个平行组中取4头刺参，设3个平行组，每个处理组共取12头刺参。将每个刺参的体腔液分别取3份，每份体积为1 mL，在4 ℃，1 000 r/min的冷冻离心机中离心30 s,去掉上清液后将剩余细胞混合，加入1 mL的Trizol试剂，立即放入液氮中速冻后转入-80 ℃超低温冰箱中保存备用。由于刺参无法耐受48 h的高温低氧(T=30 ℃, DO=2 mg/L)的协同作用而全部死亡，无法取得该实验组样品进行后续分析。

1.2 RNA文库构建

采用Trizol法提取细胞总RNA，具体步骤见李兴可等[27]。使用Nanodrop检测RNA的纯度(OD260/280比值)；使用Qubit对RNA浓度进行精确定量；使用Agilent 2100精确检测RNA的完整性。使用NEBNext® Multiplex Small RNA Library Prep Set for Illumina®(NEB, USA)生成测序文库，文库中对照组为CO1、CO2、CO3，高温组为HS1、HS2、HS3，低氧组为HY1、HY2、HY3，建库流程如图1所示。

图1 Illumina样品检测、建库和测序的流程图

1.3 miRNA测序与数据分析

采用Illumina Hiseq 2500进行miRNA测序，去除低质量的reads(质量值小于等于20 的碱基占整个reads的30%以上)；去除 N(N表示无法确定碱基信息)的比例大于10%的 reads；去除接头污染的reads；去除3’接头序列；去除polyA/T/G/C的reads。选择一定长度的clean reads用于下游的分析，miRNA的长度集中在21～22 nt。将小RNA标签映射到RepeatMasker和Rfam数据库，除去源自蛋白质编码基因、重复序列、rRNA、tRNA、snRNA和snoRNA等。用bowtie将长度筛选后的sRNA定位到刺参基因组(NCBI BioProject ID: 29023486)序列上，分析small RNA在参考序列上的分布情况。

1.4 miRNA的鉴定和差异表达分析

使用miRBase20.0作为参考，鉴定已知的miRNA；结合miRNA前体发夹结构特征，采用软件miREvo和mirdeep2预测新的miRNA。将所有sRNA标签与成熟miRNA进行比对，检测可能有碱基编辑的miRNA，已知的miRNA使用miFam.dat(http://www.mirbase.org/ftp.shtml)查找家族。将新的miRNA前体提交到Rfam库(http://rfam.sanger.ac.uk/search/)，以寻找Rfam家族。根据下面算法，应用TPM(每百万转录本)估算miRNA表达水平。

规范化表达式=(映射的读数/总读数)×1 000 000。

使用DESeq R软件包(1.8.3)对两个实验组进行差异表达分析。使用Benjamini&Hochberg方法调整P值。默认情况下，P<0.05为显著差异表达。

1.5 靶标基因的预测和分析

利用miRanda分析软件，以前期构建的三代全长转录组数据作为参考文库，从中提取3’非编码区域，预测靶标基因。将筛选参数设置为自由能小于-83.68 kJ/mol和得分阈值大于140。

1.6 GO和KEGG富集分析

基于GOseq的Wallenius非中心超几何分布调整基因长度偏差，用于GO富集分析。Pathway显著性富集分析以KEGG Pathway为单位，应用超几何检验，该分析计算公式如下：

得到的P值经Bonferroni方法校正，以校正后的P值小于0.05定义为显著富集的pathway。

1.7 qRT-PCR验证

应用SYBR®Premix Ex TaqTMII Kit(Takara，Japan，Cat No.RR820)试剂盒，分别取刺参对照组、高温组、低氧组的样品各5个平行，每个平行样品设置3个技术重复，将成熟miRNA(21～23 nt)的完整序列用作miRNA特异性5’端引物，引物序列如表1所示。

表1 miRNA qRT-PCR的引物序列

qPCR的3’端引物使用试剂盒提供的mRQ 3’通用引物，以U6作为内参基因，其引物序列为：F:GGAACGATACAGAGAAGATTAGC, R: TGGAACG-CTTCACGAATTTGCG。变性的反应条件：95 ℃ 10 s，95 ℃ 5 s，60 ℃ 20 s，40个循环。融解曲线的反应条件：95 ℃ 60 s，55 ℃ 30 s，95 ℃ 30 s。使用2-ΔΔCt法计算 miRNA 的相对表达水平，使用单因素方差分析miRNA的组间差异显著性，以P<0.05为差异显著。

2 结果与分析

2.1 测序结果分析

建库之后，利用Illumina HiSeqTM2500测序平台将原始图像数据文件经过碱基识别转化为测序序列，得到的Raw Data以FASTQ格式储存，原始数据已上传至NCBI(BioProject ID: PRJNA723369, SubmissionID: SUB9504899)。检查测序错误率异常的碱基位置分布，发现9个文库测序结果中每个位置碱基的错误率都低于0.5%，详细结果见附表1。

2.2 sRNA的测序分析结果

为了保证信息分析的质量，对得到的Raw Data数据进行过滤，得到Lean Data，筛选出sRNA进行后续的分析，详细结果见附表2。

用bowtie将长度筛选后的sRNA比对到刺参参考基因组序列上，发现每组中90%以上sRNA能够比对上参考基因组，其中HY3组能比对到参考基因组的数量最少(91.96%)，HY1组能比对上参考基因组的数量最多(98.27%)，与参考基因组方向相同链的reads数占全部clean reads数的83%以上，方向相反链的reads数占全部clean reads数的10%以下，结果如表2所示。

表2 sRNA与参考基因组比对结果统计

2.3 已知miRNA与ncRNA分析

所有样本中共比对到33个已知的miRNA，比对上的sRNA的种类有3 242种。对Novel miRNA成熟体、前体、种类和总数进行统计，详细结果见附表3。总共匹配到的前体和成熟体的数目分别为149和159，详细结果见附表4。

2.4 miRNA家族分析

对检测到的已知miRNA和新预测到的miRNA进行家族分析，结果显示，已知miRNA和新miRNA总共属于31个不同的miRNA家族，其中let-7在物种中的分布最丰富。大多数miRNA只检测到一个前体，miR-10、miR-124和miR-25检测到3个前体。miRNA家族分析部分结果见附表5。

2.5 miRNA差异表达分析

在HS组中共发掘出20个差异表达的miRNA，其中11个上调，9个下调；在HY组中共发掘出10个差异表达的miRNA，其中6个上调、4个下调(见表3)。两组中共同差异表达的miRNA有6个：novel-23、novel-43、novel-45、novel-49、miR-200-3p和miR-71。为验证结果的准确性，分别从高温组和低氧组中随机选取了10个差异miRNA，其中高温组选取miR-133、miR-71、miR-10、miR-2012、miR-200-3p和Novel-32，低氧组选取miR-2006、miR-71、miR-200-3p和Novel-43，利用荧光定量PCR进行验证。结果表明，高通量测序结果和qRT-PCR结果趋势相同(见图2)。

图2 差异表达miRNAs的qRT-PCR验证结果

表3 高温组和低氧组差异表达miRNA的详细结果

2.6 靶基因预测和功能富集分析

我们以刺参的三代全长转录组中的191 141个转录本作为参考文库，预测差异表达miRNA的靶基因。结果表明，高温组中20个差异表达miRNA共预测到14 827个靶基因，低氧组中10个差异表达的miRNA共预测到8 577个靶基因。

我们对每组的差异表达miRNA的靶基因分别进行了GO和KEGG富集分析。富集分析结果显示，在高温胁迫组中，差异表达miRNA的靶向基因在细胞组分中聚类最明显的是细胞成分，在分子功能项中聚类最明显的是转运活性和跨膜转运体活性(Transmembrane transporter activity)，在生物学进程中跨膜运输(Transmembrane transport)聚类最为明显。在低氧胁迫组中，差异表达miRNA的靶向基因在细胞组分中聚类最显著的是细胞成分和膜部分(Membrane part)；在分子功能项中聚类最显著的是转运活性，在生物学进程中内肽酶活性(Endopeptidase activity)聚类最显著，详细结果如图3所示。

图3 差异表达miRNA靶基因的GO富集分析

KEGG富集分析结果显示，高温组中显著富集的通路有36个，其中富集到基因数目较多的通路主要集中在磷脂酰肌醇激酶-丝氨酸/苏氨酸激酶信号通路(PI3K-Akt signaling pathway)、癌症通路(Pathways in cancer)、内质网蛋白质加工(Protein processing in endoplasmic reticulum)和黏着斑(Focal adhesion)。在低氧组中显著富集的通路有37个，富集靶基因数目较多的通路主要集中在黏着斑(Focal adhesion)、癌症通路、磷脂酰肌醇激酶-丝氨酸/苏氨酸激酶信号通路和丝裂原活化蛋白激酶信号通路(MAPK signaling pathway-fly)，结果如图4所示。

((a)高温组；(b)低氧组。(a)Heat stress group；(b)Hypoxia group.)

3 讨论

随着高通量测序技术的快速发展，刺参不同组织中(肠道、肌肉、体壁、管足以及体腔液等)miRNA表达谱得到广泛的构建[22]，表达谱分析表明miRNA在刺参生物学过程中具有潜在的关键调控作用，但在刺参体腔液细胞中miRNA响应高温、低氧胁迫的调控机制的研究还少有报道。

体腔液作为刺参内外环境的渗透介质，外界环境的变化可能影响体腔液细胞的活性和功能，高温和低氧是目前影响刺参生存最重要的环境因子。本研究共构建了9个刺参体腔液细胞的miRNA文库，测序结果显示大部分成熟体miRNA长度分布在21～23 nt，峰值大小为22 nt，这与大部分海洋无脊椎动物的研究结果相似，略微的差异可能是由物种差异和组织特异性导致的。

本研究从刺参的体腔液细胞中共鉴定到33个已知的miRNA和149个Novel miRNA。多数miRNA在各物种中是高度保守的，发挥了广泛而重要的调控作用。其中高温组中miR-133、miR-10、miR-2012、miR-200-3p、miR-71和低氧组中miR-1、miR-2006、miR-200-3p、miR-71都有显著差异表达。有研究表明，miR-10是刺参体腔液细胞中最丰富的miRNA[16]，且在细胞凋亡[28]、蛋白合成[29]、胚胎发育和分化[30-31]、造血干细胞调动[32]和炎症[33]等生物过程中发挥着重要的调控作用，这表明miR-10在刺参响应高温胁迫中同样可能参与调控多个生理学过程。De Lencastreetal等研究发现miR-71与细胞的衰老和DNA损伤反应有关[34]，另有研究发现miR-71同样可以增强对热休克和氧化应激的抗性[35]。Huo等最近研究发现miR-71家族成员在高温胁迫中是显著上调的[26]，这和本文的结果一致，也暗示了本研究中miR-71在刺参响应高温、低氧胁迫下的潜在功能。miR-1已被证实与凋亡相关基因如热休克蛋白基因相关，并间接调节eNOs以及细胞凋亡[36]，以往研究也发现刺参呼吸树中miR-1在低氧胁迫下显著上调[24]，这和本文的研究结果一致。miR-133家族因为在抑制人类癌症发生中的关键作用被充分研究[37]，在心肌细胞中miR-133通过靶向Hsp70抑制细胞凋亡[38]，相关研究还直接证明miR-133通过靶向AjIRAK-1参与Toll样受体(TLR)级联调节，从而促进刺参体腔细胞吞噬功能[39]。本研究中miR-133在高温胁迫下显著上调，这暗示了miR-133可能参与了高温胁迫下体腔液细胞抗凋亡过程和免疫应激。miR-200家族是哺乳动物中最著名的miRNA家族之一。作为该家族的重要一员，miR-200-3p被报道参与细胞增殖、凋亡和肿瘤发生等生物学过程[40-41]。本课题组在刺参夏眠的研究中也发现了miR-200-3p能够通过靶向AjEHHADH和AjCA参与刺参夏眠期间氧化损伤和细胞周期调控[42-43]。在本研究中，miR-200-3p在高温和低氧胁迫下均显著上调，揭示了miR-200-3p可能在刺参响应高温、低氧胁迫下参与了广泛的生物学功能调节。目前关于miR-2012的功能研究还有待开展。

为进一步发掘miRNA潜在的调控机制，我们以三代转录组数据为参考文库，预测了差异表达miRNA的候选靶基因。单个miRNA可以靶向多个靶基因，每个靶基因也可以被多个miRNA调控，因此，miRNA通过介导复杂的网络结构发挥调控功能。靶基因GO分析结果显示，大多数靶基因与细胞组分、转运以及跨膜转运活性有关。靶基因KEGG通路富集分析发现多数基因在丝裂原活化蛋白激酶信号通路、磷脂酰肌醇激酶-丝氨酸/苏氨酸激酶信号通路、寡聚化核苷酸结合结构域样受体信号通路(NOD-like receptor signaling pathway)、黏着斑等信号通路中显著富集。MAPK级联协同传递各种细胞外信号，从而控制大量不同甚至相反的细胞过程，如增殖、分化、生存、发育、应激反应和凋亡[44]，有研究发现MAPK在调节凡纳滨对虾(Litopenaeusvannamei)、香港牡蛎(Crassostreahongkongensis)等多种水生生物的先天免疫过程中发挥着重要作用[45-46]。同样在刺参体腔液高温胁迫的研究中发现p44/42MAPK可能通过靶向p90RSK参与温度诱导的代谢抑制[47]。PKB/Akt通路同样能够介导多种细胞功能，包括血管生成、代谢、生长、增殖、存活、免疫、转录和细胞凋亡[48]，在水产动物的研究中发现Akt基因可能参与介导水产动物机体的免疫防御。如香港牡蛎(Crassostreahongkongensis)在受到溶藻弧菌等病原菌感染后, 其体内Akt基因表达上调[49]。研究发现海胆中第二大免疫受体家族是NOD-like受体(NLRs)，而在哺乳动物中的NLRs通过充当支架蛋白来组装蛋白复合物，从而激活NFkB和MAPK信号通路，在免疫应答中发挥作用[6]。在刺参研究中同样发现AjNLRP10蛋白参与了对细菌感染的先天免疫应答[50]。这些都表明了本研究中miRNA通过靶向这些免疫相关通路参与刺参响应高温、低氧胁迫下的信号传导、氧化应激以及免疫反应等过程。

4 结语

本研究以高温组、低氧组和对照组的体腔液细胞为研究对象，利用Illumina测序平台构建了9个miRNA文库，有效补充了刺参体腔液细胞的miRNA表达谱，并筛选了在高温胁迫和低氧胁迫下差异表达的miRNAs。基于差异表达的miRNAs，以三代全长转录组为参考文库预测了其候选靶基因，KEGG富集分析表明靶基因主要参与了刺参应激、存活和免疫等重要生物学过程，进一步揭示了差异表达的miRNAs在刺参应对高温胁迫和低氧胁迫过程中的重要调控作用。后期将利用功能实验验证关键miRNAs及其靶基因的互作关系，进一步探讨关键差异miRNA在刺参应对高温、低氧胁迫过程中的调控作用，以期为刺参应对高温、低氧逆境关键分子标记的筛选和刺参抗逆品种与品系的选育提供理论指导。