基于轉(zhuǎn)錄組學(xué)篩選綿羊肉質(zhì)性狀相關(guān)候選基因

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

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

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

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

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

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

1 材料和方法

1.1 試驗(yàn)樣品采集

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

1.2 肉品質(zhì)指標(biāo)測(cè)定

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

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

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

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

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

1.5 基因可變剪切分析

使用rMATS(http：//rnaseq-mats.sourceforge.net/index.html)軟件對(duì)組織中基因的可變剪切結(jié)果進(jìn)行統(tǒng)計(jì)分析。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)篩選

通過(guò)StringTie分析流程計(jì)算每個(gè)轉(zhuǎn)錄本在各樣本中的表達(dá)水平，結(jié)果以FPKM(Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced)值表示。之后利用DESeq 2軟件進(jìn)行基因表達(dá)的差異顯著性分析，以padj(矯正后的Pvalue)<0.05和|log2FC|≥1作為組間差異顯著基因的篩選標(biāo)準(zhǔn)。

1.7 DEGs的GO/KEGG富集分析

為進(jìn)一步了解組間DEGs的功能和參與的代謝通路，將篩選出的DEGs進(jìn)行GO(http：//www.geneontology.org/)和KEGG(Kyoto Encyclopedia of Genes and Genomes)富集分析。

1.8 qRT-PCR驗(yàn)證

隨機(jī)選取13個(gè)差異表達(dá)基因(7個(gè)上調(diào)，6個(gè)下調(diào))進(jìn)行實(shí)時(shí)熒光定量PCR(Reverse transcription-quantitative PCR，qRT-PCR)，以驗(yàn)證轉(zhuǎn)錄組測(cè)序數(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個(gè)循環(huán)；65～95 ℃，每5 s增加0.5 ℃制作熔解曲線。每個(gè)樣品重復(fù)3次，利用2-ΔΔCt法計(jì)算各基因相對(duì)表達(dá)量，并通過(guò)Log2FoldChange標(biāo)準(zhǔn)化處理。

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

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

表1 引物信息Tab.1 Primers information

2 結(jié)果與分析

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

通過(guò)對(duì)綿羊的肉品質(zhì)指標(biāo)進(jìn)行測(cè)定，結(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 測(cè)序數(shù)據(jù)及質(zhì)量評(píng)估

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

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

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

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

2.4 基因可變剪切分析

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

2.5 DEGs的篩選和聚類分析

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

每個(gè)區(qū)域的盒形圖對(duì)應(yīng)5個(gè)統(tǒng)計(jì)量(從上而下分別為最大值、 上四分位數(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 差異表達(dá)基因的GO/KEGG富集分析

通過(guò)對(duì)DEGs進(jìn)行GO富集分析，可充分了解到DEGs在綿羊之間參與的生物學(xué)過(guò)程及生理功能。GO功能可分為三大類：生物學(xué)過(guò)程(Biological process，BP)、分子功能(Molecular function，MF)和細(xì)胞組分(Cellular component，CC)，本研究利用GOseq軟件對(duì)各組間DEGs進(jìn)行GO富集分析，結(jié)果見圖3所示：杜泊羊與灘羊組DEGs顯著富集到224個(gè)GO條目中，小尾寒羊與杜泊羊組DEGs顯著富集到517個(gè)GO條目中，小尾寒羊與灘羊組DEGs顯著富集到657個(gè)GO條目中。在生物學(xué)過(guò)程中，DEGs主要富集在代謝過(guò)程、有機(jī)物生物合成過(guò)程、 細(xì)胞對(duì)化學(xué)刺激的反應(yīng)、橫紋肌細(xì)胞分化、橫紋肌組織發(fā)育、肌細(xì)胞分化、肌肉組織發(fā)育和肌肉結(jié)構(gòu)發(fā)育等過(guò)程；在細(xì)胞組分中，主要富集在膜固有的、高分子絡(luò)合物、細(xì)胞內(nèi)非膜結(jié)合細(xì)胞器、細(xì)胞外區(qū)和核糖體等；在分子功能中，主要富集在RNA聚合酶Ⅱ調(diào)節(jié)區(qū)序列特異性DNA結(jié)合、結(jié)構(gòu)分子活性、核糖體的結(jié)構(gòu)成分等。

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

利用KOBAS軟件對(duì)組間DEGs進(jìn)行KEGG代謝途徑和信號(hào)通路分析，結(jié)果顯示，DEGs共富集到227條通路，其中杜泊羊與灘羊組DEGs富集到121條代謝通路中，小尾寒羊與杜泊羊組DEGs富集到197條代謝通路中，小尾寒羊與灘羊組DEGs富集到206條代謝通路中，3組DEGs共富集通路有104條。按P值從小到大選擇Top 20代謝通路，以散點(diǎn)圖展示(圖4)：只有1條代謝通路顯著富集(P<0.05)，為核糖體通路。剩下的19條以及其他代謝通路富集雖不顯著，但均有不同程度的DEGs富集，其中與肉品質(zhì)和風(fēng)味有關(guān)的通路及DEGs如表5所示。通過(guò)GO功能注釋、KEGG通路分析以及相關(guān)文獻(xiàn)報(bào)道，篩選到肌苷酸沉積(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.細(xì)胞對(duì)鎘離子的反應(yīng)；2.細(xì)胞對(duì)銨離子的反應(yīng)；3.含苯化合物代謝過(guò)程；4.骨吸收的正調(diào)控；5.骨重塑的正調(diào)控；6.中胚層細(xì)胞命運(yùn)承諾；7.對(duì)鎘離子的反應(yīng)；8.硝基苯代謝過(guò)程；9.氮化合物解毒；10.含氮化合物的細(xì)胞解毒；11.毒素生物合成過(guò)程；12.生物堿代謝過(guò)程；13.二苯并對(duì)二噁英代謝過(guò)程；14.氧化脫乙基；15.組織重塑的正調(diào)控；16.單萜類代謝過(guò)程；17.翻譯的晝夜調(diào)節(jié)；18.參與初級(jí)生殖層形成的細(xì)胞命運(yùn)承諾；19.中胚層細(xì)胞分化；20.異源分解代謝過(guò)程；21.AP1復(fù)合體；22.核糖體小亞基；23.胞體纖維；24.血紅蛋白復(fù)合物；25.生長(zhǎng)細(xì)胞尖端；26.新生長(zhǎng)細(xì)胞尖端；27.細(xì)胞尖端；28.精子纖維鞘；29.纖毛運(yùn)動(dòng)；30.嗜鉻顆粒膜；31.曼切特；32.嗜鉻顆粒；33.膜固有的；34.軸突丘；35.纖毛；36.核糖體亞單位；37.胞質(zhì)小核糖體亞單位；38.膜神經(jīng)元體細(xì)胞；39.氧化還原酶活性，作用于CH或CH2基團(tuán)、醌或類似化合物作為受體；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)運(yùn)體活性；53.細(xì)胞骨架的結(jié)構(gòu)成分；54.AMP脫氨酶活性；55.補(bǔ)體成分C1q結(jié)合；56.視蛋白結(jié)合；57.硫轉(zhuǎn)移酶活性；58.速激肽受體活性。B：1.翻譯；2.有機(jī)物生物合成過(guò)程；3.大分子生物合成過(guò)程；4.生物合成過(guò)程；5.細(xì)胞生物合成過(guò)程；6.細(xì)胞大分子生物合成過(guò)程；7.基因表達(dá)；8.細(xì)胞蛋白質(zhì)代謝過(guò)程；9.蛋白質(zhì)代謝過(guò)程；10.大分子代謝過(guò)程；11.代謝過(guò)程；12.蛋白質(zhì)折疊；13.核糖體；14.核糖核蛋白復(fù)合物；15.胞漿核糖體；16.核糖體亞單位；17.胞漿大核糖體亞單位；18.肌原纖維；19.收縮纖維；20.大核糖體亞單位；21.非膜結(jié)合細(xì)胞器；22.細(xì)胞內(nèi)非膜結(jié)合細(xì)胞器；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.細(xì)胞趨化性；2.細(xì)胞對(duì)化學(xué)刺激的反應(yīng)；3.蛋白質(zhì)折疊；4.橫紋肌細(xì)胞分化；5.肌細(xì)胞分化；6.呼吸爆發(fā)參與防御反應(yīng)；7. 翻譯；8.肌肉組織發(fā)育；9.肌肉結(jié)構(gòu)發(fā)育；10.橫紋肌組織發(fā)育；11.肌節(jié)；12.肌原纖維；13. 收縮纖維部分；14.收縮纖維；15.條紋；16.核糖體；17.胞漿核糖體；18.細(xì)胞外區(qū)；19.Z盤；20.胞漿部分；21.核糖體小亞基；22.未折疊蛋白結(jié)合；23.核糖體的結(jié)構(gòu)成分；24.結(jié)構(gòu)分子活性；25.rRNA結(jié)合；26.生長(zhǎng)因子活性。圖A代表D vs T組，圖B代表XHX vs D組，圖C代表XHX vs T組。圖4，5同。

A：1.膽汁分泌；2.色氨酸代謝；3.晝夜夾帶；4.卵巢類固醇生成；5.藥物代謝-細(xì)胞色素P450；6.糖胺聚糖生物合成-硫酸乙酰肝素/肝素；7.安非他明成癮；8.晝夜節(jié)律；9.cAMP信號(hào)通路；10.胰島素分泌；11.亞油酸代謝；12.雌激素信號(hào)通路；13.可卡因成癮；14.膽堿能突觸；15.ABC轉(zhuǎn)運(yùn)；16.血管平滑肌收縮；17.長(zhǎng)期抑郁癥；18.細(xì)胞色素P450對(duì)外源性物質(zhì)的代謝；19.多巴胺能突觸；20.類固醇激素生物合成。B：1.核糖體；2.雌激素信號(hào)通路；3.晝夜夾帶；4.胰島素信號(hào)通路；5.2-氧羰基酸代謝；6.膀胱癌；7.抗原處理和呈遞；8.縫隙連接；9.安非他明成癮；10.谷氨酸能突觸；11.淀粉和蔗糖代謝；12.軍團(tuán)菌病；13.慢性髓系白血??；14.酪氨酸代謝；15.p53信號(hào)通路；16.卵巢類固醇生成；17.酗酒；18.胰島素分泌；19.ErbB信號(hào)通路；20.精氨酸和脯氨酸代謝。C：1.核糖體；2.粘著；3.軍團(tuán)菌?。?.雌激素信號(hào)通路；5.膀胱癌；6.膽汁分泌；7.肌動(dòng)蛋白細(xì)胞骨架的調(diào)節(jié)；8.胃酸分泌；9.胰島素信號(hào)通路；10.谷氨酸能突觸；11.內(nèi)質(zhì)網(wǎng)中的蛋白質(zhì)加工；12.精氨酸和脯氨酸代謝；13.安非他明成癮；14.鈣信號(hào)通路；15.前列腺癌；16.ErbB信號(hào)通路；17.甲狀腺激素合成；18.瘧疾；19.嘌呤代謝；20.血小板活化。

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

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

通過(guò)qRT-PCR對(duì)隨機(jī)挑選的差異基因進(jìn)行驗(yàn)證，結(jié)果表明，各基因的定量結(jié)果與轉(zhuǎn)錄組測(cè)序結(jié)果的表達(dá)趨勢(shì)一致(圖5)，進(jìn)一步證實(shí)了轉(zhuǎn)錄組測(cè)序結(jié)果的可靠性。

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

3 結(jié)論與討論

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

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

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

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

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

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

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

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