韓信嶸,吳慧光,2,吳江鴻,王國富,3,高樹新,3
(1.內(nèi)蒙古民族大學(xué) 動物科技學(xué)院,內(nèi)蒙古 通遼 028042; 2.揚州大學(xué) 獸醫(yī)學(xué)院,江蘇 揚州 225009;3.內(nèi)蒙古民族大學(xué) 黃牛遺傳繁育研究所,內(nèi)蒙古 通遼 028000)
中國是世界上第三大牛肉生產(chǎn)國和消費國,隨著人們的生活水平不斷提高,中國居民牛肉消費持續(xù)增加。牛肉品質(zhì)的優(yōu)劣直接影響到適口性,為了迎合大眾口味,滿足大眾需求,提高牛肉品質(zhì)是當(dāng)前育種中亟待解決的問題,而研究脂肪沉積相關(guān)基因的功能和調(diào)控機制對改善肉質(zhì)至關(guān)重要。安格斯牛和中國西門塔爾牛是具有不同脂肪沉積能力的2個品種,安格斯牛的背膘厚度顯著高于中國西門塔爾牛[1],剪切力顯著低于中國西門塔爾牛[2],肌內(nèi)脂肪(Intramuscular fat, IMF)含量顯著高于中國西門塔爾牛[3]。中國西門塔爾牛于20世紀(jì)80年代引入中國,是乳、肉、役兼用的全能牛,經(jīng)育種改良的中國西門塔爾牛是科爾沁主要肉牛品種,在我國東北部、南部以及西部地區(qū)均有分布。中國西門塔爾牛肌肉豐滿,生長速度較快,胴體肉多,脂肪少而分布均勻,屠宰率可達(dá)65%,相對于黃牛、荷斯坦牛和蒙古牛來說,西門塔爾牛在屠宰率、瘦肉率及肉用性能方面具有顯著的優(yōu)勢,但其肉質(zhì)要稍遜色于世界上其他著名的肉牛品種[4]。安格斯牛是20世紀(jì)70年代引入中國,經(jīng)長期育種后,分布在新疆、內(nèi)蒙古、東北、山東等地,被毛的顏色主要為黑色,其次為黑棕色或紅色,身材矮小,耐寒能力、耐粗飼和抗病能力都很強。安格斯牛早熟、肉質(zhì)質(zhì)量好、凈肉率高、大理石花紋明顯,屠宰率達(dá)60%~65%,許多國家及地區(qū)常用安格斯牛作為母本對地方牛進(jìn)行雜交來改良品種[5]。
脂肪沉積能力是牛肉品質(zhì)的評價指標(biāo)之一,影響脂肪沉積能力的因素有很多,如遺傳因素、生長環(huán)境和營養(yǎng)水平等,其中遺傳因素占據(jù)主導(dǎo)地位。目前發(fā)現(xiàn)參與脂肪代謝的分子主要有脂肪酸結(jié)合蛋白、脂肪酸合成酶、黑色素皮質(zhì)素受體和脂聯(lián)素等。另外,許多激素和細(xì)胞因子也參與脂質(zhì)代謝,如胰島素、糖皮質(zhì)激素和成脂因子等,上述分子基因的多態(tài)性與脂肪沉積的相關(guān)性已經(jīng)通過試驗得到了證實,但目前關(guān)于安格斯牛和中國西門塔爾牛脂肪沉積的關(guān)鍵基因及調(diào)節(jié)機制的研究報道卻很少。
為了研究與脂肪沉積相關(guān)的基因,本試驗以安格斯牛和中國西門塔爾牛作為研究對象,對其背膘組織進(jìn)行轉(zhuǎn)錄組測序,通過GO富集和KEGG代謝通路的分析,結(jié)合Cytoscape等分析與脂肪沉積有關(guān)的基因,并對與脂肪沉積相關(guān)的3個基因DKKL1、ACACB和PTGS2進(jìn)行了qRT-PCR的驗證,以期找出調(diào)控脂質(zhì)代謝的候選基因,為中國西門塔爾牛和安格斯牛的后續(xù)選育提供科學(xué)依據(jù)。
安格斯牛和中國西門塔爾牛都來自內(nèi)蒙古通遼市余糧畜業(yè)育肥牛場,屠宰前牛場按照日糧標(biāo)準(zhǔn)進(jìn)行集中統(tǒng)一飼喂,且都為閹牛,育肥期為60 d。肉牛樣本是在統(tǒng)一條件下進(jìn)行采集,及時取樣,液氮保存。進(jìn)行轉(zhuǎn)錄組測序的6頭牛是通過聚類分析挑選的,能夠很好地反映總體水平。
TRIzol法提取RNA后用Nano-Drop2000測定RNA濃度。檢測合格后按照說明建庫后利用BGISEQ-500平臺對安格斯牛和中國西門塔爾牛背膘樣品的RNA測序并獲得轉(zhuǎn)錄組數(shù)據(jù)。使用SOAPnuke軟件對所得數(shù)據(jù)進(jìn)行過濾,去除低質(zhì)量reads以確保數(shù)據(jù)的準(zhǔn)確性。
與參考基因組比對之后,先使用StringTie[6]對每個樣品進(jìn)行轉(zhuǎn)錄本重構(gòu),再使用Cuffmerge(Cufflinks的工具之一)將所有樣品的重構(gòu)信息整合在一起,最后使用Cuffcompare[7]將重構(gòu)的轉(zhuǎn)錄本與參考注釋信息比較得到新轉(zhuǎn)錄本。使用CPC[8]對新轉(zhuǎn)錄本進(jìn)行蛋白編碼潛力預(yù)測,最終將預(yù)測具有蛋白編碼潛力的新轉(zhuǎn)錄本加入到參考基因序列中,得到一個完整的參考序列信息,后續(xù)將基于此參考序列進(jìn)行分析。
根據(jù)參考基因組比對結(jié)果,利用GATK[9]檢測每個樣品的SNP信息。通過比對基因組數(shù)據(jù)、標(biāo)記重復(fù)序列、分割序列和重校準(zhǔn)堿基對數(shù)據(jù)進(jìn)行預(yù)處理。然后用HaplotypeCaller檢測SNP。使用rMATS[10]檢測不同樣品之間的差異剪接基因(DSG)和樣品自身的差異剪接基因,以FDR≤0.05為判斷標(biāo)準(zhǔn),用DEGseq[11]法進(jìn)行差異基因的檢測。并將P-values矯正為qvalues。為了提高DEGs的準(zhǔn)確性,將差異倍數(shù)為2倍以上并且qvalue≤0.001的基因,篩選為顯著差異表達(dá)基因。
對轉(zhuǎn)錄組數(shù)據(jù)中的差異表達(dá)基因進(jìn)行GO(Gene Ontology)功能注釋和KEGG通路富集分析。過濾后的Clean reads保存為FASTQ格式根據(jù)差異基因檢測結(jié)果,對其GO功能進(jìn)行分類以及富集分析。
根據(jù)GO注釋結(jié)果以及官方分類,將差異基因進(jìn)行功能分類,同時使用R軟件中的phyper函數(shù)進(jìn)行富集分析。然后對P-value進(jìn)行FDR校正,通常FDR≤0.01的功能視為顯著富集。
生物體內(nèi)的生物學(xué)行為是由多種蛋白相互協(xié)調(diào)的,Pathway分析是基于公共數(shù)據(jù)庫KEGG對生物學(xué)功能的進(jìn)一步了解,可以初步了解蛋白質(zhì)的功能,還有助于驗證和闡明蛋白質(zhì)的生物調(diào)控機制[12]。KEGG提供了通過Pathway作圖過程將基因組與生命聯(lián)系起來的參考知識庫,可以將基因的基因組或轉(zhuǎn)錄組內(nèi)容物映射到KEGG參考途徑以推斷細(xì)胞或生物體的系統(tǒng)行為[13]。KEGG資源提供了將基因組與生物系統(tǒng)連接起來的參考知識庫,分類為基因組空間(KEGG GENES)和化學(xué)空間(KEGG LIGAND)中的構(gòu)建塊包括內(nèi)源分子和外源分子,以及受體互作和反應(yīng)網(wǎng)絡(luò)(Pathway)的接線圖以及計算機化的KeggBrite[12]。
熒光定量分析的主要目的是驗證轉(zhuǎn)錄組數(shù)據(jù)的差異表達(dá)基因,首先提取樣本總RNA并通過Nanodrop對RNA含量進(jìn)行檢測,然后將RNA逆轉(zhuǎn)錄成cDNA并對其進(jìn)行凝膠電泳檢測,配制25.0 μL的反應(yīng)體系:SYBR 12.5 μL,上下游引物各0.5 μL,cDNA模板2.0 μL,去離子水9.5 μL;配好后將樣品加到96孔板中,然后以4 000 r/min的轉(zhuǎn)速離心30 s,最后放到已做好設(shè)定的熒光定量PCR儀中。qRT-PCR可以對差異表達(dá)基因的表達(dá)量進(jìn)行分析。
2.1.1 測序數(shù)據(jù)過濾 過濾后的reads質(zhì)量統(tǒng)計結(jié)果見表1,通過BGISEQ-500平臺對來自中國西門塔爾和安格斯牛的背膘的6個樣品進(jìn)行轉(zhuǎn)錄組測序,過濾前reads數(shù)為120.36 Mb。安格斯牛每個樣品平均產(chǎn)出10.95 Gb,中國西門塔爾牛平均產(chǎn)出10.96 Gb。過濾后reads占比在安格斯牛為90.98%,中國西門塔爾牛為91.05%。數(shù)據(jù)過濾后,Q20(質(zhì)量大于20的堿基數(shù)目占總堿基數(shù)目的比例)安格斯牛為97.49%,中國西門塔爾牛為97.56%,都大于97.00%,Q30安格斯牛為90.22%,中國西門塔爾牛為90.36%,都大于90.00%,表明測序結(jié)果良好,符合建庫要求。
表1 過濾后的reads質(zhì)量
2.1.2 參考基因組比對 reads對基因的比對率統(tǒng)計見表2,使用HISAT[14]將Clean reads與參考基因組序列比對,結(jié)果顯示,樣品比對基因組的總比對率安格斯牛為92.19%,中國西門塔爾牛為92.47%,比對基因集的唯一比對率為安格斯牛為70.56%,中國西門塔爾牛為71.47%。同時樣品間均勻的比對率表明,樣品之間的數(shù)據(jù)具有可比性。
表2 參考基因組比對結(jié)果
2.1.3 新轉(zhuǎn)錄本預(yù)測 與參考基因組比對之后,將轉(zhuǎn)錄本重構(gòu),與參考注釋信息比較得到新轉(zhuǎn)錄本。共計得到新轉(zhuǎn)錄本19 747個,其中預(yù)測的16 576個為編碼的轉(zhuǎn)錄本數(shù),15 643個屬于已知蛋白編碼基因新的可變剪接亞型,933個屬于新的蛋白編碼基因的轉(zhuǎn)錄本,剩下的3 171個屬于LncRNA。
與參考基因組比對之后使用GATK檢測每個樣品的SNP[9],SNP統(tǒng)計類型見表3。安格斯牛的A-G和C-T變異的SNP數(shù)目158 116個,A-C、A-T、C-G 和 G-T 變異的SNP數(shù)目為60 349個,所有變異類型的總數(shù)目為218 465個。中國西門塔爾牛的A-G和C-T變異的SNP數(shù)目179 321個,A-C、A-T、C-G 和 G-T 變異的SNP數(shù)目為68 589個,所有變異類型的總數(shù)目為247 910個。安格斯牛的變異的SNP數(shù)目普遍低于中國西門塔爾牛。
表3 SNP類型統(tǒng)計
根據(jù)各個樣品基因表達(dá)水平,檢測樣品之間的差異表達(dá)基因。中國西門塔爾牛對安格斯牛的差異表達(dá)基因的上下調(diào)關(guān)系見圖1,中國西門塔爾牛對安格斯牛差異基因表達(dá)的分析中有1 296個差異表達(dá)基因,其中上調(diào)基因802個(Fold change≥2,q≤0.001),下調(diào)基因494個(Fold change≤-2,q≤0.001),即中國西門塔爾牛上調(diào)表達(dá)基因數(shù)高于安格斯牛的上調(diào)表達(dá)基因數(shù)。
圖1 DEGs的Scatter-plot分布
根據(jù)差異基因檢測結(jié)果,對所得數(shù)據(jù)進(jìn)行進(jìn)一步的分類及富集分析,由圖2可看出,對篩選的6 244個基因進(jìn)行GO功能分類注釋和富集,其中2 205個DEGs注釋到26個生物學(xué)進(jìn)程(Biologicalprocess,BP)中,3 149個DEGs注釋到15個細(xì)胞組分(Celluar component, CC)中,890個DEGs注釋到12個分子功能(Molecular function, MF)中。進(jìn)一步的分類和富集分析顯示富集比較靠前且超過300個DEGs,BP中的細(xì)胞進(jìn)程398個,占比6.37%;單細(xì)胞進(jìn)程369個,占比5.91%;CC中的細(xì)胞組分396個,占比6.34%;膜組分394個,占比6.31%;MF中的結(jié)合功能380個,占比6.09%。
圖2 差異基因GO功能分類
根據(jù)差異基因檢測結(jié)果,對其進(jìn)行KEGG生物通路分類以及富集分析,結(jié)果見圖3,通過KEGG代謝通路分析牛背膘組織的Unigene參與了6大類共44小類的代謝通路,并且一共注釋到了2 158個基因。其中基因數(shù)量排名靠前的代謝通路主要包括人類疾病通路,注釋了616個基因,占比28.54%,有機系統(tǒng)注釋了586個基因,占比27.15%,然后依次為環(huán)境信息進(jìn)程、新陳代謝通路、細(xì)胞進(jìn)程以及遺傳信息進(jìn)程,其依次占比為19.00%,12.00%,9.87%,3.43%。
圖3 差異基因Pathway分類
對差異表達(dá)基因的數(shù)據(jù)進(jìn)行了Pathway富集,由圖4可知在cAMP信號通路、神經(jīng)活性配體-受體互作以及細(xì)胞因子受體互作3個通路上基因的富集程度最高,即數(shù)量最多且富集結(jié)果顯著。
圖4 差異基因Pathway富集分析
圖5為脂肪相關(guān)基因DKKL1、ACACB及PTGS2在轉(zhuǎn)錄組數(shù)據(jù)中表達(dá)量與qRT-PCR相對表達(dá)量的對比。熒光定量結(jié)果DKKL1、ACACB在中國西門塔爾牛的背膘中的相對表達(dá)量高于安格斯牛的背膘中的相對表達(dá)量(2-(ΔΔCt)>1),PTGS2在安格斯牛背膘中相對表達(dá)量高(2-(ΔΔCt)<1),差異表達(dá)模式均與RNA-seq表達(dá)模式一致,表明轉(zhuǎn)錄組數(shù)據(jù)可靠??梢奃KKL1和ACACB在中國西門塔爾牛中表達(dá)量較高,PTGS2在安格斯牛體內(nèi)表達(dá)較高。
圖5 DKKL1、ACACB及PTGS2的qRT-PCR中相對表達(dá)量(A)及RNA-seq表達(dá)量(B)
在初期屠宰試驗中得到中國西門塔爾牛剪切力顯著高于安格斯牛,且安格斯牛背膘脂肪厚度顯著高于中國西門塔爾牛(P<0.01)[1]。IMF含量與剪切力呈負(fù)相關(guān),與嫩度呈正相關(guān)[15]。2個品種牛在脂肪沉積上有顯著差異,故候選基因的篩選具有一定的意義。
本研究通過BGISE-500對被背膘組織進(jìn)行轉(zhuǎn)錄組測序,并對測序數(shù)據(jù)進(jìn)行富集和注釋,從中找到與脂肪沉積相關(guān)功能和通路。從圖 4中可以看出注釋到cAMP信號通路中的差異基因數(shù)目最多且富集結(jié)果較為顯著。同時cAMP能促進(jìn)動物機體內(nèi)脂肪的分解代謝[16],但可降低肉中硬脂酸含量,增加不飽和脂肪酸含量,增加犢牛肉中脂肪含量[17]。線粒體內(nèi)由二氧化碳/碳酸氫鹽調(diào)節(jié)的可溶性cAMP激活蛋白激酶A(PKA)使線粒體蛋白磷酸化從而調(diào)節(jié)ATP的產(chǎn)生[18]。cAMP增加會激活cAMP依賴蛋白激酶,使激素敏感脂酶(HSL)磷酸化進(jìn)而催化甘油三酯水解為甘油和脂肪酸[19]。轉(zhuǎn)變成酰基CoA,而后進(jìn)入β-氧化和三羧酸循環(huán)(TCA)氧化產(chǎn)熱,游離脂肪酸(FFA)可使線粒體進(jìn)入解偶聯(lián)狀態(tài),底物氧化耗能增強,最終導(dǎo)致脂肪分解加強[20]。
依據(jù)|log2(Fold chage)|≥1,q≤0.001對脂肪沉積表達(dá)差異基因進(jìn)行篩選,初步篩選出3種基因DKKL1、ACACB及PTGS2。DKKL1(Dickkopf like acrosomal protein 1,Soggy)阻黑體,類似頂體蛋白1抑制物,位于牛的18號染色體,有7個外顯子,長度為985 bp,編碼233個氨基酸,屬于分泌性蛋白,是Wnt抑制劑小蛋白家族Dickkopf的成員。DKKs編碼的分泌蛋白通常通過阻止配體-受體相互作用[21]或抑制Wnt輔助受體LRP5/6來拮抗Wnt/β-catenin信號轉(zhuǎn)導(dǎo),可局部地抑制Wnt調(diào)節(jié)的過程而達(dá)到促進(jìn)脂肪沉積的功能。Mao等[22]的研究認(rèn)為DKK1通過Krm蛋白的降解,破壞Wnt蛋白配體與Fz-LRP6形成的受體復(fù)合物,從而抑制Wnt信號通路。研究顯示DKK1、DKK2、DKK3及DKK4都可以抑制Wnt通路的傳導(dǎo),促進(jìn)脂肪沉積。DKKL1雖是Wnt家族的一員,但與DKK1,DKK2和DKK4不同[23],盡管與DKK3擁有相同的Sgy結(jié)構(gòu)域,但是卻沒有2個富含Cys的LRP相互作用結(jié)構(gòu)域[24],同時也有研究表明,非特異性通路中,Wnt只與Fz作用,而不需要LRP5/6。Dakhova等[24]研究顯示DKKL1可以通過負(fù)調(diào)節(jié)類固醇生成酶來調(diào)控睪酮的生成。睪酮可以降低膽固醇和體脂,但公閹牛體內(nèi)的睪酮含量低于公牛,其作用機制和調(diào)節(jié)水平還有待進(jìn)一步確定。因此,差異表達(dá)基因DKKL1對脂肪沉積的作用及在2個品種牛中的表達(dá)機制還需要后續(xù)試驗的驗證。
ACACB(Acetyl-CoA carboxylase beta,ACC2)乙酰輔酶A羧化酶β基因,位于牛的17號染色體,有56個外顯子,全長8 924 bp,共編碼2 431個氨基酸。ACC包含2種,ACC1和ACC2,其中ACC1位于胞質(zhì),主要參與脂肪的從頭合成,ACC2主要存在于線粒體外膜,主要作用于脂肪酸的氧化。它們具有70%相同的酶活性和同源性。在AMPK信號通路當(dāng)中,AMPK會抑制ACACA和ACACB的活性,而它們的作用是催化乙酰輔酶A(CoA)羧化為丙二酰輔酶A,丙二酰輔酶A通過抑制肉毒堿棕櫚?;D(zhuǎn)移酶1(CPT1)而抑制脂酰肉堿向線粒體內(nèi)的轉(zhuǎn)運從而導(dǎo)致游離脂肪酸增加[25],其中丙二酰輔酶A通常由線粒體外膜中的ACC2產(chǎn)生[26]。ACC2位于線粒體外膜,主要參與脂肪酸細(xì)胞氧化。在不同種屬中脂肪酸調(diào)節(jié)也有差異,ACC2含量也不同。Oh等[26]認(rèn)為ACC2將成為治療肥胖和相關(guān)代謝疾病(如高脂血癥、胰島素抵抗和糖尿病)的靶基因。敲除ACACB會導(dǎo)致因為ACACB產(chǎn)生的Malonyl-COA水平的降低,進(jìn)而導(dǎo)致脂肪酸降低,脂肪含量降低[27]。高麗南等[28]研究顯示,不同飼料營養(yǎng)以及不同品種肉牛中AMPK的含量不同,其中,中國西門塔爾牛的AMPK含量較高。AMPK脂肪酸氧化中起著重要作用,即將脂肪分解為游離的脂肪酸。在肌肉中脂肪酸氧化的提高與明顯降低的丙二酸單酰輔酶A水平有關(guān),ACC2調(diào)控脂肪酸氧化,在肝臟中丙二酸單酰輔酶A水平就沒有變化,這表明ACC2必須在一個沒有ACC1產(chǎn)生的丙二酸單酰輔酶A的區(qū)域才能發(fā)揮作用。ACACB的作用機制和調(diào)節(jié)水平還有待進(jìn)一步確定,差異表達(dá)基因ACACB對脂肪沉積的作用及在2個品種牛中的表達(dá)機制還需要后續(xù)試驗的驗證。
PTGS2(Prostaglandin-endoperoxide synthase 2,COX2)前列腺素內(nèi)過氧化物合酶2,也稱環(huán)氧化酶-2。位于牛的16號染色體,含有10個外顯子,全長3 489 bp,共編碼604個氨基酸。COX2是一種催化花生四烯酸轉(zhuǎn)變成前列腺素前體(PGH2)的限速酶,是線粒體編碼蛋白,也是生物活性脂類調(diào)節(jié)物[29]。在脂肪細(xì)胞分解調(diào)節(jié)通路(Regulation of lipolysis in adipocytes)中COX2會間接促進(jìn)前列腺素E2(Prostaglandin E2,PGE2)的表達(dá),PGE2在前列腺受體EP3的作用下會促進(jìn)抑制性蛋白(Gi)的表達(dá),從而抑制腺苷酸環(huán)化酶(Adenylatecyclase,AC)的表達(dá),從而導(dǎo)致cAMP信號通路激活受阻,最終促進(jìn)脂肪沉積。PTGS2的表達(dá)以體質(zhì)量指數(shù)依賴的方式增加,并且與肥胖程度密切相關(guān)。COX2在成熟脂肪細(xì)胞中的表達(dá)受到被噻唑烷二酮特異性靶向的PPARγ2的影響而減少[30]。PPARγ可改善皮下脂肪細(xì)胞的形態(tài),并降低了促炎蛋白COX2的表達(dá)[31]。Belda等[32]在研究10e12z-cla對3t3-l1細(xì)胞分化和完全分化的影響。10e12z-cla在整個分化期都能降低3t3-l1脂肪細(xì)胞分化的程度。這些變化僅在分化過程的中途才變得明顯,并伴隨著COX2表達(dá)的增加而發(fā)生,COX2表達(dá)是NF-κB活性的標(biāo)志。10e12z-cla異構(gòu)體抑制脂肪分化并誘導(dǎo)人脂肪細(xì)胞中的NF-κB轉(zhuǎn)錄激活。在體外完全分化的人脂肪細(xì)胞中,10e12z-cla也能誘導(dǎo)環(huán)氧合酶-2[32]。10e12z-cla在3t3-l1前脂肪細(xì)胞分化過程中顯著抑制脂肪細(xì)胞表型并誘導(dǎo)炎癥標(biāo)志物COX2的表達(dá)[33]。但Yang等[34]在糖皮質(zhì)激素相關(guān)的脂肪細(xì)胞分化的研究中發(fā)現(xiàn)E4bp4E4(啟動子結(jié)合蛋白4)能夠通過反式抑制PTGS2而促進(jìn)脂肪沉積,而Fain等[35]的研究顯示,COX2與抗脂肪性前列腺素的生成具有功能耦合。該原因可能是COX2與抗脂肪性前列腺素的生成具有功能耦合,另一方面,在分化期24 h后,AA的補充導(dǎo)致了COX2轉(zhuǎn)錄水平的顯著降低,也可能是與脂肪細(xì)胞的功能調(diào)節(jié)有關(guān)。因此,PTGS2作為脂肪沉積的候選基因但還需更多試驗的驗證。
綜上,可初步確定DKKL1、ACACB及PTGS2為脂肪沉積的候選基因,對安格斯牛和中國西門塔爾牛脂肪沉積的后續(xù)研究和肉牛育種工作具有一定的積累意義和參考價值。