基于轉錄組測序挖掘奶牛乳脂代謝關鍵候選基因

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

(寧夏大學農學院，銀川 750021)

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

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

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

1 材料與方法

1.1 數據來源

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

表1 日糧成分及營養成分(以干物質為基礎)百分比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 表達譜數據的篩選及富集分析

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

1.3 細胞及組織總RNA提取、質量檢測及cDNA合成

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

1.4 差異表達基因實時熒光定量 PCR驗證

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

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

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

1.5 數據分析

用Bio-Rad CFX Manager3.1軟件收集熒光定量結果的數據，通過Microsoft Excel 2021對數據進行整理。不同組織基因的相對表達水平采用2法進行分析，結果以“平均值±標準差”來表示。最后利用SAS8.2軟件進行統計分析，<0.05認為具有生物統計學意義。

2 結 果

2.1 樣本相關性分析及差異表達基因篩選和熱圖聚類

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

A.各樣本相關性分析；B. 差異表達基因聚類熱圖；C. 差異表達基因火山圖(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 各樣本相關性分析及差異表達基因分析Fig.1 Pearson correlation between samples and analysis of differentially expressed genes

2.2 差異表達基因富集分析及脂代謝相關候選基因篩選

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

A.與脂代謝相關的差異表達基因富集分析中前30個GO條目； B.與脂代謝相關的差異表達基因富集分析中前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 差異表達基因GO及KEGG代謝通路富集分析Fig.2 Enrichment analysis of GO and KEGG metabolic pathways of differentially expressed genes

A. Apelin信號通路；B.調節脂肪細胞的脂分解信號通路A. Apelin signaling pathway; B. Regulation of lipolysis in adipocytes圖3 脂代謝相關通路Fig.3 Lipid metabolism related pathways

表4 與脂代謝相關的差異表達基因匯總Table 4 Summary of differential genes associated with lipid metabolism

2.3 蛋白互作網絡分析

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

圖4 差異表達基因蛋白網絡互作Fig.4 Differentially expressed genes protein network interaction

2.4 總RNA質量分析

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

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 乳脂代謝關鍵基因組織表達譜分析

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

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

肩標不同字母表示差異顯著(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 候選基因組織表達譜Fig.7 Tissue expression profiles of candidate genes

2.6 轉錄組數據的qRT-PCR驗證

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

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

3 討 論

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

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

4 結 論

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