基于转录组测序挖掘奶牛乳脂代谢关键候选基因

2022-10-29王川川冯小芳禹保军顾亚玲

畜牧兽医学报 2022年10期

王川川，母童，冯小芳，禹保军，张娟，顾亚玲

(宁夏大学农学院，银川 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在奶牛乳脂合成中的作用机制还不明确，对其进行乳脂合成方面的深入探究是非常有必要的。