基于高通量轉錄組測序阿爾巴斯型絨山羊絨毛生長周期差異表達基因分析

趙 濛 劉志紅* 奈日樂 謝遇春 馬麗娜 米 璐 李金泉，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月進入休止期。并且通過全年兩兩月份順序比較的差異基因的變化可以得出，絨毛毛囊啟動生長時期的基因變化劇烈，而生長期到退行期基因變化緩慢，沒有劇烈的短時間內的基因表達差異變化。