天祝白牦牛退行期毛囊细胞亚群鉴定以及特征基因生物信息学分析

郑青波 叶娜 张哓兰 包鹏甲 王福彬 任稳稳 廖月姣阎萍 潘和平

（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 材料

皮肤样品来自甘肃省武威市天祝县天祝白牦牛育种场。

主要试剂：中性树脂、二甲苯、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单细胞测序 由北京贝瑞和康生物技术有限公司完成单细胞文库的构建和测序。基于10×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个共同基因。