劉娟,王舒,左周,路暢,楊陽,蔡春波,趙燕,郭曉紅,曹果清,李步高,高鵬飛
(山西農(nóng)業(yè)大學 動物科學學院,山西 太谷030801)
豬肉是人類蛋白質(zhì)和脂肪的主要來源之一,約占世界肉類消費的30%。影響豬肉品質(zhì)的因素有遺傳、營養(yǎng)、應激等,評定肉質(zhì)的感官指標有肉色、pH值、嫩度、系水力、風味等[1]。肌纖維組織學特性包括肌纖維直徑、肌纖維密度、肌纖維類型、肌節(jié)長度、結(jié)締組織和肌內(nèi)脂肪含量,是評價肉品質(zhì)的重要經(jīng)濟性狀。肌纖維直徑(Diameter of muscle fiber,DiaMF)與肌肉嫩度相關,直徑越小,肉質(zhì)越嫩[2,3]。肌纖維密度(Density of muscle fiber,DenMF)和肌內(nèi)脂肪(Intramuscular fat,IMF)含量影響豬肉風味、嫩度和多汁性[4]。研究表明,肌內(nèi)脂肪含量低于2.5%的肉類會影響適口性,較高的肌內(nèi)脂肪含量通常被認為對肉質(zhì)特征具有積極的影響[5,6]。馬身豬是一個中國北方豬品種,先前研究分析了馬身豬和大白豬之間與生長速度和肉質(zhì)相關的表型差異。馬身豬比大白豬具有較低的生長速度和更好的肉質(zhì),馬身豬的肌內(nèi)脂肪含量高于大白豬[7]。進行馬身豬和大白豬的轉(zhuǎn)錄組分析,闡明2個豬種不同生長發(fā)育階段的分子機制,發(fā)現(xiàn)馬身豬和大白豬的主要區(qū)別在于氨基酸代謝,脂肪酸代謝和骨骼肌發(fā)育[8]。
高通量測序技術能夠準確獲得生物體或組織中特定時期的基因表達特征,近年來,國內(nèi)外研究者利用RNA-Seq技術從轉(zhuǎn)錄組水平探析了不同品種肌肉發(fā)育[9]、肌內(nèi)脂肪沉積[10]的基因表達模式,構建了相關基因調(diào)控網(wǎng)絡,發(fā)現(xiàn)并鑒定了PPARGC1B、TRIM 63、CSRP3、ACTN2、MYL 1、MYH 6等一系列基因是影響肌肉生長和脂肪沉積的關鍵基因[11],胰島素、絲裂原活化蛋白激酶(MAPK)、Wnt等信號通路能夠顯著影響肌纖維類型和肌肉發(fā)育[12]。WGCNA是一種應用于基因組學、蛋白質(zhì)組學和代謝組學中描述相應基因、蛋白或代謝物表達相關性的系統(tǒng)生物學分析方法[13]。Matthew等[14]首次將WGCNA應用于代謝組學數(shù)據(jù),基于WGCNA分析番茄代謝組學數(shù)據(jù),確定了與番茄成熟和遺傳背景相關性狀的3個主要代謝模塊。Priscila等[15]利用WGCNA揭示了miRNA在內(nèi)洛爾牛脂肪酸組成中的潛在調(diào)控作用;宋萬勤等[16]使用從歐洲生物信息研究所獲得的RNA-Seq數(shù)據(jù)鑒定與間充質(zhì)干細胞成骨分化有關的長非編碼RNA,發(fā)現(xiàn)32種新穎的LncRNA是 包 括miR-689,miR-640,miR-601和miR-544等在內(nèi)的miRNA前體。
隨機森林(Random forest,RF)[17]是一種機器學習領域的集成方法,它利用隨機重采樣技術和節(jié)點隨機分裂技術構建多棵決策樹,通過投票得到最終分類結(jié)果。程淑萍等[19]基于ncRNA和蛋白質(zhì)的序列信息提取特征,訓練包括隨機森林算法在內(nèi)的3個機器學習模型,預測ncRNA與蛋白質(zhì)的相互作用,預測的準確率較高,證明機器學習方法能夠預測ncRNA與蛋白質(zhì)是否存在相互作用。
課題組前期以馬身豬和大白豬為研究對象,對不同生長階段肌肉組織進行轉(zhuǎn)錄組測序,構建了肌肉組織中19條長非編碼RNA(Long noncoding RNA,LncRNA)539條mRNA組成的分子調(diào)控網(wǎng)絡[20]。本研究采用WGCNA分析技術并結(jié)合隨機森林算法,對2個品種肌肉轉(zhuǎn)錄組測序數(shù)據(jù)深入分析,以期獲得影響豬肌肉發(fā)育及肌內(nèi)脂肪沉積的關鍵基因。
所有動物手術都嚴格按照《世界醫(yī)學協(xié)會道德守則》(http://ec.europa.eu/environment/chemicals/lab_animals/legislation_en.htm)進行,實驗方案經(jīng)山西農(nóng)業(yè)大學動物倫理委員會批準。本研究選用馬身豬和大白豬為研究對象,飼養(yǎng)于山西省大同市種豬場,按照豬飼養(yǎng)標準(NY/T 65-2004)進行飼喂。分別在1、30、60、90、120、150、180 d隨機選擇3頭體重相近的去勢公豬屠宰,采集位于最后肋骨處的背最長肌樣品,液氮保存?zhèn)溆谩0凑毡匙铋L肌纖維方向垂直切割,采集1 cm3肌肉組織放入4%多聚甲醛固定。
肌肉組織在4%甲醛中固定,經(jīng)乙醇脫水、二甲苯透明后用石蠟包埋,6μm切片后采用蘇木精—伊紅(H&E)染色。每個采樣點平行樣品不少于9張切片,每張切片取3~5張圖像,采用馬尾松三色染色測定肌纖維的直徑和密度。在10×10倍顯微鏡下,對每個樣品作20個視野大方格內(nèi)肌纖維根數(shù)的測定(大方格面積約為0.25 mm2),換算后得到每平方毫米纖維根數(shù),得出肌纖維密度。在10×40倍顯微鏡下,采用不規(guī)則多邊形對每個樣本進行100根肌纖維橫斷面積的描畫,統(tǒng)計出所有的橫截面積,取其平均值,得到肌纖維面積。根據(jù)肌纖維面積計算,假設肌纖維呈圓柱形,由圓面積公式求得肌纖維直徑[21]。采用索氏脂肪抽提法分析肌內(nèi)脂肪含量[22]。
數(shù)據(jù)來源于本課題組前期豬肌肉轉(zhuǎn)錄組數(shù)據(jù)(NCBI登錄號為SRP068558(https://trace.ncbi.nlm.nih.gov/Traces/sra_sub/sub.cgi?login=pda)),共18個樣本,包含馬身豬和大白豬2個品種3個日齡(1、90、180 d)的測序數(shù)據(jù)。使用R包limma[23]進行差異基因分析,DEG的閾值設置為|log2(fold-change)|>0.8且P<0.05。
將18個樣本隨機分為訓練集和測試集,其中訓練集含10個樣本,測試集含有8個樣本。訓練集數(shù)據(jù)采用R-Studio 4.0.0軟件的WGCNA包[13]進行分析。首先計算所有基因的表達相關系數(shù),根據(jù)基因間相關性確定軟閾值(power),構建基因聚類樹。根據(jù)基因表達相關系數(shù)確定表達模塊,每個表達模塊最少基因數(shù)設定為30個。將每種基因表達模塊作為整體對象,計算模塊間的相關系數(shù),當模塊間相關系數(shù)高于0.85,則合并相關模塊。
基因功能注釋分析采用DAVID[24](http://david.abcc.ncifcrf.gov/)在線數(shù)據(jù)庫進行。提取GO分析和KEGG富集分析的結(jié)果,分別基于閾值P<0.05和P<0.01進行分析。
核心基因由模塊內(nèi)基因連通性決定。根據(jù)基因的表達趨勢計算皮爾遜相關系數(shù),取相關系數(shù)絕對值大于0.2的基因進行共表達分析。利用cytohubba插件[25]篩選每個模塊的核心基因,使用Cytoscape(version 3.5.1)[26]對 模 塊 構 建 互 作 網(wǎng)絡。利用R語言Random Forest模塊[17]進行隨機森林分析,計算模塊核心基因的平均基尼指數(shù)減少值(Mean Decrease in the Gini index,MDG),篩選hub基因。采用Graphpad 8.0軟件可視化作圖。
采用qPCR技術檢測hub基因在不同月齡2個品種背最長肌的發(fā)育性表達變化。在NCBI數(shù)據(jù)庫中獲得篩選得到的hub基因序列,通過Primer3plus在線軟件進行引物設計,選取18S作為檢測內(nèi)參基因,引物信息詳見表1。采用Prime ScriptTM RT reagent Kit with gDNA Eraser試劑盒(TaKaRa,大連,中國)合成cDNA,反轉(zhuǎn)錄時總RNA濃度為500 ng·μL-1。使用SYBR Green RTPCR(Taκara,Dalian,China)進行qPCR檢測。每個樣本進行3次技術重復,分別以無模板和無逆轉(zhuǎn)錄酶的體系作為陰性對照。反應程序如下:預變性階段,95℃30 s;循環(huán)階段,95℃30 s,60℃30 s,設定循環(huán)數(shù)為45;溶解階段,95℃30 s,60℃1 min,95℃30 s。
表1 關鍵基因的引物列表Table 1 The list of quantitative primers of hub genes
qPCR結(jié)果采用2-△△Ct法計算hub基因的相對表達量。運用SPSS 22.0軟件進行基因表達差異分析,顯著水平選擇0.05和0.01。采用Graphpad 8.0軟件進行基因表達量作圖。
利用R語言WGCNA包計算訓練集樣本中基因表達相關系數(shù)并聚類,表達趨勢性相似的基因合并為一組,構建基因共表達網(wǎng)絡(圖1)。訓練集樣本聚為兩大類,馬身豬和大白豬1日齡樣本聚為第一大類,90日齡和180日齡樣本聚為第二大類,每一大類下各品種分別聚一簇。樣本聚類樹下方為豬背最長肌肌纖維直徑、密度和肌內(nèi)脂肪含量的檢測熱圖,可以看出2個品種豬背最長肌肌纖維直徑隨著日齡均顯著增加,肌纖維密度逐漸降低,肌內(nèi)脂肪含量增加。同一日齡2個品種間相比,馬身豬的肌纖維密度顯著大于大白豬,肌內(nèi)脂肪含量也顯著高于大白豬。
圖1 樣本聚類圖和性狀熱圖Fig.1 Cluster diagram of samplesand heatdiagram of traits
采用平均鏈接層次聚類將表達模式相似的差異基因進行分組聚類,建立第一層模塊聚類樹狀圖(圖2A Dynamic Tree Cut),計算模塊間的相關系數(shù),建立模塊特征值聚類樹狀圖(圖2B),當模塊相關系數(shù)大于0.85,軟閾值大于5時,對模塊進行二次合并,建立第二層模塊聚類樹狀圖(圖2A Merged dynamic)。本研究中7262條差異基因首先獲得19個表達模塊,根據(jù)模塊相關性優(yōu)化后最終獲得16個模塊基因集,之后的分析按照合并后的模塊基因集進行。16個模塊平均含有450條基因,藍色(blue)模塊中含有的基因數(shù)最多(2960條),灰色60(grey60)和淺綠色(lightgreen)模塊基因最少,僅含50條。
圖2 差異表達基因基于相關系數(shù)聚類樹狀圖Fig.2 Dendrogram of all differentially expressed genes clustered based on correlation coefficient
通過計算模塊特征與表型間的相關系數(shù),獲得模塊性狀相關性熱圖(圖3)。肌內(nèi)脂肪與MEblack模塊成極顯著正相關(R=0.78,P=0.008),與MEsalmon模塊(R=-0.73,P=0.02)呈顯著負相關。肌纖維密度與MEblack模塊(R=-0.78,P=0.007)和MElightgreen模塊(R=-0.66,P=0.04)存在顯著負相關。只有MEmagenta模塊與肌纖維直徑存在極顯著正相關(R=0.9,P=0.0004)。
圖3 模塊特征基因與肌肉性狀的相關性熱圖Fig.3 Heat map of correlation between module characteristic genes and muscle traits
對MEblack、MEsalmon、MElightgreen和MEmagenta模塊基因集分別進行功能富集(圖4)。GO分析結(jié)果表明,MEblack模塊基因功能富集于生物合成、分子修飾及細胞內(nèi)各種代謝進程,MElightgreen模塊基因功能富集于核酸代謝、蛋白分子復合物組裝及細胞內(nèi)代謝進程,MEsalmon模塊基因功能富集于細胞質(zhì)、細胞器、細胞骨架、轉(zhuǎn)錄因子結(jié)合及一些細胞代謝進程,MEmagenta模塊基因顯著富集于細胞運動、細胞形態(tài)發(fā)生等進程。
圖4 GO和KEGG富集圖Fig.4 GO and KEGG enrichment plots
MEblack和MElightgreen模塊與肌纖維密度存在極顯著的負相關,這2個模塊共同富集于代謝和大分子修飾等進程,主要包含總代謝進程(GO:0008152),細胞代謝(GO:0044237)、初級 代謝(GO:0044238)等進程。MEblack模塊與肌內(nèi)脂肪含量存在極顯著的正相關,其中碳水化合物代謝(GO:0005975)、小分子代謝(GO:0044281)、氧化還原活性(GO:0016491)等生物進程可能是顯著影響肌內(nèi)脂肪沉積的關鍵分子進程。MEmagenta模塊與肌纖維直徑存在顯著的正相關,主要集中在 細 胞 骨 架(GO:0007010)、細 胞 質(zhì)(GO:0005829)和細胞器(GO:0006996)等生物進程。有趣的是MEsalmon模塊沒有與其它3個模塊共有的生物進程,該模塊與肌內(nèi)脂肪含量存在極顯著負相關,主要集中于激酶活性(GO:0016301)、離子結(jié)合(GO:0043167)、細胞及胞內(nèi)組成結(jié)構形態(tài)發(fā)生(GO:0000902)等生物進程,說明這些進程對肌內(nèi)脂肪沉積存在負調(diào)控作用。
對于KEGG通路富集分析發(fā)現(xiàn),MEblack模塊 顯 著 富 集 于cGMP-PKG(ko04022)、MAPK(ko04010)、基 質(zhì) 粘 附(ko04510)和 代 謝 途 徑(ko01100)等信號通路。lightgreen模塊基因顯著富集于配體識別(ko04371)、C型凝集素受體(ko04625)、cAMP(ko04024)以 及 生 物 合 成(ko04722)等信號通路。MEmagenta模塊基因極顯 著 富 集 于cGMP-PKG(ko04022)、Wnt(ko04310)、Rap1(ko04015)、cAMP(ko04024)、配體識別(ko04371)、甲狀腺激素(ko04919)、C型凝集素受體(ko04625)、液體剪切應力和動脈粥樣硬化(ko05418)等信號通路。salmon模塊基因顯著富集于液體剪切應力和動脈粥樣硬化(ko05418)、cAMP(ko04024)、Wnt(ko04310)等信號通路。
根據(jù)基因聚類結(jié)果,MEblack、MElightgreen、MEsalmon和MEmagenta4個模塊分別包含156條、50條、83條和113條差異基因,利用Cytoscape分別構建4個模塊基因互作網(wǎng)絡。圖5A為MEblack模塊共表達網(wǎng)絡,156個基因共形成11800條連接,根據(jù)cytohubba程序獲得各基因間的連通度,最 終 確 定CENPI、MAP3K 3、XRCC6、FAM 131B、RIC8B、PKP2和PLEKHA 5等基因為該模塊的hub基因。圖5B、圖5C和圖5D分別是MElightgreen、MEsalmon和MEmagenta模塊基因形成的表達網(wǎng)絡。與MEblack模塊類似,基于基因連通度結(jié)果,最終確定SNRNP48和HK3為MElightgreen模塊的hub基因,F(xiàn)AM 76B、TNXB、AGTRAP和CSDC2為MEsalmon模塊hub基因,AKT 3、CDK 16、RAB21、CBX7、ZNF792和PARP6為MEmagenta模塊hub基因。
圖5 模塊基因的共表達網(wǎng)絡Fig.5 The coexpression network of module genes
根據(jù)WGCNA分析結(jié)果,共獲得19個核心基因。將這些核心基因作為樣本特征屬性,建立隨機森林分類模型對訓練集樣本進行深度學習,通過測試集對模型進行檢驗,從而確定顯著影響肌肉肌纖維直徑、肌纖維密度和肌內(nèi)脂肪沉積3個性狀的關鍵基因。檢驗結(jié)果見表2,訓練集中共有10個樣本,只有1 d樣本出現(xiàn)預測錯誤,90 d和180 d的樣本全部預測準確,預測準確率為90%;測試集預測結(jié)果與訓練集相似,1 d樣品出現(xiàn)預測錯誤,其余樣本均準確預測,預測準確率為87.5%。
表2 訓練集和測試集的分類結(jié)果Table 2 Classification results of training set and test set
本研究中,隨機森林的預測準確率達到87.5%?;?qū)颖痉诸愵A測的貢獻率,可根據(jù)基因的MDG篩選(圖6)。CBX7(MDG=1.58)、PLEKHA 5(MDG=0.83)、RIC8B(MDG=0.53)、FAM 131B(MDG=0.52)和PARP6(MDG=0.29)5個基因?qū)颖痉诸愗暙I最大,最終確定這5個基因為影響肌纖維發(fā)育和肌內(nèi)脂肪沉積的hub基因。
圖6 基因MDG直方圖Fig.6 Histogram plot of MDG of Gene
為了驗證hub基因在豬生長發(fā)育過程中的變化趨勢,以馬身豬和大白豬為對象,選擇1、30、60、90、120、150、180 d共7個時間點進行采樣,采用qPCR檢測hub基表達(圖7)。大白豬中,CBX7、PARP6、FAM 131B和RIC8B四個基因均隨著日齡表達量逐漸升高。與大白豬相比,CBX7和PARP6在馬身豬中基因表達趨勢一致,隨著日齡表達量逐漸升高;FAM 131B和RIC8B基因出現(xiàn)先升高后降低的趨勢,在90 d時表達量最高。PLEKHA 5基因表達趨勢較特殊,在2個品種中均為先降低后升高的變化,90 d時表達量降到最低。
圖7 2個豬品種中5個關鍵基因的時序性表達Fig.7 Time-sequence expression line diagram of 5 hub genes in 2 pig breeds
肌肉中肌纖維的直徑、密度和肌內(nèi)脂肪含量是影響豬肉風味、嫩度、多汁性等肉質(zhì)性狀的關鍵因素[27]。我國地方品種和大白豬、長白豬等國外引進豬品種相比,肌纖維直徑小,肌內(nèi)脂肪含量高,系水力強,鮮嫩多汁,肉味香濃[27,28]。近年來隨著高通量測序技術的快速發(fā)展,國內(nèi)外學者在轉(zhuǎn)錄組水平對影響肌肉發(fā)育和脂肪沉積的關鍵基因和分子調(diào)控機制進行了大量研究,獲得了一系列關鍵基因和信號通路。
Sollero等[29]對皮奧伊州豬和長大二元豬胎兒肌肉組織轉(zhuǎn)錄本進行分析,發(fā)現(xiàn)了能夠顯著影響胎兒期肌肉發(fā)育的1300條差異表達的mRNA。Hu等[30]對沙子嶺豬和大白豬背最長肌進行轉(zhuǎn)錄組和蛋白組分析,獲得了488條顯著差異的基因轉(zhuǎn)錄本和23個差異蛋白,這些基因?qū)ωi肌肉脂肪酸代謝、糖酵解和肌肉發(fā)育等具有顯著的影響。最近課題組利用基因表達模式分析揭示了影響大白豬發(fā)育的基因(GHSR,EZR,F(xiàn)OXS1,DRD2,SH 3PXD2B,CSF1和TSHR),并 利 用GO和KEGG尋找了與脂肪酸合成相關的差異表達基因(ACACA,ACSF3,OXSM,CBR4和HSD17B8),發(fā)現(xiàn)肌內(nèi)脂肪在豬發(fā)育的90~180 d沉積[31]。本研究將WGCNA網(wǎng)絡分析和隨機森林深度學習相結(jié)合,從大數(shù)據(jù)多維角度探尋影響肌肉性狀的關鍵基因,最終篩選到5條hub基因和一些關鍵調(diào)控通路。李鑫等[32]基于網(wǎng)絡分析和隨機森林方法對肝癌分期進行研究,挑選出BUB1B、CCNA 2、CCNB1、CCNB2、CDC20、MAD2L1、MCM 4、PC-NA、RFC4和TOP2A等10條對肝細胞癌有重要影響的關鍵基因。有研究者利用WGCNA挖掘牛皮癬的關鍵控制因素和驅(qū)動因素,使用候選基因構建基于隨機森林的二元分類器,成功將驗證組中的正常樣本與牛皮癬樣本區(qū)分開,準確度為0.95~1,篩選獲得的hub基因作為潛在的牛皮癬病發(fā)進程的的治療靶標[33]。
本研究發(fā)現(xiàn)CBX7和PARP6基因歸類于MEmagenta模塊中,主要作用細胞運動、細胞骨架形成、細胞形態(tài)發(fā)生等功能,對肌纖維直徑有促進作用。研究發(fā)現(xiàn)CBX7通過基因重組蛋白(Dickkopf-1,DKK-1)介導的Wnt/β-catenin途徑能夠抑制乳腺腫瘤[34],膠質(zhì)瘤[35]等腫瘤細胞的增殖。結(jié)合網(wǎng)絡分析,CBX7在MEmagenta模塊中,MEmagenta模塊基因和脂肪沉積性狀高度相關,隨機森林分析發(fā)現(xiàn),CBX7對樣本分類貢獻最大。本研究中發(fā)現(xiàn)CBX7基因?qū)?nèi)脂肪沉積沒有產(chǎn)生顯著影響,但有研究發(fā)現(xiàn),在CBX7基因敲除小鼠的脂肪組織含量顯著高于野生型小鼠[36],說明了CBX7對脂肪沉積也具有重要作用,具體作用機制有待深入研究。因此可以推斷,CBX7基因是影響脂肪沉積性狀和肌肉發(fā)育的關鍵基因。研究表明PARP6是MAPK的信號靶點,影響肌肉發(fā)育和脂肪沉積[37],對大白豬皮下和肌內(nèi)脂肪組織進行測序發(fā)現(xiàn)MAPK信號通路中PARP6顯著減少,說明該基因在脂肪細胞分化過程中具有重要的調(diào)節(jié)作用[22]。PARP6能夠與胰島素反應性氨基肽酶(insulin-responsive amino peptidase,IRAP)和GRB14結(jié)合,調(diào)節(jié)脂肪細胞中葡萄糖轉(zhuǎn)運體4(glucose transporter 4,GLUT 4)的活性[37],猜測PARP6可能是影響肌肉發(fā)育和肌內(nèi)脂肪沉積的關鍵基因。
在本研究中,F(xiàn)AM 131B、RIC8B、PLEKHA 5均屬于MEblack模塊,能夠促進肌內(nèi)脂肪沉積并同時抑制肌纖維密度增加。FAM 131B基因參與MAPK信號通路,促進細胞增殖,機體免疫的形成。Huriye Cin等[38]發(fā)現(xiàn)FAM 131B與BRAF結(jié)合,激活MAPK信號通路,導致腫瘤發(fā)生。RIC8B基因編碼GEF,屬于G蛋白介導的信號傳導途徑,在細胞分裂過程中發(fā)揮重要作用[39]。研究發(fā)現(xiàn)在成骨細胞增殖過程中C/EBPβ調(diào)控RIC8B基因的高表達,而在成骨細胞分化過程中RIC8B基因表達量顯著下降[40]。PLEKHA 5在在脂質(zhì)和脂蛋白代謝通路中起到重要作用,對腫瘤細胞的增殖具有重要作用也有相關報道[41]。
本研究通過WGCNA分析獲得了與肌纖維發(fā)育和肌內(nèi)脂肪沉積的基因集,構建了重要模塊的共表達網(wǎng)絡,利用隨機森林算法對與性狀相關的功能基因進行深度學習與驗證,最終獲得5條hub基因,包括CBX7、PARP6、PLEKHA 5、FAM 131B和RIC8B。這些基因?qū)∪獍l(fā)育和脂肪代謝存在顯著作用,但具體的作用機制還需進一步研究。