12分Nature子刊教你分析细胞异质性还共享了数据！

Heterogeneity of human bone marrow and blood natural killer cells defined by single-cell transcriptome单细胞转录组定义的人类骨髓和血液自然杀伤细胞的异质性

一、 研究背景

自然杀伤（NK）细胞是1型先天淋巴样细胞(ILCs)，能通过发挥细胞毒性和产生细胞因子来对抗病毒感染，但其发育过程，转录组状态和异质性未有明确定义。根据CD56表达，可将NK细胞分为CD56 bright和CD56 dim两个亚群，而CD56 dim细胞基于CD57被进一步划分。CD56 bright产生细胞因子的能力更强，主要起免疫调节作用，而CD56dim主要为细胞毒作用。目前认为NK细胞是按照CD56bright → CD56dim CD57− → CD56dim CD57+这种线性发展模式发育的。但是有部分研究表明CD56bright和CD56dim来自不同的前体，此外具有GATA2突变的个体仅拥有CD56 dim细胞，这些研究均提示两者可能为独立的ILC1群体。

基于流式细胞术（CyTOF）的免疫谱分析已经揭示了基于多个表面受体表达的不同表型的NK细胞，强调了进一步定义NK群体异质性的重要性。而单细胞RNA测序（scRNA-seq）技术能够基于每个细胞的转录组分析群体的异质性，有助于进一步揭示NK细胞的异质性。

二、 分析流程

三、 结果解读

1.scRNA-Seq分析揭示了不同的人类NK亚型

首先，作者一共对来自六名健康供体的骨髓(BM)和两名健康供体的外周血的NK细胞进行了scRNA-seq。测序后，将原始数据demultiplexed来获得mRNA count，与hg19参考基因组比对，并使用10X Genomics Cell Ranger对UMI计数进行定量。然后，使用Seurat包对过滤后的barcode矩阵进行数据分析。在初次质量控制(QC)步骤中，作者过滤掉表达<200个基因或> 2500个基因的细胞，去除了线粒体转录> 5％的细胞。基于每个细胞中UMI数量和线粒体转录对基因表达值进行缩放和处理批次效应后，选择了一定的主成分PC，然后用t-SNE进行聚类分析。

由于NK祖细胞和部分未成熟的NK细胞没有CD56（NCAM1）表达，因此作者根据CD7+而不是CD56+对NK细胞进行分选。经过聚类，BM的CD7+细胞被分为9个亚群，接下来作者对这9个亚群进行检测，看其是否为NK细胞。作者确定了四个替代CD56的标记来定义NK细胞，包括有明确定义的NK细胞标记——CD94（KLRD1）和NKp80（KLRF1），以及其他研究中得出的差异基因NKG7和GNLY。然后作者发现＃8和＃9类并不是NK细胞（＃8高表达B细胞或树突状细胞标记（IGJ / MZB1 / LILRA4）；＃9中高表达T细胞特异性标记（CD3D/E/G）。（补充图2）

去除#8/9后，分别得到7个BM和5个外周血的亚群（图1a. 图2a）。根据亚群的差异基因DEG（log fc ≥0.25, min.pct = 0.1,adjusting p value <0.05，图1b. 图2b），在BM中作者将7个亚群定义为“CD56bright NK”， “Transitional NK”， “Active NK”，“Adaptive NK”，“Mature NK”，“Terminal NK”以及“Inflamed NK“。外周血的聚类也得到类似的分组，但缺乏“Adaptive NK”和“Inflamed NK“。接下来，作者将分别对每个亚群进行研究。

补充图2. BM的NK细胞初次QC后的聚类结果

图1.来自BM的人类NK细胞的无偏聚类

图2.来自外周血的人类NK细胞的无偏聚类

2.CD56 bright和CD56 dimNK细胞之间存在过渡的细胞亚群

图3.a显示的是CD56 bright细胞一些特征基因的高表达/低表达。作者在接下来的分析中，不断地结合单细胞转录组数据和CyTOF结果，增加了单细胞分析结果的可靠性。比如，BM和血液中的“CD56 bright NK”比例与流式细胞仪确定的细胞比例一致，图1b中显示的CD44高表达也与流式细胞仪的分析结果一致（图3b）。趋化因子C家族的XCL1和XCL2表达也较高。来自NK细胞的XCL1已被证实可以在肿瘤微环境中募集1型树突状细胞，对于抗肿瘤免疫至关重要。作者在蛋白质水平上通过流式胞内染色，也证实了CD56 bright细胞中的XCL1高于CD56 dim细胞（图3c）。

与此同时，作者发现了一个与CD56 bright亚群转录特征相似的亚群（图3d），该亚群下调了CD56 bright标志性基因，上调了与CD56 dim细胞相关的基因（如CD16A即FCGR3A），提示该亚群为bright和dim的过渡亚群。所以作者进一步计算了该亚群基于“CD56 bright NK”和“ Mature NK”的上调或下调的基因模块分数。“Transitional NK”均具有中等水平的分数（图3e）。接下来，作者评估了来自于其他研究的bulk RNA-seq数据中，CD56bright和CD56dimCD57+ NK细胞的DEG在本篇文章定义的所有NK亚群的表达水平，发现“Transitional NK”在BM和血液样本中均具有中等表达（图3f）。

图3.鉴定 CD56bright和CD56dim细胞的潜在过渡亚群

3.稳态下活跃的NK细胞亚群的转录组

在BM中发现的另一个亚群是“Inflamed NK”，该亚群前10个上调基因包括IFN诱导基因，TNF和趋化因子基因（图1b），提示干扰素可能会刺激NK细胞。GSEA结果也显示IFN-γ信号，Toll样受体信号和炎症反应的正调节在Inflamed NK亚群中富集（图4a）。此外，作者发现了另一个处于激活状态的亚群并将其标记为“Active NK”，“Inflamed NK”和“Active NK”的DEG较为相似（图4c）。在Active NK的DEG中，发现了一些立早基因（IEGs）（图4d），IEG在受到一系列外界刺激后能迅速并且短暂激活，不需要任何新蛋白的合成。GSEA结果显示，“Active NK”具有较高的转录活性，且多种信号通路被激活（包括KRAS-MAPK，TRAF6-NF-κB）以及各种刺激或受体（包括TNF，IL-2，TLR和GPCR）被激活，这些结果共同揭示了该亚群NK细胞的活跃状态。

接下来，作者尝试用流式细胞仪鉴定哪些细胞表面标志物可用于定义Active NK。在该亚群的DEG组中唯一的细胞表面蛋白是CD69和CXCR4（图4b,e），其中CXCR4 +NK cells在BM和血液中均可被检测到（图4f）。根据DEG（图1b,2b），转录组相似性（图3d）和模块得分分析（图3e,f），作者认为“Active NK”可能为CD56dim NK细胞。基于早期的研究，作者接下来用NKG2A，CD62L，CD57评估CD56dim CXCR4 +的BM和血液中的NK细胞成熟程度。随着NK细胞的成熟，NKG2A和CD62L的细胞表面表达会降低，而CD57通常被认为是NK细胞成熟的标志物。可以看到在未完全成熟的细胞中CXCR4+比例更高（CD57-，NKG2A+和CD62L+ 细胞中CXCR4+的百分比更高，图4g），这意味着NK细胞在功能成熟时会下调CXCR4。

因此作者得出结论，“Active NK”可能并不代表一个独特的发育阶段。作者通过评估bulk RNA-Seq中不同子集的（CD56bright, CD56dimCD57−, CD56dimCD57+）立早基因表达，发现了相似的表达水平，进一步证实了这一结论，“Active NK”可能只是来自不同发育阶段的细胞混合体，受到刺激后形成了独特的转录组谱。尽管不是一个独特的发育阶段，但该NK亚群的存在揭示了NK细胞的稳态激活，这对NK细胞的存活，增殖或分化有重要的影响。

图4. 具有独特转录组特征的Active NK细胞

4.适应性NK细胞的独特转录组谱

在BM中发现的另一个独特的亚群为“Adaptive NK”，其具有高表达的KLRC2（NKG2C）。KLRC2已被用作人类巨细胞病毒（HCMV）感染导致的适应性NK细胞的标志物。作者发现只有24岁的女性样本具有较高的KLRC2表达（图5a）。因此，接下来仅使用该样本进行了聚类分析。由于细胞数量少，只得到了四个亚群（图5b）。先前工作发现，HCMV供体的自适应NK细胞+表达更少FcεRγ（FCER1G），SYK，EAT-2（SH2D1B）和PLZF（ZBTB16），而在“Adaptive NK”中FCER1G和ZBTB16表达确实有显著降低（图5d）。

在24岁女性样本的“Adaptive NK”群中高表达的另一个基因是IL32，IL-32现已被认为是重要的促炎细胞因子。为了探讨IL-32的上调是否是与适应性NK细胞相关的独特功能有关，作者评估了所有7个BM NK亚群中IL-32的表达，发现“Adaptive NK”和“Inflamed NK”具有更多的IL-32表达。此外，IFN-γ产生的正向调节的基因也有富集，这表明了适应性NK细胞的功能特性（图5g）。

此外，与“Mature NK”相比，“Adaptive NK”具有更少的溶细胞分子表达。作者推测这些功能分子以蛋白质形式表达并存储在细胞中，因此不再检测到这些基因的主动转录。与该推测一致，GSEA发现与其他细胞相比，“Adaptive NK”内质网至高尔基体转运囊泡基因组有富集（图5g）。

图5. 24岁女性样本的适应性NK细胞

5.成熟的NK细胞表现出独特的转录谱

目前观点认为，CD57标记的成熟NK细胞具有最强的溶细胞能力，图6a可以看到“Mature NK”中与细胞毒性，细胞骨架重塑，细胞粘附和信号传导功能相关的分子表达较高。

在“Mature NK”簇中，作者发现属于膜联蛋白家族的基因表达增加（图6b），其中膜联蛋白A1对其他免疫细胞中具有免疫抑制作用，对于缓解糖皮质激素下游炎症至关重要，所以作者进一步对其进行研究。与转录数据一致，作者发现膜联蛋白A1随着NK细胞的成熟而增加（CD56bright, CD56dimCD57−, CD56dimCD57+）（图6c）。接下来，作者探讨了糖皮质激素-膜联蛋白A1轴对人NK细胞的效应功能的影响。糖皮质激素可抑制激活后NK细胞产生IFN-γ，但是膜联蛋白A1不影响PMA和离子霉素，IL-12和IL-18介导的IFN-γ产生。作者推测是由于与糖皮质激素受体（NR3C1）相比，NK细胞上膜联蛋白A1受体（FPR2）的表达较低。但是，人NK细胞胞浆中大量的膜联蛋白A1可能通过旁分泌来调节其他先天性和适应性免疫细胞的活性。为了验证这种可能性，作者发现通过刺激，NK细胞中膜联蛋白A1的细胞内水平降低，表明该蛋白可能存在分泌（图6d）。接下来，作者进一步利用酶联免疫吸附（ELISA）来测量刺激后培养上清液中膜联蛋白A1的蛋白水平，发现其含量增加（图6e）。这些数据表明，NK细胞是膜联蛋白A1的重要来源，经活化诱导的NK细胞的膜联蛋白A1分泌可能影响炎症反应。

图6. 成熟的NK细胞的转录特征

对于“ Terminal NK”，其功能分子的表达水平与“Mature NK”相似。转录因子ZEB2，以及CX3CR1和HAVCR2（TIM-3）的表达增加，表明了该亚群为成熟的NK群体（图7ab）。作者推测这两个亚群共同形成CD57+ NK细胞。

接下来，作者探索了“Mature NK”和“Terminal NK”簇之间的DEG，“ Terminal NK”小核仁RNA，C/D box 3（SNORD3）家族和各种组蛋白亚基表达较高（图7cd）。GSEA显示，核糖体基因集在“Mature NK”富集；而染色体维持，组蛋白甲基化和端粒维持基因集在“ Terminal NK”中富集（图7e）。此外，与“Mature NK”相比，“Terminal NK”的炎症反应降低，mTORC1 / STAT5信号下调以及代谢活性降低（图7f）。这些减弱的信号和代谢特征表明“Terminal NK”细胞可能处于静止状态，于是作者评估了亚群中与静止状态相关基因的表达水平，模块得分表明在BM或血液样品中“ Terminal NK”中这些基因的表达均无显著增加。尽管如此，这些数据仍表明CD57 + NK细胞的异质性，为以后的研究奠定了基础。

图7. 终末分化NK细胞的独特转录谱

6.BM和血液含有相似的发育性NK亚群

为了比较BM和血液的NK细胞亚群，作者结合了同时提供了BM和血液样本的F25y_2和F26_y的转录数据，像之前的分析一样，也发现了“CD56bright NK,” “Transitional NK,” “Active NK,” “Adaptive-like NK,” “Mature NK,”和“Terminal NK”亚群。BM和血液都具有相似的NK细胞异质性，只是组成有所不同。血液中的“Active NK”比BM多，“Adaptive-like NK”主要由来自BM。就成熟程度而言，血液具有更高的“Mature NK”，这与血液中具有更多的CD57 + NK细胞一致。

7.伪时间分析表明CD56 bright为CD56 dim NK细胞的前体

scRNA-seq分析能够根据细胞分化过程中的转录变化，模拟发育异质种群之间的差异。作者在这一步使用Monocle2包来研究NK细胞的发育过程。图8ab可以看到BM NK细胞和血液NK细胞均表现出相对线性的发展进程，几乎没有分支。“ CD56 bright NK”和“Mature/Terminal NK”在BM和血液的进展轨迹的两端占主导地位。其中“Transitional NK”均从“ CD56 bright NK”发育而来，并向“Mature NK”的方向延伸。该结果不仅证明了从CD56 bright 到CD56 dim NK细胞的过渡亚型的存在，而且也表明CD56 bright NK细胞的确是CD56 dim NK细胞的前体。

接下来，作者着重分析发展轨迹的分支。分支点A和B将细胞轨迹分为五个状态（图8c），饼图（图8d）计算了每个亚群每个状态内的百分比。State-1细胞主要由“ CD56 bright和Transitional NK”组成，而State-4细胞主要由“Mature/Terminal NK”组成。State-5细胞主要为“Adaptive NK”，其中CD3E，VIM和CCL5高表达，与图1b的基因表达对应。State-3细胞主要由“Active NK”组成。State-3和5形成的分支表明，这两个簇的转录组足迹不同于稳态成熟。

图8. 单细胞轨迹分析揭示了线性的NK细胞发育途径

8.GATA2 T354M供体中的NK细胞转录组发生改变

在GATA转录因子家族中，已知GATA2的单倍体缺失会导致人NK细胞缺乏，并特异性丧失CD56 bright亚群，这在一定程度上挑战了CD56bright NK细胞产生CD56dim NK细胞这种观点。所以作者接下来使用scRNA-seq分析具有GATA2 T354M突变的早期临床无症状女性血液样本中的NK细胞。与健康对照相比，只能发现2％的CD56 bright NK细胞，低于正常的5-10％范围（图9a）。尽管缺乏CD56 bright，但作者认为并不能排除这样的患者中也存在与健康供体CD56 bright亚群相对应的群体，只是细胞表面上缺乏CD56。如t -SNE图所示，健康供体的5个NK亚群仍能被识别出，并多了主要由突变细胞组成的＃6，（图9bc）。与流式细胞结果一致，scRNA-seq显示CD56 bright细胞中只占2％（图9d）。在突变个体和健康对照之间，其他四个亚群的组成相似。这些结果表明，GATA2突变确实导致了CD56 bright NK细胞的丢失，但不是CD56表达的丢失。

接下来，作者集中分析# 6。GSEA结果显示＃6中的细胞周期基因集，以及与细胞增殖相关的转录因子E2F和MYC的富集程度显著富集（图9e）。使用Spearman相关矩阵比较每个供体的每个簇之间的平均基因表达（图9f），突变样本的“Active, Mature,和Terminal NK”与健康对照中的相应亚群有很好的相关性。相反，“ CD56 bright和Transitional NK”相关性较弱。这一结果表明，GATA2突变使CD56 bright NK细胞转录改变比CD56 dim NK细胞更大，这可能是导致CD56 bright NK细胞特异性缺失的原因。

在GATA2 T453M供体来源的NK细胞的前10个最上调的基因中，其中四个属于GIMAP家族（图9f），这个家族对维护淋巴细胞的稳定状态发挥关键作用。GIMAP4已证明可以促细胞凋亡，而GIMAP1和GIMAP5都是抗凋亡。而流式分析也证实了突变个体中凋亡NK细胞更多（图9h）。这里要注意虽然检测到CD56 bright NK的活细胞比例比dim的更多，但当比较突变个体和正常个体时，CD56 bright NK细胞凋亡变化较大。

在GATA2 T453M中表达显著降低的基因中，作者发现“Active NK”的特征基因（包括IEG）是表达水平差异最大的基因。这表明在突变个体中NK细胞稳态激活的潜在缺陷，这可能进一步加剧了具有突变患者的发育缺陷。

图9. GATA2突变导致CD56 bright NK细胞的转录组改变更大

小结

本篇研究利用scRNA测序探索了健康的和GATA2 T354M突变的供体的BM和血液中NK细胞的异质性，还结合CyTOF保证了单细胞数据分析的可靠性。在转录组水平上证实了CD56 bright NK到 CD56 dim NK细胞的发育进程并鉴定出过渡亚群。此外，也证实了GATA2 T354M突变确实特异性地丧失了CD56 bright NK细胞，同时也具有稳态激活的改变。整篇文章的思路主要是：基本的单细胞分析（质控/批次效应/降维/聚类/注释）→ 结合大量文献的NK细胞标记对每一类细胞表征 + GSEA表征亚群特征→ CyTOF验证 + bulk RNA-Seq验证 → Monocle2伪时间分析。其中单细胞结合CyTOF值得大家去学习！

尽管本篇文章无法用独特的分子去表征不同发育阶段的亚群，但仍提供了一定的NK细胞异质性信息，在将来的应用中可以基于本篇研究构建基于机器学习的分类器算法，以表征具有各种遗传或病理状况的NK细胞的发育异质性。