不同生长速率花鲈肌肉转录组分析及生长相关基因筛选

范嗣刚，杨文燕，黄 皓,3，王鹏飞，赵 超，闫路路，张 博，邱丽华,3

（1.中国水产科学研究院南海水产研究所/广东省渔业生态环境重点实验室，广东 广州 510300；2.广东省农业技术推广中心，广东 广州 510520；3.上海海洋大学水产与生命学院，上海 201306）

除环境外，鱼类生长还受许多基因的影响和控制。目前已挖掘出一些生长相关基因。生长激素（Growth hormone，GH）-胰岛素样生长因子（Insulinlike growth factor，IGF）生长轴是鱼类生长发育的主要作用通道[1]。增加鱼类体内的GH 含量可促进生长[2]。IGFI和IGFII可促进葡萄糖进入肌肉，加速肌细胞的分裂和分化[3-4]。肌肉生长抑制素（Myo‐statin，MSTN）是转化生长因子-β (Transforming growth factor-beta，TGF-β)超家族的成员，在动物肌肉生长和发育过程中起负调控作用[5]。骨骼肌中的α-肌动蛋白基因在大、小规格虹鳟（Oncorhynchus mykiss）中差异性表达[6]，在快速生长鲇（Ictalurus punctatus）中表达量显著升高[7]。成纤维细胞生长因子（Fibroblast growth factor，FGF）可调控细胞增殖、迁移和分化[8]。转录组测序技术（RNA-sequenc‐ing，RNA-seq）可准确、有效、快速地获得组织转录组序列及表达情况，已广泛用于水产动物生长基因的筛选和鉴定等研究。孙雪等[9]对草鱼（Ctenopharyngodon idella）快长组和慢长组个体的肌肉进行RNA-seq，获得552 个差异表达基因（Differentially expressed genes，DEGs）。Zhang 等对生长快和生长慢的青鱼（Mylopharyngodon piceus）进行转录组分析，获得1 913 个DEGs，在生长和发育相关的代谢途径中富集[10]。Luo 等[11]用转录组测序技术研究生长快和生长慢的鳙（Hypophthalmichthys nobilis），鉴定出15 个关键信号途径和20 个DEGs。Cleveland等[12]通过比较虹鳟（O.mykiss）快、慢生长组的转录组数据，获得36 个DEGs。王玉梅等[13]研究2 种不同生长速率的福瑞鲤（Cyprinus carpio）转录组，获得749 个DEGs，鉴定出肌红蛋白（Myoglobin,mb）、肌球蛋白轻链2（Myosin light chain 2,myl2）等与肌肉生长相关的关键基因。

花鲈（Lateolabrax maculatus）隶属于鲈形目（Perciformes）鮨 科（Serranidae）花鲈属（Lateolabrax），俗称鲈鱼、七星鲈、花寨等，广泛分布于我国渤海至南海的近岸水域以及通海江河水域，是一种广温广盐性经济鱼类，有生长较快、适应性广、味道鲜美等特点[14-15]。目前，已有花鲈生长的研究。Liu 等[16]构建了花鲈首个高密度连锁图谱，筛选出24个与生长相关的数量性状位点（QTL），获得30个生长有关候选基因。胰岛素样生长因子结合蛋白基因（IGFBP）、IGF、生长激素受体基因（GHR）、甲状腺激素受体（Thyroid hormone receptor，TR）和瘦素A（LeptinA，LepA）等与花鲈生长有关的基因序列及其组织表达也有研究[17-18]。但尚未见对花鲈生长方面的转录组研究报道。鱼类肌肉的生长发育直接决定鱼类的总质量。本研究拟对同一养殖环境下生长快与慢的花鲈肌肉组织进行转录组测序，研究二者在转录组水平上的差异，筛选出与生长相关的差异表达基因和信号通路，为解析花鲈生长的分子机制奠定基础。

1 材料与方法

1.1 材料

花鲈同批次繁殖苗种（体长2～3 cm）1 万尾，2020 年3 月取自宁德市富发水产有限公司，在南海水产研究所珠海试验基地池塘中养殖9个月。取花鲈500 尾，其中461 尾生长快速，达到商品规格（500～1 000 g），39 尾生长缓慢，个体短小（90～150 g）。随机选取4 尾生长快[体长(34.22±2.53)cm]和4尾生长慢[体长(17.46±1.78)cm]的花鲈个体，取肌肉组织液氮保存备用。

1.2 RNA提取、cDNA文库构建及测序

用Tirzol 法（Thermofisher scientific，美国）提取花鲈肌肉组织的RNA，用DNase（Takara，日本）消化DNA。用带有Oligo (dT) 的磁珠（杭州尚炫生物科技有限公司）富集mRNA，加入打断试剂将mRNA打断成短片段。以打断后的mRNA 为模板，用六碱基随机引物合成一链cDNA，配制二链合成反应体系合成二链cDNA，并用试剂盒纯化双链cDNA。对纯化的双链cDNA 进行末端修复、加A 尾并连接测序接头、片段大小选择、PCR 扩增。构建的文库用2100 生物分析仪（Agilent，美国）质检合格后，用Illumina HiSeq X 测序仪（Illumina，美国）测序，产生150 bp的双端数据。文库构建和转录组测序均由上海生工生物工程股份有限公司完成。

1.3 测序数据组装

用Trimmomatic 软件[19]分析原始数据（Raw reads），去除reads的接头序列、低质量reads（质量值Qphred＜20）、N 碱基比例大于5%的reads 等，得到纯净读数（Clean reads）。

用Trinity[20]将干净读长从头组装成转录本（Transcript），参数为min_kmer_cov 2，其余参数为默认值。转录本去冗余后，取每个转录本聚类中最长的转录本作为unigenes，以此作为后续分析的参考序列。

1.4 基因功能注释

将获得的unigenes 在NT (NCBI nucleotide sequences)、Nr(NCBI non-redundant protein sequences)、KOG (euKaryotic Ortholog Groups)和Swiss-Prot 等数据库进行比对和注释，e值均小于1e-5。

用Blast 2 GO 软件（version 2.5）对unigenes 进行GO 功能注释（e-value ＜1e-6），用KAAS[21]软件对unigenes进行KEGG注释（e-value ＜1e-10）。

1.5 差异表达基因的鉴定与分析

用Salmon[22]计算基因的表达量。用RSEM 软件将每个样品的clean reads 比对到拼接转录本中，获得读长数量（read count）。用TMM 软件对read count 数据进行标准化处理，用DESeq2 1.12.4 软件包分析生长快和生长慢花鲈的差异表达基因。筛选条件设为：q ＜0.05 且|log2FC|＞2，FC 为差异倍数（Fold change），q为多重假设检验校正后的P值。

用topGO（version 2.24）软件包对DEGs 进行GO 富集分析，用clusterProfiler（version 3.0.5）软件包进行KEGG通路富集分析。

1.6 差异表达基因的验证

随机选取转录组中的13 个差异表达基因进行qRT-PCR验证（表1）。用与转录组测序同批次的RNA进行反转录，获得相应的cDNA。β-actin为内参基因。用SYBR Green Pro Taq HS预混型qPCR试剂盒（艾科瑞，中国）在LightCycler 480荧光定量PCR仪（Roche，美国）上进行qRT-PCR。反应体系为12.5 μL，包含cDNA（40 ng/μL）5.25 μL，上下游引物（10 μmol/L）各0.5 μL，2X SYBR Green Pro Taq HS Premix 6.25 μL。反应条件：95 ℃30 s；94 ℃5 s，60 ℃30 s，共40 个循环。反应结束后，对PCR 产物进行溶解曲线检测，以验证扩增产物的特异性。每次试验做3 个样品，每个样品重复检测3 次。用2-ΔΔCt法计算基因的相对表达量。

表1 用于实时定量PCR验证的引物Table 1 Primers used in Real-time PCR confirmation

2 结果与分析

2.1 花鲈肌肉转录组测序与组装结果

对生长快和生长慢花鲈的肌肉组织分别获得原始序列44 434 288、48 256 446 个，Clean reads 42 928 714、46 975 694 个，Q30值分别为91.65%、91.84%（表2）。测序结果已提交到NCBI SRA 数据库（PRJNA811720）。

表2 花鲈肌肉转录组数据的基本信息Table 2 Summary statistics of Lateolabrax maculatus muscle transcriptome sequencing

两个Clean reads拼接后获得137 960个unigenes；N50为637 bp，N90234 bp；总长70 483 974 bp，最大长度19 775 bp，最小长度201 bp，平均长度510.9 bp；长度≥500 bp unigenes有33 908个，长度≥1 000 bp的有13 373个（图1）。

图1 花鲈unigene长度分布Fig.1 Unigene length distribution in transcriptome of Lateolabrax maculatus

2.2 基因功能注释

unigenes 注释结果如表3 所示。有68 334 个unigenes（49.53%）在Nr 中注释，68 305 个unigenes（49.51%）在Swissprot 中注释，在GO 数据库注释到70 946 个unigenes（51.43%）。有105 310 个unigenes在至少一个数据库中获得注释信息，有5 215 个unigenes在所有数据库中均有注释。

表3 花鲈肌肉转录组测序结果统计Table 3 Statistics of Lateolabrax maculatus muscle transcriptome sequencing

2.3 差异表达基因鉴定与分析

根据q ＜0.05 且|log2FC|＞2 的标准，快速生长组与慢速生长组比较，有10 552 个差异表达基因（图2）。其中，2 064 个unigenes 表达量上调（1 661个unigenes 在Nr 中成功注释），8 488 个unigenes 表达量下调（3 793 个unigenes 在Nr 中成功注释）。在这些DEG中，有一些与肌肉生长和发育相关的基因被筛选出来，如胰岛素样生长因子结合蛋白1（Insu‐lin-like growth factor-binding protein 1，IGFBP1）、肌球蛋白重链（Myosin heavy chain，MHC）、肌球蛋白轻链（Myosin light chain，MLC）、成纤维细胞生长因子结合蛋白3（Fibroblast growth factor-binding pro‐tein 3，FGFBP3）、多表皮生长因子样结构域蛋白9（Multiple epidermal growth factor-like domain pro‐tein 9，MEGFL 9）、成纤维细胞生长因子（Fibroblast growth factor，FGF）、生长激素受体1（Growth hor‐mone receptor 1，GHR1）、肌生成抑制蛋白（Myo‐statin，MSTN）、生肌因子（Myogenic factor，MYF）、肌凝蛋白结合蛋白（Myosin binding protein，MYBP）、生长抑制蛋白（Inhibitor of growth protein，ING）等。

图2 慢速、快速生长组花鲈肌肉的差异表达基因M-A图Fig.2 M-A plot for DEGs comparison between the slowand fast-growing group of Lateolabrax maculatus

2.4 差异表达基因功能注释

对10 552个DEGs进行GO富集分析，在生物过程富集的功能分支（Term）最多，有9 493 个，其次在分子功能中，富集2 345个，在细胞成分中富集1 358个。在每个类别（Ontology）中，显著富集的前10 个功能分支如表4所示。在生物过程中显著富集的前10 个功能分支中，多为前体、多肽等物质的代谢。在细胞成分中，有6 个与核糖体有关，2 个与肌肉有关。在分子功能排名前5 的是核糖体的结构成分（Structural constituent of ribosome）、结构分子活性（Structural molecule activity）、rRNA 结 合（rRNA binding）、RNA结合（RNA binding）、细胞色素c氧化酶活性（Cytochrome-c oxidase activity）。

表4 不同规格花鲈差异表达基因在GO功能三大类别中富集前10的条目分析Table 4 Analysis of the first 10 terms about DEGs for large and small size of Lateolabrax maculatus in the three majorcategories of GO function

在KEGG 富集分析中，共富集到255 个通路。上调的DEGs 富集最多的前5 个KEGG 通路为催产素信号通路、胰高糖素信号通路、脂肪酸降解、腺苷酸活化蛋白激酶（AMPK）信号通路、碳代谢，下调的DEGs富集最多的前5个KEGG通路为核糖体、氧化磷酸化、心脏肌肉收缩、抗原加工和提呈、氨基酸生物合成（图3）。

图3 差异表达基因显著富集的前10个KEGG通路Fig.3 Top 10 enriched KEGG pathways in DEGs

2.5 差异表达基因的验证

随机挑选13个DEGs中，有8个上调基因和5个下调基因。图4 表明，所有的差异表达基因的表达模式与测序结果一致，表明测序结果的可信度高。

图4 差异表达基因的实时定量PCR验证Fig.4 qRT-PCR validation analysis of DEGs

3 讨论

本研究中，同一批花鲈苗种（体长2～3 cm）放入池塘中养殖，在养殖水体、饲料等条件一致的情况下，9 个月后获得两种规格不同的花鲈个体，可排除环境因素对生长的影响，确定某些基因和信号通路对生长的调控作用。

鱼类的肌肉生长包括肌纤维增生和肥大，受多种基因和信号通路的调控[23]。GH/IGF 轴在调节鱼类生长方面起重要作用[1]。本研究筛选出该轴上的GHR1、IGFBP1、IGFBP7和IGFBP4等 DEGs。GHR1是GH的受体，能促进细胞生长、增殖，减少蛋白降解[24]。鱼类的GHRs 是单一跨膜蛋白，包括胞外结构（包含6 或7 半胱氨酸残基）、跨膜结构和胞内结构[25]。虹鳟（O.mykiss）GHR1 主要通过信号转导及转录激活因子5（STAT5）起作用[26]。花鲈GH序列已被克隆出来[17]，但未见花鲈GHR的研究报道。本研究发现，生长快的花鲈GHR1表达量显著上调，说明GHR1能促进花鲈生长。IGFBP1受胰岛素的负调控，与鱼类的生长负相关[27]。IGFBP4在体内是一种抑制剂，通过抑制IGF的活性，阻止受体与IGF 结合[28-29]。IGFBP7，又 名IGFBP 相关蛋白1（IGFBP-related protein 1，IGFBP-rP1），与IGF 结合能力弱，与胰岛素结合能力强[30-31]。在生长快和生长慢的大口黑鲈（Micropterus salmoides）转录组比较研究中，IGFBP1也被筛选出来[32]。肌球蛋白重链和肌球蛋白轻链均为肌球蛋白的重要组成部分，参与了肌肉的生长发育过程[33]。生长抑制蛋白家族包括ING1、ING2、ING3、ING4和ING5等5 个成员[34]。本研究检测出ING3、ING4和ING5等DEGs。生长抑制蛋白能抑制细胞增殖、迁移和扩散[35]。Liu等[16]基于花鲈高密度遗传连锁图谱，发现FGFR4、FGF12a、FGF18 等3 个FGF 家族成员与花鲈生长有关。本研究DEGs 中有4 个FGF 家族成员。FGFBP3 是一种分泌分子伴侣，能影响碳水化合物和脂质代谢，在老鼠中敲除FGFBP3后，能改变脂质代谢[36]。FGF 6 的功能主要有肌源干细胞迁移，肌肉分化、增生[37]。草鱼的FGF6 表达量与肌纤维直径正相关[38]。FGFRL1 通过结合FGF，对细胞增殖起负调控作用[39]，对骨形成和细胞分化起正调控作用[40]。在草鱼饥饿再投喂实验中，FGFRL1表达量显著上调[41]。在哺乳动物中，FGF2是一种原型血管生成生长因子，能诱导血管生成和新血管形成[42]。在斑马鱼（Danio rerio）中也发现FGF2 能促进血管生成[43]。

在GO的细胞组分中，有肌原纤维和肌节等2个DEGs 显著富集，KEGG 途径中有心肌收缩的DEGs显著富集。在福瑞鲤、青鱼和鳙等鱼类中也发现有DEGs 匹配到肌原纤维和心肌收缩等GO 和KEGG富集中[10-11,13]。由此可见，在不同鱼类中，肌原纤维和心肌收缩与肌肉生长都紧密相关。

研究表明，物质和能量代谢能影响生物生长[44]。本研究也证实这一点：DEGs的GO 和KEGG 富集结果多集中在物质和能量代谢方面。在GO 富集的生物过程中，排在前10 的是蛋白、核酸等物质代谢。在生长快速的草鱼和鳙中，富集在代谢过程的DEGs表达量显著上调[11,45]。

本研究部分KEGG 分析结果与其他鱼类的研究结果有相同之处，如糖酵解/糖异生途径[32]、脂肪酸降解途径[11,13,46]、AMPK 信号通路[10]、碳代谢和脂肪酸代谢[10-11]、精氨酸和脯氨酸代谢[10-11,44]、D-谷氨酰胺和D-谷氨酸盐代谢[47]。这些KEGG 途径均与肌肉生长发育有关，如糖酵解是糖原产生ATP 的过程，糖异生是非糖有机物转变成葡萄糖的过程。虹鳟再投喂试验发现有一些基因与糖酵解/糖异生有关[48]。AMPK在调控肌肉生长和再生方面起重要作用，能促进肌肉汇总葡萄糖的吸收和脂肪酸的氧化，为肌肉收缩提供能量[49]。摄入精氨酸和脯氨酸等氨基酸能增强肌肉生长，减少过多的脂肪[50]。今后可进一步研究这些KEGG 途径与鱼类肌肉生长发育的关系。

本研究中，分别有6 个细胞组分和2 个分子功能与核糖体有关，KEGG 富集中下调最显著的是核糖体，可见核糖体对鱼类生长的重要性。事实上，核糖体是细胞合成蛋白质的主要场所，RNA 翻译到蛋白质这一过程即在核糖体中完成。核糖体蛋白和rRNA 组成两个大小不同的核糖体亚基。研究发现，核糖体蛋白能控基因转录、mRNA 翻译，影响细胞分化、增殖和凋亡，对肌肉生长有重要作用[51]。