基于高通量測(cè)序技術(shù)鑒定松遼黑豬與長(zhǎng)白豬背最長(zhǎng)肌差異基因

侍玉梅，陳少康，邢 凱，趙延輝，原佳妮，盛熙暉，齊曉龍，倪和民，郭 勇，王楚端

(1.北京農(nóng)學(xué)院 動(dòng)物科學(xué)技術(shù)學(xué)院，北京 102206；2.北京市畜牧總站，北京 100101；3.中國(guó)農(nóng)業(yè)大學(xué) 動(dòng)物科學(xué)技術(shù)學(xué)院，北京 100193)

肌肉生長(zhǎng)和脂肪沉積是養(yǎng)豬生產(chǎn)中非常重要的經(jīng)濟(jì)性狀，同時(shí)也是復(fù)雜的數(shù)量性狀，受到眾多因素的影響[1]；然而，豬肌肉生長(zhǎng)和脂肪沉積過(guò)程中的分子調(diào)控機(jī)制仍有待進(jìn)一步研究。豬的脂肪性狀直接影響豬肉品質(zhì)，而對(duì)不同品種豬的肉質(zhì)進(jìn)行比較分析是鑒定差異表達(dá)基因(DEGs)的主要方法之一。松遼黑豬作為我國(guó)培育的地方品種，具有繁殖率高、肉質(zhì)優(yōu)良、適應(yīng)性強(qiáng)、耐粗飼等特性[2]。長(zhǎng)白豬作為西方豬種，具有生長(zhǎng)速度快、飼料利用率高、瘦肉率高等特點(diǎn)，但肉質(zhì)欠佳[3]。長(zhǎng)白豬是典型的瘦肉型豬，而松遼黑豬是我國(guó)特有的脂肪型豬，2個(gè)不同品種豬的表型差異很大，是鑒別肌肉和脂肪中DEGs的良好動(dòng)物模型。

高通量測(cè)序使轉(zhuǎn)錄組研究變得容易進(jìn)行，以轉(zhuǎn)錄組測(cè)序(RNA-seq)為代表的新一代高通量測(cè)序技術(shù)能夠全面快速獲取某一物種特定器官或組織在某一狀態(tài)下的幾乎所有轉(zhuǎn)錄本信息。該技術(shù)既可以在有基因組注釋的情況下應(yīng)用以促進(jìn)轉(zhuǎn)錄本識(shí)別，也可以在沒有參考基因組的情況下應(yīng)用，其被廣泛應(yīng)用于轉(zhuǎn)錄組分析。王志秀[4]對(duì)滇南小耳豬、藏豬、長(zhǎng)白豬和杜洛克豬進(jìn)行高通量測(cè)序發(fā)現(xiàn)，在滇南小耳豬-藏豬組和長(zhǎng)白豬-杜洛克豬組共篩選出315個(gè)差異表達(dá)基因，其中，有140個(gè)上調(diào)基因，175個(gè)下調(diào)基因，并鑒定出27個(gè)與脂肪沉積相關(guān)的基因，其主要參與脂肪代謝、脂肪酸合成等過(guò)程。Xu等[5]對(duì)不同日齡梅山豬的骨骼肌進(jìn)行轉(zhuǎn)錄組研究，鑒定出338個(gè)DEGs，蛋白質(zhì)譜分析發(fā)現(xiàn)66個(gè)差異表達(dá)基因，功能分析顯示其功能主要為代謝、肌原纖維細(xì)絲形成、細(xì)胞骨架形成、收縮活動(dòng)和信號(hào)轉(zhuǎn)導(dǎo)。NOISeq是用于分析RNA-seq數(shù)據(jù)的綜合資源，是一種用于鑒別DEGs的非參數(shù)方法，不需要對(duì)數(shù)據(jù)進(jìn)行分布假設(shè)，就可以對(duì)原始計(jì)數(shù)或之前歸一化或轉(zhuǎn)換的數(shù)據(jù)集進(jìn)行差異表達(dá)分析。由于噪聲分布來(lái)自實(shí)際數(shù)據(jù)，NOISeq方法能更好地適應(yīng)數(shù)據(jù)集的大小，可更有效地控制錯(cuò)誤發(fā)現(xiàn)率[6]。

背最長(zhǎng)肌是研究豬肉品質(zhì)最常用的材料。本研究選取了6頭松遼黑豬與6頭長(zhǎng)白豬，利用其背最長(zhǎng)肌組織的轉(zhuǎn)錄本進(jìn)行RNA-seq分析，通過(guò)非參數(shù)檢驗(yàn)法NOISeq鑒定差異表達(dá)基因，并對(duì)差異表達(dá)基因的功能進(jìn)行注釋分析，以篩選出影響豬脂肪沉積的關(guān)鍵基因。

1 材料和方法

1.1 試驗(yàn)動(dòng)物

本研究選用來(lái)自天津市寧河原種豬場(chǎng)的6頭松遼黑豬與6頭長(zhǎng)白豬，所有動(dòng)物在同一條件下飼養(yǎng)至體質(zhì)量為100 kg左右時(shí)進(jìn)行屠宰，屠宰在華都陽(yáng)光食品公司進(jìn)行，屠宰標(biāo)準(zhǔn)按照GB/T 17236—2008《生豬屠宰操作規(guī)程》進(jìn)行[7]。屠宰后，采集背最長(zhǎng)肌組織樣品，并移至液氮中進(jìn)行冷凍保存，待提取RNA使用。

1.2 總RNA的提取、建庫(kù)及測(cè)序

采用TRIzol法，按照產(chǎn)品說(shuō)明提取每頭豬背最長(zhǎng)肌組織的總RNA，共提取12個(gè)樣本，然后用1%瓊脂糖凝膠電泳初步檢測(cè)總RNA是否有降解及雜質(zhì)，并用紫外分光光度計(jì)測(cè)定波長(zhǎng)，評(píng)價(jià)所得RNA質(zhì)量。cDNA文庫(kù)制備與測(cè)序在廣州基迪奧生物科技有限公司進(jìn)行。參考Illumina Tru Seq TM RNA樣品制備試劑盒操作說(shuō)明進(jìn)行構(gòu)建，最后用構(gòu)建好的文庫(kù)通過(guò)Illumina HiSeq 2500平臺(tái)進(jìn)行雙末端(Paire-end，PE)測(cè)序。

1.3 測(cè)序數(shù)據(jù)的處理

利用Trimomatic[8]對(duì)原始數(shù)據(jù)(raw reads)進(jìn)行過(guò)濾，并利用FastQC軟件[9](http：//www.bioinformatics.babraham.ac.uk/projects/fastqc/)對(duì)raw reads進(jìn)行質(zhì)控并統(tǒng)計(jì)。質(zhì)控標(biāo)準(zhǔn)為：去除含 adapter 的 reads；去除含N比例大于10%的reads；去除低質(zhì)量reads(質(zhì)量值Q≤20的堿基數(shù)占整條read的50%以上)[10]。質(zhì)控后獲得高質(zhì)量的序列數(shù)據(jù)(Clean reads)，用TopHat2[11]將其比對(duì)到豬的參考基因組Susscrofa 11.1中。Cufflinks[12]用于數(shù)據(jù)的轉(zhuǎn)錄本組裝、注釋與合并。用HT-seq軟件[13]計(jì)算每個(gè)樣本中基因的表達(dá)量。

1.4 差異表達(dá)基因的篩選

用NOISeq R數(shù)據(jù)包篩選差異表達(dá)基因，以log2|fold change|>1、P≥0.8為標(biāo)準(zhǔn)進(jìn)行篩選。NOISeq統(tǒng)計(jì)方法是基于無(wú)參估計(jì)鑒定差異表達(dá)基因的，所以并沒有常見的P值，該方法的優(yōu)點(diǎn)在于能夠根據(jù)測(cè)序數(shù)據(jù)的大小對(duì)假陽(yáng)性率進(jìn)行有效控制，所采用的閾值P=0.8意味著該基因有差異表達(dá)的可能性是無(wú)差異表達(dá)的4倍，與其他方法中(cufflink、DESeq、DEGSeq等)所用的P=0.001篩選標(biāo)準(zhǔn)相持平[14]。

1.5 基因的生物學(xué)功能分析

GO(http：//www.geneontology.org)為基因數(shù)據(jù)提供功能分類，包括生物過(guò)程(Biological process，BP)、分子功能(Molecular function，MF)和細(xì)胞成分(Cellular component，CC)3個(gè)方面。GO分析是一種廣泛使用的基因和基因產(chǎn)物注釋工具。通過(guò)GO term注釋的差異基因，計(jì)算每個(gè)term的基因列表和基因數(shù)目，應(yīng)用超幾何檢驗(yàn)，找出與整個(gè)基因組背景相比差異表達(dá)基因中顯著富集的GO條目，從而找出差異表達(dá)基因顯著相關(guān)的生物學(xué)功能。KEGG(http：//www.genome.ad.jp/kegg/)是一個(gè)網(wǎng)絡(luò)網(wǎng)站，可以分析、解釋和可視化基因功能，是通路分析的主要數(shù)據(jù)庫(kù)[15]。本研究利用基迪奧生物信息云平臺(tái)(https：//www.omicshare.com/)對(duì)差異表達(dá)基因進(jìn)行GO功能富集分析和KEGG通路富集分析，P<0.05時(shí)被鑒定為顯著富集通路和GO條目。

2 結(jié)果與分析

2.1 測(cè)序數(shù)據(jù)的質(zhì)控及比對(duì)信息

對(duì)每個(gè)樣本的測(cè)序數(shù)據(jù)質(zhì)控后獲得約1.0×107個(gè)Clean reads，與豬的參考基因組序列(Sus scrofa 11.1)比對(duì)，平均比對(duì)率約為88.8%。通過(guò)對(duì)轉(zhuǎn)錄本進(jìn)行組裝和拼接，最終得到23 471個(gè)轉(zhuǎn)錄本(FPKM≥0.1)。由測(cè)序樣本的reads數(shù)和比對(duì)率可知，樣本質(zhì)量和測(cè)序測(cè)量均達(dá)到了預(yù)期(表1)。

表1 樣本比對(duì)結(jié)果統(tǒng)計(jì)Tab.1 Statistics of sample comparison results

2.2 基因表達(dá)整體分析

為了了解reads在不同樣本中的分布，比較了試驗(yàn)中所有樣本中的reads數(shù)，發(fā)現(xiàn)基因的表達(dá)水平基本相當(dāng)(圖1-A)，只有不到35%的基因在每一個(gè)樣本中都有1個(gè)以上的CPM(圖1-B)。

A為動(dòng)態(tài)表達(dá)范圍圖，該圖顯示了每個(gè)樣本中CPM值超過(guò)0的基因的分布；B圖反映了低表達(dá)基因在總轉(zhuǎn)錄本中所占的比例。A is the dynamic expression range map，which shows the distribution of differentially expressed genes with CPM values over 0 in each sample；B map reflects the proportion of low-expressed genes in total transcripts.

2.3 生物型檢測(cè)

分別選取松遼黑豬和長(zhǎng)白豬中的1個(gè)樣本進(jìn)行生物型檢測(cè)，發(fā)現(xiàn)檢測(cè)到的大多數(shù)轉(zhuǎn)錄本為蛋白質(zhì)編碼RNA，有少部分RNA也能被檢測(cè)到，如miRNA、lincRNA和scaRNA等，檢測(cè)結(jié)果符合預(yù)期(圖2)。

圖2 生物型分布Fig.2 Biotype distribution(per sample)

2.4 差異表達(dá)基因篩選

在松遼黑豬和長(zhǎng)白豬的背最長(zhǎng)肌組織樣本中共檢測(cè)到23 471個(gè)基因表達(dá)(FPKM≥0.1)。2個(gè)不同品種豬相比，共篩選到664個(gè)差異表達(dá)基因(圖3)，其中，364個(gè)基因在松遼黑豬中高表達(dá)，300個(gè)基因在長(zhǎng)白豬中高表達(dá)。

A圖代表2個(gè)品種豬的表達(dá)值總圖，其中，紅點(diǎn)代表差異表達(dá)基因；橫坐標(biāo)代表基因在長(zhǎng)白豬中的表達(dá)水平，縱坐標(biāo)代表基因在松遼黑豬中的表達(dá)水平。B圖代表(M，D)值和差異表達(dá)基因匯總圖。橫坐標(biāo)表示值，而縱坐標(biāo)表示值，黑點(diǎn)表示(M，D)值，其中(|M|>1)，紅點(diǎn)表示差異表達(dá)基因。

2.5 差異表達(dá)基因的功能富集分析

為了更好地了解所篩選到的DEGs的功能及其與表型的關(guān)系，本研究對(duì)篩選到的DEGs進(jìn)行了功能富集分析，包括GO功能分析和KEGG信號(hào)通路分析。GO功能富集分析發(fā)現(xiàn)，差異表達(dá)基因富集于52條GO條目，且生物過(guò)程富集到的GO terms最多，由此可知，生物過(guò)程是松遼黑豬和長(zhǎng)白豬在背最長(zhǎng)肌中差異表達(dá)基因行使的主要生物學(xué)功能(圖4)。對(duì)DEGs進(jìn)行KEGG信號(hào)通路富集分析，發(fā)現(xiàn)這些基因共參與了275條KEGG信號(hào)通路，其中，有52條顯著富集通路，其通路主要為脂肪酸代謝(Fatty acid metabolism)、胰島素抵抗(Insulin resistance)、PPAR信號(hào)通路(PPAR signaling pathway)、AMPK信號(hào)通路(AMPK signaling pathway)、脂肪細(xì)胞因子信號(hào)通路(Adipocytokine signaling pathway)和胰島素信號(hào)通路(Insulin signaling pathway)等(圖5)，根據(jù)顯著性通路篩選到LPIN1、FADS1、FADS2、PLIN2、PPARGC1A、PRKAG2和ACSL1等與脂類代謝和肌肉發(fā)育相關(guān)的基因(表2)。

1.細(xì)胞過(guò)程； 2.生物調(diào)節(jié)； 3.生物過(guò)程的調(diào)節(jié)； 4.對(duì)刺激的反應(yīng)； 5.代謝過(guò)程； 6.多細(xì)胞有機(jī)體過(guò)程； 7.生物過(guò)程的正向調(diào)控； 8.信令； 9.發(fā)展過(guò)程； 10.定位； 11.細(xì)胞成分組織或生物發(fā)生； 12.生物過(guò)程負(fù)調(diào)控； 13.免疫系統(tǒng)過(guò)程； 14.運(yùn)動(dòng)； 15.多組織進(jìn)程； 16.細(xì)胞增殖； 17.生物黏附； 18.繁殖； 19.生殖過(guò)程； 20.行為； 21.成長(zhǎng)； 22.細(xì)胞殺傷； 23.有節(jié)奏的過(guò)程； 24.解毒； 25.色素沉著； 26.細(xì)胞聚集； 27.細(xì)胞； 28.細(xì)胞部分； 29.細(xì)胞器； 30.膜； 31.膜部分； 32.細(xì)胞器部分； 33.細(xì)胞外區(qū)； 34.含蛋白質(zhì)的復(fù)合物； 35.胞外區(qū)部分； 36.膜封閉內(nèi)腔； 37.突觸； 38.超分子配合物； 39.細(xì)胞連接； 40.突觸部分； 41.黏合物； 42.催化活性； 43.分子功能調(diào)節(jié)劑； 44.分子換能器活性； 45.運(yùn)輸活動(dòng)； 46.結(jié)構(gòu)分子活性； 47.轉(zhuǎn)錄調(diào)節(jié)活性； 48.抗氧化活性； 49.自噬受體活性； 50.蛋白質(zhì)標(biāo)簽； 51.分子載體活性； 52.翻譯調(diào)節(jié)活動(dòng)。

3 結(jié)論與討論

本研究選用在表型方面存在較大差異的松遼黑豬和長(zhǎng)白豬為研究對(duì)象，對(duì)其背最長(zhǎng)肌組織進(jìn)行高通量轉(zhuǎn)錄組測(cè)序，篩選與脂肪沉積相關(guān)的差異表達(dá)基因。通過(guò)非參數(shù)檢驗(yàn)法NOISeq共篩選出664個(gè)DEGs，其中，364個(gè)基因在松遼黑豬中高表達(dá)，300個(gè)基因在長(zhǎng)白豬中高表達(dá)。對(duì)其進(jìn)行生物學(xué)功能分析，鑒定出脂肪細(xì)胞因子信號(hào)通路、PPAR信號(hào)通路、AMPK信號(hào)通路、胰島素信號(hào)通路和脂肪酸代謝等與脂肪沉積相關(guān)的通路，根據(jù)顯著性通路篩選出了LPIN1、FADS1、FADS2、PLIN2、PPARGC1A、PRKAG2和ACSL1等與脂類代謝和肌肉發(fā)育相關(guān)的基因。該研究結(jié)果為深入研究松遼黑豬和長(zhǎng)白豬之間脂肪沉積和肌肉發(fā)育性狀差異的分子機(jī)制提供了堅(jiān)實(shí)的基礎(chǔ)。豬脂肪組織的沉積能力具有品種差異，且不同品種的豬脂肪沉積能力具有較大差異[16]。本研究找到了與松遼黑豬和長(zhǎng)白豬脂肪沉積相關(guān)的重要功能基因和信號(hào)通路，如mTOR信號(hào)通路(mTOR signaling pathway)、PPAR信號(hào)通路和脂肪酸代謝，與之相關(guān)的基因包括脂質(zhì)1基因(Lipin 1，LPIN1)、脂滴蛋白2基因(Perilipin-2，PLIN2)和脂肪酸去飽和酶1基

圖5 KEGG通路富集分析Fig.5 KEGG pathway enrichment analysis

表2 脂質(zhì)代謝相關(guān)KEGG通路Tab.2 Lipid metabolism-related KEGG pathway

因(Fatty acid desaturase 1，F(xiàn)ADS1)。LPIN1基因的表達(dá)是脂肪細(xì)胞分化所必需的，其作為核轉(zhuǎn)錄共激活因子與一些過(guò)氧化物酶體增殖物激活受體一起調(diào)節(jié)參與脂質(zhì)代謝其他基因的表達(dá)。He等[17]研究發(fā)現(xiàn)，在豬中LPIN1基因與肌內(nèi)脂肪沉積及肌肉品質(zhì)存在顯著相關(guān)性，可以影響肌內(nèi)脂肪含量(IMF)、瘦肉率和風(fēng)味，所以將LPIN1基因作為改善肌肉品質(zhì)的候選基因之一。PLIN2基因位于脂滴表面，稱為脂滴包被蛋白，隸屬于脂滴包被蛋白PAT家族蛋白，其主要功能是參與脂質(zhì)代謝[18]。有研究表明，PLIN2基因在肌間脂肪較高的豬種中表達(dá)水平較高[19]。FADS1基因是脂肪酸去飽和酶(FADS)基因家族的成員，去飽和酶通過(guò)在脂肪?；湹亩x碳之間引入雙鍵來(lái)調(diào)節(jié)脂肪酸的不飽和度。宋文莉等[20]研究發(fā)現(xiàn)，F(xiàn)ADS1基因?qū)Ζ?3不飽和長(zhǎng)鏈脂肪酸的代謝有重要作用，加之FADS1基因在本研究中參加了脂肪酸代謝，因此，推測(cè)FADS1基因?yàn)榕c豬脂肪沉積相關(guān)的候選基因，目前關(guān)于FADS1基因在豬上的研究還很少。另外，發(fā)現(xiàn)FADS2、過(guò)氧化物酶體增殖物激活受體γ輔助激活因子1α基因(Peroxisome proliferators activated receptor γ coactivator1α，PPARGC1A)、蛋白激酶腺嘌呤核糖核苷酸激活的非催化亞基γ2基因(Protein kinase AMP-activated non-catalytic subunit gamma 2，PRKAG2)和長(zhǎng)鏈脂酰輔酶A合成酶1基因(Long chain acyl-CoA synthetase 1，ACSL1)均同時(shí)參與了不同信號(hào)通路和代謝途徑。FADS2屬于脂肪酸去飽和酶(FADS)基因家族的成員，主要作用為通過(guò)在脂肪?；湹奶囟ㄌ荚又g引入雙鍵來(lái)調(diào)節(jié)脂肪酸的不飽和度。本研究發(fā)現(xiàn)，F(xiàn)ADS2基因在PPAR信號(hào)通路和脂肪酸代謝通路中均有參與，前人研究發(fā)現(xiàn)，F(xiàn)ADS2基因與多不飽和脂肪酸(PUFA)之間存在關(guān)聯(lián)[21]。另有研究發(fā)現(xiàn)，F(xiàn)ADS2基因SNP突變影響人體脂肪組織及血漿中PUFA的含量，從而影響人類健康，人類全基因組關(guān)聯(lián)分析證實(shí)，F(xiàn)ADS2基因?qū)χ惔x疾病有重要作用[22]；因此，推測(cè)FADS2基因可能在豬脂肪沉積中也發(fā)揮著重要作用。前人研究發(fā)現(xiàn)，PPARGC1A的表達(dá)有利于脂肪細(xì)胞的分化成熟，與脂質(zhì)代謝相關(guān)[23]，且PPARGC1A基因?qū)ι矬w內(nèi)棕色脂肪的生成具有重大影響[24]。Gandolfi等[25]研究發(fā)現(xiàn)，PPARGC1A基因與豬肉質(zhì)性狀存在顯著關(guān)聯(lián)，且由于物種間存在差異，可作為候選基因用于豬肉質(zhì)性狀的改良。本研究中，PPARGC1A基因參與了脂肪細(xì)胞因子信號(hào)通路、AMPK信號(hào)通路和胰島素抵抗等多條通路，所以，可將PPARGC1A基因作為豬背脂沉積的重要候選基因之一。PRKAG2參與了AMPK信號(hào)通路、脂肪細(xì)胞因子信號(hào)通路和胰島素信號(hào)通路等多條通路，PRKAG2基因是腺嘌呤核糖核苷酸活化蛋白激酶(AMP-activated protein kinase，AMPK)家族重要的成員，其突變與Wolff-Parkinson-White綜合征、家族性肥厚性心肌病和心臟糖原貯積病有關(guān)。其在肌肉脂肪代謝、能量代謝等方面發(fā)揮著重要作用[26]。Jing等[27]通過(guò)對(duì)高、低剩余采食量的豬骨骼肌進(jìn)行高通量測(cè)序發(fā)現(xiàn)，PRKAG2在骨骼肌能量代謝過(guò)程中起著至關(guān)重要的作用。由此推測(cè)，PRKAG2基因的表達(dá)與豬的脂肪沉積存在顯著相關(guān)性，可作為影響脂肪沉積的重要候選基因。ACSL1是長(zhǎng)鏈酯酰輔酶A合成酶家族(Long-chain fatty acyl-CoA synthetases，ACSLs)的成員之一，該家族由ACSL1、ACSL3、ACSL4、ACSL5和ACSL6組成[28]，而ACSL1是ACSLs家族中最早被發(fā)現(xiàn)的亞型，ACSL1在脂質(zhì)生物合成和脂肪酸降解中發(fā)揮關(guān)鍵作用。有研究表明，ACSL1在脂肪酸的活化、轉(zhuǎn)運(yùn)和降解以及脂類生成等過(guò)程中發(fā)揮著至關(guān)重要的作用[29]。本研究中，ACSL1基因參與了PPAR信號(hào)通路、脂肪酸代謝和脂肪細(xì)胞因子信號(hào)通路，數(shù)量性狀基因座(Quantitative trait loci，QTL)分析認(rèn)為，ACSL1作為最重要位置和功能候選基因，影響豬肉中脂肪酸組成[30]；因此，將ACSL1基因作為脂肪沉積的候選基因。

綜上所述，通過(guò)對(duì)松遼黑豬和長(zhǎng)白豬背最長(zhǎng)肌之間差異表達(dá)基因進(jìn)行鑒定與分析，篩選出多個(gè)參與肌肉發(fā)育和脂類代謝的重要候選基因，包括LPIN1、FADS1、FADS2、PLIN2、PPARGC1A、PRKAG2和ACSL1等，這些基因目前已知的功能均與脂肪沉積相關(guān)，但在豬中影響組織脂類代謝的機(jī)制還有待深入研究，未來(lái)對(duì)這些重要功能基因的深入驗(yàn)證，可能會(huì)揭示影響豬生長(zhǎng)發(fā)育和豬肉品質(zhì)的重要基因，為未來(lái)滿足社會(huì)需求的新品種遺傳育種提供理論基礎(chǔ)。