基于转录组方法鉴定贵州白山羊肌肉生长发育相关上调基因

安清明，王 星，宋兴超，孟金柱，赵园园，吴震洋

（铜仁学院，贵州省梵净山地区生物多样性保护与利用重点实验室，贵州铜仁 554300）

羊肉因其独特的风味和较高的营养价值受到消费者喜爱。近年来，我国羊肉需求量逐年增长。国家统计局统计数据显示，2020 年我国羊肉产量为492.31 万t，较2012 年的404.50 万t 增加了近90 万t，但人均年占有羊肉量只有3.52 kg，较2012 年的人均3.32 kg 几乎没有增长，表明我国羊肉产量仍存在巨大的缺口。如何提高肉羊产肉量仍是育种和养殖工作的重点关注方向。研究发现，动物肌肉组织的生长发育状况直接影响动物总体重，通常认为肌肉的重量可占哺乳动物胴体重的50% 以上[1]。研究认为，动物大多数与肌肉性状相关的因素都具有显著的遗传效应，如骨骼肌生长发育、肌内脂肪含量、嫩度等指标[2-4]，一定程度上均受相应基因的调控。如何高效准确的筛选出与动物肌肉生长发育相关的基因是分子标记育种的一个难点，但随着测序技术的不断发展，转录组测序技术被广泛应用到畜禽重要性状基因筛选和分子机制研究中。如刘瑞莉等[5]基于转录组测序技术筛选布莱凯特黑牛和鲁西黄牛背最长肌差异基因，筛选出296 个上调基因和258 个下调基因，最终推测MYL2、MYL3、MYLPF基因可能在肉牛骨骼肌生长发育过程中发挥重要作用。徐盼等[6]对苏姜猪不同部位、不同生长时期骨骼肌进行转录组测序，共发现1 655 个差异表达基因，最终筛选鉴定FOXO1、FOXO3、FOXO4、MYF6、A-FABP、H-FABP基因可作为改善猪肉质性状的候选基因。在绵羊育种中，也有研究通过转录组测序技术筛选了部分生产性状相关的候选基因，如史金平等[7]对杂交绵羊群体进行转录组测序研究，挖掘出397 个差异表达基因，筛选出ADIPOQ、IGFBP7、ACTC1基因作为肉质性状相关的候选基因。在山羊生长性能研究中，曹艳红等[8]通过对努比亚山羊和隆林山羊进行转录组测序发现ACACA、MFGE8、MYH3、ACTN2基因是山羊产肉性状的关键基因，可能在肉质脂肪含量、肌肉生长发育等方面发挥重要作用。孟宪然等[9]筛选发现ADIPOQ、PDK4、CD36基因可作为改良绒山羊肉品质的候选基因。上述研究表明，从遗传学角度分析畜禽生产相关性状已十分普遍，但目前针对贵州白山羊肉用性能的研究仍然较少，仍需进一步深入研究。

地方品种山羊相较于商品羊生产性能较低，但对当地的环境适应性好，尤其对地方特定的疾病具有一定耐受性，是发展地方畜牧业的重要种质资源。贵州白山羊是在云贵高原特殊的生态环境下，经过长期自然选择和人工选择形成的特有品种，具有耐粗饲、性成熟早、繁殖力强、适应性强、肉质细嫩、无膻味等优点，但同时存在个体较小、生长速度慢、产肉量不高等缺点。因此，本研究选用贵州白山羊为研究对象，通过屠宰实验分析不同体重山羊的屠宰情况，结合转录组测序技术挖掘影响山羊生长性状的重要基因，探索影响贵州白山羊生长发育的基因调控网络，以期为肉用山羊生长发育提供重要的参考数据。

1 材料与方法

1.1 实验动物 受试贵州白山羊饲养于贵州省铜仁市沿河县华珍牧业有限责任公司养殖试验基地，所选山羊为同等饲养水平下采用半舍饲加放牧方式进行饲养，选取同一种群中6 只2 周岁龄健康、体重不同的贵州白山羊（优势群体3 只，体重较高；劣势群体3 只，体重较轻；组内活体重相近），按照国家实验动物屠宰标准进行屠宰并测定屠宰性能（宰前活重、胴体重、屠宰率、净肉重和净肉率），采集背最长肌组织样，迅速放入液氮中，带回实验室转入-80℃冰箱保存，用于转录组测序和qRT-PCR。

1.2 总RNA 提取cDNA 文库建立 背最长肌中提取总RNA，利用Agilent Bioanalyzer 2100 检测完整性，Qubit 2.0 Flurometer 测量RNA 浓度，RNA 样品检查合格后通过Illumina Hiseq 2500 平台完成测序，测序分析由北京诺禾致源生物信息科技有限公司完成。

1.3 测序数据处理 将测序得到的原始数据进行评估，去除含有带接头和低质量的片段，利用TopHat 2 软件与已知山羊参考基因组（GCF_001704415.1）进行比对，获取测序拼装序列相应特征信息，利用Cufflinks 软件对转录本数据组装、注释和合并，利用DESeq2 软件对mRNA 进行差异表达分析，利用基因本体GO 数据库对筛选出的差异表达基因进行功能注释，利用KEGG 数据库进行DEGs 的信号通路富集分析。

1.4 引物合成与qRT-PCR 验证 以NCBI 数据库中山羊FGF11、ADIPOQ、NOTCH1、HBEGF基因的mRNA序列为参考序列，以GAPDH基因为内参基因，应用Primer Premier 5.0 软件设计相应引物，由生工生物工程（上海）股份有限公司合成，引物序列见表1。

表1 荧光定量引物基因列表

qRT-PCR 验证差异表达mRNA 相对表达水平，采用3 个样本重复，应用北京全式金生物技术有限公司生产的荧光定量试剂盒进行qPCR 扩增分析，扩增体系为20 μL，其 中cDNA模板1 μL（1 μg/L），2×qPCR super mix 10 μL，上、下游引物各0.5 μL（10 μmol/L），dd H2O 8 μL。qPCR 反应程序为：94℃预变性1 min；94℃10 s，59℃30 s，72℃10 s，46 个循环。使用2-ΔΔCT法计算各基因间的相对表达水平。

1.5 数据统计与分析 屠宰数据经Excel 整理后，采用SPSS 22.0 软件对6 只2 周岁龄健康、体重不同的贵州白山羊宰前活重、胴体重、屠宰率、净肉重和净肉率进行方差分析。表达量分析ΔΔCt=[（Ct目的基因-CtGAPDH）]实验组-[（Ct目的基因-CtGAPDH）]对照组。结果应用SPSS 22.0软件One-way ANOVA 进行显著性检验分析，P＜0.05 时，差异显著。

2 结果与分析

2.1 贵州白山羊屠宰性能差异分析 实验测定不同体重贵州白山羊2 岁龄屠宰性能（AG 为体重优势群体，DG为体重劣势群体），优势群体的活体重、胴体重和净肉重显著高于劣势群体，表明优势群体屠宰性能高于劣势群体（表2），也说明后续测序试验样品具有可靠性。

表2 贵州白山羊屠宰性能差异分析

2.2 RNA-seq 数据比对分析 由表3 可知，最终6 个样本共获得78.99 Gb 的有效数据，单个样品的有效数据均在12Gb 以上，Q20（有效数据质量不小于20 的碱基所占百分比）均在95%以上，Q30（有效数据质量不小于30 的碱基所占百分比）均在89%以上，该结果表明测序数据质量可靠，可用于进一步分析。

表3 贵州白山羊背最长肌总RNA 转录组测序数据情况

2.3 差异表达基因的筛选、功能注释及富集分析 FPKM是衡量转录组测序中基因表达水平的估算指标。在6个样本中共获得33 896 个转录本，筛选出FPKM ≥1的转录本有25 089 个。表4 列出在体重优势群体背最长肌中表达量前10 的候选基因。应用DESeq2 软件对筛选出的基因进行差异表达筛选，当设定参数FPKM ≥1 且P值小于0.05 时，在优势群体背最长肌中共筛选出563 个差异表达上调基因，其中94 个为新发现基因（转录本）。皮尔逊相关分析结果显示，2 组3个不同的重复样品的R2值基本都在0.85 以上（R2值越高，表明样品可靠性和可重复性越高），表明实验样品具有较高的可靠性（图1A），聚类分析热图显示2 个不同群体间差异表达基因重复性较好，差异基因分别聚类在一起（图1B）。

图1 贵州白山羊不同群体背最长肌样品及差异基因聚类分析

表4 优势体重群体背最长肌中表达量最高的10 个基因

2.4 差异表达基因GO 功能注释及KEGG 富集分析 GO注释系统可分为3 大类，分别为：生物学过程、细胞组分和分子功能。通过Goseq 软件对筛选出的差异表达上调基因进行GO 功能富集分析，其中36 个为显著富集（P＜0.05），这些基因在3 个分支中均有分布，包括9个生物学过程，20 个细胞组分和7 个分子功能（表5）。

表5 贵州白山羊背最长肌中上调表达基因GO 分析

应用Kobas 软件对差异表达基因进行KEGG 通路分析，结果表明注释到的差异表达上调基因共参与220条信号通路，其中FDR ≤0.05 的通路为表达基因中显著富集的信号通路，共有21 条，图2 列出了前20 个显著富集到的通路，最为富集的前5 条通路分别为代谢通路途径（含64 基因）、磷脂酰肌醇3-激酶信号通路（含25 个基因）、阿尔茨海默病（含24 个基因）、致癌的途径（含21 个基因）、帕金森病（含20 个基因）和非酒精性脂肪肝（含20 个基因）。

图2 贵州白山羊背最长肌中上调表达基因KEGG 通路富集散点图

2.5 差异表达基因qRT-PCR 验证 通过功能分析，随机选取4 个可能与贵州白山羊肌肉生长发育呈上调趋势的调控基因FGF11、NOTCH1、HBEGF和ADIPOQ进行qRT-PCR 验证，具体基因信息见表6。通过qRTPCR 验证，选取的4 个候选基因相对表达结果与转录组分析结果一致，说明测序数据准确可靠，结果见图3。

图3 候选基因在不同群体贵州白山羊背最长肌中的相对表达量

表6 可能与贵州白山羊背最长肌生长发育相关的候选上调基因

3 讨 论

贵州白山羊是在云贵高原特有的生态环境下，经过长期自然选择和人工选择而形成的地方品种，具有耐粗饲、性成熟早、繁殖力强、适应性强、肉质细嫩、无膻味等优点，但同时存在个体较小、生长速度慢、产肉量不高等缺点。近年来，由于自然和人为因素的影响，贵州白山羊个体生产性能差异较大，品种性能退化严重，并呈进一步加速趋势。因此，研究可用于辅助筛选贵州白山羊肌肉生长发育的候选基因及调控机制，对进一步提升贵州白山羊生产性能具有重要指导意义。肌肉生长发育是影响山羊生产性能的重要指标，本研究对同等饲养条件下、不同体重白山羊的背最长肌进行全转录组测序，筛选出优势体重群体相对于劣势体重群体的差异上调基因563 个，其中94 个为新发现基因（转录本），这些差异基因可能与贵州白山羊肌肉生长发育存在一定的关联性，值得进一步研究。

KEGG 富集分析得到了21 条富集的信号通路，最为富集的5 条通路中，PI3K-Akt 信号通路在细胞增殖、生长、代谢和生存中起着重要调控作用[10]，其中磷脂酰肌醇3 激酶I 类中的p85α亚基可以介导胰岛素的调控进而影响葡萄糖转运[11]。通路中蛋白激酶B（Akt/Pkb）重要的亚型Akt2 主要在骨骼肌、棕色脂肪中表达，而且研究发现，活化后的Akt 可以激活上百种底物，可以调控细胞生长、增殖、代谢、运动等功能[12-13]。吕欣等[14]也认为通过影响PI3K/Akt 信号通路中的蛋白表达、诱导通路的磷酸化，可以促进肌卫星细胞增殖分化、抑制凋亡，影响骨骼肌的再生修复。另外，筛选出基因最多的代谢途径信号通路中，通路信号众多，主要包括多聚糖生物合成代谢通路、脂类代谢通路、碳水化合物代谢通路、能量代谢通路、氨基酸代谢通路、核酸代谢通路、萜类化合物和多酮类化合物代谢通路等，均是参与能量代谢调控的通路。因此，本研究筛选出的通路可能在调解贵州白山羊肌肉生长发育中具有重要作用。

在筛选出的差异表达基因中，已有研究表明部分基因与肌肉的生长发育相关，本研究随机筛选4 个基因进行了验证。其中FGF11基因编码蛋白是一种自分泌型成纤维细胞因子，主要通过电压门控的钠离子通道相互作用调控神经系统的兴奋性[15]。已有研究表明，FGF11与毛细血管的结构稳定性及形成有关，同时还可能是肿瘤的重要调控因子[16]。Tejedor 等[17]最新研究发现FGF11基因敲除能够阻碍小鼠前肢的再生发育，表明FGF11基因可能是动物肢体发育的重要调控因子。本研究中发现，FGF11基因在体重优势群体中的表达量远高于体重劣势群体。NOTCH1基因编码蛋白是NOTCH信号通路中的重要受体，研究发现NOTCH 信号通路可以调解骨骼肌细胞的稳态、脂肪组织的相互转变、较少肥胖等[18]。亦有研究分析认为NOTCH1可以调解下游的PI3K、Akt 信号通路，而PI3K 和Akt 信号通路是参与调控肌肉细胞增殖、分化和发育的重要通路，同时该通路还参与葡萄糖的摄取与储存[19-20]。Yamaguchi 等[21]研究发现NOTCH1在小鼠体内的合成不足可以加速体内脂肪的生成。HBEGF基因编码蛋白是EGF 家族的糖蛋白成员，最早在巨噬细胞培养后的产物中获得，研究发现HBEGF可以促进平滑肌细胞、成纤维细胞和上皮细胞的增殖和迁移[22-23]。亦有研究表明，HBEGF基因是影响猪繁殖性状的重要候选基因[24-25]。本研究中发现HBEGF基因在体重优势群体中的表达量虽高于体重劣势群体，但差异不明显。因此对该基因的具体功能需进一步研究。ADIPOQ基因编码蛋白是一种对肥胖呈负相关调控的蛋白类激素，其主要参与动物机体内葡萄糖代谢和脂肪代谢。研究表明，ADIPOQ基因在脂肪组织中的表达量和在血液中的水平会显著升高，而当脂肪过度沉积时，在表达量和在血液中的水平均降低[26]。研究还发现ADIPOQ基因表达量与猪的胴体性状、肉质[27]，牛的眼肌面积和背膘厚度显著相关[28]，本研究团队也曾发现ADIPOQ基因与绵羊的生长速度和胴体性状相关[29]。罗宗刚等[30]研究发现ADIPOQ基因的表达在动物个体背最长肌中存在性别差异，表明ADIPOQ基因可能参与调解不同性别个体肌肉的生长调控。还有研究发现ADIPOQ基因在山羊脂肪组织中的表达量具有发育阶段的特异性[31]。在本研究中，不同群体间ADIPOQ基因的表达也存在显著差异，因此，ADIPOQ基因可能是肌肉生长发育的重要候选基因，其功能值得深入研究。

4 结 论

本研究基于转录组测序技术对不同体重的贵州白山羊背最长肌进行研究，在体重优势群体中发现563 个显著上调基因，其中94 个为新转录本，主要参与PI3KAkt 等信号通路的调控，筛选出的显著上调基因与转录组测序预测基因差异表达结果一致，表明差异基因结果可靠。因此，筛选出的上调基因可作为提升山羊群体肌肉生长性能的重要候选基因，值得进一步探讨其对贵州白山羊生长发育的调控机制。