鄭青波 葉娜 張嘵蘭 包鵬甲 王福彬 任穩(wěn)穩(wěn) 廖月姣閻萍 潘和平
(1.西北民族大學生命科學與工程學院,蘭州 730030;2.中國農(nóng)業(yè)科學院蘭州畜牧與獸藥研究所 甘肅省牦牛繁育工程重點實驗室,蘭州 730050)
天祝白牦牛是甘肅省天祝藏族自治縣一種獨有的牦牛地方品種,主要集中在海拔3 000 m以上的高寒草原上,有著草原“白珍珠”和祁連“雪牡丹”的美譽[1]。牦牛能在高寒高海拔地帶生存得益于獨特的混合型毛被結(jié)構(gòu),周身被毛可分為粗毛、絨毛和兩型毛。天祝白牦牛毛具有含絨率高、可紡性好以及被毛純白等特點,因而被稱為“雪絨”產(chǎn)品[2],是一種獨特的毛紡工業(yè)原料資源,具有較高的商業(yè)價值,備受國內(nèi)外市場青睞。毛囊隸屬于皮膚的亞器官是產(chǎn)生毛發(fā)的基本單位。哺乳動物出生后毛囊呈周期性生長發(fā)育,分為生長期、退行期和休止期。毛囊周期性生長使牦牛能夠在秋冬時節(jié)長出絨毛來抵御高原高寒,在春夏時節(jié)褪去絨毛給機體增加散熱,這對其適應(yīng)高寒環(huán)境具有重要意義。毛囊的周期性循環(huán)需要不同類型細胞相互協(xié)調(diào)和諸多信號通路的調(diào)控。退行期的毛發(fā)停止生長,毛囊開始皺縮,與生長期相比,毛囊形態(tài)、細胞類群及相關(guān)信號分子的表達和互作也發(fā)生巨大變化。探究退行期毛囊細胞類群和基因表達情況,對牦牛毛囊發(fā)育過程中不同細胞類型調(diào)控規(guī)律的研究具有重要意義。
單細胞轉(zhuǎn)錄組測序(single cell RNA sequencing,scRNA-seq)通過二代測序技術(shù)在單個細胞水平上對轉(zhuǎn)錄組進行測序,能夠分析單個細胞間遺傳和基因表達水平的差異,解決組織和器官發(fā)育過程中細胞異質(zhì)性難題。轉(zhuǎn)錄組測序直觀地反應(yīng)了單個細胞的基因表達情況,能更好地揭示細胞種類、亞型和狀態(tài)[3-4]。北京大學湯富酬教授2009年在Nature Methods上首次報道了單細胞轉(zhuǎn)錄組測序的研究,從而揭開了單細胞轉(zhuǎn)錄組測序研究時代的序幕[4]。近幾年scRNA-seq技術(shù)以及下游數(shù)據(jù)生物信息學的迅速發(fā)展,使其在單個試驗中可同時分析上萬個細胞的轉(zhuǎn)錄組信息,并已成為胚胎發(fā)育,器官分化,癌癥發(fā)生以及構(gòu)建細胞圖譜等研究最有力的研究手段。為探究天祝白牦牛毛囊在退行期發(fā)育過程中不同類型細胞的功能,本研究采用單細胞轉(zhuǎn)錄組技術(shù)對天祝白牦牛退行期毛囊進行異質(zhì)性分析,利用已知標記分子對細胞類群進行篩選鑒定,并對鑒定得到的細胞類群進行GO等分析,研究結(jié)果將為進一步探究牦牛毛囊發(fā)育規(guī)律以及發(fā)育調(diào)控機制提供理論基礎(chǔ)。
皮膚樣品來自甘肅省武威市天??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 皮膚組織切片與HE染色 將皮膚表面的絨毛與粗毛除凈,切成較小的組織塊,在用固定液固定24 h以上。將目的部位組織修理平整,放進脫水機內(nèi)用梯度酒精進行脫水浸蠟。再將融化的蠟放入包埋框內(nèi),在蠟凝固之前將組織取出,放-20℃條件下冷卻。用切片機將組織塊切成4 μm的薄片,把切片攤于溫水上,用載玻片將組織撈起,并在60℃條件下進行烤片。待蠟烤化后進行脫蠟、伊紅染色和中性樹膠封片,等晾干后可進行觀察。
1.2.2 免疫熒光組化分析 先將組織切片進行脫蠟,再用抗原修復緩沖液進行抗原修復,用組化筆畫圈,甩干PBS,滴加3%的BSA進行封閉。在切片上滴加稀釋好的一抗,然后放入濕盒內(nèi)4℃孵育過夜,將切片用PBS洗滌,待切片稍干后滴加二抗,進行避光孵育,最后滴加DAPI進行染核,用抗熒光淬滅封片劑封片。
1.2.3 單細胞懸液制備 選擇2歲左右健康的天祝白牦牛3只,均為母牛。采集退行期肩胛處皮膚組織,清洗消毒后置于組織保護液中帶回實驗室。在超凈工作臺繼續(xù)用PBS和75%酒精清洗消毒,用0.25%的dispase 酶消化2 h,拔出單根毛發(fā)置于培養(yǎng)皿中,將胰蛋白酶消化后的毛發(fā)放在顯微鏡下觀察細胞是否脫落。后經(jīng)反復吹打、過濾、離心和重懸獲得細胞懸液,最后進行計數(shù)和單細胞活力檢測。
1.2.4 10 ×Genomics單細胞測序 由北京貝瑞和康生物技術(shù)有限公司完成單細胞文庫的構(gòu)建和測序。基于10×Genomics ChromiumTM系統(tǒng)進行細胞標記,含有Barcode 信息的Gel Beads 與細胞的油珠混合形成GEMs。每個GEM 中,細胞破裂后釋放的mRNA被反轉(zhuǎn)錄成帶有Barcode的cDNA。然后將所有細胞的cDNA 收集到一起,擴增后按Illumina 測序文庫構(gòu)建標準流程進行測序文庫構(gòu)建。再用 Qubit 2.0進行初步定量。然后用Agilent 2100對文庫中的Insert DNA大小進行檢測。Insert DNA大小符合預期后,使用QPCR方法對文庫的有效濃度進行準確定量。庫檢合格后,使用Illumina NovaSeq 6000平臺進行PE150測序。
1.2.5 測序數(shù)據(jù)質(zhì)控與下游分析 對單細胞測序得到的原始序列進行過濾。去除N的個數(shù)大于3和含adaptor的Reads。檢測有無AT、GC分離現(xiàn)象。使用10×Genomics官方分析軟件CellRanger進行細胞檢測,用CellRanger中的STAR對 Reads 進行參考基因組剪接識別比對。通過CellRanger糾正UMI序列中的測序錯誤。CellRanger完成后就可以對scRNA-seq數(shù)據(jù)進行聚類分析。本研究通過R語言Seurat包對低質(zhì)量細胞進行過濾剔除,再用NormalizeData函數(shù)對數(shù)據(jù)進行標準化,之后用降維算法將表達矩陣降低到重要特征值,使用主成分分析將表達矩陣的維度從Cells×Genes降低到cell×M,同時將數(shù)據(jù)傳遞到具有可視化的UMAP中,在最大限度保留原始數(shù)據(jù)特征的同時降低數(shù)據(jù)維度,在質(zhì)量控制和scRNA-seq數(shù)據(jù)標準化后,得到聚類的細胞cluster結(jié)果。
1.2.6 細胞類群鑒定 細胞類群根據(jù)已知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語言將調(diào)控網(wǎng)絡(luò)的基因名稱轉(zhuǎn)化為基因的ID,應(yīng)用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,細胞相關(guān)通路頻次≥3的特征基因進行KEGG信號通路分析和基因互作分析。
通過HE染色可以直觀的觀察到毛囊結(jié)構(gòu)情況,確定毛囊發(fā)育時期。根據(jù)毛囊的形態(tài)可以清晰的觀察到天祝白牦牛毛囊IRS、Bulb、ORS、HS等基本結(jié)構(gòu)。退行期毛干部分脫落,毛囊群開始變得不完整,部分細胞已經(jīng)凋亡,毛乳頭和基質(zhì)減少,毛球萎縮變小,同時毛囊中出現(xiàn)空腔如圖1所示。
圖1 天祝白牦牛皮膚組織結(jié)構(gòu)圖Fig.1 Skin tissue structure of Tianzhu white yak
通過scRNA-seq技術(shù)分析天祝白牦牛毛囊退行期肩胛部皮膚單細胞基因表達情況,獲得退行期毛囊大約12 000個單細胞轉(zhuǎn)錄組信息。比對到基因組上的比率為93.3%,基因中位數(shù)為927,測序飽和度為70.7%,平均Reads數(shù)為62 282。數(shù)據(jù)整體質(zhì)量可靠。使用 Read 10×函數(shù)讀取細胞表達矩陣,創(chuàng)建Seruratd對象,用FilterCells函數(shù)過濾低質(zhì)量細胞,僅保留至少在3個細胞中表達的基因以及表達量不低于 200 個基因的細胞,本試驗根據(jù)樣本中每個細胞的總UMI數(shù)目、表達基因數(shù)目以及線粒體基因表達占比分布(圖2-A),對數(shù)據(jù)進行進一步篩選,保留基因數(shù)在200-4 000且線粒體基因表達量小于5%的細胞,用于后續(xù)分析。同時還分析了基因數(shù)與UMI數(shù)相關(guān)性,相關(guān)系數(shù)為0.92,說明UMI數(shù)越高,基因數(shù)也就越高,相關(guān)性符合預期(圖2-B)。
圖2 檢測到的細胞中基因數(shù)與UMI總數(shù)統(tǒng)計以及相關(guān)性分析Fig.2 Statistics and correlation analysis of gene number and total UMI in detected cells
為保證數(shù)據(jù)質(zhì)量,利用Seurat選擇前2 000個高變基因進行后續(xù)的基因表達矩陣降維和細胞聚類分析。圖3展示了所有細胞中基因表達豐度的變異情況,并挑選出在細胞間變異系數(shù)最大的前10個基因,發(fā)現(xiàn)KRT85、KRT25、PMP2、KLK7等基因在毛囊退行期變異系數(shù)較大。
圖3 基因離散度分布圖Fig.3 Gene dispersion distribution map
通過正交變換將一系列可能存在線性相關(guān)的變量轉(zhuǎn)換為一系列線性不相關(guān)的新變量,再將新變量放在更小的維度下展示數(shù)據(jù)的特征,并用散點圖的形式進行多維數(shù)據(jù)的可視化。結(jié)果如圖4-A所示,每一個點代表一個細胞所在的位置,聚集在一起的點代表這些細胞具有一定的相關(guān)性,距離越近說明其相關(guān)程度越高。挑選出已知在細胞間變異系數(shù)最大的前4個基因,觀察其在散點圖中的分布情況(圖4-B)。挑選前20個主成分進行計算,之后進行UMAP聚類分析(圖4-C)。
圖4 主成分分析Fig.4 Principal component analysis
Seurat 中的表達矩陣使用 PCA分析之后,再采用 UMAP 進行降維,得到了14個細胞cluster(圖5-A),根據(jù)已知的標記基因?qū)γ總€ cluster 的細胞類型進行鑒定(圖5-B),結(jié)果顯示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表達內(nèi)根鞘細胞(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
在對單細胞數(shù)據(jù)進行不同類型聚類之后,利用點狀圖進一步對每個cluster特異性表達的基因進行比較分析。結(jié)果如圖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
對聚類之后主要細胞類群進行富集分析,結(jié)果發(fā)現(xiàn)表皮細胞聚類得到0、3、8、9共4個cluster,特異表達基因分別為460、408、655、690個。GO富集結(jié)果顯示,4個亞群主要富集在GO:0008544表皮發(fā)育、GO:0030855上皮細胞分化、GO:0009611損傷應(yīng)答、GO:0000904調(diào)控細胞形態(tài)發(fā)生等生物過程,而circos圖也顯示,這4個亞群之間存在密集的共同信號通路(圖7-A);IFE-DC細胞聚類得到1、4、10共3個cluster,分別有582、1 061、341個特異表達基因,GO富集結(jié)果顯示:這3個亞群主要富集在GO:0008544表皮發(fā)育、GO:0097435超纖維組織發(fā)育、GO:0030155細胞黏附調(diào)節(jié)以及GO:0048729組織形態(tài)發(fā)生等生物過程中,同時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富集分析發(fā)現(xiàn),這2個亞群主要富集在GO:0008544表皮發(fā)育、GO:0034330細胞連接、GO:0000165 MAPK級聯(lián)激活以及GO:0048729組織形態(tài)發(fā)生等生物過程(圖8-A)。對黑色素細胞進行分析發(fā)現(xiàn)了752個差異基因,GO富集結(jié)果顯示黑色素細胞主要富集在GO:0030155細胞黏附的調(diào)控、GO:0030855上皮細胞分化、GO:0030335細胞遷移的正向調(diào)控以及GO:0030029肌動蛋白等生物過程(圖8-B)。
圖8 HS與黑色素細胞特征基因富集分析Fig.8 Enrichment analysis of HS and melanocytes characteristic genes
對INFU進行分析發(fā)現(xiàn),其BP共有486個條目,主要與表皮發(fā)育、上皮細胞增殖、肌動蛋白絲以及細胞-底物黏附等有關(guān);CC分析中共有82個條目,主要包括細胞-底物連接、黏著斑通路、細胞間連接以及細胞皮層;MF分析中共有27個條目,主要富集在鈣粘蛋白結(jié)合、肌動蛋白結(jié)合以及肌動蛋白絲結(jié)合等條目中(圖9-A)。對IRS進行分析發(fā)現(xiàn),其BP共有733個條目,主要與核糖核蛋白復合體的生物發(fā)生、ncRNA代謝過程和表皮發(fā)育等相關(guān);CC分析中共有104個條目,特征基因主要富集在細胞-底物連接、黏著斑通路等GO通路;MF分析中共有63個條目,主要富集于轉(zhuǎn)錄核心調(diào)節(jié)因子活性和鈣黏蛋白結(jié)合等條目中(圖9-B)。
圖9 INFU細胞與IRS細胞特征基因的GO富集分析Fig.9 GO enrichment analysis of characteristic genes in INFU and IRS cells
通過KEGG分析發(fā)現(xiàn),Epidermal特征基因主要富集于ko03010:核糖體、ko04520:黏附連接等通路;IFE-DC特征基因主要富集于ko03050:蛋白酶體等通路;HS特征基因主要富集于ko04530:緊密連接、ko04360:軸子引導等通路;IRS特征基因主要富集于ko03030:DNA復制、ko03040:剪接體等通路;INFU特征基因主要富集于ko03050:蛋白酶體、ko03013:RNA轉(zhuǎn)運、ko04110:細胞周期、ko03030:DNA復制、ko03040:剪接體等通路;Melanocyte特征基因主要富集于ko03010:核糖體、ko04360:軸子引導信號通路等(圖10-A)。運用metascape對特征基因進行互作分析,獲得不同類型細胞基因互作網(wǎng)絡(luò)圖,共同靶基因有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
利用熒光免疫組化對天祝白牦牛退行期皮膚組織進行了染色分析如圖11所示。結(jié)果顯示,SPARC陽性的細胞主要集中在HS中,其次在表皮細胞系以及內(nèi)根鞘中也存在一定的表達;KRT10主要在表皮中高表達,而在真皮細胞系中呈現(xiàn)出弱表達;此外,從皮膚組織中的表達情況中可以看出DLX3在毛干周圍的內(nèi)根鞘和外根鞘中高表達。對關(guān)鍵標記基因進行免疫組化分析得到的實驗結(jié)果與上述分析結(jié)果基本一致。
圖11 關(guān)鍵標記基因熒光免疫組化分析Fig.11 Fluorescence immunohistochemical analysis of key marker genes
毛囊作為皮膚組織中的微型器官,與微環(huán)境的其他組織、細胞共同的受到諸多信號通路調(diào)控[9]。毛囊的發(fā)育與成熟、毛囊的周期性循環(huán)均依賴于真皮與表皮的相互作用[10]。真皮細胞發(fā)出的啟動信號使相鄰部位表皮增厚,表皮細胞發(fā)出的信號促使真皮細胞在基板下聚集,誘發(fā)基板下的真皮細胞發(fā)育為毛乳頭,并由毛乳頭誘使毛囊內(nèi)表皮細胞繼續(xù)增殖分化,形成毛干并突出皮膚生長[11]。在退行期,毛乳頭分泌的調(diào)控因子減少,毛母質(zhì)細胞萎縮,毛囊再生部分退化,毛球部與外根鞘內(nèi)的部分細胞凋亡,到退行期后期基質(zhì)消失[12-13]。因此利用scRNA-seq對細胞進行生物信息學分析,可為下一步試驗提供指導。
單細胞轉(zhuǎn)錄組測序一次可以檢測產(chǎn)生上萬個細胞的轉(zhuǎn)錄信息,因此在對細胞類型進行注釋之前需要進行降維。主成分分析是一種利用投影矩陣將高維數(shù)據(jù)投影到低維空間的線性降維方法,可將多個變量轉(zhuǎn)換為少數(shù)幾個反映絕大部分信息的綜合變量即主成分[14]。tSNE是一種非線性降維方法,其巧妙的解決了集群表示過度擁擠問題。但有研究發(fā)現(xiàn)該算法比較消耗計算機資源,易丟失集群間關(guān)系而無法有效地表示數(shù)量龐大的數(shù)據(jù)集[15-16]。UMAP是建立在黎曼幾何和代數(shù)拓撲理論框架上的一種可視化降維方法,并對嵌入維數(shù)沒有計算限制。Satpathy等[17]發(fā)現(xiàn)UMAP在單細胞轉(zhuǎn)錄組測序中的可視化過程中效果理想。
Ryan等[18]發(fā)現(xiàn)Rho家族GTPases能夠通過調(diào)節(jié)鈣粘蛋白、整合素和黏附結(jié)構(gòu)來對表皮功能進行調(diào)控。Fuchs[19]發(fā)現(xiàn)表皮必須通過增殖、分化和凋亡不斷補充和維持自身穩(wěn)態(tài)。正常表皮功能的主要決定因素是細胞黏附,除了維持表皮結(jié)構(gòu)外,細胞-基質(zhì)和細胞-細胞黏附在調(diào)節(jié)信號轉(zhuǎn)導和表皮層的細胞功能方面也起著關(guān)鍵作用[20]。葛偉[21]對陜北白絨山羊毛囊發(fā)育過程中的表皮細胞進行GO富集分析發(fā)現(xiàn),表皮細胞聚類得到的5個cluster共同富集了表皮發(fā)育、上皮細胞分化等信號通路。這與本研究中牦牛表皮細胞富集結(jié)果相似,說明牦牛與其他動物毛囊發(fā)育存在相似性。葉娜[22]對聚類得到的IFE-DC細胞進行GO富集分析發(fā)現(xiàn),IFE-DC細胞主要富集在表皮發(fā)育、超纖維組織發(fā)育以及創(chuàng)傷反應(yīng)調(diào)控等通路。與本研究IFE-DC細胞富集分析結(jié)果一致。
MAPK級聯(lián)是毛干細胞主要富集通路之一。而MAPK信號通路能誘導毛囊細胞的增殖分化,促進毛囊的周期性發(fā)育。有研究顯示,Tβ4能夠影響p-catenin、VEGF以及MMP2的表達,使MAPK通路中的P38被激活,進而影響絨毛生長、毛囊分布和毛干數(shù)量[23]。在本研究中,黑色素細胞特征基因主要富集在軸突引導、細胞黏附的調(diào)控、細胞遷移的正向調(diào)控以及肌動蛋白等通路中。有研究發(fā)現(xiàn)黑色素在黑色素細胞內(nèi)合成后,通過黑色素細胞的樹突結(jié)構(gòu)以黑色素顆粒的形式轉(zhuǎn)運至周邊表皮,黑色素轉(zhuǎn)運過程需要微管驅(qū)動蛋白以及動力蛋白等的參與[24]。高澤成[25]利用HE染色發(fā)現(xiàn)天祝白牦牛毛囊毛球部黑色素細胞無法生成黑色素顆粒,導致天祝白牦牛毛纖維中不存在黑色素顆粒而呈現(xiàn)白色。在毛囊形態(tài)發(fā)生期間,黑色素細胞前體遷移到發(fā)育的毛囊中,并產(chǎn)生去分化的黑色素細胞,該黑色素細胞能夠主動產(chǎn)生色素并將色素運輸?shù)浇琴|(zhì)細胞中[26]。但黑色素細胞在軸突引導方面的生物學功能仍需進行進一步驗證。
對天祝白牦牛毛囊退行期進行分析,鑒定出IFE-DC細胞、表皮細胞以及INFU細胞等細胞類型,特征基因主要參與表皮發(fā)育、上皮細胞分化以及細胞形態(tài)學發(fā)生等生物學過程,KEGG分析顯示特征基因顯著富集于黏附連接、細胞周期以及RNA轉(zhuǎn)運等通路中,同時發(fā)現(xiàn)退行期主要細胞類型擁有GJA1、KRT1、KRT80、FGFR2等28個共同基因。