多发性硬化症（multiple sclerosis，MS）是一种中枢神经系统的自身免疫性疾病，其特征是多发病灶处炎症性脱髓鞘^［1］。免疫细胞穿过血脑屏障，引起炎症、脱髓鞘、胶质增生和轴突变性，造成感觉和视觉障碍、运动障碍、疲劳、疼痛和认知障碍^［1-3］。MS最常见的形式是复发缓解型，其特征是反复地发作和缓解，后遗症轻微。但最终80%的患者在诊断后的10~20年会继续发展为继发进行型MS^［2］。目前我国已上市的针对MS的药物共有7种，以疾病修正治疗（disease-modifying therapy，DMT）控制疾病活动^［4］。但DMT仅限复发缓解型，对进行性的亚型效果仍不理想^［5］。因此，对于MS治疗开展更深入的研究是必要的。

小胶质细胞是中枢神经系统的常驻免疫细胞，在免疫防御和维持大脑稳态方面起着至关重要的作用。作为先天免疫系统的一部分，小胶质细胞监测着中枢神经系统的潜在内部和外部异常^［6-9］。此外，小胶质细胞还具有控制神经元增殖、突触形成和消除、碎片清除、调节髓鞘生长等重要功能^［10-12］。

在MS的进展中，小胶质细胞作为大脑内的常驻免疫细胞，与MS的所有阶段有关^［13］。小胶质细胞在MS病变中的作用是复杂且具有争议的。有益特性包括其可能分泌神经营养因子^［14-16］，促进少突胶质细胞前体细胞成熟和髓鞘生成^［16］，清除髓鞘碎片^［17-19］等；有害特性则包括可能存在的轴突损伤作用和促脱髓鞘作用等^［20-23］。目前，人们认为小胶质细胞的起始反应是有益的，但慢性激活下的小胶质细胞在MS进展中起到有害的作用^［24］。MS的一大特征是在病灶局部小胶质细胞和星形胶质细胞的激活，共同促进炎症、脱髓鞘和神经退行性病变^［24-27］。小胶质细胞在MS中的作用是至关重要的。

单细胞转录组测序技术能使人们通过转录组的表达模式，在单个细胞的分辨率下，更好地解析小胶质细胞的异质性^［28-29］。ABSINTA等^［28］在snRNA-seq数据中描述了MS中发炎的小胶质细胞亚群，解剖来源为白质，并发现C1q是该亚群激活的关键。CLARK等^［29］使用狂犬病条形码和单细胞转录组测序技术发现小胶质细胞与星形胶质细胞之间存在相互作用，轴突引导分子被鉴定为小胶质细胞和星形胶质细胞相互作用的潜在介质，会加重自身免疫性脑脊髓炎小鼠模型的疾病程度。其他MS的单细胞转录组研究^［30-32］也证明了小胶质细胞在MS进展中发挥着重要的作用。但是以上研究均未从小胶质细胞的异质性出发，对来自不同解剖区域的小胶质细胞进行转录调控网络等系统性分析。

基于此，本研究从基因表达数据库（Gene Expression Omnibus，GEO）获取人脑单细胞核转录组测序数据，从灰质和白质的不同解剖部位出发，探究在MS过程中小胶质细胞的转录组转变和关键转录调控网络。

本研究使用完整的单细胞核转录组测序（single nucleus RNA sequencing，snRNA-seq）和RNA-seq数据集，下载于NCBI基因表达综合数据库（https：//www.ncbi.nlm.nih），登记编号为GSE180759^［28］、GSE227781^［32］和GSE179427^［33］。从数据集GSE227781中收集到6名MS患者和6名健康个体的snRNA-seq数据，按照解剖部位分为灰质、白质和患者特有的脱髓鞘病灶区。从数据集GSE180759中收集到5名MS患者和3名健康对照个体的snRNA-seq数据，解剖部位为白质。从数据集GSE179427中收集到41名MS患者和26名健康对照个体的RNA-seq数据，解剖部位为白质。

在R（4.2版本）中使用Seurat（4.3版，https：//satijalab.org/seurat/）^［34］读入GSE180759和GSE227781数据集的基因表达矩阵和元数据信息。使用CreateSeuratObject函数将矩阵转化为Seurat对象。质量控制方面，少于300个基因、UMI数大于50 000和线粒体基因大于10%的细胞核被去除。使用SCTransform函数^［35］进行标准化和归一化，并使用vars.to.regress参数去除线粒体基因对后续降维聚类的影响。使用3 000个高变基因进行主成分分析（PCA）。使用Harmony^［36］减轻批次效应和技术偏差。前30个主成分（principal component，PC）被选择，并使用RunUMAP函数进行均匀流形近似和投影（uniform manifold approximation and projection，UMAP）可视化。

通过Seurat软件包中的FindAllMarker函数计算每个亚群和其他细胞之间的标志基因。阈值参数设置为logfc.threshold=0.7，min.pct=0.5。只保留了校正后的P值小于0.05的基因。根据CellMarker2.0数据库（http：//bio-bigdata.hrbmu.edu.cn/CellMarker/）中现有的注释及先验知识，对细胞进行注释^［37］。

提取小胶质细胞后，按照解剖部位信息来源于灰质或白质分别进行亚群分析。亚群分析步骤同“1.2”，选取前15个PC进行UMAP可视化。

样本白质或灰质亚群细胞数与样本白质或灰质细胞总数的比值为细胞比例。MS组和对照组间的细胞比例变化分析使用Wilcoxon符号秩检验统计。

通过Seurat软件包中的FindAllMarker函数分别计算灰质和白质内MS组和对照组间的差异表达基因（differentially expressed genes，DEGs），判断上调和下调的阈值为log₂FC（log₂ Fold Change）的绝对值大于0.5，且校正后的P值小于0.05。使用clusterProfiler包的gseGO函数和gseKEGG函数，进行基因集富集分析，q-value<0.05^［38-39］。

白质小胶质细胞亚群间的DEGs及基因本体论（Gene Ontology，GO）分析使用SCP包完成，DEGs阈值选择为log₂FC大于1，且校正后的P值小于0.05。单个亚群超过50个DEGs则取log₂FC值最大的50个基因进行后续分析和绘图。

使用Monocle3软件进行拟时序分析^［40-43］。导入在Seurat软件中计算得到的降维坐标进行分析。使用cluster_cells函数将全部细胞分区，选择细胞数最多的主要分区进行拟时序分析。将细胞按照拟时序的时间顺序进行排序，使用pheatmap包进行热图绘制，并将基因分为4个模组。将每个模组的基因进行GO分析。

使用单细胞调控网络推理和聚类（single cell regulatory network inference and clustering，SCENIC）方法进行转录因子（transcription factor，TF）和靶基因之间的转录调控网络分析^［44］。SCENIC工作流程使用默认参数，使用所有样本的原始计数矩阵输入，调用python的grnboost2函数进行调控网络的构建，使用calcRSS函数计算每个调节子在不同亚群中的调节子专一性得分（regulon specificity score，RSS）。

使用GSVA包^［45］进行调节子基因集打分后使用limma包^［46］进行差异通路分析。

对基因表达量和线粒体表达量进行筛选后，共149 062个细胞核被保留。

变异量最大的3 000个基因被选取作为高变基因，用于PCA分析和UMAP可视化。将全部细胞分为10个亚群（clusters），见图1A。Cluster 1，小胶质细胞，共有12 238个细胞，占比8.21%，以DOCK8、CSF1R、P2RY12、CD74为特征；Cluster 2，少突胶质细胞，共有86 957个细胞，占比58.33%，以MBP、MOBP、PLP1为特征；Cluster 3，少突胶质细胞前体细胞（oligodendrocyte precursor cells，OPCs），共有6 245个细胞，占比4.19%，以VCAN、OLIG1、OLIG2为特征；Cluster 4，星形胶质细胞，共有16 836个细胞，占比11.29%，以AQP4、CLU、GFAP为特征；Cluster 5，室管膜细胞，共有1 198个细胞，占比0.80%，以DNAH9、CFAP44为特征；Cluster 6，兴奋性神经元，共有12 301个细胞，占比8.25%，以SYT1、ROBO2、CAMK2A、SLC17A7为特征；Cluster 7，抑制性神经元，共有7 372个细胞，占比4.95%，以SYT1、ROBO2、GAD1、GAD2为特征；Clusters 8，T淋巴细胞，共有1 008个细胞，占比0.68%，以IL7R、PTPRC为特征；Cluster 9，内皮细胞，共有2 049个细胞，占比1.37%，以FLT1、CLDN5为特征；Cluster 10，周细胞/成纤维细胞，共有2 858个细胞，占比1.92%，以PDGFRB、DCN为特征。使用UMAP图在全部细胞中展示2种小胶质细胞标志基因表达量的空间分布，以确定小胶质细胞亚群（图1B）。不同细胞的标志基因用点图展示（图1C）。

通过DOCK8、CSF1R、P2RY12、CD74等小胶质细胞标志基因，共鉴定出12 238个小胶质细胞。考虑到灰质和白质在细胞构成、代谢状态以及神经信号传递中存在本质区别，且MS主要发生在白质而非灰质中，在后续的分析中根据小胶质细胞的解剖部位来源于灰质或白质分别进行分析。

根据小胶质细胞经典标志基因和解剖部位来源，共鉴定出2 340个灰质胶质细胞。在来自灰质的小胶质细胞中进行降维聚类分析，共得到7个亚群，命名为GM_MG1‒7（Gray Matter_Microglia 1‒7）。GM_MG1亚群中包含了最多的细胞，共1 198个，占所有灰质小胶质细胞的51.20%；GM_MG7亚群中所包含的细胞数最少，占所有小胶质细胞的1.62%（图2A）。小胶质细胞的关键标志物CD74、CSF1R、DOCK8和PTPRC在所有亚群中均有表达，证明所取到的细胞中没有掺杂其他细胞（图2B）。UMAP图显示来自不同样本的细胞分布均匀，且未观察到某一样本的特异性亚群，说明样本间分布无批次效应（图2C）。MS组和对照组细胞在各个亚群的数量上未表现出明显的偏差，且未观察到对照组或MS组的特异性亚群（图2D）。

对灰质中MS组和对照组间的差异作进一步的分析。各个灰质小胶质细胞亚群在MS组与对照组间的数量比例无显著差异（图3A）。DEGs分析结果（图3B）显示：MS组对比对照组有13个基因上调（包括GRID2、DSCAM、SYNDIG1等），有20个基因下调（包括SLC11A1、DPYD、SLC16A10等）。基于表达量的富集分析表明：灰质小胶质细胞在疾病过程中内吞作用、离子稳态、脂质定位功能等GO生物学过程下调（图3C、E）；长时程抑制通路、MAPK信号通路等KEGG通路上调；核糖体通路、PPAR信号通路、溶酶体通路等KEGG通路下调（图3D、F）。该结果提示灰质小胶质细胞受到了疾病影响。

根据小胶质细胞经典标志基因和解剖部位来源，共鉴定出9 898个白质小胶质细胞。在来自白质的小胶质细胞中进行降维聚类分析，共得到6个亚群，命名为WM_MG1‒6（White Matter_Microglia 1‒6）。WM_MG1亚群中包含细胞数最多，共6 122个，占全部白质小胶质细胞的61.85%；WM_MG6亚群中包含的细胞数最少，共206个，占全部白质小胶质细胞的2.08%（图4A）。小胶质细胞的关键标志物CD74、DOCK8、CSF1R和PTPRC在各个亚群中的表达情况（图4B）表明，所获得的白质小胶质细胞中不掺杂其他杂质细胞。UMAP图显示，来自不同样本的细胞分布均匀，且未观察到某一样本的特异性亚群，说明样本间分布无批次效应（图4C）。在WM_MG2亚群和WM_MG5亚群中，观察到MS组细胞多于对照组细胞，在其他亚群中未观察到2组细胞数量的明显偏差，且未观察到对照组或MS组的特异性亚群（图4D）。

对白质中MS组和对照组间的差异作进一步的分析。WM_MG2和WM_MG6亚群的细胞比例在MS组的上调具有统计学意义，说明这2个亚群为MS相关亚群。WM_MG1、WM_MG3、WM_MG4和WM_MG5表现出上调趋势，但不具有统计学意义（图5A）。DEGs分析结果（图5B）表明，MS组对比对照组有46个基因上调（包括HSPB1、HSP90AA1、SERPINH1、FTH1等），有25个基因下调（包括KCNIP1、DLEU1、RASGEF1C等）。基于表达量的富集分析表明，白质小胶质细胞在疾病过程中蛋白质折叠、对未折叠蛋白的反应、细胞质翻译等生物学过程上调，动作电位期间膜去极化过程下调（图5C、E），核糖体通路、内质网中蛋白质加工通路等KEGG通路上调，Rap1信号通路下调（图5D、F）。对比灰质和白质小胶质细胞在MS过程中的DEGs和通路富集结果，提示白质小胶质细胞在MS过程中所受到的疾病影响大于灰质小胶质细胞。

对白质小胶质细胞内的异质性进行分析（图6）。WM_MG1亚群中富集了包括突触囊泡聚集、上皮细胞极性的建立、突触组装的正向调节等生物学过程，且WM_MG1中P2RY12的表达上调，因此推测WM_MG1为稳态的小胶质细胞亚群。WM_MG2亚群中富集了蛋白质折叠、热响应、温度刺激响应等生物学过程。值得注意的是，WM_MG2亚群中富集到的DEGs和生物学过程，与图5B中MS组相比对照组上调的生物学过程有着高度的重合，如HSPH1、HSPD1、HSP90AA1、P4HA1等。上述基因在WM_MG2中的log2FC值均大于2，因此，推测WM_MG2为MS相关的疾病亚群。WM_MG3亚群中富集了包括烟酰胺核苷回收途径的NAD生物合成、凝血纤维蛋白凝块形成、一氧化氮合酶生物合成过程的正向调控等生物学过程，推测WM_MG3为活化的小胶质细胞，具有神经保护功能。WM_MG4亚群中富集了包括细胞连接组装、蛋白质在细胞表面的定位、细胞通过细胞膜黏附分子进行细胞黏附等生物学过程，结合WM_MG4上调表达的基因，如F13A1、MRC1、CD163、STAB1、COLEC12等，推测WM_MG4可能在大脑的损伤修复、抗炎、吞噬髓鞘等方面起重要作用。WM_MG5亚群富集了细胞质翻译、核糖体生物合成、核糖体组装等生物学过程；在基因层面，该亚群富集到了TREM2、APOE、CD68等基因（log₂FC值均大于0.4，校正后的P值均小于0.05）。与核糖体相关的基因大量上调，说明WM_MG5表现出高度的代谢活跃性。WM_MG6亚群富集了细胞连接组装、突触组织、神经元迁移等生物学过程，说明WM_MG6亚群可能在促进神经系统的正常发育、神经连接的维护中发挥关键作用。DEGs分析和富集分析结果提示，白质小胶质细胞各个亚群间存在着较为清晰的DEGs和功能差异。

使用拟时序分析将细胞按照“假设时间顺序”排列，揭示MS进展中小胶质细胞内部基因表达的动态变化。首先，Monocle3将细胞分为3类（图7A），每类细胞代表着一个单独的变化轨迹。选择其中涉及细胞数最多的第一类细胞（图中以粉红色显示），根据细胞与细胞间的转录组差异，将细胞以假设时间变化进行排列（图7B、C）。在第一类细胞中共得出2个演化分支，WM_MG1被定位在演化轨迹的起始点，表明这可能是最原始的小胶质细胞状态。第一个分支由WM_MG1起始，经过WM_MG3的部分细胞，假设时间顺序结尾于WM_MG2亚群，这一分支由稳态的小胶质细胞向活化的小胶质细胞演化。WM_MG6是从WM_MG1分出的第二个分支，提示了维护神经元功能的演化。WM_MG3的部分细胞和WM_MG4被识别为处于其他单独的变化轨迹中。

将细胞在横轴上按照拟时序分析所得的假设时间进行排序，左侧为假设时间的起始，右侧为假设时间的终点。纵轴基因为拟时序轨迹高度相关基因。共分为4个基因表达的模块，热图展示4个不同模块中基因的表达量随假设时间变化的趋势（图7D）。模块1为在假设时间起始的细胞中上调的基因，随着假设时间的推移下调，包含MEGF10、ENOX1、PPP1R16B等基因。模块1中的基因富集了中枢神经系统髓鞘形成、突触前膜组织、突触组织等生物学过程，这些基因随着假设时间顺序下调，提示了小胶质细胞维持神经网络功能能力的下降。模块1中的较多基因和WM_MG1亚群中的DEGs吻合，如KCNIP1、SYNDIG1、ELMO1等基因，WM_MG1推测为稳态的小胶质细胞亚群。因此，模块1中的基因随着假设时间的推移下调，也提示了稳态小胶质细胞的这些功能在MS进程中下降。模块2为在假设时间终点的细胞中上调的基因，随着假设时间的推移上调，包含DNAJB1、FKBP4、P4HA1等基因。模块2中的基因富集了蛋白质折叠、伴侣蛋白介导的蛋白质折叠、热响应等生物学过程，与WM_MG2中上调的生物学过程是一致的；且模块2中较多基因与WM_MG2亚群中的DEGs吻合，如DNAJB1、FKBP4、SERPINH1，再次验证了小胶质细胞激活过程中上调的基因和生物学过程。模块3和模块4为在假设时间居中的细胞中上调的基因，模块3和模块4中的基因未富集到具有统计学意义的生物学过程。

用SCENIC绘制控制疾病的基因调控网络，鉴定MS进程中调节DEGs的潜在TF。对单细胞RNA测序数据进行预处理，将原始表达矩阵输入SCENIC中，仅保留在RcisTarget中可用的基因共8 318个。使用AUCell方法计算每个调节子在细胞中的活性，从而计算MS组中TF及其靶基因构成的调节子的RSS得分并进行排序，结果显示FLI1、MITF、FOXP1等调节子得分最高（图8A）。为了获得调节子与每种细胞之间的特定对应关系，计算调节子在不同小胶质细胞亚群中的RSS得分，并使用热图展示（图8B）。结果表明，在灰质小胶质细胞中，未发现明显上调或下调的TF。FLI1、MITF、FOXP1等调节子的上调主要在白质的WM_MG1亚群。同时，在WM_MG2的MS组中，还有SREBF2、NFKB1、REL等调节子上调；WM_MG4的MS组中，NKX6-2调节子上调。这些TF可能是小胶质细胞向疾病状态转变的关键。其他亚群中未发现明显上调的调节子。

为了验证SCENIC分析的结果，采用RNA测序数据结合GSVA方法计算所鉴定的调节子基因集的表达情况，并进一步对MS组与对照组之间的基因集差异表达进行定量分析（图8C）。结果显示，在SCENIC分析中识别出的得分最高的5个调节子，即FLI1、MITF、FOXP1、HIF1A和ETV6，在MS组与对照组之间均表现出了具统计学意义的上调。

MS是一种慢性中枢神经系统脱髓鞘疾病，其发病机制涉及炎症和神经退行性过程。MS确切病因尚不清楚，通常被认为是一种自身免疫性疾病，由自身反应性淋巴细胞对中枢神经系统自身抗原的反应引发^［26-27］。

本研究使用snRNA-seq数据对MS进程中小胶质细胞在灰质与白质间的差异性表达及其在疾病进展中的作用进行分析，并鉴定了疾病相关的关键转录调控网络。在灰质中聚类得到了7个小胶质细胞亚群，在白质中聚类得到了6个小胶质细胞亚群。从细胞亚群的数量和亚群间距来看，小胶质细胞在灰质中具有更高的异质性。从数量比例来看，小胶质细胞在灰质中并未发现有显著数量变化的亚群，且这些亚群也并无明显的上调或下调趋势。相反地，白质中发现了WM_MG2、WM_MG5两个疾病相关亚群，其在MS组比例上调。本研究观察到白质所受到的疾病影响大于灰质，这与MS往往发生在白质中这一现象是一致的。但灰质依然受到了MS的影响，在其他相关研究^{［30，47-48］}中也有描述。

白质小胶质细胞受到MS的影响更加严重。白质小胶质细胞内上调的基因包括HSPH1、HSPD1、HSP90AA1等，为热休克蛋白家族的成员，参与多种与蛋白质折叠相关的生物过程。研究^［49］表明，热休克蛋白家族成员的表达上调可能与多种神经退行性疾病均有关联。

使用SCENIC进行的转录调控网络分析结果表明，在疾病WM_MG1中，上调的调节子包括HIF1A、CHD2、REL、FLI1、MITF、FOXP1等。HIF1A的上调可能增强了其对炎症的反应能力。REL属于REL/NF-κB家族，调节参与细胞凋亡、炎症、免疫反应和致癌过程的基因。在WM_MG2中，SREBF2、NFKB1调节子上调。NF-κB是典型的炎症因子，在机体内是免疫应答的重要调节中心，分析结果说明疾病状态下小胶质细胞的促炎作用被激活^［50］。

本研究分析了MS中小胶质细胞的转录组病理性转变，对比了灰质和白质在疾病过程中受到的不同影响，以及在疾病中可能发挥的作用；使用拟时序分析观察了疾病进程中的基因动态变化，鉴定了疾病中的关键转录调控网络。本研究结果可为疾病诊断和治疗提供新的方向。