利用加權(quán)基因共表達網(wǎng)絡(luò)分析篩選天農(nóng)麻雞胴體性狀候選基因

趙 迪，康慧敏，譚曉冬，劉冉冉，張正芬，李 華,，趙桂蘋*

(1. 佛山科學(xué)技術(shù)學(xué)院生命科學(xué)與工程學(xué)院 廣東省動物分子設(shè)計與精準(zhǔn)育種重點實驗室，佛山 528225；2. 中國農(nóng)業(yè)科學(xué)院北京畜牧獸醫(yī)研究所，北京 100193；3.廣東天農(nóng)食品集團股份有限公司，清遠 511827)

中國是全球第二大雞肉生產(chǎn)國和消費國，雞肉產(chǎn)量約占全球雞肉產(chǎn)量的14.7%，在中國，雞肉是第二大肉類產(chǎn)品，占肉類總產(chǎn)量的20.5%。隨著生活水平的提高，人們對雞肉的需求正在從“數(shù)量”型向“質(zhì)量”型轉(zhuǎn)變。黃羽肉雞肉質(zhì)細膩，風(fēng)味鮮美，深受我國消費者喜愛。清遠麻雞是慢速型黃羽肉雞的代表性品種，也是廣東最著名、最具有代表性的優(yōu)質(zhì)地方雞種。天農(nóng)麻雞是以清遠麻雞為主要素材培育的優(yōu)質(zhì)肉雞配套系，商品代抗病力強，肉質(zhì)松軟、風(fēng)味好、雞味濃、口感佳，體型外貌符合飼養(yǎng)地區(qū)消費市場對特優(yōu)型地方雞種的要求。無論是父母代種雞還是商品代仔雞，它們在同類產(chǎn)品中性能均處于領(lǐng)先水平，具有較強的市場競爭力。

胴體性狀是肉雞的重要經(jīng)濟性狀，與肌肉、骨骼等多組織發(fā)育緊密相關(guān)，提高黃羽肉雞的生產(chǎn)效率，增加產(chǎn)肉量一直是遺傳育種工作的核心內(nèi)容。隨著高通量測序技術(shù)的快速發(fā)展，轉(zhuǎn)錄組測序(RNA-Seq)已被廣泛應(yīng)用于挖掘畜禽重要經(jīng)濟性狀的候選基因。李國輝利用金茅黑雞胸肌組織進行RNA-Seq，通過關(guān)聯(lián)分析鑒定到107基因是影響體重和屠體重的候選基因，通過調(diào)節(jié)核質(zhì)間的物質(zhì)運輸，參與生物的多種生長發(fā)育過程。Zhang等對白羽肉雞、大恒肉雞和商品代羅曼蛋雞的胸肌進行RNA-Seq，共獲得8 398個差異表達基因，發(fā)現(xiàn)與肌肉生長相關(guān)的基因15、2、3、2、-2、、等在不同品種中的表達量差異顯著。

以RNA-Seq數(shù)據(jù)為基礎(chǔ)，開發(fā)的加權(quán)基因共表達網(wǎng)絡(luò)分析(weighted gene co-expression network analysis, WGCNA)是一種高效、準(zhǔn)確的生物信息學(xué)和生物數(shù)據(jù)挖掘的方法。它可以系統(tǒng)地篩選目標(biāo)性狀關(guān)聯(lián)的共表達基因，以便于進一步了解基因相互作用的機制。同時，利用WGCNA能夠發(fā)現(xiàn)與目標(biāo)性狀高度相關(guān)的模塊，快速識別目標(biāo)性狀的關(guān)鍵候選基因，進一步通過不同的篩選標(biāo)準(zhǔn)鑒定樞紐基因。WGCNA已在動物遺傳學(xué)等生物領(lǐng)域廣泛應(yīng)用。在黃羽肉雞的研究中，本課題前期曾報道，通過對京星黃雞進行WGCNA，鑒定到轉(zhuǎn)錄因子3MBTL1和轉(zhuǎn)錄因子輔助因子1、1和6與高肌內(nèi)脂肪和低腹脂重相關(guān)。采用WGCNA探究京星黃雞肝中基因表達量與腹脂重的相關(guān)性，其中8、1、2和2基因被鑒定為腹脂重正相關(guān)模塊中的樞紐基因，表明這些基因在肝脂質(zhì)代謝中發(fā)揮著重要作用。

已有的關(guān)于慢速型黃羽肉雞胴體性狀遺傳基礎(chǔ)的研究相對較少，其候選基因仍有待進一步挖掘。本研究基于104只上市日齡天農(nóng)麻雞胸肌轉(zhuǎn)錄組數(shù)據(jù)進行WGCNA，篩選與胴體性狀顯著相關(guān)的共表達模塊及其樞紐基因，并對模塊內(nèi)基因進行富集分析，挖掘影響胴體性狀的作用通路，為解析黃羽肉雞胴體性狀的調(diào)控機制提供參考。

1 材料與方法

1.1 試驗動物

本試驗群體來源于廣東天農(nóng)食品集團股份有限公司提供的商品代天農(nóng)麻雞，該群體是慢速型優(yōu)質(zhì)黃羽肉雞，1～125日齡于農(nóng)戶家中散養(yǎng)，飼養(yǎng)期間自由飲水和采食，所有個體健康狀況良好。

1.2 樣品采集及表型測定

天農(nóng)麻雞飼養(yǎng)至126日齡，宰前禁食12 h后測定活重，并隨機選擇104只母雞進行頸動脈放血致死。統(tǒng)一采集右側(cè)胸肌相同部位樣品，置于液氮中迅速冷凍，-80 ℃保存?zhèn)溆?。參照《家禽生產(chǎn)性能名詞術(shù)語和度量統(tǒng)計方法》(NY/T823—2020)測定屠體重、全凈膛重和胸肌重，并計算屠體率、全凈膛率和胸肌率：屠體率(%)=屠體重/活重×100；全凈膛率(%)=全凈膛重/活重×100；胸肌率(%)=兩側(cè)胸肌重/全凈膛重×100。

1.3 胸肌組織RNA提取

剪取黃豆大小的胸肌樣本，使用TRIzol試劑(Invitrogen，美國)提取總RNA。使用K5500分光光度計(凱奧，中國)檢測RNA純度。同時使用生物分析儀2100系統(tǒng)的RNA Nano 6000檢測試劑盒(安捷倫科技，美國)評估RNA的完整性和濃度。并將A260/A280比值在1.8～2.0之間，RNA完整性大于7.5的樣品用于后續(xù)試驗。

1.4 RNA-Seq數(shù)據(jù)分析

首先利用FastQC軟件對下機數(shù)據(jù)進行質(zhì)量評估并利用BBMap軟件去掉read接頭，通過Fastp過濾不合格read，獲得clean data。利用HISAT2建立基因組(GRCg6A)索引，并進行比對，通過SAMtools軟件將比對后的sam文件轉(zhuǎn)換為bam文件，并進行排序、質(zhì)控、索引，使用FeatureCount軟件對基因進行定量，采用DESeq2軟件對基因表達量進行標(biāo)準(zhǔn)化。

1.5 加權(quán)基因共表達網(wǎng)絡(luò)分析

使用R語言WGCNA包的hclust函數(shù)對樣本進行聚類，剔除離群樣本。隨后使用pick-SoftThreshold(sft)函數(shù)選擇最佳軟閾值，以無標(biāo)度拓?fù)鋽M合指數(shù)()大于0.85為條件確定最佳軟閾值(sft=7)。隨后采用adjacency函數(shù)構(gòu)建拓?fù)渲丿B矩陣(topological overlap matrix, TOM)，并根據(jù)TOM矩陣的相異度，使用hclust函數(shù)對基因進行聚類。采用主成分分析法計算每個模塊的特征值(module eigengene, ME)，根據(jù)ME進行聚類，對相似度大于0.7的模塊進行合并。

將合并過的模塊重新計算ME，與表型矩陣進行關(guān)聯(lián)分析，選擇與表型最相關(guān)的模塊為目標(biāo)模塊。隨后計算目標(biāo)模塊內(nèi)的Module membership(MM)與Gene significance(GS)，以|MM| > 0.8、|GS| > 0.2為標(biāo)準(zhǔn)篩選到的基因為樞紐基因(hub gene)，采用DAVID數(shù)據(jù)庫(http://david.abcc.ncifcrf.gov/)對目標(biāo)模塊內(nèi)的基因進行KEGG富集分析，以<0.05為標(biāo)準(zhǔn)篩選顯著富集通路。

1.6 樞紐基因蛋白互作網(wǎng)絡(luò)構(gòu)建及可視化

采用String在線網(wǎng)絡(luò)(https://string-db.org/)對每一個及全部樞紐基因分別進行蛋白質(zhì)互作網(wǎng)絡(luò)(protein-protein interaction network, PPI network)的構(gòu)建，選取為參考物種。并使用Cytoscape軟件中的cytoHubba插件，根據(jù)PPI網(wǎng)絡(luò)節(jié)點中的Degree進行篩選，對排名前20的基因進行可視化。

1.7 統(tǒng)計分析

計算104個個體7種胴體表型的平均值及標(biāo)準(zhǔn)差，均以“平均值±3倍標(biāo)準(zhǔn)差”為標(biāo)準(zhǔn)剔除離群個體，利用SPSS軟件計算58個樞紐基因的表達量與7種胴體表型的Spearman相關(guān)性，使用GraphPad Prism 8軟件作圖。

2 結(jié) 果

2.1 表型數(shù)據(jù)統(tǒng)計

對104只126日齡天農(nóng)麻雞的活重、屠體重等7種胴體表型進行基本的描述性統(tǒng)計(表1)，126日齡母雞平均活重為1 541.3 g，屠體重均值為1 370.5 g，全凈膛重均值為1 040.0 g，胸肌重均值為168.9 g，屠體率均值為89.0%，全凈膛率均值為67.7%，屠宰性能良好。

表1 胴體性狀描述性統(tǒng)計Table 1 Descriptive statistics for carcass traits

2.2 基因共表達模塊的構(gòu)建

對轉(zhuǎn)錄組數(shù)據(jù)進行分析質(zhì)控后，共得到15 092個基因用于后續(xù)分析。對104個樣本進行聚類，剔除5個離群樣本。以> 0.85為標(biāo)準(zhǔn)確定最佳軟閾值sft=7(圖1a)，共得到24個模塊。根據(jù)ME進行模塊聚類，合并相似度大于0.7的模塊(圖1b)，最后得到19個模塊(圖1c)。將共表達趨勢不明顯的基因統(tǒng)一歸入Grey模塊。

a.左圖中縱軸代表對應(yīng)的網(wǎng)絡(luò)中l(wèi)og(k)與log(p(k))相關(guān)系數(shù)的平方；右圖中縱軸代表相應(yīng)的基因模塊中基因鄰接函數(shù)的均值；b.模塊聚類圖，紅色線表示相似度達0.7；c.基因動態(tài)剪切聚類樹，每個顏色代表一個模塊，第一行為首次聚類的結(jié)果，第二行為模塊合并后的結(jié)果a. The vertical axis of the left figure represents the square of the correlation coefficient between log(k) and log(p(k)) in the corresponding network; the vertical axis of the right figure represents the mean value of the gene adjacency function in the corresponding gene module; b. Module clustering graph, the height of the red line is 0.3, and the modules below the line are the modules that have greater similarity and need to be merged; c. Gene dynamic shear clustering tree, each color represents a module, the color in the first row is the result of the first clustering, the color in the second row is the result of the merging of the modules圖1 加權(quán)基因共表達網(wǎng)絡(luò)分析Fig.1 Weighted gene co-expression network analysis

2.3 目標(biāo)模塊的選擇與分析

模塊的 ME與表型矩陣進行相關(guān)性分析，得到模塊與胴體性狀間的相關(guān)系數(shù)及值(圖2)，選擇與目標(biāo)性狀相關(guān)性較強(>0.39)且<0.05的模塊為目標(biāo)模塊。最終確定4個目標(biāo)模塊，分別是“Magenta”模塊，該模塊對LW、DW、EW和BMW均呈顯著負(fù)相關(guān)(LW:=-0.55,=4×10; DW:=-0.51,=9×10; EW:=-0.54,=1×10; BMW:=-0.50,=1×10)、“Tan”模塊與DP呈顯著正相關(guān)(DP:=0.47,=9×10)、“Green”模塊與DP呈顯著負(fù)相關(guān)(DP:=-0.58,=2×10)和“Purple” 模塊與EP呈正相關(guān)(EP:=0.39,=6×10)，這4個模塊包含的基因數(shù)分別為：504、439、753、474個(表2)。根據(jù)|MM| > 0.8且|GS| > 0.2的標(biāo)準(zhǔn)確定模塊的樞紐基因展開描述。

模塊與胴體性狀的相關(guān)性分析，紅色代表正相關(guān)，藍色代表負(fù)相關(guān)，顏色越深相關(guān)性越強The correlation analysis between the module and the carcass traits, red represents positive correlation, blue represents negative correlation, the darker the color, the stronger the correlation圖2 模塊性狀相關(guān)圖Fig.2 Module-trait relationship diagram

表2 目標(biāo)模塊及樞紐基因Table 2 Target modules and hub genes

2.4 目標(biāo)模塊中基因的KEGG通路富集

KEGG富集分析表明，目標(biāo)模塊共富集到18條 通路(圖3)。Magenta模塊內(nèi)基因顯著富集于mRNA 監(jiān)測、Notch信號、黑色素生成、Wnt信號通路和TGF-β信號通路；Tan模塊內(nèi)基因顯著富集于脂肪細胞因子信號通路、肌動蛋白細胞骨架的調(diào)節(jié)和FoxO信號通路；Green模塊顯著富集的KEGG通路包括細胞周期、范可尼貧血途徑、基因復(fù)制、核苷酸切除修復(fù)、同源重組、錯配修復(fù)和嘧啶代謝通路；Purple模塊內(nèi)基因顯著富集于ECM-受體相互作用、黏著力和細胞黏附分子(CAMs)通路。富集結(jié)果表明，共有10條與肌肉生長發(fā)育相關(guān)的信號通路(表3)。

氣泡大小表示基因數(shù)量；顏色表示富集分析的顯著性The size of the bubbles indicate the number of genes; the color indicates the significance of the enrichment analysis圖3 目標(biāo)模塊KEGG富集分析Fig.3 KEGG enrichment analysis of target modules

表3 目標(biāo)模塊KEGG富集結(jié)果Table 3 KEGG enrichment results of target modules

2.5 樞紐基因PPI網(wǎng)絡(luò)構(gòu)建

對樞紐基因進行PPI網(wǎng)絡(luò)構(gòu)建，以0.4為互作標(biāo)準(zhǔn)進行篩選。結(jié)果顯示(圖4a)，前20個蛋白之間形成2個子網(wǎng)絡(luò)，8個基因的子網(wǎng)絡(luò)均來自Green模塊，12個基因的子網(wǎng)絡(luò)均來自于Purple模塊。其中12個基因的子網(wǎng)絡(luò)存在強互作，對子網(wǎng)絡(luò)中的每個基因分別構(gòu)建PPI網(wǎng)絡(luò)，發(fā)現(xiàn)COL1A2、COL6A1、COL6A2、COL6A3和SPARC的互作蛋白大多數(shù)包括在1個子網(wǎng)絡(luò)內(nèi)，表明這5個基因可能在該網(wǎng)絡(luò)中發(fā)揮核心作用，作為重要候選基因(圖4b～f)。

a.前20個樞紐基因PPI網(wǎng)絡(luò)互作圖；b～f.分別為COL1A2、COL6A1、COL6A2、COL6A3和SPARC基因PPI網(wǎng)絡(luò)互作圖a. Top 20 PPI network hub genes interaction diagram; b-f. PPI network COL1A2, COL6A1, COL6A2, COL6A3 and SPARC genes interaction diagram圖4 PPI網(wǎng)絡(luò)樞紐基因互作圖Fig.4 PPI network hub genes interaction diagram

2.6 樞紐基因相關(guān)性分析

利用SPSS軟件計算104個個體的表型數(shù)據(jù)與58個樞紐基因表達量的Spearman相關(guān)性，統(tǒng)計相關(guān)性大于0.4的基因。與PPI分析獲得的網(wǎng)絡(luò)節(jié)點基因取交集，篩選到與EP相關(guān)的1A2和基因，相關(guān)性見圖5。

a.COL1A2基因與EP的相關(guān)性；b.SPARC基因與EP的相關(guān)性a. Correlation between COL1A2 gene and EP; b. Correlation between SPARC gene and EP圖5 候選基因相關(guān)性分析圖Fig.5 Candidate gene correlation analysis diagram

3 討 論

胴體性狀是極為復(fù)雜的經(jīng)濟性狀，是優(yōu)質(zhì)肉雞產(chǎn)肉性能選育的核心指標(biāo)。屠體率和全凈膛率是衡量產(chǎn)肉性能的主要指標(biāo)，屠體率在80.0%以上及全凈膛率在60.0%以上，被認(rèn)為是產(chǎn)肉性能良好的標(biāo)志。黃得純等曾報道，138和168日齡清遠麻雞的屠體率和全凈膛率分別為92.1%、64.8%及91.7%、65.6%。本研究結(jié)果表明，126日齡天農(nóng)麻雞的屠體率和全凈膛率的均值分別達到89.0%和67.7%，產(chǎn)肉性能良好。

基于基因表達量與表型數(shù)據(jù)的WGCNA分析可以充分利用表型信息，將基因和表型之間的關(guān)聯(lián)轉(zhuǎn)化為多個基因集和表型之間的關(guān)聯(lián)，可以有效地反映基因之間的相互作用。本研究利用104個樣本的RNA-Seq數(shù)據(jù)，通過WGCNA分析挖掘胴體性狀相關(guān)候選基因，有利于鑒定到與表型相關(guān)的基因集，相比于初步的差異基因鑒定更有利于構(gòu)建調(diào)控網(wǎng)絡(luò)。結(jié)果篩選出58個樞紐基因，涉及肌動蛋白細胞骨架的調(diào)節(jié)、ECM-受體相互作用、黏著力等通路，且與肌肉生長發(fā)育、脂質(zhì)合成和能量代謝等有關(guān)。其中ECM-受體通路中細胞外基質(zhì)參與細胞的增殖、分化和遷移，調(diào)節(jié)雞的生長發(fā)育，大量研究報道，細胞外基質(zhì)在維持骨骼肌結(jié)構(gòu)和完整性及調(diào)節(jié)細胞間物質(zhì)運輸?shù)冗^程中發(fā)揮重要作用。黏著力是復(fù)雜的細胞質(zhì)膜相關(guān)的大分子組裝體，通過招募大量黏著力相關(guān)蛋白與肌動蛋白細胞骨架形成的物理連接，是信號傳導(dǎo)過程中許多信號分子組裝的平臺，與細胞的生長、發(fā)育密切相關(guān)。本研究ECM-受體相互作用和黏著力通路中檢測到1A2、3A1、6A1、6A2、6A3、4和1與表型顯著關(guān)聯(lián)。4和1是細胞外基質(zhì)糖蛋白家族的成員，是基底膜的主要非膠原成分。它們參與多種生物過程，包括細胞黏附、分化、遷移、信號傳導(dǎo)、神經(jīng)突生長和轉(zhuǎn)移等。膠原蛋白VI的形成是由6A1、6A2和6A3基因編碼的3條ɑ鏈組成的異源三聚體。研究表明，膠原蛋白VI分子組裝時產(chǎn)生的微纖維與肌肉功能的維持有關(guān)，通過各種分子相互作用在肌細胞和周圍基質(zhì)之間提供重要的聯(lián)系。本研究結(jié)果與Bordini等的研究報道一致，均鑒定到ECM-受體相互作用通路和4A1、4A2、2、4等樞紐基因，可能通過細胞外基質(zhì)組成的改變以某種方式激活級聯(lián)生物反應(yīng)，影響肌肉生長發(fā)育，并參與改變IV型膠原蛋白的結(jié)構(gòu)。此外，樞紐基因2可通過脂肪細胞因子信號通路和FoxO信號通路發(fā)揮調(diào)節(jié)作用。該基因通過部分磷酸化轉(zhuǎn)錄因子起到調(diào)節(jié)基因表達，同時參與調(diào)節(jié)脂肪酸和膽固醇的從頭合成。2基因?qū)η卮ㄅ２糠煮w尺及肉質(zhì)性狀有顯著影響，是秦川牛新品種早期選擇的候選基因。宋嬌等以愛撥益加雞和北京油雞為試驗素材，發(fā)現(xiàn)2基因表達量與肌內(nèi)脂肪沉積顯著相關(guān)。因此，推測2基因是屠體率性狀的候選基因。

本研究將以上鑒定到的天農(nóng)麻雞胴體性狀相關(guān)的58個樞紐基因進行PPI網(wǎng)絡(luò)分析，鑒定到5個與全凈膛率相關(guān)的重要候選基因(1A2、6A1、6A2、6A3和)，同時結(jié)合Spearman相關(guān)性分析結(jié)果，確定1A2和基因表達量和全凈膛率存在極顯著的正相關(guān)。在前人的報道中，是一種Ca結(jié)合糖蛋白，與形態(tài)發(fā)生、重塑、細胞遷移和增殖有關(guān)?；蚴悄z原蛋白的鈣化、細胞外基質(zhì)的合成和軟骨內(nèi)成骨過程中細胞形狀變化所必需的。在鴨胸骨骨化中同樣發(fā)揮重要作用。1A2是膠原蛋白家族的成員，膠原蛋白是動物結(jié)締組織中的一種重要蛋白質(zhì)，也是細胞外基質(zhì)的主要成分，可以為細胞生長提供支持。且與2A1和1A2之間存在相互作用。因此，1A2和是影響全凈膛率性狀的重要候選基因。1A2和可能通過協(xié)同表達影響蛋白互作，從而影響肌肉、骨骼等多組織的生長發(fā)育。

4 結(jié) 論

天農(nóng)麻雞126日齡的屠體率和全凈膛率分別達到89.0%和67.7%，產(chǎn)肉性能良好；研究利用WGCNA和通路富集分析鑒定到ECM-受體相互作用和黏著力信號通路是影響胴體性狀的主要通路，其中1A2和是影響全凈膛率性狀的重要候選基因，為獲得天農(nóng)麻雞胴體性狀分子標(biāo)記、解析復(fù)雜性狀分子調(diào)控機制奠定重要基礎(chǔ)。