1. 引言
鱼鳔是硬骨鱼特有的一个器官，它位于体腔内，在脊柱的下方，是一个主要由氧气和二氧化碳等几种气体填充的囊，通常按照有无与咽部相连接的鳔管被分为闭鳔类(physoclistous)和喉鳔类(physostomes) [1] 。在结构上，鳔主要由两层膜组成，外膜为腹膜凹陷层和胶原层，而内膜包括纤维层、平滑肌和黏膜上皮。这种结构的复杂性使鳔能够在气体分泌与重吸收、自主神经支配、辅助呼吸、调节浮力、声音产生及压力接收等众多方面发挥重要作用 [2] [3] 。Zheng等人 [4] 根据斑马鱼鳔的转录组研究表明鱼鳔主要通过上皮组织和肌肉组织发挥功能。此外，也有研究发现斑马鱼鳔的前房具有声音共振的特性，而后房则用于调节浮力 [5] 。Pelster等人 [6] 研究发现鳔会通过调节气体分泌与重吸收助力欧洲鳗鲡的长途迁徙。显而易见，不同鱼类的鳔在形态大小及功能发挥上存在很大的差异，这可能是受个体行为或生活环境的影响而产生的特异性进化。除此之外值得注意的一点是：鳔的形成和充气也是影响闭鳔鱼类人工育苗成效的重要因素之一，鳔发育异常会导致仔稚鱼死亡或脊柱前凸，育苗过程中通常以及时清除水面油膜和有机碎屑的方式助力鳔的正常充气 [7] [8] 。
河鲀，东方鲀属(Takifugu)鱼类的统称，隶属于硬骨鱼纲(Osteichthyes)，鲀形目(Tetraodontiformes)，鲀亚目(Tetraodontidei)，鲀科(Tetraodontidae)，主要分布于日本、朝鲜半岛和中国沿海，常见种类有暗纹东方鲀(Takifugu obscurus)、红鳍东方鲀(Takifugu rubripes)、菊黄东方鲀(Takifugu flavidus)、假睛东方鲀(Takifugu pseudommus)等，其中暗纹东方鲀和红鳍东方鲀分别是中国南方和北方仅有的两个可以合法养殖加工且出口占比最大的品种 [9] 。研究表明红鳍东方鲀与暗纹东方鲀亲缘关系十分密切，其基因组资源可用于暗纹东方鲀的研究 [10] 。然而有趣的是，虽然它们的亲缘关系十分密切，但其生存环境却大相径庭：红鳍东方鲀属于海养型鱼类，而暗纹东方鲀属于江海洄游型鱼类，每年10月温度下降到18℃左右时在长江中下游自然温度条件下无法正常生存，需下海在温暖洋流中越冬 [11] 。邓捷春 [12] 等研究发现呈味游离氨基酸和无机离子含量的不同是造成淡水型暗纹东方鲀与海养型红鳍东方鲀滋味差异的主要因素；Gao等 [13] 证实了由雌性暗纹东方鲀和雄性红鳍东方鲀培育出的杂交鲀相比暗纹东方鲀具有显著生长优势；Wang等 [14] 证实了由两种东方鲀皮肤提取的胶原蛋白可以作为哺乳动物来源的胶原蛋白的可行替代品。可以看出，与两种东方鲀相关的研究已全面展开，然而关于其鳔组织的研究却少之又少。与此同时，通过前文对鳔的总结我们可以发现不同鱼类的鳔通常会进化出不同的结构，发挥不同的功能，这说明鳔会因鱼类生存环境的差异而产生适应性进化。那么在亲缘关系密切且生存环境差异巨大的暗纹东方鲀与红鳍东方鲀的鳔之间又有何种不同，又是否存在功能上的差异，值得我们探究。
由此，我们利用RNA-Seq技术对暗纹东方鲀和红鳍东方鲀的鳔进行了比较。根据差异表达基因分析，在转录水平简单探究了两种东方鲀的鳔之间的功能差异，并推测这种差异可能与两种东方鲀在不同生存环境下产生的适应性进化有关。同时本研究也可为其他生物鳔组织的研究提供参考基础。
2. 材料与方法
2.1. 样品处理
动物实验得到了大连海洋大学动物保护与利用委员会的批准。
健康的暗纹东方鲀(30尾)和红鳍东方鲀(30尾)分别来自江苏中洋集团股份有限公司和大连富谷水产有限公司，组织采样前对它们分别进行了为期一周的驯化养殖，并在采样前24小时停止喂食。然后以暗纹东方鲀的鳔(TOSB)为实验组，红鳍东方鲀的鳔(TRSB)为对照组，每组各随机选取12尾鱼，用80 mg/L浓度的三卡因甲磺酸盐(MS-222, Sigma, Redmond, WA, USA)麻醉鱼后解剖，取鳔组织置于液氮中快速冷冻，−80℃保存。
2.2. 文库制备与转录组测序
采用TRIzol试剂(Invitrogen, California, CA, USA)从每个样品中提取总RNA，再将每4条鱼的鳔组织的总RNA等量混合为一组。使用紫外分光光度计(IMPLEN, California, CA, USA)检测RNA纯度，并使用Agilent 2100生物分析仪(Agilent Technologies, Santa Clara, CA, USA)精确检测RNA完整性。最终从24条鱼的鳔组织中构建了6个测序文库，分别为3个暗纹东方鲀组(TOSB1、TOSB2、TOSB3)和3个红鳍东方鲀组(TRSB1、TRSB2、TRSB3)。最后通过Illumina Novaseq 6000 (Illumina,San Diego, CA, USA)进行双端(PE)测序。
2.3. 转录组数据处理
采用FastQC对原始数据进行质控，Trimomatic v0.38被用于过滤测序引物接头、ploy-N和低质量读数的读段(reads)，从而获得干净的读段(clean reads)，所有下游分析都是基于高质量的clean reads进行的。同时，从美国国家生物技术信息中心(National Center for Biotechnology Information, NCBI)下载得到了红鳍东方鲀的参考基因组(fTakRub1.2, www.ncbi.nlm.nih.gov/genome/63)。然后使用Hisat2v2.1.0将获得的clean reads比对到参考基因组上，再通过SAMtools v1.7进行排序去重，使用HTSeq v0.11.2统计比对到每个基因上的read counts。最后采用R语言包DESeq2 v3.10进行差异表达基因分析，并以p.adj < 0.05 & |log2foldchange| ≥ 1作为筛选差异表达基因(Differentially expressed genes, DEGs)的阈值。通过R语言包pheatmap v1.0.12绘制差异基因聚类热图，通过运行python脚本进行主成分分析，此脚本已在我们以前的文章中发表 [15] 。
2.4. DEGs功能注释及富集分析
首先基于筛选出的DEGs，通过python脚本从红鳍东方鲀的基因组注释文件中获得了DEGs的蛋白序列及相关注释信息。然后采用KEGGannotation v1.0将获得的DEGs的蛋白序列比对到人类蛋白数据库，得到与DEGs的SYMBOL相对应的ENTREZID。最后使用R语言包ClusterProfilerv3.18.1进行基因本体论(Gene ontology, GO)和京都基因与基因组百科全书(Kyoto encyclopedia of genes and genomes, KEGG)富集分析。
3. 结果与分析
3.1. 转录组测序及比对结果
本实验对6个文库进行了测序并通过生物信息学手段去除了包括接头序列等在内的低质量reads，得到了足够数量的clean reads，所有下游分析都是基于高质量的clean reads进行的。我们将暗纹东方鲀和红鳍东方鲀鳔的clean reads分别与红鳍东方鲀的参考基因组进行比对，如表1所示，可以看到所有样本的测序序列与参考序列匹配以后，至少有超过89.7%的clean reads能够很好的比对到红鳍东方鲀的参考基因组上。基于SAMtool v1.7对比对结果进行排序去重后的分析见表2，所有样本中比对上的reads的重复率约为26.00%~30.60%，重复率较低，去重后也至少得到了38.0Mreads供后续分析。这些结果均表明我们的测序结果质量良好，适合用于下一步分析。所有清洁数据都已提交至NCBI (Bioproject: PRJNA910922, PRJNA645780)。
Table 1. Sequencing alignment statistics
表1. 测序数据比对率统计
Table 2. Summary of sequencing quality control results after comparison to the reference genome
表2. 比对到参考基因组后的测序质控结果汇总
3.2. 差异表达基因分析
我们通过HTSeq v0.11.2统计了read counts，然后采用R语言包DESeq2 v3.10进行了差异基因表达分析，以TOSB vs TRSB和p.adj < 0.05 & |log2foldchange| ≥ 1为标准筛选得到了DEGs。如图1所示，可以看出在两种东方鲀之间共筛选得到了4062个DEGs，其中上调表达基因2282个，下调表达基因1780个。为了更好的探索DEGs的表达模式，经HTSeq v0.11.2输出的原始readcounts被均一化，后续的分析也都是以均一化后的counts进行的。在此基础上，我们通过R语言包和python脚本分别进行了差异基因聚类分析和主成分分析。如图2(a)所示，在差异基因聚类热图中红色代表高表达量的基因，绿色代表低表达量的基因，根据表达量可以很明显的区分出TOSB与TRSB两组之间差异表达的基因，同时我们也可以看到组内基因的表达量差异很小，这说明同一处理组之间的生物重复性良好。在图2(b)主成分分析三维图中我们可以看到不同处理组的样本能够被很好的聚类到了两个区，这说明TOSB与TRSB组间差异显著，这也证明了我们实验的合理性。
Figure 1. Volcano plot of differentially expressed genes. The green dots represent down regulated genes, the red dots represent up regulated genes, and the black dots represent genes that were not significantly differentially expressed
图1. 差异表达基因火山图。绿色点代表下调的基因，红色点代表上调的基因，黑色点代表没有明显差异表达的基因
3.3. GO功能富集分析
如图3所示，我们以pcutoffer = 0.05作为截取值，对TOSB vs TRSB中筛选出的DEGs进行了GO富集分析。结果表明：4062个DEGs共富集到了1124条术语上，其中包含888条生物过程(Biological Processes, BP)术语，135条细胞组分(Cell Composition, CC)术语和101条分子功能(Molecule Function, MF)术语。在生物过程术语中，突触组织(synapse organization)、化学突触传递的调节(modulation of chemical synaptic transmission)和解剖结构大小的调节(regulation of anatomical structure size)是富集显著的前3条术语，分别注释到了111、110和121个基因；胶原细胞外基质(collagen-containing extracellular matrix)、突触膜(synaptic membrane)和突触膜的整体组成(intrinsic component of synaptic membrane)是细胞组分中富
Figure 2. (a): Heat map of differentially expressed gene clustering; (b): Principal component analysis graph
图2. (a) 差异表达基因聚类热图；(b) 主成分分析图
Figure 3. GO annotation enrichment map of differentially expressed genes. (a): biological processes; (b): cellular components; (c): molecular functions.
图3. 差异表达基因GO注释富集图。(a)：生物过程；(b)：细胞组分；(c)：分子功能
集显著的前3条术语，分别注释到了127、114和59个基因；而在分子功能中，糖胺聚糖结合(glycosaminoglycan binding)、硫化物结合(sulfur compound binding)和肝素结合(heparin binding)则是富集最显著的术语，分别注释到了77、80和57个基因。由上述结果我们可以发现，鳔中含有大量的神经轴突，突触与突触间信号传递的差异性可能是介导两种东方鲀的鳔发挥不同功能的重要途径。
除此之外，由图4我们还发现上述显著富集的GO术语可以被很好地聚类为七个大类，其中包括“adcnylatc body chemical assembly”、“alcohol import acid into”、“heart contraction circulation pressure”、“tube anatomical diameter size”及与信号传导有关的“calcium cytosolic ion concentration”、与听力有关的“amide ear hormone differentiation”、与分支形态有关的“morphogenesis axon branching oxygen”。由富集结果可以看出，两种东方鲀在信号转导及听力功能方面存在显著差异。
Figure 4. Functional grouping network diagram for GO enrichment analysis
图4. GO富集分析功能聚类网络图
3.4. KEGG通路富集分析
将TOSB vs TRSB中筛选出的DEGs注释到KEGG上，以便系统分析这些DEGs的关联性。KEGG富集结果表明：筛选出的4062个DEGs共富集到了339条信号通路上，其中有22条信号通路显著富集(p.adj < 0.05)。我们对22条显著富集的信号通路按大类进行了条理化分析，结果见图5，可以看到大部
Figure 5. Classification plot of 22 significantly enriched (p.adj < 0.05) signal pathways
图5. 22条显著富集(p.adj < 0.05)信号通路归类图
分归属于有机体系统(Organismal Systems)和环境信息处理(Environmental Information Processing)，这可能与两种东方鲀生存环境的差异性密切相关。此外如图6所示，我们按照p.adj值由低到高选择了前50条信号通路进行绘图和分析，结果显示神经活性配体-受体相互作用途径(Neuroactive ligand-receptor interaction)、钙信号传导途径(Calcium signaling pathway)、谷氨酸能突触(Glutamatergic synapse)、昼夜节律(Circadian entrainment)和cAMP信号传导途径(cAMP signaling pathway)是DEGs显著富集的前5条信号通路，分别注释到了104、71、36、32和58个基因。
Figure 6. The top 50 enriched KEGG pathways
图6. 前50条KEGG通路图
4. 讨论
鳔作为硬骨鱼特有的器官，其功能的多样性与鳔本身结构形态的差异性密切相关。我们通过阅读文献发现，通常按照有无与咽部相连的鳔管将鳔分为了喉鳔类和闭鳔类，不同类型的鳔扮演不同的角色，发挥不同的功能，甚至不同鱼类鳔腔的数量也有所不同，每个腔室亦有不同的功能 [5] 。同时我们也发现与河鲀鳔组织结构有关的研究少之又少，我们无法得到太多与暗纹东方鲀或红鳍东方鲀鳔组织结构相关的信息，只能通过组织取样时进行观察。我们发现，两种东方鲀的鳔组织只有一个腔室，整体上呈心形，内部充满气体，可以完整的剥离出来且不漏气，因此我们认为它属于闭鳔类。尽管如此，暗纹东方鲀与红鳍东方鲀亲缘关系密切但生存环境却差异巨大的特点仍然引起了我们极大的兴趣，它们也是养殖河鲀鱼中出口量占比最大的两个品种，因此我们的研究是必要的。
4.1. 基于神经活性配体受体相互作用的信号传导途径
鳔含有丰富的神经系统，不同鱼类鳔的神经类型和神经支配方式既不相同又存在进化上的保守性 [16] 。神经活性配体受体相互作用信号通路与神经功能密切相关，是质膜上所有与细胞内外信号通路相关的配体和受体的集合。其中胞外信号(配体)通常包括激素和神经递质，膜受体则以G蛋白偶联受体(G Protein-Coupled Receptors, GPCRs)为主，通常被分为A类视紫红质胺(rhodopsin-like amine)、B类分泌素样(secretin-like)、C类亲代谢型谷氨酸盐/外激素(metabotropic glutamate/pheromone)和D类离子通道/其他受体(channels/other receptors) [17] [18] 。当膜受体被激活后，以GPCRs为主的下游转导器会激活质膜上的腺苷酸环化酶(Adenylate Cyclase, AC)，将ATP转化为第二信使环磷酸腺苷(Cyclic adenosine monophosphate,cAMP)，而不同类型的胞外信号和膜受体的结合决定着cAMP通路被激活亦或是被抑制 [19] [20] [21] 。研究发现，属于A类的神经肽(NPY)作为食欲调节激素在硬骨鱼中广泛表达，其还会通过与GPCRs的结合来抑制AC的磷酸化 [22] 。同时A类的多巴胺受体D2(DRD2)也属于抑制型GPCRs，研究表明其是斑马鱼眼睛正常发挥功能所必需的 [23] 。有趣的是，多巴胺受体的另外两个亚型DRD1和DRD5却属于激活型GPCRs，多巴胺与之结合会激活AC的磷酸化。受体属于D类的催乳素(PRL)被证明是罗非鱼盐度适应中一种主要的渗透调节激素 [24] ，敲除PRL会造成斑马鱼鳔的发育性缺陷 [25] 。综上所述，种类丰富的GPCRs会通过与不同配体的结合参与调控多种信号的传导，进而发挥不同的功能。在我们的研究中，相当数量的DEGs富集在了此通路上，这说明两种东方鲀与胞内外信号传导相关的配体和受体存在明显差异性，而作为胞外信号向胞内传递的第一个环节，它的非特异性亦会导致后续信号传导的多变性，进而导致两种东方鲀的鳔在功能上的差异性。
4.2. 鳔与听觉的关系
听觉是水生生物生存和交流的一种重要方式，鱼类已经进化出了相应的听觉结构用于提高听觉的灵敏性和频率的广泛度，包括鱼鳔、内耳和韦伯氏器等 [26] 。值得注意的是，并不是所有鱼类都具有韦伯氏器，缺少韦伯氏器的鱼类通常由鳔的前方生出的盲囊与内耳连接以接收由鳔传来的振动，而不同鱼类对听觉的灵敏度各不相同，通常由是否存在听力相关结构或听力相关结构与内耳的耦合度来决定 [27] [28] 。研究发现GABA能突触(GABAergic synapse)和谷氨酸能突触(Glutamatergic synapse)在听觉系统中起重要作用，轴突引导(Axon guidance)有助于听觉的发展，cAMP信号通路(cAMP signaling pathway)则被证实是听觉受体上皮细胞增殖所必须的 [29] [30] [31] 。在我们的研究中，上述通路得到了显著富集，且在GABA能突触和谷氨酸能突触中大部分DEGs上调表达，在轴突引导和cAMP信号通路中大部分DEGs下调表达。与听觉功能相关的基因在我们的研究中也被鉴定为DEGs，如：突变后会引起听力丧失的重要基因coch [32] 、通过引起内耳基底膜毛细胞突触损伤进而介导小鼠老年性聋的kcnq2 [33] 、与听觉神经元的兴奋性突触反应密切相关的各种GABA受体基因和谷氨酸受体基因(gabrb2, gabrb3, LOC101068952, LOC101075576, LOC105416167, grin2a, grik4, gria4, LOC101064905) [34] [35] [36] 。与此同时，Yang, Wang等 [37] 通过对斑点叉尾鮰鳔和其它四种组织的转录组比较分析，筛选出了一系列与听觉相关的基因和信号通路。有趣的是我们基于两种亲缘关系密切但生存环境差异巨大的东方鲀鳔组织同样能够筛选出来一系列与听觉相关的差异基因和信号通路。这提示我们，两种东方鲀不但存在听觉功能，且听觉的灵敏度可能存在巨大差异，这可能与它们生存环境的差异性有关，或许是受到了适应性进化机制的调控。
5. 结论
本文通过对TOSB vs TRSB的转录组比较分析，以p.adj < 0.05 & |log2foldchange| ≥ 1为阈值筛选得到了4062个DEGs，其中上调表达基因2282个，下调表达基因1780个。通过GO和KEGG对DEGs进行了富集分析，富集结果显示暗纹东方鲀和红鳍东方鲀的鳔在神经活性配体受体信号传导方面存在显著差异，不同类型的配体和受体的相互作用往往介导着不同信号的传递，进而发挥不同的功能。除此之外，我们还发现了大量与听觉功能相关的基因及通路得到了显著富集。以往的研究通常是通过比对同种鱼类间鳔与其他组织的差异得出其存在听觉功能的结论，而我们通过比较两种亲缘关系密切但生存环境差异巨大的东方鲀的鳔同样发现其存在听觉功能，且听觉的灵敏度有所差异。
综上所述，可以发现生存环境的差异性决定了亲缘关系密切的暗纹东方鲀和红鳍东方鲀鳔的功能的不同，主要体现在神经活性配体受体信号传导及听觉功能方面，我们推测这可能是适应性进化的结果。
致谢
特别鸣谢大连海洋大学水产设施养殖与装备工程团队(Aquacultural Engineering R&D Team, AET)提供的养殖设备及驯养场地，也特别感谢周金旭和徐浩在取样时提供的帮助。
基金项目
本论文由国家海水鱼产业技术体系项目(CARS-47)、大连市重点领域创新团队项目(2021RT07)及大连海洋大学科技创新团队项目(B202102)资助。
参考文献
NOTES
*共同第一作者。
#通讯作者。