天祝白牦牛退行期毛囊細胞亞群鑒定以及特征基因生物信息學分析

鄭青波 葉娜 張嘵蘭 包鵬甲 王福彬 任穩穩 廖月姣閻萍 潘和平

（1.西北民族大學生命科學與工程學院，蘭州 730030；2.中國農業科學院蘭州畜牧與獸藥研究所 甘肅省牦牛繁育工程重點實驗室，蘭州 730050）

天祝白牦牛是甘肅省天祝藏族自治縣一種獨有的牦牛地方品種，主要集中在海拔3 000 m以上的高寒草原上，有著草原“白珍珠”和祁連“雪牡丹”的美譽［1］。牦牛能在高寒高海拔地帶生存得益于獨特的混合型毛被結構，周身被毛可分為粗毛、絨毛和兩型毛。天祝白牦牛毛具有含絨率高、可紡性好以及被毛純白等特點，因而被稱為“雪絨”產品［2］，是一種獨特的毛紡工業原料資源，具有較高的商業價值，備受國內外市場青睞。毛囊隸屬于皮膚的亞器官是產生毛發的基本單位。哺乳動物出生后毛囊呈周期性生長發育，分為生長期、退行期和休止期。毛囊周期性生長使牦牛能夠在秋冬時節長出絨毛來抵御高原高寒，在春夏時節褪去絨毛給機體增加散熱，這對其適應高寒環境具有重要意義。毛囊的周期性循環需要不同類型細胞相互協調和諸多信號通路的調控。退行期的毛發停止生長，毛囊開始皺縮，與生長期相比，毛囊形態、細胞類群及相關信號分子的表達和互作也發生巨大變化。探究退行期毛囊細胞類群和基因表達情況，對牦牛毛囊發育過程中不同細胞類型調控規律的研究具有重要意義。

單細胞轉錄組測序（single cell RNA sequencing，scRNA-seq）通過二代測序技術在單個細胞水平上對轉錄組進行測序，能夠分析單個細胞間遺傳和基因表達水平的差異，解決組織和器官發育過程中細胞異質性難題。轉錄組測序直觀地反應了單個細胞的基因表達情況，能更好地揭示細胞種類、亞型和狀態［3-4］。北京大學湯富酬教授2009年在Nature Methods上首次報道了單細胞轉錄組測序的研究，從而揭開了單細胞轉錄組測序研究時代的序幕［4］。近幾年scRNA-seq技術以及下游數據生物信息學的迅速發展，使其在單個試驗中可同時分析上萬個細胞的轉錄組信息，并已成為胚胎發育，器官分化，癌癥發生以及構建細胞圖譜等研究最有力的研究手段。為探究天祝白牦牛毛囊在退行期發育過程中不同類型細胞的功能，本研究采用單細胞轉錄組技術對天祝白牦牛退行期毛囊進行異質性分析，利用已知標記分子對細胞類群進行篩選鑒定，并對鑒定得到的細胞類群進行GO等分析，研究結果將為進一步探究牦牛毛囊發育規律以及發育調控機制提供理論基礎。

1 材料與方法

1.1 材料

皮膚樣品來自甘肅省武威市天?？h天祝白牦牛育種場。

主要試劑：中性樹脂、二甲苯、DAPI、BSA、DynabeadsTMMyOneTMSILANE PN-2000048、Chromium i7 MuLtiplex Kit、ChromiumTMSingle Cell 3′Library& Gel Bead Kit v3、ChromiumTMSingle Cell A Chip Kit等。

主要儀器：包埋機、病理切片機、熒光掃描顯微鏡、10×Genomics Chromium 微流控平臺（10×Genomics）；10×Genomics Chromium Controller（10×Genomics）；Illumina NovaSeq 6000；Agilent 2100 Bioanalyze（Agilent）等。

1.2 方法

1.2.1 皮膚組織切片與HE染色 將皮膚表面的絨毛與粗毛除凈，切成較小的組織塊，在用固定液固定24 h以上。將目的部位組織修理平整，放進脫水機內用梯度酒精進行脫水浸蠟。再將融化的蠟放入包埋框內，在蠟凝固之前將組織取出，放-20℃條件下冷卻。用切片機將組織塊切成4 μm的薄片，把切片攤于溫水上，用載玻片將組織撈起，并在60℃條件下進行烤片。待蠟烤化后進行脫蠟、伊紅染色和中性樹膠封片，等晾干后可進行觀察。

1.2.2 免疫熒光組化分析 先將組織切片進行脫蠟，再用抗原修復緩沖液進行抗原修復，用組化筆畫圈，甩干PBS，滴加3%的BSA進行封閉。在切片上滴加稀釋好的一抗，然后放入濕盒內4℃孵育過夜，將切片用PBS洗滌，待切片稍干后滴加二抗，進行避光孵育，最后滴加DAPI進行染核，用抗熒光淬滅封片劑封片。

1.2.3 單細胞懸液制備 選擇2歲左右健康的天祝白牦牛3只，均為母牛。采集退行期肩胛處皮膚組織，清洗消毒后置于組織保護液中帶回實驗室。在超凈工作臺繼續用PBS和75%酒精清洗消毒，用0.25%的dispase 酶消化2 h，拔出單根毛發置于培養皿中，將胰蛋白酶消化后的毛發放在顯微鏡下觀察細胞是否脫落。后經反復吹打、過濾、離心和重懸獲得細胞懸液，最后進行計數和單細胞活力檢測。

1.2.4 10 ×Genomics單細胞測序 由北京貝瑞和康生物技術有限公司完成單細胞文庫的構建和測序?；?0×Genomics ChromiumTM系統進行細胞標記，含有Barcode 信息的Gel Beads 與細胞的油珠混合形成GEMs。每個GEM 中，細胞破裂后釋放的mRNA被反轉錄成帶有Barcode的cDNA。然后將所有細胞的cDNA 收集到一起，擴增后按Illumina 測序文庫構建標準流程進行測序文庫構建。再用 Qubit 2.0進行初步定量。然后用Agilent 2100對文庫中的Insert DNA大小進行檢測。Insert DNA大小符合預期后，使用QPCR方法對文庫的有效濃度進行準確定量。庫檢合格后，使用Illumina NovaSeq 6000平臺進行PE150測序。

1.2.5 測序數據質控與下游分析 對單細胞測序得到的原始序列進行過濾。去除N的個數大于3和含adaptor的Reads。檢測有無AT、GC分離現象。使用10×Genomics官方分析軟件CellRanger進行細胞檢測，用CellRanger中的STAR對 Reads 進行參考基因組剪接識別比對。通過CellRanger糾正UMI序列中的測序錯誤。CellRanger完成后就可以對scRNA-seq數據進行聚類分析。本研究通過R語言Seurat包對低質量細胞進行過濾剔除，再用NormalizeData函數對數據進行標準化，之后用降維算法將表達矩陣降低到重要特征值，使用主成分分析將表達矩陣的維度從Cells×Genes降低到cell×M，同時將數據傳遞到具有可視化的UMAP中，在最大限度保留原始數據特征的同時降低數據維度，在質量控制和scRNA-seq數據標準化后，得到聚類的細胞cluster結果。

1.2.6 細胞類群鑒定 細胞類群根據已知marker基因進行鑒定，如表1所示。

表1 細胞類型標記分子Table 1 Cell type marker molecules

1.2.7 GO富集分析 GO 富集分析采用Metascape（http://metascape.org/gp/index.html#/main/step1）對差異基因進行 GO 富集分析，同時用R語言將調控網絡的基因名稱轉化為基因的ID，應用ClusterProfiler包、Ggplot2包、Enrichplot包從生物學途徑（biological process，BP）、細胞組成（cellular component，CC）、分子功能（molecular function，MF）3個方面對潛在通路進行GO富集分析，并將BP、CC、MF前10條進行可視化。

1.2.8 KEGG與基因互作分析 運用Metascape篩選出P-adjust值＜0.01，細胞相關通路頻次≥3的特征基因進行KEGG信號通路分析和基因互作分析。

2 結果

2.1 毛囊結構

通過HE染色可以直觀的觀察到毛囊結構情況，確定毛囊發育時期。根據毛囊的形態可以清晰的觀察到天祝白牦牛毛囊IRS、Bulb、ORS、HS等基本結構。退行期毛干部分脫落，毛囊群開始變得不完整，部分細胞已經凋亡，毛乳頭和基質減少，毛球萎縮變小，同時毛囊中出現空腔如圖1所示。

圖1 天祝白牦牛皮膚組織結構圖Fig.1 Skin tissue structure of Tianzhu white yak

2.2 單細胞數據質控

通過scRNA-seq技術分析天祝白牦牛毛囊退行期肩胛部皮膚單細胞基因表達情況，獲得退行期毛囊大約12 000個單細胞轉錄組信息。比對到基因組上的比率為93.3%，基因中位數為927，測序飽和度為70.7%，平均Reads數為62 282。數據整體質量可靠。使用 Read 10×函數讀取細胞表達矩陣，創建Seruratd對象，用FilterCells函數過濾低質量細胞，僅保留至少在3個細胞中表達的基因以及表達量不低于 200 個基因的細胞，本試驗根據樣本中每個細胞的總UMI數目、表達基因數目以及線粒體基因表達占比分布（圖2-A），對數據進行進一步篩選，保留基因數在200-4 000且線粒體基因表達量小于5%的細胞，用于后續分析。同時還分析了基因數與UMI數相關性，相關系數為0.92，說明UMI數越高，基因數也就越高，相關性符合預期（圖2-B）。

圖2 檢測到的細胞中基因數與UMI總數統計以及相關性分析Fig.2 Statistics and correlation analysis of gene number and total UMI in detected cells

2.3 高變基因分析

為保證數據質量，利用Seurat選擇前2 000個高變基因進行后續的基因表達矩陣降維和細胞聚類分析。圖3展示了所有細胞中基因表達豐度的變異情況，并挑選出在細胞間變異系數最大的前10個基因，發現KRT85、KRT25、PMP2、KLK7等基因在毛囊退行期變異系數較大。

圖3 基因離散度分布圖Fig.3 Gene dispersion distribution map

2.4 主成分分析

通過正交變換將一系列可能存在線性相關的變量轉換為一系列線性不相關的新變量，再將新變量放在更小的維度下展示數據的特征，并用散點圖的形式進行多維數據的可視化。結果如圖4-A所示，每一個點代表一個細胞所在的位置，聚集在一起的點代表這些細胞具有一定的相關性，距離越近說明其相關程度越高。挑選出已知在細胞間變異系數最大的前4個基因，觀察其在散點圖中的分布情況（圖4-B）。挑選前20個主成分進行計算，之后進行UMAP聚類分析（圖4-C）。

圖4 主成分分析Fig.4 Principal component analysis

2.5 聚類與主要細胞類型鑒定

Seurat 中的表達矩陣使用 PCA分析之后，再采用 UMAP 進行降維，得到了14個細胞cluster（圖5-A），根據已知的標記基因對每個 cluster 的細胞類型進行鑒定（圖5-B），結果顯示cluster 0，3，8，9表達表皮細胞（epidermal）的標記基因 KRT14，KRT15；cluster1，4，10表達濾泡間上皮分化細胞（interfollicular epidermis differentiated cell，IFE-DC）標記基因 KRT10，KRTDAP；cluster2，6 表達毛干（hair shaft，HS）細胞標記基因 SPARC，LHX2；cluster5表達內根鞘細胞（inner root sheath，IRS）標記基因KRT28，KRT71；cluster7表達漏斗細胞（infundibulum，INFU）標記基因TOP2A，UBE2C；cluster13表達黑素細胞標記基因PMEL，TYRP1（圖5-C）。其中特征基因 KRT71、KRT28、TOP2A、UBE2C、TYRP1、PMEL在不同細胞類型中的表達水平如圖5-D所示。

圖5 UMAP聚類與主要類型細胞鑒定Fig.5 UMAP clustering and identification of major cell types

2.6 主要細胞類型的分子特征

在對單細胞數據進行不同類型聚類之后，利用點狀圖進一步對每個cluster特異性表達的基因進行比較分析。結果如圖6所示，表皮細胞高表達 KRT14、KRT15、KRT5；IFE-DC高 表 達 KRT1、KRT10、KRTDAP、SBSN；Hair shaft高表達SPARC、LHX2、CXCL14、LGR5；IRS高 表 達 KRT28、KRT-27、KRT71、DLX3；INFU高表達TOP2A、UBE2C、CKAP2、CENPE、DLX3；黑色素細胞高表達PMEL、TYRP1、PECAM1、RGS1。

圖6 不同細胞類群差異基因比較分析Fig.6 Comparative analysis of differential genes in different cell clusters

2.7 不同細胞類型富集分析

對聚類之后主要細胞類群進行富集分析，結果發現表皮細胞聚類得到0、3、8、9共4個cluster，特異表達基因分別為460、408、655、690個。GO富集結果顯示，4個亞群主要富集在GO：0008544表皮發育、GO：0030855上皮細胞分化、GO：0009611損傷應答、GO：0000904調控細胞形態發生等生物過程，而circos圖也顯示，這4個亞群之間存在密集的共同信號通路（圖7-A）；IFE-DC細胞聚類得到1、4、10共3個cluster，分別有582、1 061、341個特異表達基因，GO富集結果顯示：這3個亞群主要富集在GO：0008544表皮發育、GO：0097435超纖維組織發育、GO：0030155細胞黏附調節以及GO：0048729組織形態發生等生物過程中，同時circos圖顯示，這3個亞群之間存在許多相同基因以及富集通路（圖7-B）。

圖7 表皮細胞系與IFE-DC細胞不同cluster特征基因富集分析Fig.7 Enrichment analysis of different clusters characteristic genes in epidermal cell lineage and IFE-DC

HS聚類得到2、6共2個cluster，分別有480、669個特異表達基因，通過GO富集分析發現，這2個亞群主要富集在GO：0008544表皮發育、GO：0034330細胞連接、GO：0000165 MAPK級聯激活以及GO：0048729組織形態發生等生物過程（圖8-A）。對黑色素細胞進行分析發現了752個差異基因，GO富集結果顯示黑色素細胞主要富集在GO：0030155細胞黏附的調控、GO：0030855上皮細胞分化、GO：0030335細胞遷移的正向調控以及GO：0030029肌動蛋白等生物過程（圖8-B）。

圖8 HS與黑色素細胞特征基因富集分析Fig.8 Enrichment analysis of HS and melanocytes characteristic genes

對INFU進行分析發現，其BP共有486個條目，主要與表皮發育、上皮細胞增殖、肌動蛋白絲以及細胞-底物黏附等有關；CC分析中共有82個條目，主要包括細胞-底物連接、黏著斑通路、細胞間連接以及細胞皮層；MF分析中共有27個條目，主要富集在鈣粘蛋白結合、肌動蛋白結合以及肌動蛋白絲結合等條目中（圖9-A）。對IRS進行分析發現，其BP共有733個條目，主要與核糖核蛋白復合體的生物發生、ncRNA代謝過程和表皮發育等相關；CC分析中共有104個條目，特征基因主要富集在細胞-底物連接、黏著斑通路等GO通路；MF分析中共有63個條目，主要富集于轉錄核心調節因子活性和鈣黏蛋白結合等條目中（圖9-B）。

圖9 INFU細胞與IRS細胞特征基因的GO富集分析Fig.9 GO enrichment analysis of characteristic genes in INFU and IRS cells

通過KEGG分析發現，Epidermal特征基因主要富集于ko03010：核糖體、ko04520：黏附連接等通路；IFE-DC特征基因主要富集于ko03050：蛋白酶體等通路；HS特征基因主要富集于ko04530：緊密連接、ko04360：軸子引導等通路；IRS特征基因主要富集于ko03030：DNA復制、ko03040：剪接體等通路；INFU特征基因主要富集于ko03050：蛋白酶體、ko03013：RNA轉運、ko04110：細胞周期、ko03030：DNA復制、ko03040：剪接體等通路；Melanocyte特征基因主要富集于ko03010：核糖體、ko04360：軸子引導信號通路等（圖10-A）。運用metascape對特征基因進行互作分析，獲得不同類型細胞基因互作網絡圖，共同靶基因有28個，分別為FKBP4、KRT1、KRT80、KRT17、KRT10、RPL15、RPL8、RPL14、NSA2、HSPA8、GJA1、CLTB、LGALS1、CST3、DSG2、PERP、PKP1、DSC2、JUP、DSG1、DSP、LGALS3、COL17A1、FGFR2、NCK1、GGT6、LY6G6C、CD109（圖10-B）。

圖10 不同細胞類群特征基因KEGG富集分析Fig.10 KEGG Enrichment analysis of characteristic gene in different cell clusters

2.8 熒光免疫組化分析

利用熒光免疫組化對天祝白牦牛退行期皮膚組織進行了染色分析如圖11所示。結果顯示，SPARC陽性的細胞主要集中在HS中，其次在表皮細胞系以及內根鞘中也存在一定的表達；KRT10主要在表皮中高表達，而在真皮細胞系中呈現出弱表達；此外，從皮膚組織中的表達情況中可以看出DLX3在毛干周圍的內根鞘和外根鞘中高表達。對關鍵標記基因進行免疫組化分析得到的實驗結果與上述分析結果基本一致。

圖11 關鍵標記基因熒光免疫組化分析Fig.11 Fluorescence immunohistochemical analysis of key marker genes

3 討論

毛囊作為皮膚組織中的微型器官，與微環境的其他組織、細胞共同的受到諸多信號通路調控［9］。毛囊的發育與成熟、毛囊的周期性循環均依賴于真皮與表皮的相互作用［10］。真皮細胞發出的啟動信號使相鄰部位表皮增厚，表皮細胞發出的信號促使真皮細胞在基板下聚集，誘發基板下的真皮細胞發育為毛乳頭，并由毛乳頭誘使毛囊內表皮細胞繼續增殖分化，形成毛干并突出皮膚生長［11］。在退行期，毛乳頭分泌的調控因子減少，毛母質細胞萎縮，毛囊再生部分退化，毛球部與外根鞘內的部分細胞凋亡，到退行期后期基質消失［12-13］。因此利用scRNA-seq對細胞進行生物信息學分析，可為下一步試驗提供指導。

單細胞轉錄組測序一次可以檢測產生上萬個細胞的轉錄信息，因此在對細胞類型進行注釋之前需要進行降維。主成分分析是一種利用投影矩陣將高維數據投影到低維空間的線性降維方法，可將多個變量轉換為少數幾個反映絕大部分信息的綜合變量即主成分［14］。tSNE是一種非線性降維方法，其巧妙的解決了集群表示過度擁擠問題。但有研究發現該算法比較消耗計算機資源，易丟失集群間關系而無法有效地表示數量龐大的數據集［15-16］。UMAP是建立在黎曼幾何和代數拓撲理論框架上的一種可視化降維方法，并對嵌入維數沒有計算限制。Satpathy等［17］發現UMAP在單細胞轉錄組測序中的可視化過程中效果理想。

Ryan等［18］發現Rho家族GTPases能夠通過調節鈣粘蛋白、整合素和黏附結構來對表皮功能進行調控。Fuchs［19］發現表皮必須通過增殖、分化和凋亡不斷補充和維持自身穩態。正常表皮功能的主要決定因素是細胞黏附，除了維持表皮結構外，細胞-基質和細胞-細胞黏附在調節信號轉導和表皮層的細胞功能方面也起著關鍵作用［20］。葛偉［21］對陜北白絨山羊毛囊發育過程中的表皮細胞進行GO富集分析發現，表皮細胞聚類得到的5個cluster共同富集了表皮發育、上皮細胞分化等信號通路。這與本研究中牦牛表皮細胞富集結果相似，說明牦牛與其他動物毛囊發育存在相似性。葉娜［22］對聚類得到的IFE-DC細胞進行GO富集分析發現，IFE-DC細胞主要富集在表皮發育、超纖維組織發育以及創傷反應調控等通路。與本研究IFE-DC細胞富集分析結果一致。

MAPK級聯是毛干細胞主要富集通路之一。而MAPK信號通路能誘導毛囊細胞的增殖分化，促進毛囊的周期性發育。有研究顯示，Tβ4能夠影響p-catenin、VEGF以及MMP2的表達，使MAPK通路中的P38被激活，進而影響絨毛生長、毛囊分布和毛干數量［23］。在本研究中，黑色素細胞特征基因主要富集在軸突引導、細胞黏附的調控、細胞遷移的正向調控以及肌動蛋白等通路中。有研究發現黑色素在黑色素細胞內合成后，通過黑色素細胞的樹突結構以黑色素顆粒的形式轉運至周邊表皮，黑色素轉運過程需要微管驅動蛋白以及動力蛋白等的參與［24］。高澤成［25］利用HE染色發現天祝白牦牛毛囊毛球部黑色素細胞無法生成黑色素顆粒，導致天祝白牦牛毛纖維中不存在黑色素顆粒而呈現白色。在毛囊形態發生期間，黑色素細胞前體遷移到發育的毛囊中，并產生去分化的黑色素細胞，該黑色素細胞能夠主動產生色素并將色素運輸到角質細胞中［26］。但黑色素細胞在軸突引導方面的生物學功能仍需進行進一步驗證。

4 結論

對天祝白牦牛毛囊退行期進行分析，鑒定出IFE-DC細胞、表皮細胞以及INFU細胞等細胞類型，特征基因主要參與表皮發育、上皮細胞分化以及細胞形態學發生等生物學過程，KEGG分析顯示特征基因顯著富集于黏附連接、細胞周期以及RNA轉運等通路中，同時發現退行期主要細胞類型擁有GJA1、KRT1、KRT80、FGFR2等28個共同基因。