细胞图谱 | Nature | 人类胚胎骨骼发育的多组学图谱

生信菜鸟团

发布于 2024-12-27 17:13:17

3820

文章被收录于专栏：生信菜鸟团生信菜鸟团

Basic Information

英文标题： A multi-omic atlas of human embryonic skeletal development
中文标题：人类胚胎骨骼发育的多组学图谱
发表日期：20 November 2024
文章类型：Article
所属期刊：Nature
文章作者：Ken To | Sarah A. Teichmann
文章链接：https://www.nature.com/articles/s41586-024-08189-z

Abstract

Para_01

人类胚胎骨骼和关节的形成是由新生骨骼中的祖细胞协调分化决定的。
这些细胞的细胞状态、表观遗传过程和关键调控因子仍然不清楚。
在这里，我们对大约336,000个核滴进行了配对转录组和表观遗传学分析，并结合空间转录组学，建立了从受精后5到11周的人类胚胎关节和颅骨发育的多组学图谱。
通过转录组和表观遗传数据的联合建模，我们描述了胚胎骨骼中区域特异性的肢体和颅骨成骨前体轨迹，并进一步描述了调控膜内成骨和软骨内成骨的调控网络。
使用新工具ISS-Patcher对我们的原位测序数据中的细胞簇进行空间定位，揭示了骨骼和关节形成过程中祖细胞分区的机制。
通过轨迹分析，我们预测了人类软骨细胞可能的非经典细胞来源，即施万细胞。
我们还介绍了SNP2Cell，这是一种将特定细胞类型的调控网络与多基因性状（如骨关节炎）联系起来的工具。
利用这里描述的成骨系细胞轨迹，我们模拟了导致单基因颅缝早闭的基因的计算机扰动，并指出了潜在的细胞状态和疾病机制。
这项工作形成了一个详细且动态的骨骼和软骨成熟的调控图谱，推进了我们对人类骨骼发育中细胞命运决定的基本理解。

Main

Para_01

人类骨骼发育始于受精后6到8周（受精后周数，PCW），在胚胎向胎儿阶段过渡期间。
在颅骨中，颅顶祖细胞通过膜内成骨分化为成骨细胞，并在出生后继续容纳成骨祖细胞。
在新生滑膜关节中，肢体芽中的中间区浓缩物在5-6 PCW时出现，并在7至8 PCW之间形成关节腔，不同关节的时间有所不同，在此过程中纤维和韧带结构发展。
软骨支架在滑膜关节的两侧形成，以促进身体平面的发展，直到它们被骨组织替代，因为从8 PCW开始的内骨化作用。
这些区域性的成骨模式支配了整个骨骼的骨发生过程。
据我们所知，它们在单细胞分辨率下的人类发育中的形成和成熟的基础细胞机制尚未完全描述。
为了解决这个问题，我们应用了单核配对RNA（snRNA）和转座酶可及染色质（snATAC）测序（seq）以及空间方法，以揭示5至11 PCW期间胚胎颅骨和四肢中不同骨形成和关节形成生态位成熟的调控景观。
通过这一过程，我们发现了成骨和软骨生成谱系中先前未描述的细胞多样性。
我们开发了ISS-Patcher工具，该工具可以从液滴数据中推断出细胞标签，应用于我们的高分辨率155重原位测序（ISS）数据集，这有助于了解胚胎滑膜关节内的空间定义生态位。
应用OrganAxis，一种新的空间转录组学注释工具，我们还定义了发育中颅骨的空间轨迹。
我们表征了颅面区域的新细胞状态，并进一步阐明了人类施万细胞和成纤维细胞发育的过程。
我们的资源和新的计算工具集，包括SNP2Cell，使得能够预测发育条件（如颅缝早闭）的机制，并允许将特定区域的间充质聚类中的基因调控网络（GRNs）与人类骨骼的老年性疾病（如骨关节炎）关联起来。

Fig. 1: Multi-omics of human embryonic skeletal development.

a, 五个主要解剖区域的解剖取样方法。顶部，颅骨内的起源细胞由前（A）-后（P）轴决定。我们根据这些经典划分进行取样。如图所示，切口用蓝色虚线表示（顶部）。附肢区域的取样如黑色虚线所示（底部）。底部，图表展示了（1）间质的初始浓缩，（2）关节内区细胞确定关节位置，以及（3）关节区域内的软骨形成。
b, 跨年龄和解剖区域取样的样本供体（n=12）概况；图例中表示了图谱模式。采样时间线上的形态事件和变化（左下）如图所示；使用从本研究生成的多组学数据进行的关键分析也显示在右侧。g1，基因；r，区域；TF，转录因子。
c, 膝关节组织切片的苏木精和伊红染色，以说明空化前后的特征。空化前存在关节间区（IZ），空化后通过关节间隙分离，并伴随韧带等软组织结构的出现。CL，交叉韧带；FEM，股骨；MEN，半月板；PAT，髌骨；TIB，胫骨。图像代表来自两个供体的切片。
d, 应用于组织样本的ISS探针检测到的标记基因，如图示意图所示。白色线条勾勒出骨骼区域。星号显示肱骨（HUM）中的软骨和锁骨（CLA）中的骨头。ACR，肩峰。图像代表来自三个供体的样本。
e, 使用从snRNA和snATAC数据计算的MultiVI潜在变量对数据集进行统一流形逼近和投影（UMAP）嵌入。颜色表示细胞簇室。
f, 各解剖位置的相对细胞类型丰度。条形图显示了每个采样解剖区域中细胞簇室的比例，表明间质室在各解剖区域中占主导地位。

Cellular taxonomy of joint development

Para_01

我们在 5 到 11 周胎龄（PCW）期间，对来自 12 名捐赠者的髋关节、膝关节和肩关节进行了配对的单核 RNA 测序（snRNA-seq）和单核 ATAC 测序（snATAC-seq，10X Genomics Multiome）（图 1a,b）。
对于发育中的头骨，据我们所知，以前从未在不同年龄阶段进行过分析，我们分别取样了颅盖和颅底的前部和后部区域，以区分假定的膜内和软骨内骨形成生态位（图 1a 和补充表 1）。
我们在 8 个共享的细胞区室中捕获了 336,162 个高质量的液滴（图 1c-e，扩展数据图 1a-e 和补充表 2）。
在各个区室之间观察到转录组和 ATAC 峰谱的高度一致性（扩展数据图 1f）。
间充质细胞在所有区域中占主导地位，而肌源性细胞在颅盖中缺失（图 1f）。
从这些细胞中，我们定义了超过 100 个细粒度的聚类（补充表 3），并与先前发表的单细胞数据相比，我们的数据中捕捉到了更多样化的软骨形成和骨形成群体（见方法；扩展数据图 2 和补充信息中的讨论）。
为了解析新生滑膜关节中骨系细胞状态及其空间组织，我们对整个完整的早期胚胎前肢（6.7 周胎龄）和后肢（5.7 周胎龄）以及晚期胚胎（7.3 周胎龄）的膝关节和肩关节区域进行了高分辨率的 155 复合 ISS（图 1d，扩展数据图 3a-c 和补充表 4）。
此外，我们对发育中的冠状缝（9 周胎龄）和额骨进行了基于测序的空间转录组学（10X Genomics Visium CytAssist）分析（扩展数据图 3d），从而能够捕捉到空间上的骨系成熟过程。
然后，我们利用这些数据系统地整理了间充质区室内细胞谱系的空间背景。
这使得我们能够在附肢区域和颅底（软骨内骨化）以及颅盖（膜内骨化）中发现骨形成细胞状态，反映了不同的成骨细胞生成机制。
尽管出生后的老鼠在缝线闭合过程中利用软骨内骨化14，但目前尚不清楚这是否是产前老鼠或人类缝线的典型情况。
我们发现颅盖中的软骨形成聚类相对较少（图 3a 和扩展数据图 4a），这与先前报道的膜内骨形成机制一致15。
我们还观察到在液滴和空间数据之间存在更大的细胞丰度差异，尤其是在较晚的发育阶段（补充图 1 和补充信息中的讨论）。

Zonation of the embryonic synovial joint

Para_01

滑膜关节部位的确定发生在四肢中的第5到6周妊娠期，在初始的间充质凝聚体中包含表达GDF5的中间区细胞。
我们对发育阶段进行了差异丰度测试，并确定了InterzoneChon是最早的关节簇（扩展数据图4a-c和补充信息中的讨论）。
InterzoneChon (GDF5+PITX2+)群体的亚聚类和RNA速度动态分析（见方法）使我们能够推断出它们的伪时间轨迹（图2a和补充信息中的讨论）。
我们应用SCENIC+来预测七个中间区亚群之间的基因程序和转录因子可及性变化（图2b和补充表5）。
早期中间区群体（PRRX1）富含与肢体间充质发育相关的转录因子（TBX18、SHOX和LHX9），RUNX2表达和目标基因可及性较低，但SOX5、SOX6和SOX9的表达和目标可及性适中，这表明其轨迹有利于软骨形成而非骨形成。
关节、纤维和GDF5hi中间区簇高度表达了GDF5，并且每个都有独特的基因特征（见补充信息中的讨论）。
我们假设GDF5hi中间区簇，显示出低水平的软骨形成和骨形成转录因子活性，是未分化的，具有维持进入形成关节的能力。
我们利用新开发的ISS-Patcher功能（见方法）来推断155-plex聚类ISS数据流形中的细胞标签（图2c）。
在胚胎后肢中，早期中间区簇弥漫分布在中间区和软骨支架区域，并被真皮中间区包围（图2d）。
关节中间区簇主要富集在即将形成的膝关节软骨区域，该区域与SOX9染色共定位（扩展数据图4c）。
相比之下，表达半月板相关（PTN）和韧带相关（POSTN和SCX）基因的纤维中间区簇（扩展数据图4b），在肱骨关节面相邻的肩部中间区区域丰富（图2d）。
后肢中纤维中间区簇的相对稀缺可能是由于肩部纤维软骨比后肢关节形成得更早。
从这些空间富集模式来看，我们展示了假定关节的分区，显示中间区中心早期软骨形成和抗骨形成转录因子的富集，以及发育中软骨支架的RUNX2富集（扩展数据图3d和4d）。
通过基因富集评分，我们的数据显示关节空洞化发生在分区之后，在第7到8周妊娠期之间发生（补充图2和补充信息中的讨论）。

Fig. 2: Formation of embryonic joints across space and time.

a, 广泛中间区软骨簇的RNA子聚类细胞状态的UMAP上的RNA速度，颜色表示每个中间区簇（顶部）。根据解剖区域着色的UMAP（底部）。
b, SCENIC+预测的转录因子表达（框的颜色）。点的大小显示目标基因的可及性（AUCell），点的阴影（灰度）显示目标基因的表达（GEX AUCell）。
c, ISS-Patcher工作流程示意图（见方法）。KNN，k近邻。
d, ISS-Patcher推测的细胞簇的空间图。显示了后肢的ISS图像，以及叠加的推测的中间区簇（左上角），还显示了膝部中间区内的各个推测的中间区簇（右上角）。还展示了前肢的ISS图像，以及肩部中间区内的推测的关节软骨细胞簇（左下角），和肩部中间区内的各个推测的中间区簇（右下角）。
e, 9孕周时颅骨的苏木精和伊红染色切片的矢状视图（左侧），以及使用10X CytAssist Visium数据剖析的相邻切片的一个区域（右侧）。图像是一个供体的代表。还展示了冠状缝合线的Cell2location结果，显示了每个体素中的簇标签富集情况（右侧）。
f, 在面板e的Visium数据的Visium体素中绘制的标准化转录因子表达。
g, SCENIC+预测的转录因子表达（框的颜色）。点的大小显示目标基因的可及性（AUCell），点的阴影（灰度）显示目标基因的表达（GEX AUCell）在缝合线祖细胞中。
h, 覆盖图显示了在膜内（颅骨间充质到HHIP+前成骨细胞）和软骨内（关节IZ到前成骨细胞）骨化过程中，具有逐渐增强的成骨表型的成骨祖细胞状态周围聚集的单细胞ATAC信号。在每个覆盖图下方，显示了SCENIC+预测的从转录起始位点到增强子之间的环（按重要性评分着色）。还展示了SCENIC+预测的一些增强子结合并调节的选定上游转录因子（左侧）。网络将成骨抑制因子（如TWIST1和LMX1B）通过总体抑制连接与促成骨基因（RUNX2和HHIP）联系起来。

Emergence of fibroblast lineages

Para_01

小鼠胚胎肢体中的成纤维细胞谱系状态源自主要的 HIC1+ 前体群体，该群体对骨软骨成分的贡献最小，并产生产后‘通用’的 PI16+ 细胞群，也在成人人体组织中被识别。
我们试图揭示第一孕期人类关节中成纤维细胞谱系的分类。
在胚胎期（小于 8 孕周）的附肢关节中，我们鉴定出成纤维细胞前体（FibroPRO1）、HIC1+ 间充质（HIC1+Mes）和真皮成纤维细胞（DermFIB1 和 DermFIB2；详见补充信息中的讨论），它们分别围绕新生关节分布，并在四肢中呈弥漫性分布（扩展数据图 5a-d）。
大约在 8 孕周时，FibroPRO2 表达 PI16 和 DPT，这些是产后健康组织中与泛组织周围相关成纤维细胞的标志物（扩展数据图 5e,f）。
BNC2 是一种与肌成纤维细胞相关的转录因子，在 FibroPRO1 和 FibroPRO2 中活性较高，这与其在血管周围组织中的产后表达一致。
此外，YBX1 是一种已被证明可以驱动小鼠胚胎成纤维细胞增殖的转录因子，也表现出富集。
基于发育时间和 RNA 速度（见方法），我们预测 HIC1+Mes 在胚胎期小于 8 孕周时作为早期成纤维细胞前体（扩展数据图 5a,b）。
这里推断它会产生腱细胞、滑膜成纤维细胞、真皮成纤维细胞和肌成纤维细胞（扩展数据图 5a,b）。
HIC1+Mes 显示出多种增殖相关转录因子如 WT1、SOX5 和 FOXC1 的高活性，这些因子与侵袭性和活化的滑膜成纤维细胞表型相关（扩展数据图 5g），并在我们的 RNA-ISS 数据中映射到胚胎肢体（扩展数据图 5h）。
HIC1+Mes 在腱发生相关转录因子如 SCX 和 MKX 中也显示出中等的可及性（扩展数据图 5i），这表明具有腱发生潜能，与小鼠的命运图谱和我们的轨迹分析一致（扩展数据图 5a-d），这一发现可能值得在未来的工作中进行功能探索。

Formation of the cranial sutures

Para_01

在颅骨中，初期缝间充质从孕周7周开始成熟，形成缝合关节。
我们确定了丰富的早期颅骨前体细胞（CranialMes、FacialMes 和 PArchMes；扩展数据图 4a 和 6a-c 以及补充信息中的讨论）和表达 RUNX2 的 SutureMes1 和 SutureMes2，它们也表达了 CTSK、SIX2 和 AXIN2，与小鼠颅骨前体一致，表明它们是膜内成骨谱系的一部分（见补充信息中的讨论）。
胎儿小鼠颅缝的经典标记物（TWIST1、ZIC1 和 ZIC4）和 THBS2 在两个 SutureMes 细胞群中富集。
使用 Cell2location，我们将这些定位到正在发育的冠状缝合关节（图 2e）。
成骨前体细胞（HHIP+PreOB）出现在 SutureMes 细胞群的额骨和顶骨边界处（图 2e 和扩展数据图 6d）。
ALX1 是一种对小鼠颅骨形成和人源细胞中神经嵴细胞（NCC）迁移所需的转录因子，以及许多 NCC 转录因子（PAX3、BMP3 和 TSHZ2）在从前部捕获的颅骨细胞中差异富集（图 1a 和扩展数据图 6b,c,e）。
尽管这些数据表明 ALX1+ 细胞可能具有神经嵴起源，但由于 NCCs 的短暂存在，我们没有捕获到真正的早期胚胎 SOX10+ NCCs。
类似于由关节间区表达的成骨抑制因子，SutureMes 表现出高活性的抗成骨转录因子，这些因子也在小鼠颅缝中富集（TWIST1、LMX1B 和 NFATC2）。
同时，成骨（SP7 和 FOXO1）转录因子活性也很高，表明成骨形成的基因调控网络已经准备好。
与胚胎膝关节的分子梯度相似（图 2f,g），我们在缝合区域观察到了 LMX1B 和 TWIST1 的表达，向两侧骨边缘扩散的同时 RUNX2 富集（图 2f 和扩展数据图 4e），这表明在四肢和颅骨中维持非骨化关节空间的机制相似。
为了揭示围绕关键成骨转录因子 RUNX2 和 HHIP 的位点增强子驱动的基因调控网络，我们结合 ATAC 覆盖度与 SCENIC+ 预测的转录因子-峰和峰-基因链接进行了可视化（图 2h 和补充图 3）。
RUNX2 和 HHIP 均被预测通过针对其位点周围增强子的中间抑制因子受到一组共同的抗成骨转录因子的抑制，包括 LMX1B、TWIST1 和 ALX4，这说明了维持成骨启动平衡的关系。
HHIP 在 HHIP+PreOB 中高度可接近，并通过 TWIST1 间接被 LMX1B 抑制。
RUNX2 在 HHIP+PreOB（膜内成骨）和前成骨细胞（软骨内成骨）中最具可接近性，并通过 TCF12 和 PRRX2 间接被相同的抑制因子靶向。
总体而言，该网络展示了关键成骨调节因子通过多条冗余路径对邻近骨组织的非骨化生态位的协调调节。

Trajectories of skeletal osteogenesis

Para_01

成骨细胞的发生从大约7到8周胚胎期开始，并在8周胚胎期的头颅中明显可见（图3a，扩展数据图4a和补充视频1和2）。
为了研究这一点，我们推断出两条主要的成骨轨迹，它们源自不同的成骨前体，在假时间过程中丰富了成骨转录因子并下调了抗成骨转录因子（图3b，c，扩展数据图6f和7a，补充表6和7以及补充信息中的讨论）。
四肢的内骨化预计起源于四肢间充质（LimbMes），这是一个取样自前肢和后肢的簇，其转录特征与胎儿人类肢体芽中的侧板中胚层（WT1）相似（图3c，扩展数据图6a-c和7a以及补充信息中的讨论）。
CranialMes 和 FacialMes 在颅盖的前部和颅底部分别表现出不同的丰度（扩展数据图4a），并表达了来自神经嵴细胞（NCC）的间充质调节因子PAX3和ALX1（参考文献35）（扩展数据图6a，e）。
我们假设这些簇构成了先前未描述的人类NCC衍生的成骨群体（见补充信息中的讨论）。
缝合发生预计从CranialMes分化成SutureMes1和/或SutureMes2（图3b和扩展数据图7a），形成一个预测的轨迹朝向HHIP+PreOB。
HHIP标记了小鼠成骨冠状缝合间充质（参考文献36），我们在本文中证明，它们在一个独特的人群后代中丰富，该人群富含TWIST1的SutureMes1和/或SutureMes2，并分布在人胎发育的骨化颅骨中（图2f）。
在缝合形成之后，定向分化的渐进波从颅缝向正在发展的骨前沿扩散（参考文献31）。
我们应用OrganAxis（见方法）定义了一个连续的成熟轴，该轴跨越冠状缝至成熟额骨区域。
基于组织学特征的分区，我们评估了沿前后轴的细胞状态映射（图3d）。
在缝合区观察到TWIST1+SutureMes1和/或SutureMes2的富集（1-3；图3e）。
在成骨前沿，成骨前体的组织学特征出现，伴随着HHIP+PreOB的富集。
成骨区的建立与抗成骨（LMX1B和TWIST1）的下调和促成骨（RUNX2，DLX5和SP7）转录因子的上调相吻合，标志着缝合区领土的空间分子开关（图2f，g和3d，e及扩展数据图6f）。
包括IRX5，SOST，SPP1，MMP9和DMP1在内的成骨细胞标志物在远离成骨区的地方达到峰值（参考文献37），并且成骨转录因子的富集与远离冠状缝的轴值对齐（图3e和扩展数据图8a）。
在四肢ISS数据中，我们应用ISS-Patcher识别了LimbMes的空间定位，这预计会在7.3周胚胎期在骨骺分化为前成骨细胞，在新形成的骨中分化为成骨细胞（图3f和扩展数据图6g）。
在更成熟的9周胚胎期颅底中，类似模式在蝶骨中被检测到，其中软骨支架的肥大软骨细胞周围是前成骨细胞（图3g）。

Fig. 3: Intramembranous and endochondral osteogenic niches.

a, 在 8.5 和 10 周胎龄（PCW）时，胚胎和颅骨的光片荧光显微镜清除图像。8.5 周胎龄胚胎的渲染图像的 3D 视图、分割的软骨颅和用于 3D 渲染的网格图像（上排）；两个样本中软骨颅的 COL2A1 染色（中排）；以及整个颅骨上的 COL2A1 和骨桥蛋白（SP7）染色叠加图像，显示在颅底 COL2A1 和 SP7 染色有明显重叠，但在顶骨中 COL2A1 表达有限（下排）。冠状缝（CSs）和矢状缝（SSs）根据标签图例着色。图像是两名捐赠者的代表性图像。比例尺，2,000 微米。MAN，下颌骨；MAX，上颌骨；OCC，枕骨；TEM，颞骨。
b, 力导向（FA）嵌入图显示了顶骨缝间充质中膜内骨化轨迹。嵌入图分别显示了 RNA 速度预测的伪时间轨迹箭头、PCW 元数据和 scFates 计算的伪时间（从左到右）。
c, 如面板 b 所示，力导向嵌入图显示了肢体间充质和前成骨细胞中软骨内骨化的轨迹。
d, 使用 OrganAxis 对 10X Visium CytAssist 图像的额骨（矢状切面）进行空间分箱。轴线跨越冠状缝至额骨（后-前）。从左到右：额骨的苏木精和伊红染色图像、来自 Cell2location 的细胞簇富集、轴值（彩虹尺度）和基于组织学特征的手动分箱。图像是一个捐赠者的三个切面的代表性图像。ANT，前部；INF，下部；OF，成骨前沿；OZ，成骨区；POST，后部；SUP，上部；SZ，缝合区。
e, 面板 d 定义的空间分箱内的富集情况，包括膜内细胞状态（左）、选定标记基因表达（中）和通路富集（右）。EC，内皮细胞。
f, 7.3 周胎龄肱骨（H；冠状切面）和推断的细胞簇显示软骨内骨化作为映射的细胞状态序列。图像是三个捐赠者切面的代表性图像。Cor，喙突。
g, 蝶骨的矢状切面和来自 Cell2location 的细胞状态富集。图像是一个捐赠者的代表性图像。源数据

Angiogenesis in the osteogenic niches

Para_01

内皮细胞萌芽的基础，这是人类膜内骨化生态位中驱动骨生成的关键因素，尚未明确。
在膜内生态位中，尖端细胞、壁细胞和毛细血管内皮细胞沿着骨生成区与成骨细胞和骨细胞逐渐共富集（图3e，扩展数据图8a-c和补充信息中的讨论）。
因此，SutureMes1 和/或 SutureMes2 以及 HHIP+PreOB 高表达 VEGFA 和 VEGFB，这两种分子是血管萌芽的调节因子，并定位于内皮细胞高表达 VEGF 受体基因（FLT1、KDR 和 NPR1）的骨生成前沿（扩展数据图9a）。
这表明 SutureMes1 和/或 SutureMes2 以及 HHIP+PreOB 可能在膜内生态位中促进血管侵入，类似于在软骨内骨化生态位中的软骨细胞。
通过 RNA 原位杂交（RNA-ISH），我们观察到 VEGFA 和 KDR–FLT1 在附肢软骨中共表达（扩展数据图10a-d）。
我们的 Visium 数据还显示，在颅底肥大软骨中 VEGFA 表达丰富，与通过 Cell2location 映射的毛细血管内皮细胞共定位（扩展数据图6h）。
通过对发芽血管生成途径的富集评分，我们观察到沿颅前骨成熟轴的血管生成的空间梯度（图3e），这表明血管生成与骨生成之间存在关联。
使用 NicheNet，我们预测了特定于膜内生态位的骨系-尖端细胞相互作用，并使用 RNA-ISH 识别了这些对的共定位（扩展数据图9b和10e-j及补充信息中的讨论）。
我们使用 CellphoneDB 预测了内皮细胞向骨生成细胞状态的信号交互作用。
尖端细胞被预测通过 NOTCH 信号传导，包括 JAG1/JAG2–NOTCH2（扩展数据图9c和补充信息中的讨论），据报道这可以促进产后周血管成骨前体细胞的分化。
壁细胞和毛细血管内皮细胞表达了配体基因 FGF2 和 RSPO3，据先前描述，这些基因可以通过 FGF2–FGFR2 和 RSPO3–LGR5 促进成骨细胞分化（扩展数据图9d和补充信息中的讨论），以及 THBS1–CD36。
内皮细胞还表达了 CCL14 和 CXCL12——编码 CCR1 和 DPP4 的配体，分别支持体外破骨细胞的募集和分化。
这些空间定义的相互作用表明，VEGFA 和 EPHB2 由骨系细胞在骨前沿招募尖端细胞。
我们推测，血管化的内皮细胞随后促进成熟的骨骼中的成骨细胞分化、骨细胞矿化和破骨细胞募集（扩展数据图6i和9e）。
其他可能调节骨生成的谱系，如神经元，在我们的液滴数据中未被捕获，未来对支配神经元的研究可能会揭示潜在的神经元-骨系相互作用（扩展数据图9f和补充信息中的讨论）。

Developmental chondrocyte heterogeneity

Para_01

发育过程中形成多种类型的软骨，包括透明软骨、纤维软骨和弹性软骨。
我们的数据显示出不同的软骨细胞簇，表现出强烈的区域特异性丰度和基因模块（图4a-c，扩展数据图11a-g，补充表3和8及补充信息中的讨论）。
除了先前描述的人群外，我们还确定了新的集群：两种区域特定的软骨前体细胞（ChondroPro1和ChondroPro2；补充信息中的讨论）和DLK1+软骨细胞（DLK1+Chon: DLK1和CD63）。
ChondroPro1和ChondroPro2分别富集于附肢关节和颅底，并表达成纤维细胞分化标志物（POSTN、COL1A1、PRRX1和TWIST1），这与小鼠早期软骨前体细胞的研究结果一致。
DLK1+Chon在核糖体基因和CD63中富集，后者已在肢体的前肥大层中被识别，而DLK1本身是增殖到前肥大表型的胚胎谱系进展的标志。
空间上，CyclingChon、DLK1+Chon和HyperChon在新生骨骼内依次组织，从干骺端延伸至骨干，即初始初级骨化中心，表明DLK1+Chon处于增殖和前肥大软骨细胞表型之间的过渡状态。
我们还表征了以前未描述的颅面群体，包括面部（FacialChon）和下颌软骨细胞（MandibularChon），它们分别高度表达PAX3和SEMA3D，暗示可能源自神经嵴。
斑马鱼和小鼠的谱系追踪研究表明，施万细胞可以在胚胎发生过程中分化为软骨细胞。
我们的数据捕捉到了多样化的施万室，并揭示了SOX9+内神经施万细胞（SOX9+ enSC），其特征是表达软骨细胞（SOX9、COL9A1、ACAN和COL2A1）和经典施万细胞标记物（MPZ和SOX10；图4d,e和扩展数据图12a-c）。
SOX9+ enSC代表轨迹分析预测的一个终点，源自施万细胞前体，并表达间充质（PRRX1、PRRX2、PDGFRA和TWIST2）和HOX（HOXA9、HOXA10、HOXA11和HOXD10）特征（扩展数据图12d,e）。
通过RNA-ISH和RNA-ISS，我们在发育后肢的髋关节软骨中观察到了广泛共定位的施万细胞标记物MPZ和SOX9转录本，以及少量共表达SOX10的细胞（图4e,f，扩展数据图12f和补充表9）。
相应地，我们的软骨细胞簇不表达SOX10或MPZ（扩展数据图12g），表明SOX9+施万群体存在于软骨中。
最近的一项研究也在发育后肢指的软骨中发现了施万细胞。
尽管需要进行谱系追踪实验以进一步调查，但我们推测，类似于小鼠，施万谱系细胞可能是人类发育中软骨细胞的非典型来源。
我们的液滴数据还捕获了一个PAX7+软骨细胞簇，该簇共表达了软骨细胞和肌细胞的标记物和基因模块，并在计算质量控制去除双峰和环境RNA后仍然存在（补充图4-8及方法和讨论在补充信息中）。
尽管RNA-ISH显示PAX7转录本在肢体软骨细胞中的潜在定位，但信号强度相对较低。
未来的工作，涉及FACS后的转录分析，需要调查这一群体的有效性。

Fig. 4: Chondrogenesis across anatomical regions.

a, 软骨细胞的 UMAP 嵌入图，按亚型着色。
b, 点图显示选定标记基因的平均表达量（点的颜色）和表达细胞的比例（点的大小）。
c, 在不同解剖区域和发育年龄（PCW）中软骨形成细胞簇的比例。
d, 施万细胞谱系状态的 UMAP 嵌入图。SC，施万细胞；SCP，施万细胞前体。
e, 点图显示跨施万细胞状态的软骨形成相关基因的平均表达量（点的颜色）和表达细胞的比例（点的大小）。蓝色基因名称表示 SOX9+ enSC 的标记。
f, 7.3 PCW 髋臼中标志基因的三重 RNA-ISH（RNAscope），显示了 SOX9+ enSC（MPZ+SOX9+）和施万细胞（MPZ+SOX10+）的标记物富集和共定位。图像代表一个供体的切片。

Developmental links to complex traits

Para_01

许多与老化骨骼相关的状况已被链接到生命胚胎阶段的关节和骨骼变化受到干扰。
值得注意的是，与成人骨关节炎相关的增强子相关变异似乎作用于特定解剖区域的调控网络，影响发育过程中的滑膜关节形态。
使用功能基因组关联研究（fGWAS；见方法）方法，我们发现膝关节骨关节炎信号在软骨形成状态中显著丰富，但不包括间区软骨（图5a和补充表10）。
相比之下，髋关节骨关节炎的富集仅观察到两个软骨细胞群（ChondroPro1 和 HypertrophicChon），但在前成骨细胞、成骨细胞和骨细胞中富集。
先前的遗传学研究表明，早期骨骼发育会影响随后发生髋关节骨关节炎的风险，可能是通过调节髋部形状和随后的机械力；而膝关节特异性发现则指向软骨形成的改变。
从涉及fGWAS的簇中获得见解，我们使用SNP2Cell来识别富含骨关节炎信号的簇特异性子-GRNs（见方法；图5b和补充信息中的讨论）。
如预期的那样，对于髋关节骨关节炎，前成骨细胞在整个成骨伪时间轨迹上显示出最大的平均富集（图5c），这与fGWAS结果一致。
我们还确定了髋关节骨关节炎和膝关节骨关节炎的子网络，优先考虑前成骨细胞（图5d）和关节软骨细胞（扩展数据图13a），分别揭示了平衡软骨形成和成骨功能的相似调控途径。
在关节软骨细胞-膝关节骨关节炎网络中，几个非转录因子基因（如COL27A1、PRKCA、SNORC和CRISPLD2），它们在软骨组成和软骨细胞分化中发挥作用，被预测由NFATC1和FOXA3调控（补充信息中的讨论）。
对于髋关节骨关节炎中的髋前成骨细胞，成骨调节因子RUNX2显示显著富集，同时多个NFAT基因（NFATC1、NFATC2、NFATC4和NFAT5）也显示富集，这些基因与额外的转录因子（ZEB1、MAF和TEAD1）一起，暗示了钙调神经磷酸酶和WNT信号通路的作用，这些通路已知在髋部形状形成和骨关节炎中起作用（补充信息中的讨论）。
通路分析还显示，在髋前成骨细胞中，对脂质的细胞反应富集（图5e和补充信息中的讨论）。
总体而言，通过应用fGWAS和我们的新工具SNP2Cell，我们在发育中的软骨形成和成骨单细胞谱型中分别鉴定了膝关节和髋关节骨关节炎GWAS信号的差异富集。

Fig. 5: Links to complex diseases of the skeleton.

a, 骨关节炎（OA）及其替代表型（膝关节和髋关节置换）的全基因组关联研究信号富集（fGWAS），展示了与髋关节相关的单核苷酸多态性与膝关节相关的单核苷酸多态性之间细胞类型富集模式的不同。富集通过每个集群特定的逻辑回归测试（Wald 测试 β1 ≠ 0，自由度 = 1，Benjamini–Hochberg (BH) 校正，错误发现率 (FDR) < 0.1）。BOER，Boer 等；logOR，对数优势比；THR，全髋关节置换术；TKR，全膝关节置换术。
b, SNP2Cell 方法的示意图；分数来自 GWAS 汇总统计数据，细胞簇标记分数被映射并整合到一个基因调控网络 (GRN) 中，突出显示预测在疾病中具有簇特异性作用的富集模块。
c, 跨骨形成簇的每个簇前十个基因的髋关节骨关节炎富集分数（热图；顶部），以及一个箱线图显示中位数富集分数，下四分位数和上四分位数，以及跨 n = 1,967 基因和峰的最大四分位数范围 0.1 的须（底部）。凹槽（中位数周围的 95% 置信区间）是基于簇之间的重叠来近似指导显著差异的。
d, 髋骨关节炎特异性和前成骨细胞特异性的最富集基因和峰的子网络。更亮的颜色对应于相对于随机排列所得分数更高的富集分数。
e, 关节软骨膝骨关节炎和前成骨细胞髋骨关节炎富集分数的基因本体生物过程术语的基因集富集分析。NES，标准化富集分数。源数据

Deciphering monogenic craniosynostosis

Para_01

颅缝早闭是一种先天性疾病，涉及胎儿和出生后发育过程中颅骨骨化和缝合线形成的障碍，导致颅骨缝合线过早融合和成骨祖细胞池耗尽（补充信息中的讨论），从而产生整体发育后果。
为了预测富含受颅缝早闭影响基因的正常发育细胞状态，我们对照了与膜内和软骨内骨化途径相关的伪时间差异表达基因与一个包含超过2,700个已知可引起人类先天性疾病的基因候选数据库（扩展数据图13b、补充表6和11及补充信息中的讨论）。
大多数富集的颅缝早闭（22个中有13个）基因在膜内伪时间轨迹的祖细胞群体中被观察到（SutureMes1、SutureMes2和HHIP+PreOB）。
这些富集基因中的大多数也具有高度可访问性，除了IHH，这表明胚胎期可能受到颅缝早闭的影响。
为了模拟候选颅缝早闭相关转录因子在正常骨形成过程中的扰动效应，我们应用了CellOracle来预测485个检测到的转录因子在膜内轨迹中的体内计算速度变化。
对于TWIST1、MSX2和LMX1B的敲除模拟预测将导致SutureMes2中的高速度变化（图6a-c）。
TWIST1和LMX1B也在冠状缝合线中显示出空间富集（图2f）。
推断的速度变化方向也与导致骨形成的敲除一致（补充信息中的讨论）。
总体而言，这些预测可能有助于了解与转录因子介导的颅缝早闭病理性特征相关的潜在转录效应。
未来功能研究将需要揭示每个基因的机制效应。
接下来，我们重建了优先考虑的转录因子（TWIST1、MSX2和LMX1B）之间的相互作用网络，以解析它们的共调控关系（图6d）。
预测显示这些转录因子之间存在编码和非编码目标的相互调节。
在连接节点中，许多转录因子（SIX1、TCF12、NFIX和ALX4）已知通过功能丧失突变与颅缝早闭相关，这表明该区域存在一个紧密调控的网络，赋予缝合线开放性。
值得注意的是，TCF12已被报道通过与TWIST1形成异二聚体来控制冠状缝合线的发展（补充信息中的讨论），并且在具有双重有害突变的小鼠中观察到了严重表型。
我们还探索了以SOST为中心的正常发育过程中的增强子介导调控，评估ECR5增强子的作用，当其发生突变时会导致Van Buchem病，这是一种骨硬化性发育异常（图6e）。
这些关于正常发育的预测可能为骨系单基因条件的转录因子和增强子驱动细胞模型提供信息。
为了探讨对胎儿骨骼发育的潜在细胞外影响，我们应用了Drug2Cell对我们的成骨簇内的致畸药物靶点基因进行富集评分（见方法；补充表12）。
这显示了在膜内祖细胞和下游成骨细胞或成骨细胞状态，例如SutureMes1和HHIP+PreOB中，已知致畸药物靶点的整体目标富集程度高于软骨内祖细胞（扩展数据图14a和补充信息中的讨论）。
我们的分析允许识别在正常发育过程中表达致畸药物靶点基因的骨系簇，这可能有助于未来功能研究的设计。

Fig. 6: Monogenic conditions affecting bone development.

a, 膜内骨化轨迹的力导向嵌入。箭头显示了使用 cytoTRACE 预测的分化方向，假设轨迹早期的细胞平均表达更多的基因。与 RNA 速度一致，预测从颅间充质细胞向骨细胞的分化。缝合线间充质细胞簇左侧指向主轨迹之外的箭头可能反映了缝合线的维持。
b, 每个细胞簇的转录因子敲除扰动评分的热图，显示每个簇中评分最高的前五名转录因子。较高的评分表明，由扰动引起的转录组变化促进了成骨作用。顶级转录因子包括几个已知在颅缝早闭中起作用的基因（用红色标记）。
c, 显示选定编码转录因子的三个基因：TWIST1、MSX2 和 LMX1B 引发的转录组变化方向的转录因子扰动向量。所有这三个基因的敲除模拟促进成骨作用，而过表达则抑制成骨作用，这表明这些转录因子调节缝合线的维持。
d, 使用 SCENIC+ 推断出的 TWIST1、MSX2 和 LMX1B 之间的相互调控增强子驱动的基因调控网络。圆圈和菱形分别代表基因和具有转录因子结合位点的区域。区域-基因链接根据峰到基因的重要性进行颜色和缩放，而转录因子-区域链接则按转录因子着色。蓝色圆圈表示受调控的目标，米黄色圆圈表示调节器，白色圆圈表示中间基因。紧密连接的网络包含多个这样的基因，并暗示其他基因可能具有类似的功能。已知在颅缝早闭中起作用的转录因子用红色文字标记。
e, 显示 SOST 预测调控的增强子驱动的基因调控网络。一个橙色菱形区域突出了在范布赫姆病中已知发生突变并影响 SOST 调控的增强子。区域-基因链接根据峰到基因的重要性进行颜色和缩放。源数据

Discussion

Para_01

我们展示了一个多组学细胞图谱，该图谱捕捉了人类滑膜和缝合关节形成过程中第一孕期的空间分辨细胞分类。
在小鼠中，GDF5+ 细胞被认为可以产生关节内的所有细胞成分，包括韧带和肌腱。
在这里，我们推断滑膜关节中的成纤维细胞和腱细胞源自胚胎期的 HIC1+ 成纤维细胞，这一群体此前在小鼠胚胎肢体中被报道过，以及一个后来的（PI16+）成纤维细胞前体，其转录与人类出生后的普遍成纤维细胞相似。
这些发现共同揭示了第一孕期发育中的人类胚胎滑膜关节的谱系起源。
我们识别出许多先前未描述的人类颅骨胚胎成骨细胞状态，包括颅骨间充质（HAND2）、面部间充质（PAX3）和前弓间充质（LHX8），它们表达的标记物与小鼠中报告的相似。
我们定义了两个空间分辨的 TWIST1+ 缝合间充质群体，以及表达 HHIP 的前成骨细胞，这些群体与胎儿小鼠缝合中的群体相呼应，并推断了它们在空间生态位中的发育轨迹。
通过使用 OrganAxis，我们展示了颅骨中成骨群体的转变，揭示了内皮细胞招募与膜内骨化之间的关联。
我们显示，肢骨软骨支架中的内软骨生态位来源于表达 ISL1 或 TBX5 的未分化前体，这标志着人类肢芽中的等效群体，并形成前成骨细胞和成骨细胞。
这些推断的轨迹为研究第一孕期人类骨骼发育提供了宝贵的参考（讨论见补充信息）。
通过我们的液滴和空间数据，这些数据捕获了人类胚胎骨骼和关节中施万细胞的发育，我们推断施万细胞可能赋予软骨生成潜力，这与之前描述的小鼠发育相类似。
未来，涉及这些人类细胞谱系分离的功能性研究将进一步了解它们的发育生物学。
利用我们的数据集，我们发现了与单基因和复杂肌肉骨骼疾病的发展联系。
我们在不同发育簇中观察到髋关节和膝关节骨关节炎 GWAS 信号的显著富集差异，暗示髋关节骨关节炎中的成骨作用和膝关节骨关节炎中的软骨生成作用。
最后，通过系统地模拟与颅缝早闭症相关的已知转录因子的体内敲除，我们确定了一个抑制缝合处成骨作用的调节网络，这为先前报道的单基因功能丧失突变如何起作用提供了机制洞察。
这种方法应用于成骨轨迹，对探索涉及胚胎成骨作用的其他疾病的基因效应具有潜在价值。
我们跨区域的多模式人类发育骨骼图谱是理解第一孕期人类软骨和骨骼发育的基础资源。
它还具有告知体外努力以分化成骨细胞和其他间充质细胞状态的潜力。

Methods

Sample acquisition and ethics

样本采集与伦理

Para_01

人类肢体和颅骨组织样本是从选择性终止妊娠中获得的，所有样本捐赠者均签署了知情同意书（东英格兰地区，经剑桥中央研究伦理委员会全面批准）。
简而言之，在解剖过程中，样本被悬挂在 PBS 中，并保持在冰上 -4°C 的温度下。
肩关节、髋关节和膝关节从肢体中整体解剖出来。
对于肩关节，切口在锁骨远端三分之一处近端进行，另一个切口在肱骨颈部远端创建。
对于胚胎期肩关节样本，如果明显的骨骼特征尚未形成，则进行估算以捕捉整个肩胛盂肱关节和肩锁关节。
对于孕周小于 8 周的颅骨样本，每个颅盖和颅底区域分别解剖了两个部分，均在额骨后缘分开。
对于大于 8 周孕周的较老颅骨样本，尽可能地解剖分离额骨、顶骨、蝶骨、筛骨、枕骨和颞骨。
样本最初嵌入最佳切割温度介质中，并在异戊烷-干冰混合物中于 -80°C 冷冻。
然后使用 Leica CM1950 冷冻切片机将冷冻切片切成 10 μm 厚的切片，并放置在 SuperFrost Plus 玻片（VWR）上用于 ISS 或 Visium CytAssist，或直接用于单核处理。
对于用于整体免疫染色的样本，样本来自妊娠终止，所有样本捐赠者均签署了知情同意书。
样本由 INSERM 的 HuDeCA 生物库提供，并符合法国法规的要求。
使用这些组织的全部授权已由法国生物医学研究机构（Agence de la Biomédecine, Saint-Denis La Plaine, France; PFS19-012）和 INSERM 伦理委员会（IRB00003888）授予。

ISS and high-resolution imaging

国际空间站和高分辨率成像

Para_01

ISS 使用 10X Genomics CARTANA HS Library Preparation Kit (1110-02，遵循协议 D025) 和 In Situ Sequencing Kit (3110-02，遵循协议 D100) 进行，这些是 HybRISS63 的商业化版本。
探针面板设计基于四肢细胞状态的倍数变化阈值（补充表 3）。
简而言之，发育中的四肢冷冻切片在 PBS 中的 3.7% 甲醛 (252549, Merck) 中固定 30 分钟，然后用 PBS 洗涤两次，每次 1 分钟，随后进行通透化处理。
切片在 37°C 下用 0.1 M 盐酸 (10325710, Fisher) 中的 0.5 mg/ml 胃蛋白酶 (P7012, Merck) 消化 15 秒（5 周胎龄）或 30 秒（6 周胎龄及以上），然后在室温下再次用 PBS 洗涤两次。
在 70% 和 100% 乙醇中各脱水 2 分钟后，将 9 毫米直径（50 微升体积）的 SecureSeal 杂交腔 (GBL621505-20EA, Grace Bio-Labs) 粘贴到每个载玻片上，并用于保持后续反应混合物。
在缓冲液 WB3 中再水化后，在缓冲液 RM1 中进行探针杂交 16 小时，温度为 37°C。
158 重探针面板包括每基因 5 个环状探针，其序列是专有的（10X Genomics CARTANA）。
切片用 PBS-T（含 0.05% Tween-20 的 PBS）洗涤两次，然后在 37°C 下用缓冲液 WB4 洗涤 30 分钟，再用 PBS-T 洗涤三次。
在 RM2 中进行探针连接 2 小时，温度为 37°C，切片用 PBS-T 洗涤三次，然后在 RM3 中进行滚环扩增 18 小时，温度为 30°C。
在 PBS-T 洗涤后，所有滚环产物 (RCPs) 在室温下与 LM（Cy5 标记混合物与 DAPI）杂交 30 分钟，切片用 PBS-T 洗涤并在 70% 和 100% 乙醇中脱水。
移除杂交腔并将载玻片用 SlowFade Gold Antifade Mountant (S36937, Thermo) 封装。
此时对 Cy5 标记的 RCPs 进行成像作为质量控制步骤，以确认 RCP（"锚点"）的生成，并在解码过程中识别位点。
使用 Perkin Elmer Opera Phenix Plus 高内涵筛选系统在共聚焦模式下进行成像，z 步长为 1 微米，使用 63×（数值孔径 1.15，像素大小 0.097 微米）水浸物镜。
对于通道：DAPI（激发波长 375 纳米，发射波长 435-480 纳米），Atto 425（激发波长 425 纳米，发射波长 463-501 纳米），Alexa Fluor 488（激发波长 488 纳米，发射波长 500-550 纳米），Cy3（激发波长 561 纳米，发射波长 570-630 纳米）和 Cy5（激发波长 640 纳米，发射波长 650-760 纳米）。
成像后，每个载玻片在 PBS 中垂直去盖片（轻轻操作，尽量减少摇动，使盖片‘掉落’以防组织受损）。
切片在 70% 和 100% 乙醇中脱水，然后在载玻片上固定新的杂交腔。
使用 100% 甲酰胺 (AM9342, Thermo) 每分钟新鲜应用 5 分钟去除前一个循环，然后用 PBS-T 洗涤。
条形码标记通过两轮杂交进行，首先是适配器探针池（AP 混合物 AP1-AP6，在后续循环中），然后是测序池（SP 混合物，定制 Atto 425），每次 1 小时，温度为 37°C，其间和之后用 PBS-T 洗涤。
切片脱水，移除腔室，并如前所述将载玻片封装和成像。
这重复了五次以生成完整的七轮数据集（锚点和六个条形码位）。

Whole-mount immunostaining, tissue clearing and image analysis

整体免疫染色、组织透明化及图像分析

Para_01

样本通过在室温下摇动孵育1周于0.5 M EDTA中进行脱钙。
孵育期间中途更换一次溶液。
样品用1X PBS洗涤两次，持续1天。
样品在室温下在H2O中的甲醇浓度逐渐增加（20%，40%，60%和80%）中脱水1小时。
然后，样品在白光（11 W 和 3,000 K°）和滚动摇动（004011000, IKA）条件下，与6%过氧化氢溶液在100%甲醇中孵育过夜。
样品在室温下在甲醇浓度逐渐减少（80%，60%，40%和20%）中再水化1小时，洗涤两次，并在0.2% PBS-明胶与0.5% Triton X-100（PBSGT）溶液中阻断1周。
样品转移至含有主要抗体（osterix，1/500；collagen2，1/500）的PBSGT溶液中，在37°C下以20 rpm摇动孵育14天。
随后，在室温下用PBSGT进行六次1小时的洗涤。
接着，将二抗稀释于PBSGT中并通过0.22-μm滤器过滤。
样品在37°C下在二抗溶液中孵育7天，并在室温下用PBSGT进行六次1小时的洗涤。

Para_02

使用 iDISCO+ 协议使样本透明化。
样本被放置在 TPP（Techno Plastic Products）管中，在旋转搅拌下（SB3，Stuart）用甲醇（20%，40%，60%，80% 和 100%（两次））脱水 1 小时。
甲醇的体积约为样本体积的五倍。
然后将样本在 67% DCM 和 33% MeOH 的溶液中孵育过夜，接着在室温下在旋转器上用 100% DCM 处理 30 分钟，随后放入 100% DBE 中。
透明化的样本使用配备 sCMOS 相机 5.5MP（2,560 × 2,160 像素）的 Blaze 光片显微镜（Miltenyi Biotec）成像，该相机由 Imspector Pro 7.5.3 采集软件（Miltenyi Biotec）控制。
光片厚度为 4 微米，由四个不同波长（488 nm、561 nm、639 nm 和 785 nm）的激光产生。
使用 1× 或 4× 物镜和不同放大率的镜头（×0.6、×1 和 ×1.66）。
样本由 Miltenyi 提供的样品架支撑，并置于充满 DBE 的水箱中，根据样本大小从一侧或两侧由激光光片照明。
采集过程中使用了 LightSpeed 模式，以便在合理的时间内以合适的分辨率获取图像。
3D 图像拼贴块以 10% 的重叠度组合在一起。
图像是以 16 位 TIFF 格式获取的。
图像最初使用 MACS iQ View 软件处理，该软件执行了拼贴块的自动对齐。
堆栈图像使用 ImarisFileConverter 转换为 imaris 文件（.ims）。
为了在 Imaris 中隔离特定结构，我们使用了带有手动选择的表面工具，然后使用表面来屏蔽图像。
图像和视频是通过 Imaris 中的快照和动画功能拍摄的。
Adobe Photoshop（v25.2）用于给缝合区域上色。

Multiplexed smFISH

多重单分子荧光原位杂交

Para_01

冰冻切片使用 Leica BOND RX 进行处理，以自动化进行 RNAscope 多重荧光试剂盒 v2 测定（Advanced Cell Diagnostics 和 Bio-Techne），按照制造商的说明操作。
探针见补充表 7。
染色前，新鲜冷冻切片在 4°C 下用含 4% 多聚甲醛的 PBS 固定 45 分钟，然后依次通过 50%、70%、100% 和 100% 的乙醇脱水，每次 5 分钟。
手动预处理后，自动处理包括用蛋白酶 III 消化 15 分钟，然后进行探针杂交。
使用 Opal 520、Opal 570 和 Opal 650（Akoya Biosciences）、TSA-生物素（TSA Plus 生物素试剂盒，Perkin Elmer）和链霉亲和素偶联的 Atto 425（Sigma-Aldrich）进行酪胺信号放大，以开发 RNAscope 探针通道。
染色切片的成像方法与上述 ISS 相同。

Flow cytometry cell sorting

流式细胞术细胞分选

Para_01

我们按照先前描述的方法对新鲜供体组织进行了全细胞解离。
在细胞提取前，样本组织（大约9周肩关节）被解剖以获取骨样，软组织则通过显微解剖移除。
所得的细胞悬液用DAPI（Invitrogen）进行活细胞染色，使用FGFR3抗体（1:50；Thermo Fisher Scientific）和TACR3抗体（1:50；Thermo Fisher Scientific），以及二级抗体进行染色。
通过BigFoot光谱细胞分选仪（Thermo Fisher Scientific）及其专有软件，对DAPI阳性的单细胞进行DAPI染色的FACS门控。
然后依次对FGFR3和TACR3进行门控，以识别双阳性细胞。
使用人外周血单核细胞作为FGFR3和TACR3的阳性对照，未染色的细胞作为阴性对照。

Image-based ISS decoding

基于图像的ISS解码

Para_01

我们使用了李等人描述的ISS解码流程。
该流程包括五个不同的步骤。
首先，我们使用Perkin Elmer提供的Acapella脚本进行图像拼接，为每个周期的所有通道生成二维最大强度投影。
接下来，我们使用Microaligner（版本1.0.0）基于DAPI信号注册所有周期，使用默认参数。
对于细胞分割，我们利用了一种可扩展的算法，该算法采用CellPose（版本3.0）作为分割方法。
预期的细胞大小设置为直径70像素，并进一步扩展10像素以模拟细胞质。
为了解码RNA分子，我们使用了PoSTcode算法（版本1.0），参数如下：rna_spot_size = 5, prob_threshold = 0.6, trackpy_percentile = 90和trackpy_separation = 2。
此外，我们使用STRtree（版本2.0.6）将解码后的RNA分子分配给分割的细胞，并根据Virshup等人的方法生成AnnData对象。
最后，只有包含超过四个RNA分子的细胞被保留用于下游分析。

Visium processing and library preparation

Visium处理和文库制备

Para_01

Visium CytAssist 空间基因表达分析（10x Genomics）按照制造商的协议进行。
根据与液滴数据相关的骨形成微环境的存在（例如，冠状缝合线）选择感兴趣的区域，并相应地与 CytAssist 机器垫片对齐。
在进行 Visium CytAssist 协议之前的后续对齐中，使用 Hamamatsu S60 幻灯片扫描仪以 40 倍放大率捕捉图像。
文库使用 SpaceRanger（10X Genomics）进行映射。

Single-nucleus isolation and library preparation

单核分离和文库制备

Para_01

从新鲜冷冻样本中通过冷冻切片后机械分离单个核，如先前工作所述。
简而言之，10-μm 切片在匀浆缓冲液（250 mM 蔗糖、25 mM 氯化钾、5 mM 氯化镁、10 mM Tris-HCl、1 mM 二硫苏糖醇、1× 蛋白酶抑制剂、0.4 U μl−1 RNaseIn、0.2 U μl−1 SUPERaseIn 和 0.1% Triton X-100 在无核酸酶水中）中使用玻璃 Dounce 组织研磨器（Merck）进行匀浆。
当组织碎片仍然存在时，样品用松散杵‘A’搅拌 10-20 次，然后用紧密杵‘B’搅拌 10 次。
所得混合物通过 50-μm 细胞过滤器，随后离心（500g，5 分钟），沉淀物重悬于 300 μl 储存缓冲液（1× PBS、4% BSA 和 0.2 U μl−1 Protector RNaseIn）中，并通过 20-μm 细胞过滤器。
核通过台盼蓝溶液染色后在显微镜下观察并评估其活力，并根据制造商的协议进一步处理用于 10X Genomics Chromium 单细胞多组学 ATAC + 基因表达分析。
每个反应加载了目标核回收量为 16,000 个液滴的核悬浮液。
对于一些核样本，来自不同供体样本的混合物在一个反应中混合，并根据基因型进行后期解复用。
cDNA 和最终文库的质量控制使用 Agilent 生物分析仪高灵敏度 DNA 分析进行。
文库使用 Illumina NovaSeq 6000 进行测序，每个液滴的最小测序深度为 20,000 读对。

Data preprocessing

数据预处理

Para_01

测序数据使用 CellRanger-ARC 软件（v2.0.0）比对到人类参考基因组（GRCh38-2020-A-2.0.0）。
来自多个样本供体的混合基因型的 10X Multiome 测序通道中的条形码通过 Souporcell（v2.0）使用 BAM 输出按基因型解复用。
随后，Souporcell 的输出按基因型聚类，以便为每个条形码分配元数据。
Visium 数据使用默认输入设置通过 SpaceRanger（v1.1.0）映射，并使用 10X Genomics LoupeBrowser（v7.0）根据捕获帧标记区域将低分辨率 CytAssist 图像与高分辨率显微镜图像对齐。
对于基因表达数据，应用 SoupX（v1.6.0）去除背景环境 RNA。
对于包含超过 16,000 个液滴的 CellRanger-ARC 调用矩阵（超出目标液滴回收数量），我们将估计的全局 rho 值增加了 0.1 以考虑可能存在的额外环境 RNA。
液滴过滤标准为超过 200 个基因且线粒体和核糖体读取比例低于 5%。
双倍体去除方法如下所述。
对于单细胞 ATAC-seq，我们应用 ArchR（v1.0.2）处理 CellRanger-ARC 的输出。
初始液滴质量控制基于独特核片段的数量、信噪比和片段大小分布进行。
此外，移除了转录起始位点富集得分 < 7 和片段数 < 1,000 的液滴。
使用默认设置丢弃了双倍体。
初始聚类在分辨率为 0.2 的条件下进行，使用迭代潜在语义索引的前 40 个维度。
然后，根据初始聚类结果为每个广泛集群按区域制作伪批量复制。
基于伪批量覆盖率，使用 MACS2（v2.2.7.1）进行峰调用（501-bp 定宽峰）。
然后，获得了细胞-峰计数矩阵并导出。
我们应用 muon（v0.1.2）进行归一化、潜在语义索引降维和聚类分析，使用 BBKNN（v1.5.1）校正解剖区域和样本供体的批次效应，以获得 ATAC 嵌入。
基于基因体周围的染色质可及性计算基因得分。
我们应用 MultiVI（通过 scVI v0.6.8）构建 snRNA-seq 和单细胞 ATAC-seq 的联合嵌入。
我们还应用 EmptyDropMultiome（v1.0.0）重复液滴调用，以识别我们的 Multiome 数据中的含核液滴，减少环境 RNA 噪声（'汤'）。
通过将 EmptyDrops 推广到多组学设置，我们使用最小的液滴创建了一个 RNA 和一个 ATAC 环境 RNA '汤' 配置文件，然后测试每个液滴与这两个配置文件的统计偏差，仅保留与汤配置文件存在统计显著差异的液滴。

Doublet detection

双重检测

Para_01

所有在RNA和ATAC模式中检测到的潜在双重物都被从我们的数据中移除。
对于RNA数据，使用Scrublet（v0.2.3）估计双重物概率，并将得分超过0.3作为阈值。
为了应用严格的双重物阈值，我们进行了改进的Scrublet工作流程，如先前所述。
简而言之，首先独立确定每个10X Multiome（基因表达）通道的CellRanger-ARC计数矩阵的每滴Scrublet分数。
然后通过标准的scanpy工作流程使用默认参数进行过度聚类，直到Leiden聚类。
每个单独的聚类进一步被聚类。
计算了每个聚类的Scrublet分数中位数。
使用以分数中位数为中心、标准差由中位绝对偏差估计的正态分布来计算每个聚类的P值。
经过Benjamini-Hochberg校正后的假发现率调整后，P > 0.65被视为高质量细胞的阈值，因为双重物是显著的异常值。
对于ATAC数据，我们首先应用ArchR中的双重物检测方法来移除疑似ATAC双重物。
此外，通过在单个snATAC-seq文库上运行AMULET（v1.1.0）来表征同型和异型双重物，移除q < 0.01的液滴。

Droplet cluster annotation

液滴簇注释

Para_01

我们采用了一种层次聚类方法，首先对全局集成的 scVI（v0.9.1；隐藏层 = 256，潜在变量 = 52，离散度 = '基因-批次'）RNA 嵌入进行 Leiden 聚类，以获得广泛的聚类。
为了验证这些聚类，我们使用 Celltypist 在胚胎肢体芽中训练模型，并将标签转移到我们的嵌入中进行检查。
我们利用这些信息以及经典标记基因来注释广泛的聚类和子谱系。
对于子谱系（软骨细胞、成纤维细胞、与成骨相关的聚类、施万细胞、免疫细胞和内皮细胞），我们进一步使用 scVI（隐藏层 = 256，潜在变量 = 52 和离散度 = '基因-批次'）对每个子集进行嵌入，并进行 Leiden 聚类（分辨率 = 0.6），随后进行差异表达基因（DEG）分析（方法 = '威克森'）以获得聚类标记。
此外，我们还利用推断的细胞状态的空间位置（如下所述）来指导注释。

Differential abundance testing

差异丰度检验

Para_01

我们应用了 MILO 包（v0.1.1）的 Python 实现进行差异丰度测试（http://github.com/emdann/milopy）。
我们使用 scVI 隐表示来创建相关隔室中液滴的 k-最近邻图，并随后应用 milopy 将液滴分配到重叠的邻域，这些液滴来自多个样本（brc_code）。
每个邻域然后基于多数投票被注释为一个聚类。
我们将前后位置和颅骨-附肢协变量的值二值化，以允许在这类变量之间进行测试。
然后，我们确定了基于 Benjamini-Hochberg 校正的差异丰度的对数倍变化值和错误发现率值。

Spatial mapping using Cell2location

使用Cell2location进行空间映射

Para_01

我们使用我们的注释多组学数据作为输入，对Visium CytAssist体素进行了Cell2location（v0.1.4）去卷积。
样本供体被用作批次变量，每个文库在回归模型中被视为协变量。
为了空间映射，我们根据组织学数据估计每个体素包含30个细胞，并将每个体素归一化的超参数检测alpha值设置为20。

ISS-Patcher

ISS-Patcher

Para_01

ISS-Patcher 是一个软件包，用于基于相关高维数据近似低维数据中未通过实验捕获的特征。
它是在本研究中开发的一种方法，用于使用匹配的单核 RNA 测序（snRNA-seq）数据作为参考，近似 ISS 数据中缺失基因的表达特征。
首先，通过子集化 ISS 池中存在的 155-158 个基因来识别两个数据集之间的共享特征空间，然后对两种模态分别进行归一化到细胞总数的中位数、对数转换和 z 分数计算。
然后，使用 Annoy Python 包为每个 ISS 细胞识别出单细胞 RNA 测序（scRNA-seq）空间中的 15 个最近邻，并将 ISS 中缺失的基因作为 scRNA-seq 邻居的平均原始计数进行插补。

Visium axis annotation using OrganAxis

使用OrganAxis进行Visium轴注释

Para_01

我们的 Visium 颅骨样本使用 TissueTag8 半自动模式进行注释，以生成一个一维成熟轴。
首先根据苏木精和伊红特征手动注释发育中的骨骼区域。
不包括骨骼形成小生境的组织区域被排除在注释之外。
存储的注释类别包括冠状缝的多个区域（0级到2级注释），从最中心的部分开始，一个成骨前沿（3级注释）具有成骨前体细胞和成骨区的组织学特征（4级到7级注释），从组织学上出现的成骨细胞开始。
所有注释均保存为 TissueTag 输出格式，该格式记录了注释分辨率、每微米像素数以及注释名称的像素值解释（例如，0 = ‘缝’）和颜色（例如，‘成骨前沿’：‘红色’）。
为了稳健且高效地将 TissueTag 注释迁移到 Visium 对象，我们首先将 TissueTag 注释从像素空间转移到高分辨率六边形网格空间（15-µm 点直径和15-µm 点间中心距离，点之间没有间隙），基于注释图像中点中心的中位像素值（半径/4）。
接下来，为了生成 Visium 数据的连续注释，对于六边形高分辨率网格中的每个点，我们测量了到每个0级注释结构的十个最近点的平均欧几里得距离，以及到1级注释结构最近点的距离。
所有注释通过注释网格与 Visium 数组中最接近的相应点的邻近性映射到 Visium 点。

GRN analysis

GRN分析

Para_01

使用 SCENIC+83（v1.0.0）管道预测转录因子、潜在靶基因以及具有结合位点的调控基因组区域。
使用 MACS2 调用的峰片段矩阵和在 ArchR73 中处理的相应 RNA 计数矩阵作为输入。
通过将液滴聚类成约 10-15 个液滴的组，并根据它们的 RNA 配置文件聚合计数和片段来创建元细胞。
该管道应用于对应于各个谱系的数据集子集：首先，应用 CisTopic（pycistopic v1.0.2）从片段计数中识别区域主题和差异可访问区域，作为转录因子结合的候选区域。
然后运行 CisTarget（pycistarget v1.0.2）扫描这些区域中的转录因子结合位点，并使用 GRNBoost2（arboreto v0.1.6）基于共表达或可访问性将转录因子和区域链接到靶基因。
与靶基因相关的区域中富集的转录因子基序用于构建转录因子-区域和转录因子-基因调控网络。
最后，根据靶基因表达和靶区域可访问性，使用 AUCell 计算调控网络活性得分，并从中得出调控网络特异性得分。
通过链接单个调控网络构建了转录因子、区域和靶基因（增强子驱动的 GRNs）的网络。
所有子集（成骨、软骨形成、纤维形成、早期关节前体、免疫和施万细胞）的转录因子-增强子-基因链接可以在补充表 9 中找到。

Trajectory analysis

轨迹分析

Para_01

为了在成骨亚室构建伪时间轨迹，非循环液滴被子集化，并使用 X_scVI 作为 palantir 的投影以获得多尺度扩散空间。
在扩散空间上使用 sc.pp.neighbors 生成邻域图，并使用前两个主成分作为初始位置，使用 sc.tl.draw_graph 创建 ForceAtlas2 嵌入。
scFates（v1.0.3）用于预测一个捕捉分化路径的原则图。
力导向嵌入和原则图导出到 R 中，并使用 monocle3（v1.0.0）基于图的测试（Moran’s I）计算沿伪时间的差异表达基因，这允许识别在伪时间任何点上调的基因。
结果通过使用 complexHeatmap（v2.6.2）和 seriation（v1.3.0）包绘制热图进行可视化，在 R 中使用平滑样条（smooth.spline；d.f. = 12）平滑基因表达后。
使用 scvelo（v0.2.3）进行了速度分析。
从未经处理的数据中使用 velocyto（v0.17.17）计算剪接和未剪接读数计数，然后使用 scvelo.pp.moments、scvelo.tl.velocity 和 scvelo.tl.velocity_graph 计算预处理液滴的速度。
使用 cytoTRACE（通过 CellRank（v2.0.2）实现）获得了另一个方向性预测，该预测独立于 RNA 速度（基于分化过程中表达基因数量减少的假设）。

Cavitation enrichment score

空化富集分数

Para_01

为了近似估算空化开始的时间，我们使用 scanpy 中的 sc.tl.score_genes() 计算了一个空化富集得分，该得分为髋关节、肩关节和膝关节间质和肌肉隔室中特定基因集的得分，包括 CD44、HAS2、ABCC5、HMMR、MSN 和 UDPGD，这些基因来源于文献和基因本体论术语，涉及透明质酸生物合成过程和透明质酸合酶活性。
我们排除了在我们的数据中表达水平较低的基因，如 HAS3。
对于通路分析，我们检索了与所有 18,640 个基因本体论术语相对应的基因集，并计算了它们的富集得分与空化富集得分之间的相关性。

In silico transcription factor perturbations

计算机模拟转录因子扰动

Para_01

使用了 CellOracle95（v0.12.0），与通过 scFates85 创建的成骨轨迹一起，以及使用 SCENIC+83 预测的相同细胞的调控子被导入 CellOracle 作为基础基因调控网络。
细胞被聚合到包含 10-15 个细胞的元细胞中，并且使用 CellOracle 按细胞簇拟合了解释转录因子来自目标基因表达的线性模型。
使用特定于细胞簇的模型推断基于调控子的转录因子扰动向量，以预测转录因子过表达和敲除的效果。
然后分别选择相应的起始点，计算膜内和软骨内骨化谱系的扩散伪时间。
使用伪时间梯度来得出基于伪时间的分化向量，并计算伪时间-扰动向量的叉积以获得扰动分数。
这些扰动分数表明，转录因子的计算机扰动是否与沿谱系（成骨）的分化一致或相反。
系统地进行了模拟，过表达和敲除基因调控网络中的所有转录因子。
对于每个转录因子和条件，扰动分数随后按细胞簇平均，并汇总在表格中，以筛选促进或抑制成骨的转录因子。

fGWAS analysis

fGWAS 分析

Para_01

fGWAS分析应用于识别疾病相关的细胞簇，具体描述见（https://github.com/natsuhiko/PHM）。
该模型利用了GWAS的完整汇总统计，将单核苷酸多态性（SNPs）与基因联系起来，并捕捉每个细胞簇中基因表达与疾病关联之间的普遍趋势。
同时，该模型还校正了连锁不平衡和其他相关因素。
我们使用了从EBI GWAS目录、开放目标数据库以及膝关节和髋关节骨关节炎和全膝关节和髋关节置换术中获得的完整GWAS汇总统计数据（https://msk.hugeamp.org/downloads.html；补充表8）。

SNP2Cell

SNP2Cell

Para_01

我们使用了一种网络传播方法，将GWAS汇总统计数据和基于细胞簇标志基因的评分整合起来，以优先确定与疾病相关且特定于细胞簇的转录因子网络子单元。
首先，从下载的汇总统计数据中计算每个SNP的得分，并根据连锁不平衡加权。
然后，将这些得分映射到一个GRN，这里是指使用SCENIC+为相应谱系计算的增强子驱动的GRN。
由于所使用的网络包含转录因子和目标基因，以及作为节点的具有转录因子结合位点的区域，因此SNP得分被映射到基因和区域，代表远端调控元件。
然后，使用带有重启的随机游走（或个性化页面排名）过程在网络中传播这些得分。
这将个体SNP的贡献整合起来，信号在相关的网络节点周围汇聚。
该程序重复了1000次，每次使用随机排列的得分来计算置换检验结果和z分数。
接下来，基于差异表达的每个细胞簇标志基因得分以相同的方式传播，从而为每个网络节点获得细胞簇特异性得分。
然后，将SNP和基于表达的得分按细胞簇组合（如参考文献99所述），通过使用每个节点的最小值。
最终得分被设定阈值，得到的连通组件作为富集子网络。
该方法已被编入一个工具，我们称之为SNP2Cell，可在https://github.com/Teichlab/snp2cell获取。

Cell–cell interactions

细胞-细胞相互作用

Para_01

配体-受体相互作用是使用 CellPhoneDB（v4.0.0）中的‘cpdb_analysis_method.call’推断的。我们包括了每个聚类中超过10%细胞表达的基因。P值大于0.001的推断相互作用被移除。我们使用NicheNet（v1.1.1）来识别软骨内和膜内微环境之间的不同相互作用。我们首先使用Seurat实现的Wilcoxon检验计算了两个微环境中成骨细胞簇和顶端细胞的差异表达基因，并使用每个聚类的最小对数倍变化来总结差异表达的配体和受体。使用排名前1,000的差异表达基因来计算配体活性。我们使用默认设置优先考虑配体-受体链接。排名前十的配体及其得分最高的受体使用热图进行了可视化。

Drug2Cell analysis

药物对细胞影响分析

Para_01

人类（Homo sapiens）的药物和靶基因信息从 ChEMBL 数据库（https://www.ebi.ac.uk/chembl/）收集。
对于致畸药物的目标，我们搜索了针对编码其报告目标的基因的临床批准分子，并从 ChEMBL 数据库中整理出一份包含 65 种临床批准药物的列表，这些药物带有致畸性警告（补充表 6）。
药物评分按先前描述的方法计算。
随后，我们根据广泛的临床用途为每种药物引入了药物类别。
Drug2Cell Python 包可在 GitHub（https://github.com/Teichlab/drug2cell）上获取。

CellHint label harmonization

CellHint 标签协调

Para_01

首先，使用 STARSOLO 将 Zhang 等人的数据集中的 fastq 文件重新映射到一个共同的基因组参考（GRCh38-2020-A-2.0.0），这一流程与 Multiome 数据处理相同。
应用 Cellbender 移除背景计数，这些计数表示为模拟的环境 RNA。
我们将这个矩阵与 Zhang 等人提供的质量控制后计数矩阵中的条形码进行交集，并使用 scVI 将其与我们的 snRNA-seq 数据整合，同时考虑了样本供体和液滴技术（细胞或核）的分类协变量，以及总计数、核糖体和线粒体计数的百分比、以及使用 scanpy 包计算的细胞周期得分‘S_score’和‘G2M_score’等连续变量。
从这些得到的潜在变量用于确定邻域，随后在 UMAP 中进行降维。
Zhang 等人的聚类标签被用作 CellHint 调和的标签，在 cellhint.harmonize() 对齐函数中使用。
使用 Cellhint.treeplot() 来检查和半自动对齐两个数据集之间的标签。
标记基因的基因表达谱用于验证两个数据集之间聚类的对齐情况。

Reporting summary

报告摘要

Para_01

关于研究设计的更多信息，请参阅与本文链接的《自然组合报告摘要》。

Data availability

Para_01

本研究中的高通量原始测序数据可从 ArrayExpress (www.ebi.ac.uk/arrayexpress) 获得，访问号为 E-MTAB-14385。
处理过的 snRNA-scATAC-seq、Visium 和 ISS 数据可供可视化，并可以从 https://developmental.cellatlas.io/skeleton-development 下载。
源数据随本文提供。