基于RNA-Seq转录组测序技术揭示南方黄牛肉质嫩度调控相关的分子机制

孟祥忍，张成龙，樊永亮，王文强，毛永江，黄必志，亏开兴，杨章平⋆（.扬州大学旅游烹饪学院，江苏扬州 57；.云南省草地与动物研究院，云南昆明650）

基于RNA-Seq转录组测序技术揭示南方黄牛肉质嫩度调控相关的分子机制

孟祥忍1，张成龙1，樊永亮1，王文强1，毛永江1，黄必志2，亏开兴2，杨章平1⋆（1.扬州大学旅游烹饪学院，江苏扬州 225127；2.云南省草地与动物研究院，云南昆明650212）

为揭示南方黄牛肉质调控的分子机制并挖掘与肌肉嫩度相关的功能基因，本研究以1周岁左右的云岭牛、文山牛和中国西门塔尔牛为试验对象，屠宰后测定背最长肌不同排酸时间（0 、1 、2 、3 、5 、7 d）下剪切力变化，揭示不同黄牛品种肌肉嫩度的差异，利用RNA-Seq技术对背最长肌进行比较转录组学分析，并结合生物信息学筛选与肌肉嫩度调控相关的信号通路及差异表达基因，最后利用qPCR进行验证。结果表明：在宰后成熟过程中，黄牛肌肉嫩度表现为明显的下降趋势，且在同一个成熟时间点下，剪切力值表现为西门塔尔牛＞文山牛＞云岭牛，表明云岭牛嫩度较好，而西门塔尔牛嫩度相对较差；云岭牛与西门塔尔牛背最长肌2组文库经测序后，筛选出3 198个差异表达基因，其中云岭牛上调基因1 750个，KEGG功能注释发现2个与肉质调控相关的重要信号通路：脂肪细胞因子信号通路和内质网蛋白加工通路，结合qPCR验证、基因表达水平与剪切力值相关性分析，最终筛选出Leptin和CAPN1 2个可能肌肉嫩度相关的功能基因。本研究筛选与鉴定出南方黄牛肌肉嫩度相关信号通路及其关键基因，为今后南方黄牛肉质遗传改良奠定理论基础。

南方黄牛；嫩度；转录组

嫩度是指肉在切割时所需的剪切力，剪切力是指测试仪器的刀具切断被测肉样时所用的力。通过测定仪器测传感器记录刀具切割肉样时的用力情况并把测定的剪切力峰值作为肉样嫩度值，结果单位以N或者kg之间切换。肉质嫩度是评价畜禽肉品质优劣的重要指标，反映了肉中各种蛋白质的结构特性，肌肉中脂肪的分布状态及肌纤维中脂肪数量等，其可以在机体宰后成熟过程中进行测定，但测定过程中存在着系统误差，且在活体上难以直接度量，迫切需要从遗传本质上寻找影响肌肉嫩度的功能基因，并利用分子标记辅助选择技术进行肉质性状改良。目前，关于肌肉嫩度相关候选基因的研究

主要集中在鸡和猪上[1-2]，对于肉牛特别是南方黄牛调控肌肉嫩度的分子机制有待进一步研究。随着高通量测序技术的发展，转录组RNA-seq测序已经广泛用于研究不同物种或者同一物种不同时期特定组织的基因表达水平差异，通过利用生物信息学（Bioinformatics）技术来揭示特定性状形成的分子机制，对今后物种优良性状的选种选育奠定理论基础。因此，本研究以云岭牛、文山牛和中国西门塔尔牛作为试验对象，通过检测不同排酸时间（0、1、2、3、5、7 d）肉质嫩度，筛选嫩度存在差异的黄牛品种进行转录组测序，并结合生物信息学分析和qPCR验证筛选肉质嫩度相关的重要调控通路和功能基因，不仅为南方黄牛肉质改良育种提供理论基础，同时也可为今后高档牛肉的生产提供重要的参考材料。

1 材料与方法

1.1 试验动物 为随机选取1周岁左右云岭牛、文山牛和中国西门塔尔牛（均来自云南省草地动物研究院），其中每个品种各5头，编号西门塔尔牛（S01、S02、S04、S03、S05）， 文 山 牛（W01、W02、W03、W04、W05）和 云 岭 牛（10513、10577、10556、10564、10568），屠宰后其中一方面取背最长肌组织立刻放入液氮中，带回实验室放入-70℃保存备用，其次将肉牛酮体放入4℃排酸库中进行6 d的排酸处理，按照不同排酸时间（0、1、2、3、5、7 d）中，测定背最长肌嫩度变化情况。

1.2 肌肉嫩度测定 采用美国F.T.C公司的TMSPro物性测试仪测定剪切力值大小，每个肉样的剪切力值以平均值±标准差表示。测量参数：探头：HDP/BSW；PPS：200.00；测试模式与选择：Measure Force in Confrerensing；测前速度：120.0 mm/min；测定速度：60.0 mm/min；测后速度：60.0 mm/min；距离20 mm；起始力0.7 N。

1.3 转录组测序及数据分析

1.3.1 RNA提取、cDNA文库建立以及上机测序 使用Trizol法提取云岭牛和中国西门塔尔牛背最长肌样品总RNA（每个3份），进一步利用Nanodrop检测RNA的纯度（OD260/280），Qubit对RNA浓度进行精确定量，Agilent 2100精确检测RNA的完整性。随后使用DNase消化DNA后，用带有Oligo（dT）的磁珠富集真核生物mRNA；加入打断试剂将mRNA打断成约200 bp短片段，以打断后的mRNA为模板，用六碱基随机引物合成一链cDNA，然后加入缓冲液、dNTPs（dTTP、dATP、dGTP和dCTP）和DNA polymeraseI合成二链cDNA，并使用试剂盒AMPure XP beads纯化双链cDNA；纯化的双链cDNA再进行末端修复、加A尾并连接测序接头，然后进行片段大小选择，最后进行PCR扩增；构建好的文库用Agilent 2100 Bioanalyzer质检合格后，使用Illumina HiSeqTM2500进行测序。

1.3.2 测序数据质量控制 由Illumina HiSeqTM2500测序所得的数据称为raw reads或raw data，随后要对raw reads进行质控（QC），以确定测序数据是否适用于后续分析。质控后，经过滤得到clean reads，用tophat/bowtie2将clean reads比对到参考序列。比对完，通过统计reads在参考序列上的分布情况、覆盖度以及Reads在不同元件的富集分析，同时利用ASprofile[3]软件分析其可变剪接事件。

1.3.3 差异表达基因筛选及聚类分析 利用RNA-seq数据比较分析2个样品中同一个Unigene是否存在差异表达时，可以选取2个标准：一是FoldChange，就是两样品中同一个Unigene表达水平的变化倍数；二是pvalue或FDR（padjust），FDR值的计算方法先要对每个Unigene进行P-value的计算，再用FDR错误控制法对P-value作多重假设检验校正。默认筛选差异的条件为P≤0.05。此 外 使 用 Cluster 3.0软 件（http://bonsai.hgc. jp/～mdehoon/software/cluster/software.htm）对不同差异表达基因RPKM值进行聚类分析。

1.3.4 GO功能和Pathway通路富集与注释 GO（Gene Ontology, http://geneontology.org/） 挖 掘差异表达的 mRNA 的主要生物学功能，包括分子功能（Molecular Function）、细胞组成（Cellular Component）、生物学过程 （Biological Process），分析软件：GOSeq[4]；KEGG（Kyoto Encyclopedia of Genes and Genomes，http://www.kegg.jp/），又称京都基因与基因组百科全书，是基因组破译方面的数据库。KEGG注释的主要目的包括2个：KO（KEGG Ortholog）注释，即将分子网络的相关信息进行跨物种注释；KEGG Pathway注释，即代谢通路注释，获得物种内分子间相互作用和反应的网络，分析软件：KOBAS[5]。

1.4 引物设计 根据GenBank数据库已公布的牛MSTN、Leptin、LEPR、STAT3、HSPA1A、CAPN1和内参基因β-actin基因序列，利用Premier 6.0软件进行跨外显子引物设计，引物由上海生工生物工程技术服务公司合成，具体引物信息如表1。

1.5 差异表达基因的qPCR验证 利用Trizol法提取各组织总RNA，以RNA为模板进行cDNA合成。cDNA合成：10 μL的反应体系中含5×Primer-Script Buffer反应液2 μL，PrimerScript RT Enzyme MixI0.5 μL，Oligo dT 0.5 μL，Random 6 mers 0.5 μL，总RNA 500 ng， RNase free H2O补足至10 μL。反应条件为37℃ 15 min，85℃ 5 s，4℃保存。荧光染料为SYBR GreenI（TaKaRa），反应体系为20 μL：cDNA 1 μL，上、下游引物各0.4 μL（10 μmol /L），ROX Reference DyeⅡ（50×）0.4 μL，SYBR Green Real-time PCR Master Mix（2×）10 μL，ddH2O 7.8 μL。Real-time PCR反应条件：95℃预变性15 s，95℃ 5 s，60℃ 34 s，40个循环，扩增结束后分析熔解曲线。

表1 引物序列信息

1.6 统计分析 相对定量的结果采用2-ΔΔCt法[6]进行处理，计算差异表达基因在云岭牛与西门塔尔牛背最长肌组织中的表达水平，并利用SPSS18.0软件对云岭牛与西门塔尔牛背最长肌中基因表达水平差异程度进行卡方独立性t检验，同时基于Pearson方法对基因表达水平与剪切值进行相关性分析。

2 结果与分析

2.1 排酸过程中肌肉嫩度变化情况 由图1可以看出，云岭牛宰后0～5 d云岭牛、西门塔尔牛和文山牛剪切力值均表现为明显的下降趋势，5 d后剪切力值下降幅度不大，趋于稳定；整体来看，同一个排酸时间点的剪切力值表现为西门塔尔牛＞文山牛＞云岭牛，表明云岭牛嫩度较好，而西门塔尔牛嫩度相对较差。

2.2 转录组测序分析 为了进一步探讨黄牛肌肉嫩度调控的分子机制，本研究选取云岭牛（嫩度较好）和西门塔尔牛（嫩度较差）2种肉牛背最长肌作为样本，云岭牛编号为Y1～Y3，西门塔尔牛编号为S1～S3，进行转录组测序分析。

2.2.1 样本RNA质量检测结果 由表2可见，经过质量检测后，选取6个背最长肌样本（S1～S3和Y1～Y3）RNA总量、RNA完整性指数（RIN）以及纯度均符合测序标准（RIN≥8.0；OD260/ OD280≥2.0；28S/18S≥1.6），可以用于下一步的RNA-Seq转录组测序。

图1 不同黄牛品种肌肉排酸过程中肌肉嫩度变化

2.2.2 测序数据质量检测结果 通过碱基质量分布分析对每个样本进行测序错误率分布检查。由表3可以看出，碱基有效比率Valid ratio和Q30均在90%以上；由碱基质量分布图（图2和图3）可以看出，测序数据质量都在平均值以上，Reads匹配到不同的物种，并且测序结果获得的转录本除了前几个碱基具有较大波动，随后所有碱基处于水平线的稳定趋势，基本没有AT和GC分离现象。此外，绝大多数reads比对到了外显子区域CDS_Exon（图4），对测序获得的数据进行Map比对到牛参考基因组上（表2、3），匹配到参考基因组的Reads数在85.43%～90.49%，综合表明测序数据非常良好，结果可利用率较高，可以用于下一步分析。

表2 样本RNA质量检测

2.2.3 基因结构分析—可变剪切事件预测 由表4可知，可变剪接分析结果（图5）表明，共发现12类型的可变剪切事件，其中主要可变剪切为可变末尾外显子（Alternative Last Exons，TTS）、可变第1外显子（Alternative First Exons，TSS）、单内含子保留（Intron Retention，IR）和可变5'或3'端剪切（Alternative Exon Ends，AE）。

表3 测序数据质量评估

图2 原始数据碱基质量分布图

图3 碱基分布图

图4 样本reads在不同原件上分布图

表4 与参考基因组比对统计结果

图5 可变剪切分析

图6 差异表达基因散点图

2.2.4 差异表达基因及聚类分析 云岭牛与西门塔尔牛背最长肌两组文库经测序后，利用P＜0.05和log2fold change＞2原则在两组之间共筛选出3 198个差异表达基因，其中云岭牛组上调基因1 750个（图6），进一步对差异表达基因进行聚类分析（图7），发现云岭牛与西门塔尔牛样本中高低表达水平的基因聚类在一起。

图7 差异表达基因聚类分析

图8 KEGG-pathway通路注释

2.2.5 KEGG通路富集分析 KEGG通路富集分析发现（图8），大部分差异表达基因主要在20种通路中富集，其中主要在脂肪细胞因子信号通路（Adipocytokine Signaling Pathway） 和 内 质 网蛋白加工通路（Protein Processing in Endoplasmic Reticulum）中高度富集，并且筛选出MSTN、Leptin、LEPR、STAT3、HSPA1A、CAPN1等几个可能与肌肉嫩度相关的差异表达基因。

图9 差异表达基因qPCR验证

图10 基因表达水平与剪切力值相关性分析

2.3 重要差异基因表达水平验证及其与肌肉剪切力值相关性分析 对6个重要差异表达基因（MSTN、Leptin、LEPR、STAT3、HSPA1A、CAPN1），利用qPCR在西门塔尔牛和云岭牛背长肌组织中进行表达检测（每组3个重复，见图9），结果表明该6个基因在两组之间表达水平均呈现显著或极显著差异（P＜0.05或P＜0.01），其表达趋势与转录组RNA-seq结果基本一致，表明测序结果可靠。此外，进一步分析以上差异基因表达水平与肌肉剪切力值相关性表明（图10），瘦素基因（Leptin）和钙蛋白酶1基因（CAPN1）分别与剪切力值表现为显著负相关（P＜0.05），相关系数分别为0.898和0.836，而其他基因无显著相关性（P＞0.05）。

3 讨 论

我国南方黄牛遗传资源分布具有品种丰富、起源独特、存栏量大以及环境适应性强等特点，然而其牛肉品质与牛肉的经济价值密切相关。本研究从肌肉嫩度入手，分析了3种不同黄牛（云岭牛、文山牛和中国西门塔尔牛）在宰后成熟过程中肌肉嫩度表现为下降趋势，由于屠宰后肌肉中仍然进行各种生化反应，但是主要进行无氧代谢，不断产生大量的乳酸等物质，从而引起肌肉中pH降低。无氧代谢使ATP合成受阻导致其浓度下降，进而使肌浆网钙泵功能失去，导致肌肉收缩连续但不可逆，收缩到最大程度时形成肌肉的宰后僵直。研究发现，肌肉嫩度在一开始屠宰后最佳，此后随着pH变化导致嫩度下降[7]。此外，本研究还比较了这3种不同黄牛肌肉嫩度变化，发现云岭牛肌肉嫩度较好，而西门塔尔牛相对较差。在此基础上，本研究利用转录组测序技术分析不同南方黄牛品种之间背最长肌差异表达基因，发现了2个重要的调控通路：脂肪细胞因子信号通路和内质网蛋白加工通路，并且筛选出MSTN、Leptin、LEPR，STAT3，HSPA1A、CAPN1等几个重要差异表达基因。肌肉生长抑制基因（Myostatin，MSTN）属于转化生长因子TGF-β家族成员之一，主要负责骨骼肌生长发育[8]，有关学者研究发现MSTN缺陷型小鼠的肌肉变大，肌群分布更加广泛，其质量也增加2倍以上[9]，同时也有研究发现牛MSTN基因多态性与双肌性状密切相关[10]。瘦素基因（Leptin）是机体脂肪细胞分泌的蛋白类激素，通过负反馈调节将肥胖信号传递给大脑，从而减少采食量，提高机体能量代谢，可以降低体重，因此Leptin在调控脂肪代谢和体重稳态中起着关键作用，同时也有研究发现Leptin基因多态性与肉质性质有关[11]。此外，瘦素受体（LEPR）是瘦素（leptin）发挥作用的关键受体物质。信号转导子与转录激活子STAT3，主要参与细胞增殖、分化与迁移等功能，在机体生长发育过程中起着重要的调控作用[12]。钙蛋白酶1（CAPN1）基因作为一种半胱氨酸疏基内肽水解酶，主要调控肌肉生长发育、细胞分化凋亡以及信号转导等生物学过程[13]。研究发现，钙蛋白酶是降解肌肉纤维的主要酶[14]，并且酶活性与肉质嫩度的变化密切相关[15]，此外也有研究发现钙蛋白酶能有效地提高牛对饲料的利用率。结合以上基因功能的相关报道，本研究分析了基因表达水平与肌肉嫩度（剪切力）的关系，发现Leptin和CAPN1均与剪切力值表现出较高的相关性，因此，本研究综合推测瘦素基因（Leptin）和钙蛋白酶1（CAPN1）基因可能在南方黄牛肉质（肌肉嫩度）形成过程中起着重要调控作用，今后需要在细胞水平上利用RNA干扰、过表达以及TALEN敲除技术来进一步验证2个关键基因的功能，此外个体水平进行遗传变异位点的检测，重点筛选影响基因表达水平的变异位点，分析其多态性与肉质性状的关系，旨在寻找南方黄牛肉质调控相关的遗传标记，为今后南方黄牛遗传改良与优质育种提高一定的理论参考。

4 结 论

3种南方黄牛（西门塔尔牛、云岭牛、文山牛）在宰后排酸过程中，嫩度表现为明显的下降趋势，并且在相同排酸时间点下，剪切力值表现为西门塔尔牛＞文山牛＞云岭牛，表明云岭牛嫩度较好，而西门塔尔牛嫩度相对较差。本研究以西门塔尔牛和云岭牛为对象，进一步利用比较转录组学揭示了南方黄牛肉质调控的分子机制，并发现脂肪细胞因子信号通路和内质网蛋白加工通路2个信号通路与肉质调控相关，筛选出Leptin和CAPN1 2个关键功能候选基因。

[1] 张瑜, 王昕陟. 钙蛋白酶系统对猪肉嫩度影响的研究进展[J]. 中国畜牧兽医, 2013, 40(10): 161‐164.

[2] 陆克龙. 肉鸡CAPN1和CAST基因遗传多态性及其与肌肉嫩度相关性[D]. 南京: 南京农业大学, 2012.

[3] Florea L, Song L, Salzberg S L. Thousands of exon skipping events differentiate among splicing patterns in sixteen human tissues[J]. F1000 Research, 2013, 2: 1‐24.

[4] Young M D, Wakefield, M J, Smyth G K, et al. Gene ontology analysis for RNA‐seq: accounting for selection bias[J]. Genome Biol, 2010, 11: R14.

[5] Mao X, Cai T, Olyarchuk J G, et al. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary[J]. Bioinformatics, 2005, 21(19): 3787‐3793.

[6] Livak K J, Schmittgen T D. Analysis of relative gene expression data using real‐time quantitative PCR and 2−ΔΔCTmethod[J]. Methods, 2001, 25:402‐408.

[7] 周光宏. 肉品学[M]. 北京: 中国农业科技出版社, 1999: 342‐344.

[8] McPherron A C, Lee S J. Double muscling in cattle due to mutations in the myostatin gene[J]. Proc Natl Acad Sci U S A, 1997, 94(23): 12457‐12461.

[9] Lee S J, Aleandra C. Regulation of myostatin activity and muscle growth[J]. PNAS, 2001, 98(16): 9306‐9311.

[10] Wiener P, Woolliams J A, Frank‐Lawale A, et al. The ef f ects of a mutation in the myostatin gene on meat and carcass quality[J]. Meat Sci, 2009, 83(1): 127‐134.

[11] Schenkel F S, Miller S P, Ye X, et al. Association of single nucleotide polymorphisms in the Leptin gene with carcass and meat quality traits of beef cattle[J]. J Anim Sci, 2005, 83(9): 2009‐2020.

[12] Herrington J, Smit LS, Schwartz J, et al. The role of STAT proteins in growth hormone signaling[J]. Oncogene, 2000, 19(21): 2585‐2597.

[13] Gitler D, Spira M E. Real time imaging of calcium‐induced localized proteolytic activity after axotomy and its relation to growth cone formation[J]. Neuron, 1998, 20(6): 1123‐1135.

[14] Koohmaraie M. Biochemical factors regulating the toughening and tenderization process of meat[J]. Meat Sci, 1996, 43: S193‐S201.

[15] Geesink G H, Koohmaraie M. Effect of calpastatin on degradation of myofibrillar proteins by calpain under postmortem conditions[J]. J Anim Sci, 1999, 77: 2685‐2692.

Molecular Mechanism Analysis of Meat Tenderness in Chinese South Cattle Based on RNA-Seq Transcriptome Sequencing Technology

MENG Xiang‐ren1, ZHANG Cheng‐long1, FAN Yong‐liang1, WANG Wen‐qiang1, MAO Yong‐jiang1, HUANG Bi‐zhi2, QU Kai‐xing2, YANG Zhang‐ping1*

( 1. Animal Science and Technology College, Yangzhou university, Jiangsu Yangzhou 225009, China; 2.Grassland and Animal Research Institute, Yunnan KunMing 650000, China)

To reveal the molecular mechanism of meat quality and excavate function genes related to tenderness in Chinese south cattle, this study selected Yunling, Wenshan and Simmental cattle at approximately 12 months of age as experimental object. Firstly, after slaughter, muscle tenderness of above cattle was measured at dif f erent stages of acid discharge (0, 1, 2, 3, 5, 7d), which aimed at revealing the dif f erences of meat quality among three varieties. Using transcriptome sequencing and bioinformatics analysis, we screened out dif f erential expression genes related to meat quality, then real‐time fl uorescence quantification PCR was further performed. Results showed, meat tenderness of three varieties (Yunling, Wenshan and Simmental) showed a decreasing trend at dif f erent stages of acid discharge. At the same stage of acid discharge, the value

of shear force showed the highest in Simmental, followed by Wenshan and Yunling, which indicated that the tenderness of Yunling was relatively good and Simmental was relatively poor. After cDNA library sequencing of longissimus dorsi between Yunling and Simmental, we screened out 3198 dif f erential expression genes, of which 1750 genes were up‐regulated in Yunling. Additionally, we found two important pathways related to meat quality, such as adipocytokine signaling pathway and protein processing in endoplasmic reticulum. Further qPCR verif i cation and correlation analysis between gene expression level and shear force value, we ultimately identif i ed Leptin and CAPN1 genes related to meat tenderness. This study screened and identif i ed important pathways and key genes regulating meat tenderness, which will provide theoretical basis for genetic improvement of meat quality in Chinese south cattle.

Chinese south cattle; Tenderness; Transcriptome

S823.2

A

10.19556/j.0258-7033.2017-09-026

2017-05-20；

2017-07-06

863国家高技术研究发展计划项目（2013AA10250 5-5）

孟祥忍（1977-），男，博士，副教授，主要从事烹饪工艺学研究，E-mail: 455609455@ qq. com

*通讯作者：杨章平（1965-），男，博士，教授，主要从事动物遗传资源评价、保护与利用的研究，E-mail: yzp@ yzu.edu.cn