文献解读|Sci Adv（13.6）：在单细胞水平上对发育中的人类大脑皮层进行多组学分析

人脑的细胞复杂性是通过整个发育过程中基因表达的动态变化而建立的，这种变化部分是由顺式调控元件(CRE)的时空活动介导的。鉴于不同神经发育障碍的发病年龄各不相同，因此检测危险因素对人类大脑发育全谱的影响至关重要。

本项研究对来自 6 个发育时期（妊娠早期中期，妊娠晚期，婴儿期，儿童期，青春期，成年期）的 12 个人类新皮质样本进行了单细胞核转录组分析(snRNA-seq)和单细胞核染色质可及性测序(snATAC-seq)分析（图1A）。

然后，他们从 53185 个单细胞核中获得了染色质可及性和基因表达的联合图谱，为了评估两种模式之间的异同，他们首先对 RNA-seq 和 ATAC-seq 数据集进行独立聚类（图 1B）。总的来说，两种模式都揭示了相同的主要新皮质细胞类型，并且 RNA-seq 和 ATAC-seq 鉴定的细胞类型的细胞身份高度一致（图1C）。

接下来，他们使用加权最近邻 (WNN) 分析对相同单细胞的配对模态进行联合聚类。WNN 分析得出 28 个聚类，包括人脑皮质中的所有主要和次要细胞类型，并进一步分为 15 种细胞类型（图 1D）。通过已知细胞类型特异性标记的基因表达和染色质可及性推断的基因活性一致证实了聚类同一性（图1F）。具体来说，他们发现了表达PAX6的神经祖细胞(NPC)，包括放射状胶质细胞 (RG)（HES5和VIM）和中间祖细胞 (IPC) (EOMES)。他们还鉴定了代表不同发育阶段的三种兴奋性神经元亚型（SATB2、SLC17A7和NEUROD2），其中一种富集于早期胎儿样本的细胞（“EN-fetal-early”；4个聚类），一组用于晚期胎儿样本（“EN-fetal-late”；2个聚类），第三类为产后样本（“EN”；2个聚类）。同样，他们发现了三种抑制性神经元亚型（GAD1和GAD2），其中两种代表内侧神经节隆起 (MGE) 衍生（“IN-MGE”;1个聚类；LHX6）和尾神经节隆起 (CGE) 衍生（“IN-CGE”；1个聚类；VIP和ADARB2）亚型在出生后样本中富集，而其余亚型在胎儿样本中富集（“IN -fetal”；1个聚类）。

此外，他们在新皮质中观察到主要胶质细胞类型的聚类，包括少突胶质细胞祖细胞（OPC；两个聚类；OLIG1和SOX10）、星形胶质细胞（3个聚类；AQP4和GFAP）、少突胶质细胞（3个聚类；MOBP和OPALIN）、小胶质细胞（4个聚类；PTPRC和CX3CR1）以及内皮细胞（1个聚类；CLDN5）、周细胞（1个聚类；PDGFRB）和血管平滑肌细胞(VSMC)（1个聚类；COL1A2）。

样品特异性细胞类型组成在不同发育阶段显著变化（图1E）。在四个胎儿样本中，神经元细胞群占细胞的绝大多数，而出生后样本的非神经元细胞比例要高得多。此外，他们发现大多数神经祖细胞，包括短暂的RG和IPC，只在两个早期胎儿样本中检测到[孕周(GW) 18-19]，这与人类大脑皮层中的大部分神经发生是在妊娠中期(GW20)发生的过程一致，并且这些祖细胞随着皮质发育的完成而开始消失或转化。

图1. 人类新皮质 RNA 表达和染色质可及性的联合单细胞分析。

(A) snRNA-seq 和 snATAC-seq 实验流程。(B) UMAP 可视化。(C) 热图显示两个聚类结果之间细胞的一致性，以 F1 分数衡量。(D) 通过使用 WNN 分析整合两种模式定义的单细胞 UMAP 可视化。(E)各年龄组细胞类型的比例。(F) 点图显示跨细胞类型的选定标记基因表达和染色质衍生基因活性。

超过80%的表达变异归因于启动子和增强子的可及性（图2A），这表明转录异质性与染色质可及性的变化广泛相关。为了将调控元件与其靶基因联系起来，他们接下来使用基于相关性的方法来检查染色质可及性和基因表达之间的成对关系，这导致了7291个显著的峰基因关联[在转录起始位点(TSS)周围±500 kb范围内]，涉及3082个独特的基因。这些链接中的大多数包含与基因表达正相关的峰。正如预期的那样，这些关联在TSS附近富集，并且相关性随着距离呈指数衰减（图2B）。然而，只有22%的峰-基因链接发生在ATAC-seq峰和最近的基因之间，这表明大多数预测的调控相互作用在线性基因组中至少跳过一个基因（图2C），提示单纯应用“最近邻基因”规则来定义调控靶点的缺点。平均而言，大多数基因的表达至少与两个不同的峰相关，而大多数峰预计与单个靶基因相互作用（图2C）。

尽管参与峰-基因连接的大多数基因与一个或两个峰相关，但有一部分基因与相对较多的峰相关，这表明靶基因活性是由作用于广泛染色质结构域的多个因子协同调节的。总之，他们确定了364个调控染色质结构域(DORC)，其中至少有5个显著的峰基因链接与同一基因相关（图2D）。DORC通常与超级增强子相关，超级增强子是在定义细胞身份中起关键调节作用的一大群增强子区域，并且在多种疾病中受到影响。

染色质可及性和基因表达的共变在RNA-seq和ATAC-seq数据中区分了鉴定的细胞类型（图2E），表明DORC基因链接具有细胞类型特异性。GO分析揭示了神经元和神经胶质细胞发育过程的强烈富集（图2F），强调了 DORC 在皮质发育过程中细胞命运决定中的重要作用。

图2. 顺式调控模式的整体和局部特征。

(A) 方差成分分析显示染色质可及性解释了基因表达的变化。(B) 每个峰到连锁基因 TSS 的距离分布。(C) 直方图显示（从左到右）每个基因显著相关的峰数量的分布。(D) 每个基因显著相关的峰数，基因按升序排序。(E) 热图显示 500 个伪bulk样本（列，按细胞类型排序）中 DORC 中关联峰-基因对的染色质可及性和基因表达（行，左：聚合峰可及性，右：关联基因表达）。(F) 与 DORC 相关的基因的前 15 个 GO 分析结果。

在鉴定了从早期胎儿皮质板到成人皮质样本的各种神经元亚型后，他们接下来使用配对的多组学单细胞核图谱来推断整个皮质发生和神经元分化过程中基因调控的发育动力学。通过关注神经元群体（包括 RG、IPC、EN-fetal-early、EN-fetal-late、EN、IN-fetal、IN-MGE 和 IN-CGE）并锚定起始点来进行伪时间轨迹分析。不同的细胞类型根据其发育阶段而拟合在推断的轨迹上（图3A），胎儿样本特异性神经元群体位于初始祖群体和出生后样本的成熟神经元之间（即EN、IN-MGE和IN-CGE）。在起始点后不久，发育轨迹分为EN谱系和IN谱系，IN谱系随后又分为IN-MGE和IN- CGE亚型（图3B）。分配给三种谱系的细胞数量分别为EN谱系(14146)、IN-MGE谱系(5728)和IN-CGE谱系(4904)。他们在神经元群体中重复了峰基因关联分析，得到了涉及930个独特基因的1638个显著关联。在整个谱系中，残差通常为正（图3C），这反映了顺式调控元件的谱系启动，因为DPRC通常在其相关基因表达开始之前就变得可及。此外，他们发现DORC的谱系启动模式随着峰值数量的增加而变得更加稳健，这表明对染色质可及性启动状态的置信度更高。

然后，它们深入研究了EN谱系上的峰值基因链接，该谱系从包括RG和IPC在内的神经元祖细胞开始，然后依次分化为胎儿早期到胎儿晚期再到出生后不同发育阶段特异性的兴奋性神经元亚型（图3A-B）。超过87%的连锁基因的表达水平沿伪时间轨迹显著变化。然后，他们使用k-均值 (km)聚类将这些基因分成4个聚类，每个聚类对应不同的发育时期（图3D），GO分析显示不同时间段发生的独特生物活性（图3E）。具体来说，在轨迹的开始（“km1”），连锁基因在与细胞命运特化、细胞分化的时间调控和神经前体细胞增殖相关的过程中富集。在接下来的胎儿早期(km2)，峰值基因相互作用与神经元迁移、形态发生、突触组织和轴突发生相关。此后，在胎儿晚期(km3)和产后阶段(km4)，兴奋性神经元获得了神经递质运输和调节的能力，表明细胞成熟。

细胞谱系发育转变过程中的动态调控活动是由转录因子的时空模式高度协调的。为了识别控制这些动态调节活动的转录因子(TF)，他们在不同的聚类中进行了 TF 基序富集分析。TF基序在细胞分化和发育的早期阶段富集，包括EN1（多巴胺能亚群的关键调节因子）和HESX1（由参与神经定型的枢纽基因编码）（图3G）。在中期阶段（包括胎儿早期和晚期），相关的峰值基因链接在神经元转录因子（如 NEUROD1、NEUROG2 和 BHLHE22）的基序中富集程度更高，表明最活跃的神经发生过程发生在这些特定的发育时期。在最后一个出生后峰值基因链接中发现富集的 TF 基序较少，包括E2F2等细胞周期调节因子。

Cut-like homeobox 2 (CUX2)是神经元特异性DORC基因，受附近最多的推测增强子调控（图3C），并作为EN谱系(km2)中第二个最早阶段的标记（图3D）。TF NEUROD1的结合基序在km2染色质可及区域中强烈富集，并且NEUROD1活性与CUX2 DORC染色质可及状态高度相关（图3F）。接下来，他们根据推测的km2阶段的假时间对单细胞进行排序，并确定了一个明确的模式，即NEUROD1的活性先于CUX2的DORC染色质状态，然后是CUX2的基因表达（图3H）。此外，随着NEUROD1活性的降低，CUX2的表达速度也相应减慢。这些结果表明，NEUROD1可能是早期神经发生过程中诱导CUX2 DORC可及进而诱导CUX2转录的关键TF。

图3. 神经元发育过程中基因调控的轨迹。

(A) 在神经元亚群内确定的轨迹。(B) 分别沿着兴奋性神经元（“EN-谱系”）和抑制性神经元（“IN-谱系”）的谱系推断伪时间。(C) 染色质可及性和基因表达之间的平均残差与神经元群体内确定的 DORC 中涉及的每个基因的显著相关峰的数量之间的关系。(D) 热图显示基因表达和峰基因链接的 DORC 染色质可及性，这些链接沿着 EN 谱系的假时间显著变化。(F) km2 峰中 TF 基序富集的P值根据 TF 基序活性与CUX2 DORC 评分的 Spearman 相关性分析。(G) EN 谱系的峰基因链接聚类中代表的峰的 TF 基序富集。(H) NEUROD1基序活性和表达的谱系动态先于EN 谱系中的 CUX2 DORC 染色质可及性和基因表达。

他们试图通过在培养的 NPC 中采用 CRISPRi，然后使用 RNAscope 直接对单细胞中的 mRNA 进行成像来验证NEUROD1和CUX2之间的预测因果关系。

在阴性对照实验中，用乱序的引导 RNA (gRNA)处理细胞，观察到CUX2广泛表达，尽管每个细胞的mRNA数量较少，而NEUROD1的表达仅限于较小的细胞子集，但每个细胞的mRNA数量范围更广，包括显示强烈转录爆发的部分（图4A）。在分化后第2周，他们定量了两个基因的每个核内荧光点的频率分布，并发现在整个群体中（通过fano因子测量）NEUROD1的表达比CUX2的变化更大（图4B）。

与对照组相比，NEUROD1失活导致NEUROD1 mRNA下调，而CUX2转录完全抑制（图4A-B）。鉴于对照 NPC 中约 80% 的细胞核不显示NEUROD1转录，NEUROD1启动子可能会长期受到严格抑制，但允许发生强烈的转录爆发。反过来，与对照相比，用 CRISPRi 灭活CUX2导致CUX2 mRNA下调和表达CUX2的细胞比例减少。总之，这些数据表明，尽管NEUROD1表达频率较低，但它是维持CUX2持续转录所必需的（图4C）。

图4. 评估 NEUROD1 和 CUX2 在区分 NPC 中的关系。

(A) RNAscope分析。(B) 所有条件下核 RNA 频率分布的小提琴图。

为了揭示与疾病相关的变异是否在特定细胞类型或发育阶段富集，他们使用连锁不平衡(LD)感知方法来评估相关全基因组范围集合之间的重叠关联研究 (GWAS) 以及谱系定义基因和染色质峰（源自本项研究的多组学单细胞数据）。通过分析 15 种细胞类型中的 9 种神经精神特征和 3 种不相关的控制特征，他们分别确定了细胞类型特异性染色质可及性和转录组数据中的 33 和 28 个显著关联（图5A）。他们观察到有意义的细胞类型- GWAS性状对有高度重叠（41对有意义的性状对中有20对是共享的），这提示两种模式都报告了可靠和有信息的关联。

与之前的研究一致，精神分裂症 (SCZ)、双向情感障碍 (BD)、注意缺陷多动障碍 (ADHD) 和重度抑郁症 (MDD) 富含多种神经元亚型。他们的结果确定了几个尚未由遗传数据描述的关联，包括少突胶质细胞与Tourette综合征(TS)、星形胶质细胞与强迫症(OCD)、 OPC与焦虑、以及胎儿神经元与神经性厌食症(AN)。对于非精神疾病、免疫相关的特征，包括AD、溃疡性结肠炎和炎症性肠病，仅观察到小胶质细胞的富集，因此进一步加强了本项结果的可信度。为了分析神经精神疾病的时间特异性，他们使用表观基因组和转录组数据比较了胎儿和成人神经元的富集，发现两种检测方法的富集具有高度一致性，因此可以将ADHD、AN和自闭症归类为与胎儿期神经元阶段的相关性更强，而对于TS、BD、SCZ和MDD，发现胎儿期和成人期神经元阶段的贡献相等（图5B）。

他们下一步的目标是为疾病相关基因座相关候选功能基因（图5D）。首先，收集了一组491个与神经精神特征相关的全基因组显著变异体，并根据高LD的存在将其扩展至16005个变异体。在这97个位点的152个基因中，7个同时与两种疾病性状相关，将假定的疾病相关变异体与证明显著峰值-基因关联的峰值重叠，从而在97个全基因组显著位点中确定至少一个受调控的基因。17个基因在神经元谱系特化的伪时间轨迹上表现出显著的表达改变(km1/2/3/4)（图5C）。虽然最初的GWAS通常为每个疾病相关基因座指定几个看似合理的基因靶点，但它们的优先级排序大多基于不精确的基于距离的注释。使用本项研究的方法，能够完善预测，并且在某些情况下，提名之前与疾病不相关的靶基因。一个重复发现的例子是DCLK3（编码神经保护激酶）与SCZ和BD的关联，这是之前通过全转录组关联研究(TWAS)和表达数量性状位点(eQTL)方法观察到的（图5E）。值得注意的是，这种关联来源于两个不同峰的假定疾病相关变异的重叠，这两个峰主要存在于成人兴奋性神经元中。

图5. 使用单细胞衍生的标记基因和峰值将与神经精神特征相关的风险变异映射到因果基因。

(A) 神经精神疾病和不相关控制特征中脑细胞类型的遗传性富集。(B) 选定神经精神疾病中胎儿和成人神经元信号的比较。(C)风险变异的候选因果基因子集。(D) 将与神经精神疾病相关的风险变异与其致病基因联系起来的总体策略示意图。(E)标准化的snATAC-seq衍生的pseudobulk轨迹展示了DCLK3基因的复杂细胞特异性调节，预测该基因为SCZ和BD GWAS风险变异的致病基因（rs75968099和rs75968099）。

 + + + + +