豹纹鳃棘鲈(Plectropomus leopardus)，俗称东星斑，是石斑鱼类中的高端养殖品种。与大多数石斑鱼类一样，豹纹鳃棘鲈的人工育苗和中间培育需要历经鱼苗分筛和饵料驯化，此过程中同一批次的仔鱼、稚鱼和幼鱼会多次出现生长差异分化的现象(图 1)，是导致其养殖周期拉长、影响养殖经济效益的主要因素之一。为解决上述问题，多数石斑鱼的育苗过程中都会进行“标粗”，并筛选生长差异显著的“鱼头”、“鱼中”和“鱼尾”，在分筛鱼苗后按照不同的体型规格分别培育(彭树锋等, 2012)，一定程度上提高了鱼苗成活率并缩短了整体养殖周期。然而，“标粗”过程不仅受人为因素影响较大，且鱼苗自身差异也影响“标粗”成效。尽管这种生长差异分化现象受环境因素的影响较大(符书源等, 2009)，但其仍存在尚未明确的潜在遗传因素。

图 1(Fig. 1)

图 1 豹纹鳃棘鲈苗种中间培育过程中的生长差异分化现象 Fig.1 Growth differentiation phenomenon during intermediate cultivation of P. leopardus juveniles

鱼类生长性能属于复杂性状，生长差异性状的分子调控和遗传作用机制解析已成为鱼类育种研究领域中的热点(Enberg et al, 2012; Guerrero-Cozar et al, 2021)。生长性状的优劣是评价鱼类养殖经济效益的重要指标，同时受遗传内因与环境外因影响，因此, 开展养殖鱼类生长相关性状的遗传机制及调控策略研究将为其品种改良和育种提供参考(苏军虎等, 2012)。全转录组测序技术因其具有分析辨别基因组结构功能、构建生物系统遗传网络以及可靠性高等优点，被广泛应用于动植物遗传机制解析研究(Zhao et al, 2021; Niu et al, 2022; Li et al, 2022)。全转录组分析在编码RNA研究的基础上，进一步涵盖了多种非编码RNA及各种类型RNA间的调控网络分析(石田培等, 2019)。Ye等(2021)对具有不同生长速度的29条6月龄草鱼(Ctenopharyngodon idella)的脑部和肝胰腺组织进行全转录组分析，筛选出多个与草鱼生长性状显著相关的mRNA、miRNA和lncRNA关键因子，为探索草鱼生长速率的分子遗传机制奠定了基础。Wang等(2022)通过全转录组测序构建了半滑舌鳎(Cynoglossus semilaevis)性别大小异形性状的lncRNA-miRNA-mRNA网络，揭示了其雌雄生长差异遗传机制。Wu等(2022)基于全转录组测序的野生型和黄色突变型虹鳟(Oncorhynchus mykiss)皮肤中lncRNA和circRNA介导的ceRNA调控网络分析，解析了二者之间先天免疫的差异，进一步扩展了对虹鳟鱼先天免疫系统的理解，并为虹鳟免疫机制研究和疾病抗性育种奠定了基础(Wu et al, 2022)。

lncRNA是一类超200 bp且不具编码功能的RNA分子，其通过与RNA、DNA或蛋白质相互作用调控靶基因的表达，广泛参与水产动物生长、肌肉品质、生殖、胚胎发育、性别分化及免疫等多种生理功能(宋飞彪等, 2020)。由于复杂性状在遗传层面通常受多个分子共同作用，同时利用mRNA和lncRNA开展调控机制研究，并运用系统生物学网络的方法对基因表达数据进行深度挖掘，构建与性状相关的基因网络已成为生物信息学分析的重要研究方向(汪涛等, 2014; 陈华峰等, 2019; Zhang et al, 2021)。加权共表达网络分析(WGCNA)因其在样本间基因关联的网络关系构建准确性高，被广泛应用于转录表达方面的研究(Li et al, 2018; Yao et al, 2023)。目前，尚未见关于豹纹鳃棘鲈生长差异分化性状共表达网络分析的报道。

本研究采用R语言的WGCNA程辑包，对豹纹鳃棘鲈苗种“标粗”过程中3个生长阶段共计18个肌肉组织样本的的转录组表达数据进行分析，构建与生长差异分化相关的lncRNA与mRNA共表达模块。通过差异表达基因分析和目标模块内基因的GO和KEGG富集分析，确定目标模块中的核心基因，深入挖掘与生长差异分化相关的基因及调控网络，以期为豹纹鳃棘鲈的生长差异分化遗传机制解析和良种选育提供依据。

本研究在豹纹鳃棘鲈生长差异分化的全转录组测序数据的基础上进行。课题组于2022年在海南省海洋与渔业科学院海研热带海水鱼类良种场开展了豹纹鳃棘鲈苗种培育工作，在鱼苗“标粗”过程中的3个生长阶段(42、70和91日龄)分别采集了“鱼头”中生长差异分化显著鱼苗的肌肉组织，每个时期大小规格各3个生物学重复，共计18个样本送至广州基迪奥生物科技有限公司进行全转录组测序。测序获得的mRNA和lncRNA原始下机数据于2023年4月返回，后续生物信息学分析在海南省海洋与渔业科学院进行。

通过去除含有接头、不确定性占比大于10%的碱基以及低质量(Q_phred≤10)的碱基数占整条read 50%以上的reads对原始下机数据(raw data)进行严格的过滤，获得高质量测序数据(clean data)以提高后续转录组分析的准确性。利用HISAT2软件(Kim et al, 2015)将过滤后的clean reads高效比对到豹纹鳃棘鲈参考基因组(https://db.cngb.org/search/project/CNP0000859)，通过StringTie软件(Pertea et al, 2015)进行基因水平的定量，获取表达谱数据。采用CPC2 (Kang et al, 2017)、CNCI (Sun et al, 2013)和LGC (Wang et al, 2019) 3个软件预测样本中的lncRNA，为降低鉴定结果的假阳性，取3个软件分析结果的交集lncRNA用于后续分析。采用FPKM (fragments per kilobase of transcript per million fragments mapped)法(Florea et al, 2013)归一化reads数和转录本长度，使不同样本间基因表达具有可比性。对苗种“标粗”过程中3个生长阶段中大规格(L)和小规格(S)样本的lncRNA和mRNA表达谱数据进行对比(L-Stage 1 vs S-Stage 1、L-Stage 2 vs S-Stage 2以及L-Stage 3 vs S-Stage 3)，根据各转录本的FPKM值，利用DESeq2软件(Love et al, 2014)进行差异表达分析，并以Fold Change≥2且FDR (false discovery rate)≤0.05为阈值标准筛选差异表达基因。

采用R语言(v4.0.5)的WGCNA程辑包(v1.72.5) (Langfelder et al, 2008)构建包括全部mRNA和lncRNA的无尺度网络。以17个(剔除1个S-Stage 2离群样本)生长差异分化转录组中标准化FPKM值的表达矩阵作为输入数据，筛选过滤表达值方差 > 0和缺失≤10%的基因构建加权共表达网络。利用WGCNA中的pickSoftThreshold函数计算无尺度拓扑拟合指数，筛选共表达网络的软阈值，并利用power- Estimate函数估计其最优值；选择最优软阈值后，使用TOMsimilarityFromExpr函数计算拓扑重叠矩阵(topological overlapmatrix, TOM)；利用拓扑相异矩阵(1-TOM)进行基因分层聚类(hclust)，并采用动态剪切算法(cutreeDynamic)划分模块(minClusterSize=50)；最后，使用mergeCloseModules合并相似度大于75%的基因模块(mergeCutHeight=0.25)，确定最终的共表达基因模块。对共表达基因模块中的所有基因表达值进行主成分分析，计算出各模块的特征向量(module eigengene, ME)，统计各模块的ME值与生长差异分化性状(不同生长规格)的相关性，按照显著性选取相关系数绝对值大于特定阈值的模块作为豹纹鳃棘鲈生长差异分化特异性功能模块进行后续分析。

采用R语言clusterProfiler程辑包(Yu et al, 2012)，基于eggNOG数据库对选取的特异性功能模块中差异表达的mRNA及lncRNA分别进行基因本体论(GO)和京都基因与基因组百科全书(KEGG)功能通路富集分析，并以P < 0.05为阈值富集获得显著性的GO term和KEGG pathway。

通过WGCNA的exportNetworkToCytoscape函数输出特异性功能基因子网络的边列表和节点列表数据，筛选基因权重(Weight)和连通度(Degree)排序靠前的mRNA和lncRNA，基于功能基因注释结果，选择前65个基因和lncRNA，利用Cytoscape(v3.10.1)软件(Shannon et al, 2003)中的CytoHubba插件估计关键性(hub)大小，并进行关键基因筛选及网络可视化。

分别从42、70及91日龄相同生长规格鱼苗中挑选的生长差异分化个体作为实验用鱼，其全长表型信息见表 1。基于全转录组测序获取mRNA和lncRNA基因表达谱，利用DESeq2软件进行不同生长阶段的差异基因表达分析，去除各生长阶段大小规格样本比对时重复出现的差异表达mRNA和lncRNA，共获得5 944个差异表达已知mRNA、302个novel mRNA以及367个差异表达lncRNA(表 2)，其中，各生长阶段差异表达上调或下调的mRNA和lncRNA数统计见图 2。3个生长阶段同时检出的差异表达mRNA占全部差异表达mRNA的17.3%，大于差异表达lncRNA的占比5.7% (图 3)。

表 1(Tab.1)

表 1 生长差异分化样本的全长表型平均数和标准差 Tab.1 Mean and standard deviation of phenotype of total length for growth differentiated samples 日龄 Days post hatching 大小规格 Size specification 平均数±标准差 Mean±SD/cm 42 L 4.00±0.10 S 3.17±0.12 70 L 5.77±0.15 S 4.30±0.17 91 L 7.40±0.10 S 4.00±0.10

表 1 生长差异分化样本的全长表型平均数和标准差 Tab.1 Mean and standard deviation of phenotype of total length for growth differentiated samples

表 2(Tab.2)

表 2 mRNA和lncRNA的差异表达统计 Tab.2 Statistics of differential expression of mRNA and lncRNA 差异表达类别 DE type L-Stage 1 vs S-Stage 1 L-Stage 2 vs S-Stage 2 L-Stage 3 vs S-Stage 1 去重复 Deduplication 差异表达已知mRNA数No. of DE known mRNAs 4 220 2 311 3 232 5 944 差异表达novel mRNA数No. of DE novel mRNAs 221 109 198 302 差异表达lncRNA数No. of DE lncRNAs 232 76 185 367

表 2 mRNA和lncRNA的差异表达统计 Tab.2 Statistics of differential expression of mRNA and lncRNA

图 1(Fig. 1)

图 2 不同生长阶段差异表达上调或下调mRNA和lncRNA数统计 Fig.2 Number of up-regulated or down-regulated mRNA and lncRNA counts at different growth stages

图 1(Fig. 1)

图 3 不同生长阶段差异表达mRNA和lncRNA韦恩图 Fig.3 Venn plot of DE mRNA and lncRNA at different growth stages

为了减少构建无尺度分布网络模块的计算复杂度和计算时间，筛选最适软阈值为22，此时，无尺度网络的拟合指数R² > 0.8，平均连通性为94.23 (附图 1)。基于相异度矩阵dissTOM进行层次聚类得到基因聚类树，利用混合动态树剪切法划分出37个基因模块(见附图 2中dynamic tree cut)，随后以相似度大于75%为标准合并相关性高的模块并获得16个共表达基因模块(见附图 2中merged dynamic)。共表达网络中各模块包含的mRNA和lncRNA及其差异表达数见表 3，其中，MEmagenta模块中的mRNA和lncRNA数最多，包含6 835个mRNA和460个lncRNA，该模块中mRNA和lncRNA差异表达数也是最多的，包含2 879个差异表达mRNA和196个差异表达lncRNA；MEyellowgreen模块中mRNA数和差异表达数最少，为66和18；MEdarkmagenta模块中不包含lncRNA，MEdarkorange和MEskyblue模块中未发现差异表达lncRNA；MEgrey模块中包含的是没有分配到其他模块中的mRNA和lncRNA集。

表 3(Tab.3)

表 3 共表达基因模块内mRNA和lncRNA及差异表达数统计 Tab.3 mRNA, lncRNA and differential expression counts in co-expression modules 基因模块 Gene module mRNA数 No. of mRNA lncRNA数 No. of lncRNA 差异表达mRNA数 No. of DE mRNA 差异表达lncRNA数 No. of DE lncRNA 基因模块 Gene module mRNA数 No. of mRNA lncRNA数 No. of lncRNA 差异表达mRNA数 No. of DE mRNA 差异表达lncRNA数 No. of DE lncRNA MEmagenta 6 835 460 2 879 196 MElightyellow 226 20 56 3 MEturquoise 4 588 115 878 20 MEdarkolivegreen 197 7 36 1 MEgrey 3 821 252 1172 74 MEdarkred 195 6 59 3 MEgreenyellow 727 46 452 37 MEdarkorange 144 7 43 0 MEblack 547 21 104 15 MEskyblue 129 1 23 0 MEdarkturquoise 469 42 121 4 MEsaddlebrown 120 9 36 2 MEsalmon 376 13 124 5 MEdarkmagenta 92 0 30 0 MEdarkgrey 254 5 201 4 MEyellowgreen 66 1 18 1

表 3 共表达基因模块内mRNA和lncRNA及差异表达数统计 Tab.3 mRNA, lncRNA and differential expression counts in co-expression modules

由豹纹鳃棘鲈不同生长阶段的生长差异分化性状与构建的模块相关关系(图 4)可知，本研究共鉴定到3个与生长差异分化特异相关的共表达模块，筛选标准为3个生长阶段都符合(0.65≤|r|≤1, P≤0.01)。其中，MEmagenta模块在3个阶段的差异分化都呈正相关，而MEgreenyellow模块都呈负相关；MEblack模块42日龄时与生长差异分化呈正相关，而70日龄和91日龄时则呈负相关。为了确定模块划分的可靠性，以上述3个共表达模块(MEmagenta、MEblack和MEgreenyellow)为特异性模块进行深入分析。

图 1(Fig. 1)

图 4 模块与性状的相关性热图 Fig.4 Heat map of the correlation between modules and traits

为了探究特异性模块中差异表达mRNA及lncRNA的主要生物学功能和参与的基因调控通路，本研究对MEmagenta、MEblack和MEgreenyellow模块中的差异表达基因进行了GO和KEGG富集分析。GO富集分析表明，MEmagenta和MEgreenyellow模块中的差异表达基因显著富集于蛋白修饰、折叠、去泛素化以及组装等肌肉蛋白形成等相关的生物学过程(图 5)，而MEblack模块中的差异表达基因则富集于横纹肌组织发育调控、主轴组织形成等生物学过程(图 5)。在KEGG富集分析结果中，MEmagenta和MEgreenyellow模块中的差异表达基因显著富集在核糖体形成、蛋白酶体、肌动蛋白细胞骨架调控及DNA修复蛋白等与生长分化相关的通路，而MEblack模块中的差异表达基因则与肌细胞的肾上腺素能信号传导、肌醇磷酸代谢和泛酸盐和辅酶a生物合成等通路相关(图 6)。

图 1(Fig. 1)

图 5 MEmagenta、MEgreenyellow和MEblack模块差异表达基因GO富集图 Fig.5 GO enrichment diagram of differential expression genes in MEmagenta, MEgreenyellow and MEblack module

图 1(Fig. 1)

图 6 MEmagenta、MEgreenyellow和模块差异表达基因KEGG富集图 Fig.6 KEGG enrichment diagram of differential expression genes in MEmagenta, MEgreenyellow and MEblack module

构建3个与豹纹鳃棘鲈生长差异分化特异性相关模块中的基因互作网络，将WGCNA分析结果中特异性模块的节点node和边edge文件导入Cytoscape软件，以权重和连通度排序靠前的65个基因(含lncRNA)进行互作网络可视化，并计算关键性(hub)值，以hub值靠前的10个基因作为模块网络中的hub基因(图 7~图 9)。MEmagenta模块中的10个hub基因分别为psmd6(26S proteasome non-ATPase regulatory subunit 6)、nmt1(Glycylpeptide N-tetradecanoyltransferase 1)、als2(alsin Rho guanine nucleotide exchange factor)、csn1(COP9 signalosome complex subunit 1)、nu107 (nuclear pore complex protein Nup107)、brcc3(Lys-63- specific deubiquitinase BRCC36)、pyrd1(pyridine nucleotide-disulfide oxidoreductase domain-containing protein 1)、rb6i2(ELKS/Rab6- interacting/CAST family member 1)、atg16 (autophagy-related protein 16)和pp2aa (serine/threonine-protein phosphatase 2A catalytic subunit alpha isoform)，另包含4个lncRNA (见图 7中三角形节点)；MEgreenyellow模块中的9个hub基因分别为kank1(KN motif and ankyrin repeat domain-containing protein 1)、bdh2(3-hydroxybutyrate dehydrogenase type 2)、cra1b(collagen alpha-1(XXVII) chain B)、pacn3(protein kinase C and casein kinase Ⅱ substrate protein 3)、leg(beta-galactoside-binding lectin)、ada12(disintegrin and metalloproteinase domain- containing protein 12)、ulk1(serine/threonine-protein kinase ULK1)、pedf (pigment epithelium-derived factor)和pp2aa(serine/threonine-protein phosphatase 2A catalytic subunit alpha isoform)，1个hub lncRNA (MSTRG.15660)和另外10个lncRNA(图 8中三角形节点)；MEblack模块中的8个hub基因分别为at131 (manganese-transporting ATPase 13A1)、atrx (transcriptional regulator ATRX)、xpo5(Exportin-5)、hnrpl (heterogeneous nuclear ribonucleoprotein L)、pyr1(pyrabactin resistance 1)、ima5(importin subunit alpha-5)、hdac3(histone deacetylase 3)和tia1 (nucleolysin TIA-1 isoform p40)，2个hub lncRNA (MSTRG.8694和MSTRG.1896)和另外4个lncRNA (图 9中三角形节点)。

图 1(Fig. 1)

图 7 MEmagenta模块基因的互作网络 Fig.7 Interaction network of genes in MEmagenta module

图 1(Fig. 1)

图 8 MEgreenyellow模块基因的互作网络 Fig.8 Interaction network of genes in MEgreenyellow module

图 1(Fig. 1)

图 9 MEbalck模块基因的互作网络 Fig.9 Interaction network of genes in MEblack module

生长差异分化是鱼类发育过程中广泛存在的一种现象，这种群体中相同生长发育阶段内的生长差异也是鱼类主要的个体差异之一，受遗传因素和环境因子共同影响。豹纹鳃棘鲈的人工育苗和中间培育过程中的个体间生长差异极其显著，决定其差异分化的遗传机制尚不清楚。目前，Zhou等(2024)开发了豹纹鳃棘鲈的20K SNP array芯片，不仅为其重要经济性状的遗传作用机制解析提供了一种高通量的检测工具，而且对体长和体质量指数性状进行了全基因组关联分析，筛选到多个与生长性状显著关联的遗传位点；Wang等(2023)利用全基因组重测序对307尾豹纹鳃棘鲈的体质量、全长、体长、体高和体厚等5个生长相关性状进行全基因组关联分析，共定位到12个与生长性状显著关联的SNP位点，为豹纹鳃棘鲈的分子辅助育种提供了理论基础；Yang等(2022)基于转录组测序揭示了豹纹鳃棘鲈的不同流速下分子应答机制；Hu等(2023)通过豹纹鳃棘鲈不同体色进行比较转录组分析解析其体色的遗传机制。上述研究中，豹纹鳃棘鲈生长性状遗传机制解析主要基于DNA层面的全基因组关联分析，而在RNA层面则主要是对抗逆和品质性状的研究。

豹纹鳃棘鲈人工育苗和中间培育过程中早期鱼苗受非遗传效应(水温、水质变化等)影响较小，且生长显著差异分化时间间隔短，适合进行遗传机制解析研究。因此，本研究选择早期鱼苗3个连续进行大小规格分筛的时间点(42、70和91日龄)内差异分化个体的肌肉组织进行全转录组测序，并以WGCNA构建了豹纹鳃棘鲈生长差异分化相关lncRNA-mRNA加权基因共表达网络，通过网络模块与生长差异分化性状的显著关联以及模块内关键基因的深入挖掘，筛选出多个与生长差异分化相关的基因和lncRNA。其中，MEmagenta模块中包含差异基因和lncRNA数最多，模块内的核心基因psmd6在斑马鱼(Danio rerio)中受volvox突变体编码，通过降解多泛素化蛋白调节细胞功能影响斑马鱼晶状体上皮细胞增殖和晶状体纤维细胞分化(Imai et al, 2010)；nmt1是哺乳动物发育必需的基因，为早期胚胎发生的主要n-肉豆蔻酰基转移酶(Yang et al, 2005)；als2基因在肌肉发育过程中发挥重要作用，该基因的突变可导致常染色体隐性肌萎缩性侧索硬化症及相关疾病(Devon et al, 2005)；brcc3作为一种细胞周期调节基因，参与磷酸化丝裂原活化蛋白激酶的产生及细胞增殖等过程(Boudreau et al, 2007)。值得注意的是，MEmagenta模块基因互作网络中有多个与psmd6同一基因家族的基因，包括psmd4、psmd7和psmd8等，表明该基因家族可能广泛参与了豹纹鳃棘鲈的生长差异分化过程。Ahongo等(2005)对虹鳟(Oncorhynchus mykiss)产卵后肉质恢复期间的肌肉进行了转录组研究，发现多个涉及蛋白酶体的psmd基因家族基因在虹鳟产卵时选择表达并于产卵后下调表达。上述研究表明，MEmagenta模块内基因主要是通过参与肌肉细胞增殖、分化、凋亡、生长以及正常条件下应激反应等生物学过程，进而影响豹纹鳃棘鲈个体间的生长差异分化。

MEgreenyellow模块内的核心基因kank1最初被认定为肾癌患者染色体9p上的一个候选肿瘤抑制基因，在肾细胞系中过表达时以抑制细胞生长，此外，kank1还通过RhoA信号传导调控肌动蛋白聚合、肌动蛋白应激纤维形成和细胞迁移(Hensley et al, 2016)；ada12(adam12)在小鼠(Rattus norvegicus)肌肉和脂肪生成中起重要作用，其在斑马鱼中合成序列和氨基酸序列上与哺乳动物保持一致，并在心血管系统、红系祖细胞、脑和颌软骨中强烈表达，且adam12基因敲除的斑马鱼在幼年期表现出体型缩小，但没有明显的形态缺陷(Tokumasu et al, 2016)。MEblack模块内的核心基因at131(atp13a1)在小鼠中对抗病毒蛋白MAVS的稳定性至关重要，at131基因敲除的小鼠表现出生长迟缓和胚胎致死等问题(Yang et al, 2024)；hdac3基因为斑马鱼肝脏形成所必需，其主要作用是抑制生长分化因子11 (gdf11)，gdf11是细胞增殖的负调节因子，也是hdac3的特异性转录靶点(Farooq et al, 2008)。2个特异性模块中核心基因ada12及at131的敲除分别导致斑马鱼和小鼠的生长缺陷，直接影响体型规格和生长速度，因此，可以将ada12与at131作为重要候选基因进行豹纹鳃棘鲈生长差异分化分子调控机制的深入研究。3个特异性共表达模块中，MEgreenyellow模块内参与了基因互作网络构建且hub值排序靠前的lncRNA最多，共有11个，表明该模块中的lncRNA相较于MEmagenta模块和MEblack模块发挥了更重要的作用。

本研究侧重分析了与豹纹鳃棘鲈生长差异分化显著相关的3个共表达模块，其他模块中也可能存在重要的基因以及lncRNA参与豹纹鳃棘鲈人工育苗和中间培育过程的生长差异分化过程。此外，本研究所分析的核心基因和lncRNA在豹纹鳃棘鲈中的生物学功能尚不明确，有待后续分子生物学实验进一步验证。

附图 1(FigureS1)

附图 1 软阈值的筛选 FigureS1 Screening of soft threshold (power)

附图 2(FigureS2)

附图 2 基因聚类树及共表达模块构建 FigureS2 Clustering dendrograms of genes and co-expression module division