苏博美利奴羊毛囊发育相关lncRNA与mRNA共表达网络的构建

陈华峰，田可川，黄锡霞，阿布来提·苏来曼，何军敏，田月珍，徐新明，付雪峰，赵冰茹，朱桦，哈尼克孜·吐拉甫

苏博美利奴羊毛囊发育相关lncRNA与mRNA共表达网络的构建

陈华峰1,2，田可川2，黄锡霞1，阿布来提·苏来曼1，何军敏3，田月珍2，徐新明2，付雪峰2，赵冰茹4，朱桦1，哈尼克孜·吐拉甫2

（1新疆农业大学动物科学学院，乌鲁木齐 830052；2新疆畜牧科学院畜牧研究所，乌鲁木齐 830011；3甘肃农业大学动物科学技术学院，兰州 730070；4中国农业大学动物科技学院，北京 100083）

【】随着高通量测序技术的发展和完善，海量的转录组测序数据随之涌现，基因网络的方法也越来越多被用于研究中；细毛羊毛囊发育受多个基因及lncRNA共同调控，单独研究某一分子并不足以发现其调控机制，本研究旨在构建毛囊发育相关lncRNA和mRNA的共表达网络，挖掘细毛羊毛囊发育的潜在候选基因。以苏博美利奴羊胎龄65 d（D65）、85 d（D85）、105 d（D105）、135 d（D135）的胎儿以及出生后7 d（A7）、30 d（A30）的羔羊肩胛部皮肤组织，每个时期有3个生物学重复，共18个样本，进行转录组测序以获得6个不同发育时期的lncRNA与mRNA表达谱数据，筛选出相邻时期的差异表达lncRNA与mRNA，运用加权基因共表达网络分析（weighted gene co-expression network analysis，WGCNA）方法构建共表达模块，使用DAVID在线工具进行GO（gene ontology）和 KEGG pathway富集分析以找到与毛囊发育相关的模块，最后从目标模块筛选高连通度的lncRNA与mRNA使用Cytoscape软件进行网络可视化。从表达谱数据中筛选得到9 070个差异表达的lncRNA与mRNA，运用WGCNA方法得到11个模块，使用DAVID富集分析发现honeydew1、paleturquoise、skyblue2模块中的基因富集在皮肤发育、毛囊发育、毛囊形态发生、Wnt信号通路负调控、细胞粘附等生物过程，以及Wnt信号通路、TGF-β信号通路、Hedgehog信号通路、紧密连接、MAPK信号通路、ECM-受体相互作用等毛囊发育相关的信号通路，并筛选这3个模块中连通度高的lncRNA和mRNA构建了子网络，得到包括、、、、、、、、在内的毛囊发育相关基因，预测出ENSOART00000029117、TCONS_00489976、TCONS_00376759等24个lncRNA的可能靶基因。首次运用WGCNA方法构建了苏博美利奴羊毛囊发育相关lncRNA与mRNA的共表达网络，鉴别出3个毛囊发育相关的共表达模块，发现多个与毛囊发育相关的潜在候选基因，并预测出24个lncRNA可能的靶基因。

苏博美利奴羊；WGCNA；mRNA；lncRNA；毛囊发育

0 引言

【研究意义】随着高通量测序技术的高速发展，尤其是基因芯片技术和第二代测序技术的发展，基因表达数据迅速增长，利用系统生物学网络的方法，对高通量测序技术下的基因表达数据进行分析和挖掘，成为了生物信息学中重要的研究方向[1]。毛囊发育及羊毛品质受多个分子共同调控，单一的个体分子研究不足以发现其机制，因此，利用网络的方法构建与毛囊发育相关的基因网络势在必行。【前人研究进展】苏博美利奴羊是在2014年培育出来的一个超细型细毛羊的新品种，该品种以澳洲美利奴羊公羊为父本，以中国美利奴羊、新吉细毛羊、敖汉细毛羊为母本进行杂交所培育的，具有着羊毛品质好、羊毛产量高、羊毛平均纤维直径为17.0—19.0 μm、羊毛平均长度为8.5—11.0 cm、抗逆性强等特点[2-3]。张立春等[4]利用PCR-SSCP方法研究了苏博美利奴羊和东北细毛羊FGF5基因多态性对毛用性状的影响；付雪峰等[5]研究发现苏博美利奴羊MYL6基因在超细型细毛羊皮肤组织表达量显著高于细型细毛羊；朱桦等[6]发现随着苏博美利奴羊毛囊的逐步成熟，DSG4基因的表达量呈上调趋势。加权基因共表达网络分析（weighted gene co-expression network analysis，WGCNA）是ZHANG等[7]于2005年提出的，LANGFELDER等[8]在2008年编写发布了WGCNA的R包。它是一种用于研究分子作用机制以及网络关系的准确性高的系统生物学方法，目前已被广泛应用于构建基因共表达网络。【本研究切入点】羊毛是纺织工业重要的原材料，研究与羊毛生长息息相关的皮肤毛囊发育相关长链非编码RNA（long noncoding RNA，lncRNA）与mRNA具有重要价值，近年来，关于毛囊的研究也逐渐成为皮肤生物学领域中的热点。毛囊是控制羊毛生长的皮肤附属器官，其形成是由初级毛囊发育至次级毛囊的一个复杂过程，受多个基因联合调控，目前关于毛囊发育的研究多关注个体分子，并不足以发现毛囊发育的分子调控机制，同时关于细毛羊毛囊发育的网络研究尚未见报道。【拟解决的关键问题】本研究采用基于R语言的WGCNA软件包，挖掘苏博美利奴羊18个皮肤样本的转录组表达数据，构建与毛囊发育相关lncRNA与mRNA的共表达模块，对所有模块进行富集分析，找到与毛囊发育相关的模块并构建了共表达网络，并确定这些模块中的核心基因。

1 材料与方法

1.1 数据来源

本研究在课题组前期的转录组测序基础上进行，根据阿布来提·苏来曼等[9]于2015年在新疆科创繁育中心进行的试验，所采集的苏博美利奴羊胎龄65 d（D65）、85 d（D85）、105 d（D105）、135 d（D135）的胎儿以及出生后7 d（A7）、30 d（A30）的羔羊肩胛部皮肤组织，每个时期有3个生物学重复，共18个样本。将采集好的皮肤样品送至上海伯豪生物技术有限公司进行转录组测序（RNA-Seq），包括lncRNA测序和mRNA测序，以获得lncRNA和mRNA的表达谱数据。后续分析于2018年在新疆畜牧科学院畜牧研究所进行。

1.2 筛选差异基因

为了使不同基因和不同样本之间基因的表达水平能够具有可比性，将每个基因的表达水平标准化为FPKM（Fragments Per Kilobase of exon model per Million mapped reads）[10]。将6个发育时期的lncRNA和mRNA的表达谱分别以后一个时期与前一个时期对比，根据各基因的FPKM值，应用edgeR包[11]进行样本间差异基因分析，当q-value≤0.05且Fold-change≥2时，认定该基因为差异表达基因。

1.3 共表达网络的构建

从表达矩阵中提取出差异表达基因的表达值，并使用WGCNA软件包进行共表达网络的分析。使用WGCNA建立检测到的差异表达基因和差异表达的lncRNA的共表达网络，步骤：

计算两两基因间的Pearson相关系数，构建表达相关矩阵，使用pick-SoftThreshold函数计算无尺度拓扑拟合指数选取合适的β值作为构建网络的软阈值，将相关矩阵转换为邻接矩阵，根据邻接矩阵建立无尺度的网络关系。

在所得邻接矩阵的基础上，使用拓扑重叠方法评估两个基因间表达模式的相关性，并构建拓扑重叠矩阵（topological overlap matrix，TOM），根据TOM矩阵计算出基因间相异度矩阵dissTOM=1-TOM。

在相异度矩阵上进行层次聚类以对具有高度相似的共表达关系的基因进行分组，设置模块至少包含30个基因，使用混合动态树剪切法进行分配，得到基因共表达模块，随后计算每个模块的模块特征值（module eigengene，ME），ME表示模块的整体表达水平，计算模块ME间的Pearson相关系数和模块ME间的平均距离，基于该平均距离，使用average-linkage层次聚类对各个模块进行聚类，合并相似度高的模块，确定最终的共表达模块。

1.4 模块的富集分析

使用注释、可视化和综合发现数据库（database for annotation，visualization and integrated discovery，DAVID）[12]分别对各个模块中的基因进行基因本体论（gene ontology，GO）和京都基因和基因组百科全书（kyoto encyclopedia of genes and genomes，KEGG）富集分析，设置＜0.05作为显著性阈值，筛选出与毛囊发育相关的共表达模块，并利用R语言的ggplot2软件包[13]对富集结果进行可视化。

1.5 网络可视化

由WGCNA生成毛囊发育相关共表达模块网络的节点和边属性文件，输出为Cytoscape文件。拓扑重叠度量（topological overlap measurement，Tom）是计算两个基因之间连通性的参数，如果两个基因间具有高的Tom值，则这两个基因与相同的基因集有高度的连接性[14]；基因的连通度是与其相连接的基因的边属性之和，与网络中度的概念类似。筛选Tom＞0.15且由Cytoscape 3.6.1软件[15]保留连通度排序前60内的mRNA和lncRNA，进行可视化和网络模块分析。

2 结果

2.1 差异基因的筛选

通过RNA-Seq测定了lncRNA和mRNA的表达谱，利用edgeR包筛选出了差异表达（differential expression，DE）lncRNA和mRNA，去除各个对比时期重复出现的差异表达lncRNA与mRNA使其保留唯一，得到140个差异表达lncRNA，4 372个差异表达已知mRNA以及4 558个差异表达novel mRNA（表1），用于构建lncRNA与mRNA共表达网络。

2.2 网络模块的构建

2.2.1 软阈值的选择 为了构建符合无尺度分布的网络模块，同时为了降低计算的复杂性并减少计算时间，选择2＞0.85时，β=18作为构建共表达模块的软阈值较为合理（图1），同时对应的平均连通度为70.5。

2.2.2 网络模块的鉴定 在选择加权系数β之后，获得基因之间的相异度矩阵dissTOM，对该相异度矩阵dissTOM进行层次聚类得到基因聚类树（图2中Cluster Dendrogram），使用混合动态树剪切法从该基因聚类树中识别出65个模块（图2中Dynamic Tree Cut），对这65个模块以0.75的相关性合并相似度高的模块后得到11个模块（图2中Merged dynamic）。最后，将9 070个lncRNA和mRNA分成11个共表达网络模块，各个网络模块所包含的基因数及lncRNA数见表2，其中grey模块是无法分配到其他模块的基因集合。

表1 差异表达lncRNA与mRNA数的统计

图1 软阈值的筛选

图2 基因聚类树与模块构建

表2 共表达网络模块中基因数与lncRNA数

2.3 模块的富集分析

为了了解各个共表达模块中基因的生物学功能，对各个模块中的基因进行富集分析，找到与毛囊发育相关的共表达模块。

利用DAVID 6.8在线分析软件分别对11个模块中的基因进行富集分析，以＜0.05作为显著性阈值，在GO富集分析中有发现honeydew1模块中的基因显著富集于皮肤发育过程（GO：0043588）和毛囊发育过程（GO：0001942）等相关的生物学过程（图3 honeydew1），paleturquoise模块中的基因显著富集于毛囊形态发生（GO：0031069）、皮肤发育（GO：0043588）以及毛囊发育（GO：0001942）等相关的生物学过程（图3 paleturquoise），而skyblue2模块中的基因显著富集于Wnt信号通路负调控（GO：0090090）、细胞粘附（GO：0007155）等生物过程（图3 skyblue2）。3个目标模块部分GO富集结果及相关基因见表3。

表3 目标模块部分GO富集

在KEGG富集分析中，发现honeydew1模块中的基因富集在Wnt信号通路（oas04310）、TGF-β信号通路（oas04350）、Hedgehog信号通路（oas04340）、粘着力（oas04510）、ECM-受体相互作用（oas04512）等毛囊发育相关的通路（图4 honeydew1），paleturquoise模块中的基因富集在紧密连接（oas04530）、MAPK信号通路（oas04010）等与毛囊发育相关的通路（图4 paleturquoise），skyblue2模块中的基因富集于Wnt信号通路（oas04310）、ECM-受体相互作用（oas04512）、粘着力（oas04510）等毛囊发育相关的通路（图4 skyblue2）。3个目标模块部分GO富集结果及相关基因见表4。

表4 目标模块部分KEGG富集结果

2.4 网络可视化

富集分析显示在honeydew1、paleturquoise、skyblue2模块中的基因与毛囊发育相关，根据功能关联原则可知这3个模块中的基因与毛囊发育相关，通过WGCNA导出3个毛囊发育相关共表达模块网络的node和边edge文件，通过筛选Tom＞0.15且由Cytoscape软件保留连通度排序前60内的mRNA和lncRNA，进行网络可视化。各模块中连通度排序前5的mRNA与lncRNA见表5。

表5 连通度前5的mRNA与lncRNA

在honeydew1模块中新基因XLOC_170356连通度最高，为该模块中的核心基因，与876个mRNA和lncRNA相连，已知基因中ANGPT1基因连通度最高，而连通度最高的lncRNA是ENSOART00000029117（表5）。通过连通度的高低，筛选出该模块中连通度高的38个mRNA和4个lncRNA构建共表达网络（图5），在该网络中发现一个有趣的基因转化生长因子-β2（Transforming Growth Factor Beta 2，），该基因是honeydew1模块子网络中唯一一个与4个lncRNA均存在着共表达关系的基因，同时该基因显著富集在毛囊发育这一生物过程中（表3）。

点的大小代表连通度高低The size of the dot represents the level of connectivity

在paleturquoise模块中TEC基因共与1 462个mRNA和lncRNA相连接，连通度最高，为该模块中的核心基因，而连通度最高的lncRNA是TCONS_00489976（表5）。筛选出paleturquoise模块中连通度高的19个lncRNA和41个mRNA，进行该模块子网络的可视化（图6），该子网络中的19个lncRNA与41个mRNA均存在着共表达关系。

在skyblue2模块中新基因XLOC_001041共与 1 462个mRNA和lncRNA相连接，连通度最高，为该模块中的核心基因，连通度最高的已知基因为CDH11基因，而连通度最高的lncRNA是TCONS_00041508（表5）。筛选skyblue2模块中连通度高的18个lncRNA和21个mRNA，进行该模块子网络的可视化（图7），发现该子网络中的有3个lncRNA（TCONS_ 00376759、TCONS_00571105、TCONS_00479936）与21个mRNA均存在着共表达关系。

LncRNA Trans作用靶基因预测认为lncRNA的功能与其共表达的基因相关，因此在本研究中lncRNA与mRNA的WGCNA分析同时预测出lncRNA的靶基因。在honeydew1模块子网络中发现ENSOART00000029117和TCONS_00236664的38个靶基因，在paleturquoise模块子网络中鉴定出TCONS_00489976等19个lncRNA的41个靶基因，在skyblue2模块子网络中发现3个lncRNA（TCONS_00376759、TCONS_00571105、TCONS_00479936）的21个靶基因。

3 讨论

基因共表达网络分析是基于系统生物学方法，将表达模式相似的基因进行聚类，以网络模块的形式表现基因的共调控特性，使用模块作为基因分类和功能识别的框架，使用功能关联原则以及模块内研究较为透彻的基因，借助其与其他基因之间的网络关系，推测其他基因的功能，有利于快速鉴定新的目标基因[16]。

点的大小代表连通度高低The size of the dot represents the level of connectivity

点的大小代表连通度高低 The size of the dot represents the level of connectivity

WGCNA是基因共表达网络的常用算法，它是一种高效、准确的生物信息学、生物数据挖掘方法。与基于相关系数的共表达网络算法相比，WGCNA特点是在阈值选择上采用的是软阈值的方法代替之前的硬阈值。之前的基因共表达网络分析方法存在一个很明显的局限，如将阈值定为0.85时，甚至是相关系数为0.84的基因对也将被分类至“无关”的组中，而WGCNA使用基于软阈值的决策方法来避免此问题。

WGCNA兴起于人类医学的相关研究中，而近年来在畜禽研究中也已有所应用，例如：MIAO等[17]为了鉴定与绵羊繁殖力有关的新lncRNA并更好地了解它们的分子机制，在鉴定差异表达mRNA和lncRNA后，利用WGCNA对重要模块中的基因和lncRNA进行功能富集分析，确定了与绵羊繁殖力相关的几个基因模块，以及由hub基因和lncRNA组成的网络，这些网络可能通过调节与TGF-β和催产素信号传导途径有关的目标mRNA来促成绵羊多产性；XU等[18]为了研究在两个绵羊品种的不同胎儿阶段背最长肌有关的基因表达变化的动态，使用时间序列表达分析在各个胎儿阶段鉴定了1 472个差异表达的基因，并利用WGCNA进行分析，识别与胎儿体重显着相关的基因模块，并使用具有高度重要性的基因构建网络，鉴定了，，等与胎儿体重相关的基因；PACHOLEWSKA等[19]使用RNA-seq分析了LPS诱导的马外周血单核细胞中基因表达的差异，并利用加权基因共表达网络分析，识别了10个共同表达的基因模块受LPS体外刺激调节。但是应用WGCNA构建细毛羊毛囊相关的研究尚未见报道。

本研究使用苏博美利奴羊不同发育时期的转录组数据，利用WGCNA方法得到了11个模块，使用DAVID进行GO、KEGG富集分析确定了3个与毛囊发育相关的共表达模块，其中honeydew1模块包含1 853个mRNA和14个lncRNA，paleturquoise模块中包含1 748个mRNA和38个lncRNA，skyblue2模块中含有1 870个mRNA和23个lncRNA，分别从这3个模块中筛选连通度高的基因构建了子网络。

子网络分析发现在3个模块中均存在着与毛囊发育相关的基因。在honeydew1模块中，显著富集在毛囊发育生物学过程上，并且已有研究发现在小鼠中能够诱导静止毛囊干细胞增殖并进入组织再生模式[20]。在paleturquoise模块中，发现41个mRNA与19个lncRNA均存在着共表达关系，其中组织蛋白酶B（cathepsin B，）据报道是小鼠毛囊中潜在的半胱氨酸蛋白酶抑制剂M/E的靶标[21]；14-3-3σ（stratifin，）基因是一种细胞周期调节因子，在小鼠中发现对于毛囊发育是至关重要的，特别是对于毛干的形成[22]，而卜然等的研究表明基因可能参与了细毛羊羊毛生长的调控[23]；丝氨酸肽酶抑制剂，Kunitz 1型（serine peptidase inhibitor, Kunitz type 1，）与肝细胞生长因子激活剂抑制剂1型（hepatocyte growth factor activator inhibitor type 1，）形成一种膜结合的丝氨酸蛋白酶抑制剂，在小鼠中发现/在表皮角质化和毛发发育中具有重要作用[24]；序列相似性家族83成员G（family with sequence similarity 83 member G，）在狗和小鼠的研究中发现参与了毛发形态发生，毛囊分化和毛囊周期[25]；任亮在藏系绵羊中的研究结果表明Gasdermin A（）基因可能通过Wnt信号通路影响毛囊的发育以及周期变化[26]，杨雪的研究表明牦牛GSDMA基因在毛生长周期中参与了毛囊细胞凋亡与增值的调控[27]，而LEI等的研究发现家族成员GSDMA3基因突变可以阻止小鼠毛囊周期中退行期的发生[28]；LEIVA等的研究表明髓磷脂蛋白零样3（myelin protein zero like 3，）在控制皮肤发育，毛发生长和脂肪细胞功能中的复杂作用[29]。在skyblue2模块中，发现18个lncRNA中有3个lncRNA与21个mRNA均存在着共表达关系，其中波形蛋白（vimentin，）在内蒙古绒山羊的研究中发现可能通过影响外根鞘参与毛囊生长周期的调节[30]；Collins等的研究发现细胞维甲酸结合蛋白1（cellular retinoic acid binding protein 1，）在发育形态发生，成人毛囊周期和皮肤肿瘤形成期间被动态调节[31]，表明CRABP1基因与毛囊相关。

子网络分析中找到了部分与毛囊相关的基因，根据功能关联原则，与这些基因存在着共表达关系的基因以及lncRNA也与毛囊发育相关。子网络的可视化结果根据连通度的高低仅显示了排名靠前的部分基因，所筛选的3个模块中还包含着许多已报道的与毛囊发育相关的基因，根据这些已知功能的基因，可以构建相应的子网络，快速找到与之共表达的基因，并推测未知基因的功能。

在lncRNA Trans作用靶基因预测的基本原理是lncRNA的功能与其共表达的蛋白编码基因相关[32]，当lncRNA与一些距离较远的基因在表达量上存在正相关或负相关的情况时，能够通过WGCNA共表达分析来预测lncRNA的靶基因。因此本研究的网络分析中预测出ENSOART00000029117、TCONS_00489976、TCONS_00376759等24个lncRNA的靶基因，为后续进行功能验证提供了参考依据。

本研究构建了苏博美利奴羊毛囊发育相关lncRNA与mRNA的共表达网络，找到了lncRNA与mRNA间的共表达关系，但其具体调控方式还需要进一步深入研究。

4 结论

在本研究中，利用苏博美利奴羊皮肤组织的RNA-seq数据，通过使用加权基因共表达网络分析方法构建共表达网络模块，将筛选出的差异基因划分为11个模块，利用注释、可视化和综合发现数据库进行富集分析发现honeydew1、paleturquoise、skyblue2这3个与毛囊发育相关的模块，鉴定出许多与毛囊发育相关的潜在候选基因及lncRNA，并预测出24个lncRNA可能的靶基因。

[1] 汪涛, 蒋庆华, 彭佳杰, 王亚东. 基因共表达网络的构建及分析方法研究综述. 智能计算机与应用, 2014, 4(6): 47-50.

WANG T, JIANG Q H, PENG J J, WANG Y D. A review of the construction and analysis of gene co-expression network., 2014, 4(6): 47-50.(in Chinese)

[2] ZHANG L C, SUN F L, JIN H G, DALRYMPLE B P, CAO Y, WEI T, VUOCOLO T, ZHANG M X, PIAO Q L, INGHAM A B. A comparison of transcriptomic patterns measured in the skin of Chinese fine and coarse wool sheep breeds.2017, 7(1): 14301-14312.

[3] 柳楠, 田可川, 石国庆, 贺建宁, 刘积凤, 狄江, 杨永林. 苏博美利奴羊核心群级进杂交不同世代对产毛性状影响的研究. 中国畜牧杂志, 2015, 51(15): 6-10.

LIU N, TIAN K C, SHI G Q, HE J N, LIU J F, DI J, YANG Y L. Effect of different generations on wool traits during upgrading cross stages in subo merino nuleus herd.2015, 51(15): 6-10.(in Chinese)

[4] 张立春, 尹峰, 柳俭强, 李梦姝, 王春昕, 曹阳, 朴庆林, 金海国, 张明新. 绵羊FGF5基因Exon 3多态性对毛用性状的影响. 中国草食动物科学, 2018, 38(06): 1-5.

ZHANG L C, YIN F, LIU J Q, LI M S, WANG C X, CAO Y, PIAO Q L, JIN H G, ZHANG M X. Relationship between polymorphisms of FGF5 gene Exon 3 and wool traits., 2018, 38(06): 1-5.(in Chinese)

[5] 付雪峰, 黄锡霞, 刘春洁, 徐新明, 张艳花, 吴伟伟, 于丽娟, 石刚, 田可川. 不同羊毛纤维直径细毛羊MYL6基因克隆及差异表达研究. 中国畜牧杂志, 2016, 52(13): 5-11.

FU X F, HUANG X X, LIU C J, XU X M, ZHANG Y H, WU W W, YU L J, SHI G, TIAN K C. Study on clone and differential expression of MYL6 gene in fine-wool sheep with different fiber diameter.2016, 52(13): 5-11.(in Chinese)

[6] 朱桦, 潘家华, 徐新明, 付雪峰, 田月珍, 阿布来提·苏来曼, 黄锡霞, 田可川. Dsg4基因在苏博美利奴羊毛囊发育不同时期及不同组织中的表达差异. 中国畜牧杂志, 2018, 54(10): 47-51.

ZHU H, PAN J H, XU X M, FU X F, TIAN Y Z, ABLAT S, HUANG X X, CHUANTIAN K. Differential expression of Dsg4 gene in different stages of hair follicles development and different tissues in subo merino sheep., 2018, 54(10): 47-51.(in Chinese)

[7] ZHANG B, HORVATH S. A general framework for weighted gene co-expression network analysis., 2005, 4(1): 1-37.

[8] LANGFELDER P, HORVATH S. WGCNA: an R package for weighted correlation network analysis., 2008, 9(1): 559-571.

[9] 阿布来提·苏来曼, 田月珍, 许琴, 何军敏, 赵冰茹, 徐新明, 付雪峰, 热娜古丽·艾尼瓦尔江, 黄锡霞, 田可川. 苏博美利奴羊胎儿皮肤毛囊结构及形态发育. 中国农业科学, 2017, 50(16): 3226-3235.

ABLAT S, TIAN Y Z, XU Q, HE J M, ZHAO B R, XU X M, FU X F, RANAGUL A, HUANG X X, TIAN K C. Study on fetal skin hair follicle structure and morphological development of Subo Merino sheep., 2017, 50(16): 3226-3235.(in Chinese)

[10] MORTAZAVI A, WILLIAMS B A, MCCUE K, SCHAEFFER L, WOLD B. Mapping and quantifying mammalian transcriptomes by RNA-Seq.2008, 5(7): 621-628.

[11] ROBINSON M D, MCCARTHY D J, SMYTH G K. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data., 2010, 26(1): 139-140.

[12] HUANG D W, SHERMAN B T, LEMPICKI R A. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources.2009, 4(1): 44-57.

[13] 哈德利・威克姆. ggplot2: 数据分析与图形艺术. 西安: 西安交通大学出版社, 2013.

HADLEY W. ggplot2.. Xi'an: Xi'an Jiaotong University Press, 2013.(in Chinese)

[14] CHEN Y C, GUO Y F, HE H, LIN X, WANG X F, ZHOU R, LI W T, PAN D Y, SHEN J, DENG H W. Integrative analysis of genomics and transcriptome data to identify potential functional genes of BMDs in females., 2016, 31(5): 1041-1049.

[15] SHANNON P, MARKIEL A, OZIER O, BALIGA N S, WANG J T, RAMAGE D, AMIN N, SCHWIKOWSKI B, IDEKER T. Cytoscape: a software environment for integrated models of biomolecular interaction networks., 2003, 13(11): 2498-2504.

[16] 魏凯, 张婷婷, 马磊. 猪基因共表达网络模块的构建及功能分析. 畜牧兽医学报, 2017, 48(11): 2205-2215.

WEI K, ZHANG T T, MA L. Construction and functional analysis of gene co-expression network modules.2017, 48(11): 2205-2215.(in Chinese)

[17] MIAO X Y, LUO Q M, ZHAO H J, QIN X Y. Co-expression analysis and identification of fecundity-related long non-coding RNAs in sheep ovaries.2016, 6(1): 39398-39407.

[18] XU L Y, ZHAO F P, REN H X, LI L, LU J, LIU J S, ZHANG S F, LIU G E, SONG J Z, ZHANG L, WEI C H, DU L X. Co-expression analysis of fetal weight-related genes in ovine skeletal muscle during mid and late fetal development stages., 2014, 10(9): 1039-1050.

[19] PACHOLEWSKA A, MARTI E, LEEB T, JAGANNATHAN V, GERBER V. LPS-induced modules of co-expressed genes in equine peripheral blood mononuclear cells., 2017, 18(1): 34-45.

[20] OSHIMORI N, FUCHS E. Paracrine TGF-beta signaling counterbalances BMP-mediated repression in hair follicle stem cell activation., 2012, 10(1): 63-75.

[21] OORTVELD M A W, VAN VLIJMEN WILLEMS I M J J, KERSTEN F F J, CHENG T, VERDOES M, VAN ERP P E J, VERBEEK S, REINHECKEL T, HENDRIKS W J A J, SCHALKWIJK J, ZEEUWEN P L J M. Cathepsin B as a potential cystatin M/E target in the mouse hair follicle.2017, 31(10): 4286-4294.

[22] HAMMOND N L, HEADON D J, DIXON M J. The cell cycle regulator protein 14-3-3sigma is essential for hair follicle integrity and epidermal homeostasis.2012, 132(6): 1543-1553.

[23] 卜然, 李晶, 刘开东, 刘积凤, 柳楠. 细毛羊肩颈部与腹部皮肤组织SFN基因差异表达的研究. 黑龙江畜牧兽医, 2014(19): 80-82.

BU R, LI J, LIU K D, LIU J F, LIU N. Study on the differential expression of the SFN gene in the skin tissues of the shoulder /neck and abdomen in fine-wool sheep.2014(19): 80-82.(in Chinese)

[24] NAGAIKE K, KAWAGUCHI M, TAKEDA N, FUKUSHIMA T, SAWAGUCHI A, KOHAMA K, SETOYAMA M, KATAOKA H. Defect of hepatocyte growth factor activator inhibitor type 1/serine protease inhibitor, Kunitz type 1 (Hai-1/Spint1) leads to ichthyosis- like condition and abnormal hair development in mice., 2008, 173(5): 1464-1475.

[25] BALMER P, FELLAY A K, SAYAR B S, HARITON W V J, WIENER D J, GALICHET A, MULLER E J, ROOSJE P J. FAM83G/ Fam83g genetic variants affect canine and murine hair formation., 2018, 28(4): 13729-13738.

[26] 任亮. 藏系绵羊GSDMA基因皮肤特异表达载体的构建 [D]. 成都: 西南民族大学, 2015.

REN L. Construction of skin-specific expression vector for GSDMA gene of Tibetan sheep [D]. Chengdu: Southwest Minzu University, 2015.(in Chinese)

[27] 杨雪. 牦牛皮肤毛囊周期性结构变化规律及其相关调控因子的研究 [D]. 兰州: 甘肃农业大学, 2017.

YANG X. Study on histological characteristics of hair follicle and expression pattern of related factors in the skin during hair cycle in yak [D]. Lanzhou: Gansu Agricultural University, 2017.(in Chinese)

[28] LEI M X, XIANG G, LI Y, TIAN Y, LIAN X H. Gsdma3 gene is needed for the induction of apoptosis-driven catagen during mouse hair follicle cycle.2011, 136(3): 335-343.

[29] LEIVA A G, CHEN A L, DEVARAJAN P, CHEN Z, DAMANPOUR S, HALL J A, BIANCO A C, LI J, BADIAVAS E V, ZAIAS J, MITEVA M, ROMANELLI P, NOURI K, CAO WIKRAMANAYAKE T. Loss of Mpzl3 function causes various skin abnormalities and greatly reduced adipose depots.2014, 134(7): 1817-1827.

[30] RILE N, LIU Z H, GAO L X, QI J K, ZHAO M, XIE Y C, SU R, ZHANG Y J, WANG R J, LI J, XIAO H M, LI J Q. Expression of vimentin in hair follicle growth cycle of inner Mongolian cashmere goats.2018, 19(1): 38-45.

[31] COLLINS C A, WATT F M. Dynamic regulation of retinoic acid-binding proteins in developing, adult and neoplastic skin reveals roles for beta-catenin and Notch signalling.2008, 324(1): 55-67.

[32] YAN P X, LUO S, LU J Y Y, SHEN X H. Cis- and trans-acting lncRNAs in pluripotency and reprogramming.2017, 46: 170-178.

Construction of Co-expression Network of lncRNA and mRNA Related to Hair Follicle Development of Subo Merino Sheep

CHEN HuaFeng1,2, TIAN KeChuan2, HUANG XiXia1, Ablat Sulayman1, HE JunMin3, TIAN YueZhen2, XU XinMing2, FU XueFeng2, ZHAO BingRu4, ZHU Hua1, Hanikezi Tulafu2

(1College of Animal Science, Xinjiang Agricultural University, Urumqi 830052;2Institute of Animal Science, Xinjiang Academy of Animal Science, Urumqi 830011;3College of Animal Science and Technology, Gansu Agricultural University, Lanzhou 730070;4College of Animal Science and Technology, China Agricultural University, Beijing 100083)

【】With the development and improvement of high-throughput sequencing technology, massive transcriptome sequencing data has emerged, and more and more methods of gene network have been used in research; the hair follicles development of fine wool sheep is regulated by multiple genes and lncRNAs. It is not enough to study a certain molecule alone to discover its regulatory mechanism. The aim of this study was to construct a co-expression network of lncRNAs and mRNAs related to hair follicle development, and to explore potential candidate genes for hair follicle development. 【】The study was performed on fetal skin tissues of the 65th, 85th, 105th, and 135th days and lambs skin tissues of the 7th and 30th days from Subo Merino sheep, with 3 biological replicates in each period, a total of 18 samples were sequenced by transcriptome to obtain lncRNAs and mRNAs expression profiles of 6 different developmental stages, and the differential expressed lncRNAs and mRNAs in adjacent periods was screened to construct a co-expression module by Weighted Gene Co-expression Network Analysis (WGCNA) method, and GO (Gene Ontology) and KEGG pathway enrichment analysis of DAVID online tool was used to find hair follicle related modules. Finally, the high-interconnect lncRNAs and mRNAs were screened from the target module for network visualization using Cytoscape software.【】From the expression profile data, 9070 differentially expressed lncRNAs and mRNAs were screened, and 11 modules were obtained by WGCNA method. The DAVID enrichment analysis revealed that the genes in the honeydew1 module, paleturquoise module, and skyblue2 module were involved in biological processes such as skin development, hair follicle development, hair follicle morphogenesis, negative regulation of Wnt signaling pathway, and cell adhesion, et al; and the signaling pathways involved in hair follicle development such as Wnt signaling pathway, TGF-β signaling pathway, Hedgehog signaling pathway, tight junction, MAPK signaling pathway, ECM-receptor interaction, et al. The sub-networks of lncRNAs and mRNAs with high connectivity were screened in these 3 modules, and the hair follicle development related genes including,,,,,,,andwere obtained, and possible target genes of 24 lncRNAs including ENSOART00000029117, TCONS_00489976, TCONS_00376759 were predicted.【】In this study, we first constructed a co-expression network of lncRNAs and mRNAs related to the hair follicles development of Subo merino sheep using WGCNA method, identified 3 co-expression modules related to hair follicle development, and found several potential candidate genes related to hair follicle development,and predicted possible target genes of 24 lncRNAs.

Subo Merino sheep; WGCNA; mRNA; lncRNA; hair follicle development

10.3864/j.issn.0578-1752.2019.19.016

2018-12-19；

2019-07-03

国家自然科学基金（31460593，31360543）、国家现代农业产业技术体系建设项目（CARS-40）、新疆绒毛用羊遗传育种与繁殖自治区重点实验室项目（XJYS1105）、新疆农业大学研究生科研创新项目（XJAUGRI-2017013）

陈华峰，E-mail：1504159731@qq.com。

黄锡霞，E-mail：au-huangxixia@163.com。通信作者田可川，E-mail：tiankechuan@163.com

（责任编辑 林鉴非）