基于轉(zhuǎn)錄組測序挖掘奶牛乳脂代謝關(guān)鍵候選基因

王川川，母 童，馮小芳，禹保軍，張 娟，顧亞玲

(寧夏大學(xué)農(nóng)學(xué)院，銀川 750021)

乳脂是牛奶中主要營養(yǎng)成分之一，也是影響乳品質(zhì)的重要指標(biāo)，它含有人體所必需的脂肪酸及多種脂溶性維生素，是一種優(yōu)質(zhì)的油脂資源，為人體提供大量的營養(yǎng)物質(zhì)與能量。在生產(chǎn)中，乳脂主要用于制作黃油和酸奶。乳脂合成是一個復(fù)雜的生物學(xué)過程，在反芻動物乳腺中，乳脂合成的1/2是來自奶牛每天的飼喂飲食，其余是由乳腺細(xì)胞自身合成。在合成與分泌過程中，不同種類與數(shù)量的脂肪酸會形成不同分子量和飽和度的甘油三酯，而用于甘油三酯合成的脂肪酸的來源主要有兩個途徑：一種是乳腺上皮細(xì)胞的從頭合成，另一種是從血液中吸收而來。乳脂率的高低受到多種因素的影響，如遺傳、飼養(yǎng)管理和營養(yǎng)等，其中參與乳脂合成的酶和轉(zhuǎn)錄因子的表達(dá)對乳脂合成具有極其重要的作用。此外，許多重要基因也參與到脂質(zhì)合成中，并形成一個復(fù)雜的基因調(diào)控網(wǎng)絡(luò)。隨著分子生物學(xué)技術(shù)的發(fā)展，對乳脂合成的深入研究已經(jīng)揭示了脂肪酸的從頭合成、血液中脂肪酸的攝取、轉(zhuǎn)運與去飽和、脂滴的形成、甘油三酯的合成與運輸?shù)冗^程中關(guān)鍵基因及信號通路的調(diào)控作用。如今，牛奶中乳脂的含量不僅是乳品核心競爭力的重要指標(biāo)，也是奶牛育種的目標(biāo)性狀，如何提高奶牛乳脂含量，一直是研究人員關(guān)注的熱點。

自從2006年完成牛的全基因組測序工作至今，國內(nèi)外許多學(xué)者在與奶牛泌乳相關(guān)的轉(zhuǎn)錄組學(xué)、蛋白組學(xué)和代謝組學(xué)等方面開展了眾多研究并取得了重要成果。轉(zhuǎn)錄組測序(RNA-seq)是利用二代高通量測序技術(shù)更全面、快速地獲取某一物種特定器官、組織或細(xì)胞在某一狀態(tài)下的幾乎所有的轉(zhuǎn)錄本，包括mRNA和一些非編碼RNA。自從發(fā)現(xiàn)RNA是基因組和蛋白質(zhì)組之間的關(guān)鍵中間體以來，轉(zhuǎn)錄本鑒定和基因表達(dá)的量化一直是分子生物學(xué)的核心活動。目前，高通量測序技術(shù)已廣泛應(yīng)用于各種家畜組織差異表達(dá)基因的分析，例如在牛的胚胎、體細(xì)胞以及豬、羊等各種動物組織。同時，轉(zhuǎn)錄組技術(shù)已廣泛應(yīng)用于奶牛乳腺組織的研究，Dai等使用RNA-seq轉(zhuǎn)錄組學(xué)分析鑒定了來自泌乳期和非泌乳期奶牛乳腺組織中的差異表達(dá)基因，Gao等和Lin等也利用轉(zhuǎn)錄組測序技術(shù)對奶牛乳腺組織進(jìn)行了轉(zhuǎn)錄組分析。但有關(guān)于影響奶牛高、低乳脂率的泌乳機制方面的轉(zhuǎn)錄組研究相對較少。因此，通過RNA-seq闡明牛乳腺轉(zhuǎn)錄組對于挖掘影響奶牛乳脂相關(guān)的候選基因至關(guān)重要。

本研究對課題組前期基于高、低乳脂率(MFP)奶牛乳腺上皮細(xì)胞(BMECs)Illumina PE150測序結(jié)果中的mRNA表達(dá)譜數(shù)據(jù)進(jìn)行綜合分析，探索與MFP相關(guān)的標(biāo)志性基因和重要功能富集途徑，為理解奶牛乳脂合成過程中復(fù)雜的生物學(xué)過程和機制提供理論依據(jù)。

1 材料與方法

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

本研究中14 561條mRNA表達(dá)譜數(shù)據(jù)來源于課題組前期基于高、低MFP BMECs Illumina PE150的測序結(jié)果。試驗樣本來自寧夏農(nóng)墾賀蘭山奶業(yè)茂盛牧場，試驗牛飼喂相同的平衡全混合日糧(表1)。選取飼養(yǎng)背景一致，月齡相近，處于泌乳中后期的荷斯坦奶牛245頭。采集每頭牛早、中、晚班次的乳樣進(jìn)行奶牛生產(chǎn)性能測定(DHI)，篩選體細(xì)胞數(shù)在100 000個·mL以內(nèi)且MFP具有極端差異的荷斯坦奶牛各4頭(表2)。無菌采集新鮮乳液分離BMECs，待細(xì)胞鑒定完成后由北京諾禾致源科技股份有限公司進(jìn)行細(xì)胞總RNA提取、文庫構(gòu)建和轉(zhuǎn)錄組測序。

表1 日糧成分及營養(yǎng)成分(以干物質(zhì)為基礎(chǔ))百分比Table 1 Ingredient and nutrient composition in diet (dry matter basis) %

表2 荷斯坦牛高低乳脂率和SCCTable 2 High and low MFP and SCC of Holstein cattle

1.2 表達(dá)譜數(shù)據(jù)的篩選及富集分析

基于課題組前期14 561條mRNA表達(dá)譜數(shù)據(jù)，將<0.05和|logFoldChange|≥1.5設(shè)置為差異表達(dá)基因篩選的閾值，使用在線網(wǎng)站KOBAS(http://kobas.cbi.pku.edu.cn)數(shù)據(jù)庫對所有差異表達(dá)基因進(jìn)行GO功能注釋和KEGG通路富集分析，再以<0.05為顯著性閾值篩選顯著富集的GO條目和KEGG通路。利用北京諾禾致源科技股份有限公司提供的在線售后工具平臺(https://magic.novogene.com/customer/ main#/ home/2 d9 dc26 d1e059b93 1b9ac5364)對mRNA表達(dá)譜數(shù)據(jù)進(jìn)行相關(guān)性分析、表達(dá)量熱圖聚類、火山圖分析等。使用String數(shù)據(jù)庫進(jìn)行差異表達(dá)基因的蛋白互作網(wǎng)絡(luò)分析。

1.3 細(xì)胞及組織總RNA提取、質(zhì)量檢測及cDNA合成

使用課題組前期保存的高、低乳脂組奶牛BMECs的 cDNAs驗證轉(zhuǎn)錄組數(shù)據(jù)的準(zhǔn)確性。奶牛的小腸、肝、乳腺、心、子宮、腎、卵巢組織均為課題組前期采集(3頭泌乳中期的荷斯坦奶牛)。按照Trizol試劑盒(TaKaRa，大連)使用說明提取奶牛不同組織總RNA。使用紫外可見光分光光度計(Eppendorf BioSpectrometer basic)進(jìn)行濃度和純度檢測，總RNA完整性利用1%瓊脂糖凝膠電泳進(jìn)行檢測。挑取完整性及純度較好的RNA利用諾唯贊生物公司反轉(zhuǎn)錄試劑盒將其反轉(zhuǎn)錄為cDNAs，反轉(zhuǎn)錄完成后置于-20 ℃冰箱保存?zhèn)溆谩?/p>

1.4 差異表達(dá)基因?qū)崟r熒光定量 PCR驗證

1.4.1 引物設(shè)計與合成 隨機選擇8個差異表達(dá)基因(、、4、2G16、、、4D、1)及篩選出的脂代謝相關(guān)關(guān)鍵候選差異表達(dá)基因(2、2、4、5)，以牛的基因為內(nèi)參進(jìn)行qRT-PCR。根據(jù)NCBI上查找到的基因序列，運用Primer5.0設(shè)計熒光定量引物(表3)，由陜西中科羽瞳生物公司代為合成。

表3 引物序列及退火溫度Table 3 Primer sequence and annealing temperature

1.4.2 熒光定量PCR 采用實時熒光定量PCR驗證基因的表達(dá)水平。每個樣本設(shè)置3個技術(shù)重復(fù)，按照艾科瑞生物公司的產(chǎn)品SYBR Green Pro Taq HS預(yù)混型qPCR試劑盒進(jìn)行檢測，qRT-PCR反應(yīng)體系：ddHO 8 μL，2×SYBR Green PCR Master Mix 10 μL，cDNA 1 μL，上、下游引物各0.5 μL。qRT-PCR反應(yīng)程序：95 ℃ 預(yù)變性30 s；95 ℃變性5 s，退火30 s，60 ℃延伸5 s，40個循環(huán)。

1.5 數(shù)據(jù)分析

用Bio-Rad CFX Manager3.1軟件收集熒光定量結(jié)果的數(shù)據(jù)，通過Microsoft Excel 2021對數(shù)據(jù)進(jìn)行整理。不同組織基因的相對表達(dá)水平采用2法進(jìn)行分析，結(jié)果以“平均值±標(biāo)準(zhǔn)差”來表示。最后利用SAS8.2軟件進(jìn)行統(tǒng)計分析，<0.05認(rèn)為具有生物統(tǒng)計學(xué)意義。

2 結(jié) 果

2.1 樣本相關(guān)性分析及差異表達(dá)基因篩選和熱圖聚類

樣品間基因表達(dá)水平相關(guān)性是檢驗試驗可靠性和樣本選擇與否的重要指標(biāo)。本研究根據(jù)各樣本所有基因的FPKM值，計算組內(nèi)及組間樣本的相關(guān)性系數(shù)，并繪制成熱圖(圖1A)，圖中可直觀看出組內(nèi)相關(guān)系數(shù)均在0.90以上，表明樣品之間表達(dá)模式的相似度較高，同時組間樣本也存在差異。本研究以|logFoldChange|≥1.5、<0.05作為差異表達(dá)基因的篩選標(biāo)準(zhǔn)，結(jié)果如圖1C所示，高乳脂率組和低乳脂率組中共有578個差異表達(dá)mRNA，其中高乳脂組(試驗組)有332個基因顯著上調(diào)(低乳脂組為對照組)，246個基因顯著下調(diào)。聚類熱圖是差異表達(dá)基因的另一種展示方式，如圖1B所示，表達(dá)模式相近的基因被聚在一起，這些基因可能具有共同的功能或參與到共同的代謝途徑。

A.各樣本相關(guān)性分析；B. 差異表達(dá)基因聚類熱圖；C. 差異表達(dá)基因火山圖(HMF_VS_LMF)A. Pearson correlation between samples；B. Heat map of clustering for differentially expressed genes；C. Differentially expressed genes volcano map(HMF_VS_LMF)圖1 各樣本相關(guān)性分析及差異表達(dá)基因分析Fig.1 Pearson correlation between samples and analysis of differentially expressed genes

2.2 差異表達(dá)基因富集分析及脂代謝相關(guān)候選基因篩選

對578個差異表達(dá)基因進(jìn)行GO和KEGG富集分析，結(jié)果共顯著富集到366個GO條目，包含224個生物學(xué)過程(BP)，58個細(xì)胞組分(CC)及84個分子功能(MF)。圖2A中列出了顯著富集到與脂代謝相關(guān)的前30個GO條目，其中與乳脂代謝關(guān)系最為密切的條目是中長鏈脂肪酸的運輸、乳腺肺泡發(fā)育、花生四烯酸的結(jié)合和磷脂酰肌醇-4,5-二磷酸結(jié)合。利用在線數(shù)據(jù)庫KOBAS(http://kobas.cbi.pku.edu.cn)對差異表達(dá)基因進(jìn)行KEGG富集分析，結(jié)果發(fā)現(xiàn)有47個信號通路被顯著富集(<0.05)，圖2B中列出了與脂代謝相關(guān)的前15的信號通路，其中與脂代謝相關(guān)的代謝通路為Apelin信號通路、Hippo信號通路、脂肪細(xì)胞中脂解的調(diào)節(jié)及MAPK信號通路等。之后結(jié)合GO和KEGG富集分析結(jié)果，最終篩選出可能調(diào)控乳脂代謝的4個候選差異表達(dá)基因，即2、2、4和5。其中5、2基因表達(dá)下調(diào)，2、4基因表達(dá)上調(diào)。值得關(guān)注的是，2和4基因顯著富集至脂代謝相關(guān)BP和MF，同時富集在脂代謝相關(guān)KEGG通路中(圖3A、B)，提示其參與調(diào)控乳脂代謝的可能性更大(表4)。

A.與脂代謝相關(guān)的差異表達(dá)基因富集分析中前30個GO條目； B.與脂代謝相關(guān)的差異表達(dá)基因富集分析中前15條KEGG通路圖A. The top 30 GO terms in the enrichment analysis of differentially expressed genes in lipid metabolism; B. The top 15 KEGG pathways in the enrichment analysis of differentially expressed genes in lipid metabolism圖2 差異表達(dá)基因GO及KEGG代謝通路富集分析Fig.2 Enrichment analysis of GO and KEGG metabolic pathways of differentially expressed genes

A. Apelin信號通路；B.調(diào)節(jié)脂肪細(xì)胞的脂分解信號通路A. Apelin signaling pathway; B. Regulation of lipolysis in adipocytes圖3 脂代謝相關(guān)通路Fig.3 Lipid metabolism related pathways

表4 與脂代謝相關(guān)的差異表達(dá)基因匯總Table 4 Summary of differential genes associated with lipid metabolism

2.3 蛋白互作網(wǎng)絡(luò)分析

對篩選出可能調(diào)控乳脂代謝的差異表達(dá)基因2、2、4、5分別進(jìn)行蛋白互作分析，結(jié)果如圖4所示，蛋白互作網(wǎng)絡(luò)圖中各有10個節(jié)點，置信度為0.7(高等)，線的粗細(xì)表示了數(shù)據(jù)支撐的強度。從圖中可以看出，2、2、4基因與多個脂代謝相關(guān)基因存在較強的互作關(guān)系，在調(diào)控乳脂代謝方面它們可能共同扮演著重要的角色，5基因與調(diào)控乳脂代謝方面的基因互作較少，數(shù)據(jù)支撐相對較少，但仍然具有一定的參考價值。

圖4 差異表達(dá)基因蛋白網(wǎng)絡(luò)互作Fig.4 Differentially expressed genes protein network interaction

2.4 總RNA質(zhì)量分析

各組織樣本總RNA的提取質(zhì)量會影響到后續(xù)實時熒光定量的結(jié)果，經(jīng)1%的瓊脂糖凝膠電泳檢測后，如圖5所示，1～7分別是奶牛的小腸、肝、乳腺、心、子宮、腎、卵巢7個組織RNA凝膠電泳檢測，可以清晰看到兩條較明顯的條帶分別為28S、18S，并無明顯降解，表明樣品提取的總RNA的完整性較好。使用核酸濃度檢測儀檢測總RNA的OD/OD值均在1.8～2.0之間，說明RNA純度較高可滿足后續(xù)試驗的要求。

1～7. 奶牛小腸、肝、乳腺、心、子宮、腎、卵巢不同組織的RNA條帶1-7. RNA bands extracted from small intestine, liver, breast, heart, uterus, kidney and ovary of dairy cows, respectively 圖5 不同組織RNA提取電泳圖Fig.5 Electrophoresis of RNA extracted from different tissues

2.5 乳脂代謝關(guān)鍵基因組織表達(dá)譜分析

對2、2、4、5基因分別進(jìn)行不同組織的qRT-PCR擴增，定量結(jié)果如圖6所示，所有基因的擴增曲線均呈現(xiàn)S型且具有指數(shù)增長期和平臺期，均在15～30個循環(huán)之間起峰，表明擴增效率較高。基因的熔解曲線均為單峰，且起峰迅速形狀狹窄，說明不存在非特異性片段，擴增特異性較好。qRT-PCR結(jié)果發(fā)現(xiàn)，以小腸作為參照，2基因在7個組織中都有表達(dá)，表達(dá)水平從高到低依次為卵巢、乳腺、子宮、小腸、腎、心及肝。其中，在卵巢組織中表達(dá)水平顯著高于其他組織，在肝中的表達(dá)水平最低。2基因在小腸中表達(dá)水平最高，與心和腎差異不顯著，與其余組織表達(dá)水平差異顯著(<0.05)，4在乳腺組織中表達(dá)水平最高，與其他組織差異顯著(<0.05)，其他組織之間的表達(dá)量差異不顯著。5基因在心中表達(dá)水平最高，乳腺組織次之，與其余組織之間表達(dá)水平差異顯著(<0.05,圖7)。

圖6 Real-time PCR的熔解曲線和擴增曲線Fig.6 Melting curves and amplification curves of real-time PCR

肩標(biāo)不同字母表示差異顯著(P<0.05)，相同小寫字母表示差異不顯著(P>0.05)Different letters on the shoulders means significant difference between the treatments(P<0.05), same letter on the shoulders indicate not significant difference between treatments(P>0.05)圖7 候選基因組織表達(dá)譜Fig.7 Tissue expression profiles of candidate genes

2.6 轉(zhuǎn)錄組數(shù)據(jù)的qRT-PCR驗證

本研究隨機選擇了8個差異表達(dá)基因，在高、低組試驗樣品中(3個重復(fù))運用qRT-PCR分析了它們的相對表達(dá)水平，通過logFoldChange對差異倍數(shù)進(jìn)行轉(zhuǎn)換。結(jié)果如圖8所示，3個基因(4、4D、)表達(dá)上調(diào)，5個基因(、、2G16、1、)表達(dá)下調(diào)。qRT-PCR試驗結(jié)果與RNA-seq測序結(jié)果均表現(xiàn)一致，證實了轉(zhuǎn)錄組測序結(jié)果的可靠性。

圖8 qRT-PCR驗證差異表達(dá)基因Fig.8 The differentially expressed genes confirmed by qRT-PCR

3 討 論

乳脂是乳品質(zhì)調(diào)控的主要目標(biāo)性狀，也是反映奶產(chǎn)品行業(yè)經(jīng)濟(jì)效益的主要指標(biāo)之一。本研究以<0.05和|logFoldChange|≥1.5為閾值，在荷斯坦奶牛高、低乳脂組之間共發(fā)現(xiàn)578個差異表達(dá)基因，其中上調(diào)基因332個，下調(diào)基因246個。功能富集分析發(fā)現(xiàn)，差異表達(dá)基因顯著富集于長鏈脂肪酸的運輸、磷脂酰肌醇介導(dǎo)的信號及乳腺肺泡的發(fā)育等GO條目。此外，與脂代謝密切相關(guān)的通路有15條，包括脂肪細(xì)胞脂解的調(diào)節(jié)、磷脂酶D信號通路、河馬信號通路及Apelin信號通路等。在Apelin信號通路中，5基因通過對cAMP的調(diào)控間接對抑制脂質(zhì)分解起到促進(jìn)的作用。2基因通過AMPK信號通路間接促進(jìn)脂肪的分解及褐色脂肪的形成。4參與甘油酯代謝，生成脂肪酸。這說明5、2及4基因的表達(dá)與脂代謝密切相關(guān)。此外，利用String蛋白質(zhì)互作數(shù)據(jù)庫對候選基因蛋白進(jìn)行蛋白互作分析，2、2、4基因與多個脂代謝相關(guān)基因存在較強的互作關(guān)系，它們在調(diào)控乳脂代謝方面發(fā)揮著相似的功能。因此，本研究選取2、2、4、5基因作為乳脂代謝關(guān)鍵功能候選基因。

2目前在牛上未見報道，Liu等研究發(fā)現(xiàn)2是鑒定DCIS(導(dǎo)管原位癌)非侵襲性乳腺癌起始的關(guān)鍵基因，在本研究中其在奶牛乳腺中的表達(dá)水平較高，顯著高于小腸、肝、心和腎。2基因?qū)η卮ㄅ２糠煮w尺及肉質(zhì)性狀有顯著影響，可作為秦川肉牛新品種早期選擇的候選基因和分子輔助標(biāo)記。Khan等研究發(fā)現(xiàn)，在奶牛泌乳早期，AMPK的完全激活需要1和2，2表達(dá)的變化表明其參與了脂解的過程，脂肪生成的速率降低，從而刺激脂解的作用增加。此外，在牦牛泌乳周期中，2在整個哺乳期內(nèi)上調(diào)。蛋白互作分析進(jìn)一步發(fā)現(xiàn)，2與、、1、2基因均有互作，這些都是與脂代謝密切相關(guān)的基因。本研究發(fā)現(xiàn)，2在乳腺中的表達(dá)水平處于中間水平，但其在乳脂合成中的作用是不可忽視的。Cai等的研究表明，奶牛產(chǎn)奶量與生育力之間存在著遺傳負(fù)相關(guān)性，使得在奶牛的育種中難以同時提高產(chǎn)奶量和生育力，奶牛產(chǎn)奶性狀的快速改善伴隨著奶牛生育能力的下降，這種相關(guān)性的遺傳基礎(chǔ)仍然知之甚少。5作為奶牛生育能力的候選基因，有助于剖析母牛生育能力的遺傳學(xué)。脂肪酸結(jié)合蛋白(fatty acid binding protein，F(xiàn)ABP)可與長鏈脂肪酸結(jié)合，并將脂肪酸從質(zhì)膜轉(zhuǎn)運到脂肪酸氧化位點及其與甘油三酯和磷酸的合成位點，調(diào)節(jié)細(xì)胞內(nèi)脂肪的合成與轉(zhuǎn)運，并和其他脂質(zhì)類調(diào)節(jié)因子調(diào)節(jié)細(xì)胞內(nèi)信號通路。在奶牛的乳腺組織中，F(xiàn)ABP 家族的3、4及5表達(dá)量最高。其中，脂肪酸結(jié)合蛋白4(FABP4)最早被Kacey等在脂肪細(xì)胞中發(fā)現(xiàn)，是畜禽脂肪沉積與代謝的熱門候選基因。本研究發(fā)現(xiàn)，4在乳腺組織中的表達(dá)尤為突出，顯著高于其它組織，但目前4在奶牛乳脂合成中的作用機制還不明確，對其進(jìn)行乳脂合成方面的深入探究是非常有必要的。

4 結(jié) 論

本研究共篩選到578個差異表達(dá)基因，其中上調(diào)基因有332個，下調(diào)基因有246個。通過功能富集分析和組織表達(dá)譜分析得到了4個乳脂代謝相關(guān)的關(guān)鍵候選基因，分別為2、2、4、5。本研究為闡明奶牛乳脂代謝的分子調(diào)控機制和調(diào)控網(wǎng)絡(luò)提供了基礎(chǔ)資料，也為荷斯坦奶牛分子選育工作提供了素材。