影響藏豬與大約克夏豬肌肉品質(zhì)的lncRNA的篩選及功能分析

宋明坤，薛明明，張力戈，顏 鐸，夏 寧，商 鵬，李新建，韓雪蕾，喬瑞敏，李秀領(lǐng)，李 明，王克君*

(1. 河南農(nóng)業(yè)大學(xué)動(dòng)物科技學(xué)院，鄭州 450002； 2. 西藏農(nóng)牧學(xué)院動(dòng)物科學(xué)學(xué)院，林芝 860000)

在生豬產(chǎn)業(yè)中，豬的肉質(zhì)性狀是重要的經(jīng)濟(jì)性狀，主要包括肉色、滴水損失、大理石花紋、肌內(nèi)脂肪含量等指標(biāo)。豬肉品質(zhì)與肌肉生長發(fā)育和脂肪沉積密切相關(guān)，在很大程度上直接影響消費(fèi)者對(duì)豬肉的選擇[1-2]。藏豬是我國寶貴的地方種質(zhì)資源，主要起源于西藏高原。由于長期在惡劣的環(huán)境中經(jīng)歷自然選擇[3-5]，藏豬形成了耐粗飼、耐低氧、抗病力較強(qiáng)、脂肪沉積能力強(qiáng)等品種特性[6-7]。藏豬肉因其肉質(zhì)鮮美，富含蛋白質(zhì)和氨基酸等特點(diǎn)深受廣大消費(fèi)者的喜愛[8]。大約克夏豬在世界范圍內(nèi)廣泛分布，是世界著名的瘦肉豬品種，具有生長速度快、瘦肉率高等優(yōu)良特點(diǎn)。已有研究表明，藏豬和大約克夏豬在肌肉品質(zhì)上存在較大差異[9]，因此兩者是比較豬肌肉品質(zhì)的理想試驗(yàn)?zāi)Ｐ汀?/p>

長鏈非編碼RNA(long non-coding RNA，lncRNA)是一類長度大于200 bp，幾乎不具有編碼蛋白能力的非編碼RNA。大量研究表明，lncRNAs在細(xì)胞增殖、分化、發(fā)育、疾病等生命活動(dòng)中發(fā)揮著重要的調(diào)節(jié)功能[10-11]。隨著高通量測序技術(shù)的不斷發(fā)展與應(yīng)用，發(fā)現(xiàn)lncRNAs在家養(yǎng)動(dòng)物肌肉生長發(fā)育和脂肪沉積過程中發(fā)揮著重要作用[12-13]。目前已有一些關(guān)于不同豬種背最長肌轉(zhuǎn)錄組的比較分析。Wang等[14]通過對(duì)大約克夏豬和淮南豬胚胎期背最長肌的組織進(jìn)行比較分析，發(fā)現(xiàn)大量差異表達(dá)的lncRNAs可能在肌肉生長發(fā)育和脂肪沉積過程中發(fā)揮著重要作用，并對(duì)候選差異表達(dá)lncRNA進(jìn)行功能驗(yàn)證，發(fā)現(xiàn)IMFlnc1可抑制脂肪的生成。Li等[15]通過分析圩豬和大約克夏豬背最長肌的轉(zhuǎn)錄組數(shù)據(jù)，發(fā)現(xiàn)許多差異表達(dá)lncRNAs可能通過調(diào)節(jié)其潛在靶基因從而影響肌內(nèi)脂肪的發(fā)育過程。然而，關(guān)于藏豬和大約克夏豬背最長肌組織轉(zhuǎn)錄組lnc-RNAs的比較分析還鮮有報(bào)道。

本研究基于藏豬和大約克夏豬的背最長肌組織轉(zhuǎn)錄組測序，對(duì)不同豬種中的lncRNAs和基因進(jìn)行生物信息學(xué)分析，以確定影響藏豬和大約克夏豬肉品質(zhì)的候選lncRNAs和基因。這些研究結(jié)果將有助于從肌肉生長發(fā)育和脂肪沉積的角度解析藏豬和大約克夏豬肉質(zhì)性狀形成的分子機(jī)制，為豬的遺傳改良提供理論參考。

1 材料與方法

1.1 肌肉品質(zhì)表型和RNA-seq數(shù)據(jù)的收集

本研究試驗(yàn)動(dòng)物為180日齡健康藏公豬(n=3)和大約克夏公豬(n=3)，均飼養(yǎng)于相同飼養(yǎng)環(huán)境，且保證飼養(yǎng)條件一致[16]。背最長肌肉品質(zhì)測定數(shù)據(jù)來源于四川省畜牧科學(xué)院[16]。肉質(zhì)指標(biāo)主要包括滴水損失、pH45 min、pH24 h、肉色、大理石花紋、肌內(nèi)脂肪含量、眼肌面積、肌纖維橫截面積和維生素B1(VB1)含量。RNA-seq數(shù)據(jù)來源于同一樣本藏豬和大約克夏豬的背最長肌組織，均收集于NCBI數(shù)據(jù)庫(登錄號(hào)：SRP090525)[16]。

1.2 lncRNAs的鑒定流程

測序數(shù)據(jù)包含帶接頭和低質(zhì)量的reads，故使用FastQC軟件對(duì)測序數(shù)據(jù)進(jìn)行質(zhì)量控制。過濾后的高質(zhì)量序列(clean reads)用于后續(xù)分析。使用HISAT2(v2.0.0-beta)軟件將FASTQ格式的clean reads比對(duì)到豬的參考基因組Sus scrofa 11.1。

使用StringTie(v1.3.0)軟件將比對(duì)后的數(shù)據(jù)用于轉(zhuǎn)錄本的組裝。利用gffcompare(v0.10.1)軟件將組裝好的轉(zhuǎn)錄本與Ensemble基因集(Sscrofa 11.1.91)相比較，將具有“i”，“u”，“o”和“x”這類代碼的轉(zhuǎn)錄本進(jìn)行后續(xù)lncRNA的鑒定。為了避免不完全組裝和假陽性率，對(duì)轉(zhuǎn)錄本以長度>200 bp和外顯子數(shù)>2的標(biāo)準(zhǔn)進(jìn)行過濾。使用CNCI[17]和CPC[18]軟件對(duì)過濾后的轉(zhuǎn)錄本進(jìn)行蛋白編碼潛力的預(yù)測，CNCI和CPC得分均低于0的轉(zhuǎn)錄本被認(rèn)為是lncRNAs。

1.3 lncRNAs和基因的差異表達(dá)分析

利用HTseq軟件[19]計(jì)算比對(duì)到每個(gè)基因的reads數(shù)。使用DEseq2軟件[20]鑒定差異表達(dá)lncRNAs和基因，篩選標(biāo)準(zhǔn)為差異表達(dá)倍數(shù)|log2(fold change)|≥1且顯著性P-adj<0.05。

1.4 lncRNAs的靶DEGs預(yù)測及功能富集分析

lncRNAs的順式調(diào)控作用被定義為對(duì)鄰近基因的影響[21]。Bedtools軟件[22]用于鑒定位于差異表達(dá)lncRNAs上游或下游100 kb的順式調(diào)控靶DEGs。利用皮爾遜相關(guān)性檢驗(yàn)計(jì)算lncRNAs和DEGs之間的表達(dá)相關(guān)性來預(yù)測差異表達(dá)lncRNAs反式調(diào)控的靶DEGs。

基于Metascape數(shù)據(jù)庫[23]對(duì)差異表達(dá)lncRNAs的靶DEGs進(jìn)行GO生物學(xué)過程(biological process)和KEGG信號(hào)通路功能富集分析。Cytoscape(v3.8.2)軟件[24]用于構(gòu)建與肌肉生長發(fā)育和脂肪沉積相關(guān)的候選差異表達(dá)lncRNAs、靶DEGs和生物學(xué)過程三者間的調(diào)控網(wǎng)絡(luò)。

1.5 lncRNAs的靶DEGs與肉質(zhì)性狀的關(guān)聯(lián)分析

基于R語言環(huán)境下使用“Corrplot”和“Psych”軟件包對(duì)差異表達(dá)lncRNAs的靶DEGs的表達(dá)水平和肉質(zhì)性狀的表型數(shù)據(jù)進(jìn)行關(guān)聯(lián)分析。皮爾遜相關(guān)系數(shù)用于度量基因與表型間的相關(guān)程度。FDR校正法用于確定每個(gè)相關(guān)事件的顯著性水平。

1.6 數(shù)據(jù)統(tǒng)計(jì)分析

使用EXCEL 2019對(duì)試驗(yàn)表型數(shù)據(jù)進(jìn)行初步整理，利用SPSS 22.0統(tǒng)計(jì)軟件進(jìn)行獨(dú)立樣本t檢驗(yàn)。結(jié)果使用“平均值±標(biāo)準(zhǔn)差”表示，P<0.05為差異顯著判斷標(biāo)準(zhǔn)，P<0.01為差異極顯著判斷標(biāo)準(zhǔn)。

2 結(jié) 果

2.1 藏豬和大約克夏豬肌肉品質(zhì)的比較

180日齡藏豬和大約克夏豬的背最長肌肌肉品質(zhì)測定結(jié)果顯示，大約克夏豬組的肌纖維橫截面積和眼肌面積均極顯著高于藏豬組(P<0.01)；藏豬組的肌肉pH45 min和肉色均顯著高于大約克夏豬組(P<0.05)，且肌內(nèi)脂肪含量、VB1含量極顯著高于大約克夏豬組(P<0.01)；而肌肉滴水損失、pH24 h和大理石花紋在兩組間不存在顯著性差異(P>0.05，圖1)。

pH45 min和pH24 h分別代表屠宰后45 min和24 h肌肉的pH。下圖同pH45 min and pH24 h represent the pH of the muscle at 45 min and 24 h after slaughtering, respectively. The same as below圖1 藏豬和大約克夏豬肌肉品質(zhì)的比較分析Fig.1 Comparative analysis of meat quality between Tibetan and Yorkshire pigs

2.2 轉(zhuǎn)錄組數(shù)據(jù)評(píng)估和lncRNAs的鑒定

本研究測序數(shù)據(jù)質(zhì)量控制后共獲得227 087 843個(gè)clean reads，包括藏豬組的114 520 669個(gè)和大約克夏豬組的112 567 174個(gè)。clean reads比對(duì)到豬參考基因組Sscrofa 11.1的比對(duì)率超過82%，其中超過85.5%的reads比對(duì)到CDS、5′UTR和3′UTR區(qū)域(圖2A)。Reads進(jìn)行轉(zhuǎn)錄本組裝后共獲得21 765個(gè) 轉(zhuǎn)錄本。對(duì)這些轉(zhuǎn)錄本進(jìn)行鑒定分析，共篩選出20 219個(gè)mRNAs和1 546個(gè)lncRNAs。藏豬組和大約克夏豬組共有的mRNAs有15 329個(gè)，在藏豬組和大約克夏豬組中特異性表達(dá)的mRNAs分別有2 839和2 051個(gè)(圖2B)；藏豬組和大約克夏豬組中共有的lncRNAs有1 314個(gè)，在藏豬組和大約克夏豬組中特異性表達(dá)的lncRNAs分別有145和87個(gè) (圖2C)。對(duì)lncRNAs進(jìn)行分類，如圖2D所示，1 546個(gè)lncRNAs中基因間lncRNAs(lincRNAs)有654個(gè)，反義lncRNAs(lncNATs)有556個(gè)，內(nèi)含子lncRNAs(ilncRNAs)有336個(gè)。

2.3 差異表達(dá)lncRNAs和DEGs的篩選

根據(jù)|log2(fold change)|≥1且P-adj<0.05的標(biāo)準(zhǔn)篩選差異表達(dá)基因，結(jié)果表明在兩組中共篩選出49個(gè)差異表達(dá)lncRNAs和154個(gè)DEGs。與藏豬組相比，大約克夏豬組中有28個(gè)上調(diào)lncRNAs和21個(gè)下調(diào)lncRNAs(圖3A)，74個(gè)上調(diào)DEGs和80個(gè) 下調(diào)DEGs(圖3B)。

2.4 差異表達(dá)lncRNAs靶基因的預(yù)測

對(duì)差異表達(dá)lncRNAs進(jìn)行順式作用靶基因預(yù)測，共獲得了8對(duì)具有潛在順式調(diào)控作用的差異表達(dá)lncRNAs和靶DEGs(表1)。其中MSTRG.42387.1上下游100 kb共存在3個(gè)潛在順式作用靶DEGs，包括 ATP合酶亞基8(ATP synthase F0 subunit 8，ATP8)、NADH脫氫酶亞單位1(NADH dehydrogenase subunit 1，ND1)和NADH脫氫酶亞單位2(NADH dehydrogenase subunit 1，ND2)。MSTRG.27311.5、MSTRG.42191.1、MSTRG.43545.1分別靶向ENSSSCG00000031583、β-防御素-1(beta defensin 1，DEFB1)和diaphanous相關(guān)成蛋白2(diaphanous related formin 2，DIAPH2)基因。MSTRG.8729.4和MSTRG.8729.8同時(shí)靶向肌球蛋白重鏈13(myosin heavy chain 13，MYH13)。

根據(jù)皮爾遜相關(guān)性系數(shù)|correlation|≥0.95且P<0.05 的篩選標(biāo)準(zhǔn)發(fā)現(xiàn)48個(gè)差異表達(dá)lncRNAs與138個(gè)DEGs存在493個(gè)反式作用關(guān)系對(duì)。在大約克夏豬組和藏豬組中，差異表達(dá)倍數(shù)前5的lncRNAs 及其反式作用的靶DEGs如表2所示。值得注意的是，不同的lncRNAs其潛在的反式作用靶DEGs的數(shù)量差別很大，MSTRG.42621.1反式作用調(diào)控25個(gè)DEGs，MSTRG.151.1反式作用調(diào)控18個(gè)DEGs，而MSTRG.8729.3和MSTRG.34836.20分別只有2個(gè)反式作用靶DEGs。

A. 6個(gè)樣本的reads分布統(tǒng)計(jì)圖:TA代表藏豬，YA代表大約克夏豬；B. 藏豬與大約克夏豬中鑒定出的mRNAs數(shù)量；C. 藏豬與大約克夏豬鑒定出的lncRNAs數(shù)量；D. lncRNAs的分類統(tǒng)計(jì)圖:lincRNAs代表基因間lncRNAs，lncNATs代表反義lncRNAs，ilncRNAs代表內(nèi)含子lncRNAsA. Reads distribution of 6 samples： TA is Tibetan pigs, and YA is Yorkshire pigs; B. The number of mRNAs identified in Tibetan and Yorkshire pigs; C. The number of lncRNAs identified in Tibetan and Yorkshire pigs; D. The classification and statistics of lncRNAs: lincRNAs represent intergenic lncRNAs, lncNATs represent antisense lncRNAs, ilncRNAs represent introns lncRNAs圖2 lncRNAs和mRNAs的鑒定Fig.2 Identification of lncRNAs and mRNAs

紅色代表上調(diào)基因；灰色代表表達(dá)無顯著變化的基因；綠色代表下調(diào)基因Red represents up-regulated genes; gray represents genes with no significant change in expression; green represents down-regulated genes圖3 差異表達(dá)lncRNAs和DEGs的火山圖Fig.3 Volcano plots of differentially expressed lncRNAs and DEGs

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

為了揭示差異表達(dá)lncRNAs潛在的生物學(xué)功能，本研究對(duì)差異表達(dá)lncRNAs的反式作用DEGs進(jìn)行GO生物學(xué)過程和KEGG通路功能富集分析，結(jié)果顯示，靶DEGs顯著富集到了循環(huán)系統(tǒng)過程(circulatory system process)、骨骼系統(tǒng)發(fā)育(skeletal system development)、心臟發(fā)育(heart development)、磷脂酰肌醇3-激酶信號(hào)的正調(diào)控(positive regulation of phosphatidylinositol 3-kinase signaling)、脂肪細(xì)胞分化(fat cell differentiation)等生物學(xué)過程和一個(gè)糖尿病并發(fā)癥的AGE-RACE信號(hào)通路(AGE-RAGE signaling pathway in diabetic complications)(圖4)。為進(jìn)一步挖掘與豬肌肉生長發(fā)育和脂肪沉積密切相關(guān)的候選差異表達(dá)lncRNAs和靶DEGs，本研究構(gòu)建了差異表達(dá)lncRNAs及其反式作用DEGs共表達(dá)調(diào)控網(wǎng)絡(luò)，其中共有39個(gè)差異表達(dá)lncRNAs與23個(gè)差異表達(dá)基因存在較強(qiáng)的共表達(dá)關(guān)系且與肌肉生長發(fā)育和脂肪沉積相關(guān)(圖5)，這些靶DEGs顯著富集在磷脂酰肌醇3-激酶信號(hào)的正調(diào)控、肌肉結(jié)構(gòu)發(fā)育(muscle structure develop-ment)、脂肪細(xì)胞分化(fat cell differentiation)、JAK-STAT級(jí)聯(lián)的正調(diào)控(positive regulation of JAK-STAT cascade)、MAPK級(jí)聯(lián)的正調(diào)控(positive regulation of MAPK cascade)等生物學(xué)過程(表3)。lncRNAs-DEGs-生物學(xué)過程調(diào)控網(wǎng)絡(luò)顯示,MSTRG.34836.10、MSTRG.151.1調(diào)控的靶DEGs最多，其次是MSTRG.21406.1和MSTRG.26709.1。與lncRNAs存在靶向關(guān)系最多的基因?yàn)門GFB2，其次是JAK2，這些靶基因共同被MSTRG.12095.1、MSTRG.27311.5、MSTRG.33150.6、MSTRG.42387.11、MSTRG.42387.18和MSTRG.42387.2所調(diào)控，并共同參與到磷脂酰肌醇3-激酶信號(hào)的正調(diào)控、MAPK級(jí)聯(lián)的正調(diào)控和MAPK級(jí)聯(lián)的調(diào)控的生物學(xué)過程中。

2.6 靶DEGs與肉質(zhì)性狀的關(guān)聯(lián)分析

為了更好地解析差異表達(dá)lncRNAs反式作用靶DEGs的潛在功能，本研究將上述篩選的23個(gè)與肌肉生長發(fā)育和脂肪沉積生物學(xué)過程相關(guān)的靶DEGs與藏豬和大約克夏豬的肉質(zhì)性狀進(jìn)行關(guān)聯(lián)分析。結(jié)果表明，有11個(gè)差異表達(dá)靶基因與肉質(zhì)性狀顯著相關(guān)(圖6)。CCN1(cellular communication network factor 1)與滴水損失顯著相關(guān)(P<0.05)；肽酶抑制蛋白16(peptidase inhibitor 16，PI16)和原癌基因FOS(Fos proto-oncogene，F(xiàn)OS)與肌肉pH24 h顯著相關(guān)(P<0.05)；popeye結(jié)構(gòu)域2(popeye domain containing 2，POPDC2)、PI16、FOS和生物鐘基因PER2(period circadian regulator 2)共4個(gè)基因與肉色(CLR)顯著相關(guān)(P<0.05)；PI16、FOS、鈣調(diào)磷酸酶B同源蛋白2(RB1 inducible coiled-coil 1，RB1CC1)、腱蛋白(chordin，CHRND)和受體活性修飾蛋白3(receptor activity modifying protein 3，RAMP3)共5個(gè)基因與大理石花紋(MBL)顯著相關(guān)(P<0.05)；POPDC2、PER2和P-選擇素(selectin P，SELP)共3個(gè)基因與肌內(nèi)脂肪含量顯著相關(guān)(P<0.05)；轉(zhuǎn)化生長因子β2(transforming growth factor beta 2，TGFB2)和Janus激酶2基因(Janus kinase 2，JAK2)與眼肌面積(LMA)顯著相關(guān);POPDC2、PI16、FOS和PER2共4個(gè)基因與維生素B1含量(VB1)顯著相關(guān)(P<0.05)；關(guān)聯(lián)分析沒有發(fā)現(xiàn)任何基因與pH45 min和肌纖維橫截面積(CSAF)顯著相關(guān)(P>0.05)。有趣的是，PI16和FOS基因與pH24 h、CLR、MBL和VB1含量均顯著相關(guān)(P<0.05)；POPDC2和PER2基因與CLR、IMF和VB1含量均顯著相關(guān)(P<0.05)。

3 討 論

肉質(zhì)性狀屬于微效多基因控制的數(shù)量性狀，了解影響肉質(zhì)性狀的潛在分子機(jī)制對(duì)提高豬肉產(chǎn)量和肉質(zhì)改良具有重要意義。先前的研究表明，哺乳動(dòng)物基因組中含有大量的lncRNAs，一些lncRNAs已被發(fā)現(xiàn)在動(dòng)物生長發(fā)育等生命過程中發(fā)揮著非常重要的調(diào)節(jié)作用[12-13]。目前，lncRNAs的研究主要集中在人和小鼠等模式動(dòng)物[25-26]，在豬上的研究相對(duì)較少，對(duì)豬肌肉品質(zhì)相關(guān)lncRNAs的認(rèn)識(shí)更是不足。藏豬和大約克夏豬的肌肉品質(zhì)存在較大差異,是研究豬肌肉品質(zhì)理想的試驗(yàn)?zāi)Ｐ蚚27]。與藏豬組相比，大約克夏豬組具有較高的眼肌面積和較低的肌內(nèi)脂肪含量。本研究通過分析藏豬和大約克夏豬的背最長肌組織的轉(zhuǎn)錄組測序數(shù)據(jù)，篩選潛在影響肌肉生長發(fā)育和脂肪沉積相關(guān)的候選lncRNAs和候選基因，以探究lncRNAs在藏豬和大約克夏豬肌肉品質(zhì)調(diào)控中的作用。

左下圖所富集的聚類根據(jù)聚類的ID添加顏色。右上圖代表聚類的P值The enriched clusters in the left bottom figure were colored by the cluster ID. The color of clusters in the upper right figure represents the P-value of the clusters圖4 差異表達(dá)lncRNAs反式作用靶DEGs的功能富集分析圖Fig.4 Functional enrichment analysis of trans-acting target DEGs of differentially expressed lncRNAs

菱形代表差異表達(dá)lncRNAs，圓代表靶DEGs，正方形代表與肌肉生長和脂肪沉積相關(guān)的生物學(xué)過程，紅色代表基因表達(dá)上調(diào)，綠色代表基因表達(dá)下調(diào)The diamonds represent the differentially expressed lncRNAs, the circles represent the target DEGs, and the squares represent the biological processes associated with muscle development and fat deposition. The red edge represents up-regulated gene expression, and the green edge represents down-regulated gene expression圖5 差異表達(dá)lncRNAs和與肌肉生長和脂肪沉積相關(guān)的生物學(xué)過程相關(guān)靶DEGs的調(diào)控網(wǎng)絡(luò)Fig.5 Regulation network of differentially expressed lncRNAs and target DEGs enriched in biological processes related to muscle growth and fat deposition

3.1 差異表達(dá) lncRNA可能通過順式作用參與豬肌肉生長發(fā)育過程

lncRNAs一般不具有編碼蛋白的功能，可通過順式調(diào)控模式直接結(jié)合到臨近位置的DNA調(diào)節(jié)基因表達(dá)或以反式調(diào)控的模式調(diào)控lncRNA遠(yuǎn)端位置的基因表達(dá)[28]。Chen等[29]通過識(shí)別肌肉特異性表達(dá)lncRNAs臨近位置的蛋白質(zhì)編碼基因(<100 kb)以預(yù)測lncRNA的功能。本研究通過預(yù)測差異表達(dá)lncRNAs的順式調(diào)控靶基因，發(fā)現(xiàn)共有8對(duì)具有潛在順式調(diào)控作用的差異表達(dá)lncRNAs-靶DEGs對(duì)，其中l(wèi)ncRNA MSTRG.43545.1與順式調(diào)控靶基因DIAPH2存在較強(qiáng)的相關(guān)性。與藏豬相比，大約克夏豬的MSTRG.43545.1和DIAPH2的表達(dá)水平均顯著提高。有研究報(bào)道，DIAPH2在果蠅成肌細(xì)胞融合過程中起著至關(guān)重要的作用，可誘導(dǎo)成肌細(xì)胞融合成核并促進(jìn)肌動(dòng)蛋白絲的延長[30]。MSTRG.43545.1可能通過上調(diào)DIAPH2的表達(dá)促進(jìn)肌肉的生長發(fā)育過程。

3.2 差異表達(dá) lncRNAs可能通過反式作用調(diào)控豬肌肉生長和肌內(nèi)脂肪沉積過程

Li等[31]通過lncRNA-mRNA共表達(dá)關(guān)系分析了lncRNAs在骨骼肌纖維形成的反式調(diào)控作用。在本研究中，共確定了48個(gè)差異表達(dá)lncRNAs與138個(gè)DEGs存在493個(gè)反式作用關(guān)系對(duì)。為了明確這些差異表達(dá)lncRNAs的生物學(xué)功能，本研究對(duì)其反式作用的靶DEGs進(jìn)行富集分析，結(jié)果表明大量靶基因顯著富集到與肌肉生長發(fā)育和脂肪沉積相關(guān)的生物學(xué)過程。其中Notch信號(hào)通路已被證明在骨骼肌的發(fā)育和再生過程中起著關(guān)鍵作用，在骨骼肌衛(wèi)星細(xì)胞再生期間激活Notch信號(hào)通路可使肌肉干細(xì)胞保持靜息狀態(tài)[32]。此外，Notch信號(hào)還參與了調(diào)節(jié)白細(xì)胞和棕色脂肪細(xì)胞祖細(xì)胞的增殖和脂肪分化的過程[33]。磷脂酰肌醇3-激酶信號(hào)通路作為經(jīng)典的胰島素信號(hào)通路，在調(diào)節(jié)脂肪代謝中起著關(guān)鍵作用[34]，還可以通過激活A(yù)KT、mTOR等下游靶蛋白參與到骨骼肌生長發(fā)育的過程中[35]。MAPK級(jí)聯(lián)[36-37]和JAK-STAT級(jí)聯(lián)[38-39]也均已被證明在肌肉生長發(fā)育和脂肪沉積過程中發(fā)揮著重要作用。這些生物學(xué)過程可能是導(dǎo)致藏豬和大約克夏豬肌肉生長發(fā)育和脂肪沉積性狀差異顯著的重要因素。

為進(jìn)一步挖掘潛在影響肌肉生長發(fā)育和脂肪沉積的基因的互作關(guān)系，本研究構(gòu)建了lncRNAs-DEGs-生物學(xué)過程調(diào)控網(wǎng)絡(luò)，并將這些DEGs與藏豬和大約克夏豬的肉質(zhì)性狀進(jìn)行關(guān)聯(lián)分析。調(diào)控網(wǎng)絡(luò)顯示，MSTRG.12095.1、MSTRG.14538.2、MSTRG.25815.1、MSTRG.26232.1等共11個(gè)lncRNAs對(duì)TGFB2存在反式調(diào)控作用。TGFB2作為TGF-β/SMAD信號(hào)通路的成員之一，具有抑制骨骼肌衛(wèi)星細(xì)胞增殖和分化，促進(jìn)凋亡的作用[40]。本研究中，TGFB2 與眼肌面積顯著相關(guān)，且顯著富集到Notch信號(hào)通路、肌肉結(jié)構(gòu)發(fā)育等肌肉發(fā)育相關(guān)生物學(xué)過程。由此推測，在大量lncRNAs的調(diào)控作用下TGFB2可能 通過多種途徑共同參與到肌肉的生長發(fā)育過程中。在肌肉發(fā)育和再生過程中，HES1作為Notch信號(hào)通路下游的重要靶基因可通過抑制MYOD和MYOG的表達(dá)控制激活的肌肉干細(xì)胞增殖和分化之間的平衡[41]。與藏豬相比，MSTRG.151.1、MSTRG.21406.1、MSTRG.26709.1和HES1在大約克夏豬的表達(dá)顯著下調(diào)，推測MSTRG.151.1、MSTRG.21406.1和MSTRG.26709.1可能通過正向調(diào)控HES1的表達(dá)進(jìn)而參與肌肉的生長發(fā)育過程。

表3 差異表達(dá)靶基因顯著富集的與肌肉生長發(fā)育和脂肪沉積相關(guān)的GO生物學(xué)過程

差異表達(dá)lncRNAs調(diào)節(jié)的靶DEGs還通過參與脂肪沉積的相關(guān)生物學(xué)過程導(dǎo)致藏豬與大約克夏豬肌肉品質(zhì)差異。與藏豬相比，MSTRG.42191.1、MSTRG.11899.1和LEP在大約克夏豬中顯著下調(diào)，且MSTRG.42191.1是差異倍數(shù)第二大的lncRNAs。LEP基因多態(tài)性與豬[42]和綿羊[43]的脂肪特征相關(guān)。LEP基因是白色脂肪的標(biāo)志基因，調(diào)節(jié)脂肪沉積的重要因子，可促進(jìn)脂肪的生成和脂滴的形成[44]。LEP基因被顯著富集到了磷脂酰肌醇3-激酶信號(hào)的正調(diào)控、JAK-STAT級(jí)聯(lián)和MAPK級(jí)聯(lián)的調(diào)控和脂肪細(xì)胞分化等GO生物學(xué)過程中。據(jù)此推測，MSTRG.42191.1和MSTRG.11899.1可能通過靶向調(diào)控LEP基因的表達(dá)促進(jìn)藏豬的肌內(nèi)脂肪沉積，MSTRG.42191.1可能在肌內(nèi)脂肪沉積中發(fā)揮著更為重要的作用。PER2是一個(gè)節(jié)律基因，可通過直接調(diào)節(jié)PPARγ基因影響脂肪細(xì)胞的分化，PER2基因缺陷的小鼠表現(xiàn)出脂質(zhì)代謝的改變以及總?cè)８视秃头酋セ舅峒眲p少[45]。本研究中，PER2與肌內(nèi)脂肪含量顯著相關(guān)，并顯著富集到了脂肪細(xì)胞分化生物學(xué)過程。與藏豬相比，MSTRG.3334.1、MSTRG.8729.8在大約克夏豬中顯著上調(diào)，而PER2在大約克夏豬中顯著下調(diào)，推測MSTRG.3334.1、MSTRG.8729.8可能通過下調(diào)PER2的表達(dá)進(jìn)而阻礙脂肪細(xì)胞分化過程，從而導(dǎo)致了藏豬和大約克夏豬的肌內(nèi)脂肪含量的顯著差異。本研究篩選到許多與肌肉生長發(fā)育和脂肪沉積相關(guān)的新lncRNAs，為深入解析肌肉品質(zhì)的形成提供了重要基礎(chǔ)依據(jù)，其與靶基因的調(diào)控作用和具體的作用機(jī)制還有待進(jìn)一步探究。

圖中的肉質(zhì)性狀均使用原始數(shù)據(jù)除以對(duì)應(yīng)個(gè)體胴體重的方法進(jìn)行了歸一化校正。圓顏色和圓大小代表基因表達(dá)水平和肉質(zhì)性狀的相關(guān)系數(shù)，即圓顏色越深和圓越大代表相關(guān)系數(shù)越高，藍(lán)色代表正相關(guān)，紅色代表負(fù)相關(guān)。采用FDR校正法進(jìn)行顯著性檢驗(yàn)。*. P<0.05The meat quality traits in the figure are normalized and corrected by dividing the original data by the corresponding individual’s carcass weight. The color and size of the circles represent the correlation coefficients between gene expression levels and meat quality traits, that is, the darker the circle color and the larger the circle represent the higher the correlation coefficient, blue represents positive correlation, and red represents negative correlation. FDR correction method was used for significance test. *. P<0.05圖6 lncRNAs的反式差異表達(dá)靶基因的表達(dá)水平與肉質(zhì)性狀的關(guān)聯(lián)分析Fig.6 Correlation analysis between expression of trans differentially expressed target genes of lncRNAs and meat quality traits

4 結(jié) 論

藏豬和大約克夏豬背最長肌中部分lncRNAs和基因存在顯著性差異，推測這些lncRNAs在肌肉品質(zhì)的調(diào)節(jié)中發(fā)揮重要作用。MSTRG.12095.1、MSTRG.151.1、MSTRG.42191.1、MSTRG.3334.1等一些lncRNAs可能通過調(diào)節(jié)TGFB2、HES1、LEP和PER2等基因的表達(dá)來影響肌肉生長發(fā)育和肌內(nèi)脂肪沉積的過程，從而導(dǎo)致藏豬和大約克夏豬在肌肉生長發(fā)育和肌內(nèi)脂肪含量方面的差異。本研究從lncRNA的角度為豬肉質(zhì)改良工作中候選功能基因的研究提供了理論依據(jù)。