不同月龄滩羊背最长肌组织的基因表达谱分析

马丽娜，蒋秋斐，马 青，额尔和花*

（1.宁夏农林科学院动物科学研究所，宁夏银川 750002；2.宁夏回族自治区畜牧工作站，宁夏银川 750001）

滩羊是主要分布于我国宁夏、甘肃等西部省份的优良裘皮用绵羊品种，除生产闻名于世的二毛皮之外，其肉质具有细嫩、不膻不腻等特点，受到消费者的广泛青睐。但滩羊与杜泊羊等国外优秀肉用绵羊品种相比，存在生长发育速度较慢的劣势，难以满足日益增长的消费需求。近年来，转录组学技术（RNA-seq）是解析家畜肌肉发育机理、筛选相关候选基因的常用技术手段，许多学者利用RNA-seq 技术对不同发育阶段的山羊背最长肌的基因表达谱进行了分析，并广泛应用于滩羊被毛卷曲、长脂尾等特定性状密切相关的候选基因及信号通路的研究，但是利用该技术筛选与滩羊肌肉发育相关候选基因的研究还相对较少。以CRISPR/Cas9为代表的基因编辑技术日趋成熟，显著加快了猪、牛、羊等家养动物新品种的精准培育进程。因此，筛选滩羊肌肉发育相关的候选基因，解析相关基因调控滩羊肌肉发育的分子机理，是利用基因编辑技术培育生长速度快、肉质优良的滩羊新品种的基础性工作。研究发现，34～36 kg 体重阶段为育肥滩羊的最佳屠宰体重，滩羊体重在32～34 kg 时羊肉品质最佳，通过优质滩羊胴体标准化生产技术的推广应用，出栏活重达到41～45 kg，9 月龄（体重33.89 kg±1.80 kg）滩羊肉嫩度、熟肉率、粗脂肪等指标最优。基于前期团队研究成果，本研究选取3 种不同体重阶段（体重4.83 kg±0.61 kg、34.71 kg±3.77 kg、体重42.62 kg±5.71 kg）的滩羊背最长肌组织，采用RNA-seq 技术对基因表达谱进行分析，旨在筛选到与滩羊肌肉发育相关的可靠候选基因用于后续功能验证，进而为通过基因编辑技术快速培育生长发育速度快的滩羊新品种提供备选基因。

1 材料和方法

1.1 实验动物 本实验选取9 只健康状况良好的盐池滩羊，按照出生日期和体重分为3 组：A 组（1 月龄，体重4.83 kg±0.61 kg）、B 组（6 月龄，体重34.71 kg±3.77 kg）和C 组（10 月龄，体重42.62 kg±5.71 kg），每组各3 只，且体重相近。该实验羊饲喂于宁夏朔牧盐池县滩羊有限公司（位于宁夏吴忠市盐池县），按照不同阶段滩羊相关饲养标准进行饲养。屠宰后采集背最长肌组织样品约20 g，经无菌PBS 冲洗2 次后，迅速分成若干小块，分装入冻存管中置于液氮中保存备用。

1.2 背最长肌组织总RNA 提取及建库测序 采用Trizol法提取滩羊背最长肌组织总RNA，具体操作步骤如下：把约100 mg 背最长肌组织放入已预冷的研钵中进行研磨，至组织研磨成粉末状转移至1.5 mL 离心管中；加入1 mL Trizol 试剂，颠倒混匀，室温孵育2 min；加氯仿0.2 mL，剧烈上下颠倒离心管若干次，充分混匀，室温下静置3 min；4℃ 12 000 r/min 离心15 min，吸取上层水相至另外一个新离心管中，加入等体积异丙醇，充分颠倒混匀，静置10 min；4℃ 12 000 r/min 离心15 min后小心弃掉上清液，加入无酶水配置的75% 乙醇，吹打若干次，4℃ 12 000 r/min 离心5 min；小心弃掉液体，室温下放置5 min，加入50 μL 无酶水溶解白色沉淀；最后，采用NanoDrop 分光光度计检测RNA 纯度和浓度，采用Agilent 2100 检测RNA 完整度，凝胶电泳检测RNA 降解度。

待所有样品检测合格后，各取2 μL RNA 用于建库测序。具体建库流程如下：使用带有Oligo（dT）的磁珠富集mRNA；加入片段化缓冲液将mRNA 打断成短片段；用六碱基随机引物以mRNA 为模板进行反转录合成一链cDNA，再加入缓冲液、dNTPs 和DNA 聚合酶I 合成第2 链cDNA；利用AMPure XP 磁珠纯化双链cDNA；对纯化后的双链cDNA 进行末端修复、加A和接头；通过AMPureXP 磁珠对双链cDNA 进行片段大小选择，选择大小在200～250 bp 的片段进行PCR 扩增以构建cDNA 文库；最后，文库质检合格后，采用Illumina Hiseq4000 测序平台PE150（paired-end）测序策略上机测序。

1.3 测序数据质控与参考基因组比对 将上述测序所得原始数据即Raw data（Raw Reads）通过北京奥维森基因科技有限公司开发的Perl 脚本进行数据过滤，具体内容如下：使用Trimmomatic（v0.33）软件去除带有测序接头的Reads、不确定碱基含量比例大于10% 的reads 以及低质量碱基即Q ≤20 含量大于50% 的reads，最终得到净数据即Clean Data（Clean Reads）。对所得净数据进行统计，将Q20>95%、Q30>89%（Q20、Q30：Phred 质量值大于20、30 的碱基数占总碱基数的比例）定义为测序质量好的数据，用于后续比对及数据分析。采用Tophat 2（v2.1.0）软件将上述质控所得净数据与最新版的绵羊参考基因组序列（序列号：GCF_002742125.1）进行比对。

1.4 转录本的组装及基因表达水平分析 采用Cufflinks（v2.1.1）软件对上述比对到参考基因组上的净读数进行组装，利用Cuffcompare 程序和已知的参考基因序列进行比较，以发现未注释的新基因或者已知基因新的外显子，同时可以优化已知基因的起始和终止位置。采用HTSeq（v0.5.4p3）软件对不同月龄滩羊背最长肌组织进行基因表达水平分析，以FPKM（Fragments Per Kilobase Of Exon Model Per Million Mapped Reads）即每百万Reads 中来自某一基因每千碱基长度的Fragments 数目以校正测序深度和基因长度对测序所得Reads 计数的影响，从而客观衡量基因表达水平。本实验选择以FPKM 值1 作为基因表达与否的标准，后续的分析只选择FPKM ≥1 的基因。

1.5 不同样本间的差异表达基因分析 采用上述步骤中基因表达水平分析所得到的FPKM 值来比较基因在不同组别样本间的表达，使用软件为DESeq（1.10.1），所得差异表达分析结果值使用Benjamini 和Hochberg方法控制FDR（False Discovery Rate）。差异基因筛选的标准一般为：值<0.05。

1.6 差异表达基因GO 和KEGG 富集分析 采用GOseq（v1.22）软件对不同组别背最长肌样本间的差异表达基因进行GO 富集分析，准确计算基因所富集分子功能（Molecular Function）、生物过程（Biological Process）和细胞 组成（Cellular Component）3 个GO分类条目及子条目的概率；采用KOBAS3.0对差异表达基因进行KEGG 通路分析，确定差异表达基因所参与的最主要生化途径和信号转导途径。

1.7 实时荧光定量PCR 验证结果 为验证测序结果的准确性，挑选和3 个在不同阶段表达逐步上调的基因和和3个在不同阶段表达逐步下调的基因进行验证，以-作为内参基因，采用实时荧光定量PCR（qRT-PCR）进行验证，引物信息见表1。试剂采用北京康为世纪公司生产的FastSYBR Mixture（CW0955S），反应体系为20 μL：上、下游引物（10 µmol/L）各0.8 µL，2×Fast SYBR Mixture 10 µL，模板cDNA（20 μg/μL）1 µL，ddHO 补足至20 µL。PCR 采取试剂说明书所推荐两步法：95℃预变性20 s；95℃变性3 s，60℃退火延伸30 s，共40 个循环，采集荧光；65～95℃进行溶解曲线分析。采用2法计算基因的相对表达量，以1 月龄样本作为对照组。

表1 引物信息

2 结果与分析

2.1 组织所提取总RNA 质量检测 由表2 可知，所有样本所提取总RNA 的OD和OD比值均在1.8～2.0 之间，RIN（RNA Integrity Number）值在7.5～8.4之间，符合RNA-seq 测序实验对RNA 质量的一般要求。上述结果说明所有样本提取RNA 质量良好，可用于下游的RNA-seq 实验。

表2 背最长肌组织RNA 提取质量检测

2.2 测序数据质量评估及参考基因组比对 由表3 可知，各个样本测序数据的原始读数在39 440 656～56 556 972 之间，原始碱基数在5.9～8.48G 之间；净读数在38 855 940～55 525 176 之间，净碱基数在5.82～8.32G 之间；测序碱基的错误率均为0.02%，Q20 值在98.09%～98.81%之间，Q30 值在94.11%～96.09%之间，GC 含量在50.81%～51.88%之间。上述各项指标均符合精准转录组测序建库的标准，可进行下一步分析。

表3 背最长肌组织原始测序数据质量评估

采用Tophat2（v2.1.0）软件将各个样本建库获得的净读数匹配到绵羊参考基因组上，结果发现（表4）：各个样本93.5%～94.59% 的净读数都能匹配到参考基因组上；在匹配到参考基因组上的净读数中，72.08%～76.94%的读数匹配到基因外显子区域，7.20%～10.02%的读数匹配到基因内含子区域，而15.63%～18.07% 的读数匹配到基因间区域。该现象符合匹配到基因外显子区域的读数在匹配读数占比最高的正常情况，同时匹配到内含子区域的读数来源于前体mRNA 的残留及在前体mRNA 发生可变剪切形成成熟mRNA 过程中的内含子保留事件，而定位到基因间区域的读数可能是由于目前参考基因组注释不完备引起的。

表4 背最长肌组织净读数的参考基因组比对情况

2.3 基因表达定量分析 在对匹配到参考基因组上的净读数进行组装后，采用FPKM 值表示基因的表达量并对所有基因的表达量进行统计分析，结果如图1 所示。可见，不同组别样本小提琴图具有类似的形状，表明不同组别样本的基因表达在整体上高度类似。

图1 不同样本基因表达定量分析

2.4 差异表达基因筛选 由图2 可知，6 月龄组和1 月龄组相比，差异表达基因共664 个，上调基因236 个，下调基因428 个；10 月龄组和1 月龄组相比，差异表达基因共2 131 个，上调基因927 个，下调基因1 204个；10 月龄组和6 月龄组相比，差异表达基因共442个，表达上调基因188 个，表达下调基因254 个。对不同比较组之间上调的基因取交集发现，共有的上调基因10 个，包括等已知促进不同物种肌肉发育相关的基因；共有的下调基因16 个，包括等已知的抑制肌肉发育相关的基因。相关基因在不同组别样本间的表达水平（以FPKM 表示）见表5。

表5 3 个比较组别间共有的差异表达基因

图2 不同组别差异表达基因火山图

2.5 差异表达基因的GO 分析 采用GOseq 软件对不同组别间差异表达基因进行GO 富集分析，将同基因归类到GO 三大分类及所属子条目中，结果如图4 所示。可知，1 月龄组和6 月龄组差异基因所富集到的生物过程包括刺激响应、免疫响应和免疫系统过程等，富集到的细胞组分包括外切脱氧核糖核酸酶7 复合物等，而富集到的分子功能包括核苷结合、细胞骨架蛋白结合和GTP 结合等；1 月龄组和10 月龄组差异基因所富集到的生物过程包括生物合成过程、小分子代谢过程和含氧酸代谢过程等，富集到的细胞组分包括胶原蛋白三聚物，而富集到的分子功能包括催化活性、氧化还原酶活性和焦磷酸酶活性等；6 月龄组和10 月龄组差异基因所富集到的生物过程包括免疫响应、免疫系统过程和多细胞生物发育等，细胞组分包括内质网、三磷酸甘油脱氢酶复合物和RNA 聚合酶III 复合物等，而分子功能包括蛋白结合、信号受体结合和信号受体激活子活性等。

图3 不同组别差异表达基因GO 分析

图4 不同组别差异基因KEGG 分析

2.6 差异表达基因的KEGG 分析 由图4 可知，1 月龄组和6 月龄组差异基因共富集到263 条KEGG 通路，其中前20 条显著富集的信号通路包括病毒性心肌炎、心肌细胞肾上腺激素信号通路、心肌收缩和扩张性心肌病等；1 月龄组和10 月龄组差异基因共富集到315 条KEGG 通路，其中前20 条显著富集的信号通路包括胰岛素分泌、胰岛素抗性、胰高血糖素信号通路、ECM-受体互作、心肌细胞肾上腺激素信号通路和AMPK 信号通路等；6 月龄组和10 月龄组差异基因共富集到233条KEGG 通路，其中前20 条显著富集的信号通路包括癌症中的转录失调、MAPK 信号通路、FoxO 信号通路和细胞凋亡等。

2.7 测序结果的qRT-PCR 验证 如图5 所示，通过qRT-PCR 检 测和基因在不同月龄组织样本间的表达趋势同RNA-seq 完全一致，而和基因的表达趋势同RNA-seq 略有差异，具体表现为6 月龄组和10 月龄组表达量均显著低于或者高于1 月龄组，同RNA-seq 结果相符，但前2 组之间基因表达趋势同RNA-seq 结果不同。上述结果表明，和可以作为候选基因直接进行下一步的功能验证，而和还需要进一步的分析。

图5 qRT-PCR 验证测序结果

3 讨论

本研究利用RNA-seq 技术构建了不同月龄滩羊背最长肌组织的基因表达谱并通过生信分析，成功筛选到了甲基转移酶样1（Methyltransferase Like 1，METTL1）、泛素特异性蛋白酶43（Ubiquitin Specific Peptidase 43，USP43）等调控滩羊肌肉发育的候选基因。现有研究证实：动物在出生后主要依赖肌纤维体积的增加而实现肌肉的生长，而该过程与肌卫星细胞密切相关。肌卫星细胞是肌肉中通常处于静息状态的肌源性干细胞，在接受其所附属肌纤维、周围微环境等来源的信号后被激活，持续增殖和分化形成肌细胞、肌管，同已有肌纤维融合，最终增加肌纤维的体积和长度。因此，推测本研究中筛选到的候选基因和信号通路可能同滩羊出生后肌肉中肌卫星细胞激活、增殖和分化密切相关。

本研究发现，等基因在出生后滩羊背最长肌组织的表达量逐步上升，表明它们是可能促进滩羊肌卫星细胞增殖和分化的关键基因。基因已被发现是维持小鼠胚胎干细胞自我更新的必需基因，而自我更新是肌卫星细胞通过细胞分裂产生相同细胞以维持干细胞池中细胞数目稳定的主要方式。芳香烃受体核转位因子样蛋白（）基因亦被发现通过刺激小鼠肌卫星细胞的增殖和分化而参与肌肉再生过程，表明该基因是动物出生后肌肉发育的重要参与者，因此可能通过刺激出生后滩羊肌肉中肌卫星细胞的增殖和分化而促进滩羊肌肉组织的发育。在不同阶段滩羊背最长肌肌肉组织中表达量较高且也呈逐步上升的趋势，其已被证实参与鸡肌卫星细胞分化成为肌管的过程而不影响肌卫星细胞增殖的过程。同时，Miao 等发现在不同生长发育速度的陶赛特和小尾寒羊背最长肌组织表达量存在显著差异，暗示该基因可能同绵羊肌肉组织生长发育调控有关。此外，基因的单核苷酸多态性也被发现同猪、牛等家畜肉质性状相关联。已经被发现同人肌卫星细胞的晚期增殖和早期分化有关，且在处于静息期的肌卫星细胞中不表达，说明的表达与肌卫星细胞由静息期向增殖分化状态转变相关。同时，也有研究证实与小鼠锻炼诱导的肌卫星细胞增殖和分化导致的骨骼肌肉再生有关。结合上述基因在出生后滩羊背最长肌中的表达模式和在其他动物上的已有研究推测，上述基因可能通过调节肌卫星细胞的自我更新、增殖和分化等来促进出生后滩羊肌肉组织的发育。

同时，本研究也发现一些在滩羊肌肉生长发育过程中表达量逐步下调的基因，说明这些基因可能在滩羊肌肉生长发育过程中起着负向调控的作用。其中，在小鼠中的过表达会导致肌肉发育迟缓，而这种现象同维持肌卫星细胞处于静息状态而抑制其增殖分化有关。也是参与骨骼肌能量代谢的基因之一，其被发现特异性地表达于猪骨骼肌中且在肌肉发育过程中表达水平持续上调，同时在约克夏猪背最长肌组织中表达量高于梅山猪，而且基因单核苷酸多态性同眼肌面积等肉质性状相关。但是本研究发现该基因在滩羊肌肉组织发育过程中表达量持续下调，暗示该基因可能是滩羊肌肉发育的负调控因子，这种不一致可能与物种特异性有关。考虑到与肌肉发育和肉质的密切相关性，未来的工作需要进一步探索其在滩羊肌肉组织发育的功能和相应的分子机理。

对差异表达的基因进行GO 和KEGG 功能富集发现，免疫响应、胰岛素分泌和MAPK 等生物学过程和信号通路都得到了显著富集。免疫细胞是肌肉的重要组成部分，它们可以通过分泌TNF-等因子调节肌卫星细胞增殖和分化的过程。胰岛素可以刺激绵羊肌卫星细胞的增殖，而MAPK 信号通路活化与静息期肌卫星细胞的激活密切相关。上述信号通路的显著富集不仅说明它们在滩羊肌肉组织生长发育过程中也起着至关重要的调控作用，而且也暗示调节相关信号通路可能是促进滩羊肌卫星细胞增殖和分化、加速肌肉发育的潜在手段。

4 结论

本研究通过RNA-seq 技术对1 月龄、6 月龄和10月龄滩羊背最长肌组织进行基因表达谱分析，通过对不同阶段差异表达基因取交集，成功筛选到了等10 个促进滩羊肌肉发育的候选基因及等16个起负向调控作用的基因；通过对差异表达基因进行功能富集，发现免疫响应、胰岛素分泌和MAPK 等信号通路是参与滩羊肌肉生长发育的关键生物学过程和信号通路。上述候选基因和信号通路的成功筛选为通过基因编辑等技术手段培育生长发育速度快的滩羊新品种提供了可靠信息。