不同品种肉牛肌内脂肪沉积相关基因比较研究

施 安 孙文阳 陈志龙 李 博 侯鹏霞* 梁小军*

(1.宁夏农林科学院动物科学研究所,银川 750002;2.宁夏农林科学院固原分院,固原 756099;3.西北农林科技大学动物科技学院,杨凌 712100)

西门塔尔(Simmental)和安格斯(Angus)是我国具有优势特色和发展潜力的主要肉牛品种。研究发现,安格斯牛肉有呈现低饱和脂肪酸(saturated fatty acid,SFA)、高亚麻酸(C18∶3n3)和不饱和脂肪酸(unsaturated fatty acid,UFA)含量的趋势,并且其肌内脂肪(intramuscular fat,IMF)含量较高[1-2]。此外,安格斯牛肉在多汁、嫩度和风味上均优于中国西门塔尔牛肉[3-4],是打造高端牛肉市场上的主推品种。

IMF直接影响牛肉的品质和风味[5-6],其由肌原性干细胞和前脂肪细胞分化转移形成[7],也可以从成纤维细胞样前脂肪细胞到成熟脂肪细胞的多步骤分化形成[8];其中,前脂肪细胞还通过旁分泌因子的细胞间通讯、细胞外基质(extracellular matrix,ECM)和细胞-细胞黏附调节与周围结缔组织的增殖、分化和脂肪生成[9]。前脂肪细胞向成熟脂肪细胞的转化是由一系列转录因子控制的,其中过氧化物酶体增殖物激活受体γ(peroxisome proliferator-activated receptor γ,PPARγ)和CCAAT/增强子结合蛋白α(CCAAT/enhancer binding protein α,C/EBPα)被认为是关键的成脂因子,它们共同发挥多效转录激活因子的作用,包括脂肪酸结合蛋白4(fatty acid binding protein 4,FABP4)、脂蛋白脂肪酶(lipoprotein lipase,LPL)、脂肪酸转位酶(fatty acid translocase,FAT)和产生脂肪细胞表型的周围脂素[10-11]。在PPARγ和C/EBPα诱导之前,还有其他转录因子驱动成脂程序,包括CCAAT/增强子结合蛋白β(CCAAT/enhancer binding protein β,C/EBPβ)、CCAAT/增强子结合蛋白δ(CCAAT/enhancer binding protein δ,C/EBPδ)、Krüpple样转录因子(Krüppel-like factor,KLF)和成纤维细胞生长因子(fibroblast growth factor,FGF)等早期成脂因子[12-13]。

转录组测序(transcriptome sequencing,RNA-Seq)可精确获得某一物种在特定条件下所有转录产物的表达状况。利用高通量测序计算mRNA表达量,结合生物信息分析可高效获得性状相关基因的表达模式与功能差异[14-15]。Zheng等[3]对安格斯牛和中国西门塔尔牛的粪便和背最长肌进行多组学分析发现,安格斯牛的肠道微生物多样性明显高于中国西门塔尔牛,中国西门塔尔牛蛋白上调的差异表达基因(differentially expressed genes,DEGs)主要参与免疫和炎症反应,而安格斯牛蛋白上调的DEGs主要参与肌肉系统和肌原纤维的调节,且该研究鉴定出了17个与肌肉和脂肪代谢相关的基因,包括肌肉生长抑制素(myostatin,MSTN)、肌球蛋白轻链磷酸化快肌(myosin light chain, phosphorylatable, fast skeletal muscle,MYLPF)、快速骨骼肌钙蛋白T3(troponin T3, fast skeletal type,TNNT3)和脂肪酸结合蛋白3(FABP3)/FABP4等。Zhang等[16]对云岭牛和中国西门塔尔牛的背肌转录组检测出1 347个DEGs,且被差异调节的主要代谢途径是脂肪分解和糖代谢。Ueda等[17]通过RNA-Seq鉴定出与日本和牛皮下和IMF形成的8个关键基因,包括羧肽酶E(carboxypeptidase E,CPE)、生腱蛋白C(tenascin C,TNC)、转凝蛋白(transgelin,TAGLN)、Ⅳ型胶原α5(collagen type IV α5,COL4A5)、半胱氨酸和甘氨酸富集蛋白2(cysteine and glycine-rich protein 2,CSP2)、PDZ和LIM域3(PDZ and LIM domain 3,PDLIM3)、磷酸酶1调节抑制因子14A亚基(phosphatase 1 regulatory inhibitor subunit 14A,PPP1R14A)和钙调神经磷酸酶调节因子2(regulator of calcineurin 2,RCAN2)。本试验对宁夏地区西门塔尔和安格斯肉牛的背最长肌肉质特性进行了对比研究及RNA-Seq,结合生物信息鉴定出可能影响IMF沉积的DEGs,并对其表达特性进行分析,为阐明肉牛IMF沉积的调控机制和改善肉品质提供一定的理论借鉴。

1 材料与方法

1.1 试验动物和样品采集

选用来自宁夏永宁某养殖场亲缘关系和月龄相近、饲粮结构相同以及饲养管理条件一致的西门塔尔和安格斯公牛各3头,平均体重为(700±10) kg,组间差异不显著(P>0.05)。

试验牛于2022年5月统一屠宰,采集胴体第12根肋骨处背最长肌样品。取500 g样品于4 ℃熟化24 h,用于熟肉率和系水力(参考NY/T 1333—2007方法)、剪切力(参考NY/T 1180—2006方法)以及氨基酸(参考GB 5009.124—2016方法)和脂肪酸(参考GB 5009.168—2016方法)含量的测定;另取500 g样品于-20 ℃保存,用于蛋白质和IMF含量的测定(剔除表面筋膜和脂肪);采集5 g左右样品组织于专用固定液保存,用于肌纤维组织特性的测定[苏木精-伊红(HE)染色];取1 cm3左右样品锡纸包裹投入液氮于-80 ℃保存,用于RNA提取。

1.2 RNA的提取与建库测序

采用TRIzol(Invitrogen,美国)试剂盒从背最长肌样品中提取RNA样本,通过Agilent 2100生物分析仪(Agilent Technologies,美国)检测RNA样本的完整性和纯度。筛选出符合28S rRNA∶18S rRNA≥1.0且RNA完整度(RIN)≥7条件的样本,通过Oligo(dt)磁珠吸附法构建mRNA文库,库检合格后,不同文库按照有效浓度及目标上机数据量的需求进行Illumina NovaSeq 6000测序。

1.3 测序数据分析

对测序获取的原始数据筛选过滤后得到高质量数据(clean reads),利用HISAT2 2.0.5软件与牛(Bostaurus)的参考基因组比对(参考基因组版本为USDA ARS-UCD1.2)[18]。采用FeatureCounts 1.5.0统计每个基因的每千个碱基的转录每百万映射读取的片段(fragments per kilobase of transcript per million mapped reads,FPKM),再结合StringTie 1.3.3b对基因进行定量[19]。采用DESeq2 1.20.0软件进行2个肉牛品种之间的差异基因表达分析[20],显著差异表达阈值为校正后的P值(Padj)<0.05及|log2[差异倍数(fold change,FC)]|>1,对DEGs进行主成分分析(PCA),并可视化每个个体的基因表达谱。

1.4 DEGs富集分析

采用clusterProfile 3.8.1 R语言包对DEGs集进行基因本体(gene ontology,GO)功能和京都基因与基因组百科全书(Kyoto encyclopedia of genes and genomes,KEGG)通路富集分析[21],均以Padj<0.05作为显著性富集的阈值,根据DEGs的显著性水平和肌纤维、脂质发育关联注释信息筛选出与IMF沉积相关的主控基因和调控通路。

1.5 实时荧光定量PCR(RT-qPCR)分析

为验证RNA测序数据,采用RT-qPCR检测DEGs。背最长肌组织分离的总RNA样品采用PrimeScriptTMRT试剂盒(TaKaRa,日本)反转录成cDNA,去除基因组DNA。采用SYBR Green试剂盒(QIAGEN,美国)进行RT-qPCR分析。将所选DEGs的表达水平与内参基因β-肌动蛋白(β-actin)进行归一化处理。PCR反应条件为:95 ℃反应10 min,40个循环,95 ℃反应15 s,60 ℃反应60 s。采用2-ΔΔCt法计算DEGs的相对表达量。所有供试基因引物信息见表1。

表1 引物信息

1.6 数据统计分析

采用SPSS 19.0统计软件中的独立样品t检验进行差异显著性分析,结果用“平均值±标准差”表示,P<0.05为差异显著,P<0.01为差异极显著,P>0.05为差异不显著。

2 结果与分析

2.1 西门塔尔牛和安格斯牛肉品质及肌纤维特性分析

由表2可知,安格斯牛IMF含量和系水力极显著高于西门塔尔牛(P<0.01),肌纤维直径显著低于西门塔尔牛(P<0.05);安格斯牛肌肉粗蛋白质含量、熟肉率和肌纤维密度都高于西门塔尔牛,但差异均不显著(P>0.05)。

2.2 西门塔尔牛和安格斯牛肌肉氨基酸和脂肪酸含量分析

由表3可知,试验测定了非必需氨基酸(non essential amino acid,NEAA)和必需氨基酸(essential amino acid,EAA)各8种,2个品种牛肌肉的氨基酸组成成分一致,但含量存在差异。在NEAA中,安格斯牛肌肉酪氨酸含量极显著高于西门塔尔牛(P<0.01),肌肉精氨酸含量极显著低于西门塔尔牛(P<0.01)。在EAA中,2个品种牛肌肉各氨基酸含量差异均不显著(P>0.05)。2个品种牛肌肉NEAA和EAA的总含量差异均不显著(P>0.05),且EAA/TAA值均接近40%(理想值),即2个品种牛肌肉蛋白质含量丰富。

肌肉中的脂肪酸可分为SFA、单不饱和脂肪酸(monounsaturated fatty acid,MUFA)和多不饱和脂肪酸(polyunsaturated fatty acid,PUFA)。由表4可知,2个品种牛分别检测到4种SFA、3种MUFA和4种PUFA,且各分类中的脂肪酸含量差异均不显著(P>0.05),但西门塔尔牛肌肉SFA总含量高于安格斯牛,肌肉MUFA总含量低于安格斯牛,且肌肉PUFA含量2个品种接近,脂肪酸组成整体差异不大。

表2 西门塔尔牛和安格斯牛肉品质及肌纤维特性分析

表3 西门塔尔牛和安格斯牛肌肉氨基酸含量分析

表4 西门塔尔牛和安格斯牛肌肉脂肪酸含量分析

2.3 西门塔尔牛和安格斯牛肌肉RNA-Seq结果分析

由表5可知,采用Illumina NovaSeq 6000平台对西门塔尔牛和安格斯牛背最长肌进行RNA-Seq,结果显示,每组3个重复平均共产生45 677 571条cleanreads,150bp clean reads的碱基质量值≥30的碱基所占百分比(Q30)平均值为93.96%,而每个样品的平均鸟嘌呤和胞嘧啶(GC)含量为51.02%,表明测序数据的结果可靠,质量稳定。此外,每个样本约96%被定位到参考基因组,有93%左右被定位到唯一参考基因组的基因区域。

表5 西门塔尔牛和安格斯牛肌肉RNA-Seq结果分析

2.4 西门塔尔牛和安格斯牛肌肉DEGs分析

由图1可知,通过盒形图展示西门塔尔牛和安格斯牛之间基因表达值FPKM分布的对比情况发现,西门塔尔牛和安格斯牛肌肉DEGs的表达模式相似。PCA结果(图2)表明,主成分1(PC1)和主成分2(PC2)分别为42.88%和20.57%。由图3可知,安格斯牛3个重复样本聚合紧密,重复性较好;而西门塔尔牛组内分散,重复性较差,但2个品种之间比较分散,说明所用样本分组合理。由图4可知,与安格斯牛相比,西门塔尔牛共有1 214个候选基因在肌肉组织中差异表达,其中上调基因736个,下调基因478个。

图1 西门塔尔牛和安格斯牛肌肉基因

图2 西门塔尔牛和安格斯牛肌肉FPKM PCA结果

2.5 西门塔尔牛和安格斯牛肌肉DEGs GO功能富集分析

DEGs注释到GO数据库中共显著富集在了43个GO条目,在生物过程(biological process,BP)、细胞组分(cellular component,CC)和分子功能(molecular function,MF)分类中,包含的亚类分别为20、15和8个。从每个分类中各选取最显著的前10个(分子功能为8个)条目绘制柱状图(图5)。在生物过程分类中,参与生物合成和代谢的DEGs最多,数量分别为103和121个;在细胞组分分类中,细胞质部分和细胞内无膜细胞器富集的DEGs最多,数量分别为123和117个;在分子功能分类中,DEGs主要参与核糖体结构成分、结构分子活性和磷酸酯水解酶活性,数量分别为71、82和50个。

2.6 西门塔尔牛和安格斯牛肌肉DEGs KEGG通路富集分析

从KEGG通路富集分析结果中,选取最显著的20个注释绘制散点图(图6),结合GeneCards在线查询基因功能,进一步筛选与IMF沉积相关的DEGs,发现其主要富集在了心肌收缩(cardiac muscle contraction)、非酒精性脂肪肝病(non-alcoholic fatty liver disease)和逆行内源性大麻素信号(retrograde endocannabinoid signaling)通路中,数量分别为42和47个;同时,共筛选出3个可能与IMF沉积相关的DEGs。与安格斯牛相比,西门塔尔牛肌肉淀粉样前体蛋白(APP)为下调基因,肌肉原肌球蛋白2(TPM2)和钙电压门控通道亚基α1 S(CACNA1S)为上调基因,这3个调控基因的关键信息见表6。

2.7 RT-qPCR验证结果

由图7可知,在2种肉牛品种背最长肌中,安格斯牛肌肉CACNA1SRNA-Seq表达量极显著高于西门塔尔牛(P<0.01),肌肉TPM2 RNA-Seq表达量也显著高于西门塔尔牛(P<0.05),表明这2个基因对IMF沉积有促进作用;而肌肉APPRNA-Seq表达量在2种牛之间差异不显著(P>0.05)。RT-qPCR验证结果与RNA-Seq结果相一致。

3 讨 论

牛肉IMF含量与肉质的细嫩多汁和口感风味等特性密切相关[22]。通常来说,IMF与肌肉嫩度[23-25]、系水力[26]、多汁性[27]以及肌纤维密度[28]呈正相关,与剪切力[29]和肌纤维直径[28]呈负相关。本试验测定肉质相关性状发现,安格斯牛的IMF含量和系水力极显著高于西门塔尔牛,肌纤维直径显著低于西门塔尔牛,肌肉粗蛋白质含量、熟肉率和肌纤维密度都高于西门塔尔牛,而剪切力低于西门塔尔牛。这与上述IMF与肉品质的相关性规律一致,也与赵伟明等[30]和韩信嵘[2]研究发现安格斯牛IMF含量高、肉品质好的结果一致。此外,2个品种肉牛肌肉氨基酸组成的分析结果发现,两者EAA/TAA值都接近40%理想值,表明肌肉中都具有丰富的优质蛋白质。不过,与西门塔尔牛相比,安格斯牛肌肉的SFA总含量较低,而肌肉的MUFA总含量较高,但差异均不显著,且两者的PUFA总含量接近,总体而言,2个品种肌肉脂肪酸组成差异不大。通常来说,过高的肌肉SFA含量不利于人体健康,将会提升血液胆固醇水平[31],而不饱和脂肪酸不仅可以阻止动脉硬化[32],还能预防高血脂和冠心病的发生[33]。Chambaz等[4]在背最长肌平均IMF含量为3.25%的条件下,比较安格斯、西门塔尔、夏罗莱和利木赞阉牛的肉质和大理石花纹特性发现,安格斯牛的脂肪含量得分最高,IMF沉积最好,同时安格斯和利木赞牛比西门塔尔牛肉质更嫩,这进一步说明了安格斯牛肉的优秀品质。为了进一步确定影响这些表型差异的候选基因和通路,本试验采用Illumina NovaSeq 6000平台对安格斯和西门塔尔2个肉牛品种的背最长肌组织进行RNA-Seq,并通过对DEGs GO功能和KEGG通路进行富集分析,鉴定出与IMF沉积有关的功能基因,结果发现,一些DEGs参与了脂肪形成的生物过程,一些DEGs参与了脂肪组织发育的重要信号通路,结合GeneCards和RT-qPCR对相关功能基因进一步验证筛选,得出了APP、TPM2和CACNA1S这3个基因,可为肉牛的脂肪沉积机制研究提供参考依据。

图3 西门塔尔牛和安格斯牛肌肉DEGs热图及聚类分析

图4 西门塔尔牛和安格斯牛肌肉DEGs火山图

APP是一类具有胞外N端区、跨膜结构域和胞质C端尾短的大蛋白分子[34]。APP与阿尔茨海默病(Alzheimer disease,AD)斑块相关,已知存在于大脑突触和成人神经肌肉接头中[35],还存在于小鼠骨骼肌和培养的肌源细胞中,成年小鼠骨骼肌中APPmRNA下调,会抑制脂肪细胞的成脂分化[36-37]。此外,Wu等[38]研究发现,当APP在小鼠肌肉中过度表达时,会引起肌无力和萎缩,表明APP低表达有利于肌肉及脂肪的适度沉积;通过荧光素酶报告基因检测证明,APP是miR-101-1的潜在靶基因,而miR-101-1是一种广泛表达的微小RNA(miRNA),主要在牛骨骼肌中富集,其丰度在胎牛、犊牛和成年牛的骨骼肌中存在差异;在C2C12成肌细胞细胞系分化过程中,内源性miR-101-1会上调,转染miR-101-1会抑制C2C12成肌细胞分化,同时显著降低肌原分化(myogenic differentiation,MyOD)、肌细胞生成素(myogenin,MyOG)和肌球蛋白重链(myosin heavy chain,MyHC)的mRNA表达。APP与肉质性状有强相关,可作为猪肉剪切力和嫩度的数量性状位点(QTL)[39]。而Lee等[40]在人体脂肪组织中研究发现,APP参与肥胖相关胰岛素抵抗和脂肪组织炎症,在脂肪组织和肥胖群体中均高表达,其表达水平与脂肪沉积水平呈正相关,这与本试验结果呈相反趋势,可能与研究对象不同有关,其脂肪组织的细胞组成及成脂机制存在较大差异,具体原因需进一步研究。

图5 西门塔尔牛和安格斯牛肌肉DEGs GO功能富集分析

CACNA1S是编码骨骼肌细胞中缓慢失活的L型电压依赖性钙通道的5个亚基之一,基因突变与骨骼肌代谢生长和脂肪堆积相关[41]。骨骼肌电压传感器Cav1.1也是一个电压依赖性钙离子(Ca2+)通道,允许少量Ca2+内流进入哺乳动物骨骼肌纤维[42]。Georgiou等[43]研究发现,骨骼肌中发现了一条通路Cav1.1→钙离子/钙调蛋白依赖激酶Ⅱ(CaMKⅡ)→一氧化氮和酶(NOS),能调节脂肪酸转运蛋白CD36在细胞内的分布从而改变脂肪酸代谢。阻断这一通路途径可减少线粒体β氧化和能量消耗,促进肌纤维密度增大,由此说明CACNA1S能通过影响肌肉线粒体功能调节骨骼肌的能量消耗降低,从而导致IMF的沉积。Xue等[44]研究发现,CACNA1S还涉及胰岛素分泌、抵抗和信号通路相关的表观遗传调控网络,通过调控长链非编码RNA(long noncoding RNA,lncRNA)的变化和甲基化影响能量代谢、信号通路和细胞增殖相关过程,是脂质代谢的上游调控途径,还与其他靶基因共同调控肌内沉积过程。Ma等[45]研究结果也证实了上述结果,CACNA1S在启动子区发生甲基化,与修饰基因表达呈负相关。Xue等[44]和Ma等[45]分别在蒙古羊和牦牛中均发现CACNA1S在启动子区发生甲基化,参与脂质代谢的上游调控途径,调控动物肌肉的发育和肉质。方晓敏等[46]推测CACNA1S是一个适合作为影响猪肉品质性状QTL的候选基因。有研究发现,基本螺旋-环-螺旋芳香烃受体核转位蛋白样1(BMAL1)通过反向定位c-erbA-1α(reverse orientation the c-erbA-1alpha,REV-ERBα)控制CACNA1S表达来调节钙信号,从而决定肌肉合成代谢或分解状态的信号通路,可能是肌肉活动的分子传感器[47]。

图6 西门塔尔牛和安格斯牛肌肉DEGs前20 KEGG通路富集气泡图

表6 IMF沉积相关DEGs

TPM2编码β-原肌球蛋白,是肌动蛋白纤维结合蛋白家族的一员,主要表达在慢速Ⅰ型肌纤维中,该基因的突变可改变其他肌原肌球蛋白蛋白的表达。TPM2在肌肉和非肌肉细胞中与肌动蛋白纤维结合,与肌钙蛋白复合物一起,在脊椎动物横纹肌收缩的钙依赖调节中起核心作用[48]。Cho等[49]通过RNA和蛋白定量验证了原肌球蛋白家族成员在韩牛、阉牛和公牛的肌内脂肪组织(intramuscular adipose tissue,IMAT)和大血管脂肪组织(macrovascular adipose tissue,MAT)中参与肌肉发育密切相关的调控通路,与奶牛或阉牛相比,IMAT中TPM2 mRNA和蛋白表达水平均较低,表明它们与大理石花纹评分和质量等级呈正相关。本试验中,DEGs富集注释结果显示,TPM2与肌原纤维、收缩纤维部分和肌动蛋白细胞骨架结构的KEGG通路相关,这与之前的研究结果[50-51]也一致,这为TPM2直接或间接参与脂肪沉积和脂肪酸代谢提供了证据。

*表示差异显著(P<0.05),**表示差异极显著(P<0.01),ns表示差异不显著(P>0.05)。

在本研究中,在RT-qPCR验证DEGs表达结果中,安格斯牛肌肉CACNA1S和TPM2 RT-qPCR相对表达量显著高于西门塔尔牛,表明其对IMF沉积有促进作用,且与RNA-Seq的结果趋势一致;而2种肉牛之间肌肉APP表达量差异不显著,但相比于西门塔尔牛,安格斯牛肌肉APP表达量呈降低趋势,表明其对脂肪沉积有负调控作用,这与RNA-Seq结果也是一致的,这也进一步支持了APP、CACNA1S和TPM2能够参与调节IMF沉积关键通路,具有调节脂肪细胞分化和脂肪酸新陈代谢的分子功能的结论。

4 结 论

安格斯牛IMF沉积能力较强,肉品质较好。RNA-Seq共筛选出了3个(APP、CACNA1S和TPM2)与肉牛IMF沉积相关的候选基因,这为研究肉牛IMF沉积调控机制提供了参考。