0 引言
海洋无脊椎动物体内维持适量的锰，对其机体的酶活性和代谢过程都至关重要[6]。锰作为超氧化物歧化酶等关键酶的辅助因子，参与多种调控过程，并参与能量生成和生物分子合成等核心代谢途径的酶促反应，对维持无脊椎动物生理稳态、骨骼形成及整体健康具有关键基础性作用。而环境中过量的锰，可能对动物的早期发育、生长、代谢及生理功能产生不利影响[7-8]。目前已有研究针对甲壳类和软体动物开展了锰的生物毒性效应相关实验[9-11]。转录组学技术是解析水生生物逆境响应分子机制的重要手段[12-13]，然而现有与锰相关的研究主要集中于锰对海洋无脊椎动物的毒性和生理效应，关于在锰生理剂量与毒性剂量下生物体产生差异响应的转录调控机制仍鲜有报道[14-15]。
1 材料与方法
1.1 实验设计
实验用小刺蛇尾于2023年采集自中国青岛海域（湿重0.61±0.24 g，盘径4.12±0.71 cm），在实验室水族箱中暂养3~5 d以适应环境。将无水氯化锰（分析纯）溶解于超纯水中制备氯化锰储备溶液（Mn2+质量浓度为10 g/L），并用人工海水配制成Mn2+质量浓度梯度为0、10、20、30、40、50、60、70、80、90和100 mg/L的海水溶液（NaCl质量浓度为35 g/L， pH为8.0）。在各梯度海水中分别放入小刺蛇尾6只，培养96 h。实验期间，使用直径为3 mm的小塑料管持续泵入空气，水温控制在（13±2） ℃，避光，并根据水分蒸发量补充相应的去离子水以保持海水的盐度和Mn2+质量浓度不变。在暴露实验开始后的48 h和96 h两个时间点，分别统计各质量浓度组小刺蛇尾的死亡数量（以对机械刺激无反应作为死亡判定标准），并采用线性回归法计算48 h半数致死质量浓度（48 h-LC50） 和96 h半数致死质量浓度（96 h-LC50）。通过比较两个时间点的LC50变化趋势，评估Mn2+毒性效应随暴露时间延长而增强的程度。
为探究亚致死质量浓度下小刺蛇尾对锰离子的分子响应，开展转录组实验，旨在通过测序分析对照组与暴露组个体的基因表达差异，以揭示锰胁迫下的分子响应机制。参考急性毒性实验结果，选择10 mg/L作为暴露质量浓度——该质量浓度远低于96 h-LC50，但已位于显著生理响应区间，可确保在个体不死亡的前提下获得稳定的转录水平变化。暴露组与对照组各设3个生物学重复样，编号分别为T1~T3和C1~C3，暴露条件与急性毒性实验一致。实验结束后，收集存活个体并立即用液氮固定，用于后续RNA提取与测序分析。
1.2 RNA提取、cDNA文库构建及Illumina HiSeq测序
RNA提取、cDNA文库构建及Illumina HiSeq测序的实验操作按照文献[18]中的方法进行。首先，使用TRIzol®试剂（Invitrogen）按照说明提取组织总RNA，之后采用DNase I（TaKaRa）去除基因组DNA。随后，使用2100 Bioanalyser（Agilent Technologies）和NanoDrop-2000（NanoDrop Technologies）进行RNA质量检测和定量，选取高质量RNA样本（OD260/OD280=1.8~2.2，OD260/OD230≥2.0，RIN≥6.5，28S∶18S≥1.0，总量>10 μg）用于后续分析。
取1 μg总RNA，使用TruSeqTM RNA试剂盒（Illumina，美国）构建cDNA文库。首先，采用oligo（dT）磁珠富集mRNA，随后用片段化缓冲液纯化并片段化。之后，按照Illumina标准流程进行cDNA合成、末端修复、加A尾以及连接Illumina索引接头。在2% 低分子量超琼脂糖凝胶（Bio-Rad，美国）上筛选200~300 bp的cDNA目标片段，并使用Phusion DNA聚合酶（New England Biolabs，美国）进行15个循环的PCR扩增。对文库进行TBS380定量后，所有样本交由杭州明科生物技术有限公司，使用Illumina HiSeq平台进行双端测序（PE 2×150 bp）。
1.3 从头组装与注释
使用Trimmomatic对原始双端测序数据进行修剪和质量控制[18]。去除包含未知核苷酸（N含量超过5%）或低质量（Phred值<28）的读段（reads），修剪过程中设置最小读长阈值为75 bp。采用Trinity软件（http://trinityrnaseq.sourceforge.net/）对高质量读段进行从头转录组组装。所有组装获得的转录本均通过Blastx v2.7.1进行注释（E值阈值=1.0×10-5），注释数据库包括NCBI非冗余蛋白数据库（NCBI Non-Redundant Protein Sequence Database，NR）、基因本体论（Gene Ontology，GO）、直系同源蛋白簇（Clusters of Orthologous Groups of proteins，COG）、Swiss-Prot蛋白数据库以及KEGG通路数据库（Kyoto Encyclopedia of Genes and Genomes，KEGG， http://www.genome.jp/kegg/）。
1.4 差异表达基因分析与功能富集
为鉴定对照组与暴露组间的差异表达基因（differentially expressed genes， DEGs），采用FPKM（Fragments Per Kilobase of exon model per Million mapped reads）标准方法计算各转录本的表达量[18]。具体而言，通过RSEM软件（http://deweylab.github.io/RSEM/）定量基因及异构体表达水平，并基于其定量结果，计算基因表达差异倍数（Fold Change）及统计显著性， p值<0.05且|log2（差异倍数）|≥1的基因定义为显著差异表达基因。进一步采用Goatools（https://github.com/tanghaibao/Goatools）和KOBAS（http://kobas.cbi.pku.edu.cn/home.do）分别对差异表达基因进行GO功能富集分析（以下简称GO分析）和KEGG通路分析。富集分析以全转录组基因集为背景，多重假设检验校正后p值<0.05的功能类别或通量被判定为显著富集。
1.5 基于实时定量PCR（qRT-PCR）的差异表达基因验证
为验证RNA测序（RNA-Seq）结果的可靠性，采用StepOne Plus实时PCR系统（Applied Biosystems，美国）和SYBR Green Master mix（Promega，美国）对筛选出的差异表达基因进行定量PCR分析。以β-actin作为内参基因[18]。PCR反应体系总体积为20 μL，包含：10 μL 2×SYBR Green I Master Mix、各1.0 μL 10 μmol/L基因特异性正反向引物（正向：5'-TGAGTCTCTTTCGGCAGTTCG-3'；反向：5'-GGTCC-AGCAATCGTCCATC-3'）、6 μL DEPC水和2 μL cDNA模板。扩增程序设置为：95 ℃预变性2 min；随后进行45个循环，每个循环包括95 ℃变性10 s、58 ℃退火10 s，72 ℃延伸20 s三个步骤。反应结束后，采用2-ΔΔCt法计算目标基因的相对表达量。所有样本均设置3个技术重复[19]。
2 结果
2.1 小刺蛇尾行为与生理变化
在实验条件下，无Mn2+海水中的小刺蛇尾表现出较好的活力。在Mn2+质量浓度大于10 mg/L的环境中，小刺蛇尾出现不同程度的反应迟钝、运动能力丧失，并伴有变色、蜷缩、断腕等症状。大部分小刺蛇尾在Mn2+质量浓度为50~70 mg/L的海水中暴露48 h内出现翻身困难，当Mn2+质量浓度达到90 mg/L以上时，所有实验小刺蛇尾在48 h内丧失翻身能力。当Mn2+质量浓度大于50 mg/L时，小刺蛇尾在48 h内即出现死亡。分析结果表明，Mn2+对小刺蛇尾的48 h半数致死质量浓度为61.6 mg/L（R2=0.94），96 h半数致死质量浓度为58.2 mg/L（R2=0.89）。
2.2 测序数据与组装结果
6个实验个体转录组测序共获得278 649 764条原始读段，质控后得到267 862 286条有效读段（clean reads），其中对照组有133 175 912条，锰暴露组有134 686 374条，各组样本测序数据详见表1。质量控制显示，对照组与锰暴露组的平均错误率均为0.015%。本研究中产生的原始转录组测序数据已保存于NCBI SRA数据库（登录号：SRX19924348， SRX19924349，SRX19924350，SRX19924351，SRX19924352，SRX19924353）。
表1 转录组测序数据统计Tab.1 Summary of transcriptome sequencing datasets |
| 样品号 | 原始读段数/条 | 原始碱基数/bp | 有效读段数/条 | 有效碱基数/bp | 有效碱基比例/% | Q20/% | Q30/% | GC含量/% |
|---|---|---|---|---|---|---|---|---|
| T1 | 47 008 152 | 7 051 222 800 | 45 208 780 | 6 753 508 368 | 95.78 | 98.93 | 96.39 | 36.26 |
| T2 | 44 955 096 | 6 743 264 400 | 43 261 720 | 6 461 804 011 | 95.83 | 98.85 | 96.07 | 36.45 |
| T3 | 48 034 842 | 7 205 226 300 | 46 215 874 | 6 903 036 744 | 95.81 | 98.87 | 96.14 | 36.68 |
| C1 | 44 488 222 | 6 673 233 300 | 42 734 038 | 6 382 807 089 | 95.65 | 98.86 | 96.10 | 37.10 |
| C2 | 45 459 166 | 6 818 874 900 | 43 701 244 | 6 528 372 113 | 95.74 | 98.92 | 96.39 | 36.96 |
| C3 | 48 704 286 | 7 305 642 900 | 46 740 630 | 6 982 389 487 | 95.58 | 98.97 | 96.46 | 37.10 |
注：Q20表示碱基质量值大于等于20的碱基所占百分比，Q30表示碱基质量值大于等于30的碱基所占百分比。 |
如表2所示，转录组从头组装结果显示：测序读段共组装获得324 097条转录本（transcripts），总长度为216 359 856个碱基；得到125 817条单基因，总长度为77 768 182个碱基。单基因的最大长度为17 697 bp，平均长度为618.11 bp。N50值可用于评估组装连续性（其数值越高，通常表明组装完整性与连续性越好），在本研究中单基因的N50值为612 bp。在所有单基因中，分别有20 857条（16.58%）、11 852条（9.42%）、17 828条（14.17%）、11 810条（9.39%）和13 886条（11.04%）成功在NR、GO、COG、KEGG和Swiss-Prot数据库中获得注释，其中，有21 757条（占单基因总数的17.29%）被至少一个数据库注释（表3）。
表2 转录本与单基因的质量及长度统计Tab.2 Quality and length statistics of the transcripts and unigenes |
| 项目 | 值 |
|---|---|
| 单基因总数/条 | 125 817 |
| 单基因总长度/bp | 77 768 182 |
| 转录本总数/条 | 324 097 |
| 转录本总长度/bp | 216 359 856 |
| 单基因平均长度/bp | 618.11 |
| 最长单基因长度/bp | 17 697 |
| N50/bp | 612 |
注：N50是指将转录本按长度从长到短排序后，当转录本长度累加达到总长度的50%时，所对应的单个转录本长度值。 |
表3 单基因在不同数据库中的注释结果Tab.3 Annotation of unigenes in different databases |
| 数据库 | 注释单基因数 | 占单基因总数的比例/% |
|---|---|---|
| NR | 20 857 | 16.58 |
| GO | 11 852 | 9.42 |
| COG | 17 828 | 14.17 |
| KEGG | 11 810 | 9.39 |
| Swiss-Prot | 13 886 | 11.04 |
| 全部数据库共有 | 7 602 | 6.04 |
| 至少一个数据库注释 | 21 757 | 17.29 |
2.3 差异表达基因注释
锰暴露组与对照组间共鉴定到3 756个差异表达基因（DEGs）。锰暴露组（Mn2+质量浓度为10 mg/L）与对照组之间差异表达基因的火山图显示，Mn2+暴露后小刺蛇尾下调基因数量约为上调基因数量的1.8倍（图1）。
图2为锰暴露组与对照组差异表达基因的GO功能注释柱形图。上方横坐标（0%~100%）表示注释到某一GO条目的差异表达基因数占所有被注释的差异表达基因总数的比例；下方横坐标以分数形式表示具体数量，其中分子为注释到该GO条目的差异表达基因数，分母为用于富集分析的背景基因总数（即该物种所有被注释的基因数）；纵坐标为按照GO三大类别（生物过程、细胞组分、分子功能）排列的具体功能条目。GO富集分析结果显示，在生物过程中，注释基因集中于单生物体过程（65条）、细胞过程（64条）、代谢过程（54条）、多细胞生物体过程（52条）、生物调控（49条）等多个条目；在细胞组分中，注释基因主要集中于细胞组分（62条）、细胞（62条）、细胞器（51条）、细胞器组分（40条）、大分子复合物（31条）等条目；在分子功能中，注释基因主要集中于结合（37条）、催化活性（28条）、结构分子活性（17条）、分子转导活性（12条）、信号转导活性（11条）等条目。
2.4 差异表达基因的GO富集分析
在设定的显著性阈值下（p<0.05），共有3 756个差异表达基因被注释到630个不同的GO条目中。这些条目可进一步划分为：458个生物过程条目、92个细胞组分条目和80个分子功能条目。图3展示了富集最显著的30个GO条目的气泡图，图中横坐标为富集因子（rich factor），即注释到某GO条目的差异表达基因数与该GO条目中总注释基因数的比值；纵坐标为GO条目名称；气泡颜色代表p值（颜色越红表示富集越显著）；气泡面积代表注释到该GO条目的差异表达基因数量。在排名前30的GO条目中，未出现分子功能类别的条目，它们全部属于生物过程和细胞组分类别。
在生物过程类别中，富集最显著的条目包括：细胞组分解聚（差异表达基因数为54）、病毒生命周期（51）、蛋白质复合体解聚（51）、翻译终止（49）以及共翻译蛋白质靶向膜（48）等。其中“细胞组分解聚”的差异表达基因数量最多，但其富集因子相对较低（约为0.11），提示该条目在背景基因组中本身规模较大，但差异表达基因的富集程度并不特异。在细胞组分类别中，前30个条目中，共有7个显著富集的条目，包括：胞质核糖体大亚基（43）、内吞小泡（36）、胞质核糖体小亚基（32）、外部包被结构（31）、细胞壁（31）、吞噬小泡（30）和前核糖体（15）。由图3可以看出，富集显著性与富集因子之间并非严格对应。p值受GO条目总基因数的影响较大，总基因数较多的条目，即使差异表达基因占比不高，也易获得极显著的p值。相比较而言，富集因子更能反映某一GO条目对Mn2+胁迫的特异性响应程度。部分低丰度通路虽然差异表达基因的绝对数量较少，但对胁迫表现出了高度的特异性响应。
2.5 KEGG通路富集分析
在显著性阈值以下（p<0.05），差异表达基因共富集到229条KEGG通路（图4展示了富集最显著的前30条通路）。这些通路可系统归类为：代谢通路、遗传信息处理、环境信息处理、细胞过程、有机体系统及疾病相关通路等六大类。
在细胞稳态层面，线粒体自噬调控通路与萜类骨架（辅酶Q10前体）生物合成通路的显著富集提示线粒体功能受损可能为锰离子暴露下的早期核心事件。在信号转导层面，Hedgehog、Jak-STAT、Hippo等多条发育相关通路被扰动，提示细胞通讯与组织模式调控受到干扰。在代谢层面，原核生物碳固定、甘油脂、脂肪酸、2-氧羧酸及硫胺素代谢等基础能量与物质代谢通路发生广泛重塑，同时，硒化合物、牛磺酸、氰基氨基酸及多种生物碱（异喹啉、莨菪烷类）代谢通路显著富集，反映出小刺蛇尾激活多层次抗氧化与解毒机制。群体感应、细菌感染以及外源物降解（氯烷、香叶醇、芳香族化合物）通路显著富集，暗示了小刺蛇尾微生物防御屏障的破坏及广谱外源物代谢应激的诱导。其他免疫、感觉、内分泌及排泄系统相关通路的广泛富集，表明锰暴露引发的是全身性的、多器官的复杂应激状态。
部分锰离子暴露所涉及的信号通路途径不属于富集最显著的通路，但这些通路的生物学意义不容忽视，这些通路包括胰岛素信号通路（差异表达基因数为8）、神经营养因子信号通路（7）、GnRH信号通路（6）、MAPK信号通路（5）、PI3K-Akt信号通路（4）、HIF-1信号通路（2）和FoxO信号通路（1）等。与神经退行性疾病相关的通路以及线粒体自噬、氧化磷酸化等与线粒体功能相关的通路中部分基因在不同物种（如人、小鼠、大鼠、果蝇等）中高度保守，如本研究中表达水平显著上调的PINK1基因，与帕金森病的发生高度相关。
2.6 实时定量PCR验证差异表达基因
3 讨论
已有研究表明，环境中过量的锰离子会对海洋生物的运动能力和感知能力造成影响，如高浓度的锰会损害龙虾的触角摆动能力，导致其迷失方向和无法找到食物[23]。在本实验中，锰离子对小刺蛇尾的影响首先体现在运动能力的丧失以及感知和反应能力大幅削弱。
研究结果表明，锰离子对小刺蛇尾的生理代谢具有多重复杂影响。与对照组相比，锰暴露组中大部分差异表达基因呈现下调趋势（2 409/3 756）。GO富集结果中的条目可划分为外部防御结构、内部隔离与解毒、能量分配与生长、蛋白处理与应激、遗传信息调控等不同功能模块，表明小刺蛇尾对锰暴露的响应是多系统、多层次的复杂生理过程。在外部结构防御层面，个体可能通过改变细胞壁、几丁质结构或加速蜕皮，将重金属离子阻挡在体外或直接摆脱已吸附重金属的体表结构；在内部隔离与解毒上，则可能通过内吞和吞噬作用形成囊泡，降低细胞内游离锰离子质量浓度以减轻毒性；而“解组装”过程既可能是细胞器受损的表现，也可能是一种主动回收利用受损组件的适应性过程；在能量分配上，个体可能通过减少核糖体生成和蛋白质合成以减少能量消耗，从而将有限的能量重新分配用于激活解毒、修复和维持基本生存的代谢活动。胰岛素样信号通路是无脊椎动物中重要且古老的调控通路之一，其核心功能涉及代谢、能量分配和生长调控。本研究中小刺蛇尾在锰暴露下胰岛素样信号通路的富集，可能表明个体正在重新配置其能量代谢的模式以应对胁迫。
在人类与啮齿类的研究中已证实，锰与运动障碍和认知功能障碍密切相关。锰易在脑部多巴胺富集区域蓄积，过量摄入可引发类帕金森氏症运动障碍[24]。NICOLAI等[25]在线虫（Caenorhabditis elegans）中发现帕金森病相关同源基因与锰中毒存在关联。本研究发现帕金森病相关基因PINK1显著上调。PINK1 基因在进化上极为古老，其感应线粒体损伤、启动线粒体自噬的核心功能在所有后生动物中高度保守。在人类中，线粒体损伤会激活PINK1，而该基因的突变会导致线粒体功能紊乱，特别是在多巴胺神经元中，最终引发早发性帕金森病。本研究结果提示，小刺蛇尾在锰暴露下受到了广泛的线粒体损伤，从而激活了由PINK1上调介导的、这一高度保守的细胞防御机制。小刺蛇尾在线粒体受损、蛋白质合成抑制、氧化应激和炎症反应的多重作用下，最终表现出运动能力的大幅下降。需要特别指出的是，深海采矿过程中另一关键环境因素是沉积物羽流导致的再悬浮与覆盖效应。若深海动物因锰暴露导致运动能力下降，将显著降低其逃离此类覆盖层的可能。因此，在评估深海采矿的环境影响时，需综合考虑多种因素的协同作用。
小刺蛇尾在锰暴露下的线粒体功能损伤是深刻而多层次的。转录组分析显示，除了直接与线粒体自噬调控相关的PINK1外，还包括大量与能量代谢、抗氧化防御、细胞凋亡以及线粒体生物发生和动力学相关的DEGs。例如，注释到电子载体活性的基因发生显著变化，包括细胞色素C、辅酶Q及线粒体电子传递链复合物的多个亚基。这些与线粒体功能相关的DEGs的显著变化表明，Mn2+会影响小刺蛇尾的呼吸电子传递、ATP合成及产热过程。ALAIMO等[26]和SARKAR等[27]的研究均指出，锰可激活线粒体凋亡途径、导致线粒体形态蛋白表达失衡、破坏线粒体网络结构、损害ATP生成能力，并通过线粒体缺陷驱动炎症小体信号放大。然而，当前关于锰离子暴露所引发的炎症与应激反应的分子机制研究仍较有限。
过量锰离子可破坏小刺蛇尾的神经递质系统，GO分析结果显示，与突触结构和功能相关的基因集发生显著变化，包括“突触组分”和“参与突触传递的突触前过程”等条目，表明锰离子毒性直接靶向了神经突触这一关键结构，可能直接干扰了突触囊泡的胞吐作用以及神经递质受体的功能， 从而导致信号传递效率的急剧下降。此外，注释到转运蛋白活性的基因表达变化，提示其可能影响了突触间隙内神经递质（如谷氨酸、GABA）的再摄取效率，破坏递质作用的时空精确性。无脊椎动物的核心神经递质系统，如胆碱能、谷氨酸能、GABA能、多巴胺能系统等在进化上高度保守，锰离子可能通过直接毒性和间接的细胞能量剥夺，破坏神经递质系统的正常功能，这或许是导致其在实验中对外界刺激反应能力下降及出现翻正反射障碍等运动能力丧失表现的重要原因。
4 结论与展望
锰作为重要的工业原料虽具有广泛应用价值，但其潜在的环境危害不容忽视。目前，针对锰污染生态风险的认识仍存在明显不足。本研究通过分析小刺蛇尾在Mn2+暴露下的转录组响应，发现实验组与对照组间存在3 756个差异表达基因（DEGs），其中1 347个基因上调、2 409个基因下调。这些DEGs涉及多种功能类别，揭示了过量锰离子对小刺蛇尾多层次的复杂影响：一方面，小刺蛇尾运动协调性下降会直接影响个体的摄食效率和面对捕食者的躲避能力；另一方面，神经递质系统的紊乱可能改变其昼夜节律和繁殖行为。值得注意的是，深海环境中锰离子往往与其它金属污染物共存，其协同毒性效应可能比单一金属的毒性效应更为显著。因此，我们建议后续研究从生理、代谢及基因表达等多维度深入探究锰暴露对海洋生物的影响机制，着重关注长期低剂量锰暴露的累积效应、锰与其他污染物（如钴、镍）的联合毒性及个体水平行为变化与种群动态之间的关联机制。此外，应推动建立基于转录组标志物的生物监测体系，为构建相关环境阈值标准提供科学依据。