荷斯坦奶牛肝脏组织中与泌乳时期及繁殖力相关的基因共表达网络构建

张志飞,唐雪颖,,闵 力,童 雄,陈卫东,巨向红,李大刚*

(1.广东省农业科学院动物科学研究所 广东省畜禽育种与营养研究重点实验室,广州 510000;2.广东海洋大学滨海农业学院,湛江 524000)

荷斯坦牛(Holstein,Bostaurus)是全球最重要的奶牛品种之一,据记载起源于至少2 000年前,原产于荷兰北部的西弗里斯兰省(Friesland)和北荷兰省(North Holland),以及德国北部的荷斯坦(Holstein)地区[1]。中国在十九世纪中叶从美国引进荷斯坦牛,从20世纪50年代起有计划开展品种培育工作,1992年经原农业部批准命名为“中国荷斯坦牛”[1]。荷斯坦奶牛具有产奶性能好、适应性强和易于饲养等特点,因而被广泛应用于奶业生产[2]。荷斯坦奶牛来源的鲜奶、奶粉、酸奶等产品是人类饮食中优质蛋白的重要来源,为改善人们的健康状况做出了重要贡献[3]。营养不良、热应激、乳房炎、围产应激等生理应激或环境应激往往对奶牛的生产和繁殖性能造成不良影响,是目前阻碍奶牛产业健康发展的关键问题[4-6]。

肝脏参与奶牛生理代谢与调控相关的众多关键过程,包括葡萄糖代谢、能量稳态、营养因子的转运、免疫系统调节、生长发育刺激及有害毒素降解等[5]。因此,肝脏在调控不同生理时期奶牛生产及繁殖性能等方面起着关键作用[7]。在奶牛的围产期,胎儿生长对营养物质需求激增,但干物质采食量受胃肠道系统影响显著降低,极易发生能量负平衡(NEB)[8]。此时,肝脏中葡萄糖合成、胆固醇代谢及脂肪生产过程增强,同时对脂肪和氨基酸的动员增加,维持能量平衡的同时产生大量酮体,损伤免疫系统,诱发产褥热、酮病、胎衣不下、流产等代谢性疾病[9]。在奶牛的产奶高峰期,奶牛由于大量合成乳汁与泌乳,同样经历着能量负平衡,当大量脂肪被消耗时,脂肪组织合成的瘦素水平降低,使体内促性腺激素释放激素分娩减少,抑制卵泡发育,引起高产奶牛低繁殖力问题[10]。

加权基因共表达网络分析(weighted gene co-expression network analysis, WGCNA)是一种系统生物学方法,可用于对表达特性上高度相关的基因进行聚类(模块),通过与样本的性状进行关联分析,可用来筛选出影响或调控生理性状的特征模块及Hub基因,进而可结合其它生物学技术鉴定性状标记基因或靶点基因[11-12]]。王子渲等[13]的研究利用RNA-seq和WGCNA分析,挖掘了肉鸡脾脏中与热应激相关性状显著关联的Hub基因。李晓波等[14]的研究利用WGCNA和GSEA方法鉴定了与中卫山羊羊毛弯曲相关的Hub基因。王丽敏等[15]利用WGCNA和PPI网络技术鉴定了山羊金黄色葡萄球菌型乳腺炎关键应答基因。由此可见,WGCNA技术在畜禽关键性状调控或标记基因的挖掘方面具有广阔的应用前景。

因此,本研究旨在利用WGCNA结合PPI分析技术构建荷斯坦奶牛肝脏中与泌乳时期和繁殖能力性状相关的基因共表达网络,并分析和鉴定调控网络中的关键基因和调控因子。这将有助于我们深入了解荷斯坦奶牛生产和繁殖性能的调控机理,从而为改善荷斯坦奶牛的繁殖能力和产奶效率提供理论依据和基础数据。

1 材料与方法

1.1 数据来源

本试验的数据来源于GEO数据库(https:∥www.ncbi.nlm.nih.gov/geo/)的数据集GSE62159。该数据集共包含48头健康荷斯坦奶牛的肝脏转录组测序数据,按照泌乳时期划分:妊娠末期(分娩前18 d)16头,泌乳早期(分娩后1 d)16头,泌乳中期(分娩后147 d)16头;按照繁殖力(基于生产间隔及世代信息的育种值估计值)划分:高繁殖力奶牛24头,低繁殖力奶牛24头[16]。数据集中的荷斯坦奶牛饲养于爱尔兰(55°10′N 8°16′W),其系谱信息来自于爱尔兰牛产业联盟(the Irish Cattle Breeding Federation,ICBF)的国家奶牛数据库,试验用牛的系谱记录开始于2007年,2008年获得首个产奶性能、繁殖性状(产犊间隔)的育种值(EBV),测序样品采集于2011年,其中高繁殖力组奶牛属于产犊间隔(85.6 d)育种值排名前20%的家系,低繁殖力组奶牛属于产犊间隔(113.8 d)育种值排名后5%的家系[16-17]。

1.2 构建基因共表达网络

使用R语言中的WGCNA包[18],根据hclust函数,对样本数据进行分析将基因表达模式类似的样本进行聚类,并且除去离群样本;利用pick-SoftThreshold函数选择出合适的软阈值;利用adjacency函数构建拓扑重叠矩阵(TOM)并对基因进行聚类,得到共表达模块;将共表达模块与表型信息矩阵进行关联分析,依据相关系数r值与P值(P<0.05 表示差异显著)选择与不同泌乳时期或繁殖力相关的模块为目标模块;选择模块内连接度排名前30的基因作为Hub基因,构建基因共表达网络。

1.3 模块的富集分析

对每个目标模块的Hub基因使用KOBAS 3.0(http:∥kobas.cbi.pku.edu.cn/kobas3/)进行KEGG富集分析;使用DAVID网站(https:∥david.ncifcrf.gov)进行GO富集分析,并使用R语言将富集结果可视化。

1.4 构建蛋白互作网络和确认目标基因

对所有目标模块中的基因使用String 11.0网站(http:∥string-db.org/)进行蛋白互作网络(protein-protein interaction networks,PPI)的构建。使用Cytoscape软件根据PPI网络中节点的度值degree选出前30个基因作为核心基因。将模块中的Hub基因与所得的核心基因取交集,即得到与不同泌乳时期或繁殖调控相关的关键基因。

2 结 果

2.1 基因共表达网络的构建

对48个样本进行聚类分析,发现3个离群样本,对离群样本进行剔除后选择45个样本进行分析,样本树状与性状热图如图1A所示。根据数据自由度和连通度分析确定最佳软阈值数值取12,进行无尺度网络的构建。对模块特征基因进行聚类,对红线以下(高度<0.05)的模块进行合并(图1B),共得到14个模块(图1C)。对模块特征值与表型信息值进行相关性分析,选择与泌乳早期(EL)相关性最强tan模块(r=0.54,P=1×10-4),与泌乳中期(ML)相关性最强greenyellow模块(r=0.58,P=3×10-5),和与繁殖力相关性最强black模块(r=-0.57,P=4×10-5)为目标模块(图1D)。

A. 样本聚类树图;B.模板特征基因聚类树;C.模块聚类树;D.模板-性状相关性图A. Sample dendrogram and trait heatmap; B. Template feature gene clustering tree; C. Gene dynamic splicing clustering tree; D. Template and phenotype association analysis

2.2 泌乳早期肝脏中相关基因的共表达网络及关键基因分析

对tan模块进行GS-MM相关性验证,结果表明tan模块中基因与泌乳早期相关性显著(P=1.5×10-9,图2A),模块内连接度排名前30的Hub基因见表1。

表1 各模块中的hub基因

A.GS-MM分析散点图;B. tan模块KEGG富集分析;C. tan模块GO富集分析;D. PPI网络核心基因互作图A. Scatter plots of GS-MM analysis; B. Tan module KEGG enrichment analysis; C. Tan module GO enrichment analysis; D. PPI network core genes interaction diagram

KEGG分析结果显示,与tan模块关联的主要功能信号通路包括:内质网蛋白质加工、代谢信号通路、N-聚糖生物合成、蛋白质分泌、不同亚型N-聚糖生物合成、核苷酸糖生物合成、补体和凝血级联、丙酸代谢、丙酮酸代谢、半胱氨酸和蛋氨酸代谢等(图2B)。GO分析结果表明,与tan模块关联的主要生物学过程有:信号肽加工、内质网应激的响应、蛋白质聚合、蛋白N连接糖基化、血小板激活、内源性肽酶活性的负调控、纤维蛋白溶解、介导囊泡由内质网到高尔基体的转运、羧酸代谢过程和急性期反应;与tan模块关联的主要细胞组分包括:信号肽酶复合物、低聚糖转移酶复合物、高密度脂蛋白颗粒、纤维蛋白原复合物、细胞外隙、内质网膜、内质网腔等;与tan模块关联的主要分子功能包括:结构分子活性、丝氨酸型内肽酶抑制剂活性、受体结合、蛋白质二硫氧化还原酶活性、蛋白质二硫异构酶活性、L-乳酸脱氢酶活性、整合素结合、同蛋白结合、酶激活剂活性和核糖核酸内酶活性等(图2C)。

Tan模块中Hub基因与PPI网络分析核心基因取交集得到的关键基因有12个,分别为:RPN1、SEC61A1、SEC61B、SEC61G、SSR1、SSR3、STT3A、DAD1、DDOST、ERLEC1、HM13和OSTC(图2D)。

2.3 泌乳中期与妊娠后期肝脏中相关基因的共表达网络及关键基因分析

对greenyellow模块进行GS-MM相关性验证,结果表明greenyellow模块中基因与泌乳中期相关性显著(P=0.002 5,图3A)。模块内连接度排名前30的Hub基因见表1。

A. GS-MM分析散点图;B. Greenyellow模块KEGG富集分析;C. Greenyellow模块GO富集分析;D. PPI网络核心基因互作图A. Scatter plots of GS-MM analysis; B. Greenyellow module KEGG enrichment analysis; C. Greenyellow module GO enrichment analysis; D. PPI network core genes interaction diagram

KEGG分析结果显示,与greenyellow模块关联的主要功能信号通路包括:病毒性心肌炎、Ⅰ型糖尿病、肺结核、弓形体病、金黄色葡萄球菌感染、类风湿性关节炎、吞噬小体、百日咳、溶酶体、白细胞跨内皮细胞迁移、利什曼病、炎症性肠病、人类T细胞白血病病毒Ⅰ型感染、造血细胞谱系、移植物抗宿主病、EB病毒感染、补体和凝血级联、细胞黏附分子、抗原处理和呈递以及同种异体移植排斥等(图3B)。GO分析结果表明,与greenyellow模块关联的主要生物学过程有:细胞形状调控、磷脂酰肌醇3-激酶信号的正向调控、血管生成的正向调控、炎症反应、免疫反应、细胞迁移、通过MHC II类抗原和外源性肽抗原的呈递、抗原加工呈递和血管生成等;与greenyellow模块关联的主要细胞组分包括:受体复合体、MHCⅡ类蛋白复合体、膜筏、溶酶体、免疫突触、粘着斑、细胞外间隙、质膜外侧、细胞表面和基底外侧质膜等;与greenyellow模块关联的主要分子功能包括:结构跨膜受体蛋白酪氨酸激酶活性、受体结合、蛋白质同二聚体活性、蛋白结合、整合素结合、肝素结合、钙依赖蛋白结合、钙粘蛋白结合、淀粉样蛋白结合和肌动蛋白丝结合等(图3C)。

Greenyellow模块中Hub基因与PPI网络分析核心基因取交集得到的关键基因有6个,分别为ITGAL、ITGB2、LAPTM5、PTPRC、C3AR1和CTSS(图3D)。

2.4 繁殖力相关基因在肝脏中的共表达网络及目标基因分析

对black模块进行GS-MM相关性验证,结果表明black模块中基因与低繁殖力相关性显著(P=3.8×10-21,图4A)。模块内连接度排名前30的Hub基因见表1。

A. GS-MM分析散点图;B. black模块KEGG富集分析;C. black模块GO富集分析;D.PPI网络核心基因互作图A. Scatter plots of GS-MM analysis; B. Black module KEGG enrichment analysis; C. Black module GO enrichment analysis; D. PPI network core genes interaction diagram

KEGG分析结果显示,与black模块关联的主要功能信号通路包括:淀粉和蔗糖代谢、胰岛素抵抗、粘着斑、肌动蛋白细胞骨架的调节、EGFR酪氨酸激酶抑制剂耐药性、胰岛素信号通路、癌症中的蛋白聚糖、癌症中的中枢碳代谢、Ⅱ型糖尿病、甲状腺激素信号通路、黏附连接、胃癌、黑素瘤、非小细胞肺癌、mTOR信号通路、FoxO信号通路、胰高血糖素信号通路、C型凝集素受体信号通路、癌症通路以及子宫内膜癌等(图4B)。GO分析结果显示,与black模块关联的主要生物学过程有:蛋白质磷酸化、蛋白去磷酸化、葡萄糖输入的正向调节、肽基-苏氨酸磷酸化、应力纤维组装负调控、Rho蛋白信号转导负调控、细胞迁移的负调节、眼睛的胚胎视网膜形态发生、胚胎眼形态发生以及细胞对紫外线的反应;与black模块关联的主要细胞组分包括:受体复合体、核浆、核斑点、高尔基体、初级内体、胞浆、细胞质、细胞表面、细胞-细胞连接和陷窝;与black模块关联的主要分子功能包括:锌离子结合、小GTPase结合、核糖核酸结合、蛋白丝氨酸/苏氨酸激酶活性、核小体依赖性ATP酶活性、甘露糖寡糖1,2-α-甘露糖苷酶活性、胰岛素受体结合、胰岛素结合、水解酶活性,水解O-糖基化合物以及ATP结合(图4C)。

Black模块中Hub基因与PPI网络分析核心基因取交集得到的关键基因有4个,分别为PDS5A、ROCK1、AQR、和LTN1(图4D)。

2.5 各模块筛选到的关键基因mRNA表达水平分析

对tan、greenyellow及balck模块中筛选到的肝脏组织中与泌乳时期或繁殖力相关的关键基因进行表达水平分析,并绘制基因在各样本中的表达水平热图(图5),结果显示:tan模块中的关键基因在妊娠后期表达水平最高、在泌乳中期表达水平最低;Greenyellow模块中的关键基因在泌乳早期表达水平较高、在泌乳中期表达水平较低;balck模块中的基因在高繁殖力组中的表达水平高于低繁殖力组。

3 讨 论

本研究通过构建荷斯坦奶牛肝脏组织中与泌乳时期以及繁殖力相关的基因共表达网络,得到了肝脏中与不同泌乳时期相关的tan模块、greenyellow模块以及与繁殖力正相关的black模块。且分别筛选出12、6及4个与目标性状相关的关键基因。

在泌乳早期,荷斯坦奶牛的肝脏需要加强对葡萄糖、脂肪酸和氨基酸等物质的代谢和利用,以保证能量和营养物质的供应[19]。此外,肝脏在泌乳早期还需要清除乳酸、亚硝酸和其它代谢废物,维持内环境的稳定性[20]。蛋白质糖基化是一种常见的蛋白质翻译后修饰类型,糖类物质在糖基转移酶催化下与蛋白质上特定氨基酸残基相连形成糖苷键[21]。与糖基分子互作的蛋白称为糖结合蛋白(glycan-binding protein,GBP),参与调控细胞信号识别与传递、细胞内吞、细胞生长、分化、凋亡等生物学过程[22]。肝脏组织发生病变时,GBP的结构及其次级产物与正常肝脏组织存在显著差异,因此能够作为判断特异性肝脏疾病的标志物[23]。本研究筛选出的与泌乳初期荷斯坦奶牛肝脏代谢关键基因中:RNP1编码一种核糖体结合蛋白,在肝细胞中参与蛋白质高效转运,与细胞凋亡和应激响应有关[24];SEC61A1、SEC61B、SEC61 G共同编码转录因子Sec61,在蛋白质合成与转运过程发挥重要作用[25];SSR1、SSR3编码两种糖蛋白,能够与Sec61发生蛋白互作,增强内质网活性、促进蛋白质运输[26];SST3A、DDOST均是蛋白质N-糖基化的催化酶基因[27];ERLC1、OSTC编码的蛋白均是内质网糖蛋白组分,参与糖蛋白修饰、内质网质量控制等过程[28-29];HM13参与糖酵解过程,是能量代谢稳态和氧化应激响应的重要标志物之一[30]。综合发现认为,在泌乳初期荷斯坦奶牛肝脏中,蛋白质合成与转运相关代谢过程旺盛,蛋白质糖基化这一蛋白质翻译后修饰类型可能是参与调控这一时期肝脏蛋白合成、能量代谢与应激响应等作用的关键生物学过程。泌乳初期肝脏糖基化相关酶及基因编码蛋白的表达及功能研究,是潜在的泌乳初期荷斯坦奶牛代谢稳态失衡的有效监测标准或重要调控途径之一。

泌乳中期与妊娠末期都是荷斯坦奶牛重要的生理时期,它们的肝脏代谢也存在许多异同点。在这两个阶段,肝脏组织均扮演着能量代谢、脂肪代谢及葡萄糖代谢的稳态调节者,同时也是机体免疫功能的重要参与者[31];不同点在于,奶牛在泌乳中期需要更多能量维持泌乳、同时需要调节血糖、血脂水平稳定以支持乳汁合成,而妊娠末期则需要更多营养以支持胎儿生长发育、能量供应不足时需要分解储存的糖原供能[32]。本研究发现,荷斯坦奶牛泌乳中期及妊娠末期与greenyellow模块均存在显著相关,而greenyellow模块中基因的功能主要富集在疾病相关通路,如金黄色葡萄球菌感染、机体的炎症反应、免疫反应以及急性期等。本模块中筛选到的6个关键基因中:ITGAL和ITGB2基因编码整合素,能够保护肝脏免受异常代谢或病变的损伤[33];LAPTM5编码一种在免疫细胞中优先表达的蛋白质,它与泛素连接酶的Nedd4家族相互作用,被鉴定为炎症信号通路的正调节剂以及预测高血压患者LVH的潜在生物标志物[34];PTPRC基因编码CD45蛋白,是免疫细胞的标记基因之一[35];C3AR1和CTSS基因均与炎症反应和免疫系统功能相关[36-37]。由此推测,无论处于泌乳中期或妊娠末期,奶牛机体均处于疾病易感状态;外界环境或内在因素诱导的炎症反应、免疫系统功能损伤或营养代谢性疾病作为威胁奶牛健康的主要原因,对奶牛肝脏功能调控造成巨大挑战。

目前认为,对产奶性状的高强度选择导致了奶牛繁殖性能衰退,但营养因素、疾病和健康状况、生殖管理等因素同样对奶牛繁殖性能具有重要影响[38]。本研究所使用的数据集中,评判奶牛个体繁殖力高低的主要依据是经过长期选择、不同家系的荷斯坦奶牛产犊间隔性状的估计育种值。可以看到,肝脏组织中与繁殖力相关的black模块中基因所富集到的多数信号通路包括淀粉和蔗糖代谢、胰岛素抵抗、GFR酪氨酸激酶抑制剂耐药性、胰岛素信号通路、Ⅱ型糖尿病、甲状腺激素信号通路、黏附连接、mTOR信号通路、FoxO信号通路、胰高血糖素信号通路等,与肝脏自身发挥的主要代谢功能高度符合,表明肝脏代谢能力与荷斯坦奶牛的繁殖能力高度相关,这也部分解释了肝脏代谢旺盛会促使奶牛体内与发情、排卵相关激素分解,进而减弱了奶牛繁殖能力的说法[39-40]。本研究通过分析筛选到的与繁殖力相关的4个关键基因PDS5A、ROCK1、AQR和LTN1,它们在肝脏代谢及奶牛繁殖性能的协同调控方面的功能仍待深入验证。

4 结 论

本研究使用WGCNA、PPI、基因功能富集等生物信息学方法,鉴定得到了不同泌乳时期荷斯坦奶牛肝脏代谢的关键基因,其中泌乳早期12个(RPN1、SEC61A1、SEC61B、SEC61G、SSR1、SSR3、STT3A、DAD1、DDOST、ERLEC1、HM13和OSTC)、泌乳中期及妊娠后期6个(ITGAL、ITGB2、LAPTM5、PTPRC、C3AR1和CTSS),鉴定得到繁殖力相关的荷斯坦奶牛肝脏代谢关键基因4个(PDS5A、ROCK1、AQR和LTN1),为高繁殖力及高产奶牛培育方向的科研工作积累了理论资料。