基于转录组测序分析绿鳍马面鲀脑组织响应密度胁迫的分子机制

段晓晨，程起群

(1.中国水产科学研究院东海水产研究所，上海 200090；2.上海海洋大学水产与生命学院，上海 201306)

绿鳍马面鲀(Thamnaconusseptentrionalis)，隶属于鲀形目(Tetraodontiformes)，单角鲀科(Monacanthidae),马面鲀属，又称马面鱼、剥皮鱼等，是我国沿海较为常见的经济鱼类[1]。经过近40年的高强度开发，绿鳍马面鲀资源量迅速衰减，目前已不能充分满足市场需求[2]。以“养”代“捕”，是缓解野生资源压力、保障绿鳍马面鲀市场供给的有效途径。2011年，绿鳍马面鲀的人工繁育技术取得突破[3]，之后开始应用于工厂化养殖。

密度过高是水产养殖活动中的重要限制因子[4]，密度胁迫会使种群内对空间、食物和氧气等资源的竞争加剧，基因表达和适应能力受到影响[5]。如SUN等[6]的研究表明，拥挤可以显著诱导鱼的免疫反应，而长期处于密度胁迫下会损害鱼类的免疫能力；QIANG等[7]的研究显示，高密度条件会改变肌肉脂肪酸的组成，刺激脂质代谢酶显著上调。LIU等[8]发现，高密度养殖导致多种消化酶的活性显著降低，进而影响生物的免疫能力。目前，关于绿鳍马面鲀脑组织对密度胁迫的分子响应机制尚不明确。

转录组学是揭示生物体内分子机制的重要手段[9]，应用转录组测序技术可以对水产动物进行高通量测序，通过测序结果进一步实现基因表达情况分析和基因功能注释，明确特定条件下相关基因的表达情况[10]，目前转录组分析被普遍应用到水产动物的生物学研究中。NORIKO等[11]发现，背景颜色不仅影响平均红细胞血红蛋白量(mean corpuscular hemoglobin，MCH)水平，而且影响大脑中的促性腺激素释放激素II(chicken-type gonadotropin-releasing hormone II，cGnRH-II)水平，并表明cGnRH-II可能在调节比目鱼大脑中的 MCH 神经功能、食物摄入方面发挥作用。周江等[12]研究发现，星斑川鲽(Platichthysstellatus)脑组织中除促卵泡激素(follicle-stimulating hormone, FSH)外，促性腺激素(gonadotropins hormone, GnH)、促性腺激素释放激素(gonadotropin-releasing hormone, GnRH)以及催乳素释放肽受体(prolactin-releasing peptide-receptor, PrRPR)等均在卵巢退化期时期个体的脑组织中上调表达。刘怡南等[13]发现，igdcc3、prlh、mch、pitx2等基因在七彩神仙鱼(Symphysodonharaldi)的脑组织中呈性别差异表达，可能在该鱼脑组织类固醇激素形成及配子发生的调节过程中发挥重要作用。

本研究利用 Illumina 平台对绿鳍马面鲀的脑组织进行转录组测序，再对测序所得数据进行质控、比对、差异性分析，然后再进行GO(Gene Ontology)富集分析、KEGG(Kyoto Encyclopedia of Genes and Genomes)富集分析、趋势分析，以期找到绿鳍马面鲀脑组织响应密度胁迫的差异表达基因和关键通路,为了解绿鳍马面鲀响应密度胁迫的分子机制及绿鳍马面鲀的集约化养殖提供参考依据。

1 材料与方法

1.1 实验材料

实验用绿鳍马面鲀取自中国水产科学研究院东海水产研究所赣榆研究中心，系人工自繁。选取体质健壮、规格较为一致且平均体质量为(9.5± 0.5)g的个体进行实验。

实验前，将绿鳍马面鲀放在水容量40 m3的室内水泥池中暂养一周，密度为50尾·m-3，将此池样本作为对照组；实验时，从大池将鱼分至水容量1 m3的聚乙烯微流水养殖桶中，设置2个密度试验组，即100尾·m-3的中密度组和500尾·m-3的高密度组，每个密度设置3个重复。

实验过程中，水温保持(17±1)℃，pH为8.20±0.90，采用微流水饱和充气养殖，溶氧保持在(6.7±0.4)mg·L-1，换水率120%。每日于8∶00、12∶00、15∶00分别投饵1次，总投饵量为鱼平均体质量的3%。实验周期50 d。

1.2 组织样本采集

首先，取对照组大池中的绿鳍马面鲀脑组织；然后，在实验开始后的第25天和第50天时，分别在两个密度试验组的养殖桶中取3尾鱼，用丁香酚对实验鱼进行麻醉、致晕、致死，然后将鱼捞起放置于干净的解剖盘中，用酒精擦拭鱼体表面，随后用解剖剪剪开鱼脑上方的骨骼，用镊子取出脑组织。脑组织取出后，快速切成厚度<0.5 cm的任意片状，迅速浸入10倍体积 RNAwait 液(Solarbio)中，然后按照RNAwait液的使用说明，将其4 ℃暂存一昼夜，再放置于-20 ℃冰箱冷冻保存，带回实验室进行后期实验。为便于后续的叙述，各组代号的含义如下：1)CB指对照组的脑组织；2)D1B指第25天时100尾·m-3的脑组织；3)D2B指第25天时500尾·m-3的脑组织；4)d1B指第50天时100 尾·m-3的脑组织；5)d2B指第50天时500尾·m-3的脑组织。

1.3 RNA的提取与检测、文库构建及测序

取出冷冻的绿鳍马面鲀脑组织，采用Trizol法对脑组织进行RNA提取，检测合格后得到该鱼脑组织的总RNA。通过带有Oligo(dT)的磁珠富集具有polyA尾巴的真核mRNA后，用缓冲液把mRNA打断。以片段化的mRNA为模版、随机寡核苷酸为引物，在M-MuLV逆转录酶体系中合成cDNA第一条链，随后用RNaseH降解RNA链，并在DNA polymerase I 体系下,以dNTPs为原料合成cDNA第二条链。纯化后的双链cDNA经过末端修复、加A尾并连接测序接头，用AMPure XP beads筛选200 bp左右的cDNA，进行PCR扩增并再次使用AMPure XP beads纯化PCR产物，最终获得文库。构建好的文库利用Illumina HiSeq平台进行测序。

1.4 测序数据组装及基因功能注释

将测序获得的原始数据(raw reads)利用fastp[14]进行质控、去除含接头reads、去除含N比例大于10%的reads、去除全部都是A碱基的reads、去除低质量reads(质量值Q≤20的碱基数占整条read的50%以上)。得到高质量干净序列(clean reads)进行后续分析。将质控后的过滤数据利用HISAT2软件与绿鳍马面鲀的参考基因组进行比对[15]，以便进行后续的基因表达和差异分析。

1.5 基因表达和差异分析

本研究采用 FPKM(expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced)方法计算基因的表达量[16]。再用DESeq2[17]软件包筛选差异表达基因, 筛选阀值为Fold change>2和Q值<0.05，且Q值越小基因表达量差异越显著[18]。最后，进行差异表达基因的GO和 KEGG富集分析(false discovery rate, FDR< 0.05)。

2 结果与分析

2.1 脑转录组测序数据

运用 Illumina HiSeq 高通量测序平台分别对不同密度组绿鳍马面鲀脑组织进行了转录组测序，共得到100.3 Gb RawData，过滤后得到98.9 Gb CleanData。碱基质量及组成分析显示，每组GC含量区间为48.46%～50.38%，各样品Q20的碱基质量值比例均不小于96.67%，各样品Q30的碱基质量值比例均不小于91.36%(表1)。

表1 绿鳍马面鲀脑转录组测序数据统计

2.2 基因表达量总体分布

对照组和密度胁迫组绿鳍马面鲀脑组织间比较能反映它们在基因表达水平上的基本差异(图1-a～d)。图1中所示每个点代表一个基因，横坐标表示两个分组间的差异倍数对数值，纵坐标表示两个分组差异的FDR的-Lg值，红色(group_2相对于group_1表达量上调)和蓝色(表达量下调)的点表示基因的表达量有差异(判断标准为 FDR <0.05, 且差异倍数在两倍以上)，黑色的点为没有差异。其中，越靠近两端的基因，差异程度越大。在CB-vs-D1B中，有197个差异基因上调，418个差异基因下调；在CB-vs-D2B中，有114个差异基因上调，301个差异基因下调；在CB-vs-d1B中，有317个差异基因上调，491个差异基因下调；在CB-vs-d2B中，有311个差异基因上调，608个差异基因下调(图1-e)。

图1 密度胁迫组和对照组差异基因表达状况

2.3 GO功能富集分析

将转录组信息与 GO(Gene Ontology)数据库作比对，分别注释到了生物过程(biological process, BP)、分子功能(molecular function, MF)、细胞组分(cellular component, CC)3个亚类(图2)，其中红色代表上调表达，蓝色代表下调表达。在生物过程亚类中，得到最显著注释的条目依次为：细胞过程、单组织过程、代谢过程、生物调控、生物过程的调节；在分子功能亚类中, 得到最显著注释的条目依次为：结合、催化活性；在细胞组分亚类中, 得到最显著注释的条目依次为：细胞组分、细胞、细胞膜。但在CB-vs-D2B中，显著差异注释的条目相较其他3组更少。

图2 差异表达基因的 GO 功能富集分析

2.4 KEGG 功能富集分析

用KEGG数据库富集出在密度胁迫后表现出显著差异的代谢途径绘制KEGG富集气泡图(图3)，其中Q值为多重假设检验校正之后的P值。利用Q值最小的前20个pathway来作图，纵坐标为pathway，横坐标为富集因子(该pathway中差异基因数量除以该pathway中所有数量)，大小表示数量多少，颜色越红Q值越小。富集分析的结果表明：在CB-vs-D1B中，差异表达的基因富集到267条通路上，其中富集程度最高的通路依次是：氧化磷酸化、昼夜节律、产热、病毒蛋白与细胞因子和细胞因子受体的相互作用(图3-a)。在CB-vs-D2B中，差异表达的基因富集到230条通路上，其中富集程度最高的通路是：昼夜节律(图3-b)。在CB-vs-d1B中，差异表达的基因富集到278条通路上，其中富集程度最高的通路是：昼夜节律、蛋白质消化吸收(图3-c)。在CB-vs-d2B中，差异表达的基因富集到269条通路上，其中富集程度最高的通路是：昼夜节律、赖氨酸降解、MAPK信号通路、神经活性配体-受体相互作用(图3-d)。

2.5 趋势分析

趋势分析(trend analysis)是针对多个连续型样本(至少3个)的特点(样本间包含特定的时间、空间或处理剂量大小顺序)而对基因的表达模式(在多阶段中表达曲线的形状)进行聚类的方法。然后从聚类结果中挑选符合一定生物学特性(如表达量持续上升)的基因集。

2.5.1 100尾·m-3组的趋势分析

来自CB_D1B_d1B的1 520个基因被分配到8个趋势中(图4)，其中 Profile 2、Profile 1、Profile 0和Profile 6为显著趋势。在Profile 2中，有361个基因表达先下调后恢复到胁迫前的状态，KEGG通路分析显示其中最主要富集的通路是白细胞跨内皮迁移和细胞凋亡。在Profile 1中，有273个基因表达先下调后趋平，KEGG通路分析显示其中最主要富集的通路是昼夜节律。在Profile 0中，有271个基因表达持续性下调，KEGG通路分析显示其中最主要富集的通路是MAPK信号通路、昼夜节律、皮质醇合成与分泌、钙信号通路、ECM-受体相互作用。在Profile 6中，有234个基因表达情况先上调后趋平，KEGG通路分析显示其中最主要富集的通路是氧化磷酸化、产热、代谢途径、糖酵解/糖异生、碳代谢、氨基酸的生物合成、核糖体和RNA降解。

图4 CB_D1B_d1B所有趋势总图

2.5.2 500尾·m-3组的趋势分析

来自CB_D2B_d2B的1 547个基因被分配到8个趋势中(图5)，其中Profile 2、Profile 0、Profile 4和Profile1为显著趋势。在Profile 2中，有363个基因表达先下调后恢复到胁迫前的状态，KEGG通路分析显示其中最主要富集的通路是Toll样受体信号通路。在Profile 0中，有349个基因表达持续性下调，KEGG通路分析显示其中最主要富集的通路是MAPK信号通路、cAMP 信号通路、多巴胺能突触、钙信号通路和谷氨酸能突触。在Profile 4中，有178个基因表达先不变再上调，KEGG通路分析显示其中最主要富集的通路是PPAR信号通路、过氧化物酶体、氮代谢和类固醇生物合成。在Profile 1中，有168个基因表达先下调再趋平，KEGG通路分析显示其中最主要富集的通路是昼夜节律。

图5 CB_D2B_d2B所有趋势总图

3 讨论

为了满足人类社会需求，水产动物高密度养殖逐渐成为降低养殖成本、提高养殖产量的重要养殖模式，然而过度追求高密度养殖会造成种群内对资源的掠夺加剧，最终可能造成水产动物生存能力的降低[19]。当密度发生改变时，水产动物体内并非单个的基因或者信号通路发生变化，而是整个调控网络均会发生改变[20]。目前，关于绿鳍马面鲀的的研究多集中在生长发育、摄食习性、养殖技术、营养成分和遗传结构分析[21-25]。近些年，转录组技术在多种水产动物研究中均得到广泛应用，其在生物进化研究、免疫应答反应和毒理学等方面取得了大量的研究进展[26]。为了阐明绿鳍马面鲀脑组织的基因表达模式，探究明显富集通路的关键基因，本研究首次构建了绿鳍马面鲀脑组织的cDNA文库，并通过Illumina平台进行了高通量测序分析。

3.1 脑组织差异表达基因数目和GO富集

对绿鳍马面鲀脑组织的不同组进行转录组测序、组装、注释，发现在CB-vs-D1B中，有197个差异基因上调，418个差异基因下调；在CB-vs-D2B中，有114个差异基因上调，301个差异基因下调；在CB-vs-d1B中，有317个差异基因上调，491个差异基因下调；在CB-vs-d2B中，有311个差异基因上调，608个差异基因下调。表明绿鳍马面鲀在受到密度胁迫后可能激活了细胞的代谢，以应答应激环境下对绿鳍马面鲀体内造成的伤害。

在GO功能富集分析中，注释得到的最显著的条目有细胞过程、单组织过程、代谢过程、结合和催化活性等部分，表明绿鳍马面鲀脑组织的此类生理过程受到了显著影响。但与其他组相比，CB-vs-D2B中的显著差异注释的条目相较于其他3组相对较少，可能是相对其他组来说，D2B在密度胁迫下受到的影响更小。

3.2 各组中均得到富集的通路

在对绿鳍马面鲀脑组织进行KEGG分析后发现，各个组中的昼夜节律通路都得到了显著富集。昼夜节律又称日夜节律，是生物体内的生物钟面对外界信号影响时所作出的一种内源性调控应答反应机制，从而使生物的行为节律以及内分泌系统与外界信号同步, 因此是机体对外界环境信息的一种适应性应答[27]。clock基因编码的蛋白质在调节昼夜节律中起到核心作用。该蛋白具有组蛋白乙酰转移酶活性，是碱性-螺旋-环-螺旋家族的转录因子，在DNA损伤、细胞凋亡等过程具有重要作用[28]。cry是一种蛋白质编码基因，与cry相关的疾病包括延迟睡眠阶段障碍和睡眠障碍，其相关通路包括BMAL1-CLOCK、NPAS2，激活昼夜节律基因的表达与该基因相关的GO注释包括核苷酸结合和转录因子结合[29]。这些与生物体昼夜节律相关的基因异常表达可能代表着密度胁迫影响了绿鳍马面鲀脑组织对于生物钟的调控，且随着时间改变未恢复到正常水平。

3.3 趋势分析中的关键通路

在绿鳍马面鲀脑组织的趋势分析中，发现MAPK信号通路在CB_D1B_d1B和CB_D2B_d2B均持续性下调表达，是真核细胞介导细胞外信号到细胞内反应的重要信号转导系统，调节细胞的生长、分化、凋亡和死亡等多种生理过程。其中bdnf基因是编码神经生长因子蛋白质家族的成员，选择性剪接产生多个转录变体，其中至少一个编码前体蛋白，该前体蛋白经过蛋白水解加工产生成熟蛋白[30]。rps6ka3是一种蛋白质编码基因，该基因是编码丝氨酸/苏氨酸激酶的 RSK家族的成员[31]。il1r1基因编码属于白细胞介素-1 受体家族的细胞因子受体，它是参与许多细胞因子诱导的免疫和炎症反应的重要介质[32]。由此看来，MAPK信号通路的持续性下调表达对绿鳍马面鲀的蛋白质编码造成了不利影响，并可能造成该鱼的抗病免疫功能降低。

另外，钙离子信号通路也在CB_D1B_d1B和CB_D2B_d2B中持续性下调表达，钙离子信号通路可以调节有氧代谢[33]，其中htr5a是一种蛋白质编码基因，起神经递质、有丝分裂原和激素的作用，G蛋白介导这种受体的活性，它可能通过调节细胞内Ca2+水平而发挥作用，可能与精神疾病有关[34]。ryr2基因编码主要位于内质网膜上的RYR-2蛋白，RYR-2作为同四聚体的一个亚基，构成了调节细胞内钙离子浓度的钙通道，其在心脏和大脑中广泛分布和高表达[35]。pln基因可逆地抑制心脏肌质网中ATP2A2的活性，通过对ATP2A2的影响，调节心肌对生理刺激的收缩性，在肌肉放松期间调节钙的再摄取，并在心肌中的钙稳态中发挥重要作用[36]。

除去以上几个通路出现显著富集外，笔者还关注到PPAR信号通路、Toll样受体信号通路。PPAR信号通路和动物免疫防御和抗炎方面有关，PPARα调节因子影响脂蛋白和长链脂肪酸代谢[37]。推测本研究中绿鳍马面鲀在密度胁迫下能量消耗变大，随后糖酵解/糖异生增强，长时间的密度胁迫导致该鱼免疫抗病能力的降低。Toll样受体可识别不同的病原体相关分子模式，并在先天性免疫应答中扮演着不可或缺的角色。它们是抵御病原体入侵的第一道防线，在炎症、免疫细胞调控、存活和增殖方面发挥着关键作用[38]。

致谢：绿鳍马面鲀养殖实验过程中，得到本所吴艳庆、刘威等同志的帮助;绿鳍马面鲀全基因组序列由水科院黄海所陈四清研究员、边力博士提供；谨致谢忱。