文献解读|Sci Adv（13.6）：跨越成年恒河猴大脑的单细胞多组学图谱

复杂的人类思想和行为的细胞和分子起源在很大程度上仍然是一个谜。“人脑的相对大小、细胞数量惊人、或较大的皮质表面积和厚度”等数据无法单独解释人类许多独特的能力，也无法解释伴随人类神经发育、神经精神和神经退行性疾病而来的极其多样和复杂的损伤。

因此，灵长类动物大脑的多样化细胞结构对于理解人类的认知、行为和疾病至关重要。

研究者团队从皮质、皮质下、小脑和脑干的30个不同区域生成了单细胞核转录组分析 (snRNA-seq)数据，并使用三水平单细胞组合索引转录组分析 (sci-RNA-seq3)（图1A），确定了17种分子上不同的细胞类型，称为“细胞类型”（图1B-C），分别注释为：(i)神经元细胞，包括皮质谷氨酸能神经元(CAMK2A)、皮质GABA能神经元（GAD1和GAD2）、篮子细胞（GRID2和SORCS3）、其他小脑神经元（主要是颗粒细胞；GRM4）、中等棘神经元（DACH1、PPP1R1B和BCL11B）、血清素能神经元(TPH2)、多巴胺能神经元（TH和DBH）；(ii) 非神经元细胞，包括小胶质细胞(DOCK2)、少突胶质细胞前体细胞（OPC和VCAN）、星形胶质细胞 （ALDH1A1和GFAP）、少突胶质细胞 （MOG和MBP）、血管细胞 (CFH) 和室管膜细胞 (FOXJ1) 。同时还捕获了4种稀有的细胞群（尚未在其它研究中发现）：三个RBFOX3 + (NeuN +) 神经元样群（标记基因：APOA2、F5、KIR2DL1/2）和一个RBFOX3-小胶质细胞样群体（标记基因：KIR3DL1/2+和P2RY12+/PTPRC+/ENTPD1+）。

通过对同一个体的广泛解剖区域进行采样，他们能够分析30个不同大脑区域的细胞组成——这是迄今为止区域范围最广泛的非人灵长类单细胞脑图谱（图1D-E）。主要细胞类别的分布在性别和半球之间是平衡的，但在不同区域之间差异很大，反映了区域特异性功能背后的细胞组成（图1E）。根据细胞类别组成对脑区进行的无监督分层聚类在很大程度上符合更广泛的解剖学分类，皮质、皮质下、脑干和小脑区域通常聚在一起，基于基因表达的前50个主成分对区域进行聚类时也是如此。在这四个广泛的区域类别中，有两个主要由单个细胞类别组成：在皮层，谷氨酸能神经元是最富集的细胞类型，并且数量超过GABA能神经元几乎4倍（图1E），而小脑几乎全部由小脑神经元组成。

相比之下，皮质下和脑干在细胞组成方面的异质性更大，这些区域的样本包含大致相同比例的谷氨酸能神经元，GABA能神经元和少突胶质细胞。他们根据神经解剖学分组将皮质和皮质下样本进一步细分为“区域亚类”，在这些亚类中，细胞组成的变化更有限（图1E）。例如，在皮层下，中棘神经元 (MSN) 约占基底神经节细胞的一半[伏隔核(NAc)平均= 44.7%;尾状核(CN)平均= 60.0%]，而丘脑则富集GABA能神经元[外侧膝状体(LGN)平均= 55.7%；丘脑内侧背侧核(mdTN)平均= 43.8%;丘脑腹外侧核(vlTN)平均= 28.6%]。

图1. 猴大脑图谱 snRNA-seq 数据集的实验方案和摘要。

(A) 用于 sci-RNA-seq3 实验的活检大脑区域示意图。(B) UMAP 可视化所有 snRNA-seq 分析的细胞。(C)条形图显示了每种细胞类型的log2转化细胞计数（左）、区域特异性评分（中）和区域组成（右）。(D)所有snRNA-seq细胞按细胞类型着色的UMAP可视化。(E)条形图显示各区域细胞类型组成 （左）、谷氨酸能神经元与gaba能神经元的log2转化比（中）、神经元与胶质细胞的log2转化比（右）。

为了表征细胞类别内的异质性，他们对数据集进行了分区，并对 17 个细胞类别中的每一个类别分别重复预处理和聚类。总的来说，使用较粗粒度的聚类标准鉴定了 112 个不同的聚类，并使用更精细粒度的标准鉴定了 397 个不同的聚类，捕获了灵长类动物大脑中的神经元和非神经元多样性（图2A），将较粗粒度的聚类称为“细胞亚型”，将较细粒度的聚类称为“细胞亚聚类”。他们发现谷氨酸能（39 种亚型）和 GABA 能（20 种亚型）神经元主要存在于皮质和皮质下的某些区域[例如海马 (HIP) 和丘脑]，而来自其他非皮质大脑区域（例如小脑和纹状体）的神经元在这些区域内的转录水平不同的并且相对一致（图2A）。这是由于在这些区域中存在大量的特化神经元，包括小脑中的颗粒细胞和浦肯野细胞，以及基底神经节中的中棘神经元。

许多细胞亚型，包括神经元型和非神经元型，具有高度的区域特异性。例如，少突胶质细胞亚型8，这最罕见的少突胶质细胞亚型，绝大多数来源于高髓鞘胼胝体（图2A）。在皮质的特定区域内，细胞亚型组成的差异变得更加细微，需要重点量化。对于皮层中每一个足够丰富的中间神经元和星形胶质细胞亚型（> 100个细胞），他们计算一个区域中细胞亚型构成的log2转换比率，并与所有皮层区域中该亚型的平均构成进行比较（图2B）。

在五种最丰富的中间神经元亚型中，他们不仅注意到所有皮质区域的总体平衡，而且还观察到聚类2 (PVALB+)在枕叶 [初级视觉皮层 (V1)] 中的相对富集，在颞叶内区域减少以及在V1 中的聚类5 (ADARB/PAX6+) 的减少。在颞上沟 (STS) 和颞中视区 (MT)，星形胶质细胞亚型 3 (LUZP2/GPC5+)显著减少，但亚型 6 (KCNIP4/RBFOX1+)富集。

中间神经元是通过释放GABA（γ-氨基丁酸）来抑制控制的主要驱动因素，从而强烈影响神经回路。GABA能中间神经元的不适当发育和随后抑制调节的丧失有助于神经发育障碍，包括癫痫和自闭症。他们专注于区域特异性中间神经元亚型的基因标记。11 种中间神经元亚型主要存在于皮质中，可分为小鼠和人脑之间保守的 4 个主要中间神经元组，以SST、PVALB、VIP和LAMP5表达为标志。与皮质相比，脑干和丘脑的中间神经元亚型分布独特（图2C）。丘脑中间神经元利用前馈抑制来传递和调节丘脑皮层神经元的视觉输入，表达高水平的NTNG1 和 RNF220（图 2D），这表明丘脑一级中继核中存在远距离中间神经元。

接下来，对于脑区域子类水平的星形胶质细胞，他们计算了每个细胞的11个lochNESS评分，每个区域亚类1个，每个评分都量化了给定区域亚类在细胞转录附近的富集情况。然后，确定了细胞邻近区域中最富集的区域亚类，并研究了与聚类分配的亚型特征无关的区域异质性（图2E）。他们还扩展了lochNESS，以确定可以通过lochNESS分数预测给定区域表达的基因。他们使用广义线性回归将基因表达建模为每个细胞中每个区域的lochNESS评分的函数。使用这种方法，他们鉴定了特定区域（例如，枕叶中的TCAF2和FRK）和区域组合（例如，脑干、基底神经节和丘脑中的PGD ）中星形胶质细胞的标记（图2F）。这一策略有助于鉴定更复杂的区域特异性基因表达模式。例如EMID1是星形胶质细胞样NG2细胞亚群的标志，它在皮质的星形胶质细胞中表达较高，但在丘脑、脑干或小脑中不表达。相反，参与保护机体免受RNA病毒感染的ADAP2对丘脑内的一个星形胶质细胞亚群具有高度特异性（图2F）。

因此，与传统的聚类相比，LochNESS可以提供一种更细致的方法来识别区域偏向的细胞亚型和基因表达。虽然本研究中他们关注的是星形胶质细胞，但lochNESS可以迭代地应用于每个细胞类别的一个亚类内的区域，例如，适用于所有皮质区域的所有谷氨酸能神经元或所有皮质下区域的少突胶质细胞。

图2. 整个大脑的细胞亚型分布和变异。

(A)条形图显示区域特异性得分和细胞亚型的组成。(B) 热图显示皮质区域内 GABA 能神经元和星形胶质细胞亚型组成。(C) 按细胞亚型（左）和区域亚类（右）着色的 GABA 能神经元的 UMAP 可视化。(D) GABA 能神经元的 UMAP 可视化，由细胞亚型标记基因表达着色。(E) 星形胶质细胞的 UMAP 可视化。( F ) UMAP 星形胶质细胞的可视化。

为了补充转录组数据集和识别大脑细胞中的关键调控基因组区域，他们应用了三水平单细胞组合索引ATAC-seq (sci-ATAC-seq3)来分析单细胞转座酶可及染色质测序(snATAC-seq)结果和表观基因组数据。

他们首先进行了UMAP可视化分析，鉴定出了42个聚类（图3A-B）。为了整合数据集，他们使用了图形链接统一嵌入(graph-linked unified embedding, GLUE)方法，并生成了4171847个细胞核的统一转录组和表观基因组嵌入（图3C-D）。从snATAC-seq和snRNA-seq数据中捕获的细胞类别的区域分布在区域内（图3E）和总体上（图3F）都高度一致，这表明本项研究中的均质化和核分离方案在相同区域中通过两种模式捕获了相同的异质性细胞群。

他们利用基于 snRNA 的细胞类别注释（图 3G）来探索大脑中细胞类型特异性基因调控的异质性。为此，按预测的细胞类别对所有独特的 snATAC-seq 读取进行分区（图 3G）。

图3. 猕猴脑图谱 sci-ATAC-seq 数据集的生成和细胞类别的识别。

(A) 所有 snATAC-seq 细胞的 UMAP 可视化。(B) 细胞标记物启动子可及性评分。(C) 条形图显示 snRNA-seq、snATAC-seq 和集成数据集的大脑区域的细胞核计数。(D) 集成多模态数据的 UMAP 可视化。(E) Spearman 等级相关系数分析。(F) 散点图显示整个 snRNA-seq 和 snATAC-seq 数据集中的细胞类别比例之间的相关性。(G) 在 (A) 中所示的相同 snATAC-seq UMAP 可视化上可视化的集成衍生细胞类注释。

在具有 snATAC-seq 指定细胞核的 11 个细胞类别中，他们平均识别出每个细胞类别 210572 个峰，范围从小胶质细胞中的 99323 个到皮质 GABA 能神经元中的 425738 个（图4A）。平均而言，对于任何给定的细胞类别，这些峰覆盖了基因组的 7.7%，并且发现 28.8% 的峰距离最近的基因或启动子>2 kb（图4A）。

为了识别每个细胞类别中的候选反式作用调控网络，他们对每组细胞类别特异性峰进行了 转录因子(TF) 结合基序富集分析，将这些峰定义为细胞类别顺式调控元件(cCRE) 的子集，不与其他细胞类别中任何峰重叠。细胞类特异性cCRE高度富集许多可能参与细胞特异性基因调控的TF结合基序（图4B），包括许多先前涉及的基序（图4C），例如，小胶质细胞cCRE包含的核TF SPI1（也称为PU.1）的结合位点比预期多6.6倍（图4B-C）。除了这些典型的例子外，他们还发现了许多区分相对相似细胞类别的基序。例如，来自核呼吸因子(NRF) TF家族的NFE2的TF结合基序在中棘神经元和篮状细胞的cCRE中富集最多，而NEUROD1的结合基序在篮状细胞的cCRE中富集最多，其中该TF是篮状细胞末端分化所必需的，因此，轴突生长和抑制回路的形成。

页面比较了 4 个区域亚类的细胞类别中 369 个 TF 结合基序的可及性及其相应基因的表达，其中 189 个 TF 在基因表达和同源基序的可及性之间显示正相关性，180 个显示负相关（图4D）。在Pearson正相关值或负相关值最大的tf中，有很强的细胞类别特异性激活因子和抑制因子。例如，SPI1通过各种功能遗传学方法确定为阿尔茨海默病的候选基因，SPI1基因在小胶质细胞中高表达，SPI1结合基序高可及性，具有很强的激活作用。相反，NFATC2在小胶质细胞和血管细胞中具有抑制作用，在这些细胞类型中，NFATC2基因的高表达与较低的NFATC2基序结合相关。

图4. 候选调控元件中 TF 结合位点基序的富集。

(A)条形图显示了分别调用的峰值集的汇总统计数据。(B) 热图显示细胞类别之间 TF 结合基序的富集。(C) 六个示例细胞类别最富集的 TF 基序的位置权重矩阵。(D) 散点图显示 TF 结合基序的 snATAC-seq 可及性与四个示例 TF 的区域类别中细胞类别内相应 TF 基因的 snRNA-seq 基因表达之间的相关性。

与 GLUE 调节评分相反，逻辑回归分析使他们能够区分峰和基因之间的正调节相互作用和负调节相互作用。如果两项分析的Padj < 0.05，则认为峰基因对为调节性的（图 5A）。对于每个单元格类，他们还使用正则化逻辑回归和t检验扫描差异可访问峰，检测给定单元格类中的可及性与所有其他单元格类中的可及性，考虑具有不同高可及性的cCRE作为细胞类型特异性基因的候选调节因子（图5A）。

他们将分析重点放在 snRNA-seq 数据集中的 6000 个变异最大的基因上，并测试了位于基因启动子 150 kb 范围内的所有 snATAC-seq 峰（定义为 TSS 向上游延伸 2 kb）。绝大多数候选调控相互作用是正相关的（即，在元细胞逻辑回归中具有正效应大小）——无论峰值是在上游、下游还是与基因的转录起始位点(TSS)重叠。当峰靠近TSS时，GLUE调节得分最高（表明峰与基因特征嵌入之间的相似性较高）（图5B）。对于TSS下游的峰值，GLUE调节分数在所有距离上都保持很高，只有在离TSS更远的地方略有下降（图5B）。

为了说明假定的相互作用图如何有助于研究单个位点水平的调控景观，他们重点关注髓磷脂碱性蛋白(MBP)基因（图5C），该基因编码中枢神经系统髓鞘中最丰富的蛋白质之一，MBP具有一系列剪接同种异构体，并且是少突胶质细胞的典型标记物。MBP位于恒河猴(Mmul_10)基因组的18号染色体上，有8个注释的mRNA亚型(Ensembl)。在人类中，典型的MBP亚型3 (18.5 kDa)在成人髓磷脂中占主导地位。在整体峰集中（所有细胞），94个峰落在MBP启动子的150kb以内，并纳入本项研究得分析。在这些峰中，83个确定为MBP的候选调节因子(crMBP)，其中38个与MBP表达呈正相关。在所有的crMBP中，只有一个不位于MBP基因边界内——然而，它位于可能的启动子区域上游不到2kb处。

通过检测成年少突胶质细胞中 crMBP 和主要MBP转录本之间的基因组距离关系，他们发现所有 16 个与亚型 TSS 重叠或下游的 crMBP 均与MBP表达呈正相关。在位于TSS上游的67个crMBP中，22个与MBP表达呈正相关，而45个与MBP表达呈负相关。其中一些负相关的 crMBP 与其他细胞类型中的 sci-ATAC-seq3 峰相对应，特别是 OPC 和小胶质细胞（图 5C）。

图5. 猕猴大脑图谱中顺式调节相互作用的景观。(A) 概述 cCRE 识别标准的示意图。(B) 根据峰和基因 TSS 之间的最小符号距离划分的基因峰 GLUE 调节分数的分布。(C) 候选调控元件从上到下显示为连接峰与MBP表达的推断调控链接的强度和方向。

最后，他们使用 cCRE 来识别可能增加多基因疾病风险的细胞类型相关调控网络。在将恒河猴 cCRE 提升到人类基因组坐标后，他们使用连锁不平衡评分回归 (LDSC) 工具 检测了疾病性状遗传力的富集。检测了总共 53 种与神经系统疾病、紊乱、综合征、行为或其他特征相关的表型，研究了11个细胞类别中分别为cCRE的细胞类别的富集情况，大致概括了几种已知的细胞类别在神经系统疾病中的作用（图6）。例如，与心脏栓塞性中风或缺血性中风相关的位点仅在血管细胞中富集，而血管细胞在形成和维持血脑屏障中起着至关重要的作用。他们还发现，阿尔茨海默病相关位点仅在小胶质细胞中富集——使用三个独立全基因组关联研究(GWAS)的位点重复了这一结果——这与小胶质细胞增殖和激活在阿尔茨海默病中的突出作用一致。最后，与注意力缺陷/多动障碍(ADHD)相关的遗传位点仅在中等多刺神经元的开放区域中富集。中等棘神经元可能是未来多动症治疗的一个有希望的目标，值得进一步研究。

图6. 候选调控元件中遗传性疾病相关位点的富集。

热图显示了所测试疾病、综合征和表型的细胞类 snATAC-seq 峰中疾病的遗传性富集。

 + + + + +