基于轉錄組數據挖掘藏羊立毛肌發生的關鍵基因

楚金雨，李紹梅，楊 戈，牟春燕

(華中農業大學動物科學技術學院 農業動物遺傳育種與繁殖教育部重點實驗室，武漢 430070)

皮膚是人體最大的器官，由表皮、真皮和皮下組織構成，并含有毛囊、汗腺、皮脂腺、指甲等附屬物，具有保護身體、維持水分、感覺冷暖等功能[1]。由于羊毛囊的形態發生類似于人的頭部毛囊，而綿羊的軀干皮膚所含有的毛囊和汗腺與人類腋窩皮膚的解剖結構相似[2]，所以，綿羊可成為研究人類毛發生長的優良動物模型。

毛囊是由上皮細胞和間充質之間通過信號分子介導的一系列相互作用形成的[3]。成熟的毛囊包括真皮乳頭、外根鞘、內根鞘、毛球、毛干等結構[4]，且在發育形成的過程中通常伴隨著汗腺、皮脂腺和立毛肌的出現。立毛肌是一類平滑肌，由皮膚結締組織的成纖維細胞分化形成[5]，位于表皮下方，始于毛囊底部隆起處，終止于真皮細胞外基質[6]。在發育過程中，立毛肌與毛囊和交感神經的關系十分密切。立毛肌附著于毛囊干細胞所在的毛囊隆起區域，交感神經使立毛肌平滑肌細胞束化，交感神經與立毛肌和毛囊一起形成三系單元(tri-lineage unit)[7-8]。Shwartz等[9]對小鼠皮膚的研究進一步闡述了三系單元的作用機制：在發育過程中，毛囊干細胞會分泌音猬因子(Sonic hedgehog)指導立毛肌交感神經小生境的形成，進而控制成體毛囊再生。Torkamani等[10]研究發現，脫發患者的毛囊附著的立毛肌結構殘缺，表明立毛肌可以影響成熟毛囊的完整性。也有研究表明，皮膚組織中立毛肌肌動蛋白的存在有利于毛囊調控毛干的運動[11]。除此之外，Sato等[12]研究表明，立毛肌在毛發移植中與毛囊的再生有關。可見立毛肌對皮膚其他附屬器官的發育過程起重要作用。

藏羊是我國特有的地方綿羊品種，近些年來，已經有研究人員報道了關于胚胎期藏羊皮膚中毛囊、汗腺等皮膚附屬物發育的分子調控機制[13-15]，但目前還沒有關于調控藏羊胚胎期立毛肌發育分子機制的報道。本研究初步探索藏羊立毛肌發生及發育的過程，通過形態學變化明確立毛肌發生的關鍵時期，并通過轉錄組數據分析該時期的整體基因表達變化，探索影響立毛肌發生過程的動態基因調控網絡和潛在的候選基因。該結果將豐富有關立毛肌形態發生過程中組織學和分子變化的研究，為了解人類毛發發育的分子調控機制提供參考，對了解皮膚生物學復雜性具有重要的參考意義。

1 材料與方法

1.1 試驗材料

藏羊(西藏羊)取自青海當地屠宰場，采集30只75～110日齡健康情況良好的藏羊胚胎并解剖背部皮膚，將樣品平均分為兩部分，一部分立即放在4 ℃用4%的多聚甲醛固定，另一部分凍存在液氮中[13-14]。

1.2 綿羊背部皮膚石蠟組織切片及HE染色

為了確定藏羊皮膚中立毛肌的發生階段，將藏羊不同胚齡背部皮膚樣品用梯度酒精脫水，二甲苯透明，石蠟包埋并切成5 μm組織切片。然后將背部皮膚切片進行脫蠟和HE(蘇木精和伊紅)染色，并用中性樹脂固定。對染色的皮膚切片進行顯微鏡拍照記錄，根據形態學初步推斷樣品胚胎日齡。

1.3 綿羊背部皮膚組織轉錄組測序

挑選藏羊胚齡大約在75和85 d(E75和E85)各3個背部皮膚組織樣品，使用TRIzol試劑(美國Invitrogen公司)分別從6個樣品中提取總RNA，并使用帶有Bioanalyzer 2100系統(美國加利福尼亞州安捷倫科技公司)的RNA Nano 6000分析試劑盒評估RNA完整性。按照制造商的程序，使用用于Illumina的NEBNext9 UltraTM RNA文庫制備試劑盒(美國NEB)，在Novogene(中國北京)構建測序文庫，利用Agilent Bioanalyzer 2100系統評估文庫質量，在Illumina Hiseq平臺(Hiseq X ten)上進行150個堿基對配對末端序列的測序。

1.4 綿羊背部皮膚組織轉錄組測序數據分析

原始數據的上游處理均在服務器中進行，利用fastqc和trimglore軟件進行數據的質控和過濾；通過Hisat2[16](版本:2.1.0)軟件將clean reads映射到藏羊參考基因組(版本：Oarv3.1)；用Samtools軟件對比對結果進行格式轉換、排序；使用Htseq[17]軟件對clean reads 進行計數；用Stringtie[18](版本：1.3.5)軟件計算每個基因的FPKM值。使用R語言(版本：3.6.3)進行轉錄組下游數據分析；用DESeq2軟件包進行基因表達量差異分析[19]。對于生物學重復，將P<0.05和|log2(fold change)|≥1設置為差異表達基因的閾值，獲得差異表達的基因。使用bioDBnet網站進行基因ID的轉換。因為藏羊是非模式生物，先用AnnotationHub R包獲取綿羊注釋文件，再使用 clusterProfiler[20]R包對DEGs進行GO富集分析，差異表達基因顯著富集于校正后P<0.05的GO條目，使用KOBAS[21]網站進行KEGG通路中差異表達基因的富集分析。使用Cytoscape StringApp軟件[22]實現蛋白質相互作用網絡的可視化。

1.5 差異表達基因實時熒光定量PCR驗證

1.5.1 引物的設計與合成 隨機選出6個差異表達基因，以羊的GAPDH為內參基因進行實時熒光定量PCR(RT-qPCR)驗證。利用NCBI網站中的primer-BLAST在線工具設計引物(表1)，引物由生工生物工程(上海)股份有限公司合成。

表1 RT-qPCR引物序列信息

1.5.2 熒光定量PCR 每個時期的樣本設置3個技術性重復(n=3)。按照Vazyme公司的產品ChamQ Universal SYBR qPCR Master Mix說明書對目的基因進行RT-qPCR檢測。反應體系為：cDNA 6.6 μL，酶Mix 7.5 μL，上、下游引物各0.45 μL；反應程序：95 ℃預變性30 s；95 ℃變性10 s，62 ℃退火30 s，72 ℃延申30 s，40個循環；65～95 ℃，每5 s增加0.5 ℃，5 min制作熔解曲線。用Bio-Rad CFX Manager 3.1軟件收集處理數據，采用2-ΔΔCT方法對mRNA進行定量數據分析，并用t檢驗進行統計分析。

2 結 果

2.1 藏羊背部皮膚中立毛肌的形態發生及發育

本試驗選擇西藏地毯羊(一種典型的粗羊毛綿羊)為研究對象，詳細研究皮膚中立毛肌的發生過程。通過組織學HE染色觀察到立毛肌從無到有直至結構完整的3個時期(圖1A～C)。根據Rogers[2]的描述確定TF2a(約E75，圖1A，T代表藏羊)和TF3a(約E85，圖1B)為立毛肌開始形成的關鍵時期。

A.TF2a(E75)皮膚組織切片圖：沒有觀察到明顯的立毛肌；B.TF3a(E85)皮膚組織切片圖：顯示形成早期立毛肌；C.TF4(E100)皮膚組織切片圖：顯示立毛肌的延伸。B、C圖紅色虛線內標識部位為立毛肌A. TF2a (E75) skin tissue section image: no obvious arrector pili muscle; B. TF3a (E85) skin tissue section image: the early formation of arrector pili muscle; C. TF4 (E100) skin tissue section image： the down-growth of arrector pili muscle. In B, C figures,the parts in the red dashed line are the arrector pili muscle圖1 藏羊胚胎期皮膚組織切片顯示立毛肌的早期形態發生過程圖Fig.1 The histological identification of arrector pili muscle development in Tibetan sheep skin during embryo period

2.2 藏羊背部皮膚組織轉錄組測序質量分析

取兩個時期(TF2a、TF3a)，每個時期3個樣品，共6個藏羊背部皮膚組織樣品進行RNA測序。將數據上傳于服務器中進行質控、比對、計數。每個樣品產生超過12G的原始數據，測序結果如表2所示，其中Q30數據基本都在90%以上，表明測序數據滿足后續分析的條件。

表2 轉錄組測序數據質量檢測分析

2.3 藏羊背部皮膚組織轉錄組數據差異基因表達分析

利用R軟件包DESeq2進行差異分析得到初步差異分析的結果。設置篩選差異基因的標準為|log2(fold change)|≥1 和P<0.05，得到的差異基因結果如圖2所示，共獲得1 159個差異表達基因，其中900個基因上調表達，259個基因下調表達。

虛線左上側代表顯著下調基因，虛線右上側代表顯著上調基因。TF2a代表E75；TF3a代表E85The upper left side of the dotted line represents significantly down-regulated genes, and the upper right side of the dotted line represents significantly up-regulated genes. TF2a stands for E75； TF3a stands for E85圖2 立毛肌早期形態發生期藏羊皮膚中差異表達基因火山圖(立毛肌雛形 vs. 無立毛肌)Fig.2 Volcano map of differentially expressed genes in Tibetan sheep skin during early morphogenesis of arrector pili muscle (the germ of arrector pili muscle vs. no arrector pili muscle)

2.4 差異表達基因的GO和KEGG富集分析

對1 159個差異基因進行GO和KEGG功能富集分析，一共顯著富集到71個GO條目，包括29個生物學過程(BP)，27個細胞組分(CC)，15個生物學功能(MF)。前30條顯著富集到的GO條目如圖3A所示，在生物學過程中，差異基因顯著富集到肽生物合成過程(peptide biosynthetic process，35個基因)、酰胺生物合成過程(amide biosynthetic process，38個基因)、翻譯(translation，34個基因)。在生物學功能歸類中，差異基因顯著富集到結構分子活性(structural molecule activity，40個基因)。所有富集條目大部分都與轉錄和翻譯過程緊密相關，表明細胞在此時期的活動，如細胞增殖和分化等非常活躍，此時期可能為器官發育形成的重要時期。

在KEGG富集分析中，共顯著富集到28條信號通路，其中圖3B列出了前20個顯著富集到的通路，顯著富集到的通路為代謝途徑(metabolic pathways，77個基因)、神經營養蛋白信號通路(neurotrophin signaling pathway，11個基因)、ECM-受體相互作用(ECM-receptor interaction，9個基因)等。

2.5 RT-qPCR驗證

對隨機選取的6個差異基因進行RT-qPCR驗證，結果如圖4所示，2個基因(MATN4、LAMB3)上調表達，4個基因(JUN、CNTFR、NOTCH3、COL4A2)下調表達，6個基因在立毛肌發育的2個時期(TF3avs. TF2a)的藏羊皮膚組織中表達趨勢與轉錄組測序結果一致。進一步表明測序結果具有可靠性。

2.6 調控立毛肌發生的候選基因篩選

立毛肌發生的初期伴隨著毛囊的成熟和汗腺的出現，這些器官的協同變化表明，一些參與調控毛囊或腺體發育的轉錄因子也可能會對立毛肌的發育有影響。本研究富集到了與管發育、腺體發育、肌肉發育、干細胞發育以及皮膚發育相關的GO條目中共23個非冗余差異基因，見表3。這些條目均屬于生物學過程類別。

KEGG富集分析得到了可能與立毛肌發育有關的通路基因(表4)。這些通路包括粘著斑(focal adhesion)信號通路、Wnt信號通路(Wnt signaling pathway)、Notch信號通路(Notch signaling pathway)、TGF-β信號通路(TGF-beta signaling pathway)和鈣離子信號通路(calcium signaling pathway)。

通過篩選GO和KGEE富集到的條目以及信號通路，共得到47個非冗余的差異基因(圖5A)，并對其進行了蛋白網絡互作分析。共有29個基因之間有互作關系，其中有16個上調基因，13個下調基因。這些基因構成復雜的調控網絡，例如，Wnt信號通路中的LEF1基因參與到肌肉結構發育和腺體發育等生物學過程；粘著斑信號通路中的JUN基因參與到腺體發育以及管發育等生物學過程(圖5B)。

3 討 論

組織切片結果展示出立毛肌發育過程中3個比較重要的發育時期的動態變化。根據已有文獻的描述[2]以及試驗觀察[13-14]，推測藏羊TF2a和TF3a兩個時期分別對應E75和E85左右，TF4則對應E100左右的胚齡。其中，E75和E85左右時期立毛肌分別處于未發育和發育初期階段，在E100左右時立毛肌已經發育較為成熟。因此，E75和E85藏羊立毛肌組織形態學的明顯差異為后續的轉錄組數據分析提供了基礎。

KEGG富集分析得到了候選信號通路。粘著斑通路(focal adhesion)可以與肌動蛋白細胞骨架途徑協同調節平滑肌的遷移[23]，粘著斑是連接細胞骨架與細胞外基質的一種細胞結構，是信號傳導過程中許多信號分子組裝的平臺[24]，差異基因在該信號通路富集表明，在立毛肌細胞形成后向真皮外基質方向擴展的過程中，該信號通路中某些基因的表達可能促進了立毛肌細胞的分化與形態的變化。Wnt分泌型糖脂蛋白家族介導的信號傳導通路在胚胎發育和組織穩態過程中起著調節細胞增殖、分化、遷移和極性等重要作用[25]，它被證明是調節癌癥的重要信號通路之一[26-27]，Wnt信號通路還廣泛參與到皮膚組織與毛囊等器官的發育過程并調節該程中的細胞穩態[28-29]。Notch信號通路是一種進化上保守的細胞內途徑，參與細胞的增殖、分化、遷移和細胞命運等過程，該途徑不僅參與皮膚組織器官的發育，而且與人類皮膚疾病的發生密切相關[30]，另外它也參與了腸道平滑肌的發育[31]。轉化生長因子β(TGF-β) 信號通路在動物中高度保守，被認為是多細胞動物進化初期最早出現的信號通路之一，其在胚胎期可以介導組織特異性的分化和增殖等活動[32]，TGF-β 信號通路在調節肌肉的生長和重塑方面也起著重要作用[33]。Ca2+信號通路則是促進肌肉形成、肌肉動態平衡和再生的信號傳導網絡的重要組成部分[34]。上述信號通路被富集表明，TF2a到TF3a時期這一發育過程中與肌肉發育相關的基因活躍表達，各種信號通路之間的相互通訊參與了肌肉細胞發育的調節，立毛肌展現出的顯著形態變化也證實了這一點。

BP、MF、CC分別代表生物學過程、生物學功能、細胞組分BP, MF, and CC represent biological process, molecular function, and cellular component, respectively圖3 立毛肌早期形態發生期差異表達基因富集分析中前30個GO條目(A)與前20條KEGG通路圖(B)(TF3a vs. TF2a)Fig.3 Representation of top 30 GO terms(A) and top 20 KEGG pathways(B) in the enrichment analysis of differentially expressed genes in early morphogenesis of arrector pili muscle(TF3a vs. TF2a)

*. P<0.05圖4 RT-qPCR檢測6個差異表達基因在藏羊皮膚組織中的表達模式(TF3a vs. TF2a)Fig.4 RT-qPCR verification of expression tendency of 6 differentially expressed genes in Tibetan sheep skin(TF3a vs. TF2a)

表3 參與立毛肌發育的候選GO條目及基因

表4 參與立毛肌發育的候選信號通路及基因

在47個候選差異表達基因中，有研究表明一些基因與肌肉發育相關。其中，SPP1在立毛肌發生期間高表達，有文獻報道SPP1與人類肌肉炎癥以及肌肉再生有關[35]，它也可能參與調控雞成肌細胞的增殖發育過程[36]，表明SPP1可能會在立毛肌發生早期誘導肌肉發育。BAMBI基因在立毛肌發生期間表達量顯著上調，Yao等[37]的研究也證明，BAMBI蛋白在肌肉生長和再生期間富集在細胞膜上，這表明BAMBI介導的重要信號通路可能是肌肉生長和再生的重要組成部分。同樣，作為表達量顯著上調的TNNI3基因，在動物心肌形成過程中發揮著重要的作用，是影響肌纖維類型的一種肌小節基因，它的突變常常會引發肥厚型心肌病的發生[38]，該基因的高表達也暗示了其可能與立毛肌的發育密切相關。HOX基因家族是胚胎發生過程中調控干細胞和組織模式的關鍵調節因子，本研究結果表明，HOXA9基因在E85天時顯著下調，這與Schw?rer等[39]的試驗結果一致，敲低HOXA9基因的小鼠呈現肌纖維顯著增加，通過抑制HOXA9基因的表達可以使肌肉再生，表明HOXA9基因的下調促進了肌肉干細胞的分化。SOX15基因在E85天時顯著上調，Lee等[40]的研究表明，SOX15基因對肌肉的再生有著重要的作用，Savage等[41]的研究顯示，SOX15基因可能是早期肌肉前體細胞的調節因子，它的表達對肌肉發生是正向調控作用，表明SOX15在E85天時的上調可能促進立毛肌的發生過程。因此，推測SPP1、BAMBI、TNNI3、HOXA9、SOX15為影響藏羊立毛肌發生的候選基因。

黑色代表表達量上調的基因；淺灰色代表表達量下調的基因Black signal represents genes with up-regulated expression; Light gray signal represents genes with down-regulated expression圖5 差異表達基因聚類熱圖(A)及網絡互作圖(B)Fig.5 Heat map of clustering (A) and protein-protein network interaction (B)of differentially expressed genes

4 結 論

本研究利用組織切片和HE染色技術推測出藏羊胚胎期立毛肌發生的關鍵時期為E75～E85，通過皮膚轉錄組數據分析，得到1 159個差異表達基因，其中有47個可能參與到肌肉發育過程，進一步篩選得到5個與立毛肌發育顯著相關的候選基因，分別為SPP1、BAMBI、TNNI3、HOXA9、SOX15。該結果為進一步闡述藏羊立毛肌發生的分子作用機制以及皮膚調控分子網絡奠定了基礎。