基于基因调控网络对鸡胚皮肤发育的枢纽基因和关键通路分析

胡斯乐 勿都巴拉 陈宇杰 丽 春 霍红雁 白海花 吴江鸿*

(1.内蒙古民族大学 生命科学与食品学院，内蒙古 通辽 028000；2.内蒙古民族大学 动物科学技术学院，内蒙古 通辽 028000；3.内蒙古民族大学 生物信息学重点实验室，内蒙古 通辽 028000)

皮肤是动物最大的器官，它覆盖了全身。脊椎动物皮肤具有感知、排汗、保护、免疫和屏障等重要作用，这些作用是通过进化出特定的皮肤类型及其附属物来实现[1]。动物皮肤发育过程可分为胚胎期体内形态发生及成体期体外发育2个阶段。鸟类属于卵生动物，其受精卵在母体内发育24 h后大部分组织器官的发生均在体外完成。鸡属于早成鸟，羽毛发育起始于胚胎发育第5～6天[2]，雏鸡出壳时体表已有稠密的绒羽覆盖。鸡羽毛的发育过程完全是在壳内进行的，因此鸡胚是研究脊椎动物皮肤及皮肤附属物发育的理想模型。

禽羽毛的形态发生过程遵循严格发育顺序模式，分别为形成羽基板、毛乳头分化和形成Wnt信号梯度，最后形成区域特异性皮肤附属物[3-4]。此过程有许多信号分子参与并发挥重要作用，如FGF10可激活羽基板的形成[5]，SHH的表达可促进羽芽的发育[6]，BMP2和BMP4可通过抑制羽芽命运来调节羽迹大小[7]。而Wnt/β-catenin信号通路贯穿了羽芽的形成到拓扑排列全过程[3]。一直以来，由于伦理及妊娠期母体影响，胚胎期组织发育机制的研究受到了很大的限制及阻碍。到目前为止，对鸡胚皮肤发育的分子机制相关研究较多，如对早期鸡胚(45 h前)的杂合基因组激活过程及细胞分化不对称等研究[8]，通过对鸡胚发育E6～E10皮肤中JAK-STAT信号通路研究[9]，发现甲硫氨酸可激活Wnt/β-catenin信号通路进而促进鸡胚羽毛发育[10]。但这些研究多数集中在其中某一阶段的分子调控过程上，对整个鸡胚皮肤调控网络变化的研究相对较少，而胚胎期皮肤发育过程是按照严格的顺序模式进行，皮肤及其附属物的形态发生在分子水平上受严格的调控[11]。因此，本研究基于NCBI的SRA数据子库中已有的鸡胚皮肤转录组数据，分析鸡胚连续发育的生长模式，构建与发育阶段对应的基因调控网络，并根据网络拓扑属性及蛋白互作关系深入挖掘了每个调控网络中的关键基因，旨在为鸟类皮肤形态发生以及胚胎发育机制的研究奠定基础。

1 数据来源与研究方法

利用来源于NCBI数据库中的SRA数据子库中胚胎发育第6～21天的皮肤组织RNA-seq数据进行聚类分析，在基因共表达网络、信号通路和蛋白互作网络3个层面基于网络拓扑属性逐步挖掘核心基因集，并构建鸡胚在不同发育阶段的特征网络，整体的分析流程见图1。

图1 数据挖掘分析流程Fig.1 Pipeline of data mining

1.1 数据来源

从NCBI数据库中的SRA数据子库下载鸡胚皮肤转录本表达数据集[12]，该数据是本课题组前期研究结果，前期研究结合组织学试验初步将鸡胚的皮肤及羽毛发育分成了2个较大的形态发生阶段，为了分析核心基因集在形态发生过程中的调控细节。本研究进一步细分鸡胚发育阶段，提取在特定发育阶段的高表达基因，并将基因表达模块与功能模块整合分析。

1.2 研究方法

1.2.1基因表达模式分析

通过Gene Cluster3.0软件进行基因表达模式分析，寻找协同表达基因的共表达模块(Module)，其中对基因(gene)及不同天数(Array)进行聚类，聚类算法采用非加权分组平均法进行。利用Java treeview对分析结果进行可视化分析。每个聚类模块命名为module1～5。同时利用pvclust(R package)对聚类结果进行验证，其中检验方式采用步长检验，检验次数为500次。

1.2.2基因共表达网络的构建及HGecgn基因的筛选

本研究采用Cytoscape3.8.1软件的插件Expression Correlation构建表达趋势相同的共表达网络，利用NetworkAnalyzer插件计算网络节点的3个拓扑属性，分别为最短路径长度(Shortest path length)、介数中心性(Betweenness centrality)和节点度属性(Node degree)。它们从不同层面反映了节点在该网络中的重要程度。通过权重排序评分，计算每个节点的中心性值(图2)。本研究中从基因共表达网络中筛选得到的前25%(TOP25) Hub基因称之为HGecgn基因(Hub gene for expression correlation gene network，基因共表达网络核心基因)。

图2 基于网络属性的节点重要性打分过程 (以节点度(Node)为例)Fig.2 Ranking process of gene importance based on network attribute (node degree is used as an example)

计算中心性值公式如下[13]：

SN(v)=(Cd(v)+Cb(v)+Cspl)/3

(1)

式中：v表示任一个基因，S表示中心性值core的首字母。d是(degree)节点度属性的首字母,spl为(shortest path length)最短路径的首字母，b为(betweenness)介数中心性的首字母。SN(v)为某一节点基因v的中心性值；Cd(v)为该节点的度评分函数；Cb(v)为介数中心性打分函数；Cspl为该节点的最短路径长度打分函数。这3个函数对中心性值的影响都很重要，所以在公式中取了均值。根据中心性值选取排名前25%的节点为HGecgn基因，作为候选基因进一步分析。

1.2.3功能注释分析及通路融合分析

功能注释是利用Cytoscape 3.8.1中的ClueGO插件及DAVID数据库将HGecgn基因富集到KEGG(Kyoto encyclopedia of genes and genomes)信号通路数据库中，其中物种选择鸡，并利用Cytoscape 3.8.1里的CluePedia插件[14]将信号通路与信号通路间，信号通路与基因间的关系进行了可视化分析。发现一部分HGecgn基因参与了多个(Npathway≥2)信号通路，而这些信号通路间在功能上存在明显相关性。因此，将同时参与多个信号通路的候选基因命名为HGkegg基因(Hub gene for KEGG，KEGG信号通路核心基因)，作为后续蛋白互作网络分析的节点基因。

1.2.4蛋白互作网络的构建及基因功能注释分析

利用Cytoscape 3.8.1中的String APP插件对上一步得到的HGkegg基因构建了蛋白互作网络，互作关系中包括了已知互作关系(包括来源于实验数据及数据库)、预测互作关系(包括基因内生性、基因融合和基因协同)及其他蛋白互作关系(文献挖掘、共表达和同源蛋白)等，进一步筛选处于核心位置的具有互作关联性的HGpin基因(Hub gene for protein interaction network，蛋白互作网络核心基因)。最后利用Cytoscape 3.8.1软件展示皮肤发育阶段性相关的特征网络，每个网络中包含信号通路与基因间关系、信号通路与信号通路间关系及蛋白互作关系等信息。

2 结果与分析

2.1 基因表达数据

NCBI数据库中的SRA数据子库中输入Accession序号(PRJNA397795)下载原始数据(Raw sequence reads)，其中包含鸡胚皮肤数据16例的RNA-seq数据。通过edgeR对上述数据库进行预处理后筛选差异基因。表达差异基因筛选条件为：P≤0.05，|log2FC|≥1。在不同天数的鸡胚皮肤组织中，共鉴定得到了5 830个转录本，3 400个差异基因。利用RPKM (ReadsPer Kilobase per Million)算法对转录本的表达量进行标准化。大多数基因的表达水平处于RPKM值1～500，仅有395个转录本具有非常高的表达水平(RPKM>500)。

2.2 基因表达模式分析

聚类分析结果表明，从第6天皮肤开始发育，胚胎期鸡皮肤发育过程大致可以分成5个阶段(图3(a))，第6～9天为第1阶段、第10～12天为第2阶段、第13～14天为第3阶段、第15～17天为第4阶段、第18～21天为第5阶段。从pvclust验证结果上看，虽第6聚类结果有细微出入，但整体分类结果与Gene Cluster分析结果一致，步长检验可信度较高，聚类结果较为可信(图3(b))。

从基因共表达模式分析，5 830个差异表达的转录本对应每个阶段表现出了一种典型模式，将其分别命名为module1～module5(图3(a))。提取了每个module中所包含的基因集，module1至module5中分别包含1 646、380、416、413和1 338个差异表达基因。

2.3 构建基因共表达网络及筛选HGecgn基因

由图4可知，module1和module5的网络节点数量最多，分别为1 260和1 248个。而module2、module3和module4形成的是相对较小的网络，节点数分别为191、250和307个(图4)。从聚类系数中可以看出module1、module4和module5的成簇性较高，聚类系数分别为0.610、0.624和0.688，表明在这3个模块的基因在鸡胚皮肤发育过程调控网络较为复杂，相关性更强。从网络直径看，5个module中直径最短的是module5，直径为9，直径最长的是module4，直径为17。从网络密度上看，module4的网络密度最小，为0.115(表1)。

表1 基因共表达网络网络拓扑属性Table 1 Topological properties of gene co expression network

为探索每个module中发挥关键调控作用的枢纽基因，利用NetworkAnalyzer对基因共表达网络的节点拓扑属性进行分析，并参考每个网络中基因节点的最短路径、介数中心性值和度属性等网络拓扑属性，在本网络中进行排序打分，筛选TOP25的基因定义为HGecgn基因。从筛选结果中可以看出，module1～module5的基因共表达网络中HGecgn基因数量分别为268、34、49、66和254个。

(a)利用Gene Cluster分析鸡胚皮肤基因表达模式；(b)利用pvclust分析鸡胚皮肤基因表达模式。(a) Gene expression pattern of chicken embryo skin by Gene Cluster; (b) Gene expression pattern of chicken embryo skin by pvclust.图3 鸡胚皮肤基因表达模式Fig.3 Gene expression pattern of chicken embryo skin

(a) module1; (b) module2; (c) module3; (d) module4; (e) module5图4 鸡胚皮肤基因共表达网络Fig.4 Co-expression network of chicken embryo skin genes

2.4 基因功能富集分析及HGkegg基因的筛选

从富集分析图5可以看出，module1中的HGecgn基因主要参与细胞的基础代谢，例如糖代谢，如丙酮酸代谢(Pyruvate metabolism)、氨基酸的生物合成(Biosynthesis of amino acids)、糖酵解/糖异生作用(Glycolysis/Gluconeogenesis)、碳代谢(Carbon metabolism)、β-丙氨酸代谢(Beta-Alanine metabolism)、乙醛酸和二羧酸代谢(Glyoxylate and dicarboxylate metabolism)等。除此之外，还包括一些细胞周期发育相关的信号通路，例如细胞周期信号通路(Cell cycle)、DNA复制信号通路(DNA replication)、卵母细胞减数分裂(Oocyte meiosis)等。module2中的HGecgn基因主要参与一些与皮肤及皮肤附属物发育相关的信号通路，例如黏附链接(Adherens junction)、黑色素生成(Melanogenesis)、干细胞多能性调节信号通路(Signaling pathways regulating pluripotency of stem cells)、Wnt信号通路(Wnt signaling pathway)、Hippo信号通路(Hippo signaling pathway)、Rap1信号通路(Rap1 signaling pathway)等。而module3中HGecgn基因主要参与Wnt信号通路，module4中的HGecgn基因主要参与紧密连接信号通路(Tight junction)，module5中的HGecgn基因主要参与ECM受体互作(ECM-receptor interaction)及GnRH信号通路(GnRH signaling pathway)。

为分析不同信号通路的相互作用，将每个模块中同时参与2个KEGG信号通路的基因定义为HGkegg基因。module1中包含的HGkegg基因数为21个，为ARG2、GOT2、SHMT1、PCK2、ACO1、ALDOB、PCNA、MCM6、MCM5、MCM3、YWHAZ、YWHAE、PLK1、CCNB2、DPYS、FBP1、PKM、PC、LDHB、ALDH9A1和GRHPR，它们主要参与了糖类等碳水化合物的代谢及细胞周期两类与基础代谢相关的信号通路。module2中包含的HGkegg基因数为6个，为WNT11、TIAM1、MMP2、LEF1、FGFR1和CTNNB1，它们主要参与了Wnt信号通路及黑色素合成等与皮肤毛囊发育相关的信号传导通路。由于module3、module4和module5中没有筛选到HGkegg基因，因此将基于HGecgn基因的富集结果构建蛋白互作网络。

(a)module1基因共表达网络中HGecgn基因功能富集；(b)module2基因共表达网络中HGecgn基因功能富集。(a) The functional enrichment of HGecgn gene in module1 gene co-expression network; (b) The functional enrichment of HGecgn gene in module2 gene co-expression network.图5 HGecgn基因富集分析Fig.5 Gene enrichment analysis of HGecgn

(a) module1; (b) module2; (c) module3; (d) module4; (e) module5图6 不同module构建的蛋白互作网络Fig.6 Protein interaction networks constructed by different modules

2.5 蛋白互作网络的构建

通过构建蛋白互作网络，发现在module1中HGkegg基因被分成了两个网络,包含了参与的信号通路、关键基因/蛋白及蛋白互作等多个调控关系。第一个网络由参与Pyruvate metabolism、Glycolysis/Gluconeogenesis和Biosynthesis of antibiotics信号通路的PCK2、LDHB和PKM构成，三者存在共表达现象，而PCK2与LDHB间有相似的系统发育谱，结构上同源。另一个网络由参与Cell cycle、DNA replication和Oocyte meiosis信号通路的MCM3、MCM5、MCM6、PCNA、PLK1和CCNB2构成。这6个基因间有相似的系统发育谱，存在共表达趋势。CCNB2在结构上与其他5个基因具有同源区域，因此CCNB2有可能是该网络的枢纽基因。module2网络由参与Adherens junction、Thyroid cancer、Melanogenesis、Signaling pathways regulating pluripotency of stem cells、Wnt signaling pathway、Rap1 signaling pathway的FGF1、LEF1、CTNNB1和WNT11构成。其中CTNNB1同时参与6个信号通路，有可能是该网络中的枢纽基因，且CTNNB1与LEF1存在共表达。module3只富集到了1个信号通路，该网络由参与Wnt signaling pathway的WNT16和TCF7组成，经文献挖掘发现两者间存在基因融合现象。module4网络由组成Tight junction信号通路的MYH1G、ACTN2、MYH1F、MYH1C、MYH1B构成，存在共表达。module5被分成了2个网络。第1个网络由参与ECM-receptor interaction的COMP、ITGB5、LAMB3、ITGB4、LAMC2、CD36和CD44组成，这些蛋白结构上同源，也存在共表达现象。另1个网络由参与GnRH signaling pathway的LOC420419、MAP2K3与MAPK13组成，其中MAP2K3与MAPK13，MAP2K3与LOC420419结构上同源。MAP2K3与MAPK13有相似的系统发育谱，MAPK13与LOC420419，MAP2K3与MAPK13间存在共表达现象(图6)。

3 讨 论

脊椎动物不同区域的皮肤发育模式早在胚胎发育期就已被决定[15]，而鸡胚是脊椎动物胚胎期皮肤发育研究的理想模型。随着测序技术的快速发展，成本的不断降低，公共数据库中可利用的共享数据也与日俱增，如何实现这些数据的联合分析及深度挖掘已成为研究学者们所关注的问题。本研究发现与发育过程实时相关的5个模块。5个模块在基因共表达、蛋白互作和信号通路等不同水平上挖掘了关键的枢纽基因，最后构建了信号通路-基因-信号通路的基因调控网络，为鸟类及其他脊椎动物的皮肤发育模式的研究提供一些新的视角和分析方法。

3.1 鸡胚早期皮肤发育模式

已有研究表明，在早期皮肤发育(第6～9天)中，上皮细胞具有高度的可塑性，可分化成羽毛或鳞片[16]。但随着皮肤的进一步发育，上皮细胞的可塑性逐渐降低，发生区域特化[15]。本研究表明，发育第6～9天的基因表达模式更为相近，有1 646个基因在这一阶段表达量增加，之后逐渐下降。基因共表达调控网络分析结果显示，这些基因间在表达模式上成簇性很高，暗示着基因间紧密相连，相关性很强。对该模式下基因进行功能注释发现，它们主要参与了糖代谢及与细胞周期性发育密切相关的信号通路。HGkegg基因中的CCNB2是细胞周期蛋白B型家族成员之一，为细胞周期蛋白B2。B型细胞周期蛋白包括B1和B2，其中细胞周期蛋白B2可与转录因子βRⅡ相结合，从而在转录生长因子β介导的细胞周期控制中发挥作用。该蛋白在细胞有丝分裂的G2/M期的细胞周期转变中发挥重要作用[17]。微小染色体维持蛋白(Mini-chromosome maintenance proteins MCM)家族由MCM基因编码，高度保守，主要成员包括MCM1-10。其中MCM2-7主要在复制型DNA解螺旋酶活性及真核细胞基因组复制的启动中发挥作用[18]，可以多倍形式复合物发挥作用[19]。研究表明MCM基因与CCNB2协同表达调控细胞周期[20-22]。MCM家族成员与CCNB2在胚胎发育早期协同表达，暗示胚胎发育第九天前鸡胚皮肤处于分化的初级阶段，可塑性强，尚未在分子水平上开启皮肤附属物发育分子调控机制，这与Hughes等[16]研究结果一致。

3.2 鸡胚中期皮肤发育模式

第10～17天的鸡胚发育模式被分成了为连续的3个小网络，模块基因分别在第10～12天、13～14天和15～17天高表达，这些网络基因参与毛囊发育密切相关的信号通路，例如Wnt信号通路[23-24]、黑色素生成[25]和Hippo信号通路等[26]。鸟类胚胎期毛囊的形成要经历诱导、器官发生和分化等3个阶段。毛囊起始阶段始于表皮和真皮成纤维细胞间的相互作用[27]。Wnt信号通路与皮肤毛囊发育密切相关，本研究发现Wnt信号通路中的CTNNB1同时参与6个信号通路，有可能是该网络中的枢纽基因。在毛囊生长早期CTNNB1基因编码的β-catenin蛋白存在于未分化的间充质细胞中，起到毛芽前体形成的第一真皮信号作用，从而诱导上皮细胞聚集，并形成上皮基板[28]。同时结果也显示Wnt11作为候选枢纽基因与CTNNB1具有相互作用，Wnt11是非经典Wnt/PCP信号通路的配体[29]。Gabriela等[30]发现在鸡胚模型中，Wnt11在形成致密的真皮组织及真皮附属物的过程中必不可少，Wnt11缺失小鼠表现出了背部真皮组织发育不全的现象。而在module3中的Wnt16是经典Wnt信号通路配体，可通过经典Wnt/β-catenin信号来促使小鼠皮肤细胞角质化的发生，同时可促进膜间表皮细胞的增殖[31]。因此进一步推测，第10～17天为鸡胚皮肤附属物形成阶段，而Wnt经典信号通路与非经典通路在10～14 d轮流发挥调控作用，为毛囊形态发生关键期，该结果与我们前期的研究结果相一致。在此阶段Wnt信号通路贯穿其中，但不同的Wnt信号通路很可能通过不同的Wnt配体与膜受体来发挥不同的分子调控作用。

3.3 鸡胚晚期皮肤发育模式

鸡胚皮肤发育最后4天(第18～21天)的基因调控网络节点较多，由1 338个节点组成，这些基因主要参与ECM信号通路与GnRH信号通路。细胞外基质(Extracellular matrix, ECM)是由细胞合成并分泌到胞外的一些大分子物质，皮肤ECM可为皮肤细胞生存提供微环境，对皮肤细胞起到支撑保护及提供营养的作用[32]。本研究筛到的ECM受体相互作用网路中的CD44是皮肤中角质化细胞表面透明质酸的主要受体[33]，CD44-/-敲除小鼠表现为角质化细胞表面透明质酸的显著减少，并丢失了角化细胞的增殖、分化及脂类合成功能，导致屏障功能的降低[34]。CD36编码血小板表面糖蛋白，可与长链脂肪酸相结合，促进脂肪酸在质膜上的跨膜运动[35]。LAMB3和LAMC2编码层黏连蛋白，此蛋白属于基底膜蛋白家族之一。层黏连蛋白α亚基、β亚基与γ亚基可共同形成层粘连蛋白5，并与其他细胞外基质相互作用，介导细胞在胚胎发育过程中的细胞黏附、迁移和组织过程[13]。ITGB4和ITGB5共同编码整合素β亚基，与整合素α亚基组成异二聚体，属于非卵形相关跨膜糖蛋白受体。α与β亚基的不同组合方式可与不同的配体相结合，形成有功能的复合物，介导细胞粘附，传递细胞表达信号[36]。整合素α6/β4可作为层粘连蛋白受体，在上皮细胞半桥粒中起到重要作用，并可调节角质细胞极性和运动[37]。COMP编码软骨寡聚基质蛋白，这是一种非胶质ECM蛋白，可与其他细胞外基质蛋白如胶原及纤维连接蛋白相互作用，在组织结构完整性中发挥作用[38]。本研究发现其中CD44、CD36、LAMC2、LAMB3、ITGB4、ITGB5和COMP等节点共同组成了ECM信号通路网络，这也暗示着在出壳前的鸡胚皮肤发生了进一步分化及角化，保证了细胞的功能及结构的完整性，为出壳做准备。

4 结 论

本研究利用前期的研究数据，进一步细分了鸡胚皮肤发育阶段，提取在特定发育阶段的高表达基因，并将基因表达模块与功能模块整合分析。根据网路拓扑属性结合蛋白互作数据库提出了参与不同发育模块的关键枢纽基因。本研究将为脊椎动物的胚胎期皮肤阶段性发育模式的分析提供一种新的分析流程，对转录组测序数据的深入分析提供了一种可选方案。