基于转录组学筛选绵羊肉质性状相关候选基因

梁 鹏，张 稳，冯登侦，强 浩，荣 轩，孟 科

(宁夏大学 农学院，宁夏 银川 750021)

我国是绵羊养殖和羊肉消费大国，绵羊种类占世界绵羊品种的10%左右[1]。经过绵羊育种工作者多年的努力和研究，我国的绵羊良种繁育体系也在不断完善。同时，随着经济的快速发展和人民生活水平的提高，人们对肉类消费的理念已经从单纯追求数量转变为数量和质量并重。因此，人们也把绵羊育种的目标逐渐转向为培育具有多优良性状的肉羊新品种(系)，以满足市场需求的供应。研究表明，品种是影响畜禽肉品质的主要因素[2]，也是种质基因库的重要保存单位，一些具有特色的畜禽品种便成为理想的研究对象。

滩羊是在宁夏地区经过漫长的自然和人工选择形成的优良特色畜种，其肉质以鲜美细嫩、营养丰富和风味独特而闻名，是世界认可的优质羊肉之一[3]，但滩羊个体小、早期生长发育慢、产肉力不足以及繁殖力低下，羊肉总产量只能靠多饲养、高出栏为主，这使得滩羊养殖的成本增加，且养殖效益低，也不符合当前的产业发展要求。因此，为了达到对滩羊这一种质资源的创新与利用，选用产肉性能好的杜泊羊和繁殖性能高的小尾寒羊与本地滩羊进行三元杂交，以期利用各亲本的优良性状培育出适合宁夏地区资源环境特点的肉羊新品种(系)。

利用传统方法对肉品质和风味的研究只能通过检测肉嫩度、脂肪酸和氨基酸含量等进行分析，而转录组测序技术能够深入解析与这些表型性状有关的调控因子，从分子水平上对调控肉品质和风味性状的相关基因进行筛选和鉴定。目前，转录组测序技术已经广泛应用于畜禽功能性状基因的挖掘与育种研究中，Wang等[4]利用转录组学技术研究了影响晋南牛与西门塔尔牛肌内脂肪(IMF)含量差异的相关调控基因，结果筛选到的这些差异基因(CYP21A2、PC、ACACB、APOA1和FADS2)可能通过调节脂质代谢来影响IMF的沉积。Sun等[5]利用转录组学技术研究了千华肉用美利奴羊和小尾寒羊之间与肌肉生长发育有关的转录调节因子，结果显示，960个基因在品种间差异表达，通过代谢途径鉴定到MRF、GXP1和STAC3等差异基因在肌肉生长发育过程中起着重要作用。韩信嵘等[6]以背膘组织为研究对象，利用转录组学技术分析安格斯牛和西门塔尔牛与脂肪沉积有关的差异表达基因，结果筛选到DKKL1、ACACB和PTGS2基因通过cAMP通路调控脂肪沉积。

本研究以宁夏优质肉羊培育中的3个亲本为研究对象，旨在通过转录组学技术挖掘影响不同品种绵羊肉品质和风味性状差异的相关调控基因，为宁夏优质肉羊新品种(系)培育提供理论依据。

1 材料和方法

1.1 试验样品采集

选择相同饲养环境下的8月龄杜泊羊、滩羊和小尾寒羊各8只(公母各1/2)作为研究对象，试验动物均饲养于宁夏宇泊科技有限公司威震肉羊繁育场，屠宰前禁食24 h、停水2 h，于颈静脉放血屠宰后，迅速采集第3～12肋骨之间的背最长肌组织约500 g，剔除表面筋膜后装入自封袋中并标记，用于肉品质指标测定。另取少部分样品装入无RNA酶的冻存管中，标号后快速置于液氮罐中，待试验结束后将样品带回实验室，于-80 ℃冰箱保存备用。

1.2 肉品质指标测定

采用石蜡切片法观察4%多聚甲醛固定样品的肌纤维直径和密度；使用pH计测定屠宰24 h后的pH值；使用嫩度仪和压力法分别测定剪切力和系水力；按照GB/T9695.7—2008《肉及肉制品中脂肪含量测定》测定肌肉中的脂肪含量；采用高效液相色谱法测定肌肉中肌苷酸和肌苷含量。

1.3 文库构建及转录组测序

随机选取杜泊羊、滩羊和小尾寒羊母羊各3只，提取背最长肌组织总RNA，利用琼脂糖凝胶电泳、Nanodrop和Agilent 2100 bioanalyzer检测RNA样品的浓度、纯度及完整性。经质检合格后采用去除核糖体RNA的方法构建链特异性文库，文库构建完成后，先使用Qubit 2.0进行初步定量，稀释文库至1.5 ng/μL，随后使用Agilent 2100对文库的插入片段大小(Insert size)进行检测，符合预期后，使用qRT-PCR方法对文库有效浓度(高于3 nmol/L)进行准确定量，以保证文库质量。文库经检验合格后，利用Illumina HiSeqTM4000平台进行PE150双端测序。

1.4 转录本拼接、筛选及预测

对测序获得的原始数据进行过滤、错误率及GC含量分布检查，获得后续分析使用的高质量数据(Clean reads)。利用Hisat 2(http：//ccb.jhu.edu/software/hisat2)软件将Clean reads数据比对到绵羊参考基因组上，然后使用Stringtie软件在基于比对到基因组上的结果，将reads拼接成转录本进行定量。使用Cuffmerge软件对各样本拼接得到的转录本进行合并，去掉其中链方向不确定和转录本长度不超过200 nt的转录本，之后利用Cuffcompare软件和已知数据库进行比较，过滤掉数据库中已知的转录本，最后对筛选到的新转录本进行编码潜能预测，预测结果取CPC、PFAM和CNCI的交集。

1.5 基因可变剪切分析

使用rMATS(http：//rnaseq-mats.sourceforge.net/index.html)软件对组织中基因的可变剪切结果进行统计分析。rMATS分析的可变剪切事件有外显子跳跃(Skipped exon，SE)、外显子互斥(Mutually exclusive exon，MXE)、5′端可变剪切(Alternative 5′ splice site，A5SS)、3′端可变剪切(Alternative 3′ splice site，A3SS)和内含子滞留(Retained intron，RI)这5种。

1.6 组间差异基因(DEGs)筛选

通过StringTie分析流程计算每个转录本在各样本中的表达水平，结果以FPKM(Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced)值表示。之后利用DESeq 2软件进行基因表达的差异显著性分析，以padj(矫正后的Pvalue)<0.05和|log2FC|≥1作为组间差异显著基因的筛选标准。

1.7 DEGs的GO/KEGG富集分析

为进一步了解组间DEGs的功能和参与的代谢通路，将筛选出的DEGs进行GO(http：//www.geneontology.org/)和KEGG(Kyoto Encyclopedia of Genes and Genomes)富集分析。

1.8 qRT-PCR验证

随机选取13个差异表达基因(7个上调，6个下调)进行实时荧光定量PCR(Reverse transcription-quantitative PCR，qRT-PCR)，以验证转录组测序数据的可靠性。利用FastKing RT Kit(With gDNase)试剂盒(天根)将原始样品的RNA反转录得到cDNA。内参基因选用甘油醛-3-磷酸脱氢酶基因(Glyceraldehyde-3-phosphate dehydrogenase gene，GAPDH)，引物均由中科羽瞳(陕西)生物科技有限公司合成，序列信息见表1。荧光定量PCR的反应体系为20 μL：TB Green Premix Ex Taq Ⅱ 10 μL，ROX Reference Dye Ⅱ 0.4 μL，上、下游引物(10 μmol/L)各0.8 μL，cDNA 2 μL，ddH2O 6 μL。反应条件为：95 ℃ 30 s； 95 ℃ 10 s，退火30 s，40个循环；65～95 ℃，每5 s增加0.5 ℃制作熔解曲线。每个样品重复3次，利用2-ΔΔCt法计算各基因相对表达量，并通过Log2FoldChange标准化处理。

1.9 数据处理

利用SPSS 22.0软件对肉品质指标进行单因素方差分析，用Duncan法进行多重比较，数据结果以“平均值±标准误”表示。以P<0.05表示组间差异显著。

表1 引物信息Tab.1 Primers information

2 结果与分析

2.1 肉品质差异分析

通过对绵羊的肉品质指标进行测定，结果显示(表2)：杜泊羊肌纤维密度显著低于滩羊和小尾寒羊(P<0.05)，肌纤维直径显著高于滩羊(P<0.05)。小尾寒羊的剪切力和失水率均显著高于杜泊羊和滩羊(P<0.05)，而脂肪含量显著低于滩羊(P<0.05)，与杜泊羊差异不显著(P>0.05)。滩羊的pH24h值显著低于小尾寒羊和杜泊羊(P<0.05)。杜泊羊的肌苷酸含量显著低于滩羊和小尾寒羊(P<0.05)，而肌苷含量显著高于滩羊(P<0.05)，与小尾寒羊差异不显著(P>0.05)。

表2 绵羊肉品质差异分析Tab.2 Analysis of sheep meat quality difference

2.2 测序数据及质量评估

通过建立9个背最长肌组织样品的文库，上机测序后数据经过过滤、错误率和GC含量分布检查，共获得808 043 288条有效序列(表3)，有效碱基达到121.21 G。各样品的GC含量在47.12%～52.04%，Q20≥97.63%，Q30≥93.49%，错误率≤0.03%；测序的有效数据通过与参考基因组比对，结果显示比对率均在89.25%以上。这些结果都表明，测序数据质量可靠，可用于后续的分析。

表3 测序数据质量评估Tab.3 Sequencing data quality assessment

2.3 背最长肌组织转录本筛选与预测

基于比对结果，将reads拼接成转录本，共得到952 263个转录本。通过将各样本拼接的转录本进行合并，且除去不符合要求的转录本后，将剩余转录本与已知数据库进行比较，过滤掉数据库中已知的转录本，得到新转录本8 083个，其中位于“+”链的新转录本3 992个，位于“-”链的新转录本4 091个。新转录本的长度为201～71 939 bp，外显子个数为2～35个，开放阅读框(Open reading frame，ORF)在29～1 807 bp。经编码潜能预测，结果显示，新转录本中280个具有编码潜能，可作为Novel-mRNA数据集。

2.4 基因可变剪切分析

通过rMATS软件对各组织样品中的可变剪切进行定量，结果见表4。由表4可知，各比较组中检测到的可变剪切事件分别为18 239，18 258，19 222个，说明许多基因都存在多种剪切方式，造成基因在表达水平上发生变化，进而使生物体能够表达多种蛋白，引起物种的表型性状产生差异。在5种可变剪切事件中，各比较组合均以SE事件最多，MXE事件次之，A5SS事件最少。

2.5 DEGs的筛选和聚类分析

通过FPKM方法计算各样本中基因的表达水平，结果以盒形图展示，由图1可知，不同品种绵羊背最长肌组织样品中基因的表达水平有差异。通过计算各样品间皮尔逊相关系数(Pearson′s correlation coefficient)，可以说明重复样品的相关性强弱，r越接近1，表明样品相关性越强。本研究中各样品的r值均大于0.95，说明试验的可靠性和组内样品的重复性较高。使用DESeq 2软件对定量分析完成的样品进行组间基因表达差异显著性分析，结果显示，共鉴定出820个DEGs。杜泊羊与滩羊组鉴定出99个DEGs，其中51个上调，48个下调；小尾寒羊与杜泊羊组鉴定出436个DEGs，其中270个上调，166个下调；小尾寒羊与滩羊组鉴定出552个DEGs，其中320个上调，232个下调。3个组共表达差异显著基因有4个(MAOB、SMAD1、AMPD3和ENAH)。对筛选出来的DEGS进行聚类分析(图2)，结果显示，同一品种的3个生物学重复样本能够很好地聚到一起，说明本研究所用样本的生物学重复较好。小尾寒羊与杜泊羊和滩羊的基因表达模式出现明显分离，说明组间差异较大，而杜泊羊和滩羊重叠区域较多，说明组间差异较小。

每个区域的盒形图对应5个统计量(从上而下分别为最大值、 上四分位数、中值、下四分位数和最小值)。 The box plot for each region corresponds to five statistics(maximum， upper quartile，median，lower quartile and minimum from top to bottom).

2.6 差异表达基因的GO/KEGG富集分析

通过对DEGs进行GO富集分析，可充分了解到DEGs在绵羊之间参与的生物学过程及生理功能。GO功能可分为三大类：生物学过程(Biological process，BP)、分子功能(Molecular function，MF)和细胞组分(Cellular component，CC)，本研究利用GOseq软件对各组间DEGs进行GO富集分析，结果见图3所示：杜泊羊与滩羊组DEGs显著富集到224个GO条目中，小尾寒羊与杜泊羊组DEGs显著富集到517个GO条目中，小尾寒羊与滩羊组DEGs显著富集到657个GO条目中。在生物学过程中，DEGs主要富集在代谢过程、有机物生物合成过程、 细胞对化学刺激的反应、横纹肌细胞分化、横纹肌组织发育、肌细胞分化、肌肉组织发育和肌肉结构发育等过程；在细胞组分中，主要富集在膜固有的、高分子络合物、细胞内非膜结合细胞器、细胞外区和核糖体等；在分子功能中，主要富集在RNA聚合酶Ⅱ调节区序列特异性DNA结合、结构分子活性、核糖体的结构成分等。

左侧根据表达相似程度对基因进行聚类，上方根据表达谱的相似程度对每个样本进行聚类，由蓝至红表达量逐渐上调，数字为均一化后的相对表达量。

利用KOBAS软件对组间DEGs进行KEGG代谢途径和信号通路分析，结果显示，DEGs共富集到227条通路，其中杜泊羊与滩羊组DEGs富集到121条代谢通路中，小尾寒羊与杜泊羊组DEGs富集到197条代谢通路中，小尾寒羊与滩羊组DEGs富集到206条代谢通路中，3组DEGs共富集通路有104条。按P值从小到大选择Top 20代谢通路，以散点图展示(图4)：只有1条代谢通路显著富集(P<0.05)，为核糖体通路。剩下的19条以及其他代谢通路富集虽不显著，但均有不同程度的DEGs富集，其中与肉品质和风味有关的通路及DEGs如表5所示。通过GO功能注释、KEGG通路分析以及相关文献报道，筛选到肌苷酸沉积(AMPD1、AMPD3、NT5C1A)、脂肪沉积(LPIN1、FABP4、ACOT7、DUSP1、GOT1、CKMT2、CAPN2、SLC2A4、SCOS3、RXRG、PPARGC1A)、肌纤维调控(FOS、ACTC1、GYS1、CS)、肌肉糖酵解/糖异生(PFKM、PKM、GPI、PGAM1、FBP2)、肌肉发育(CKM、SMAD1)和脂质代谢(PIK3C2G、ABCA1)等与肉品质和风味有关的基因。

A：1.细胞对镉离子的反应；2.细胞对铵离子的反应；3.含苯化合物代谢过程；4.骨吸收的正调控；5.骨重塑的正调控；6.中胚层细胞命运承诺；7.对镉离子的反应；8.硝基苯代谢过程；9.氮化合物解毒；10.含氮化合物的细胞解毒；11.毒素生物合成过程；12.生物碱代谢过程；13.二苯并对二噁英代谢过程；14.氧化脱乙基；15.组织重塑的正调控；16.单萜类代谢过程；17.翻译的昼夜调节；18.参与初级生殖层形成的细胞命运承诺；19.中胚层细胞分化；20.异源分解代谢过程；21.AP1复合体；22.核糖体小亚基；23.胞体纤维；24.血红蛋白复合物；25.生长细胞尖端；26.新生长细胞尖端；27.细胞尖端；28.精子纤维鞘；29.纤毛运动；30.嗜铬颗粒膜；31.曼切特；32.嗜铬颗粒；33.膜固有的；34.轴突丘；35.纤毛；36.核糖体亚单位；37.胞质小核糖体亚单位；38.膜神经元体细胞；39.氧化还原酶活性，作用于CH或CH2基团、醌或类似化合物作为受体；40.咖啡因氧化酶活性；41.芳香化酶活性；42.转铁蛋白受体活性；43.吡喃结构域结合；44.货物受体活性；45.犬尿氨酸3-单加氧酶活性；46. UDP-N-乙酰氨基葡萄糖-溶酶体-酶N-乙酰氨基葡萄糖磷酸转移酶活性；47.[硫酸乙酰肝素]-氨基葡萄糖3-磺基转移酶1活性；48. P物质受体活性；49.钙和钙调素反应性腺苷酸环化酶活性；50.RNA聚合酶Ⅱ调节区序列特异性DNA结合；51.RNA聚合酶Ⅱ调节区DNA结合；52.氧转运体活性；53.细胞骨架的结构成分；54.AMP脱氨酶活性；55.补体成分C1q结合；56.视蛋白结合；57.硫转移酶活性；58.速激肽受体活性。B：1.翻译；2.有机物生物合成过程；3.大分子生物合成过程；4.生物合成过程；5.细胞生物合成过程；6.细胞大分子生物合成过程；7.基因表达；8.细胞蛋白质代谢过程；9.蛋白质代谢过程；10.大分子代谢过程；11.代谢过程；12.蛋白质折叠；13.核糖体；14.核糖核蛋白复合物；15.胞浆核糖体；16.核糖体亚单位；17.胞浆大核糖体亚单位；18.肌原纤维；19.收缩纤维；20.大核糖体亚单位；21.非膜结合细胞器；22.细胞内非膜结合细胞器；23.肌节；24.胞浆部分；25.高分子络合物；26.收缩纤维部分；27.胞质小核糖体亚单位；28.多体核糖体；29.核糖体小亚基；30.条纹；31.核糖体的结构成分；32.结构分子活性；33.未折叠蛋白结合；34.rRNA结合。C：1.细胞趋化性；2.细胞对化学刺激的反应；3.蛋白质折叠；4.横纹肌细胞分化；5.肌细胞分化；6.呼吸爆发参与防御反应；7. 翻译；8.肌肉组织发育；9.肌肉结构发育；10.横纹肌组织发育；11.肌节；12.肌原纤维；13. 收缩纤维部分；14.收缩纤维；15.条纹；16.核糖体；17.胞浆核糖体；18.细胞外区；19.Z盘；20.胞浆部分；21.核糖体小亚基；22.未折叠蛋白结合；23.核糖体的结构成分；24.结构分子活性；25.rRNA结合；26.生长因子活性。图A代表D vs T组，图B代表XHX vs D组，图C代表XHX vs T组。图4，5同。

A：1.胆汁分泌；2.色氨酸代谢；3.昼夜夹带；4.卵巢类固醇生成；5.药物代谢-细胞色素P450；6.糖胺聚糖生物合成-硫酸乙酰肝素/肝素；7.安非他明成瘾；8.昼夜节律；9.cAMP信号通路；10.胰岛素分泌；11.亚油酸代谢；12.雌激素信号通路；13.可卡因成瘾；14.胆碱能突触；15.ABC转运；16.血管平滑肌收缩；17.长期抑郁症；18.细胞色素P450对外源性物质的代谢；19.多巴胺能突触；20.类固醇激素生物合成。B：1.核糖体；2.雌激素信号通路；3.昼夜夹带；4.胰岛素信号通路；5.2-氧羰基酸代谢；6.膀胱癌；7.抗原处理和呈递；8.缝隙连接；9.安非他明成瘾；10.谷氨酸能突触；11.淀粉和蔗糖代谢；12.军团菌病；13.慢性髓系白血病；14.酪氨酸代谢；15.p53信号通路；16.卵巢类固醇生成；17.酗酒；18.胰岛素分泌；19.ErbB信号通路；20.精氨酸和脯氨酸代谢。C：1.核糖体；2.粘着；3.军团菌病；4.雌激素信号通路；5.膀胱癌；6.胆汁分泌；7.肌动蛋白细胞骨架的调节；8.胃酸分泌；9.胰岛素信号通路；10.谷氨酸能突触；11.内质网中的蛋白质加工；12.精氨酸和脯氨酸代谢；13.安非他明成瘾；14.钙信号通路；15.前列腺癌；16.ErbB信号通路；17.甲状腺激素合成；18.疟疾；19.嘌呤代谢；20.血小板活化。

表5 部分KEGG通路与DEGs富集结果Tab.5 Results of partial enrichment of KEGG pathway with DEGs

2.7 qRT-PCR验证结果

通过qRT-PCR对随机挑选的差异基因进行验证，结果表明，各基因的定量结果与转录组测序结果的表达趋势一致(图5)，进一步证实了转录组测序结果的可靠性。

图5 DEGs的qRT-PCRFig.5 qRT-PCR results for DEGs

3 结论与讨论

肉品质和风味是消费者评价畜禽肉质好坏的重要标准，也是畜禽育种改良的重要经济指标之一。肉的风味主要与脂肪、氨基酸、核苷酸和有机酸等前体物质密切相关，在肌肉中它们含量的不同会导致肉品质和风味产生差异[7-8]。通过对杜泊羊、滩羊和小尾寒羊的肉品质指标进行测定，结果显示，品种间存在不同程度的显著差异。因此，为了解造成这些肉品质指标产生差异的潜在遗传机制，进一步采用转录组学技术对3个亲本群体进行高通量测序，筛选出组间差异表达基因及富集的通路进一步比较物质沉积能力的不同，有助于后期的绵羊分子辅助育种，加速宁夏优质肉羊新品种(系)的育成。

本研究通过对3个亲本绵羊群体的背最长肌组织进行转录组测序，以|log2FC|≥1和P<0.05为标准筛选组间DEGs，结果显示，共有820个DEGs可用于后续的生物信息学分析，其中杜泊羊与滩羊组99个，小尾寒羊与杜泊羊组436个，小尾寒羊与滩羊组552个，说明杜泊羊与滩羊差异较小，与小尾寒羊差异较大，这与检测肉品质指标分析的结果相似，进一步说明测序结果的准确性。

通过对DEGs进行GO功能注释和KEGG通路富集分析，能够深入了解基因的生物学功能和参与的代谢通路。本研究筛选出的DEGs富集到与肉品质调控和风味物质代谢有关的通路，如cAMP、TGF-β、MAPK和AMPK等经典信号通路以及嘌呤代谢、甘油磷脂代谢、氨基酸和脂肪酸代谢、肌动蛋白细胞骨架的调节和糖酵解/糖异生等通路。研究表明，cAMP信号通路通过与其下游的反应元件(CREB)和调控转录的辅助激活因子(CRTC)发挥作用，当肌肉中过表达CRTC2时会增加肌内脂肪含量和肌纤维的横截面积[9]；当过表达CRTC3时会促进甘油酰基转移酶的活性，进而增加肌肉中甘油三酯的含量[10]；在敲除小鼠的CRTC3基因后，其分解脂肪的能力增强[11]。在3T3-L1细胞分化过程中，cAMP通路通过活化CREB和C/EBPβ来促进脂肪的生成[12]。本研究结果表明，cAMP信号通路差异基因较多，其中FOS在调节细胞增殖、分化和凋亡过程中具有重要作用，能够作为启动基因转录的枢纽，通过细胞的刺激和应激调节基因表达[13]。成志敏等[14]研究表明，EGR1基因能够通过促进FOS基因的表达，引起大白猪肌细胞增长，导致肌纤维直径增大，肌纤维密度降低，本研究结果也发现，EGR1和FOS基因在小尾寒羊中均高表达，且EGR1基因显著高于杜泊羊和滩羊，而FOS基因在杜泊羊中的表达显著低于小尾寒羊和滩羊，表型测定结果显示，杜泊羊肌纤维密度显著低于滩羊和小尾寒羊，说明EGR1和FOS基因在调控杜泊羊肌纤维特性方面具有重要作用。同时，FOS基因也富集在MAPK信号通路，而MAPK通路介导细胞对外界做出反应的重要信号系统之一，不仅对脂肪细胞的增殖和分化具有重要作用[15]，在肌纤维发育过程也具有重要作用，通过调控肌纤维类型来影响肉质性状[16-17]。魏立民等[18]在研究五指山猪和长白猪背最长肌组织差异基因时，发现FOS基因通过MAPK信号通路调控肌肉发育。

AMPK通路是机体能量代谢的重要信号通路，能够为肌肉发育提供能量。过氧化物酶体增殖物激活受体γ辅激活因子1A(PPARGC1A)可以与不同的转录因子结合，参与多个代谢过程，如肌纤维转化和脂肪酸代谢等[19]。研究表明，PPARGC1A基因在猪肌肉中的表达与肌纤维面积和脂肪细胞体积呈负相关，而与脂肪组织中的表达呈正相关[20]。刘玉芳等[21]研究表明，PPARGC1A的表达量在脂肪沉积较高的金华猪背脂中显著低于瘦肉型大白猪，进一步推测该基因表达与脂肪沉积呈负相关。本研究中该基因在脂肪沉积较高的滩羊中表达量高于杜泊羊和小尾寒羊，表明该基因与脂肪沉积呈正相关，这与上述研究结果相反，可能是与研究的物种不同有关。而刘远等[22]通过转录组学技术研究表明，PPARGC1A基因在脂肪沉积高的福清山羊中表达量高于低脂肪的努比亚黑山羊，说明该基因的表达具有物种特异性。

甘油磷脂代谢通路中产生的脂肪酸和溶血磷脂等可进一步合成磷脂，而磷脂是肉风味形成的重要前体物质。该通路中的GPAM可催化动物脂质代谢中甘油酰基的合成过程，进而影响三酰甘油和磷脂的合成[23]。氨基酸代谢、脂肪酸代谢和肌醇磷酸代谢通路均与动物脂质代谢有关，这几个通路中的基因如PLA2G4E[24]、LPIN1[25]、CS[26]、SOCS3[27]等均报道与脂代谢有关，可能在调控绵羊肉品质和风味方面发挥作用，在今后还需进一步研究。

肌苷酸(IMP)又被称作次黄嘌呤核苷酸，是反映肉质鲜味的重要指标，它是机体在嘌呤代谢过程中合成的物质。研究表明，腺苷单磷酸脱氨酶(AMPD)是影响嘌呤核苷酸代谢的限速酶，与肌肉中的IMP代谢有关，它包括AMPD1、AMPD2和AMPD3共 3个亚型[28]。IMP在体内合成包括从头合成和补救合成2条途径，从头合成过程有10种酶参与反应，而AMPD1是催化AMP水解脱氨生成IMP的酶，主要在骨骼肌中表达[29]。多项研究也表明，AMPD1基因表达量与肌肉中IMP含量密切相关[30]，主要在禽类动物上报道较多。张会丰等[31]研究了城口山地鸡、大宁河鸡和青脚麻鸡胸肌中IMP含量与AMPD1基因表达量的关系，结果显示二者呈正相关；而陈星等[32]通过转录组学技术表明，AMPD1在肌苷酸含量低的沔阳麻鸭中表达量显著高于连城白鸭，二者差异倍数达到3.73倍，qRT-PCR进一步验证的结果也显示差异倍数达到2.11倍。本研究结果显示，AMPD1在肌苷酸含量最高的小尾寒羊中表达量最低，显著低于肌苷酸含量最低的杜泊羊，这与罗玉龙等[33]在苏尼特羊中的研究结果相反，推测AMPD1有可能与IMP沉积呈负相关，这与所研究的物种或品种不同有关，也可能与其他因素有关，具体还有待进一步探究。AMPD3和NT5C1A也富集在嘌呤代谢通路，有研究报道，这2个基因也与IMP代谢途径有关，因此可用于绵羊IMP沉积的候选基因进行研究[34-35]。

糖酵解/异生代谢途径是影响畜禽肉品质的重要调节过程。糖酵解能力(Glycolytic potential，GP)是影响动物宰后肌肉pH变化的主要因素，这与宰前动物肌肉内的糖原含量、糖酵解酶活性以及糖酵解底物等有关，pH变化越大，说明糖酵解能力越强[36]。本研究结果显示，富集在糖酵解途径的差异基因PFKM、GPI、PKM和FBP2在小尾寒羊中的表达量显著低于滩羊和杜泊羊，说明小尾寒羊肉的pH值显著高于滩羊和杜泊羊，这与表型测定的结果相一致。李小金等[37]通过RNA-Seq技术研究发现，差异基因GPI、PGK1、PFKM和PGAM2在pH值较低的大约克猪中表达量显著高于pH值高的皖南花猪，说明这几个基因在调控肌肉糖酵解过程具有重要作用。段艳宇等[38]研究发现，猪肉的GP与pH24h呈显著负相关，与脂肪沉积呈正相关，这与检测的滩羊脂肪显著高于小尾寒羊，而pH24h显著低于小尾寒羊的结果相似。上述的这些分析结果都为初步揭示杜泊羊、滩羊和小尾寒羊肉质性状差异的遗传机理提供理论参考。

本研究利用转录组学技术对杜泊羊、滩羊和小尾寒羊背最长肌组织进行测序和分析，结果显示，共获得820个DEGs，其中杜泊羊与滩羊之间99个，杜泊羊与小尾寒羊之间436个，小尾寒羊与滩羊之间552个。通过对DEGs进行KEGG富集分析，发现主要富集在cAMP、MAPK、AMPK、嘌呤代谢、甘油磷脂代谢、氨基酸和脂肪酸代谢以及糖酵解/异生等信号通路，进一步筛选出FOS、PLA2G4E、LPIN1、AMPD1、AMPD3、NT5C1A、GPI、PFKM和PKM等与肉品质调控和风味物质代谢有关的候选基因，为宁夏优质肉羊新品种(系)在分子辅助育种方面提供基础资料。