基于高通量转录组测序阿尔巴斯型绒山羊绒毛生长周期差异表达基因分析

赵 濛 刘志红* 奈日乐 谢遇春 马丽娜 米 璐 李金泉，2 李 婕

(1.内蒙古农业大学 动物科学学院，呼和浩特 010018；2.内蒙古自治区动物遗传育种与繁殖重点实验室，呼和浩特 010018；3.东北农业大学 动物科学技术学院，哈尔滨 150030)

山羊绒是重要的毛纺织原料，具有纤细、柔软、轻薄等优点。内蒙古阿尔巴斯型绒山羊是我国重要的绒肉兼用型绒山羊品种，主要分布于内蒙古自治区伊克昭盟鄂托克旗境内，是经过长期自然选择和人工选育而成的优良品种，该品种可适应西北地区寒冷干旱和多风的气候环境。绒毛是来源于绒山羊下层毛被的一种优良的纺织原料，绒山羊毛被具有2种截然不同的纤维结构，位于毛被上层的是粗且硬的硬质粗毛，位于毛被下层的是细且软的绒毛。绒毛纤维来源于生于皮肤中的次级毛囊结构，粗毛来源于生于皮肤中的初级毛囊结构[1-3]。

出生后的哺乳动物皮肤毛囊生长处在动态变化当中，其毛囊生长变化具有往复循环规律性，也就是出生后毛囊生长具有周期性，其1个生长周期大致可分为生长期、退行期和休止期[4-6]。近年来，科学研究针对控制毛囊周期生长规律的分子调控机制进行了大量的探索，发现某些信号通路在毛囊周期生长中具有重要的作用[7-11]。毛囊生长周期内不同阶段的转换是由于不同的信号通路间相互促进以及抑制作用的结果。规律性的调控机制作用导致了毛囊生长的规律性以及生长周期间的转换。绒山羊皮肤次级毛囊同样具有类似的生长周期性，并且,1个生长周期具有生长期、退行期以及休止期之分，但是针对绒山羊次级毛囊生长周期性的分子调控机制研究甚少。所以，从分子生物学角度研究绒毛周期生长中不同阶段的基因表达情况以及不同阶段的基因表达差异性对了解调控其生长分子机制具有十分重要的作用。

Jiang等[12]通过构建抑制消减杂交文库的方法对绒山羊绒毛生长的生长期和休止期皮肤样本的差异表达基因进行研究发现，生长期和休止期绒山羊皮肤具有差异表达基因45个，其中23个具有明确的功能。但这种方法对研究差异基因具有随机挑选的特性，并不能对基因表达的全谱进行深入研究。随着研究技术的进步，第二代高通量测序技术应用于差异基因的筛选研究方面，使得对基因表达全谱的研究成为可能。RNA-Seq技术是一种高通量的测序技术，可以针对基因转录本进行深入研究，发掘低丰度转录本和新转录本以及鉴定不同样本间的差异表达转录本[13-15]。本试验采用第二代高通量测序技术针对绒毛生长不同阶段的皮肤样本进行转录组测序，旨在研究绒毛生长不同阶段的皮肤差异表达基因以及基因差异表达对调控绒毛周期转换的相关性，这些绒毛生长不同阶段的差异基因所具有的生物学功能对深入了解绒毛生长的调控机理具有重要的作用，并对研究绒毛生长调控机理奠定理论基础。

1 材料与方法

1.1 试验动物及皮肤组织样品采集

阿尔巴斯型内蒙古绒山羊是绒肉兼用的优良品种，本试验选取内蒙古伊克昭盟阿尔巴斯绒山羊种羊场的阿尔巴斯型内蒙古绒山羊3只，试验用羊选取条件为：选自同一放牧环境下、体重体尺相当、周龄相同、无亲缘关系、生长状况良好的周岁母羊，采样周期为期1年，起于2011年12月，止于2012年12月，采样周期内试验用羊的饲养管理制度与当地饲养管理一致，且试验用羊均未产羔。皮肤样本采集时间为每月初。皮肤样本采集部位为体侧体中线部距离肩胛骨10～15 cm处，采取手术直接剥离大小为3 cm直径的皮肤样本，皮肤样本手术剥离时间控制在5 min之内，采集后样本进行PBS冲洗后快速液氮冷冻，然后放入液氮罐带回实验室存储于-80 ℃ 冰箱中备用。

1.2 试验方法

1.2.1皮肤总RNA提取

采用RNAiso Plus 试剂盒(Trizol法)操作说明分别提取3只绒山羊体侧皮肤组织的总RNA，无菌无酶水溶解，分别使用微量紫外分光光度计以及2100 bioanalyzer对所提取皮肤总RNA进行浓度、纯度以及完整度的检测。总RNA置于-80 ℃冰箱保存。

1.2.2测序文库构建

根据RNA样本浓度，分别将3只绒山羊皮肤的总RNA等量混池，为构建12个月的cDNA文库做准备。转录组测序cDNA文库构建参考Illumina公司(美国)的TruSeqTMRNA样本制备试剂盒操作说明进行。利用oligo dT磁珠法将mRNA从总RNA中富集纯化出来，利用Illumina特有的片段化缓冲液将mRNA片段化成为大小在100～400 bp之间的小片段mRNA，利用Invitrogen公司的cDNA合成试剂盒以片段化的mRNA为模板合成双链cDNA，通过核酸外切酶和聚合酶的作用将双链cDNA片段的末端补平，并将双链cDNA磷酸化连接测序接头和加Poly A尾，利用Bio-Rad公司的certified low range ultra agarose试剂盒回收大小在200～300 bp的cDNA，利用Phusion DNA聚合酶(NEB)对cDNA进行PCR扩增，PCR循环次数为15次，获得测序文库，随后使用TBS380对文库cDNA进行质检。

1.2.3转录组测序

使用Illumina HiSeqTM2000测序平台对cDNA进行双末端测序，进行2×100 bp 测序试验，由上海美吉生物医药科技有限公司测序。

1.2.4测序数据除杂、拼接及拼接结果注释

测序得到的原始fastq数据，运用SeqPrep和condetri_v2.0.pl软件进行过滤除杂，去除测序质量较低的reads(测序质量值小于25)，去除reads中的adapter序列，去除含有较多N的reads，去除小于25 nt的短序列，得到高质量的序列数据。对除杂过滤后的序列运用Trinity软件进行无参考基因组信息de novo拼接，并将组装结果利用Trinity软件进行ORF预测，并运用BLAST(Version 2.2.25)软件将预测出ORF的蛋白序列与未预测出ORF的序列分别注释，将预测出的ORF的蛋白序列分别于NR、string、gene数据库进行blastp数据库比对，剩余未预测出ORF的序列与NR、string、gene数据库进行blastx比对。

1.2.5序列注释结果功能分类

利用blast2go软件将注释后的序列，与GO数据库进行比对，将注释序列按照它们参与的生物学过程、构成细胞的组分以及实现的分子功能进行分类，GO注释分类可以帮助理解基因背后所代表的生物学意义。利用blastx(2.2.24)软件将注释后的序列与string 9.0数据库进行比对，该数据库选取66株已完成的基因组的蛋白序列，根据系统进化关系分类构建而成，通过与该数据库进行比对分类，可以进行基因的功能注释、归类以及蛋白进化分析。运用BLAST算法(blastx/blastp2.2.24)将注释序列与KEGG的基因数据库(GENES)进行比对，根据比对到的KO编号去查找具体的生物学通路，可以提供所分析基因可能参与的所有生物学通路。

1.2.6表达差异分析

某一个基因的表达水平的直接体现就是其转录本的丰度，转录本丰度越高，则基因的表达水平越高。在转录组测序数据分析中，通过对定位到基因组区域的测序序列(Clean reads)的计数来估计基因的表达水平。使用RESM计算表达量，由于质量剪切后，如果一个测序片段map到序列上确实只记一个count，另外还存在2种情况就是测序序列只有一部分map到了参考序列上，或者序列map到了参考序列的多个位置上，RESM会用最大似然法(Expectation-Maximization)估计一个count值[16-17]。表达量的计算使用FPKM法(Fragments per kb per million fragment)[18]，其计算公式：

式中：C为比对到isogene A的片段数，N为比对到所有isogene总片段数，L为isogene A的碱基数。最终对所有基因/转录本在各组样本中的表达进行差异显著性分析，运用RSEM和edgeR软件找出相对差异表达的基因/转录本，并对其进行可视化分析。并运用edgeR软件计算基因的表达差异[19-20]，其计算步骤可以分为以下几步：

1)使用TMM(Trimmed mean ofMvalues)对序列的count数进行均一化。假设基因g在样本k中的表达量用Ygk表示，对于2个样本k与k′可以定义：

①Log-fold_change：(基因g在2个样本表达量的差值)可以由下式定义

②Absolute expression levels由下式定义：

使用TMM均一化，是指除去所有基因中M值(Log-fold_change)过高以及过低的基因共30%，以及A值(Absolute expression levels)过高以及过低的基因共5%后对剩余基因的M值进行均一化，使其平均值为0。

2)离散度估计(Dispersion estimation)。假设基因i在样本j中的表达量可以用Yij满足负二项分布：

Yij～NB(uij，φ)

其中：uij为该分布的均值；φ表示该分布的离散度。φ与方差的关系如下：

Var(Yij)=uij(1+uijφ)

基因g离散度的对数似然只可以用下式估计

3)参数检验。使用假设检验计算负二项分布中该离散度出现的显著性P值，通常使用fisher精确检验法计算。

2 结果与分析

2.1 高质量测序数据统计

得到原始的FASTQ数据后，对其进行过滤除杂，去除测序质量较低的数据、去除序列中的接头序列、去除含有N较多的序列、去除小于25 nt的短序列后得到高质量的序列数据见表1。

表1 高质量数据统计结果Table 1 Statistics of high quality raw data

如表所示，测序所获得原始数据其中总reads数约为3.09亿条，测序总量超过30.7 GB。对测序所得原始数据进行质量控制过滤掉低质量序列后，获得高质量的reads 2.3亿条。

2.2 测序数据拼接结果统计

转录组测序数据denovo拼接产生总基因数为55 541条，可变剪切体为105 854条，获得拼接序列总读长为207 066 099 bp，平均序列长度为1 956.15 bp，其中最长片段为23 311 bp，最短片段为351 bp(表2)。

表2 拼接结果主要长度分布统计

从转录组数据拼接结果长度分布来看，denove拼接的转录本大部分长度在4 000 bp以下，主要长度集中在401～600 bp，随后呈现随转录本长度的增加，其所占总数的比例递减的趋势，其中长度在1～400 bp的转录本为7 737条，占总数比例7.31%，401～600 bp长度的转录本19 281条，占总数比例最大，所占比例为18.21%。601～800 bp长度的转录本11 147条，所占比例第二，比例为10.53%，801～1 000 bp长度的转录本7 796条，占比例为7.36%，1 001～1 200 bp长度的转录本6 054条，占总数比例5.72%，1201～1 400长度的转录本5 041条，占总数比例4.76%；1 401～1 600 bp长度的转录本4545条，占总数比例4.29%；1 601～1 800 bp长度的转录本4 165条，占总数百分比为3.93%；1 801～2 000 bp长度转录本3 755条，占总数百分比为3.54%；长度在2 000 bp以下的转录本占到总数比例的65.56%；4 600 bp长度以下的转录本占到总数比例的90.79%。

2.3 基因功能注释分类

将表达基因进行GO分析(图1)，可以用来研究皮肤转录组表达基因在不同时期的生理学功能以及参与的生物学过程。

通过GO分类统计数可以看出皮肤表达基因主要分为三大类中，分别为生物学功能(Biological process)、细胞组分(Cellular component)和分子功能(Molecular function)，共有51 078条转录本进行了GO注释。

在生物学功能中，细胞学功能(Cellular process)、代谢功能(Metabolic process)、生物学调控(Biological regulation)、生物功能调控(Regulation of biological process)和单一生物功能(Single-organism process)以及刺激应答(Response to stimulus)等方面，注释到的基因百分比最大。注释在Cellular process方面的转录本最多，为39 318条，注释在Metabolic process方面的转录本为32 353条，注释在Biological regulation方面的转录本25 423条，注释在Regulation of biological process方面的转录本为24 507条，注释在Single-organism process方面的转录本为19 432条，注释在Response to stimulus方面的转录本为17 863条。

在细胞组分中，细胞结构(Cell)、细胞器(Organelle)、细胞部分(Cell part)以及细胞膜(Membrane)等方面注释到的基因百分比最大。在细胞组分这一大类中注释到Cell和Cell part方面的转录本为41 112条，注释到Organelle方面的转录本为21 058条，注释到Membrane方面的转录本为17 852条。

在分子功能方面，结合(Binding)、催化作用(Catalytic activity)等方面注释到的基因百分比最大。在分子功能这一大类中注释到Binding方面的转录本共有35 643条，注释到Catalytic activity方面的转录本共有20 037条。

横坐标表示GO的二级分类描述，左边纵坐标表示百分比，右边纵坐标表示数量，3个颜色表示三大分类。X-axis means secondary category description of GO,the Y-axis on the left means percent of genes,the Y-axis on the right means number of genes,the different color annotations mean different classifications.图1 GO分类统计Fig.1 Gene ontology annotation

将拼接得到的基因利用COG数据库，进行功能注释、归类以及蛋白进化分析，共分析得到注释转录本21 189条(图2)。

图2 COG分类统计Fig.2 COG function classification

通过COG分类注释可以看出皮肤表达基因的同源蛋白主要功能集中在基因转录、翻译后修饰、一般功能和信号传导机制方面，其中注释到一般功能的COG数据为3 741条，KOG数据为3 503条。值得注意的是注释到未知功能方面的COG数据为398条，KOG数据竟达到1 379条。

然而，在同工酶转运和代谢、细胞动力、二级代谢生物合成、转运和催化、核结构等方面注释到极少的基因。

除了对表达基因进行GO以及COG分类注释分析外，还利用KEGG数据库对转录组拼接到的基因进行表达基因的通路分析。将所有预测基因与KEGG的基因数据库进行比对，共有26 201条预测基因KEGG数据库比对成功，33 990条比对失败。此次皮肤转录组拼接到的基因共富集到288条通路中去。其中与皮肤毛囊相关信号通路，例如，MAPK signaling pathway、p53 signaling pathway、WNT signaling pathway、NOTCH signaling pathway、HEDGEHOG signaling pathway、TGF-β signaling pathway、JAK-STAT signaling pathway和FOCAL signaling pathway等均有大量转录本预测基因比对到上述信号通路中。

2.4 基因表达差异分析

将绒山羊12个月的皮肤转录组分别进行顺序两两比较，将基因表达差异进行可视化分析见图3所示。顺序两两月份间差异基因统计如图4所示，1和 2月差异表达基因19个，2和3月差异表达基因陡增到1 059个，3和4月差异表达基因731个，4和5月差异表达基因13个，5和6月差异表达基因1个，6和7月差异表达基因418个，7和8月差异表达基因26个，8和9月差异表达基因2个，9和10月差异表达基因1个，10和11月差异表达基因47个，11和12月差异表达基因8个，12和1月差异表达基因18个。

红色点表示在2个样本中都表达的差异基因，黑色点表示在2个样本中都表达的非差异基因，黄色部分表示仅在一个样本中表达的基因。(a)为1和2月份相比较，(b)为2和3月份相比较，以此类推。下同。The red dot means DEGs between months,the black dot means non-differential genes between months,the yellow dot means gene that expressed only in one months.(a) is January compares February,(b) is February compares March,and so on.The same as below.图3 1—12月份可视化显著差异表达基因Fig.3 Visualization figure of significantly DEGs between neighbor months

图4 1—12月份皮肤差异表达基因统计Fig.4 Histogram of differentially expressed gene statistics between neighbor months

通过对1年当中12个月两两比较所找到的差异基因数量，进行简单统计后发现，皮肤基因在一年当中共发生了4次显著的基因差异变化，第1次差异基因数量显著变化发生在2—3月期间，差异基因数量达到1 059个；第2次差异基因的显著差异变化发生在3—4月期间，差异基因数量为731个；因此3月为皮肤基因变化显著且相对独立的一个时期；皮肤基因第3次出现显著差异变化是在6—7月，差异表达基因为418个； 7月之后直至10月，皮肤差异基因变化不明显；第4次皮肤差异基因变化显著的时间点发生在10—11月，皮肤差异基因为47个；之后11月直至次年2月，皮肤基因差异变化几乎为零。

通过转录组不同月份间差异基因的统计可以对绒山羊皮肤基因表达变化规律进行总结，可以得出，绒山羊皮肤基因在不同月份间具有显著的差异性，全年共出现4次皮肤基因的显著差异变化，处于绒毛启动生长时期的基因差异变化较为剧烈，而绒毛生长进入退行期和休止期时期的皮肤基因差异变化不明显，为一个缓慢的渐变过程。随后，将全年各月份间转录组数据进行两两比较，进行各月份间的差异基因比较，进行整理后得到皮肤全年月份间差异基因数据如表3。

其中发现3月份较全年其他月份的皮肤转录组具有明显的差异性(图5)。

如图5所示，发现3月份的皮肤样本转录组数据与全年其他月份皮肤样本转录组数据相比较所具有的差异表达基因数最多，说明3月份皮肤基因活动为全年最活跃时期。

并且，从表3全年各月份间皮肤差异基因统计表中发现3月份皮肤与7、8、9和10月这4个月的差异基因变化数是最多的，将3月份与7、8、9和10月这4个月均具有显著差异的基因进行筛选，从拼接产生的55 541个基因中筛选出3 740个基因。其中包括角蛋白(Keratin)基因家族相关基因以及大量的锌指蛋白基因家族(Zinc finger protein)。

表3 各月份间皮肤差异基因统计Table 3 Statistics of differentially expressed gene between each two months

图5 3月份较全年其他月份间差异基因可视化图Fig.5 Visualization figure of DEGs between March and other months

2.5 绒毛生长相关信号通路部分基因的差异表达规律

某一个基因表达水平的直接体现就是其转录本的丰度，转录本丰度越高，则基因表达水平越高。在RNA-seq分析中，通过对定位到基因组区域的测序序列(Clean reads)的计数来估计基因的表达水平，一般采用每百万条读段中定位到转录本每一千碱基上的片段数(Fragment per kilobase of transcript per million mapped reads,FPKM)来表达。对于绒毛生长相关的信号通路进行相关基因表达量的整合分析，其中同时富集到Wnt，TGF-β，Cell cycle信号通路中的SMAD家族成员SMAD4和SMAD3基因的基因表达FPKM值见图6，同时富集到MAPK,TGF-β，Cell cycle信号通路中的转化生长因子β(Transforming growth factor beta，TGFB)家族的TGFB1，TGFB2和TGFB3的表达FPKM值(图7)，作为毛囊生长周期的Marker基因的淋巴增强结和因子1(Lymphoid enhancer-binding factor 1，LEF1)富集到wnt信号通路中，其基因表达变化FPKM值(图8)，富集到TGF-β和HEDGEHOG信号通路的骨形态蛋白4(Bone morphogenetic protein 4，BMP4)基因表达变化FPKM值见图9。

由图6 可知SMAD3与SMAD4基因在3月出现基因的波动，并且SMAD4基因上调、SMAD3基因下调，且在随后的月份SMAD4基因表达量上升，SMAD3基因表达量下降。3月为皮肤基因最为活跃的月份，并且绒毛生长处在启动前期，推测SMAD基因在绒毛启动生长时期的基因调控中具有特定的功能作用。由图7可知，TGFB基因家族的TGFB1、TGFB2、TGFB3基因在全年皮肤中的表达量具有显著差异，TGFB1与TGFB3基因的表达趋势相同，在绒毛生长的休止期表达量较高，在绒毛生长期表达量下降。而TGFB2基因的表达趋势相反，在绒毛生长的休止期表达量较低，在绒毛生长期表达量上升。推测TGFB基因家族内不同的基因在绒毛生长调控中的作用不同，甚至可能相反。由图8可知，LEF1基因在全年皮肤中的表达量具有差异性，在绒毛的生长期，该基因的表达量显著上调。由图9可知，BMP4基因在3月皮肤中表达量显著下降，并且其表达量在绒毛生长期呈现上升趋势。

图6 SMAD基因表达FPKM值Fig.6 FPKM value of SMAD family member 3 and member 4 gene

图7 TGFB基因表达FPKM值Fig.7 FPKM value of transforming growth factor beta family gene

图8 LEF1基因表达FPKM值Fig.8 FPKM value of lymphoid enhancer-binding factor 1 gene

图9 BMP4基因表达FPKM值Fig.9 FPKM value of bone morphogenetic protein 4 gene

3 讨 论

3.1 皮肤转录组注释结果

此次阿尔巴斯型内蒙古绒山羊皮肤转录组测序数据denovo拼接产生总基因数为55 541条，isogene(Unigene)为105 854条。其中共有51 078条转录本进行了GO注释，21 189条转录本在COG数据库分析得到注释。将所有预测基因与KEGG基因数据库进行比对，共有26 201条预测基因比对成功，33 990条比对失败。阿尔巴斯型内蒙古绒山羊是我国重要的绒山羊品种，具有重要的经济价值以及物种遗传资源价值。然而，目前对于这种具有重要经济以及遗传资源价值的动物纤维生长的分子机理却知之甚少，此次针对阿尔巴斯型绒山羊皮肤绒毛生长全年时期的转录组测序获得的转录本数据对了解绒毛生长相关基因以及基因表达变化在绒毛生长时期间的差异性具有重要作用。通过皮肤转录组测序序列的功能注释以及功能分类来看，发现皮肤基因GO注释在生物学功能分类中，细胞过程、代谢过程、生物学调控、生物功能调控以及刺激应答方面注释到的转录本比例最大，说明皮肤表达基因在调控以及次级应答方面较为丰富；在分子功能分类中，结合以及催化活性作用方面注释到的转录本比例最大，说明皮肤表达基因中与结合和催化活性相关的功能基因较为丰富；而通过COG数据库的比对注释，看到皮肤表达基因的同源蛋白功能主要集中在基因转录、翻译后修饰和信号传导机制方面，反而在同工酶转运和代谢、细胞动力、二级代谢生物合成、转运和催化以及核结构方面注释到的基因极少，说明皮肤功能基因表达蛋白在信号传导等方面较为丰富；通过将皮肤转录组数据与KEGG数据库比对，得到皮肤相关转录本较为富集的几条通路分别为MAPK signaling pathway[21],Wnt signaling pathway[22],Notch signaling pathway[23],Hedgehog signaling pathway[24],TGF-β signaling pathway[25],JAK-STAT signaling pathway[26]以及Focal signaling pathway[27]均与绒毛生长基因及绒毛品质相关基因的信号通路相关。

3.2 皮肤全年各月份差异表达基因变化规律

毛发对于哺乳动物是至关重要的，在自然界动物的毛发普遍存在周期性生长的规律。毛发来源于皮肤组织中生长的毛囊结构，毛囊具有复杂的生理结构并且由多种细胞构成，是存在于皮肤中的一种微小的可再生器官。动物出生后毛囊的生长就处在不断变化当中，毛囊具有自我更新和周期性生长的规律特点，毛囊的生长具有周期性，分为生长期(Anagen)、退行期(Catagen)和休止期(Telogen)，并且3个时期的交替变化构成毛囊活动的1个周期。所以研究绒山羊毛囊的变化规律和其调控基因的表达差异，对提高绒山羊的经济价值具有重要的指导作用。

针对阿尔巴斯绒山羊次级毛囊的组织切片研究发现，1—3月的皮肤厚度，初次级毛囊的长和深，初次级毛球宽度、密度和活性等数据变化不明显，基本没有变化，各性状数据统计值较低；而从4月开始，毛囊根部细胞分裂加速并向真皮层延伸，各形态数据也开始变高，直至7月绒毛长出体表，8、9月大部分绒毛长出体表，此时毛囊各性状统计值大多也达到了全年最大值，此时期为绒毛生长的高峰时期；而10月开始毛囊毛球细胞开始逐渐衰老死亡，毛乳头细胞也凯式萎缩，此时毛囊各形态学数据统计值也开始下降变小；到达12月的时候，毛囊根部上升到皮脂腺附近，毛囊各统计数值也达到全年最低值，直至次年2月；通过皮肤毛囊的切片观察认为绒山羊毛囊的生长期为4—9月，为期6个月；退行期为10—11月，为期2个月；休止期为12月—第2年3月，为期4个月[28-29]。

通过对绒山羊全年12个月的皮肤进行转录组测序，发现绒山羊皮肤表达基因在全年具有明显的差异性。通过将月份间转录组数据进行两两比较发现，1月处于毛囊生长的退行期，由1月相较于其他月份的差异基因变化来看，1月与3月的差异基因达到1 059个，随后4月开始直至11月各月份相较于1月的差异表达基因数平稳且差异数量不大，说明，3月出现由基因调控的绒毛生长时期的转换，也就是由休止期向生长期转换，并且调控基因对毛囊时期转换的调控具有时限性，由休止期到生长期调控基因启动生长后，维持毛囊生长期的基因与退行期的基因差异数不大；2月处于绒毛退行期末期，由2月与全年其他各月份相比较来看，同样在2月与3月之间出现显著的差异表达基因数量剧增，且随后在2月与7、8和9月的比较中表达差异基因数也出现数量上的增加，说明2月与3月间出现的基因表达的剧烈变化为毛囊周期转的调控时期，随后在7、8和9月，毛囊生长较2月再次出现基因差异的波动，即毛囊由生长前期转换为生长后期也就是快速生长期；3月毛囊处于休止后期，从3月较全年其他各月份的差异基因数目上来看，3月与7、8、9和10月的差异基因数目较大，说明这4个月毛囊处于快速生长期；4月与全年其他月份相比较，与3月差异基因数为731个，出现1个基因波动，随后在8、9、10月出现差异基因波动，而相反5、6和7月差异基因平稳且不明显，说明4、5、6和7月为生长前期，基因表达稳定，随后与8、9和10月出现差异基因波动，说明生长期由生长前期进入生长后期，也就是快速生长期；随后的5、6和7月较全年其他各月份的差异基因变化与4月相仿，仅有7月与8、9和10月的差异基因不再明显，说明6、7月的差异基因波动，已使得绒毛生长期由前期向快速生长后期进行了转换；8、9和10月与全年各月份的差异基因变化趋势基本一致，即与3月出现剧烈的基因波动，而与4、5和6月的差异基因变化波动趋于剧烈，也就是说8月与6月的基因差异数较8月与4月的基因差异数大，说明在绒毛生长期，越接近快速生长期，基因活动越剧烈；其中10月较为特殊，10月与生长期的7、8和9月差异基因不明显，与11、12月退行期的差异基因也不明显，说明10月为绒毛周期由生长期向退行期的过度时期；而11、12 月处于毛囊的退行期，这2个月与全年其他月份的差异基因变化趋势一样，即3月出现基因波动后，较7、8、9月差异基因变化剧烈，说明休止期与快速生长期的差异基因具有较大波动；通过对全年12个月皮肤转录组差异表达基因的数据将绒山羊毛囊生长期的3个时期分为，生长期为4—10月，其中4—7为生长前期，8、9、10月为生长后期也就是快速生长期；10月为生长期向退行期过渡时期，11、12月为退行期，1—3月为休止期，其中3月为休止末期，3月基因表达最为活跃。

说明毛囊的周期间转换与基因间的差异变动具有相关性，基因间构成的复杂网络调控毛囊的周期间转换。并且在绒毛启动生长时的基因差异变化剧烈，也就是休止期到生长期基因水平的表达存在剧烈波动，而绒毛停止生长的基因差异变化则相对缓慢，绒毛停止生长的基因水平的变化不会出现短时间内的剧烈波动，而是1个循序渐进的过程。所以针对绒毛启动生长的相关基因研究对于发现调控绒毛周期变化以及绒毛生长等影响绒山羊产绒性状的基因具有潜在研究价值，并且在基因水平休止期到生长期的基因变化较为剧烈，差异较为明显，针对基因表达研究具有优势。

4 结 论

通过转录组不同月份间差异基因的统计得出皮肤差异基因变化规律与绒毛生长周期的相关性，每一次皮肤差异基因剧烈变化都伴随着绒毛生长周期的转换，2—3月出现第1次差异基因剧烈变化，绒山羊次级毛囊生长期进入休止末期，为下1个生长期的到来准备蓄能。随后3月与4月之间出现第2次剧烈变化，次级毛囊由休止期进入生长期，直至11月。6月与7月之间出现第3次差异基因剧烈变化，将生长期分为生长前期与生长后期，10月与11月之间出现第4次差异基因剧烈变化，绒毛生长周期由生长期转换为退行期，直至次年3月进入休止期。并且通过全年两两月份顺序比较的差异基因的变化可以得出，绒毛毛囊启动生长时期的基因变化剧烈，而生长期到退行期基因变化缓慢，没有剧烈的短时间内的基因表达差异变化。