榮麗云，黃 寧，王孝義，李明麗，陳 強 ，魯紹雄

(云南農(nóng)業(yè)大學(xué) 動物科學(xué)技術(shù)學(xué)院，云南 昆明 650201)

肉質(zhì)性狀是養(yǎng)豬生產(chǎn)中最重要的經(jīng)濟(jì)性狀之一，對養(yǎng)豬生產(chǎn)效益有著重要影響。提高豬肉品質(zhì)一直是豬育種工作的重要目標(biāo)，但由于肉質(zhì)性狀活體難以準(zhǔn)確測定，且包括肌內(nèi)脂肪含量、大理石紋、嫩度、pH 值和滴水損失等諸多相互關(guān)聯(lián)的子性狀，常規(guī)育種方法改良難度大，難以獲得理想的遺傳進(jìn)展。隨著分子生物學(xué)及基因組學(xué)相關(guān)技術(shù)的發(fā)展，人們發(fā)掘了大量與豬肉質(zhì)性狀相關(guān)的基因[1-4]。盡管如此，關(guān)于豬肉質(zhì)性狀的分子機制迄今仍不甚清楚。

基因共表達(dá)網(wǎng)絡(luò)分析(gene co-expression network analysis，GCNA)是基于對同一性狀的基因在表達(dá)模式上具有共表達(dá)特征這一設(shè)定，從構(gòu)建的表達(dá)調(diào)控網(wǎng)絡(luò)中發(fā)現(xiàn)并識別目標(biāo)性狀的相關(guān)基因。該方法對目標(biāo)性狀關(guān)鍵新基因的發(fā)掘及其功能預(yù)測具有重要作用。加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(weighted gene co-expression network analysis，WGCNA)是目前開展GCNA 最主要的方法[5]，人們利用該方法鑒定了大量與復(fù)雜性狀相關(guān)的基因共表達(dá)模塊和關(guān)鍵基因[6-9]。

近年來，利用GCNA 方法探索畜禽重要經(jīng)濟(jì)性狀分子機制的研究逐漸增多，有關(guān)豬肉質(zhì)性狀的基因共表達(dá)模塊研究也有陸續(xù)報道[10-13]。如ZAPPATERRA 等[2]對意大利大白豬半膜肌進(jìn)行研究，識別了4 個與肌內(nèi)脂肪含量相關(guān)的基因共表達(dá)模塊；ZHAO 等[14]利用WGCNA 法識別了與豬肉滴水損失相關(guān)的AASS、BCKDHB、ALDH6A1、MUT和MCCC1等5 個基因。這些研究為進(jìn)一步揭示豬肉質(zhì)性狀形成的基因表達(dá)機制提供了有益探索，但目前關(guān)于肉質(zhì)性狀相關(guān)的基因共表達(dá)研究及可利用的基因共表達(dá)信息仍較少，尚需進(jìn)一步深入研究。

相較于國外瘦肉型豬種，中國地方豬種肉質(zhì)優(yōu)良的遺傳特性突出。由于兩類豬種在肉質(zhì)性狀等方面存在較大差異，基于兩類豬種的比較分析是揭示肉質(zhì)性狀分子機制的一種有效途徑。本研究基于公共GEO 數(shù)據(jù)庫中大蒲蓮豬(中國地方品種)和長白豬(國外瘦肉型品種)背最長肌基因表達(dá)數(shù)據(jù)集，通過WGCNA、蛋白質(zhì)與蛋白質(zhì)互作(protein-protein interaction，PPI)網(wǎng)絡(luò)及中心性分析等方法，挖掘與肉質(zhì)性狀相關(guān)的關(guān)鍵基因共表達(dá)模塊和關(guān)鍵基因，并對識別的關(guān)鍵基因表達(dá)進(jìn)行qPCR 驗證，以探索豬肉質(zhì)性狀形成的分子機制，為聯(lián)合運用經(jīng)典育種和分子育種手段改良豬肉質(zhì)性狀奠定基礎(chǔ)并提供科學(xué)依據(jù)。

1 材料與方法

1.1 數(shù)據(jù)來源

數(shù)據(jù)來源于美國生物技術(shù)信息中心(NCBI)基因表達(dá)數(shù)據(jù)庫(GEO，https://www.ncbi.nlm.nih.gov/geoprofiles/)中的數(shù)據(jù)集GSE119368，包括8 頭大蒲蓮豬(中國地方品種)和8 頭長白豬(國外瘦肉型品種) 6 月齡背最長肌基因表達(dá)數(shù)據(jù)。豬群飼料營養(yǎng)相同、飼養(yǎng)管理條件相近?；虮磉_(dá)數(shù)據(jù)由Illumina HiSeq 2500 測序平臺產(chǎn)生。

1.2 差異表達(dá)基因分析

差異表達(dá)基因分析利用基于R 軟件的limma包[15](version 3.44.3)進(jìn)行，以矯正的P＜0.05 和|log2fold change|＞1 作為基因差異表達(dá)的顯著性閾值。差異表達(dá)基因(differentially expressed genes,DEGs)分布的火山圖和表達(dá)熱圖采用ggplot2 軟件包(https://cran.r-project.org/web/packages/ggplot2/,version 3.3.3)繪制。

1.3 基因共表達(dá)模塊分析

基因共表達(dá)模塊分析采用WGCNA 方法[5]，使用基于R 環(huán)境的WGCNA 軟件包(https://cran.r-project.org/web/packages/WGCNA/，version 1.70-3)進(jìn)行。以大蒲蓮豬為處理組、長白豬為對照組，質(zhì)控后獲得的12 494 個基因按如下流程進(jìn)行分析：

(1)計算兩兩基因間的皮爾遜相關(guān)系數(shù)(r)，并據(jù)此構(gòu)建基因表達(dá)相似矩陣；(2)將基因表達(dá)相似矩陣轉(zhuǎn)換成鄰接矩陣，并通過軟閾值β (β=10，無尺度拓?fù)錁?biāo)準(zhǔn)R2＞0.85)產(chǎn)生增強鄰接矩陣；(3)使用皮爾遜相關(guān)分析構(gòu)建增強鄰接矩陣的無監(jiān)督基因共表達(dá)關(guān)系矩陣，并基于基因間的關(guān)聯(lián)程度指標(biāo)拓?fù)渲丿B尺度(topological overlap measure，TOM)將其轉(zhuǎn)化為拓?fù)渚仃嚕?4)基于基因間相異度(1-TOM)，使用平均連鎖層次聚類法將具有相似基因表達(dá)模式的基因歸為同一個基因模塊，采用動態(tài)剪切算法從系統(tǒng)聚類樹中識別基因共模塊，并將相似度高于75%的基因模塊進(jìn)行合并。

1.4 基因功能分析

差異表達(dá)基因及關(guān)鍵模塊中基因的功能采用基因本體論(gene ontology，GO)富集及京都基因和基因組百科全書(Kyoto encyclopedia of genes and gnomes，KEGG)通路富集方法，利用在線的DAVID 軟件[16](https://david.ncifcrf.gov/，version 6.8)進(jìn)行分析，通路顯著性富集閾值為P＜0.05。

1.5 PPI 網(wǎng)絡(luò)構(gòu)建

關(guān)鍵模塊中編碼蛋白基因間的互作關(guān)系采用PPI 網(wǎng)絡(luò)進(jìn)行分析，互作關(guān)系信息來源于在線數(shù)據(jù)庫STRING[17](http://string-db.org,version 11.5)，篩選標(biāo)準(zhǔn)為互作得分大于0.7。PPI 網(wǎng)絡(luò)采用Cytoscape 軟件[18](https://cytoscape.org/,version 3.8.0)構(gòu)建。

1.6 高度關(guān)聯(lián)模塊及關(guān)鍵基因識別

PPI 網(wǎng)絡(luò)中與肉質(zhì)性狀高度關(guān)聯(lián)模塊的分析采用分子復(fù)合檢測(molecular complex detection，MCODE)算法，并利用Cytoscape 軟件中的MCODE 插件[19](version 1.5.1)進(jìn)行，相應(yīng)參數(shù)設(shè)置為Degree Cutoff=2、Max Depth=100、K-Core=2 以及Node Score Cutoff=0.2。關(guān)鍵基因的識別采用Subgraph、Eigenvector、Information、Closeness、Betweenness、LAC 和Network 等7 種中心性方法，并利用Cytoscape 軟件中的CytoNCA 插件[20](version 2.1.6)進(jìn)行，其中，7 種方法得分最高的前10 個基因即為關(guān)鍵基因。

1.7 關(guān)鍵基因驗證

選擇彼此無親緣關(guān)系的地方品種撒壩豬和引進(jìn)瘦肉型品種大白豬各3 頭，在相同飼養(yǎng)管理條件下飼養(yǎng)至體質(zhì)量為90～100 kg 時屠宰，采集背最長肌組織用于肉質(zhì)關(guān)鍵基因驗證。從每個識別的肉質(zhì)性狀高度關(guān)聯(lián)模塊中各選擇3 個關(guān)鍵基因，采用qPCR 方法檢測關(guān)鍵基因在撒壩豬和大白豬背最長肌組織中的表達(dá)水平，并采用非配對試驗的t檢驗進(jìn)行差異分析，差異顯著和極顯著閾值分別為P＜0.05 和P＜0.01。

2 結(jié)果與分析

2.1 DEGs 識別及功能分析

差異基因表達(dá)分析結(jié)果如圖1 所示。在大蒲蓮豬和長白豬的背最長肌組織中共識別出390 個DEGs，其中，大蒲蓮豬中顯著高表達(dá)基因189 個、顯著低表達(dá)基因201 個。

圖1 大蒲蓮和長白豬差異表達(dá)基因分布及表達(dá)水平Fig.1 Distribution and expression level of differentially expressed genes in Dapulian and Landrace pigs

GO 分析結(jié)果(表1)顯示：189 個高表達(dá)基因極顯著富集在3 個GO 生物學(xué)過程、3 個GO 細(xì)胞組分和2 個GO 分子功能上，201 個低表達(dá)基因極顯著富集在12 個GO 生物學(xué)過程、5 個GO細(xì)胞組分和2 個GO 分子功能上。高表達(dá)基因主要與膠原纖維組織、成肌細(xì)胞融合和PI3K-Akt信號傳導(dǎo)等生物學(xué)過程相關(guān)，低表達(dá)基因主要與骨骼肌細(xì)胞分化、蛋白質(zhì)折疊及神經(jīng)元凋亡過程調(diào)控等生物學(xué)過程相關(guān)。

表1 差異表達(dá)基因極顯著富集的GO 條目Tab.1 GO terms with extremely significant enrichment of differentially expressed genes

KEGG 通路分析結(jié)果(表2)顯示：189 個高表達(dá)基因和201 個低表達(dá)基因分別極顯著富集在5 個和3 個KEGG 通路上。其中，高表達(dá)基因主要與蛋白質(zhì)消化吸收、醛固酮調(diào)節(jié)的鈉重吸收以及心肌細(xì)胞的腎上腺素能信號等通路相關(guān)，低表達(dá)基因主要與MAPK 信號通路、內(nèi)質(zhì)網(wǎng)蛋白質(zhì)加工以及雌激素信號通路等相關(guān)。

表2 差異表達(dá)基因極顯著富集的KEGG 通路Tab.2 KEGG pathways with extremely significant enrichment of differentially expressed genes

2.2 基因共表達(dá)模塊分析

聚類結(jié)果顯示：大蒲蓮豬和長白豬樣本均各自聚成一類(圖2a)。根據(jù)計算的基因間相異度，識別了6 個基因共表達(dá)模塊(圖2b)；合并相似性高于75%的模塊后，獲得紅色、藍(lán)色和綠色3 個基因共表達(dá)模塊(圖2c)。其中，紅色模塊(r=0.91，P=7E-7)和藍(lán)色模塊(r=0.89，P=5E-6)與肉質(zhì)性狀的關(guān)聯(lián)度最高(圖2d)，兩模塊中基因表達(dá)與肉質(zhì)均呈極顯著相關(guān)(r=0.93，P＜1E-200；r=0.90，P＜1E-200)(圖2e 和2f)。

圖2 WGCNA 分析結(jié)果Fig.2 The results of WGCNA analysis

2.3 與肉質(zhì)性狀最顯著關(guān)聯(lián)的模塊基因功能分析

與肉質(zhì)性狀最顯著關(guān)聯(lián)的紅色模塊660 個基因極顯著(P＜0.01)富集到17 個GO 生物學(xué)過程和12 個KEGG 通路上(顯著性最高的前5 個生物學(xué)過程和KEGG 通路見表3)，主要與蛋白質(zhì)折疊、骨骼肌細(xì)胞分化、糖原代謝過程和雌激素信號、蛋白質(zhì)加工以及蛋白激酶等信號通路相關(guān)；藍(lán)色模塊527 個基因極顯著(P＜0.01)富集到12個GO 生物學(xué)過程和10 個KEGG 通路上(顯著性最高的前5 個生物學(xué)過程和KEGG 通路見表4)，主要與肌肉組織形態(tài)形成、心臟肌肉收縮、快慢纖維轉(zhuǎn)換過程和蛋白質(zhì)消化吸收、ECM-受體相互作用以及黏著斑等信號通路相關(guān)。

表3 紅色模塊基因極顯著富集的前5 個GO 條目和KEGG 通路Tab.3 Top 5 GO terms and KEGG pathways extremely significantly enriched by the genes in the red module

表4 藍(lán)色模塊基因極顯著富集的前5 個GO 條目和KEGG 通路Tab.4 Top 5 GO terms and KEGG pathways extremely significantly enriched by the genes in the blue module

2.4 PPI 網(wǎng)絡(luò)構(gòu)建及關(guān)鍵基因識別

所構(gòu)建的紅色模塊基因PPI 網(wǎng)絡(luò)包括163 個基因(節(jié)點)和371 條邊(圖3a)，最高得分的子網(wǎng)絡(luò)包含11 個基因(FBXL4、RNF115、FBXO40、ASB4、ASB15、UBE2L6、RNF41、ASB5、HERC5、SOCS1、ASB11)和55 條邊(圖3b)，11 個基因的中心性綜合得分相同(表5)，為肉質(zhì)性狀相關(guān)的關(guān)鍵基因。

表5 紅色模塊基因PPI 子網(wǎng)絡(luò)基因中心性分析Tab.5 Centrality analysis of genes in the PPI sub-network of the red module

圖3 基因互作網(wǎng)絡(luò)與高關(guān)聯(lián)度模塊識別Fig.3 Gene interaction network and highly related module identification

所構(gòu)建的藍(lán)色模塊基因PPI 網(wǎng)絡(luò)包含128 個基因和341 條邊(圖3c)，最高得分的子網(wǎng)絡(luò)包含22 個基因和112 條邊(圖3d)，中心性綜合得分最高的前10 個基因為COL3A1、COL5A1、COL1A2、COL14A1、CRTAP、COL5A3、LEPREL1、DEPTOR、COL4A2和COL6A2(表6)，為肉質(zhì)性狀相關(guān)的關(guān)鍵基因。

表6 藍(lán)色模塊基因PPI 子網(wǎng)絡(luò)基因中心性分析Tab.6 Centrality analysis of genes in the PPI sub-network of the blue module

2.5 關(guān)鍵基因表達(dá)分析及qPCR 驗證

由圖4 可知：2 個模塊各挑選的3 個肉質(zhì)性狀關(guān)鍵基因在大蒲蓮豬和長白豬背最長肌中的表達(dá)水平均呈極顯著差異，除RNF115(P=0.001 84)外，COL3A1(P=0.000 16)、COL5A1(P=0.001 65)、COL1A2(P=0.000 073)、FBXL4(P=0.000 97)和FBXO40(P=0.000 54) 5 個基因的表達(dá)水平均呈現(xiàn)為長白豬高于大蒲蓮豬。

圖4 大蒲蓮和長白豬背最長肌中6 個關(guān)鍵基因的表達(dá)水平Fig.4 The expression levels of six key genes in longissimus dorsi muscle of Dapulian and Landrace pigs

qPCR 結(jié)果(圖5)顯示：藍(lán)色模塊中的COL-3A1(P＜0.000 1)、COL5A1(P=0.000 6)和COL1A2(P＜0.000 1)在撒壩豬和大白豬背最長肌中的表達(dá)水平呈極顯著差異，而紅色模塊中3 個基因(RNF115、FBXL4和FBXO40)的表達(dá)水平在2 個品種間則無顯著差異(P＞0.05)。

圖5 6 個關(guān)鍵基因在撒壩和大白豬背最長肌中的表達(dá)水平Fig.5 The expression of six genes in longissimus dorsi muscle of Saba and Large White pigs

3 討論

本研究基于WGCNA 方法識別了2 個與豬肉質(zhì)性狀顯著關(guān)聯(lián)的基因共表達(dá)模塊，并通過構(gòu)建的模塊基因間互作關(guān)系網(wǎng)絡(luò)分別識別出10 個和11 個關(guān)鍵基因。

在藍(lán)色模塊識別的10 個豬肉質(zhì)性狀關(guān)鍵基因中，有7 個基因(COL3A1、COL5A1、COL1A2、COL14A1、COL5A3、COL4A2和COL6A2)為膠原蛋白家族相關(guān)基因，該家族基因編碼一類膠原蛋白，主要存在于動物結(jié)締組織中。有研究表明：膠原蛋白對動物肌內(nèi)脂肪形成、肌肉嫩度、系水力及滴水損失等具有重要影響[21-24]。NAKAJIMA等[25]研究了牛胸最長肌的肌內(nèi)脂肪，發(fā)現(xiàn)膠原蛋白V 和VI 型對肌內(nèi)脂肪的形成具有正向調(diào)控作用；曾勇慶等[26]對萊蕪豬背最長肌的膠原蛋白含量及性質(zhì)進(jìn)行分析后發(fā)現(xiàn)：隨著膠原蛋白含量的逐漸增加，豬肉的大理石紋和持水性能逐漸增加。

已有研究[27-29]表明：COL3A1、COL5A1、COL1A2和COL14A1等4 個膠原蛋白家族基因與豬肉質(zhì)性狀相關(guān)。LIM 等[27]對巴克夏豬背最長肌轉(zhuǎn)錄組的研究顯示：COL1A2、COL5A1和COL14-A1與豬肌肉肌內(nèi)脂肪含量顯著相關(guān)；BAO 等[28]對大白豬雜交群體的研究發(fā)現(xiàn)：COL3A1基因的表達(dá)上調(diào)有利于豬肌纖維的生成；FERNáNDEZ-BARROSO 等[29]對伊比利亞豬的研究發(fā)現(xiàn)：COL14A1基因?qū)ωi肉嫩度有顯著影響?？梢姡z原蛋白家族基因可能主要與豬肉的肌內(nèi)脂肪含量和嫩度等相關(guān)。本研究的qPCR 驗證結(jié)果表明：COL3A1、COL5A1和COL1A2基因在肉質(zhì)差異較大的地方品種撒壩豬和引進(jìn)品種大白豬背最長肌中的表達(dá)水平差異極顯著，進(jìn)一步提示上述基因?qū)θ赓|(zhì)性狀有著潛在重要影響。而有關(guān)COL-5A3、COL4A2和COL6A2基因?qū)ωi肉質(zhì)性狀的影響迄今未見相關(guān)研究報道，其確切功能尚需進(jìn)一步研究。此外，本研究識別的另外3 個關(guān)鍵基因(CRTAP、LEPREL1和DEPTOR)中，CRTAP是軟骨關(guān)聯(lián)蛋白編碼基因，主要與骨質(zhì)發(fā)育相關(guān)，有研究表明該基因參與膠原形成和胞外基質(zhì)降解[30]，而LEPREL1和DEPTOR基因?qū)θ赓|(zhì)性狀的影響目前仍未見報道，2 個基因與肉質(zhì)的確切相關(guān)性仍需進(jìn)一步通過功能試驗確認(rèn)。

在紅色模塊識別的11 個關(guān)鍵基因中，已有研究[31-37]顯示其中的一些基因與肉質(zhì)性狀相關(guān)，如FBXL4、FBXO40、ASB4和ASB15等。FBXO40為F-box 蛋白家族成員之一，該基因的敲除可顯著促進(jìn)豬骨骼肌纖維增粗[31-32]。FBXL4是一個與線粒體功能控制有關(guān)的核基因，與FBXO40同屬于F-box 蛋白家族，當(dāng)FBXL4基因敲除時小鼠出現(xiàn)線粒體功能障礙及體質(zhì)量減輕等現(xiàn)象[33]。ASB4、ASB15和ASB5同屬ASB 基因家族成員，在骨骼肌和心肌中呈高表達(dá)，被認(rèn)為是豬肉質(zhì)性狀遺傳改良的潛在候選基因[34-35]。UBE2L6在動物脂肪細(xì)胞中可作為介導(dǎo)脂肪分解的負(fù)調(diào)節(jié)因子促進(jìn)高脂飲食誘導(dǎo)產(chǎn)生肥胖，敲除該基因能夠抑制脂肪生成和減少脂肪細(xì)胞數(shù)量[36-37]。目前，關(guān)于RNF41、RNF115、HERC5和SOCS1基因與肉質(zhì)性狀的相關(guān)性尚未見報道。

盡管撒壩豬和大白豬背最長肌中FBXL4和FBXO40基因的表達(dá)差異與大蒲蓮和長白豬的趨勢一致，但RNF115基因的表達(dá)差異則呈相反趨勢，且這3 個關(guān)鍵基因在撒壩豬和大白豬中的表達(dá)水平差異均不顯著。究其原因，可能在于撒壩豬和大蒲蓮豬雖同為肉質(zhì)優(yōu)良的中國地方豬種，長白豬和大白豬均為國外瘦肉型豬種，但同類豬種在部分肉質(zhì)指標(biāo)上仍存較大差異，如體質(zhì)量為100 kg 的撒壩豬肌內(nèi)脂肪含量和肉色評分分別為(6.53±1.65)%和(3.19±0.31)分[38]，而大蒲蓮豬則分別為(4.63±0.88)%和(4.25±0.14)分[39]，可見其肉質(zhì)形成的分子機制可能存在一定差異，這可能也是造成驗證結(jié)果不一致的原因之一。如能在不同豬種中進(jìn)一步開展表達(dá)水平驗證，重點比較不同地方品種(如撒壩豬和大蒲蓮豬)的基因表達(dá)差異情況，深入揭示這些基因影響的具體性狀及其效應(yīng)，將有助于豬肉質(zhì)性狀的基因發(fā)掘并為精準(zhǔn)育種提供更有價值的分子標(biāo)記。

4 結(jié)論

本研究基于WGCNA 方法識別了與豬肉質(zhì)性狀相關(guān)的2 個關(guān)鍵基因共表達(dá)模塊和21 個關(guān)鍵基因，其中7 個膠原蛋白基因家族相關(guān)基因COL3A1、COL5A1、COL1A2、COL14A1、COL5A3、COL4A2和COL6A2在豬肉質(zhì)性狀形成中可能有重要作用。大蒲蓮豬和撒壩豬雖同為地方豬種，但其肉質(zhì)性狀形成的分子機制可能存在較大差異。