梁 鵬，張 穩(wěn)，馮登偵，強 浩，榮 軒，孟 科

(寧夏大學(xué) 農(nóng)學(xué)院，寧夏 銀川 750021)

我國是綿羊養(yǎng)殖和羊肉消費大國，綿羊種類占世界綿羊品種的10%左右[1]。經(jīng)過綿羊育種工作者多年的努力和研究，我國的綿羊良種繁育體系也在不斷完善。同時，隨著經(jīng)濟的快速發(fā)展和人民生活水平的提高，人們對肉類消費的理念已經(jīng)從單純追求數(shù)量轉(zhuǎn)變?yōu)閿?shù)量和質(zhì)量并重。因此，人們也把綿羊育種的目標逐漸轉(zhuǎn)向為培育具有多優(yōu)良性狀的肉羊新品種(系)，以滿足市場需求的供應(yīng)。研究表明，品種是影響畜禽肉品質(zhì)的主要因素[2]，也是種質(zhì)基因庫的重要保存單位，一些具有特色的畜禽品種便成為理想的研究對象。

灘羊是在寧夏地區(qū)經(jīng)過漫長的自然和人工選擇形成的優(yōu)良特色畜種，其肉質(zhì)以鮮美細嫩、營養(yǎng)豐富和風(fēng)味獨特而聞名，是世界認可的優(yōu)質(zhì)羊肉之一[3]，但灘羊個體小、早期生長發(fā)育慢、產(chǎn)肉力不足以及繁殖力低下，羊肉總產(chǎn)量只能靠多飼養(yǎng)、高出欄為主，這使得灘羊養(yǎng)殖的成本增加，且養(yǎng)殖效益低，也不符合當前的產(chǎn)業(yè)發(fā)展要求。因此，為了達到對灘羊這一種質(zhì)資源的創(chuàng)新與利用，選用產(chǎn)肉性能好的杜泊羊和繁殖性能高的小尾寒羊與本地灘羊進行三元雜交，以期利用各親本的優(yōu)良性狀培育出適合寧夏地區(qū)資源環(huán)境特點的肉羊新品種(系)。

利用傳統(tǒng)方法對肉品質(zhì)和風(fēng)味的研究只能通過檢測肉嫩度、脂肪酸和氨基酸含量等進行分析，而轉(zhuǎn)錄組測序技術(shù)能夠深入解析與這些表型性狀有關(guān)的調(diào)控因子，從分子水平上對調(diào)控肉品質(zhì)和風(fēng)味性狀的相關(guān)基因進行篩選和鑒定。目前，轉(zhuǎn)錄組測序技術(shù)已經(jīng)廣泛應(yīng)用于畜禽功能性狀基因的挖掘與育種研究中，Wang等[4]利用轉(zhuǎn)錄組學(xué)技術(shù)研究了影響晉南牛與西門塔爾牛肌內(nèi)脂肪(IMF)含量差異的相關(guān)調(diào)控基因，結(jié)果篩選到的這些差異基因(CYP21A2、PC、ACACB、APOA1和FADS2)可能通過調(diào)節(jié)脂質(zhì)代謝來影響IMF的沉積。Sun等[5]利用轉(zhuǎn)錄組學(xué)技術(shù)研究了千華肉用美利奴羊和小尾寒羊之間與肌肉生長發(fā)育有關(guān)的轉(zhuǎn)錄調(diào)節(jié)因子，結(jié)果顯示，960個基因在品種間差異表達，通過代謝途徑鑒定到MRF、GXP1和STAC3等差異基因在肌肉生長發(fā)育過程中起著重要作用。韓信嶸等[6]以背膘組織為研究對象，利用轉(zhuǎn)錄組學(xué)技術(shù)分析安格斯牛和西門塔爾牛與脂肪沉積有關(guān)的差異表達基因，結(jié)果篩選到DKKL1、ACACB和PTGS2基因通過cAMP通路調(diào)控脂肪沉積。

本研究以寧夏優(yōu)質(zhì)肉羊培育中的3個親本為研究對象，旨在通過轉(zhuǎn)錄組學(xué)技術(shù)挖掘影響不同品種綿羊肉品質(zhì)和風(fēng)味性狀差異的相關(guān)調(diào)控基因，為寧夏優(yōu)質(zhì)肉羊新品種(系)培育提供理論依據(jù)。

1 材料和方法

1.1 試驗樣品采集

選擇相同飼養(yǎng)環(huán)境下的8月齡杜泊羊、灘羊和小尾寒羊各8只(公母各1/2)作為研究對象，試驗動物均飼養(yǎng)于寧夏宇泊科技有限公司威震肉羊繁育場，屠宰前禁食24 h、停水2 h，于頸靜脈放血屠宰后，迅速采集第3～12肋骨之間的背最長肌組織約500 g，剔除表面筋膜后裝入自封袋中并標記，用于肉品質(zhì)指標測定。另取少部分樣品裝入無RNA酶的凍存管中，標號后快速置于液氮罐中，待試驗結(jié)束后將樣品帶回實驗室，于-80 ℃冰箱保存?zhèn)溆谩?/p>

1.2 肉品質(zhì)指標測定

采用石蠟切片法觀察4%多聚甲醛固定樣品的肌纖維直徑和密度；使用pH計測定屠宰24 h后的pH值；使用嫩度儀和壓力法分別測定剪切力和系水力；按照GB/T9695.7—2008《肉及肉制品中脂肪含量測定》測定肌肉中的脂肪含量；采用高效液相色譜法測定肌肉中肌苷酸和肌苷含量。

1.3 文庫構(gòu)建及轉(zhuǎn)錄組測序

隨機選取杜泊羊、灘羊和小尾寒羊母羊各3只，提取背最長肌組織總RNA，利用瓊脂糖凝膠電泳、Nanodrop和Agilent 2100 bioanalyzer檢測RNA樣品的濃度、純度及完整性。經(jīng)質(zhì)檢合格后采用去除核糖體RNA的方法構(gòu)建鏈特異性文庫，文庫構(gòu)建完成后，先使用Qubit 2.0進行初步定量，稀釋文庫至1.5 ng/μL，隨后使用Agilent 2100對文庫的插入片段大小(Insert size)進行檢測，符合預(yù)期后，使用qRT-PCR方法對文庫有效濃度(高于3 nmol/L)進行準確定量，以保證文庫質(zhì)量。文庫經(jīng)檢驗合格后，利用Illumina HiSeqTM4000平臺進行PE150雙端測序。

1.4 轉(zhuǎn)錄本拼接、篩選及預(yù)測

對測序獲得的原始數(shù)據(jù)進行過濾、錯誤率及GC含量分布檢查，獲得后續(xù)分析使用的高質(zhì)量數(shù)據(jù)(Clean reads)。利用Hisat 2(http：//ccb.jhu.edu/software/hisat2)軟件將Clean reads數(shù)據(jù)比對到綿羊參考基因組上，然后使用Stringtie軟件在基于比對到基因組上的結(jié)果，將reads拼接成轉(zhuǎn)錄本進行定量。使用Cuffmerge軟件對各樣本拼接得到的轉(zhuǎn)錄本進行合并，去掉其中鏈方向不確定和轉(zhuǎn)錄本長度不超過200 nt的轉(zhuǎn)錄本，之后利用Cuffcompare軟件和已知數(shù)據(jù)庫進行比較，過濾掉數(shù)據(jù)庫中已知的轉(zhuǎn)錄本，最后對篩選到的新轉(zhuǎn)錄本進行編碼潛能預(yù)測，預(yù)測結(jié)果取CPC、PFAM和CNCI的交集。

1.5 基因可變剪切分析

使用rMATS(http：//rnaseq-mats.sourceforge.net/index.html)軟件對組織中基因的可變剪切結(jié)果進行統(tǒng)計分析。rMATS分析的可變剪切事件有外顯子跳躍(Skipped exon，SE)、外顯子互斥(Mutually exclusive exon，MXE)、5′端可變剪切(Alternative 5′ splice site，A5SS)、3′端可變剪切(Alternative 3′ splice site，A3SS)和內(nèi)含子滯留(Retained intron，RI)這5種。

1.6 組間差異基因(DEGs)篩選

通過StringTie分析流程計算每個轉(zhuǎn)錄本在各樣本中的表達水平，結(jié)果以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個上調(diào)，6個下調(diào))進行實時熒光定量PCR(Reverse transcription-quantitative PCR，qRT-PCR)，以驗證轉(zhuǎn)錄組測序數(shù)據(jù)的可靠性。利用FastKing RT Kit(With gDNase)試劑盒(天根)將原始樣品的RNA反轉(zhuǎn)錄得到cDNA。內(nèi)參基因選用甘油醛-3-磷酸脫氫酶基因(Glyceraldehyde-3-phosphate dehydrogenase gene，GAPDH)，引物均由中科羽瞳(陜西)生物科技有限公司合成，序列信息見表1。熒光定量PCR的反應(yīng)體系為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。反應(yīng)條件為：95 ℃ 30 s； 95 ℃ 10 s，退火30 s，40個循環(huán)；65～95 ℃，每5 s增加0.5 ℃制作熔解曲線。每個樣品重復(fù)3次，利用2-ΔΔCt法計算各基因相對表達量，并通過Log2FoldChange標準化處理。

1.9 數(shù)據(jù)處理

利用SPSS 22.0軟件對肉品質(zhì)指標進行單因素方差分析，用Duncan法進行多重比較，數(shù)據(jù)結(jié)果以“平均值±標準誤”表示。以P<0.05表示組間差異顯著。

表1 引物信息Tab.1 Primers information

2 結(jié)果與分析

2.1 肉品質(zhì)差異分析

通過對綿羊的肉品質(zhì)指標進行測定，結(jié)果顯示(表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 綿羊肉品質(zhì)差異分析Tab.2 Analysis of sheep meat quality difference

2.2 測序數(shù)據(jù)及質(zhì)量評估

通過建立9個背最長肌組織樣品的文庫，上機測序后數(shù)據(jù)經(jīng)過過濾、錯誤率和GC含量分布檢查，共獲得808 043 288條有效序列(表3)，有效堿基達到121.21 G。各樣品的GC含量在47.12%～52.04%，Q20≥97.63%，Q30≥93.49%，錯誤率≤0.03%；測序的有效數(shù)據(jù)通過與參考基因組比對，結(jié)果顯示比對率均在89.25%以上。這些結(jié)果都表明，測序數(shù)據(jù)質(zhì)量可靠，可用于后續(xù)的分析。

表3 測序數(shù)據(jù)質(zhì)量評估Tab.3 Sequencing data quality assessment

2.3 背最長肌組織轉(zhuǎn)錄本篩選與預(yù)測

基于比對結(jié)果，將reads拼接成轉(zhuǎn)錄本，共得到952 263個轉(zhuǎn)錄本。通過將各樣本拼接的轉(zhuǎn)錄本進行合并，且除去不符合要求的轉(zhuǎn)錄本后，將剩余轉(zhuǎn)錄本與已知數(shù)據(jù)庫進行比較，過濾掉數(shù)據(jù)庫中已知的轉(zhuǎn)錄本，得到新轉(zhuǎn)錄本8 083個，其中位于“+”鏈的新轉(zhuǎn)錄本3 992個，位于“-”鏈的新轉(zhuǎn)錄本4 091個。新轉(zhuǎn)錄本的長度為201～71 939 bp，外顯子個數(shù)為2～35個，開放閱讀框(Open reading frame，ORF)在29～1 807 bp。經(jīng)編碼潛能預(yù)測，結(jié)果顯示，新轉(zhuǎn)錄本中280個具有編碼潛能，可作為Novel-mRNA數(shù)據(jù)集。

2.4 基因可變剪切分析

通過rMATS軟件對各組織樣品中的可變剪切進行定量，結(jié)果見表4。由表4可知，各比較組中檢測到的可變剪切事件分別為18 239，18 258，19 222個，說明許多基因都存在多種剪切方式，造成基因在表達水平上發(fā)生變化，進而使生物體能夠表達多種蛋白，引起物種的表型性狀產(chǎn)生差異。在5種可變剪切事件中，各比較組合均以SE事件最多，MXE事件次之，A5SS事件最少。

2.5 DEGs的篩選和聚類分析

通過FPKM方法計算各樣本中基因的表達水平，結(jié)果以盒形圖展示，由圖1可知，不同品種綿羊背最長肌組織樣品中基因的表達水平有差異。通過計算各樣品間皮爾遜相關(guān)系數(shù)(Pearson′s correlation coefficient)，可以說明重復(fù)樣品的相關(guān)性強弱，r越接近1，表明樣品相關(guān)性越強。本研究中各樣品的r值均大于0.95，說明試驗的可靠性和組內(nèi)樣品的重復(fù)性較高。使用DESeq 2軟件對定量分析完成的樣品進行組間基因表達差異顯著性分析，結(jié)果顯示，共鑒定出820個DEGs。杜泊羊與灘羊組鑒定出99個DEGs，其中51個上調(diào)，48個下調(diào)；小尾寒羊與杜泊羊組鑒定出436個DEGs，其中270個上調(diào)，166個下調(diào)；小尾寒羊與灘羊組鑒定出552個DEGs，其中320個上調(diào)，232個下調(diào)。3個組共表達差異顯著基因有4個(MAOB、SMAD1、AMPD3和ENAH)。對篩選出來的DEGS進行聚類分析(圖2)，結(jié)果顯示，同一品種的3個生物學(xué)重復(fù)樣本能夠很好地聚到一起，說明本研究所用樣本的生物學(xué)重復(fù)較好。小尾寒羊與杜泊羊和灘羊的基因表達模式出現(xiàn)明顯分離，說明組間差異較大，而杜泊羊和灘羊重疊區(qū)域較多，說明組間差異較小。

每個區(qū)域的盒形圖對應(yīng)5個統(tǒng)計量(從上而下分別為最大值、 上四分位數(shù)、中值、下四分位數(shù)和最小值)。 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在綿羊之間參與的生物學(xué)過程及生理功能。GO功能可分為三大類：生物學(xué)過程(Biological process，BP)、分子功能(Molecular function，MF)和細胞組分(Cellular component，CC)，本研究利用GOseq軟件對各組間DEGs進行GO富集分析，結(jié)果見圖3所示：杜泊羊與灘羊組DEGs顯著富集到224個GO條目中，小尾寒羊與杜泊羊組DEGs顯著富集到517個GO條目中，小尾寒羊與灘羊組DEGs顯著富集到657個GO條目中。在生物學(xué)過程中，DEGs主要富集在代謝過程、有機物生物合成過程、 細胞對化學(xué)刺激的反應(yīng)、橫紋肌細胞分化、橫紋肌組織發(fā)育、肌細胞分化、肌肉組織發(fā)育和肌肉結(jié)構(gòu)發(fā)育等過程；在細胞組分中，主要富集在膜固有的、高分子絡(luò)合物、細胞內(nèi)非膜結(jié)合細胞器、細胞外區(qū)和核糖體等；在分子功能中，主要富集在RNA聚合酶Ⅱ調(diào)節(jié)區(qū)序列特異性DNA結(jié)合、結(jié)構(gòu)分子活性、核糖體的結(jié)構(gòu)成分等。

左側(cè)根據(jù)表達相似程度對基因進行聚類，上方根據(jù)表達譜的相似程度對每個樣本進行聚類，由藍至紅表達量逐漸上調(diào)，數(shù)字為均一化后的相對表達量。

利用KOBAS軟件對組間DEGs進行KEGG代謝途徑和信號通路分析，結(jié)果顯示，DEGs共富集到227條通路，其中杜泊羊與灘羊組DEGs富集到121條代謝通路中，小尾寒羊與杜泊羊組DEGs富集到197條代謝通路中，小尾寒羊與灘羊組DEGs富集到206條代謝通路中，3組DEGs共富集通路有104條。按P值從小到大選擇Top 20代謝通路，以散點圖展示(圖4)：只有1條代謝通路顯著富集(P<0.05)，為核糖體通路。剩下的19條以及其他代謝通路富集雖不顯著，但均有不同程度的DEGs富集，其中與肉品質(zhì)和風(fēng)味有關(guān)的通路及DEGs如表5所示。通過GO功能注釋、KEGG通路分析以及相關(guān)文獻報道，篩選到肌苷酸沉積(AMPD1、AMPD3、NT5C1A)、脂肪沉積(LPIN1、FABP4、ACOT7、DUSP1、GOT1、CKMT2、CAPN2、SLC2A4、SCOS3、RXRG、PPARGC1A)、肌纖維調(diào)控(FOS、ACTC1、GYS1、CS)、肌肉糖酵解/糖異生(PFKM、PKM、GPI、PGAM1、FBP2)、肌肉發(fā)育(CKM、SMAD1)和脂質(zhì)代謝(PIK3C2G、ABCA1)等與肉品質(zhì)和風(fēng)味有關(guān)的基因。

A：1.細胞對鎘離子的反應(yīng)；2.細胞對銨離子的反應(yīng)；3.含苯化合物代謝過程；4.骨吸收的正調(diào)控；5.骨重塑的正調(diào)控；6.中胚層細胞命運承諾；7.對鎘離子的反應(yīng)；8.硝基苯代謝過程；9.氮化合物解毒；10.含氮化合物的細胞解毒；11.毒素生物合成過程；12.生物堿代謝過程；13.二苯并對二噁英代謝過程；14.氧化脫乙基；15.組織重塑的正調(diào)控；16.單萜類代謝過程；17.翻譯的晝夜調(diào)節(jié)；18.參與初級生殖層形成的細胞命運承諾；19.中胚層細胞分化；20.異源分解代謝過程；21.AP1復(fù)合體；22.核糖體小亞基；23.胞體纖維；24.血紅蛋白復(fù)合物；25.生長細胞尖端；26.新生長細胞尖端；27.細胞尖端；28.精子纖維鞘；29.纖毛運動；30.嗜鉻顆粒膜；31.曼切特；32.嗜鉻顆粒；33.膜固有的；34.軸突丘；35.纖毛；36.核糖體亞單位；37.胞質(zhì)小核糖體亞單位；38.膜神經(jīng)元體細胞；39.氧化還原酶活性，作用于CH或CH2基團、醌或類似化合物作為受體；40.咖啡因氧化酶活性；41.芳香化酶活性；42.轉(zhuǎn)鐵蛋白受體活性；43.吡喃結(jié)構(gòu)域結(jié)合；44.貨物受體活性；45.犬尿氨酸3-單加氧酶活性；46. UDP-N-乙酰氨基葡萄糖-溶酶體-酶N-乙酰氨基葡萄糖磷酸轉(zhuǎn)移酶活性；47.[硫酸乙酰肝素]-氨基葡萄糖3-磺基轉(zhuǎn)移酶1活性；48. P物質(zhì)受體活性；49.鈣和鈣調(diào)素反應(yīng)性腺苷酸環(huán)化酶活性；50.RNA聚合酶Ⅱ調(diào)節(jié)區(qū)序列特異性DNA結(jié)合；51.RNA聚合酶Ⅱ調(diào)節(jié)區(qū)DNA結(jié)合；52.氧轉(zhuǎn)運體活性；53.細胞骨架的結(jié)構(gòu)成分；54.AMP脫氨酶活性；55.補體成分C1q結(jié)合；56.視蛋白結(jié)合；57.硫轉(zhuǎn)移酶活性；58.速激肽受體活性。B：1.翻譯；2.有機物生物合成過程；3.大分子生物合成過程；4.生物合成過程；5.細胞生物合成過程；6.細胞大分子生物合成過程；7.基因表達；8.細胞蛋白質(zhì)代謝過程；9.蛋白質(zhì)代謝過程；10.大分子代謝過程；11.代謝過程；12.蛋白質(zhì)折疊；13.核糖體；14.核糖核蛋白復(fù)合物；15.胞漿核糖體；16.核糖體亞單位；17.胞漿大核糖體亞單位；18.肌原纖維；19.收縮纖維；20.大核糖體亞單位；21.非膜結(jié)合細胞器；22.細胞內(nèi)非膜結(jié)合細胞器；23.肌節(jié)；24.胞漿部分；25.高分子絡(luò)合物；26.收縮纖維部分；27.胞質(zhì)小核糖體亞單位；28.多體核糖體；29.核糖體小亞基；30.條紋；31.核糖體的結(jié)構(gòu)成分；32.結(jié)構(gòu)分子活性；33.未折疊蛋白結(jié)合；34.rRNA結(jié)合。C：1.細胞趨化性；2.細胞對化學(xué)刺激的反應(yīng)；3.蛋白質(zhì)折疊；4.橫紋肌細胞分化；5.肌細胞分化；6.呼吸爆發(fā)參與防御反應(yīng)；7. 翻譯；8.肌肉組織發(fā)育；9.肌肉結(jié)構(gòu)發(fā)育；10.橫紋肌組織發(fā)育；11.肌節(jié)；12.肌原纖維；13. 收縮纖維部分；14.收縮纖維；15.條紋；16.核糖體；17.胞漿核糖體；18.細胞外區(qū)；19.Z盤；20.胞漿部分；21.核糖體小亞基；22.未折疊蛋白結(jié)合；23.核糖體的結(jié)構(gòu)成分；24.結(jié)構(gòu)分子活性；25.rRNA結(jié)合；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.晝夜節(jié)律；9.cAMP信號通路；10.胰島素分泌；11.亞油酸代謝；12.雌激素信號通路；13.可卡因成癮；14.膽堿能突觸；15.ABC轉(zhuǎn)運；16.血管平滑肌收縮；17.長期抑郁癥；18.細胞色素P450對外源性物質(zhì)的代謝；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.肌動蛋白細胞骨架的調(diào)節(jié)；8.胃酸分泌；9.胰島素信號通路；10.谷氨酸能突觸；11.內(nèi)質(zhì)網(wǎng)中的蛋白質(zhì)加工；12.精氨酸和脯氨酸代謝；13.安非他明成癮；14.鈣信號通路；15.前列腺癌；16.ErbB信號通路；17.甲狀腺激素合成；18.瘧疾；19.嘌呤代謝；20.血小板活化。

表5 部分KEGG通路與DEGs富集結(jié)果Tab.5 Results of partial enrichment of KEGG pathway with DEGs

2.7 qRT-PCR驗證結(jié)果

通過qRT-PCR對隨機挑選的差異基因進行驗證，結(jié)果表明，各基因的定量結(jié)果與轉(zhuǎn)錄組測序結(jié)果的表達趨勢一致(圖5)，進一步證實了轉(zhuǎn)錄組測序結(jié)果的可靠性。

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

3 結(jié)論與討論

肉品質(zhì)和風(fēng)味是消費者評價畜禽肉質(zhì)好壞的重要標準，也是畜禽育種改良的重要經(jīng)濟指標之一。肉的風(fēng)味主要與脂肪、氨基酸、核苷酸和有機酸等前體物質(zhì)密切相關(guān)，在肌肉中它們含量的不同會導(dǎo)致肉品質(zhì)和風(fēng)味產(chǎn)生差異[7-8]。通過對杜泊羊、灘羊和小尾寒羊的肉品質(zhì)指標進行測定，結(jié)果顯示，品種間存在不同程度的顯著差異。因此，為了解造成這些肉品質(zhì)指標產(chǎn)生差異的潛在遺傳機制，進一步采用轉(zhuǎn)錄組學(xué)技術(shù)對3個親本群體進行高通量測序，篩選出組間差異表達基因及富集的通路進一步比較物質(zhì)沉積能力的不同，有助于后期的綿羊分子輔助育種，加速寧夏優(yōu)質(zhì)肉羊新品種(系)的育成。

本研究通過對3個親本綿羊群體的背最長肌組織進行轉(zhuǎn)錄組測序，以|log2FC|≥1和P<0.05為標準篩選組間DEGs，結(jié)果顯示，共有820個DEGs可用于后續(xù)的生物信息學(xué)分析，其中杜泊羊與灘羊組99個，小尾寒羊與杜泊羊組436個，小尾寒羊與灘羊組552個，說明杜泊羊與灘羊差異較小，與小尾寒羊差異較大，這與檢測肉品質(zhì)指標分析的結(jié)果相似，進一步說明測序結(jié)果的準確性。

通過對DEGs進行GO功能注釋和KEGG通路富集分析，能夠深入了解基因的生物學(xué)功能和參與的代謝通路。本研究篩選出的DEGs富集到與肉品質(zhì)調(diào)控和風(fēng)味物質(zhì)代謝有關(guān)的通路，如cAMP、TGF-β、MAPK和AMPK等經(jīng)典信號通路以及嘌呤代謝、甘油磷脂代謝、氨基酸和脂肪酸代謝、肌動蛋白細胞骨架的調(diào)節(jié)和糖酵解/糖異生等通路。研究表明，cAMP信號通路通過與其下游的反應(yīng)元件(CREB)和調(diào)控轉(zhuǎn)錄的輔助激活因子(CRTC)發(fā)揮作用，當肌肉中過表達CRTC2時會增加肌內(nèi)脂肪含量和肌纖維的橫截面積[9]；當過表達CRTC3時會促進甘油?；D(zhuǎn)移酶的活性，進而增加肌肉中甘油三酯的含量[10]；在敲除小鼠的CRTC3基因后，其分解脂肪的能力增強[11]。在3T3-L1細胞分化過程中，cAMP通路通過活化CREB和C/EBPβ來促進脂肪的生成[12]。本研究結(jié)果表明，cAMP信號通路差異基因較多，其中FOS在調(diào)節(jié)細胞增殖、分化和凋亡過程中具有重要作用，能夠作為啟動基因轉(zhuǎn)錄的樞紐，通過細胞的刺激和應(yīng)激調(diào)節(jié)基因表達[13]。成志敏等[14]研究表明，EGR1基因能夠通過促進FOS基因的表達，引起大白豬肌細胞增長，導(dǎo)致肌纖維直徑增大，肌纖維密度降低，本研究結(jié)果也發(fā)現(xiàn)，EGR1和FOS基因在小尾寒羊中均高表達，且EGR1基因顯著高于杜泊羊和灘羊，而FOS基因在杜泊羊中的表達顯著低于小尾寒羊和灘羊，表型測定結(jié)果顯示，杜泊羊肌纖維密度顯著低于灘羊和小尾寒羊，說明EGR1和FOS基因在調(diào)控杜泊羊肌纖維特性方面具有重要作用。同時，F(xiàn)OS基因也富集在MAPK信號通路，而MAPK通路介導(dǎo)細胞對外界做出反應(yīng)的重要信號系統(tǒng)之一，不僅對脂肪細胞的增殖和分化具有重要作用[15]，在肌纖維發(fā)育過程也具有重要作用，通過調(diào)控肌纖維類型來影響肉質(zhì)性狀[16-17]。魏立民等[18]在研究五指山豬和長白豬背最長肌組織差異基因時，發(fā)現(xiàn)FOS基因通過MAPK信號通路調(diào)控肌肉發(fā)育。

AMPK通路是機體能量代謝的重要信號通路，能夠為肌肉發(fā)育提供能量。過氧化物酶體增殖物激活受體γ輔激活因子1A(PPARGC1A)可以與不同的轉(zhuǎn)錄因子結(jié)合，參與多個代謝過程，如肌纖維轉(zhuǎn)化和脂肪酸代謝等[19]。研究表明，PPARGC1A基因在豬肌肉中的表達與肌纖維面積和脂肪細胞體積呈負相關(guān)，而與脂肪組織中的表達呈正相關(guān)[20]。劉玉芳等[21]研究表明，PPARGC1A的表達量在脂肪沉積較高的金華豬背脂中顯著低于瘦肉型大白豬，進一步推測該基因表達與脂肪沉積呈負相關(guān)。本研究中該基因在脂肪沉積較高的灘羊中表達量高于杜泊羊和小尾寒羊，表明該基因與脂肪沉積呈正相關(guān)，這與上述研究結(jié)果相反，可能是與研究的物種不同有關(guān)。而劉遠等[22]通過轉(zhuǎn)錄組學(xué)技術(shù)研究表明，PPARGC1A基因在脂肪沉積高的福清山羊中表達量高于低脂肪的努比亞黑山羊，說明該基因的表達具有物種特異性。

甘油磷脂代謝通路中產(chǎn)生的脂肪酸和溶血磷脂等可進一步合成磷脂，而磷脂是肉風(fēng)味形成的重要前體物質(zhì)。該通路中的GPAM可催化動物脂質(zhì)代謝中甘油酰基的合成過程，進而影響三酰甘油和磷脂的合成[23]。氨基酸代謝、脂肪酸代謝和肌醇磷酸代謝通路均與動物脂質(zhì)代謝有關(guān)，這幾個通路中的基因如PLA2G4E[24]、LPIN1[25]、CS[26]、SOCS3[27]等均報道與脂代謝有關(guān)，可能在調(diào)控綿羊肉品質(zhì)和風(fēng)味方面發(fā)揮作用，在今后還需進一步研究。

肌苷酸(IMP)又被稱作次黃嘌呤核苷酸，是反映肉質(zhì)鮮味的重要指標，它是機體在嘌呤代謝過程中合成的物質(zhì)。研究表明，腺苷單磷酸脫氨酶(AMPD)是影響嘌呤核苷酸代謝的限速酶，與肌肉中的IMP代謝有關(guān)，它包括AMPD1、AMPD2和AMPD3共 3個亞型[28]。IMP在體內(nèi)合成包括從頭合成和補救合成2條途徑，從頭合成過程有10種酶參與反應(yīng)，而AMPD1是催化AMP水解脫氨生成IMP的酶，主要在骨骼肌中表達[29]。多項研究也表明，AMPD1基因表達量與肌肉中IMP含量密切相關(guān)[30]，主要在禽類動物上報道較多。張會豐等[31]研究了城口山地雞、大寧河雞和青腳麻雞胸肌中IMP含量與AMPD1基因表達量的關(guān)系，結(jié)果顯示二者呈正相關(guān)；而陳星等[32]通過轉(zhuǎn)錄組學(xué)技術(shù)表明，AMPD1在肌苷酸含量低的沔陽麻鴨中表達量顯著高于連城白鴨，二者差異倍數(shù)達到3.73倍，qRT-PCR進一步驗證的結(jié)果也顯示差異倍數(shù)達到2.11倍。本研究結(jié)果顯示，AMPD1在肌苷酸含量最高的小尾寒羊中表達量最低，顯著低于肌苷酸含量最低的杜泊羊，這與羅玉龍等[33]在蘇尼特羊中的研究結(jié)果相反，推測AMPD1有可能與IMP沉積呈負相關(guān)，這與所研究的物種或品種不同有關(guān)，也可能與其他因素有關(guān)，具體還有待進一步探究。AMPD3和NT5C1A也富集在嘌呤代謝通路，有研究報道，這2個基因也與IMP代謝途徑有關(guān)，因此可用于綿羊IMP沉積的候選基因進行研究[34-35]。

糖酵解/異生代謝途徑是影響畜禽肉品質(zhì)的重要調(diào)節(jié)過程。糖酵解能力(Glycolytic potential，GP)是影響動物宰后肌肉pH變化的主要因素，這與宰前動物肌肉內(nèi)的糖原含量、糖酵解酶活性以及糖酵解底物等有關(guān)，pH變化越大，說明糖酵解能力越強[36]。本研究結(jié)果顯示，富集在糖酵解途徑的差異基因PFKM、GPI、PKM和FBP2在小尾寒羊中的表達量顯著低于灘羊和杜泊羊，說明小尾寒羊肉的pH值顯著高于灘羊和杜泊羊，這與表型測定的結(jié)果相一致。李小金等[37]通過RNA-Seq技術(shù)研究發(fā)現(xiàn)，差異基因GPI、PGK1、PFKM和PGAM2在pH值較低的大約克豬中表達量顯著高于pH值高的皖南花豬，說明這幾個基因在調(diào)控肌肉糖酵解過程具有重要作用。段艷宇等[38]研究發(fā)現(xiàn)，豬肉的GP與pH24h呈顯著負相關(guān)，與脂肪沉積呈正相關(guān)，這與檢測的灘羊脂肪顯著高于小尾寒羊，而pH24h顯著低于小尾寒羊的結(jié)果相似。上述的這些分析結(jié)果都為初步揭示杜泊羊、灘羊和小尾寒羊肉質(zhì)性狀差異的遺傳機理提供理論參考。

本研究利用轉(zhuǎn)錄組學(xué)技術(shù)對杜泊羊、灘羊和小尾寒羊背最長肌組織進行測序和分析，結(jié)果顯示，共獲得820個DEGs，其中杜泊羊與灘羊之間99個，杜泊羊與小尾寒羊之間436個，小尾寒羊與灘羊之間552個。通過對DEGs進行KEGG富集分析，發(fā)現(xiàn)主要富集在cAMP、MAPK、AMPK、嘌呤代謝、甘油磷脂代謝、氨基酸和脂肪酸代謝以及糖酵解/異生等信號通路，進一步篩選出FOS、PLA2G4E、LPIN1、AMPD1、AMPD3、NT5C1A、GPI、PFKM和PKM等與肉品質(zhì)調(diào)控和風(fēng)味物質(zhì)代謝有關(guān)的候選基因，為寧夏優(yōu)質(zhì)肉羊新品種(系)在分子輔助育種方面提供基礎(chǔ)資料。