基于WGCNA鑒定山羊金黃色葡萄球菌型乳腺炎關(guān)鍵應(yīng)答基因

王麗敏，靳麗娜，劉亞芳，閆曉茹，鄭惠玲

(西北農(nóng)林科技大學(xué) 動(dòng)物科技學(xué)院，陜西 楊凌 712100)

乳腺炎是山羊的常見疾病，具有發(fā)病率高、治愈率低、高產(chǎn)奶山羊泌乳盛期常見的特點(diǎn)。該病導(dǎo)致母羊泌乳量和乳品質(zhì)下降，使母羊乳腺機(jī)能減退、使用年限縮短，養(yǎng)殖生產(chǎn)經(jīng)濟(jì)損失嚴(yán)重[1]。山羊乳腺炎多由病原微生物入侵引起，金黃色葡萄球菌(Staphylococcusaureus)是主要菌種之一[2]，其所致乳房炎常為慢性，有時(shí)也可是急性[3]。乳腺組織由乳腺上皮細(xì)胞、基底層結(jié)締組織、乳腺腺泡腔、中性粒細(xì)胞、巨噬細(xì)胞、淋巴細(xì)胞等多種細(xì)胞組成。當(dāng)乳腺組織受到金黃色葡萄球菌感染時(shí)，機(jī)體的免疫系統(tǒng)使炎癥組織中白細(xì)胞增多[4]，同時(shí)損傷脫落的上皮細(xì)胞進(jìn)入乳中[5]，最終使乳中體細(xì)胞數(shù)增加。

隨著轉(zhuǎn)錄組測序技術(shù)的發(fā)展，以基因芯片為代表的基因表達(dá)數(shù)據(jù)已成廣大研究者的研究熱點(diǎn)之一，而利用加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(weighted gene co-expression network analysis, WGCNA)等生物信息學(xué)方法分析基因芯片測序數(shù)據(jù)，可在短時(shí)間內(nèi)篩選出與疾病相關(guān)的關(guān)鍵基因，對(duì)解析疾病發(fā)生發(fā)展的內(nèi)在機(jī)理發(fā)揮著重要作用[6-7]。因此，本研究采用WGCNA對(duì)乳腺感染金黃色葡萄球菌24 h山羊的乳汁中體細(xì)胞的轉(zhuǎn)錄組數(shù)據(jù)集綜合分析，并結(jié)合基因本體(gene ontology, GO)富集、蛋白互作分析(protein-protein interaction, PPI)等方法，鑒定金黃色葡萄球菌感染山羊乳腺后的關(guān)鍵應(yīng)答基因，為揭示山羊金黃色葡萄球菌型乳腺炎發(fā)生發(fā)展的分子機(jī)制，為尋找新的山羊臨床乳腺炎的防治方法提供新的理論參考。

1 材料與方法

1.1 數(shù)據(jù)的獲取及預(yù)處理

本研究所使用的原始數(shù)據(jù)來源于NCBI的GEO數(shù)據(jù)庫中GSE33894[8]數(shù)據(jù)集(https://www.ncbi.nlm.nih.gov/geo)。該數(shù)據(jù)集是以金黃色葡萄球菌侵染10只健康法國高山山羊的乳房，分別于感染前及感染后24 h采集乳汁體細(xì)胞樣本的微陣列轉(zhuǎn)錄組測序數(shù)據(jù)。基于平臺(tái)文件GPL14856對(duì)原始數(shù)據(jù)集進(jìn)行注釋，將探針名轉(zhuǎn)換為基因名，對(duì)于2個(gè)探針對(duì)應(yīng)同1個(gè)基因名的情況取平均值，對(duì)于缺失值采用K臨近(k-nearest neighbor，K-NN)算法進(jìn)行補(bǔ)充。

1.2 共表達(dá)模塊的構(gòu)建

使用R語言WGCNA包[9]的hclust函數(shù)對(duì)樣本進(jìn)行聚類，剔除離群樣本。使用pick-Soft Threshold函數(shù)計(jì)算，以無標(biāo)度拓?fù)鋽M合指數(shù)(R2)大于0.8且網(wǎng)絡(luò)具有較好的平均連通性為條件選擇適當(dāng)?shù)能涢撝?。使用adjacency函數(shù)將原始矩陣轉(zhuǎn)換為拓?fù)渲丿B矩陣(TOM)，依據(jù)TOM矩陣相異度(dissTOM)，使用hclust函數(shù)對(duì)基因進(jìn)行聚類，得到多個(gè)基因共表達(dá)模塊。采用主成分分析法計(jì)算每個(gè)模塊的特征值(module eigengene，ME)，根據(jù)ME對(duì)模塊聚類，對(duì)相似度大于0.75的共表達(dá)模塊進(jìn)行合并。

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

將合并后的模塊與表型信息矩陣進(jìn)行關(guān)聯(lián)分析，選擇與感染階段相關(guān)的模塊為目標(biāo)模塊(依據(jù)相關(guān)系數(shù)r及P值進(jìn)行選擇，P<0.05為差異顯著)。計(jì)算目標(biāo)模塊內(nèi)每個(gè)基因的模塊身份(moudle membership，MM)和基因顯著性(gene significance，GS)的關(guān)系，驗(yàn)證所選模塊的可靠性。選擇目標(biāo)模塊內(nèi)基因連接度排名前30的基因?yàn)闃屑~基因(hubgene)。通過DAVID[10]網(wǎng)站對(duì)目標(biāo)模塊基因進(jìn)行GO功能富集分析。

1.4 關(guān)鍵應(yīng)答基因的選擇

通過String在線網(wǎng)站[11]對(duì)目標(biāo)模塊基因構(gòu)建PPI網(wǎng)絡(luò)，利用Cytoscape[12]軟件的CytoHubba插件計(jì)算PPI網(wǎng)絡(luò)中節(jié)點(diǎn)的度值(degree)，選擇排名前30的基因?yàn)楹诵幕?core gene)，并與模塊內(nèi)的樞紐基因取交集得到關(guān)鍵應(yīng)答基因。

2 結(jié)果與分析

2.1 共表達(dá)模塊的構(gòu)建

數(shù)據(jù)預(yù)處理后得到8 032個(gè)基因用于后續(xù)分析。對(duì)20個(gè)樣本進(jìn)行聚類分析，未發(fā)現(xiàn)離群樣本。以R2>0.8且具有一定連通性為標(biāo)準(zhǔn)選擇軟閾值為5，根據(jù)ME進(jìn)行模塊聚類(圖1)，圖中紅色線的高度是0.25，在紅色線條以下的模塊是相似度較大需要合并的模塊。以相似度大于0.75合并共表達(dá)模塊，得到11個(gè)模塊(圖2)，圖中每個(gè)顏色代表1個(gè)模塊，第1行的顏色表示首次聚類的結(jié)果，第2行的顏色是模塊經(jīng)過合并后的結(jié)果。其中灰色模塊的基因?qū)儆诠脖磉_(dá)趨勢不明顯的基因，單獨(dú)形成1個(gè)模塊。

圖1 模塊特征基因的聚類

圖2 基因動(dòng)態(tài)剪切聚類樹

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

依據(jù)每個(gè)模塊的特征值(ME值)，將模塊與表型信息矩陣進(jìn)行關(guān)聯(lián)分析，得到每個(gè)模塊與不同感染階段的相關(guān)性及P值(圖3)，選擇與感染24 h相關(guān)性最強(qiáng)且P<0.05的blue模塊(r=0.69，P=7×10-4)和green模塊(r=0.68，P=0.001)為目標(biāo)模塊。

圖3 模塊與感染階段關(guān)聯(lián)分析

GO富集是對(duì)基因在不同維度和不同層次的描述，包括3部分：生物過程(biological process，BP)、細(xì)胞組分(cellular component，CC)和分子功能(molecular function，MF)。分別表示基因參與的生物學(xué)過程、基因存在的位置及基因在分子層面的功能。通過對(duì)所選目標(biāo)模塊的基因進(jìn)行GO(gene ontology)富集(每部分富集結(jié)果取基因count值前5展示，不足5個(gè)的全部展示)。

結(jié)果顯示，金黃色葡萄球菌感染24 h的blue模塊基因主要具有蛋白同源二聚化、相同蛋白識(shí)別、ATP結(jié)合等分子功能；green模塊基因主要參與細(xì)胞內(nèi)蛋白水解、泛素依賴的蛋白酶體降解過程，主要具有泛素結(jié)合、SMAD蛋白結(jié)合、分子伴侶HSP90蛋白結(jié)合等分子功能。

圖4 Blue模塊GO富集

圖5 green模塊GO富集

2.3 關(guān)鍵應(yīng)答基因的篩選

利用Cytoscape軟件構(gòu)建核心基因的PPI網(wǎng)絡(luò)見圖6、圖7，圖中紅色越深表示基因連接度越高。依據(jù)目標(biāo)模塊基因連接度篩選出的樞紐基因表示基因在轉(zhuǎn)錄水平發(fā)揮關(guān)鍵作用，依據(jù)PPI網(wǎng)絡(luò)篩選出的核心基因代表基因在蛋白水平發(fā)揮關(guān)鍵作用[13]。而二者的交集基因(圖中菱形表示的基因)即為本研究所選擇的關(guān)鍵應(yīng)答基因。因此，山羊乳腺感染金黃色葡萄球菌24 h的4個(gè)關(guān)鍵應(yīng)答基因?yàn)镹FKB1、MRPL16、TMED10和ERLEC1。

圖6 Blue模塊PPI網(wǎng)絡(luò)核心基因互作

圖7 Green模塊PPI網(wǎng)絡(luò)核心基因互作

3 討 論

金黃色葡萄球菌感染山羊乳腺組織的乳汁體細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)來源于GSE33894數(shù)據(jù)集，Cremonesi等[8]于2012年通過利用Limma函數(shù)對(duì)數(shù)據(jù)進(jìn)行差異表達(dá)分析，發(fā)現(xiàn)金黃色葡萄球菌感染山羊乳腺組織24 h后乳汁體細(xì)胞中有300個(gè)基因呈現(xiàn)差異表達(dá)，其中251個(gè)基因表達(dá)上調(diào)，49個(gè)基因表達(dá)下調(diào)。本研究以樣本數(shù)據(jù)基因表達(dá)量為基礎(chǔ)，利用WGCNA技術(shù)對(duì)數(shù)據(jù)構(gòu)建基因共表達(dá)網(wǎng)絡(luò)，使表達(dá)趨勢相同的基因聚為同一模塊，表達(dá)趨勢不同的基因聚為不同模塊，將模塊與表型(感染時(shí)間)相關(guān)聯(lián)后，結(jié)合PPI方法，依次獲得模塊核心基因和模塊關(guān)鍵基因，最終交集獲得山羊乳腺感染金黃色葡萄球菌24 h的4個(gè)關(guān)鍵應(yīng)答基因(NFKB1、MRPL16、TMED10和ERLEC1)。

NF-κB是真核生物細(xì)胞調(diào)節(jié)炎癥和對(duì)外來刺激免疫應(yīng)答的重要響應(yīng)蛋白[14]，對(duì)維持免疫系統(tǒng)的穩(wěn)態(tài)十分重要。NF-κB1和RelA是NF-κB家族的部分成員[15]。IKBKE是NF-κB激活因子，其可通過解除IκB對(duì)NF-κB的抑制，使NF-κB入核并激活下游靶基因的轉(zhuǎn)錄[16]。王晉鵬等[17]研究表明乳腺炎發(fā)生后，細(xì)胞通過激活NF-κB信號(hào)通路啟動(dòng)機(jī)體的先天免疫防御。NF-κB1是NF-κB信號(hào)通路的關(guān)鍵節(jié)點(diǎn)，是炎癥反應(yīng)和免疫反應(yīng)的熱門靶點(diǎn)[18]。李蕊等[19]研究證實(shí)NF-κB1是奶牛金黃色葡萄球菌型乳腺炎的關(guān)鍵抗性候選蛋白。楊星哲等[20]研究發(fā)現(xiàn)NF-κB1是丙酸桿菌誘導(dǎo)的兔耳痤瘡炎癥模型的關(guān)鍵應(yīng)答因子。他們的研究結(jié)果與本研究結(jié)果相符。提示雖然感染菌種屬和受菌物種存在差異，但NF-κB1在多種炎癥免疫模型中具有重要作用。

蛋白質(zhì)在細(xì)胞生長發(fā)育、增殖分化、物質(zhì)代謝等方面有極其復(fù)雜的調(diào)控過程。本研究中g(shù)reen和blue 2個(gè)目標(biāo)模塊是與金黃色葡萄球菌感染24 h最為相關(guān)的2個(gè)模塊，GO富集分子功能中，大多數(shù)模塊基因聚焦于蛋白同源二聚化活性、同質(zhì)蛋白質(zhì)結(jié)合、ATP結(jié)合、泛素結(jié)合、Hsp90蛋白結(jié)合、SMAD結(jié)合等分子功能。這正說明金黃色葡萄球菌侵染乳腺后，細(xì)胞調(diào)動(dòng)并激活了多種蛋白信號(hào)傳導(dǎo)途徑以啟動(dòng)免疫防御。蛋白質(zhì)泛素化通過將泛素分子遞續(xù)連接到目標(biāo)蛋白上，對(duì)蛋白進(jìn)行降解，調(diào)節(jié)細(xì)胞增殖、凋亡和免疫信號(hào)傳導(dǎo)等過程。Ebstein等[21]通過炎癥誘導(dǎo)物刺激樹突狀細(xì)胞后，許多泛素相關(guān)基因表達(dá)上升。而高慧杰[22]的研究中，以生物信息學(xué)挖掘到泛素化相關(guān)的基因(UBE2L6、UBA7、USP18、ISG15)，是奶牛金黃色葡萄球菌乳房炎的關(guān)鍵應(yīng)答基因。本研究中一些泛素相關(guān)基因(USP5、UBE2M、UBE2D1、UBE2S)，在金黃色葡萄球菌感染山羊乳腺24 h的green模塊中也被鑒定到。同時(shí)green模塊GO富集到的生物過程集中于細(xì)胞蛋白水解過程、泛素依賴的蛋白酶解過程。表明細(xì)胞中蛋白質(zhì)降解系統(tǒng)的激活與重塑在炎癥反應(yīng)中發(fā)揮了重要作用。

鈣結(jié)合蛋白S100A8在炎癥相關(guān)研究中已廣見報(bào)道，通常被認(rèn)為是啟動(dòng)炎癥通路、促進(jìn)炎癥發(fā)展進(jìn)程的關(guān)鍵炎癥介質(zhì)[23]。在Cremonesi等[8]的研究中，S100A8被定義為金黃色葡萄球菌感染山羊乳腺后的抵抗因子和促炎介質(zhì)。同時(shí)其也是本研究中感染金黃色葡萄球菌24 h最相關(guān)的blue模塊的樞紐基因，但該基因并未出現(xiàn)在PPI網(wǎng)絡(luò)中，可能是因?yàn)樵谘装Y條件下S100A8mRNA表達(dá)和蛋白表達(dá)水平存在差異。關(guān)鍵基因MRP16編碼線粒體核糖體蛋白，在惡性腫瘤中通常高表達(dá)，具有促進(jìn)癌細(xì)胞增殖、正向調(diào)控惡性腫瘤的功能[24]。但MRP16在乳腺炎相關(guān)研究中鮮有報(bào)道。關(guān)鍵基因TMED10編碼跨膜轉(zhuǎn)運(yùn)蛋白，通過內(nèi)質(zhì)網(wǎng)-高爾基體系統(tǒng)參與蛋白質(zhì)的囊泡運(yùn)輸，在早期胚胎發(fā)育及免疫反應(yīng)中發(fā)揮重要作用[25]。關(guān)鍵基因ERLEC1編碼常駐內(nèi)質(zhì)網(wǎng)蛋白，具有識(shí)別n-多聚糖的功能，在乳腺炎相關(guān)研究尚未見報(bào)道。

乳腺炎的發(fā)生發(fā)展涉及多通路的調(diào)控網(wǎng)絡(luò)，其內(nèi)在機(jī)理繁多復(fù)雜。目前，關(guān)于乳腺炎機(jī)制的研究報(bào)道較多，但本次鑒定出的一些關(guān)鍵基因在乳腺炎領(lǐng)域鮮有報(bào)道，其在炎癥發(fā)展進(jìn)程的分子作用機(jī)制尚待探究。

4 結(jié) 論

本研究通過利用WGCNA、GO富集、PPI方法對(duì)樣本數(shù)據(jù)進(jìn)行綜合分析發(fā)現(xiàn)，金黃色葡萄球菌感染山羊乳腺24 h，多種蛋白參與細(xì)胞信號(hào)轉(zhuǎn)導(dǎo)(包括蛋白泛素化、分子伴侶蛋白HSP90結(jié)合、受體信號(hào)轉(zhuǎn)導(dǎo)蛋白Smad結(jié)合、蛋白二聚化等)；4個(gè)關(guān)鍵基因NFKB1、MRPL16、TMED10、ERLEC1在炎癥發(fā)展進(jìn)程中發(fā)揮核心調(diào)控作用。