斑马鱼软骨内骨化过程中骨骼的转录组测序分析

斑马鱼(Danio rerio)是研究骨骼发育与疾病的良好模型[1]。以往骨骼疾病研究通常是在小鼠中进行，但小鼠繁育周期长且产仔量少，加上动物试验3R原则的出现，促使科研工作者的视线聚焦到了斑马鱼上[1]。斑马鱼骨骼与哺乳动物骨骼细胞类型相同，包含骨细胞、成骨细胞和破骨细胞3种细胞类群；成骨方式相似，具有膜内骨化和软骨内骨化两种成骨方式；骨骼稳态的保持也相似，依赖骨吸收和骨形成的协同调控[2]。斑马鱼骨骼发育的分子调控机制也非常保守，哺乳动物中骨骼发育的关键调控因子runx2、osx和osn等在斑马鱼中高度同源且功能相似[3]。哺乳动物中已知的与骨骼发育紧密相关的信号通路如Wnt/β-catenin、Hedgehog、TGF-β等在斑马鱼中也发挥同样作用[4-6]。研究者已构建了多种斑马鱼骨骼疾病模型，如斑马鱼骨质疏松症模型，也为研究骨骼疾病提供了新的思路[7]。

骨骼不仅能够给动物机体提供支撑，与骨骼肌一起协助机体完成复杂的运功，还具有强大的内分泌功能。骨骼能分泌多种生物活性因子，以旁分泌或自分泌的方式调控骨骼的发育和代谢过程[8]，其分泌的多种生物活性物质还可经循环系统参与调节机体内多种组织器官的代谢和功能。如成骨细胞分泌的骨钙素可直接作用于胰岛β细胞，调节其糖代谢过程[9]，也能增加脂质的利用从而减少脂肪的堆积。此外，成骨细胞还是雌激素的主要靶标，雌激素水平降低会间接激活破骨细胞分化，使破骨细胞主导的骨吸收大于成骨细胞主导的骨形成过程，最终导致骨质疏松[10]。骨骼系统的稳态受到内分泌和旁分泌的协同调控，一些内分泌调控因子，如甲状腺激素、降钙素、性激素等，以及一些旁分泌调控因子如前列腺素和生长因子等的异常都会破坏骨重建稳态，从而诱发骨骼疾病[8]。由于骨骼疾病受机体复杂调控网络的影响，骨骼疾病的治疗始终是个难题，再加上现有研究大多致力于如何降低骨吸收过程，缺少对骨形成过程的利用，因此，对活体骨骼进行详尽研究就显得至关重要[11]。

目前，以斑马鱼作为骨骼疾病研究模型的相关研究较多，但一些关键时期骨骼发育模式和调控机制及调控网络还有待深入探究。斑马鱼骨骼从受精后3 d(day post fertilization，dpf)开始发育，在30 dpf左右时脊椎骨骼基本形成，70 dpf左右时所有头骨骨骼基本形成，一直到3～4个月时斑马鱼骨骼发育完全成熟。在30～70 dpf时，斑马鱼的骨骼发育正在经历软骨内骨化的过程，而这一过程的分子调控机制尚不清楚。在NCBI数据库中已有在30～45 dpf对斑马鱼进行转录组测序的结果。在45～70 dpf时，骨骼仍在经历软骨内骨化过程，但缺乏转录组测序分析的报道。本研究中，选择斑马鱼两个时间点(50 dpf和65 dpf)的骨骼进行转录组测序，以期挖掘更多骨骼发育的关键基因，丰富斑马鱼软骨内骨化过程中骨发育的分子调控网络，为骨骼疾病的治疗提供新思路。

1 材料与方法

1.1 材料

试验用斑马鱼为AB品系，购于中国科学院上海生命科学研究院生物化学与细胞生物学研究所，并根据上海海洋大学动物伦理委员会的要求进行处理，按照斑马鱼养殖标准饲养。试验用鱼均等密度养殖于循环养殖系统(Aquaneering Inc，San Diego，CA，USA)，水温控制在27～28 ℃，盐度为0.25～0.75，pH为7.0～8.0，氨氮及硝态氮质量浓度不超过200 mg/L，光照条件处理为14 h照明交替10 h黑暗。

1.2 方法

1.2.1 阿尔辛蓝-茜素红双重染色取50 dpf和65 dpf斑马鱼样品，用体积分数为4%的多聚甲醛(PFA)在4 ℃冰箱固定48 h。将PFA吸出，样品用1×PBST清洗两次，每次5 min,在1%(质量分数)KOH和5%(质量分数)H2O2的混合液(二者体积比为1∶1)中,自然光下漂白2 h，直到鱼眼呈淡黄色。加入体积分数为0.05%的阿尔辛蓝避光染色24 h，移去染色液，用体积分数为100%、70%、50%、30%的乙醇和1×PBST清洗样品，每次清洗30 min。用质量分数为0.1%的胰蛋白酶消化过夜，移去消化液，用1×PBST清洗两次，每次30 min。再用0.5%的KOH配制成0.5%的茜素红溶液，染色 1 h，待硬骨染成红紫色时，用体积分数为25%、50%、75%的甘油保存，待样品透明后，再置换于体积分数为100%的甘油溶液中保存。

1.2.2 斑马鱼骨组织解剖与获取从斑马鱼全身骨骼中获取测序使用的骨组织，包括头骨、脊椎骨和鱼鳍等。具体如下：分别取50 dpf和65 dpf的斑马鱼，使用1×MS222进行麻醉，将其固定在泡沫板上剥去鱼皮，用手术刀和镊子切除内脏，再将头部和身体分离，置于解剖镜下将肌肉组织剔除，再将取出的斑马鱼骨骼组织置于质量分数为0.1%的胰蛋白酶溶液中，37 ℃下水浴消化残存肌肉，并于解剖镜下确认肌肉消化程度，以确保获取组织仅为斑马鱼骨组织。反应结束后，使用1×PBS缓冲液反复润洗骨组织。将采集的骨组织样本立即放入液氮中速冻，于-80 ℃超低温冰箱中保存待用。

1.2.3 骨组织总RNA提取及转录组测序取适量的骨组织样品，用试剂盒(mirVanaTM miRNA ISOlation Kit，Ambion-1561)进行RNA提取，用NanoDrop 2000检测RNA浓度，用Agilent 2100 Bioanalyzer生物分析仪检测RNA质量(RIN≥7且28S/18S≥0.7)，样品总RNA使用DNase消化DNA后，再用带有Oligo(dT)的磁珠富集mRNA，加入打断试剂将mRNA打断成短片段，用打断后的mRNA为模板合成一链cDNA，然后继续合成二链cDNA，并使用试剂盒纯化双链cDNA，纯化的双链cDNA进行末端修复再加A尾并连接测序接头，然后进行片段大小选择，最后进行PCR扩增。将构建好的文库经Agilent 2100 Bioanalyzer质检合格后，使用Illumina HiSeq X Ten测序平台进行测序，产生150 bp的双端数据。

1.2.4 转录组数据分析通过Illumina平台得到大量的样本双端测序数据，测序中产生大量原始数据(raw reads)。为获得可用于后续分析的高质量reads，需要对raw reads进行过滤。首先使用Trimmomatic软件进行质量检测并去除接头，过滤掉低质量碱基，获得高质量的clean reads[12]。再使用Hisat 2软件将clean reads比对到斑马鱼的参考基因组，并通过基因组比对率来评估样本的情况[13]。使用Cufflinks软件定量基因FPKM表达量值[14]，在计算基因的表达量差异时，通过Htseqcount软件获取落到各个样本中基因的reads数目[15]，利用DESeq(2012)R Package对数据进行标准化，并用NbinomTest函数计算差异比较的P值和fold change值，使用DESeq(2012)R软件进行差异表达分析[16]，将|log2(fold change)|≥1且P<0.05设置为显著差异表达的阈值。对差异表达基因(DEG)进行分层聚类分析，以证明不同组样本中基因的表达模式。基于超几何分布分别对DEGs进行GO(gene ontology)富集和KEGG(kyoto encyclopedia of genes and genomes)途径富集分析，以分析差异基因主要影响的生物学功能或通路[17]。

1.2.5 荧光定量PCR(qRT-PCR) 分别取50 dpf和65 dpf的斑马鱼各3尾提取骨骼组织总RNA。使用试剂盒(TaKaRa，RR047A)将RNA(1 μg)反转录成cDNA，并保存于-80 ℃超低温冰箱中。qRT-PCR所用引物在NCBI网站设计(表1)，所有引物扩增效率均为90%～105%。以cDNA为模板、β-actin为内参进行荧光定量检测。反应条件：95 ℃下变性10 s，60 ℃下退火30 s，72 ℃下延伸1 min，共进行40个循环。荧光定量数据采用2-ΔΔCt法进行计算，最后使用Prism 6.0软件作图。

2 结果与分析

2.1 50 dpf和65 dpf斑马鱼软骨内骨化进程分析

利用阿尔辛蓝-茜素红双重染色(阿尔辛蓝染色标记软骨结构，茜素红染色标记硬骨结构)，对50 dpf和65 dpf两个时期的斑马鱼骨骼进行染色。结果发现，50 dpf斑马鱼头骨基本上全由软骨构成，而65 dpf斑马鱼头骨已经发育分化为硬骨结构，仅剩少部分骨骼仍处在软骨内骨化进程中(图1)。这表明，与50 dpf相比，65 dpf的斑马鱼骨骼发育已较为成熟，大部分骨骼已经发育完善，但仍有部分骨骼正在经历软骨内骨化阶段。

2.2 转录组序列质量分析

为挖掘两个骨骼发育时期调控通路的变化，对50 dpf和65 dpf野生型斑马鱼骨骼共6个样本进行转录组测序分析，共获得41.44×109个clean bases，各样本的有效数据量为(6.60×109)～(7.17×109)，Q30碱基分布为93.46%～93.83%，平均GC含量为46.63%(表2)。说明测序质量较好，均可用于后续分析。

2.3 差异表达基因

基因表达数据相关性分析表明，65 dpf和50 dpf野生型斑马鱼，同一组样本间呈现高相关性(图2(a))。比较分析后发现，65 dpf斑马鱼相比50 dpf斑马鱼共有610个差异表达基因，其中，237个基因表达显著上调，373个基因表达显著下调(图2(b)、(c))。

表达差异最显著的14个基因全部是显著上调的基因，即β晶体蛋白基因(crybb1l3)、γ晶体蛋白基因(crygm2d1)、麦芽糖酶-葡萄糖淀粉酶基因(si:ch73-111m19.2)、E3泛素蛋白连接酶基因(LOC558300)、GTPase IMAP家族成员(si:ch1073-365p7.2)、未知蛋白的编码基因(LOC- 110438923、zgc:173545、zgc:173544、LOC11043- 8920、LOC100002043、si:ch211-146l10.8)、晶状体纤维b的主要内在蛋白基因(mipb)、肽酶抑制剂分支a基因(serpina10a)和透明带精子结合蛋白2.2基因(zp2.2)(表3)。

2.4 差异表达基因的GO富集分析

对差异表达基因进行GO富集分析发现，本次共有848个unigene获得注释，主要包括生物学过程(biologolical process)、分子功能(molecular function)和细胞组分(cellular compent)3大类。其中，注释到生物学过程的基因最多，涉及479个基因，其次有271个基因注释到分子功能，注释到细胞组分的基因有100个。

GO富集显示：差异最显著的生物学过程主要涉及代谢和生殖调节等过程，如高密度脂蛋白颗粒重塑、精子与透明带的结合、脂质运输、卵外膜形成、胆固醇逆转运、顶体反应的正向调节、脂蛋白的代谢过程、鳍再生、胆酰酯化的正向调控、磷脂外排、极低密度脂蛋白颗粒重塑、甘油三酯分解过程和胆固醇平衡等；在分子功能方面，主要涉及胆固醇结合、胆固醇运输活动、转录抑制因子活性、碳水化合物结合、脂质结合、GTP结合和细胞外基质结构成分等；在细胞组分层面，主要涉及基底膜中间纤维、细胞连接、胶原蛋白三聚物、基底膜和胞外区细胞外空隙等(图3)。这些差异说明，相较50 dpf斑马鱼，65 dpf斑马鱼骨骼组织的脂质利用与转运过程变得较为活跃，脂质代谢途径占据较大部分比例，这表明斑马鱼骨骼的成熟与脂质代谢活动存在密切关系。

2.5 差异表达基因的KEGG通路富集分析

对检测到的所有的差异表达基因进行KEGG分析，共富集到89个Pathway，这些Pathway主要涉及代谢途径、遗传信息处理、生物系统、细胞过程、环境信息处理和人类疾病6大类。其中，代谢途径占据较大部分比例，占总富集通路数的42.7%，其次是生物系统和环境信息处理，人类疾病占比最少，只占总富集通路数的4%。这表明，差异表达基因主要与代谢通路、蛋白质的加工合成和细胞间相互作用的信号通路相关(图4)，主要包括甘油磷脂代谢、细胞衰老、黏着斑、内质网蛋白质加工、脂肪细胞因子信号通路、FoXO信号通路、亚油酸代谢、叶酸生物合成、α-亚麻酸代谢、细胞因子及其受体相互作用、精氨酸与脯氨酸代谢。其中，富集程度最高的是FoXO信号通路，FoXO属于转录调控因子家族，广泛表达于正在发育或成熟的组织器官中，FoXO1家族成员还是成骨细胞增殖分化的关键转录因子。

2.6 荧光定量PCR验证

随机选取6个差异表达基因进行qRT-PCR验证，结果显示，相较于50 dpf斑马鱼，65 dpf斑马鱼的ucmab和comp基因表达极显著下调(P<0.01)，而sox9b、ak1、zp2.2和si:ch1073-365p7.2等基因表达显著上调(P<0.05)(图5)，其变化趋势与转录组测序分析结果一致，证明转录组测序结果真实可靠。

3 讨论

3.1 软骨内骨化过程骨骼发育相关基因及通路分析

与哺乳动物一样,斑马鱼骨骼包括软骨细胞、成骨细胞和破骨细胞3种细胞类型，且同样具有膜内骨化和软骨内骨化两种成骨方式，故被当作良好的骨骼研究模型。然而,目前有关斑马鱼软骨内骨化过程分子调控机制的研究相对较少。本研究中，利用Illumina转录组测序技术对不同时期的斑马鱼骨骼进行研究，相较50 dpf斑马鱼，65 dpf斑马鱼成骨细胞标记基因(runx2、sp7)的表达无显著性差异，但一些关于生长板上部区域和软骨基质相关基因(comp)[18]、软骨发育相关基因(dlx2a、foxd3、sec23b)[19]、碱性磷酸酶基因(alpi.2)、骨矿化调节基因(ptk2bb)、骨骼细胞分化正调控基因(dlx2a)、颅骨系统发育基因(cyldb、gpatch3)[20]和软骨发育与软骨内骨形态发生基因(foxd3)[21]等表达量均下调，还有一些调控骨骼生长分化及修复的信号通路，如Notch信号通路、BMP信号通路和经典的Wnt信号通路等也下调。结合两个时间节点的阿尔辛蓝-茜素红双重染色结果，说明在这两个时间点斑马鱼硬骨骨骼结构发育已较为成熟，骨细胞分化进程平稳，但调控软骨形成和软骨内骨化的相关基因表达均下调，表明软骨形成和软骨内骨化进程减缓，65 dpf左右软骨内骨化接近尾声。

3.2 斑马鱼骨骼发育过程代谢通路分析

研究表明，哺乳动物的骨骼发育受到代谢通路的调控，如载脂蛋白E基因(ApoE)是间充质祖细胞骨形态发生蛋白的调控基因，在小鼠胚胎发育过程中高表达，且ApoE还在小鼠成骨细胞发育后期表达[22]。进一步的研究还发现，细胞内胆固醇和脂质利用过程可以调控软骨细胞分化和动态平衡[23]。然而，目前还不了解斑马鱼骨骼发育是否受到代谢相关基因及通路的调控。本研究中，相比50 dpf斑马鱼，65 dpf斑马鱼糖代谢和脂质代谢途径均发生了显著性的改变，表明骨骼和代谢过程间存在复杂的相互调控机制。其中，调控脂质分解及合成的基因，如磷脂酶C基因(plcd1a)、胆固醇25-羟化酶样1基因(ch25hl1.1)、磷脂转移蛋白基因(pltp)、烷基甘油单加氧酶基因(agmo)、载脂蛋白Ea基因(apoea)、脂肪酸结合蛋白基因(fabp1b.1)和载脂蛋白A-IV b基因(apoa4b.1)等，这些差异表达基因主要富集到类固醇生物合成途径、脂肪细胞因子信号通路、甘油酯代谢途径、乙醚脂质代谢途径和亚油酸代谢途径；调控碳水化合物及糖代谢的基因，如胰岛素受体底物基因(irs2b)、葡聚糖谷氨酸连接酶基因(glula)、甲壳素酶基因(chia.1)、磷酸烯醇丙酮酸激酶基因(pck1、pck2)、乳糖酶基因(lct)、血糖皮质激素调节酶基因(sgk1、sgk2b)和凝集素-半乳糖苷结合基因(lgals2b)等，这些基因主要富集到胰岛素信号通路、糖酵解途径、磷酸戊糖途径、果糖和甘露糖代谢途径。这些结果表明，斑马鱼骨骼细胞分化过程同样需要大量的糖代谢和脂肪分解过程来提供能量支撑，同时骨骼分泌的一些骨源性因子也同样调控着机体的代谢途径[24]。

3.3 斑马鱼骨骼发育与生殖功能相关通路分析

哺乳动物中生殖功能的损伤会影响骨骼重建过程，并且成骨细胞分泌的活性因子也会调控机体的生殖功能。本研究中，对50 dpf和65 dpf两个时期斑马鱼骨组织的差异表达基因进行GO功能注释和KEGG富集分析发现，除了与机体生长发育相关的基因差异显著外，与生殖过程相关的基因也产生显著性差异，且呈上调趋势，如凝集素-甘露糖结合基因(lman1)、雌酮磺基转移酶基因(sult1st2、sult2st3)、透明带精子结合蛋白2.2基因(zp2.2)[25]、卵母细胞调控因子基因(buc、star)和苯丙氨酸羟化酶基因(pah)等，这些基因会刺激卵母细胞成熟，并可作用于雌激素代谢过程，还可以加强细胞对雌激素的刺激反应。有研究表明，雌激素介导的途径是在骨骼和软骨形成(包括颅面骨骼发生)中起作用的分子信号之一，而凝集素-甘露糖结合基因(lman1)可能对雌激素信号做出反应，并在颅面发育中发挥作用[26]。由此推测，骨骼与生殖功能之间可能存在紧密的相互作用关系。

3.4 斑马鱼骨骼、代谢和生殖功能三者间的相互调控

哺乳动物的新陈代谢、骨骼和生殖功能三者间存在紧密联系，受到某些信号通路的共同调控。如脂肪细胞分泌的瘦素分子具有广泛的生物学功能，瘦素可以调节身体代谢，提高脂肪利用率，避免肥胖[27]。哺乳动物缺乏瘦素，会导致严重肥胖，并且不育，表明瘦素不仅可以调节身体的能量利用过程，维持身体的能量平衡，而且可以通过下丘脑-垂体-性腺轴调节机体的生殖功能[28]。另外，瘦素受体存在于所有卵泡细胞中，可与促黄体生成素一起刺激卵泡膜细胞和卵母细胞的成熟[28]。对瘦素和瘦素受体敲除小鼠模型研究表明，瘦素能够影响骨骼发育，缺乏瘦素或瘦素受体会导致骨量增加，而注射瘦素会导致骨质流失[29]。综上所述，瘦素可作为一种重要的内分泌调节分子，调节新陈代谢、生殖和骨量。然而，目前尚不清楚这种调控模式是否存在于鱼类中。本研究中，结合GO富集和KEGG通路富集分析发现，大部分差异表达基因主要富集到脂质、糖代谢和生殖功能调节等途径，如类固醇生物合成途径、脂肪细胞因子信号通路、甘油酯代谢途径、乙醚脂质代谢途径、胆固醇逆转运、胆固醇平衡、脂质运输、脂蛋白代谢过程、胆酰酯化的正向调控、甘露糖代谢途径、糖酵解途径、胰岛素信号通路、糖异生、对葡萄糖的反应、卵外膜形成、卵母细胞成熟途径和顶体反应的正向调节等。在斑马鱼50～65 dpf的骨骼发育进程中，随着斑马鱼骨骼发育的成熟，除调控骨骼生长发育的相关基因及通路发生改变外，一些糖代谢、脂质合成与转运及生殖功能调节过程等相关途径均产生了显著差异。综合转录组测序结果及早期对哺乳动物的研究，作者推测，斑马鱼骨骼发育可能与哺乳动物一样，与代谢网络和生殖过程存在紧密联系，且三者间存在相互调控的关系。

4 结论

1)通过对50 dpf和65 dpf斑马鱼进行阿尔辛蓝-茜素红双重染色，发现50 dpf斑马鱼头骨正在经历软骨内骨化，而65 dpf斑马鱼头骨软骨内骨化过程基本完成。

2)50 dpf和65 dpf斑马鱼骨组织的差异基因主要富集在脂质代谢、糖代谢及生殖过程相关的多个通路。

3)本研究中转录组测序数据首次揭示了鱼类的代谢、骨骼和生殖功能三者间可能存在相互调控关系，这为研究骨骼疾病提供了一个新的切入点，为进一步研究鱼类骨骼发育，阐明骨骼发育过程中的调控模式，建立复杂的疾病模型提供了新思路。

[1] WACHSMUTH L,MENSEN A,BARCA C,et al.Contribution of preclinical MRI to responsible animal research:living up to the 3R principle[J].Magnetic Resonance Materials in Physics,Biology and Medicine,2021,34(4):469-474.

[2] HUANG W,YANG S Y,SHAO J Z,et al.Signaling and transcriptional regulation in osteoblast commitment and differentiation[J].Frontiers in Bioscience,2007,12:3068-3092.

[3] LLERAS-FORERO L,WINKLER C,SCHULTE-MERKER S.Zebrafish and medaka as models for biomedical research of bone diseases[J].Developmental Biology,2020,457(2):191-205.

[4] BENSIMON-BRITO A,RAMKUMAR S,BOEZIO G L M,et al.TGF-β signaling promotes tissue formation during cardiac valve regeneration in adult zebrafish[J].Developmental Cell,2020,52(1):9-20.

[5] SCHWEND T,AHLGREN S C.Zebrafish con/disp1 reveals multiple spatiotemporal requirements for Hedgehog-signaling in craniofacial development[J].BMC Developmental Biology,2009,9:59.

[6] TANG D M,HE Y Z,LI W Y,et al.Wnt/β-catenin interacts with the FGF pathway to promote proliferation and regenerative cell proliferation in the zebrafish lateral line neuromast[J].Experimental & Molecular Medicine,2019,51(5):1-16.

[7] CARNOVALI M,BANFI G,MARIOTTI M.Zebrafish models of human skeletal disorders:embryo and adult swimming together[J].BioMed Research International,2019(20):1253710.

[8] WANG H,ZHENG X Q,ZHANG Y,et al.The endocrine role of bone:novel functions of bone-derived cytokines[J].Biochemical Pharmacology,2021,183:114308.

[9] DOGA M,BONADONNA S,GIUSTINA A.Glucocorticoids and bone:cellular,metabolic and endocrine effects[J].Hormones,2004,3(3):184-190.

[10] KHOSLA S,MELTON L J,RIGGS B L.Estrogens and bone health in men[J].Calcified Tissue International,2001,69(4):189-192.

[11] LIN C X,CHEN Z,GUO D,et al.Increased expression of osteopontin in subchondral bone promotes bone turnover and remodeling,and accelerates the progression of OA in a mouse model[J].Aging,2022,14(1):253-271.

[12] BOLGER A M,LOHSE M,USADEL B.Trimmomatic:a flexible trimmer for Illumina sequence data[J].Bioinformatics,2014,30(15):2114-2120.

[13] KIM D,LANGMEAD B,SALZBERG S L.HISAT:a fast spliced aligner with low memory requirements[J].Nature Methods,2015,12(4):357-360.

[14] 魏亚丽,周艳,黄思婕,等.高温胁迫下尼罗罗非鱼肝脏组织的转录组分析[J].大连海洋大学学报,2021,36(2):222-228.

WEI Y L,ZHOU Y,HUANG S J,et al.Transcriptome analysis of liver tissue of Nile tilapia Oreochromis niloticus exposed to high temperature stress[J].Journal of Dalian Ocean University,2021,36(2):222-228.(in Chinese)

[15] TRAPNELL C,WILLIAMS B A,PERTEA G,et al.Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation[J].Nature Biotechnology,2010,28(5):511-515.

[16] LEE J,JI Y,LIANG S D,et al.Bayesian hierarchical model for differential gene expression using RNA-seq data[J].Statistics in Biosciences,2015,7(1):48-67.

[17] KANEHISA M,ARAKI M,GOTO S,et al.KEGG for linking genomes to life and the environment[J].Nucleic Acids Research,2008,36:D480-D484.

[18] 白伦浩，王瑜，张田鹏.骨性关节炎生物学标记物软骨寡聚基质蛋白的研究进展[J].中国修复重建外科杂志, 2010,24(6):714-717.

BAI L H,WANG Y,ZHANG T P.Research progress of cartilage oligomeric matrix protein[J].Chinese Journal of Reparative and Reconstructive Surgery,2010,24(6):714-717.(in Chinese)

[19] SPERBER S M,SAXENA V,HATCH G,et al.Zebrafish dlx2a contributes to hindbrain neural crest survival,is necessary for differentiation of sensory ganglia and functions with dlx1a in maturation of the arch cartilage elements[J].Developmental Biology,2008,314(1):59-70.

[20] FERRE-FERN width=7,height=9,dpi=110

NDEZ J J,AROCA-AGUILAR J D,MEDINA-TRILLO C,et al.Whole-exome sequencing of congenital glaucoma patients reveals hypermorphic variants in GPATCH3,a new gene involved in ocular and craniofacial development[J].Scientific Reports,2017,7:46175.

[21] KAMALDINOV T,ERNDT-MARINO J,DIAZ-RODRIGUEZ P,et al.Tuning forkhead box D3 overexpression to promote specific osteogenic differentiation of human embryonic stem cells while reducing pluripotency in a three-dimensional culture system[J].Journal of Tissue Engineering and Regenerative Medicine,2018,12(12):2256-2265.

[22] B

CHNER D,SCHRÖDER D,BETAT N,et al.Apolipoprotein E (ApoE),a Bmp-2 (bone morphogenetic protein) upregulated gene in mesenchymal progenitors (C3H10T1/2),is highly expressed in murine embryonic development[J].BioFactors,1999,9(1):11-17.

[23] TSUSHIMA H,TANG Y J,PUVIINDRAN V,et al.Intracellular biosynthesis of lipids and cholesterol by Scap and Insig in mesenchymal cells regulates long bone growth and chondrocyte homeostasis[J].Development,2018,145(13):dev162396.

[24] HAN Y J,YOU X L,XING W H,et al.Paracrine and endocrine actions of bone-the functions of secretory proteins from osteoblasts,osteocytes,and osteoclasts[J].Bone Research,2018,6:16.

[25] MOLD D E,DINITZ A E,SAMBANDAN D R.Regulation of zebrafish zona Pellucida gene activity in developing oocytes[J].Biology of Reproduction,2009,81(1):101-110.

[26] PASHAY AHI E,WALKER B S,LASSITER C S,et al.Investigation of the effects of estrogen on skeletal gene expression during zebrafish larval head development[J].PeerJ,2016,4:e1878.

[27] ZHOU Y J,RUI L Y.Leptin signaling and leptin resistance[J].Frontiers of Medicine,2013,7(2):207-222.

[28] TSAKOUMIS E,AHI E P,SCHMITZ M.Impaired leptin signaling causes subfertility in female zebrafish[J].Molecular and Cellular Endocrinology,2022,546:111595.

[29] REID I R,BALDOCK P A,CORNISH J.Effects of leptin on the skeleton[J].Endocrine Reviews,2018,39(6):938-959.