绒山羊次级毛囊生长调控基因的筛选

李晓燕，张 璐，王乐乐，邬 娇，宋慧子，涂藩镜，张燕军

（内蒙古农业大学动物科学学院，动物遗传育种与繁殖自治区重点实验室，呼和浩特010018）

绒山羊作为中国特色地方畜种，其绒毛具有较高经济价值，是重要出口产品。研究绒毛生长调控机制，为绒产量与绒品质提升提供理论基础，对国民经济的发展具有重要意义。绒毛是毛囊的衍生物，毛囊则是形态和结构较为复杂的皮肤附属结构，由内根鞘、外根鞘、真皮鞘、毛球和毛干几个部分组成［1］。毛囊分为初级毛囊和次级毛囊，初级毛囊形成髓质发达的粗毛，而次级毛囊产生无髓质的绒毛。毛囊密度和次级毛囊与初级毛囊比值（S／P）随羊的品种不同而存在很大差异［2］。绒山羊体表既生长初级毛囊也生长次级毛囊，而奶山羊的体表一般只有粗毛即只生长初级毛囊，不生长次级毛囊［3］。此外，哺乳动物的毛囊呈周期性生长，一般在一个生物年内毛囊要经历生长期、退行期和休止期3个阶段［4］。研究发现，绒山羊次级毛囊生长周期为一年，在4-11月份为生长期，9月份进入兴盛期，12-1月份逐渐进入退行期，2-3月份进入休止期，但周期中各个时期界限并不绝对［5］。近年来，许多学者从分子水平上对次级毛囊的生长发育机理进行探索，筛选出一些与次级毛囊生长发育相关的信号通路，发现基因对皮肤及毛囊发育的调控是一个多方位、多层次、多角度、多基因微效的时空动态过程。研究发现，Wnt信号通路在毛囊的发育和分化中起重要的作用［6］；SH H信号可能参与上皮细胞与间充质细胞之间的信息交换，影响表皮分化和毛囊发育［7-8］；MAPK信号通路能够调控细胞生长、分化、凋亡与炎症反应，EGF、FGF、PDGF、BDNF等为该途径成员，以上成员均被证明参与毛囊生长发育调控［9-13］；TGFβ1、TGF-β2、BMP2、BMP4等均为TGFβ超家族成员，被证明能够影响毛囊周期性生长，调控毛囊生长周期［14-16］；此外，TNF信号途径影响细胞的凋亡，参与毛囊的退行调控［17-18］；Notch通路能够调控毛囊分化且与其他通路之间存在着紧密的联系［19-20］。但是，绒山羊毛囊生长发育与周期性变化的机制尚不明确，需要进一步深入研究。

本研究选取遗传背景相似，以3、6、9、12月份绒山羊和奶山羊皮肤样品为材料，运用RNA-Seq转录组测序技术，结合KEGG等生物信息学方法，分析绒山羊与奶山羊差异表达基因及绒山羊4个时期之间差异表达基因，以期获得调控次级毛囊周期性生长相关调控基因，为次级毛囊生长调控机制的阐明及绒纤维产量和品质的提升提供理论基础。

1 材料与方法

1.1 试验材料

根据绒山羊次级毛囊周期生长变化特点，分别于3月份（休止期）、6月份（生长期）、9月份（兴盛期）、12月份（退行期）的中旬，在4只内蒙古绒山羊成年母羊和4只萨能奶山羊成年母羊体侧部位割取1 cm2左右大小的皮肤样品，-80℃冰箱保存备用。

1.2 研究方法

用TRIZOL进行样品总RNA提取，所得总RNA经1%的琼脂糖凝胶电泳检验其纯度，采用分光光度计（BioRad）检测RNA的浓度和OD260nm／OD280nm的比值。纯化RNA进行mRNA捕获，将捕获的mRNA片段化，加入接头进行反转录再扩增。本试验将同一物种同一时期的4个体混为一个样本混样建库，建好的测序文库采用Illumina Hiseq2000进行双端测序。测序获取原始序列，经fastx-toolkit软件去掉含有两个N的reads，去掉序列中的adaptor，截掉低质量值的reads，剩余的长度仍在16 nt及其以上的reads为clean reads。将有效RNA-seq序列用Tophat2（http：／／goat.kiz.ac.cn／GGD／download.htm）以每个reads容1 nt的错配比对到山羊的参考基因组上，采用RPKM值表征每个基因的表达丰度，以消除mRNA长度和样本测序丰度差异带来的偏差。选用软件edgeR做基因的差异表达分析，认定fold change≥2或≤0.5，P-value≤0.01为差异表达基因。使用WEGO（http：／／wego.genomics.org.cn）与KEGG（http：／／www.genome.jp／kegg／pathway.html）在线平台对这些基因进行功能及pathway分析。不同组织样品的总RNA利用Prime Script RT reagent Kit（Ta KaRa公司）反转录试剂盒进行反转录，用SYBR®Premix Ex TaqTMⅡ试剂盒（TaKaRa公司）在MX3000P荧光定量PCR仪（Agilent公司）上采用SYBR green荧光染料法进行Real time-PCR，以GAPDH基因为内参，采用相对定量2-ΔCt法进行结果分析。所有实时定量试验均进行3次重复，所得数据以“平均值±标准差”表示，Excel绘制基因表达量柱形图，并用SAS7.0软件进行方差分析，P＜ 0.01为差异极显著。利用Primer 5.0软件设计基因的引物序列（表1），送上海生工生物工程有限公司合成。

表1 RT-qPCR引物Table1 List of primers used in RT-qPCR analysis

2 结 果

2.1 绒山羊及奶山羊皮肤转录组组装

测序获取原始数据去掉低质量序列，获得clean reads，本研究获得有效clean reads均在28 435 878以上，为后续研究提供丰富的原始数据。将有效clean reads用Tophat2比对到山羊的参考基因组上，77%～82%的reads都能map上，其中90%～92%的reads为单位置比对，结果见表2。

表2 有效RNA-seq序列及山羊参考基因组定位Table 2 High-quality clean reads and mapping on the goat’s genome

2.2 差异表达基因筛选

采用RPKM值（Reads per kilo bases per million reads）表征每个基因的表达丰度，以消除mRNA长度和样本测序丰度差异带来的偏差。通过差异表达基因（fold change≥2或≤0.5，P-value≤ 0.01）分析发现，绒山羊与奶山羊相同时期共有2 409个显著差异表达基因，其中上调差异表达基因1 175个，下调差异表达基因1 234个，绒山羊不同生长时期显著差异表达基因共550个，其中上调差异表达基因302个，下调差异表达基因248个（表3）。

表3 各样本差异表达基因统计结果Table 3 Total number of differentially expressed genes

2.3 差异基因的GO功能分类

对绒山羊与奶山羊同时期差异基因进行GO注释，从图1可以看出，差异基因根据GO功能注释分为细胞组分、分子功能和生物学过程3大类，在细胞组分中细胞、细胞部分与细胞器部分富集较显著；在生物过程中细胞过程、生物调控以及代谢过程是最大的3个功能结点；而在分子功能中，结合所占比例最大，其次是催化，说明转录组序列与结合和催化相关的功能基因非常丰富。绒山羊不同时期差异基因GO功能注释结果表明，差异表达基因再次显著富集在以上功能，同绒山羊与奶山羊同时期差异基因GO注释结果一致。

2.4 差异基因Pathway分析

利用KEGG平台对DEG进行pathway分析，更深入地了解基因功能及其之间的互相联系，共42个基因富集在毛囊相关WNT、SH H、TGFβ、TNF、MAPK、Notch信号途径中（表4），作为毛囊生长发育相关候选功能基因，后期进行功能验证研究。

2.5 荧光实时定量PCR分析

选取10个富集在WNT信号途径中的差异表达基因WNT 11、WNT 16、FZD10、CTN NB1、VANGL 2、FOSL 1、MMP7、SFRP2、SFRP4和WIF 1进行real-time PCR分析，检测其在绒山羊4个发育时期mRNA表达情况（图2）。与RNA-Seq测序结果比对分析发现，CTNNB1、SFRP4、WNT 16、MMP7、WIF 1、WNT 11、FZD10等7个基因表达模式与测序结果一致，说明测序所产生的差异表达基因生物学重复好，结果准确可靠。

3 讨 论

图1 绒山羊较奶山羊差异基因GO注释Fig.1 Gene ontology of cashmere goat vs dairy goat differential expression genes

表4 富集于毛囊调控通路中的差异基因Table 4 DEG in hair follicle growth related pathway

本次测序共获得38G数据，平均每个样品4.7G，说明测序深度足够，建库质量较好。在转录组数据中，单位置比对多数来自成熟的mRNA，多位置比对则以r RNA和t RNA为主，本研究获得了90%～92%的单位置比对率，奠定了后续分析的准确性。从8个样品中检出的山羊基因总数为20 780，而山羊参考基因组基因总数为21 389，推算覆盖度达到97.15%，说明基因信息代表性充足。绒山羊与奶山羊相同时期共有2 409个显著差异表达基因，其中上调差异表达基因1 175个，下调差异表达基因1 234个，由于奶山羊并不具有次级毛囊，所以推测差异基因可能与次级毛囊生长调控相关。绒山羊不同生长时期显著差异表达基因共550个，鉴于绒山羊毛囊呈周期性生长，推测该550个差异基因可能调控次级毛囊的周期性生长。差异表达基因GO分类结果与许腾发表的绒山羊皮肤毛囊生长期与静息期基因差异表达的研究成果相一致，表明差异表达基因可能通过结合与催化功能，参与生物调控、细胞过程、代谢过程等生物过程，调控毛囊生长，但具体机制有待解析［21］。利用KEGG平台对DEG进行pathway分析，大量基因富集在代谢、疾病、细胞因子和细胞因子受体相互作用等途径中，但以上途径与毛囊生长的关系，有待深入研究。本研究结合前人研究成果，从已报道与毛囊生长相关的WNT、SH H、TGFβ、TNF、MAPK、Notch等信号途径中筛选出42个差异表达基因，推测其为毛囊生长发育相关候选功能基因。

图2 绒山羊不同时期部分WNT通路基因mRNA相对表达量Fig.2 Part of WNT pathway gene relative expression levels in different period of cashmere

WNT信号通路是调控毛囊的起始发育和成纤维细胞的增殖的重要通路，Wnt分子能够与位于细胞膜表面跨膜的Frizzle受体结合，激活胞质dishevelled（Dvl），Dvl是组织细胞内广泛存在的胞浆蛋白，能够阻止胞质内游离的β连环蛋白（β-catenin）降解，β-catenin在胞质内积累到一定量后进入细胞核，与核内淋巴样增强因子及T细胞因子（LEF／TCF）形成复合物，激活下游FOSL 1、MMP7等基因的转录。本研究中有10个差异表达基因富集在该信号通路。分别是WNT 11、WNT 16、FZD10、CTNNB1、VANGL 2、FOSL 1、MMP7、SFRP2、SFRP4和WIF 1。WNT 11、WNT 16是Wnt家族成员，实时定量结果显示较其他基因Wnt11与Wnt16的表达较低。FZD10是Frizzle卷曲蛋白家族成员，9月份有极显著的上调（P＜0.01）与次级毛囊9月份进入生长兴盛期表型相符。β-catenin基因在4个时期中的表达量要高于其他基因，在9月份表达量最高，且较3月份上调极显著，6月份与12月份的表达量也显著高于3月份（P＜0.05），同样符合次级毛囊周期生长规律。然而MMP7表达模式则相反，3月份表达量最高，之后各时期较前一时期表达量均下调。VANGL 2能够指导Dvl在细胞中的定位，其作用模式与Frizzle导致Dvl激活类似，可能通过招募Dvl竞争或促进WNT信号通路［22］，呈现3、6、12月份的表达量均很低，9月份却有显著上调的表达模式。SFRP2及SFRP4均为FRP家族成员，因其结构与Frizzle相似，能够竞争性抑制Wnt因子与Frizzle受体结合，从而影响下游基因激活［23］，SFRP2基因表达量在3月份较12月份有显著上调，符合次级毛囊在3月份处于休止期，生长确实受到抑制的现实。SFRP4表达量的变化同FOSL1一致，在3个时期中的表达量都较低，只有在6月份有显著的上调。WIF 1是Wnt的抑制因子，可以直接与Wnt分子结合，抑制Wnt分子与受体细胞膜上的Frizzled受体及辅助受体LRP5／6形成三聚体复合物，从而抑制了信号传导［24］。以上10个基因虽然在WNT信号通路中具有十分重要的作用，但除β-catenin之外，其余基因在毛囊生长调控中的研究极少。

4 结 论

本研究以有次级毛囊的绒山羊和无次级毛囊的奶山羊为试验动物，利用Illumina Hiseq 2000对绒山羊与奶山羊皮肤样品进行了RNA-Seq测序，共获得38G数据，检出20 780个山羊基因。通过差异表达基因分析发现，绒山羊与奶山羊相同时期共有2 409个差异表达基因，绒山羊不同生长时期差异表达基因共550个。GO功能注释分类显示，差异表达基因主要在细胞与细胞部分、细胞器部分、结合、催化、生物调控、细胞过程、代谢过程等功能富集。KEGG pathway分析发现，42个差异表达基因富集在与毛囊生长相关的WNT、SH H、TGFβ、TNF、MAPK和Notch信号通路中，由此推测其可能参与次级毛囊生长调控过程。

（

）：

［1］ ALONSO L，FUCHS E.The hair cycle［J］.J Cell Sci，2006，119：391-393.

［2］ YANG J B，GAN S Q，YANG Y L.ILK gene cloning and expression of hair follicle growth period of sheep［J］.Genetic，2012，34（6）：179-125.

［3］ ZHU B，XU T，YUAN J L，et al.Transcriptome sequencing reveals differences between primary and secondary hair follicle-derived dermal papilla cells of the Cashmere goat（Capra hircus）［J］.PLoS ONE，2013，8（9）：e76282.

［4］ 李金泉，尹 俊，周欢敏.绒山羊绒毛生长调控机理研究进展［J］.中国草食动物，2001，3（3）：41-44.LI J Q，YIN J，ZHOU H M.The research progress of cashmere goat growth regulation mechanism［J］.Herbivores of China，2001，3（3）：41-44.（in Chinese）

［5］ 张春兰.内蒙古阿尔巴斯白绒山羊毛囊生长周期性变化与KAP6 cDNA的特征分析［D］.呼和浩特：内蒙古农业大学，2003 ZHANG C L.Inner Mongolia Albath White Cashmere goat hair follicle growth periodic and KAP6 cDNA characteristics analysis［D］.Hohhot：Inner Mongolia Agricultural University，2003.（in Chinese）

［6］ ENSHELL SEIJFFERS D，LINDON C，KASHIWAG M，et al.beta-catenin activity in the dermal papilla regulates morphogenesis and regeneration of hair［J］.Dev Cell，2010，18（4）：633-642.

［7］ WOO W M，ZHEN H，ORO A.Shh maintains dermal papilla identity and hair morphogenesis via a Noggin-Shh regulatory loop［J］.Genes Dev，2012（26）：1235-1246.

［8］ TEICHERT A E，ELALIEH H，ELIAS P M，et al.Overexpression of Hedgehog signaling is associated with epidermal tumor formation invitamin d receptornull mice［J］.Sci Int Dermatol，2011：2289-2297.

［9］ MOU C Y，JACKSON B，SCHNEIDER P，et al.Generation of the primary hair follicle pattern［J］.PNAS，2006（24）：9075-9080.

［10］ UEKI M K，YUKO O，JUNKO O，et al.Hair cycle resting phase is regulated by cyclic epithelial FGF18 signaling［J］.J Int Dermatol，2012，132（5）：1338-1345.

［11］ HUH S H，NÄRHI K，LINDFORS P H，et al，Fgf20 governs formation of primary and secondary dermal condensations in developing hair follicles［J］.Genes Dev，2013，27（4）：450-458.

［12］ KARLSSON L，BONDJERS C，BETSHOLTZ C.Roles for PDGF-A and sonic hedgehog in development of mesenchymal components of the hair follicle［J］.Development，1999，126：2611-2621

［13］ ULLRICH R，PAUS R.Molecular principles of hair follicle induction and morphogenesis［J］.Biol Essays，2005，27（3）：247-261.

［14］ SOMA T，OGO M.Analysis of apoptotic cell death in human hair follicles in vivo and in vitro［J］.J Int Dermatol，1998，111：948-954.

［15］ OSHIMORI1 N，FUCHS E.Paracrine TGF-βsignaling counterbalances BMP-mediated repression in hair follicle stem cell activation［J］，Cell Stem Cell，2012（10）：63-75.

［16］ MYLLYLA R M，HAAPASAARI K M，PALATSI R，et al.Multiple miliary osteoma cutis is adistinct disease entity：four case reports and review of the literature［J］.Br J Dermatol，2011，164：544-552.

［17］ MIKKOLA M L.TNF superfamily in skin appendage development［J］.Cytokine Growth Factor Rev，2008（19）：219-230.

［18］ LEI M X，BAI X F，YANG T，et al.Gsdma3 is a new factor needed for TNF-αmediated apoptosis signal pathway in mouse skin keratinocytes［J］.Histochem Cell Biol，2012，138：385-396.

［19］ CROWE R，HENRIQUE D，HOROWICZ D，et al.A new role for Notch and Delta in cell fate decisions pat-terning the feather array［J］.Development，1998，125（4）：767-775.

［20］ ZHANG Y，YU J，SHI C Y，et al.Lef1 contributes to the differentiation of bulge stem cells by nuclear translocation and cross-talk with the notch signaling pathway［J］.Int J Med Sci，2013 10（6）：738-746.

［21］ 许 腾.绒山羊皮肤绒毛生长期与静息期基因差异表达的研究［D］.呼和浩特：内蒙古大学，2013.XU T.Research of differentially expressed genes of cashmere goat hair follicles anagen and telogen［D］.Hohhot：Inner Mongolia University，2013.（in Chinese）

［22］ 王 娟，杨雪艳，王红艳，等.Vangl2及其相关基因与神经管畸形的关系［J］.生命的化学，2010，30（5）：730-733.WANG J，YANG X Y，WANG H Y，et al.Progress of Vangl2 and related genes in NTD［J］.Chemistry of Life，2010，30（5）：730-733.（in Chinese）

［23］ 罗丽敏，王 俐.Wnt通路抑制因子Sfrp的研究进展［J］.医学综述，2007，13（10）：729-731.LUO L M，WANG L.The research progress of Wnt pathways inhibitory factor Sfrp［J］.Medical Recapitulate，2007，13（10）：729-731.（in Chinese）

［24］ 胡毓安，赵春杰.Wif1在神经系统发育中的研究进展［J］.浙江大学学报（医学版），2010，39（1）：93-96.HU Y A，ZHAO C J.The research progress of Wif1 in nervous system development［J］.Journal of Zhejiang University（Medical Sciences），2010，39（1）：93-96.（in Chinese）