鸡胚卵巢基因共表达网络模块的构建及功能分析

任一帆，孟 婕，杨 晴，张 贝，李可强，张传生*，耿立英，赵雪艳，李祥龙

(1 河北科技师范学院 河北省特色动物种质资源挖掘与创新重点实验室，河北 秦皇岛，066600；2 山东省农业科学院畜牧兽医研究所)

卵巢作为鸡繁殖活动的基础，近年来受到研究人员的广泛关注，特别是对鸡胚卵巢转录组的研究。例如：曹顶国[1，2]对高、低产汶上芦花鸡的卵巢组织转录组分析，分别获得有功能注释的基因4 675个和4 472个，并以低产汶上芦花鸡为对照，在高产个体的卵巢组织中获得了158个差异表达基因。前人在鸡胚卵巢后期生长发育的转录组研究取得了阶段性研究成果，但对鸡胚胎期的卵巢发育研究还很匮乏，尤其缺乏利用系统生物学网络分析的方法对鸡胚胎期卵巢组织转录组基因的研究。此外，鸡卵巢发育受多个分子共同调控，单一的个体分子研究不足以发现其机制，而加权基因共表达网络分析作为一种系统生物学方法，目前在多个领域被广泛应用[3～6]。通过WGCNA可系统的分析基因之间的相互关系，发现共表达基因的调控中心，并将基因模块与性状相联系，从而找出目标性状的关键基因[7～10]。为此，笔者以不同发育阶段的转录组测定的卵巢组织鸡的RNA-Seq数据为基础，利用WGCNA分析鸡卵巢发育过程中的基因表达谱，挖掘出与鸡卵巢发育高度相关的基因，并试图筛选出关键的调控基因，为鸡卵巢发育具体机制的研究提供帮助。

1 材料与方法

1.1 数据获取

鸡卵巢在不同发育阶段的转录组测序数据来自文献[11]，样品的RNA-Seq数据在NCBI的登记号码为：SRP081829。该数据以白来航鸡为材料，分别选取孵化第6天(E6)、E12和出壳第1天(D1)的左右侧性腺，进行高通量测序，每组包括2个生物学重复。所有基因的表达量通过FPKM标准化，并去除标准差小于5的基因。

1.2 加权基因共表达网络构建

通过R软件(R version 3.6.2)的WGCNA包(版本1.6.9)提供的函数构建加权基因共表达网络[12，13]。为了尽可能的满足无尺度网络分布要求，调用WGCNA中的函数pick Soft Threshold，设置合适的软阈值β，并采用WGCNA自带的动态切割法划分基因模块，设定模块的最小基因总数为100，切割高度为0.85，其他参数参考函数的默认值。

1.3 组织特异性模块鉴定

计算出模块基因的模块特征向量(module eigengene, ME)，模块的特征向量与不同组织性状之间的相关性系数r>0时表示模块表达量与性状呈正相关[14]。绘出模块与性状的相关性热图，本次研究定义相关系数r>0.65且显著性P<0.01的模块为组织特异性模块。

1.4 模块的富集分析

通过R软件提供的clusterProfiler，enrichplot等软件包，对所有8个模块进行GO(gene ontology，GO)和KEGG (Kyoto Encyclopedia of Genes and Genomes)分析。GO分析主要是将每个模块的基因分为生物过程(Biological Process，BP)、分子功能(Molecular Function，MF)以及细胞组分(Cellular Component，CC)等3个部分，经过GO分析可以明确每个模块基因具有的生物学作用。KEGG分析通过信号通路的显著性富集，确定差异表达基因参与的主要生化代谢途径和信号转导途径。设置阈值p<0.1作为显著性阈值，筛选出与鸡卵巢发育相关的基因，并利用R语言的ggplot2软件包[15]对富集结果进行可视化。

1.5 hub gene的筛选

根据模块中各基因的模块连通性和临床性状关系选择候选hub基因。模块连通性定义为基因间皮尔逊相关性的绝对值(模块成员，即MM)，临床性状关系定义为各基因与性状的皮尔逊相关绝对值(基因显著性,即GS)。设置候选hub基因的MM>0.8，GS>0.6。同时，将hub模块中的候选基因，使用搜索工具For the Retrieve of Interactive gene(String；https://string-db.org/)数据库构建PPI网络，对基因的功能进行注释。将候选基因导入Cytoscape筛选连接值前30的基因，并呈现网络进行下一步分析。

2 结果与分析

2.1 筛选表达谱基因和样本

将所有的基因通过FPKM标准化并去除低表达量的数值之后，得到了5 557个表达基因的FPKM值。通过对样本进行聚类分析，结果表明，12个样本分为2组，一组包含8个样本，另一组为4个样本，本次研究数据的样本聚集度相对较好，没有进行离群值的删除(图1)。

图1 鸡胚卵巢基因样本聚类分析

2.2 软阈值的确定

在确定无尺度网络拟合指数R2>0.85时，软阈值β确定为16比较合适，在β=16时所对应的平均连通度为257，最大的连接系数为882，中位数为171(图2)。

注：图中的横轴均代表软阈值(β)。 A: 纵坐标对应的是无尺度网络模型指数; B: 每一个软阈值对应的网络平均连接度。 图2 鸡胚卵巢加权基因共表达网络分析软阈值(β)的确定

2.3 基因共表达网络构建

在软阈值为16的条件下，采用动态剪切法构建加权基因共表达网络，共得到8个模块(图3)，每个模块都代表一组相关度比较高的基因，灰色模块代表不能被纳入其他模块的基因。各个模块都用不同的颜色表示，各模块包含的基因数量相差较大，最少的是黑色模块，包含109个基因；最多的是青色模块，包含2 118个基因，其余6个模块包含的基因数分别为：534(黄色),602(绿色),854(黑色)和143(红色)。

注：基因聚类树的每一个颜色对应一个模块分支1～7代表不同的模块，1：黑色模块，2：蓝色模块，3：红色模块，4：绿色模块，5：黄色模块，6：棕色模块，7：青色模块，灰色代表灰色模块代表不能被纳入其他模块的基因图3 基于拓扑重叠构建的鸡胚卵巢基因聚类树

注：横轴表示不同的性状, 纵轴表示每一个模块的特征向量。红色的格子代表性状与模块具有正相关性, 绿色的格子代表性状与模块具有负相关性。矩形框里的数字代表模块与性状之间的相关系数及相应p值。(E6L：孵化6天鸡胚左侧卵巢；E6R：孵化6天鸡胚右侧卵巢；E12L：孵化12天鸡胚左侧卵巢；E12R孵化12天鸡胚右侧卵巢；D1L：出壳后第1天雏鸡左侧卵巢；D1R：出壳后第1天雏鸡右侧卵巢)图 4 鸡胚卵巢基因共表达网络模块与性状相关性

2.4 组织特异性模块鉴定

在8个模块中，除去灰色模块共有6个模块与组织存在高度特异性(图4)，分别为：黑色(r=0.76，P=4×10-3)、蓝色(r=0.95，P=2×10-6)、棕色(r=0.97，P=3×10-7)、青色(r=0.88，P=2×10-4)、红色(r=0.81，P=10-3)与黄色(r=0.96，P=10-6)。大多数模块都有与其高度相关联的组织,例如棕色模块与D1L卵巢组织存在高度相关性,并且连接度也最高。而绿色模块没有与任何一个组织存在高度相关性，通过筛选最终得到6个模块进行下一步的富集分析以及目标模块的鉴定。

2.5 鸡卵巢发育相关模块的筛选

为进一步了解各模块的生物学功能，本次研究对筛选的6个模块的KEGG通路和GO富集分析进行深度解析。

在孵化期的第6天，鸡卵巢特征显现，只有右侧卵巢组织有与之显著相关的模块：黑色，该模块富集到(GO:0002443)白细胞介导的免疫反应以及跟免疫相关的94条GO 条目，并富集到糖异生，丙酮酸代谢等9条KEGG通路。在孵化进行到第12天时，与左侧卵巢组织显著相关的青色模块可以富集的GO条目有：(GO:0008585)雌性性腺发育，(GO:0061458)生殖系统发育，(GO:0030111) Wnt信号通路的调控等共328条，而KEGG富集的通路有：细胞周期，卵母细胞减数分裂，p53信号通路，TGF-β信号通路等共25条；与右侧卵巢高度相关的模块有2个,分别是红色模块和黄色模块，红色模块显著富集到：(GO:0006342)染色质沉默，(GO:0016458)基因沉默，(GO:0004857)酶抑制剂的活动等122条GO条目，富集到的KEGG通路有：氨基酸的生物合成等共17条。黄色模块共富集到 151 个 GO 条目, 其中多个GO条目与体液免疫反应以及防御反应相关，KEGG共富集到39条通路，大多是酪氨酸代谢和精氨酸生物合成等与氨基酸代谢相关的通路。

出壳第1天，左侧卵巢组织与棕色模块存在最高的相关性，棕色模块富集到(GO:0001838)胚胎上皮管形成等166条GO条目，并且显著富集到孕酮介导的卵母细胞成熟和GnRH信号通路等9条KEGG通路。而右侧卵巢组织与蓝色模块相关性最高，蓝色模块主要富集于(GO:0055085)跨膜转运等248条GO条目，以及氨基酸合成与代谢等21条KEGG信号通路。

虽然6个模块各自都有相关的性状，但根据富集分析结果发现，与鸡卵巢发育相关的通路主要集中在青色模块，并且该模块富集到本试验领域已经报道的Wnt通路。右侧卵巢的发育停滞以及萎缩主要与红色模块的基因相关，所以选择青色模块与红色模块作为特征性模块进行下一步的分析。青色模块与红色模块的GO富集分析气泡见图5。

图5 鸡胚卵巢基因共表达网络模块中青色模块与红色模块的GO富集分析

2.6 枢纽基因的筛选

模块内基因的连通性(Connectivity)代表该基因与其他基因之间的调控关系，连通性越高，说明该基因在模块中的调控作用越大，而枢纽基因的连通性通常情况下是在该模块中连通性较大的基因。参考WGCNA 的标准，对特征性进行枢纽基因筛选，确定了2个模块的候选枢纽基因数量，分别为：1 863(青色)，442(红色)。每个模块MM值排名前10的候选枢纽基因见附表1。将候选枢纽基因利用STRING计算出其互作关系，输出为Cytoscape文件，并借助STRING网站注释了基因的功能(附表2)。

附表1 鸡胚卵巢基因共表达网络模块中的6个模块前10个枢纽基因

附表2 鸡胚卵巢基因共表达网络模块候选枢纽基因的功能注释

将特征性模块的候选基因利用Cytoscape软件的cytoHubba插件根据程度值筛选出连接度前30的基因作为中心网络，候选基因的连线分别为：青色模块205条，红色模块328条。

使用Cytoscape将结果可视化，青色模块和红色模块的中心网络见图6。

注：在网络中的每一个基因代表一个节点，颜色由深到浅代表连接度从大到小。A:青绿色模块的中心网络，B：红色模块的中心网络。图6 鸡胚卵巢基因共表达网络模块中青色模块和红色模块的中心网络结构

青色模块的核心基因CDC5L与红色模块的核心基因PLG在卵巢的表达量变化见图7。可以看出CDC5L基因与PLG基因都在胚胎左侧卵巢发育的第12天表达量激增，不同的是CDC5L在左侧卵巢表达量增加，PLG基因在右侧卵巢表达量增加。此外，PLG基因在胚胎期第12天右侧卵巢的表达量远大于其他时期，但CDC5L基因在同时期卵巢左侧的表达量与出壳第1天的数值相差不大。

A:青色模块的核心基因表达量，B：红色模块的核心基因表达量图7 鸡胚卵巢基因共表达网络模块中青色模块和红色模块的核心基因表达量

3 讨 论

3.1 WGCNA 分析的重要意义

传统的生物学研究侧重于从分子水平阐释单个功能元件(如DNA, mRNA和蛋白质)对生命活动的影响,只能局部地解释某一生命活动发生的原因。而新兴的共表达网络分析方法可以将复杂的数据进行归纳和整理，高效研究基因整体表达规律。相比于其他调控网络,WGCNA可以特异地筛选出与性状相关的基因，进行模块化分类,得到具有高度生物学意义的共表达模块,已经被证明是一种高效的数据挖掘手段[16]。

3.2 鸡胚卵巢发育相关共表达模块的鉴定

本次研究根据鸡卵巢不同发育时期的转录组数据，利用 WGCNA 方法进行分析，最终获得了8个共表达模块，通过对6个组织特异性模块进行富集分析,发现这些模块均可以得到具有生物学意义的调控通路。例如：青色模块富集到(GO:0007049)细胞周期等细胞过程，蓝色模块可富集到(GO:0015267)通道活性，(GO:0055085)跨膜转运等通路。富集结果提示了参与鸡卵巢发育生理活动的差异基因除了集中于细胞内，还广泛分布于细胞膜上，该结论与兰道亮等[17]对牦牛卵巢转录组研究获得的结果相一致，这说明细胞膜组分在卵巢生理生殖活动中可能也扮演了重要角色。

3.3 鸡胚卵巢发育节点特征

鸡胚的原始生殖腺位于中肾嵴腹内侧，由中胚层发育而来，在孵化期第3.5天左右出现[18]；性腺分化在孵化期第5.5天出现，到孵化期第6天时通过形态学方法能够分辨出性别[19]，性腺分化在孵化期第7天时差异更显著；原始生殖细胞在孵化期第11天左右完全分化[20]。卵巢皮质部在孵化期第11～12天时开始出现原始卵泡，并且数量不断增加。卵原细胞(oogonia)减数分裂生成初级卵母细胞出现在孵化期第15.5天[21，22]。到了孵化期第16～18天，雌性左侧卵巢继续发育，外观似腔卵泡样结构；而右侧性腺退化，外观呈睾丸样结构。雌鸡性腺一般情况下为左侧正常发育，右侧退化并逐渐萎缩[23，24]。但有学者研究发现，在雌鸡性腺发育过程中会出现双侧功能卵巢闭锁，并证实该现象并非巧合，原因之一是等位基因的组合配对导致出现这样的表型[25]。

3.4 鸡胚卵巢发育基因的筛选

本次研究根据WGCNA枢纽基因的筛选标准：GS>0.6，MM>0.8，并结合富集分析结果，筛选出2个特定模块的候选枢纽基因；并将候选基因利用在线网络工具：STRING，根据已有的研究数据建立基因之间的连通性，最后将基因间的连通性信息导入到 Cytoscape 软件中，筛选每个模块连接度前30的基因，并根据STRING对每个基因的注释，进一步识别出每个模块的枢纽基因。

在孵化期第12天，卵巢显著相关的模块所参与的生物调控通路差别较大。左侧卵巢与青色模块高度相关，参照富集分析结果，该模块主要调控卵母细胞分裂并调控卵巢的发育；模块的核心基因CDC5L基因大量存在于细胞核，在mRNA剪接前的第2个催化步骤起重要作用；CDC5L在癌症以及肿瘤方面的研究成果较多，研究表明：CDC5L对于猪卵母细胞减数分裂和早期胚胎发育至关重要，但是该基因在鸡胚发育方面研究较少，本次试验结果可能提示CDC5L在鸡胚胎卵巢发育早期起到重要作用。右侧卵巢与红色模块显著相关，富集分析显示红色模块主要调控基因以及染色质的沉默，并参与酶抑制剂的合成活动，右侧性腺在该时期处于停滞发育并出现萎缩的时期，富集分析结果与该时期节点特征相一致。红色模块的核心基因PLG在胚胎右侧卵巢发育的第12天表达量激增，该基因编码蛋白为纤维溶酶原，在包括胚胎发育，组织重塑在内的多种过程中充当蛋白水解因子，但目前在鸡胚卵巢发育方面对PLG的研究资料很少，本次试验结果可能提示该蛋白在胚胎右侧卵巢发育停滞时期起到关键作用。此外，PLG大量存在于细胞外，这个结果提示核心基因不仅在细胞内，在细胞外也同样存在。

4 结 论

利用权重基因共表达网络分析(WGCNA)的方法将5 558个差异基因划分成8个共表达模块，将6个有意义的模块完成了模块内基因的功能富集分析，揭示了模块所蕴含的生物学意义，得到不同发育阶段卵巢组织差异表达基因的功能和代谢通路，并筛选出2个特征模块的枢纽基因：CDC5L和PLG。