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

趙 濛 劉志紅* 奈日樂 謝遇春 馬麗娜 米 璐 李金泉，2 李 婕

(1.內(nèi)蒙古農(nóng)業(yè)大學(xué) 動物科學(xué)學(xué)院，呼和浩特 010018；2.內(nèi)蒙古自治區(qū)動物遺傳育種與繁殖重點實驗室，呼和浩特 010018；3.東北農(nóng)業(yè)大學(xué) 動物科學(xué)技術(shù)學(xué)院，哈爾濱 150030)

山羊絨是重要的毛紡織原料，具有纖細、柔軟、輕薄等優(yōu)點。內(nèi)蒙古阿爾巴斯型絨山羊是我國重要的絨肉兼用型絨山羊品種，主要分布于內(nèi)蒙古自治區(qū)伊克昭盟鄂托克旗境內(nèi)，是經(jīng)過長期自然選擇和人工選育而成的優(yōu)良品種，該品種可適應(yīng)西北地區(qū)寒冷干旱和多風的氣候環(huán)境。絨毛是來源于絨山羊下層毛被的一種優(yōu)良的紡織原料，絨山羊毛被具有2種截然不同的纖維結(jié)構(gòu)，位于毛被上層的是粗且硬的硬質(zhì)粗毛，位于毛被下層的是細且軟的絨毛。絨毛纖維來源于生于皮膚中的次級毛囊結(jié)構(gòu)，粗毛來源于生于皮膚中的初級毛囊結(jié)構(gòu)[1-3]。

出生后的哺乳動物皮膚毛囊生長處在動態(tài)變化當中，其毛囊生長變化具有往復(fù)循環(huán)規(guī)律性，也就是出生后毛囊生長具有周期性，其1個生長周期大致可分為生長期、退行期和休止期[4-6]。近年來，科學(xué)研究針對控制毛囊周期生長規(guī)律的分子調(diào)控機制進行了大量的探索，發(fā)現(xiàn)某些信號通路在毛囊周期生長中具有重要的作用[7-11]。毛囊生長周期內(nèi)不同階段的轉(zhuǎn)換是由于不同的信號通路間相互促進以及抑制作用的結(jié)果。規(guī)律性的調(diào)控機制作用導(dǎo)致了毛囊生長的規(guī)律性以及生長周期間的轉(zhuǎn)換。絨山羊皮膚次級毛囊同樣具有類似的生長周期性，并且,1個生長周期具有生長期、退行期以及休止期之分，但是針對絨山羊次級毛囊生長周期性的分子調(diào)控機制研究甚少。所以，從分子生物學(xué)角度研究絨毛周期生長中不同階段的基因表達情況以及不同階段的基因表達差異性對了解調(diào)控其生長分子機制具有十分重要的作用。

Jiang等[12]通過構(gòu)建抑制消減雜交文庫的方法對絨山羊絨毛生長的生長期和休止期皮膚樣本的差異表達基因進行研究發(fā)現(xiàn)，生長期和休止期絨山羊皮膚具有差異表達基因45個，其中23個具有明確的功能。但這種方法對研究差異基因具有隨機挑選的特性，并不能對基因表達的全譜進行深入研究。隨著研究技術(shù)的進步，第二代高通量測序技術(shù)應(yīng)用于差異基因的篩選研究方面，使得對基因表達全譜的研究成為可能。RNA-Seq技術(shù)是一種高通量的測序技術(shù)，可以針對基因轉(zhuǎn)錄本進行深入研究，發(fā)掘低豐度轉(zhuǎn)錄本和新轉(zhuǎn)錄本以及鑒定不同樣本間的差異表達轉(zhuǎn)錄本[13-15]。本試驗采用第二代高通量測序技術(shù)針對絨毛生長不同階段的皮膚樣本進行轉(zhuǎn)錄組測序，旨在研究絨毛生長不同階段的皮膚差異表達基因以及基因差異表達對調(diào)控絨毛周期轉(zhuǎn)換的相關(guān)性，這些絨毛生長不同階段的差異基因所具有的生物學(xué)功能對深入了解絨毛生長的調(diào)控機理具有重要的作用，并對研究絨毛生長調(diào)控機理奠定理論基礎(chǔ)。

1 材料與方法

1.1 試驗動物及皮膚組織樣品采集

阿爾巴斯型內(nèi)蒙古絨山羊是絨肉兼用的優(yōu)良品種，本試驗選取內(nèi)蒙古伊克昭盟阿爾巴斯絨山羊種羊場的阿爾巴斯型內(nèi)蒙古絨山羊3只，試驗用羊選取條件為：選自同一放牧環(huán)境下、體重體尺相當、周齡相同、無親緣關(guān)系、生長狀況良好的周歲母羊，采樣周期為期1年，起于2011年12月，止于2012年12月，采樣周期內(nèi)試驗用羊的飼養(yǎng)管理制度與當?shù)仫曫B(yǎng)管理一致，且試驗用羊均未產(chǎn)羔。皮膚樣本采集時間為每月初。皮膚樣本采集部位為體側(cè)體中線部距離肩胛骨10～15 cm處，采取手術(shù)直接剝離大小為3 cm直徑的皮膚樣本，皮膚樣本手術(shù)剝離時間控制在5 min之內(nèi)，采集后樣本進行PBS沖洗后快速液氮冷凍，然后放入液氮罐帶回實驗室存儲于-80 ℃ 冰箱中備用。

1.2 試驗方法

1.2.1皮膚總RNA提取

采用RNAiso Plus 試劑盒(Trizol法)操作說明分別提取3只絨山羊體側(cè)皮膚組織的總RNA，無菌無酶水溶解，分別使用微量紫外分光光度計以及2100 bioanalyzer對所提取皮膚總RNA進行濃度、純度以及完整度的檢測。總RNA置于-80 ℃冰箱保存。

1.2.2測序文庫構(gòu)建

根據(jù)RNA樣本濃度，分別將3只絨山羊皮膚的總RNA等量混池，為構(gòu)建12個月的cDNA文庫做準備。轉(zhuǎn)錄組測序cDNA文庫構(gòu)建參考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循環(huán)次數(shù)為15次，獲得測序文庫，隨后使用TBS380對文庫cDNA進行質(zhì)檢。

1.2.3轉(zhuǎn)錄組測序

使用Illumina HiSeqTM2000測序平臺對cDNA進行雙末端測序，進行2×100 bp 測序試驗，由上海美吉生物醫(yī)藥科技有限公司測序。

1.2.4測序數(shù)據(jù)除雜、拼接及拼接結(jié)果注釋

測序得到的原始fastq數(shù)據(jù)，運用SeqPrep和condetri_v2.0.pl軟件進行過濾除雜，去除測序質(zhì)量較低的reads(測序質(zhì)量值小于25)，去除reads中的adapter序列，去除含有較多N的reads，去除小于25 nt的短序列，得到高質(zhì)量的序列數(shù)據(jù)。對除雜過濾后的序列運用Trinity軟件進行無參考基因組信息de novo拼接，并將組裝結(jié)果利用Trinity軟件進行ORF預(yù)測，并運用BLAST(Version 2.2.25)軟件將預(yù)測出ORF的蛋白序列與未預(yù)測出ORF的序列分別注釋，將預(yù)測出的ORF的蛋白序列分別于NR、string、gene數(shù)據(jù)庫進行blastp數(shù)據(jù)庫比對，剩余未預(yù)測出ORF的序列與NR、string、gene數(shù)據(jù)庫進行blastx比對。

1.2.5序列注釋結(jié)果功能分類

利用blast2go軟件將注釋后的序列，與GO數(shù)據(jù)庫進行比對，將注釋序列按照它們參與的生物學(xué)過程、構(gòu)成細胞的組分以及實現(xiàn)的分子功能進行分類，GO注釋分類可以幫助理解基因背后所代表的生物學(xué)意義。利用blastx(2.2.24)軟件將注釋后的序列與string 9.0數(shù)據(jù)庫進行比對，該數(shù)據(jù)庫選取66株已完成的基因組的蛋白序列，根據(jù)系統(tǒng)進化關(guān)系分類構(gòu)建而成，通過與該數(shù)據(jù)庫進行比對分類，可以進行基因的功能注釋、歸類以及蛋白進化分析。運用BLAST算法(blastx/blastp2.2.24)將注釋序列與KEGG的基因數(shù)據(jù)庫(GENES)進行比對，根據(jù)比對到的KO編號去查找具體的生物學(xué)通路，可以提供所分析基因可能參與的所有生物學(xué)通路。

1.2.6表達差異分析

某一個基因的表達水平的直接體現(xiàn)就是其轉(zhuǎn)錄本的豐度，轉(zhuǎn)錄本豐度越高，則基因的表達水平越高。在轉(zhuǎn)錄組測序數(shù)據(jù)分析中，通過對定位到基因組區(qū)域的測序序列(Clean reads)的計數(shù)來估計基因的表達水平。使用RESM計算表達量，由于質(zhì)量剪切后，如果一個測序片段map到序列上確實只記一個count，另外還存在2種情況就是測序序列只有一部分map到了參考序列上，或者序列map到了參考序列的多個位置上，RESM會用最大似然法(Expectation-Maximization)估計一個count值[16-17]。表達量的計算使用FPKM法(Fragments per kb per million fragment)[18]，其計算公式：

式中：C為比對到isogene A的片段數(shù)，N為比對到所有isogene總片段數(shù)，L為isogene A的堿基數(shù)。最終對所有基因/轉(zhuǎn)錄本在各組樣本中的表達進行差異顯著性分析，運用RSEM和edgeR軟件找出相對差異表達的基因/轉(zhuǎn)錄本，并對其進行可視化分析。并運用edgeR軟件計算基因的表達差異[19-20]，其計算步驟可以分為以下幾步：

1)使用TMM(Trimmed mean ofMvalues)對序列的count數(shù)進行均一化。假設(shè)基因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)。假設(shè)基因i在樣本j中的表達量可以用Yij滿足負二項分布：

Yij～NB(uij，φ)

其中：uij為該分布的均值；φ表示該分布的離散度。φ與方差的關(guān)系如下：

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

基因g離散度的對數(shù)似然只可以用下式估計

3)參數(shù)檢驗。使用假設(shè)檢驗計算負二項分布中該離散度出現(xiàn)的顯著性P值，通常使用fisher精確檢驗法計算。

2 結(jié)果與分析

2.1 高質(zhì)量測序數(shù)據(jù)統(tǒng)計

得到原始的FASTQ數(shù)據(jù)后，對其進行過濾除雜，去除測序質(zhì)量較低的數(shù)據(jù)、去除序列中的接頭序列、去除含有N較多的序列、去除小于25 nt的短序列后得到高質(zhì)量的序列數(shù)據(jù)見表1。

表1 高質(zhì)量數(shù)據(jù)統(tǒng)計結(jié)果Table 1 Statistics of high quality raw data

如表所示，測序所獲得原始數(shù)據(jù)其中總reads數(shù)約為3.09億條，測序總量超過30.7 GB。對測序所得原始數(shù)據(jù)進行質(zhì)量控制過濾掉低質(zhì)量序列后，獲得高質(zhì)量的reads 2.3億條。

2.2 測序數(shù)據(jù)拼接結(jié)果統(tǒng)計

轉(zhuǎn)錄組測序數(shù)據(jù)denovo拼接產(chǎn)生總基因數(shù)為55 541條，可變剪切體為105 854條，獲得拼接序列總讀長為207 066 099 bp，平均序列長度為1 956.15 bp，其中最長片段為23 311 bp，最短片段為351 bp(表2)。

表2 拼接結(jié)果主要長度分布統(tǒng)計

從轉(zhuǎn)錄組數(shù)據(jù)拼接結(jié)果長度分布來看，denove拼接的轉(zhuǎn)錄本大部分長度在4 000 bp以下，主要長度集中在401～600 bp，隨后呈現(xiàn)隨轉(zhuǎn)錄本長度的增加，其所占總數(shù)的比例遞減的趨勢，其中長度在1～400 bp的轉(zhuǎn)錄本為7 737條，占總數(shù)比例7.31%，401～600 bp長度的轉(zhuǎn)錄本19 281條，占總數(shù)比例最大，所占比例為18.21%。601～800 bp長度的轉(zhuǎn)錄本11 147條，所占比例第二，比例為10.53%，801～1 000 bp長度的轉(zhuǎn)錄本7 796條，占比例為7.36%，1 001～1 200 bp長度的轉(zhuǎn)錄本6 054條，占總數(shù)比例5.72%，1201～1 400長度的轉(zhuǎn)錄本5 041條，占總數(shù)比例4.76%；1 401～1 600 bp長度的轉(zhuǎn)錄本4545條，占總數(shù)比例4.29%；1 601～1 800 bp長度的轉(zhuǎn)錄本4 165條，占總數(shù)百分比為3.93%；1 801～2 000 bp長度轉(zhuǎn)錄本3 755條，占總數(shù)百分比為3.54%；長度在2 000 bp以下的轉(zhuǎn)錄本占到總數(shù)比例的65.56%；4 600 bp長度以下的轉(zhuǎn)錄本占到總數(shù)比例的90.79%。

2.3 基因功能注釋分類

將表達基因進行GO分析(圖1)，可以用來研究皮膚轉(zhuǎn)錄組表達基因在不同時期的生理學(xué)功能以及參與的生物學(xué)過程。

通過GO分類統(tǒng)計數(shù)可以看出皮膚表達基因主要分為三大類中，分別為生物學(xué)功能(Biological process)、細胞組分(Cellular component)和分子功能(Molecular function)，共有51 078條轉(zhuǎn)錄本進行了GO注釋。

在生物學(xué)功能中，細胞學(xué)功能(Cellular process)、代謝功能(Metabolic process)、生物學(xué)調(diào)控(Biological regulation)、生物功能調(diào)控(Regulation of biological process)和單一生物功能(Single-organism process)以及刺激應(yīng)答(Response to stimulus)等方面，注釋到的基因百分比最大。注釋在Cellular process方面的轉(zhuǎn)錄本最多，為39 318條，注釋在Metabolic process方面的轉(zhuǎn)錄本為32 353條，注釋在Biological regulation方面的轉(zhuǎn)錄本25 423條，注釋在Regulation of biological process方面的轉(zhuǎn)錄本為24 507條，注釋在Single-organism process方面的轉(zhuǎn)錄本為19 432條，注釋在Response to stimulus方面的轉(zhuǎn)錄本為17 863條。

在細胞組分中，細胞結(jié)構(gòu)(Cell)、細胞器(Organelle)、細胞部分(Cell part)以及細胞膜(Membrane)等方面注釋到的基因百分比最大。在細胞組分這一大類中注釋到Cell和Cell part方面的轉(zhuǎn)錄本為41 112條，注釋到Organelle方面的轉(zhuǎn)錄本為21 058條，注釋到Membrane方面的轉(zhuǎn)錄本為17 852條。

在分子功能方面，結(jié)合(Binding)、催化作用(Catalytic activity)等方面注釋到的基因百分比最大。在分子功能這一大類中注釋到Binding方面的轉(zhuǎn)錄本共有35 643條，注釋到Catalytic activity方面的轉(zhuǎn)錄本共有20 037條。

橫坐標表示GO的二級分類描述，左邊縱坐標表示百分比，右邊縱坐標表示數(shù)量，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分類統(tǒng)計Fig.1 Gene ontology annotation

將拼接得到的基因利用COG數(shù)據(jù)庫，進行功能注釋、歸類以及蛋白進化分析，共分析得到注釋轉(zhuǎn)錄本21 189條(圖2)。

圖2 COG分類統(tǒng)計Fig.2 COG function classification

通過COG分類注釋可以看出皮膚表達基因的同源蛋白主要功能集中在基因轉(zhuǎn)錄、翻譯后修飾、一般功能和信號傳導(dǎo)機制方面，其中注釋到一般功能的COG數(shù)據(jù)為3 741條，KOG數(shù)據(jù)為3 503條。值得注意的是注釋到未知功能方面的COG數(shù)據(jù)為398條，KOG數(shù)據(jù)竟達到1 379條。

然而，在同工酶轉(zhuǎn)運和代謝、細胞動力、二級代謝生物合成、轉(zhuǎn)運和催化、核結(jié)構(gòu)等方面注釋到極少的基因。

除了對表達基因進行GO以及COG分類注釋分析外，還利用KEGG數(shù)據(jù)庫對轉(zhuǎn)錄組拼接到的基因進行表達基因的通路分析。將所有預(yù)測基因與KEGG的基因數(shù)據(jù)庫進行比對，共有26 201條預(yù)測基因KEGG數(shù)據(jù)庫比對成功，33 990條比對失敗。此次皮膚轉(zhuǎn)錄組拼接到的基因共富集到288條通路中去。其中與皮膚毛囊相關(guān)信號通路，例如，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等均有大量轉(zhuǎn)錄本預(yù)測基因比對到上述信號通路中。

2.4 基因表達差異分析

將絨山羊12個月的皮膚轉(zhuǎn)錄組分別進行順序兩兩比較，將基因表達差異進行可視化分析見圖3所示。順序兩兩月份間差異基因統(tǒng)計如圖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月份皮膚差異表達基因統(tǒng)計Fig.4 Histogram of differentially expressed gene statistics between neighbor months

通過對1年當中12個月兩兩比較所找到的差異基因數(shù)量，進行簡單統(tǒng)計后發(fā)現(xiàn)，皮膚基因在一年當中共發(fā)生了4次顯著的基因差異變化，第1次差異基因數(shù)量顯著變化發(fā)生在2—3月期間，差異基因數(shù)量達到1 059個；第2次差異基因的顯著差異變化發(fā)生在3—4月期間，差異基因數(shù)量為731個；因此3月為皮膚基因變化顯著且相對獨立的一個時期；皮膚基因第3次出現(xiàn)顯著差異變化是在6—7月，差異表達基因為418個； 7月之后直至10月，皮膚差異基因變化不明顯；第4次皮膚差異基因變化顯著的時間點發(fā)生在10—11月，皮膚差異基因為47個；之后11月直至次年2月，皮膚基因差異變化幾乎為零。

通過轉(zhuǎn)錄組不同月份間差異基因的統(tǒng)計可以對絨山羊皮膚基因表達變化規(guī)律進行總結(jié)，可以得出，絨山羊皮膚基因在不同月份間具有顯著的差異性，全年共出現(xiàn)4次皮膚基因的顯著差異變化，處于絨毛啟動生長時期的基因差異變化較為劇烈，而絨毛生長進入退行期和休止期時期的皮膚基因差異變化不明顯，為一個緩慢的漸變過程。隨后，將全年各月份間轉(zhuǎn)錄組數(shù)據(jù)進行兩兩比較，進行各月份間的差異基因比較，進行整理后得到皮膚全年月份間差異基因數(shù)據(jù)如表3。

其中發(fā)現(xiàn)3月份較全年其他月份的皮膚轉(zhuǎn)錄組具有明顯的差異性(圖5)。

如圖5所示，發(fā)現(xiàn)3月份的皮膚樣本轉(zhuǎn)錄組數(shù)據(jù)與全年其他月份皮膚樣本轉(zhuǎn)錄組數(shù)據(jù)相比較所具有的差異表達基因數(shù)最多，說明3月份皮膚基因活動為全年最活躍時期。

并且，從表3全年各月份間皮膚差異基因統(tǒng)計表中發(fā)現(xiàn)3月份皮膚與7、8、9和10月這4個月的差異基因變化數(shù)是最多的，將3月份與7、8、9和10月這4個月均具有顯著差異的基因進行篩選，從拼接產(chǎn)生的55 541個基因中篩選出3 740個基因。其中包括角蛋白(Keratin)基因家族相關(guān)基因以及大量的鋅指蛋白基因家族(Zinc finger protein)。

表3 各月份間皮膚差異基因統(tǒng)計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 絨毛生長相關(guān)信號通路部分基因的差異表達規(guī)律

某一個基因表達水平的直接體現(xiàn)就是其轉(zhuǎn)錄本的豐度，轉(zhuǎn)錄本豐度越高，則基因表達水平越高。在RNA-seq分析中，通過對定位到基因組區(qū)域的測序序列(Clean reads)的計數(shù)來估計基因的表達水平，一般采用每百萬條讀段中定位到轉(zhuǎn)錄本每一千堿基上的片段數(shù)(Fragment per kilobase of transcript per million mapped reads,FPKM)來表達。對于絨毛生長相關(guān)的信號通路進行相關(guān)基因表達量的整合分析，其中同時富集到Wnt，TGF-β，Cell cycle信號通路中的SMAD家族成員SMAD4和SMAD3基因的基因表達FPKM值見圖6，同時富集到MAPK,TGF-β，Cell cycle信號通路中的轉(zhuǎn)化生長因子β(Transforming growth factor beta，TGFB)家族的TGFB1，TGFB2和TGFB3的表達FPKM值(圖7)，作為毛囊生長周期的Marker基因的淋巴增強結(jié)和因子1(Lymphoid enhancer-binding factor 1，LEF1)富集到wnt信號通路中，其基因表達變化FPKM值(圖8)，富集到TGF-β和HEDGEHOG信號通路的骨形態(tài)蛋白4(Bone morphogenetic protein 4，BMP4)基因表達變化FPKM值見圖9。

由圖6 可知SMAD3與SMAD4基因在3月出現(xiàn)基因的波動，并且SMAD4基因上調(diào)、SMAD3基因下調(diào)，且在隨后的月份SMAD4基因表達量上升，SMAD3基因表達量下降。3月為皮膚基因最為活躍的月份，并且絨毛生長處在啟動前期，推測SMAD基因在絨毛啟動生長時期的基因調(diào)控中具有特定的功能作用。由圖7可知，TGFB基因家族的TGFB1、TGFB2、TGFB3基因在全年皮膚中的表達量具有顯著差異，TGFB1與TGFB3基因的表達趨勢相同，在絨毛生長的休止期表達量較高，在絨毛生長期表達量下降。而TGFB2基因的表達趨勢相反，在絨毛生長的休止期表達量較低，在絨毛生長期表達量上升。推測TGFB基因家族內(nèi)不同的基因在絨毛生長調(diào)控中的作用不同，甚至可能相反。由圖8可知，LEF1基因在全年皮膚中的表達量具有差異性，在絨毛的生長期，該基因的表達量顯著上調(diào)。由圖9可知，BMP4基因在3月皮膚中表達量顯著下降，并且其表達量在絨毛生長期呈現(xiàn)上升趨勢。

圖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 皮膚轉(zhuǎn)錄組注釋結(jié)果

此次阿爾巴斯型內(nèi)蒙古絨山羊皮膚轉(zhuǎn)錄組測序數(shù)據(jù)denovo拼接產(chǎn)生總基因數(shù)為55 541條，isogene(Unigene)為105 854條。其中共有51 078條轉(zhuǎn)錄本進行了GO注釋，21 189條轉(zhuǎn)錄本在COG數(shù)據(jù)庫分析得到注釋。將所有預(yù)測基因與KEGG基因數(shù)據(jù)庫進行比對，共有26 201條預(yù)測基因比對成功，33 990條比對失敗。阿爾巴斯型內(nèi)蒙古絨山羊是我國重要的絨山羊品種，具有重要的經(jīng)濟價值以及物種遺傳資源價值。然而，目前對于這種具有重要經(jīng)濟以及遺傳資源價值的動物纖維生長的分子機理卻知之甚少，此次針對阿爾巴斯型絨山羊皮膚絨毛生長全年時期的轉(zhuǎn)錄組測序獲得的轉(zhuǎn)錄本數(shù)據(jù)對了解絨毛生長相關(guān)基因以及基因表達變化在絨毛生長時期間的差異性具有重要作用。通過皮膚轉(zhuǎn)錄組測序序列的功能注釋以及功能分類來看，發(fā)現(xiàn)皮膚基因GO注釋在生物學(xué)功能分類中，細胞過程、代謝過程、生物學(xué)調(diào)控、生物功能調(diào)控以及刺激應(yīng)答方面注釋到的轉(zhuǎn)錄本比例最大，說明皮膚表達基因在調(diào)控以及次級應(yīng)答方面較為豐富；在分子功能分類中，結(jié)合以及催化活性作用方面注釋到的轉(zhuǎn)錄本比例最大，說明皮膚表達基因中與結(jié)合和催化活性相關(guān)的功能基因較為豐富；而通過COG數(shù)據(jù)庫的比對注釋，看到皮膚表達基因的同源蛋白功能主要集中在基因轉(zhuǎn)錄、翻譯后修飾和信號傳導(dǎo)機制方面，反而在同工酶轉(zhuǎn)運和代謝、細胞動力、二級代謝生物合成、轉(zhuǎn)運和催化以及核結(jié)構(gòu)方面注釋到的基因極少，說明皮膚功能基因表達蛋白在信號傳導(dǎo)等方面較為豐富；通過將皮膚轉(zhuǎn)錄組數(shù)據(jù)與KEGG數(shù)據(jù)庫比對，得到皮膚相關(guān)轉(zhuǎn)錄本較為富集的幾條通路分別為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]均與絨毛生長基因及絨毛品質(zhì)相關(guān)基因的信號通路相關(guān)。

3.2 皮膚全年各月份差異表達基因變化規(guī)律

毛發(fā)對于哺乳動物是至關(guān)重要的，在自然界動物的毛發(fā)普遍存在周期性生長的規(guī)律。毛發(fā)來源于皮膚組織中生長的毛囊結(jié)構(gòu)，毛囊具有復(fù)雜的生理結(jié)構(gòu)并且由多種細胞構(gòu)成，是存在于皮膚中的一種微小的可再生器官。動物出生后毛囊的生長就處在不斷變化當中，毛囊具有自我更新和周期性生長的規(guī)律特點，毛囊的生長具有周期性，分為生長期(Anagen)、退行期(Catagen)和休止期(Telogen)，并且3個時期的交替變化構(gòu)成毛囊活動的1個周期。所以研究絨山羊毛囊的變化規(guī)律和其調(diào)控基因的表達差異，對提高絨山羊的經(jīng)濟價值具有重要的指導(dǎo)作用。

針對阿爾巴斯絨山羊次級毛囊的組織切片研究發(fā)現(xiàn)，1—3月的皮膚厚度，初次級毛囊的長和深，初次級毛球?qū)挾?、密度和活性等?shù)據(jù)變化不明顯，基本沒有變化，各性狀數(shù)據(jù)統(tǒng)計值較低；而從4月開始，毛囊根部細胞分裂加速并向真皮層延伸，各形態(tài)數(shù)據(jù)也開始變高，直至7月絨毛長出體表，8、9月大部分絨毛長出體表，此時毛囊各性狀統(tǒng)計值大多也達到了全年最大值，此時期為絨毛生長的高峰時期；而10月開始毛囊毛球細胞開始逐漸衰老死亡，毛乳頭細胞也凱式萎縮，此時毛囊各形態(tài)學(xué)數(shù)據(jù)統(tǒng)計值也開始下降變?。坏竭_12月的時候，毛囊根部上升到皮脂腺附近，毛囊各統(tǒng)計數(shù)值也達到全年最低值，直至次年2月；通過皮膚毛囊的切片觀察認為絨山羊毛囊的生長期為4—9月，為期6個月；退行期為10—11月，為期2個月；休止期為12月—第2年3月，為期4個月[28-29]。

通過對絨山羊全年12個月的皮膚進行轉(zhuǎn)錄組測序，發(fā)現(xiàn)絨山羊皮膚表達基因在全年具有明顯的差異性。通過將月份間轉(zhuǎn)錄組數(shù)據(jù)進行兩兩比較發(fā)現(xiàn)，1月處于毛囊生長的退行期，由1月相較于其他月份的差異基因變化來看，1月與3月的差異基因達到1 059個，隨后4月開始直至11月各月份相較于1月的差異表達基因數(shù)平穩(wěn)且差異數(shù)量不大，說明，3月出現(xiàn)由基因調(diào)控的絨毛生長時期的轉(zhuǎn)換，也就是由休止期向生長期轉(zhuǎn)換，并且調(diào)控基因?qū)γ視r期轉(zhuǎn)換的調(diào)控具有時限性，由休止期到生長期調(diào)控基因啟動生長后，維持毛囊生長期的基因與退行期的基因差異數(shù)不大；2月處于絨毛退行期末期，由2月與全年其他各月份相比較來看，同樣在2月與3月之間出現(xiàn)顯著的差異表達基因數(shù)量劇增，且隨后在2月與7、8和9月的比較中表達差異基因數(shù)也出現(xiàn)數(shù)量上的增加，說明2月與3月間出現(xiàn)的基因表達的劇烈變化為毛囊周期轉(zhuǎn)的調(diào)控時期，隨后在7、8和9月，毛囊生長較2月再次出現(xiàn)基因差異的波動，即毛囊由生長前期轉(zhuǎn)換為生長后期也就是快速生長期；3月毛囊處于休止后期，從3月較全年其他各月份的差異基因數(shù)目上來看，3月與7、8、9和10月的差異基因數(shù)目較大，說明這4個月毛囊處于快速生長期；4月與全年其他月份相比較，與3月差異基因數(shù)為731個，出現(xiàn)1個基因波動，隨后在8、9、10月出現(xiàn)差異基因波動，而相反5、6和7月差異基因平穩(wěn)且不明顯，說明4、5、6和7月為生長前期，基因表達穩(wěn)定，隨后與8、9和10月出現(xiàn)差異基因波動，說明生長期由生長前期進入生長后期，也就是快速生長期；隨后的5、6和7月較全年其他各月份的差異基因變化與4月相仿，僅有7月與8、9和10月的差異基因不再明顯，說明6、7月的差異基因波動，已使得絨毛生長期由前期向快速生長后期進行了轉(zhuǎn)換；8、9和10月與全年各月份的差異基因變化趨勢基本一致，即與3月出現(xiàn)劇烈的基因波動，而與4、5和6月的差異基因變化波動趨于劇烈，也就是說8月與6月的基因差異數(shù)較8月與4月的基因差異數(shù)大，說明在絨毛生長期，越接近快速生長期，基因活動越劇烈；其中10月較為特殊，10月與生長期的7、8和9月差異基因不明顯，與11、12月退行期的差異基因也不明顯，說明10月為絨毛周期由生長期向退行期的過度時期；而11、12 月處于毛囊的退行期，這2個月與全年其他月份的差異基因變化趨勢一樣，即3月出現(xiàn)基因波動后，較7、8、9月差異基因變化劇烈，說明休止期與快速生長期的差異基因具有較大波動；通過對全年12個月皮膚轉(zhuǎn)錄組差異表達基因的數(shù)據(jù)將絨山羊毛囊生長期的3個時期分為，生長期為4—10月，其中4—7為生長前期，8、9、10月為生長后期也就是快速生長期；10月為生長期向退行期過渡時期，11、12月為退行期，1—3月為休止期，其中3月為休止末期，3月基因表達最為活躍。

說明毛囊的周期間轉(zhuǎn)換與基因間的差異變動具有相關(guān)性，基因間構(gòu)成的復(fù)雜網(wǎng)絡(luò)調(diào)控毛囊的周期間轉(zhuǎn)換。并且在絨毛啟動生長時的基因差異變化劇烈，也就是休止期到生長期基因水平的表達存在劇烈波動，而絨毛停止生長的基因差異變化則相對緩慢，絨毛停止生長的基因水平的變化不會出現(xiàn)短時間內(nèi)的劇烈波動，而是1個循序漸進的過程。所以針對絨毛啟動生長的相關(guān)基因研究對于發(fā)現(xiàn)調(diào)控絨毛周期變化以及絨毛生長等影響絨山羊產(chǎn)絨性狀的基因具有潛在研究價值，并且在基因水平休止期到生長期的基因變化較為劇烈，差異較為明顯，針對基因表達研究具有優(yōu)勢。

4 結(jié) 論

通過轉(zhuǎn)錄組不同月份間差異基因的統(tǒng)計得出皮膚差異基因變化規(guī)律與絨毛生長周期的相關(guān)性，每一次皮膚差異基因劇烈變化都伴隨著絨毛生長周期的轉(zhuǎn)換，2—3月出現(xiàn)第1次差異基因劇烈變化，絨山羊次級毛囊生長期進入休止末期，為下1個生長期的到來準備蓄能。隨后3月與4月之間出現(xiàn)第2次劇烈變化，次級毛囊由休止期進入生長期，直至11月。6月與7月之間出現(xiàn)第3次差異基因劇烈變化，將生長期分為生長前期與生長后期，10月與11月之間出現(xiàn)第4次差異基因劇烈變化，絨毛生長周期由生長期轉(zhuǎn)換為退行期，直至次年3月進入休止期。并且通過全年兩兩月份順序比較的差異基因的變化可以得出，絨毛毛囊啟動生長時期的基因變化劇烈，而生長期到退行期基因變化緩慢，沒有劇烈的短時間內(nèi)的基因表達差異變化。