Gene expression dynamics of human and mouse craniofacial development at the single-cell level

Domain: chemistry
Published: 2026-03-09
Editorial outcome: major_revision

Abstract
Craniofacial development is a complex process that involves the specification of diverse and
transient cell types. However, our understanding of these processes and the cell-types present
during human craniofacial developmental remains limited. We address this gap in knowledge
through single-nucleus RNA sequencing of human craniofacial development spanning 4 to 8 post-
conception weeks. This resource identifies multiple subtypes of mesenchyme, epithelium, and
cranial neural crest, among other functionally distinct cell types. Extensive comparisons to single-
nucleus gene expression and spatial transcriptomics of comparable mouse developmental stages
reveal functional conservation of most cell types and identification of anatomically distinct cell
subtypes. We find distinct contributions of cellular subtypes to normal facial morphology and
common risk factors for orofacial clefts. Additionally, we find specific ectodermal and epithelial
subtypes whose gene expression networks are most significantly affected by rare de novo
protein-altering variants from orofacial cleft patients. Together our data provide an extensive
resource for understanding human craniofacial development at the cellular level.

Introduction
Craniofacial development orchestrates the formation of the human face through the
interplay of multiple cell types, including mesenchyme, ectoderm, endothelium, and cranial
neural crest cells (CNCCs), which differentiate into a diverse array of tissues such as bone,
cartilage, muscle, skin, and vasculature1-3. Together, these cells and tissues perform essential
functions of the face like respiration, mastication, communication, and sensory perception4-7,
Disruptions to craniofacial developmental processes rank amongst the most common causes of
human congenital anomalies, with orofacial clefts representing a significant portion of global
birth defects8-11. Thus, there exists significant need to understand the molecular, genetic, and
cellular mechanisms underlying human craniofacial development.

Studies utilizing model organisms, particularly mice, have offered key insights into
craniofacial development and abnormalities12-16. However, significant differences exist between
mouse and human craniofacial development, including variations in timing, cellular
contributions, and gene regulatory networks1,13,17-22. Furthermore, human craniofacial features
exhibit evolutionary adaptations that distinguish them from other mammals and primates,
underscoring the necessity for human-specific studies 5,13,21,22.

Advances in single-cell and single-nucleus RNA sequencing (scRNA-seq and snRNA-seq
respectively) technologies have enabled de tailed characterization of cellular diversity and gene
expression during development15,23-27. These tools are particularly valuable for resolving the
dynamics of rare or transient cell populations, such as CNCCs, that play critical roles in
craniofacial formation13. While previous efforts have developed single-cell atlases for murine
craniofacial tissues, corresponding human datasets are limited15,23. Few studies have examined
bulk gene expression during human craniofacial development13,14,28-31, with only two datasets
sampled during the embryonic period of human development13. While mapping of human
genetics findings to mouse craniofacial cell types has indicated potential disease-causing
subtypes32, limited replicates and differences between human and mouse craniofacial
development preclude confident interpretation.

To address these gaps in knowledge, we constructed a time-resolved gene expression
atlas of human craniofacial development by collecting craniofacial tissues from 24 individual
human embryos encompassing six key time points from 4 to 8 post-conception weeks33,34. We
profiled over 95,000 nuclei and identified eight major cell types, including mesenchyme,
ectoderm, endothelium, muscle progenitors, SOX10+ CNCCs and SOX2+ early progenitors.
Comparative analysis with newly generated and previously published23 murine craniofacial
datasets highlighted significant conservation of major cell types and their gene expression
programs. Integration with human spatial transcriptomics further validated the localization of
these subtypes within the developing human face.

Additionally, this study explores the genetic factors shaping human craniofacial features.
By integrating genome-wide association studies (GWAS) with our atlas, we identified cell type-
specific enrichments for genetic variants associated with normal facial variation. We found that
specific subtypes of ectoderm and mesenchyme contribute to different aspects of facial
appearance and shape. We also examined rare variants associated with congenital craniofacial
disorders and found that de novo protein damaging variants identified in orofacial clefting trios
were enriched in genes that specify distinct ectodermal and epithelial subtypes.

This comprehensive atlas provides a high-resolution view of the cellular and molecular
landscape of human craniofacial development, integrating gene expression, spatial mapping,
and evolutionary genomics. Our work enhances our understanding of human craniofacial
biology and establishes a framework for future studies aimed at uncovering therapeutic targets
and evolutionary insights into craniofacial biology. These data can be explored through an
interactive web application accessible via https://cotneyshiny.research.chop.edu/shiny-
apps/craniofacial_all_snRNA_v2/ and at the Chan-Zuckerberg CellXGene Discover resource35.

Results
Single nuclei gene-expression atlas of the developing human face
To characterize the cellular landscape of human craniofacial development we
performed snRNA-seq analysis on the entire craniofacial prominence of 24 individual human
embryos across 6 distinct time points, encompassing major milestones of human craniofacial
development from 4 to 8 post conception weeks (Fig. 1A). Following multiple levels of quality
control and filtering (Methods and Supplementary Fig. 1) we obtained gene expression profiles
from 95,892 nuclei with an average of 3,996 nuclei from each sample (Supplementary Fig. 2A)
and a median of 6,890 counts from 2,988 genes per nucleus (Supplementary Fig. 2B,C). To
compare to our previously published bulk gene expression analysis, we combined data from all
nuclei of a specific stage into pseudo-bulk gene expression profiles. This revealed the greatest
similarity between samples separated by the least developmental time (Supplementary Fig. 3).
Principal component and differential expression analyses of these pseudo-bulk profiles revealed
very similar results to those previously obtained from bulk tissues13. Specifically, samples were
ordered by developmental stage along the first principal component (Fig. 1B, Supplementary
Fig. 4A) and the greatest number of differentially expressed genes were observed between the
earliest and latest comparable timepoints (CS13 vs CS17) (Supplementary Fig. 4B,C).

Clustering and projection of these single nuclei expression profiles in two dimensions
using Uniform Manifold Approximation and Projection (UMAP)36 revealed eight distinct
populations (Fig. 1C). The clusters were contributed to by samples from each of the stages in
very similar proportions (Fig. 1D,E). To characterize the eight major clusters, we first examined
expression of four genes defined by Li et al12, ALX1, EPCAM, CDH5, and FCERG1 as markers of
mesenchyme, ectoderm, endothelium, and immune cells, respectively (Fig. 2A). ALX1 was most
strongly expressed in cluster 0 (mesenchyme, a tissue with dual origin from the mesoderm and
cranial neural crest during development), EPCAM in cluster 1 (ectoderm), CDH5 in cluster 4
(endothelium), and FCERG1 in cluster 6 (immune). We also examined markers of erythrocytes
including SPTA1, ALAS2, RHAG, and HEGMN, which were strongly enriched in cluster 3
(Supplementary Fig. 5). Clusters 2, 5 and 7 did not show signal for any of these genes.

To characterize the remaining clusters, we identified the top 5 genes that were most
strongly differentially expressed between the clusters (Fig. 2B, Supplementary Data 1). In
agreement with the initial cluster identities suggested by the markers of Li et al12, cluster 0
(mesenchyme) was marked by SATB237,38 and RSPO239,40, cluster 1 (ectoderm) was identified by
GRHL241 and TP63, cluster 3 (erythrocytes) was marked by SPTA142,43 and RHAG44,45, cluster 4
(endothelium) showed highly biased expression for FLT146 and MSYCT147, and cluster 6 (immune) 
was marked by MPEG148, CD8649, and CD16349-52. Cluster 2 showed highly biased expression of
MYOG, MYH3, and MYBPH, all genes related to muscle specification and function53-55, while
cluster 5 was marked by FOXD3, a developmental transcription factor which has been linked to
pluripotency maintenance in stem cells and specification of neural crest in multiple species56-59.
This cluster was also marked by high expression of SOX10 but not SOX2 (Fig. 2A). Thus, we
labeled this cluster as SOX2-/SOX10+ putative neural crest. Cluster 7 showed highly biased
expression of SOX260,61, PAX662, and POUE3F263, which are all markers of early progenitor cells
(Fig. 2A). However, these cells did not express SOX10 and thus were labeled as SOX2+/SOX10-
early progenitors. We inferred RNA velocity and transition states to identify trajectories in the
data across the newly identified main cell types. We found that the trajectories emanate from
the SOX2+/SOX10- population, confirming that these cells likely represent an early progenitor
population (Fig. 2C).

To determine biological and disease relevance of each of these clusters, we analyzed the
top 100 marker genes from each cluster for gene ontology, pathway, and disease enrichments.
The genes that defined the mesenchyme were enriched for processes related to extracellular
matrix development and skeletal, cartilage, and roof of mouth development as well as diseases
including malar flattening and isolated cleft palate (Fig. 2D, Supplementary Fig. 6,
Supplementary Data 2). Genes most strongly expressed in ectoderm were enriched for
processes related to tight junctions, cell adhesion, and the plasma membrane and diseases
including Epidermolysis Bullosa, Keratoderma, and Palate fistula (Fig. 2D, Supplementary Fig. 6,
Supplementary Data 2). The marker genes identified in the SOX2-/SOX10+ plausible cranial
neural crest cell (CNCC) population were enriched for processes related to nervous system
development, pathways related to EGR2 and SOX10-mediated initiation of Schwann cell
myelination and diseases related to neuropathies, primary brain tumors, and neural crest
tumors (Fig. 2D, Supplementary Fig. 6, Supplementary Data 2). Marker genes for the
erythrocyte, immune, endothelial, muscle progenitor, and SOX2+/SOX10- early progenitor
clusters also showed enrichment of relevant cell type-specific processes and disease ontologies
(Supplementary Fig. 6, Fig. 2D, Supplementary Data 2).

Characterization of putative human cranial neural crest
To further characterize the SOX2-/SOX10+ presumptive cranial neural crest (CNCC)
cluster, we performed additional analyses using known cranial neural crest markers including
EDNRB, ERBB3, PAX3, SOX10, SPP1, TFAP2B, and ZEB2, all genes well known to be involved in
various aspects of neural crest specification and migration64. While these genes are biased
toward the SOX2-/SOX10+ putative CNCC cluster, they are not exclusively expressed in cells
found in this cluster (Supplementary Fig. 7). Since no single gene was capable of confirming the
identity of this population, we generated a module score of genes from regulatory networks
identified in cultured human and chimpanzee CNCCs 65. We found significantly higher
expression of this module in the SOX2-/SOX10+ putative CNCC cell population (Fig. 2E)
confirming our initial labeling.

Since few studies have been able to identify populations of CNCCs from primary human
tissue66-68, we sought to subcluster this putative CNCC population (n = 3,294) to better
understand the gene expression programs active in these cells. We identified 19 distinct
subclusters from this original population, some of which were biased towards early vs. late time
points (Fig. 3A,B; Supplementary Fig. 8A-C) . We analyzed the genes enriched in each of these
subclusters to label based on functionality. For instance, CNCC cluster 16 was enriched for
melanocyte-specific gene markers MITF and PMEL (Supplementary Fig. 8D), while two other
clusters (clusters 10 and 14) were significantly enriched for the G2M stage of the cell cycle and
labeled as cycling (Supplementary Fig. 8E). Since SOX10 is also a known marker for neural crest
derived Schwann cells69, we generated a Schwann cell module score using marker genes from
human fetal development24. We identified five clusters that were enriched for the Schwann cell
module score (Supplementary Fig. 8F) and were biased towards the later time points profiled
(Supplementary Fig. 8C). To identify additional subclusters we generated additional module
scores from markers of E9.5 mouse neural crest70 labelled by Wnt1Cre/R26RTomato , which consist
of various functionally specific, fate-committed CNCC subtypes. Using this methodology, we
identified one cluster of mesenchymal-committed CNCCs and one cluster of sensory-committed
CNCCs (Supplementary Fig. 8G). We also identified four more punctate clusters, which we
annotated as neural crest like (ncl). The resulting subclusters were not differentially enriched
for any of the previously computed module scores and were labeled according to stage-specific
bias, i.e., early (early.CNCC) vs. intermediate (int.CNCC) (Fig. 3B, Supplementary Fig. 8H). While
SOX10, NR2F1 and NR2F2 were observed in all subclusters, TFAP2B expression was lower in the
intermediate CNCCs and TFAP2A was absent from the int.CNCC3 cluster. Interestingly, SOX9
was specifically expressed in the mesenchymal committed CNCCs (Fig. 3C).

We identified marker genes for each of the putative CNCC subtypes and found 5,279
genes that were differentially expressed (Supplementary Data 3). We highlighted the top five
transcription factor marker genes that distinguished each subtype including MEIS2, MECOM
and PAX3 in early CNCCs; DACH1 and BCL11B in the int.CNCC.1 subcluster; FLI1 and ALX1 in
mesenchymal committed CNCCs, and NEUROG2 in sensory committed CNCCs (Fig. 3D,
Supplementary Fig. 8I, Supplementary Data 3). Identification of enriched gene ontology
categories for each subtype revealed distinct functions (Supplementary Fig. 9, Supplementary
Data 4). Early CNCCs were enriched for early pattern specification and axon guidance, and
intermediate CNCCs were enriched for ATP synthesis and oxidative phosphorylation.
Presumptive melanocytes were enriched for melanin biosynthetic processes and
developmental pigmentation, while Schwann cells were enriched for synapse organization and
sympathetic nervous system development. The ncl1 cluster was enriched for several categories
shared with early CNCCs and Schwann cells suggesting this was an early neural crest type with
possible downstream fate commitment towards a Schwann cell phenotype. The ncl2 cluster
was enriched for processes overlapping with the mesenchymal committed CNCC’s including
hair follicle development and skin development. The ncl3 cluster was enriched for synaptic
signaling, while the ncl4 was enriched for gas transport processes. We also examined the
previously generated CNCC module score (Fig. 2E) across subtypes identified from ectoderm
and mesenchyme clusters, which we describe in further detail below, and found that all
subtypes labeled as SOX2-/SOX10+ putative CNCC have high module scores, while those from
other main cell type clusters are considerably lower (Supplementary Fig. 8J).

To determine if these subclusters have relevant spatial patterns of expression, we
examined recently published spatial transcriptomics on two sections of a human CS13 embyro25
(Methods and Supplementary Fig. 10). We found expression of marker genes that identified
CNCCs versus the other major craniofacial cell types had both broad and restricted patterns of
expression. For example, NR2F1 was broadly expressed across the embryo whereas PAX3,
TFAP2A, and TFAP2B were more regionally restricted to the head and neural tube regions (Fig.
3E, Supplementary Fig. 11A). Genes identified as markers of CNCC subtypes showed unique
patterns of expression with MECOM, a marker of early.CNCC.1 and TWIST1, a marker of
mesenchymal CNCCs, expressed in the craniofacial region as well as putative somites. ALX4 was
generally restricted to the head region and putative frontal nasal process region, while HAND2
was expressed in the putative pharyngeal arch regions, heart, and limb (Supplementary Fig.
11B). We then calculated module scores on these spatial data using the marker genes from
each of the CNCC subtypes. We found that each of these marker genes sets were biased toward
the neural tube region (Supplementary Fig. 11C) consistent with their initial site of specification.

In summary while the SOX2-/SOX10+ cluster is clearly enriched for gene expression
signatures of CNCCs, not all subclusters are CNCCs. Instead, some subclusters represent more
terminally differentiated cell states including Schwann cells and melanocytes. Moreover, since
the neural crest have been somewhat elusive to detect in previous single cell atlases it is
possible other CNCC states are present in the developing embryo that we are unable to sample
here.

Conservation of cell types in craniofacial development
Next, we sought to determine if these subtypes are present in mouse and if similar sets
of genes specify all the main cell types we found. To allow direct comparisons to our human
data, we generated snRNA-seq data from mouse craniofacial tissues harvested from multiple
biological replicates of E10.5 to E12.5, which reflect the major morphological landmarks of the
human tissue profiled. We integrated these data with recently published single cell gene
expression data from E13.5 and E15.523, resulting in 79402 expression profiles (Supplementary
Fig. 12, Methods). When we clustered these data using approaches identical to the human
data, we obtained seven main clusters with remarkably similar cluster ratios and organization in
the UMAP projection space (Fig. 4A, Supplementary Data 5). When we examined gene
expression of the same major markers profiled in human (Fig. 2A), we readily identified the
same major cell types including the muscle and Sox2-/Sox10+ putative CNCC cell types (Fig. 4B).
To determine if these cell types were specified by the same sets of genes, we compared marker
gene identities and found significant sharing of marker genes between the orthologous major
cell types (Fig. 4C). The highest degree of sharing was found between endothelium, followed by
mesenchyme and erythrocytes. While still significant, the lowest degree of marker gene
conservation was observed between the Sox2-/Sox10+ putative CNCC cell type cluster of each
species (Fig. 4C). Despite these substantial similarities, we did not recover a separate Sox2+
progenitor cluster in mouse craniofacial tissue.

Since the majority of main cell types were functionally conserved between human and
mouse, we wondered if we could identify conserved subtypes and infer spatial identities for
human subtypes. To address this question, we leveraged the substantial annotation resources
that have been recently generated for mouse craniofacial development12,14,23. We first
independently subclustered the mouse CNCCs in an identical manner to that of human,
resulting in 12 distinct subtypes (Fig. 4D, Supplementary Fig. 13). When examining the
functional enrichment of marker genes of these subclusters we found similar results as in
human, including a clear population of melanocytes (Supplementary Fig. 14, Supplementary
Data 6). Examination of the same neural crest markers as in human CNCC subtypes revealed
very similar patterns of expression with Sox10 and Nr2f2 expressed across all the subtypes as
well as a similar trend in expression of Foxd3, Pax3, Tfap2a, and Tfap2b across subtypes
(Supplementary Fig. 15A). When we inspected the transcription factor markers for each of
these subtypes, we found many similarities to human including SOX5 in early CNCC’s, NR2F1
and NR2F2 in intermediate CNCC’s, TWIST2 in mesenchymal CNCC’s, and NEUROD1 and
NEUROG1 in sensory CNCC’s (Supplementary Fig. 15B, Supplementary Data 7). We reprocessed
mouse E11.5 spatial transcriptomic data71 (Methods) and found similar expression patterns for
many of the human CNCC markers in mouse tissues (Supplementary Fig. 16A). Examination of
module scores calculated from the mouse CNCC subtypes also revealed similar patterns across
the mouse embryo as observed for human (Supplementary Fig. 16B). To identify the
orthologous CNCC subtypes across species we compared sharing of orthologous marker genes
as we did with the main cell types, as well as an orthogonal method, SAMap72, which compares
cell types with shared expression programs across species. Both approaches showed highest
similarity amongst the human and mouse mesenchymal committed putative CNCC’s as well as
human and mouse melanocytes (Supplementary Fig. 17A,D).

We then subclustered the large mouse mesenchyme cluster and identified 19
subclusters across the mouse timeseries (Supplementary Fig. 18A). Using a combination of gene
ontology enrichments of marker genes and previous annotations of mouse craniofacial single
cell and spatial transcriptomics, we assigned functional and/or positional labels to each cluster
(Fig. 4E, Supplementary Data 8 and 9). For example, the well-established lateral nasal process
(LNP) marker Pax773-75 and the osteoblast marker Sp776,77 were used to define respective
clusters (Supplementary Fig. 18B). We observed two clusters expressing similar osteoblast
markers but were biased in cells from different stages of development, thus we further refined
these as early and late osteoblasts (Supplementary Fig. 18B,C). When we examined mouse
E11.5 spatial transcriptomics data, we found good concordance between marker gene
expression and generalized localization in the embryo (Supplementary Fig. 19A). Furthermore,
projection of modules scores for mouse mesenchymal clusters readily identified specific regions
of the developing craniofacial structures that corresponded well to their assigned labels
(Supplementary Fig. 19B).

We performed identical analyses for the ectodermal cluster revealing an additional 19
subclusters (Fig. 4F, Supplementary Fig. 20A). Applying the same analysis of marker genes from
the literature, gene ontology enrichments, and expression in mouse single cell transcriptomics
data we annotated each of these clusters with functional and spatial labels (Supplementary
Data 10 and 11). We identified highly specific ectodermal populations like periderm marked by
Gabrp78,79 and Grhl380, cells that will form structures of the inner ear marked by Oc9081,82,
palate ectoderm identified by Foxe183, and the putative pituitary progenitors marked with Lhx3
and Lhx484,85 among others (Supplementary Fig. 20B). Several of the ectodermal subtypes
showed stage-specific biases (Supplementary Fig. 20C). As with the mesenchyme, we found the
spatially resolved expression of marker genes corresponded well to expected positions of the
mouse embryo (Supplementary Fig. 21A). Module score calculations for each subtype resulted
in refined spatial identification of subtypes that matched the labels and expected positions well
(Supplementary Fig. 21B). Mouse mesenchyme and ectoderm subtype labels were further
confirmed and validated using multiple databases (Methods; Supplementary Fig. 22).

We next sought to predict the spatial localization of main cell types identified in both
human and mouse by projecting them onto previously published spatial transcriptomics
datasets of the whole mouse embryo from E9.5 to E16.586 and the E15.5 mouse head23,87. We
found expected patterns of cell type localization for human cell types projected onto the mouse
that reflect known anatomical structures (Fig. 4G). Additionally, we projected these cell types
on our recent analysis of spatial gene expression in mouse E15.5 craniofacial sections23 and
similarly found expected patterns of cell type localization (Supplementary Fig. 23A,B). Overall,
the analyses performed here confirmed the identities of multiple cell types and subtypes across
the development of mouse craniofacial tissues. Moreover, the demonstration of conserved
marker genes provides a framework for transferring cell type labels to subclusters identified in
human data, as well as putative spatial inferences from data that originally lacked that
information which we explore below.

Characterization of human mesenchymal cell subtypes
We then subclustered the human mesenchymal population and identified 25 subtype
clusters (Supplementary Fig. 24A). We applied the same methods of assigning cell types and/or
functional labels to each of the human clusters using the gene marker approach based on
orthologous gene expression in addition to our orthogonal method using SAMap (Fig. 5A,
Supplementary Fig. 17B,E). In some cases, multiple human clusters correlated well with a single
mouse cluster and were labelled as separate populations (e.g., mouse palatalShelf2 and human
palatalShelf2-2.4). The most abundant cell types were obtained from the lateral nasal process
(LNP) and the fusion mesenchyme (fusion.mes2) and were well represented from each of the
timepoints. Some of the transient structures like the muscle mesenchyme (muscle.mes) and
cells labeled as palatal.fusion were biased towards early timepoints, while later forming cell
types and structures such as cartilage and palatal shelves were dominated by cells derived from
CS20 samples (Fig. 5B, Supplementary Fig. 24B). When we examined the transcription factor
marker genes identifying each of these clusters we found MSX2 in the fusion mesenchyme 2
cluster (fusion.mes2), TFAP2B in the lateral nasal process population (LNP), SPX in the palatal
fusion zone population (palatal.fusion); HOX genes in the posterior lateral nasal process 2
population (pLNP2); MKX in palatal shelf population 1 (palatal.shelf1); and MYF5 in the muscle
mesenchyme population (muscle.mes), among many others (Fig. 5C, Supplementary Fig. 24C,
Supplementary Data 12).

Gene ontology analysis revealed subtype relevant enrichment such as bone
morphogenesis and ossification in osteoblasts, cartilage development in cartilage, and
embryonic organ development and regionalization in fusion.mes2 populations (Supplementary
Fig. 25, Supplementary Data 13). Examination of disease enrichments across cluster marker
genes revealed some tissue-specific disorders like Osteogenesis imperfecta in late osteoblasts
and Chondrodysplasia in cartilage. Enrichment for genes related to isolated cleft palate were
found in several clusters including cartilage, fusion.mes2, LNP, MandibularArch.1, palatalShelf1,
palatalShelf2, and palatalShelf2.2 (Supplementary Fig. 25F).

We next examined some of the markers that defined the mesenchyme versus other cell
types, such as TWIST1 and SATB2 in the CS13 human spatial transcriptomics data. We found
fairly broad expression of TWIST1 across the embryo with some bias toward the craniofacial
region while, SATB2 was much more regionally restricted to the pharyngeal arches (Fig. 5D,
Supplementary Fig. 26A). We examined subtype marker genes including BARX1, a marker of
MxP.aLNP, and DLX5, a marker of the MandibularArch, were rather specifically localized in the
general region of the pharyngeal arches, while BARX1 was also expressed in the developing
stomach. EMX2, a marker of the MandibularArch, and DMRTA2, a marker of muscle.mes, were
both more restricted to the head region of the embryo at this stage (Supplementary Fig. 26B).

As was observed in mouse, we found much more regionalized signals from module
scores for each subtype. The fusion.mes2 cluster were clearly enriched in the pharyngeal arch
region of the CS13 embryo and biased toward the more anterior portion of this region (Fig. 5E).
The LNP cluster was enriched in distinct areas of the head near the frontonasal region, while
MxP2 and other palatal shelf clusters showed good spatial concordance with the labels that
had been assigned (Supplementary Fig. 26C). When we projected these subtypes onto mouse
spatial data, we found good concordance for cell subtype labels and known anatomical features
(Fig. 4H). Additionally, we projected these subtypes on our recent analysis of spatial gene
expression in mouse E15.5 craniofacial sections23 and similarly found expected patterns of
localization (Supplementary Fig. 23C,D).

Characterization of human ectodermal cell subtypes
Using the same approach as described for the mesenchyme, we identified 22 distinct
ectodermal clusters (Fig. 6A, Supplementary Fig. 27A). Transferring of mouse labels
(Supplementary Fig. 17C,F) revealed cells that give rise to specific ectodermal-derived organs
like the pituitary and thyroid, structures of the inner ear (auditory), and surfaces of several
structures including periderm (Fig. 6A). As with mesenchymal subclusters, many of the
ectodermal subclusters displayed stage-specific bias (Fig. 6B, Supplementary Fig. 27B). Amongst
transcription factor marker genes of ectodermal subtypes, GATA6 was enriched in the palate.1
subtype, SP8 and NR2E1 in the nasal placode (NaP), SALL3, TBX2, and LMX1A in the auditory
subtypes; PAX1 in the fusion.zone cluster, among others (Fig. 6C, Supplementary Fig. 27C,
Supplementary Data 14).

Consistent with our findings for CNCCs and mesenchyme, gene ontology analyses
revealed subtype relevant enrichment including inner ear morphogenesis and development for
auditory subtypes, structural components of the lens in the eye cluster, cornified envelope and
skin development in the periderm cluster, thyroid hormone synthesis in the thyroid progenitor
cluster, and epithelial to mesenchymal transition in the EMT subcluster (Supplementary Fig. 28,
Supplementary Data 15). Ectodermal surface clusters were enriched for a variety of categories
with palate.surface.3 marker genes biased for oxidative phosphorylation and salivary gland
development (Supplementary Fig. 28). Examination of enriched human diseases revealed many
tissue- or region-specific disorders including microcornea and ectopia lentis in the eye cluster,
nonsyndromic deafness in auditory clusters (auditory and auditory.1), thyroid agenesis for the
thyroid progenitor cluster, anterior pituitary hypoplasia for the pituitary progenitor cluster
(Supplementary Fig. 28F).

Examination of ectodermal markers GRHL2 and ESRP1, revealed relatively restricted
expression to various surfaces in the human spatial transcriptomics data, while OC90 was
strongly expressed in putative inner ear (Fig. 6D, Supplementary Fig. 29A). Subtype markers
also showed generally restricted expression particularly for SIX2, NR2E1, SIX6, and NKX2-1
(Supplementary Fig. 29B). When we inspected module scores of each subtype, we observed
exquisitely specific localization for some clusters like eye and auditory. Other clusters were
generally enriched at surfaces of the pharyngeal arches and the putative esophagus
(Supplementary Fig. 29C). As we observed for mesenchyme subclusters, projection onto mouse
spatial data revealed known anatomical features (Fig. 4H).

Cell type specific contributions to craniofacial phenotypes
The analyses above demonstrated strong concordance between human and mouse cell
types and showed strong support of our cell type labelling, giving us the opportunity to
interrogate the cell type-specific expression profiles for enrichment of craniofacial related
genetic signals. The genetic contributions of common variants to many aspects of craniofacial
variation have been studied in multiple populations based on frontal and profile
photographs88,89. However, the cell types and embryonic landmarks that drive these differences
are currently unknown. Following uniform processing of GWAS summary statistics, assigning
variants to genes, and calculating expression weighted cell type enrichments (Methods), we
found that measures likely to be influenced by bone or cartilage structure such as jaw, chin, and
brow protrusion as well the positioning of the eyes relative to the base of the left nostril (EnL-
AIL) and the positioning of the base of the left nostril to the base of the right nostril (AIL-AIR)
were enriched primarily amongst mesenchymal subtypes (Fig. 7A, Supplementary Fig. 30).
For profile landmark measures, we found metrics related to soft tissues including multiple
measures of lip thickness and shape, ear size, and nose shape were enriched primarily in
ectodermal subtypes (Supplementary Fig. 31). Interestingly, while CNCCs give rise to many
downstream cell types and tissues, we found relatively few shape associations for putative
CNCC subtypes. Overall, this analysis suggests specific cell subtypes contribute differentially to
facial variation and suggest these effects begin to manifest early in human development.

Next, we sought to determine if these cell subtypes might be linked to human disease.
We posited that integrating continuous expression patterns in an unbiased fashion might yield
phenotypic associations reflective of their anatomical location or function. To address this, we
employed a systematic examination of the entire Human Phenotype Ontology (HPO)90-92 (Fig.
7B). The ectoderm cluster was enriched for abnormality of the integument (158 of 686
phenotypes) and abnormality of head or neck (98 of 1056 phenotypes) (Supplementary Fig. 32).
Amongst subtypes we found the eye.ect subcluster was strongly enriched for phenotypes
related to abnormalities of the eye (78 of 731), the periderm was most enriched for
abnormalities of the integument, and the ect.GDNF subtype was significantly associated with
abnormalities of the voice. Among the mesenchymal subclusters, many were enriched for
abnormalities of the head and neck as well as abnormalities of the musculoskeletal system. The
cartilage cluster showed the most diverse enrichments including phenotypes related to growth
abnormalities and abnormalities of the musculoskeletal system and limbs. The SOX10+ putative
CNCC cluster was most enriched for abnormalities of the nervous system, while the SOX2+ early
progenitor cluster was most enriched for abnormalities of the eye. Together these results
suggest that some subtypes show phenotypic specificity, while others are more generic.

We next turned to studies of the genetic underpinnings of craniofacial abnormalities,
particularly, the genome-wide association studies that have been systematically performed on
a large cohort of Finnish ancestry93 to determine the cell types that potentially influence risk for
orofacial clefting. When we analyzed these GWAS using the same approach as for facial
variation we found similar partitioning of enrichment between specific classes of cell subtypes
(Fig. 7C). The immune cluster was most significantly enriched for Crohn’s and Systemic Lupus
Erythematosus (SLE), serving as a good control for the specificity of this analysis
(Supplementary Fig. 33). We found that cell subtypes were enriched for ankyloglossia, other
congenital malformation of the tongue and mouth, and cleft palate (Fig. 7C, Supplementary Fig.
33). Specifically, we identified ectoderm subtype annotated as palate.3 as the most significantly
enriched for associations with cleft palate. For mesenchymal subtypes we found significant
associations with Cleft Lip or Lip and Palate in mesenchymal subtypes related to general fusion,
palatal fusion, and the maxillary process surface (Fig. 7C). A few of these subtypes, particularly
the palatal fusion subtypes, were associated with both craniofacial disease and normal facial
variation suggesting their dual contribution to facial shape as well as risk for clefting. When we
examined the spatial mapping of the risk associations of Cleft Palate and Cleft Lip and Cleft
Palate using gsMap94 (Methods), we found significant associations in the pharyngeal regions of
the CS13 human embryo Slice 1 section (Fig. 7D,E) agreeing well with the localization of the
subtype module scores (Fig. 5E, Supplementary Fig. 29C). One of the strongest correlations with
specificity of expression driving the spatial association of the Cleft Lip and Cleft Palate was
MAB21L1 (Fig. 7F), a gene we previously implicated as a putative craniofacial disease-associated
gene13. As a control, we identified localization of the narrow left ventricular outflow tract
obstruction (LVOTO Narrow) trait and highly correlated B4GALNT3 gene specificity to the
putative heart region (Supplementary Fig. 34).

Cell type specific enrichment of disease relevant gene lists
Thus far our analyses of the craniofacial cell types have only leveraged annotated
ontology categories and common variant associations but do not include other gene lists that
may be of use to the craniofacial field. To address this we assembled multiple gene lists
relevant for orofacial clefting including those compiled by CleftGeneDB95, genes co-expressed in
important gene modules or prioritized for craniofacial disease in our recent work1313, and genes
with distinct classes of de novo mutations (synonymous vs protein-altering) in orofacial cleft
trios sequenced as part of the Gabriella Miller Kids First program 96,97 and CPSeq Studies98. We
also curated genes at the extremes of tolerance to loss of function mutations in otherwise
healthy populations that have been suggested to be enriched or depleted of disease-relevant
genes99. Lastly given our findings for common facial variation across humans, we wondered
whether genes potentially regulated by Neanderthal derived sequences might have craniofacial
cell type specific enrichments. As a control for this evolutionary analysis, we included genes
near human accelerated regions (HARs), which have been reported to be enriched in neuronal
related functions and expression100-103.

We again employed the EWCE approach and found that the CleftGeneDB, craniofacial
black co-expression modules, and our prioritized gene lists showed similar patterns of
significant enrichments in mesenchymal subtypes including multiple clusters related to the
maxillary process, palatal shelves, and lateral nasal process (Fig. 8A). Additionally, the genes
identified by gnomAD to have the least tolerance for loss of function mutations (LOUEF decile
1) were also significantly enriched in these mesenchymal subtypes. This contrasted with those
genes identified by gnomAD with the most tolerance for loss of function mutations (LOUEF
decile 9) that showed few enrichments and were generally non-overlapping with the LOUEF
decile 1 enrichments. Interestingly, relatively few ectodermal and CNCC subtypes were
enriched for these gene lists.

When we analyzed the genes near Neanderthal-derived sequences, we found patterns
of cell type enrichment distinct from the more disease-focused lists described above. The
specialized ectodermal EMT and ect.GDNF clusters as well as both mesenchymal maxillary
process surface clusters (MxP.surface and MxP.surface.1), multiple palatal shelf clusters, and
both auditory subtypes were enriched in this analysis. This was contrasted by only a single
progenitor cell subtype enriched in the HAR-associated gene list, consistent with their
previously published association with brain cell types and neuronal function 100-102. We found
no significant enrichments from any of our randomly selected gene lists (Methods) across cell
subtypes and only found enrichment for the gnomAD decile 9 genes in erythrocytes and
immune cell types (Fig. 8A). We then turned to recently identified de novo variants from
orofacial clefting trios from the Gabriella Miller Kids First program96,97 and CPSeq studies98. We
found no enrichment in any cell types for genes affected by de novo synonymous variants (Fig.
8A). In contrast, we identified multiple ectodermal cell subtypes that strongly express genes
with de novo protein altering variants including palate.1 and NaP, as well as multiple palatal
shelf mesenchymal subtypes.

We finally explored whether the total number of de novo variants observed in genes
might reveal additional disease associations. We applied a computational framework that
examines gene lists for excess de novo mutational load104,105 and identified 29 clusters for which
at least one phenotype was significantly enriched for protein-altering variants using a Benjamin-
Hochberg false discovery rate of <10%106 (Fig. 8B and Supplementary Data 16). These include 27
subtypes for all trios with orofacial clefts (OFCs), 21 subtypes for trios with CLP, and 7 subtypes
for trios with CP. We identified 16 significant enrichments across ectodermal cell types (n=22),
5 enrichments from mesenchymal cell types (n=25), and 8 enrichments from the early
progenitor population (n=19). Only 3 clusters were significantly enriched in all three categories
(palate.surface.3, EMT, and NeuralProg.5), none of which were mesenchymal subtypes,
whereas there were 16 shared between all OFCs and CL/P and 4 shared between all OFCs and
CP. We also found 2 clusters that were only significant in the CL/P group (pituitary and
MxP.surface.1) and 4 clusters that were only significant in the full cohort (NeuralProg.4,
MesenchymalLikeProg.1, NeuralProg.2, auditory). We identified strongest enrichment for de
novo variants identified in the whole cohort and those probands with cleft lip and/or cleft
palate (CL/P) in the ectodermal subtypes: palate.surface.3, EMT, palatal.surface.1, NaP,
palate.3, surface3, among others. Interestingly several subtypes were biased toward significant
enrichment related to CP vs CL/P. For instance, palate.1, Prog.1, Prog.2, and cartilage were
enriched for the former while both dental clusters, and multiple palate and palate surface-
related clusters were enriched for the latter (Fig. 8B).

For de novo variants identified in probands with CL/P from significantly enriched
ectodermal subtypes (n=14), the main driver genes were ESRP1, GRHL1, GRHL2 and TPD52.
Variants identified in probands with CP from significantly enriched ectodermal subtypes (n=3),
additionally included CDH1, and LLGL2. These genes have fairly restricted expression in the
head region and presumptive somites of the CS13 embryo (Supplementary Fig. 35).
This enrichment highlights the importance of these cells in OFC etiology, but the difference in signal
drivers may provide insight into the heterogeneity of the genetic architecture between CP and
CL/P.

Methods
Human tissue samples
The use of human embryonic tissue was reviewed and approved by the Human Subjects Protection
Program at UConn Health (UCHC 710-2-13-14-03) and Children’s Hospital of Philadelphia (IRB 24-
022258). Human embryonic craniofacial tissues were collected via the Joint MRC/Wellcome Trust
Human Developmental Biology Resource (HDBR) under-informed ethical consent with Research
Tissue Bank ethical approval (18/LO/0822 and 18/NE/0290, project 200225). Donors are asked to
give explicit written consent for the fetal material to be collected, and only after they have been
counseled about the termination of their pregnancy. Donors do not get anything in return for their
donation and receive the same hospital care as those who decide not to donate. Further
documentation of all policies and ethical approvals for HDBR sample collection can be found
at https://www.hdbr.org/ethical-approvals. Tissues were flash-frozen upon collection and stored at
−80 °C. Upon thawing, the samples were quickly inspected for intactness of the general craniofacial
prominences and the whole sample was processed for single nucleus multiomics. As these are not
entire embryos, we lack many of the additional anatomical landmarks that are available in intact
mouse embryos. Thus, we expect that additional tissues and/or cell types will be present in human
samples that can only be removed through computational filtering.

Mouse embryonic tissue samples
The use of mouse embryonic tissues was reviewed and approved by the UConn Health Institutional
Animal Care and Use Committee (Protocol AP-2000061-0723). Eight-week-old wild-type male and
female C57BL6/J mice were obtained from Jackson Laboratory. Mice were housed according to
recommendations by Jackson Laboratory with 12 h light:dark cycle beginning at 7 a.m. The ambient
temperature was maintained between 20 and 22 °C and humidity was maintained at 40–60%. Mice
were given ad libitum access to food and water. Timed matings were established by the identification
of vaginal plugs the morning following the housing of a single male with multiple female mice.
Pregnant mothers were euthanized by carbon dioxide narcosis at mid-day either 10, 11, or 12 days
after identification of the vaginal plug and embryos were harvested. The staging was confirmed by
counting somites and comparing overall morphology to the Theiler Staging Criteria132. All embryos
from a given litter were combined for individual biological replicates, and at least three biological
replicates were collected and processed for each stage. Craniofacial prominences were collected in a
very similar fashion to human samples, flash frozen, and stored at -80C until prepared for single
nucleus multiomics.

Single-nucleus multiomics
Primary human craniofacial prominences excluding the eyes and more anterior tissue from CS12,
CS13, CS14, CS16, CS17 and CS20, each stage represented by a minimum of 3 replicates, were
obtained from HDBR. Tissue from each embryo was mechanically broken into single-cell suspensions
and cells were checked for viability and counted using Trypan blue staining following the 10X
Genomics protocol for single-cell multiome sequencing using the ChromiumX controller. Samples
were sequenced on multiple Illumina NovaSeq runs according to 10X Genomics recommendations.
Raw fastqs were processed using CellRanger ARC133 (v2.0.2) using hg38 genome and gene
annotations provided by 10X Genomics. For purpose of this paper, we have focused on gene
expression to identify cell types in depth and understand their biological relevance. The chromatin
accessibility data was of generally high quality based on reports from CellRanger ARC. This data forms
the foundation of a separate manuscript that will leverage the cell-type and subtype labels generated
and validated in this study.

Primary mouse craniofacial prominences excluding the eyes and more anterior tissue from E10.5-
E12.5 from multiple (3-5 depending on stage) mixed sex C57BL/6J Mus Musculus embryos (Jackson
Laboratories) were pooled. Animals were raised and sacrificed in compliance with UConn Health
IACUC approval (protocol AP-200061-0723). Samples were mechanically broken into single-cell
suspensions, processed for multiome using the ChromiumX controller, and sequenced in the same
fashion as for human samples above. Raw fastqs were processed using CellRanger ARC (v2.0.2) using
mm10 genome and gene annotations provided by 10X Genomics.

Processing of snRNA and identification of major cell types.
Raw barcode matrices from each human samples generated by CellRanger ARC were utilized for
downstream processing to maximize total nuclei. Raw barcode matrices for each sample were
individually loaded with Read10X_h5 command in Seurat 134 (v. 5.2.1). Subsequently, the emptyDrops
function within the DropletUtils R package (v. 1.24.0) was used to identify nuclei containing ambient
RNA for each sample (emptyDrops FDR < 0.001). Individual sample Seurat objects were then merged
into a single object. Percentage of mitochondrial reads were calculated for each cell and filtering was
performed to only retain cells with less than ten percent mitochondrial derived. Further filtering was
performed based on number of counts per cell (500 < x < 25000) and number of genes detected per
cell (500 < x < 7000), resulting in a total of 130,427 nuclei. Filtered data were normalized with default
values and cell cycle scores were calculated using Seurat. Data was scaled based on S and G2M score
regression and dimensionality reduction with principal component analysis (PCA) were performed
using respective commands in Seurat. The top 2000 variable features were identified, and sample-
level batch correction was performed using the harmony R package135 (v. 1.2.3). Nearest neighbors
based on harmony corrected embeddings were calculated with up to 30 dimensions and clusters
were identified with multiple resolutions from 0.1 to 1 in Seurat. Resolution parameters were
selected for clustering based on generated elbow plots and identification of the elbow point to
determine the optimal number of clusters. We then performed Uniform Manifold Approximation and
Projection (UMAP)36 dimensionality reduction using harmony corrected embeddings in Seurat
(dimensions = 30, minimum distance = 0.3). Resulting clusters were inspected for expression of
multiple craniofacial markers from Li et al 201912 and marker genes were identified for each cluster.
Module scores were generated using the AddModuleScore function in Seurat for neuronal markers
and radial glia markers from Eze et al 2021136. Cells from clusters identified with high neuronal or
radial glia module were removed (See Supplementary Note), resulting in 95,892 nuclei. The process
of normalization, harmonization, and clustering was repeated with remaining cells from all samples.
We observed that early samples had consistently higher mitochondrial reads (Supplementary Fig.
2D), however, this was not due to increased apoptosis as tested by generating an apoptosis module
score using gene markers from the Hallmark Apoptosis gene set within the human MSigDB Core
Collection (Supplementary Fig. 2D). The sex of each embryo was also identified using known sex
determination genes for males (UTY, DDX3Y, KDM 5D, and NLGN4Y) and females (XIST and TSIX)
(Supplementary Fig. 2E), and each of the main cell type clusters contained roughly the same number
of male and female samples (Supplementary Fig. 2F). Marker genes for major cell types were
identified using FindAllMarkers (logfc.threshold = 0.25, min.pct = 0.1, test.use = “wilcox”,
min.cells.feature = 3, mi.cells.group = 3, pseudocount.use =1, and return.thresh = 0.01). The top 100
marker genes for each cluster (p < 0.05, ranked by log2fold change versus all other clusters) were
then analyzed for gene and disease ontology enrichments using compareCluster in clusterProfiler R
package (v. 4.12.6). Cranial neural crest genes were compiled based on markers identified by
regulatory network construction in human cultured CNCC and craniofacial tissue data65. Major cell
type labels were applied to each cluster.

RNA velocity analyses
RNA velocity analyses were performed on the main cell types to assess their transcriptional
dynamics. We used the “run” command within the velocyto python package137 (v. 0.17.17) with the
sorted bam files from each sample and the human gencode v32 gtf file as input. This resulted in the
computation of a loom file for each sample containing two count matrices: one containing unspliced
counts and another containing spliced counts. The loom files for each of the samples were merged
and subset to exclude the previously filtered out cells. We then used the scVelo package (v. 0.3.3) in
python to further filter, normalize and process the data. We used the filter_and _normalize function
with paramters min_shared_counts set to 20 and n_top_genes set to 2000. Next, we computed the
moments using moments function with n_pcs and n_neighbors set to 30. Then, the
recover_dynamics function was run with n_jobs set to 30. Using the CellRank package138 (v. 2.0.6) in
python, we set up the VelocityKernel containing the scVelo-computed velocities to compute a
transition matrix using the compute_transition_matrix function. To increase the robustness of this
calculation, we combined the VelocityKernel with the similarity-based ConnectivityKernel. To
accomplish this, we also used the compute_transition_matrix function with the ConnectivityKernel
and then combined the two kernels using the following function: combined_kernel = 0.8 *
VelocityKernel + 0.2 * ConnectivityKernel. We visualized the transition matrix on the UMAP of all
main cell types using the plot_projection function.

For mouse data sets, filtered barcode matrices from each mouse E10.5 to 12.5 samples generated by
CellRanger ARC were individually loaded with Read10X_h5 command in Seurat134 and merged with
E13.5 to E15.5 data from Pina et al 202323 (GSE205448). Calculation of percent mitochondrial reads
and filtering were performed with similar thresholds to human data above. Subsequent
harmonization, dimensionality reduction, and clustering were performed identically to those for
human data above. Less significant contamination of neuronal cell types was observed in mouse
data, which was identified and filtered as described for human. Identification of marker genes and
gene ontology enrichments, and CNCC modules scores were performed as above for the human data.
Major cell type labels were applied to each mouse cluster.

Marker gene comparisons across species
Lists of all marker genes for each of the eight main subtypes for each species (p < 0.05) were
compiled and orthology based on HGNC symbol =was obtained using the convert_orthologs
command in the orthogene R package139 (v. 1.10.0) with non121_strategy set to
“drop_both_species” and method set to “gprofiler”. Only genes that had one ortholog in each
species and a HGNC symbol were retained (n = 15,699). Significant overlaps between all orthologous
marker gene lists were determined using the newGOM command in GeneOverlap R package140 (v
1.40.0). Conserved and species-specific genes were determined based on HGNC symbol and the
intersection matrix obtained by getMatrix in GeneOverlap. Gene and disease ontology enrichments
were calculated using clusterProfiler.

Final Seurat objects were prepared for display in an interactive webapp using the ShinyCell R
Package141.

Subclustering of major cell types
For major cell types labelled as mesenchymal and ectodermal, further subclustering was first
performed on mouse data. For each major cell type, normalization, scaling with regressed cell cycle
impacts, harmonization, and subclusters were identified using the same procedure as described
above. Marker genes were identified for each cluster and functional enrichments were determined
using clusterProfiler. Annotations for each cluster were manually assigned based on those originally
described13-16. Mouse cell subtype assignments were also confirmed with ToppGene142 using the
toppFun function from the scToppR package143. Mouse main and subtype annotations were further
confirmed by projection on mouse E15.5 spatial transcriptomics data23 (GSE245469). Links between
snRNA and spatial data were determined using FindTransferAnchors and transferred using
TransferData in Seurat.

Following annotation of subtypes and for comparison with human data, an intermediate data set was
created where mouse genes were reduced and converted to those to those with one-to-one fig.
orthology with human genes using the orthogene R package (v. 1.10.0) as described above. The
intermediate dataset was used to transfer mouse subtype annotations to human subtypes by first
identifying shared features across clusters using FindTransferAnchors in Seurat with log
normalization and canonical correlation analysis (cca) or reciprocal PCA projection (rpca). Predicted
subtype labels were transferred to human subtypes using TransferData in Seurat and further
confirmed with a confusion matrix. As an orthogonal approach to label transferring in Seurat, we
used SAMap72 (v. 1.0.15), an algorithm implemented in python that uses self-assembling manifolds to
map homologous cell types with shared expression programs across species. Briefly, we used the
map_genes function within the SAMAp package, which performs a BLAST search between species,
with t1 and t2 both set to ‘nucl’, tr1 set to the gencode v32 human fasta file, tr2 set to the the
gencode v10 mouse fasta file, n1 set to ‘hu’ for human, and n2 set to ‘mo’ for mouse. This generates
a “maps” directory with the transcriptome mapping BLAST results. Next, we used the SAMAP, run,
and get_mapping_scores functions within the SAMap package to calculate alignment scores between
cell types and subtypes across species. To finalize cell subtype annotations, we used a combinatorial
approach considering the three annotation methods (label transfer by Seurat using cca, label transfer
by Seurat using rpca, and label mapping using SAMap). Final annotations were confirmed using
contingency matrices, significant marker gene overlap between human and mouse subtypes using
the GeneOverlap R package140, and functional enrichments using the compareCluster function in the
clusterProfiler R package144.

The putative CNCC and early progenitor clusters were annotated by generating functional and/or
regional module scores using marker genes from various publicly available databases and resources.
To label the subclusters of the putative CNCC population, we generated a Schwann cell module score
using the marker genes from Cao et al. 202024. Additionally, we used the gene markers from various
functionally specific, fate committed CNCCs in Soldatov et al. 201970 and identified one cluster of
mesenchymal committed CNCCs and one cluster of sensory committed CNCCs. The remaining
putative CNCC subclusters were labeled according to Carnegie stage-specific biases. We randomly
downsampled the number of cells to 200 cells per stage and quantified the percentage of cells in
each cluster; those that were biased towards CS12 and CS13 were labeled as early (early.CNCC) and
those that were biased towards CS13-C17 were labeled as intermediate (int.CNCC). When
subclustering the early progenitor population, we generated module scores for neural progenitors
using markers genes from the ScType database145, neural stem cells using marker genes from Wang
et al. 2024146, mesenchymal progenitors using markers from Wang et al. 2024146 and epithelial
progenitors using markers from the human CellMarker2.0 database147.

Marker genes for human subtypes were identified as performed for major cell types and functional
enrichments were characterized with compareCluster in clusterProfiler. Final Seurat main objects and
subtype objects were prepared for display in an interactive webapp using the ShinyCell R Package141.
Seurat objects were also converted to h5 files using the SaveH5Seurat and Convert functions in the
SeuratDisk R package (v. 0.0.0.9021). These files are compatible with python packages including
scanpy148 and anndata149 and were used for hosting at the Chan Zuckerberg Initiative CELLxGENE
Discover resource150.

Processing of spatial transcriptomics
Spatial transcriptomics data for two sagittal sections of a human CS13 embryo25 were retrieved from
https://heoa.shinyapps.io/code/. Slice 1 is more laterally localized while Slice 2 is more medially
localized. Thus, representative plots for ectoderm are shown in Slice 1 and for putative CNCC’s and
mesenchyme are shown in Slice 2.

Raw sequence count matrices were loaded using the Read10x command of Seurat134 and converted
to HDF5 format. These counts were then combined with spot coordinates and section images using
CreateSeuratObject. Data from both slices were merged and variable features were determined
using Seurat. The percentage of mitochondrial reads was determined for each cell and was used to
transform all data in the merged object using SCTransform151 from Seurat. Data was clustered using
UMAP and plotted which revealed a strong batch effect between the two spatial objects. Data was
further normalized using Harmony (v. 1.2.3), projected using UMAP, and clustered with a resolution
of 0.8. Marker genes for the 22 clusters were identified using FindAllMarkers in Seurat. The top 100
marker genes for each cluster (p < 0.05, ranked by log2fold change versus all other clusters) were
then analyzed for gene and disease ontology enrichments using compareCluster in clusterProfiler R
package (v. 4.12.6). Enrichments and spatial localization were compared to previous annotations by
Xu et al 202325 and labelled accordingly. The top 100 marker genes from each of the subclusters
identified in human craniofacial data were used to calculate module scores across the merged spatial
object and plotted using SpatialFeaturePlot in Seurat.

E9.5-16.5 spatial transcriptomics datasets were retrieved from the MOSTA resource71 and loaded
into Seurat as for human data above. For representative images in the supplement, we chose to use
the mouse E11.5 spatial transcriptomics data as it is most morphologically similar to CS13 human
embryos. For the E9.5-E16.5 combined stages object, we converted the mouse genes to those with
one-to-one orthology with human genes using the orthogene R package (v. 1.10.0). Following, we
used FindTransferAnchors in Seurat with log normalization and canonical correlation analysis (cca) to
transfer the human main cell type and subtype annotations to the mouse spatial transcriptomics
dataset using only the highly variable genes identified using FindVariableFeatures in Seurat.
Predicted labels were then transferred to mouse spatial transcriptomics data using TransferData in
Seurat. For each E11.5 section, module scores of top 100 marker genes for each main cluster or
subcluster were calculated. Gene spatial feature plots for selected genes and module scores were
then generated with Seurat.

Facial variation and congenital abnormality GWAS enrichments
We retrieved summary statistics for facial variation88,89 and all congenital abnormality GWAS
summary statistics from FinnGenn93. Raw summary stats were further processed and standardized
with hg38 coordinates with MungeSumstats R package152 (v. 1.10.1). Variants were mapped to genes
+/- 100kb using MAGMA153. Frontal facial measures and FinnGen data were processed based on
1000 genome European population while profile facial measures were processed with the 1000
genome Middle/South American population all obtained from the MAGMA website
(https://cncr.nl/research/magma/). We converted the Seurat snRNA-seq expression data to a
CellTypeData set with the Expression Weighted Celltype Enrichment (EWCE) R package154 (v. 1.12.0)
and then assessed each study trait for a linear positive correlation of cell type gene expression
specificity and gene-level genetic association using MAGMA.Celltyping155 (v. 2.0.14). Plots were
generated using tidyheatmaps 156 (v. 0.2.1) and ggplot2 (v. 3.5.1) in R.

Spatial GWAS enrichments
The human CS13 spatial transcriptomics Seurat object was subset to include only Slice 1 and
then converted to h5ad format using the Convert function from the SeuratDisk package in R.
GWAS summary statistics of facial variation and congenital abnormalities from FinnGen were
then mapped to the human CS13 spatial transcriptomics data using gsMap94. The
format_sumstats function was first used to convert the FinnGenn GWAS summary statistics to a
format compatible by gsMap. A sumstats configuration file was then generated with
information for all GWAS traits analyzed and then gsmap quick_run was run with default
parameters. The gsMap run_cauchy_combination and run_report functions were run for each
GWAS trait to visualize the results. Genetic spatial mapping and gene specificity score plots
were generated using Seurat’s SpatialFeaturePlot function in R.

Gene list enrichments per cell type
The CellTypeData-formatted human craniofacial snRNA-seq objects were generated using
generate_celltype_data in EWCE (v. 1.12.0). Mean and specificity metrics for several marker genes
(SOX10, TP63, and MSX1) were inspected across main cell types and subtypes using plot_ctd in
EWCE. Gene lists were compiled from multiple resources including gnomAD (v4.1), CleftGeneDB,
prioritized genes and black module from craniofacial WGCNA13, and genes affected by de novo
variation in orofacial clefting probands32,98. For Neanderthal introgressed regions and human
accelerated regions, coordinates were obtained from respective publications157,158 and assigned
single nearest gene using rGREAT (v. 2.6.0) 159 with “oneClosest” association rule. Each gene list was
then tested for linear association using bootstrap enrichment test in EWCE (reps = 10,000;
geneSizeControl = TRUE). As an additional negative control, we generated random gene lists based
on the mean number of all the lists tested here. Results from all gene lists were then merged,
corrected for the total number of gene lists and cell types tested using the Benjamini-Hochberg
approach, and then plotted using the pheatmap package160 (v. 1.0.12) in R.

Phenotype-cell type association tests
To map the relationships between cell types and phenotypes, we ran pairwise association tests
between all combinations of cell types in our snRNA-seq-derived CellTypeData and phenotypes
across the Human Phenotype Ontology (HPO; release v2025-08-11)91 using the run_phenomix
function from MSTExplorer (v. 1.0.9). In contrast to the gene list-based approaches (e.g. EWCE) this
function reframes the problem as a series of generalized linear regressions by leveraging continuous
scores that summarize the current strength of evidence for a causal relationship between each gene-
phenotype pair (using additional data from the Gene Curation Coalition)90,161. The continuous nature
of these data allows us to more accurately capture phenotype-cell type relationships, especially for
phenotypes with large gene lists where only some genes have strong causal evidence for the
phenotype. The gene signature vectors for each phenotype were previously merged and shared as a
single precomputed gene (5,180 unique gene symbols) x phenotype (11,606 unique HPO
phenotypes) association matrix. Next, a series of linear regressions tests were performed between
the gene specificity vectors of each cell subtype (n=89 vectors) and the gene association vectors of
each phenotype (n=11,606 vectors). Finally, multiple-testing correction was applied using Benjamini-
Hochberg False Discovery Rate106 (at FDR<5% significance).

For the purposes of summarization and visualization, the number of significantly associated
phenotypes per cell type were then computed within each major HPO branch (Fig. 7C). Here, we
define HPO branches as groups of related phenotypes that can be labeled according to their shared
ancestral term, e.g., ‘Abnormality of the immune system’. Next, we sought to determine whether
some cell types were disproportionately more often associated with phenotypes of a particular HPO
branch. To accomplish this, we performed a series of proportion tests comparing the proportion of
total phenotypes that a given cell type was significantly associated with within a target HPO branches
relative to all other HPO branches. In practice, we computed 2x2 contingency tables (number of
significant phenotype association vs. number of non-significant phenotype associations x target
branch vs. non-target branches) for each cell type within each HPO branch, which were then used as
inputs to the prop_test function within the rstatix R package162 (v. 0.7.2). This test appropriately
takes into account the different number of phenotypes across HPO branches. Only one-sided tests
were performed to test whether the target HPO branch was greater than all other (non-target)
branches (set with the alternative = "greater" parameter). All proportion tests were then corrected
for multiple testing at FDR<5%.

Orofacial Clefting de novo variant analysis
We used the R package DenovolyzeR104 (v. 0.2.0) to test enrichment of de novo variants (DNs) in
a dataset of OFC case-parent trios. Enrichment is calculated by comparing the expected number
of variants, as determined by mutation models described by Samocha, et al105, to the observed
number of variants in a given gene or group of genes using the ‘DenovolyzeByClass’ and
‘includeGenes’ functions. Using our dataset of 2031 DNs in 1171 genes identified in 1676 trios
with OFCs, we first compared this list of genes to those with calculated mutational rates in the
R package ‘DenovolyzeR’ using the ‘viewProbabilityTable()’ function. Broken down by subtype,
there were 1180 cleft lip with or without cleft palate (CL/P; 226 cleft lip (CL), 966 cleft lip and
palate (CLP)), and 484 cleft palate (CP) trios. We then tested enrichment of all OFC trios and by
subtype within the top 20% of genes by log2FC derived from single nucleus RNA sequencing of
human craniofacial tissue at CS20.

Statistics and Reproducibility
Detailed statistical analyses are described the Methods section above. Biological replicates for
the human embryonic face samples were collected in the following way: four at CS12, six at
CS13, three at CS14, three at CS16, three at CS17, and five at CS20 for single nucleus
multiomics. Biological replicates of pooled mouse embryonic face samples were collected in the
following way: three at E10.5, three at E11.5, and three at E12.5 for single nucleus multiomics.
Appropriate negative controls were used throughout the study for GWAS enrichment across
main cell types and subtypes. The investigators were not blinded to allocation during
experiments and outcome assessment. No statistical method was used to pre-determine
sample size. No data that passed quality control criteria for experiments were excluded from
the analyses. The experiments were not randomized. Unless otherwise stated, default
parameter settings were employed for any software tool that was used in the analyses.
Whenever a p-value is reported in the text, the statistical test is also indicated. All statistics
were calculated, and plots were generated using R v4.4.0.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript “Gene expression patterns of the developing human face at single cell resolution reveal cell type contributions to normal facial variation and disease risk” by Khouri-Farah et al. describes a new resource providing high-resolution single-cell transcriptome data for multiple stages of human embryonic development. Importantly, these stages cover critical morphogenetic processes relevant to human facial variation and craniofacial birth defects.\n\nThere is currently a sparsity of high-quality data from actual human tissues covering these processes. Considering differences in craniofacial development and regulatory sequences between human and mouse (a model for which various datasets are available), the data provided in this study will be a rich resource for craniofacial researchers interested in studying general developmental processes, as well as the impact of human genetic variation on these processes in the context of birth defects. The authors provide an easily usable interface for researchers to explore the gene expression profiles in specific cell-type clusters, which will facilitate the utilization of their data by others.\n\nWhile the data provided through this study are without doubt valuable and deserve to be adequately described in a publication, the manuscript itself has significant room for improvement. The text is much longer and detailed than needed and many of the figures are convoluted and contain unnecessary detail. Also, the rationale for some of the analyses performed in this study is not clearly described.\n\nOverall, the generation of human single-cell RNA sequencing data across multiple timepoints is an important addition to the craniofacial field, but the current manuscript would require major revisions to make this valuable resource easily accessible and provide the advance and insight that would make it a strong candidate for Nature Communications.\n\nMajor points\n1. All figures need major revisions so that they can be easily interpreted. All plots, text, legends, axis labels need to be readable. The authors should think carefully about what is the key takeaway from each figure and which plots are absolutely required to make that point. Keeping a consistent color theme throughout would also improve the accessibility of the figures.\n2. The main text is far too long and needs to be at least 50% reduced. The readability of this manuscript also needs improvement. There are many different analyses presented but the rationale for their inclusion and key conclusions for the reader are not always clear. I would recommend adding a Supplementary Information section where more details can be added to improve the structure of the main text and make it easier to follow.\n3. The human single-cell data from this study was produced using 10X Multiome but only the RNA-seq component of the data is presented here. The ATAC component of the data is not needed for the purposes of this manuscript (which focuses on transcriptomes) but it would be good to mention for the reader why it is not reported. Was it of insufficient quality to provide useful results? Will it be described in a complementary manuscript? Some other reason?\n4. It’s surprising that Sox10 does not come up in any of the human CNCC subtype clusters, considering that it is a key marker in the mouse CNCC and given its known role in cranial neural cell development and differentiation. Why do the authors think that is? Given the overlap in markers with neural cells, is there any concern that human nuclei of CNCC were accidentally filtered out when removing neuronal cell types? E.g., TUBB3 is expressed in premigratory cranial neural crest cells in birds and maintained in neural crest derived tissues/cell types. See: doi: 10.1016/j.gep.2019.119067.\n\nMinor comments\n1. The PDF files are very large (causing poor display/printing issues on some systems) because dot plots are pasted in vector format with large numbers of individual datapoints. It would be helpful to rasterize plots so that the Figures can render more quickly and make it easier to review.\n2. The sentence in the introduction about Neanderthal-introgressed sequences seems a little out of place. Perhaps move this up to the section about normal facial variation and give some clearer rationale for why this analysis was done/included.\n3. Typo on 105 - doublets.\n4. Line 113: Which clusters were problematic? It’s not easy to tell from S1b as the colors merge together.\n5. Line 117: Please specify that nuclei are filtered specifically based on TUBB3 and MAP2 expression and not the other genes shown in S1.\n6. Figure 2B: Perhaps only show top 5 gene markers per cluster so that the y-axis is more easily read.\n7. Figure 2C - it’s not possible to read the terms in the smaller font. Simplify to top 3-5 enriched terms per cluster that can all be shown and read.\n8. The Venn Diagram in FigureS4B may not be necessary. It's not easy to interpret and the upset plot clearly shows the intersections.\n9. Line 136-138: This statement is somewhat unclear as at this point in the manuscript no enrichments for craniofacial biology and developmental abnormalities has been shown. I think the conclusion of this paragraph can be kept as simply “Overall, these results suggest that our single-nucleus expressions data closely resembles the bulk gene expression data.”\n10. Reword sentence at 154-155. It’s not necessary to state “some of the timepoints”. I suggest the authors rephrase so that it’s clear that at CS12 all samples were derived from female embryos and at CS13 all samples were derived from male embryos, but that this did not introduce any significant bias.\n11. Typo Line 155 - main cluster -> main clusters.\n12. Figure S5A - y-axis text is overlapping and not readable.\n13. Line 258 sub-cluster cnl2 is not described - what is this cluster?\n14. Line 271: clusters cnl1, cln2 and cnl3 need descriptions.\n15. Figure 3C: Do the different colors represent the Carnegie stages in this plot? Can this plot be simplified to highlight the genes referred to in the text? Or perhaps a heatmap would be easier to interpret.\n16. Figure 3D: Challenging to read the y-axis.\n17. Figure 3E: I suggest changing the color scale for these plots so that cream or a light color represents lowly expressed genes and a darker color (e.g., red/blue) represents highly expressed genes\n18. Line 300: Muscle cells clusters have been observed in the single-cell data from https://doi.org/10.1038/s41467-024-46396-4 and https://doi.org/10.1038/s41467-023-41349-9\n19. Line 301: It seems unlikely that sampling differences and/or timepoints would be the reason these cells are not observed but I would imagine it’s due to the lack of specificity in resolved clusters. E.g., Sox10, one the main markers used by the authors for specifying CNCC cells, is also a specific marker for neural-crest derived Schwann cells. Is it also possible that CNCCs had not been resolved in preexisting single-cell datasets because they were instead labelled broadly as neural cell types?\n20. What is the rationale for selecting E10.5-E12.5 in mouse for resolving muscle cells?\n21. Figure 4 needs to be simplified. It has too many panels and it is difficult to understand what the key takeaway is meant to be. What are the absolute key plots to show? And can the rest be put in supplements? For example, Figure 4K and Figure 4N are not readable.\n22. Figure 4G: Please specify what the subtype labels represent.\n23. When comparing single-cell data between species the authors could consider joint-embedding with SAMap (https://doi.org/10.7554/eLife.66747) to identify homologous cell types.\n24. For the cross-species cell type comparisons, what gene list was used as background for the enrichment analyses? Is it all human/mouse genes or genes that have 1-to-1 orthologs between mouse and human?\n25. Figure 7: It is difficult to understand the main takeaway from these plots. The heatmap legend is not labelled and the values are not easily interpreted. Perhaps parring these heatmaps down to only those that are significant for the main figure could help. The x-axes are not readable so it’s impossible to know which sub-cell types have specific variant enrichments. In Figure 7C, do the different colors represent the different human stages?\n26. Line 621: What is this referring to? The HARs or overall for all gene lists?\n27. Line 713: From Supplementary Table 7, Sox10 does not appear in the CNCC subtypes for human data.", "reviewer_2": "The manuscript by Khouri-Farah et al describes a novel dataset derived from single-nucleus RNA-seq analysis of human embryonic face and head development across weeks 4-8 after conception. They define several broad classes of cell types, including epithelial and mesenchymal, and a population defined as neural crest, and delve deeper into the cell identities by sub-clustering. Further, they perform a similar analysis using approximately equivalent stages of mouse embryogenesis and compare and contrast the data obtained between the two species. The authors next used a published human spatial transcriptomics dataset as a means to validate several of their findings concerning their human assignments. The human single-nucleus RNA-seq data are also used in conjunction with several public databases to make predictions concerning disease associations, human facial shape variation, connections between human and Neanderthal facial morphology, and clefting genetics. Overall, the dataset is a valuable additional resource for people interested in human embryonic head development that adds to other bulk, single-cell, and spatial datasets focused on the entire embryo. However, there are several significant issues with the manuscript, particularly with the analysis and interpretation of the datasets, that considerably reduce its impact and utility.\n\nMajor issues\n1. The manuscript is too long, and there are many sections that are redundant or repetitive. The journal’s instructions for authors state that ideally the main text should be 5,000 words. This manuscript is well over 10,000 words, with the discussion alone about 3,000 words. Much of the data from the results section is reiterated in the discussion and there are paragraphs dealing with disease correlation analyses of red blood cells and immune cells that are not pertinent for an analysis focused on the craniofacial region. For example, the paragraph beginning on line 191 concerning ontologies associated with red blood cells including malaria or hemolytic anemia or the paragraph beginning on line 203 concerning immune cell related diseases. The length of the manuscript and its multiple individual correlations with other datasets means that the focus of the paper is lost, and the key points are difficult to extract from the study. Overall, the manuscript would be more impactful if it was considerably shortened and focused on key findings concerning the craniofacial region.\n2. A major flaw in the manuscript is a lack of clarity concerning the regions under analysis in the sn-RNA-seq study. For human, only the following statements are provided: Results “we profiled the entire craniofacial prominence” and Methods “Primary human craniofacial tissues were obtained”. For mouse: Results “we generated single-nucleus gene expression from mouse craniofacial tissue. These samples reflected the major morphological landmarks of the human tissue profiled allowing a more direct comparison of cell types”. Methods “primary mouse craniofacial tissues … were pooled”. Unfortunately, these extremely limited descriptions do not provide sufficient information to understand or review the manuscript with any degree of certainty. One illustration of how this lack of detail affects the manuscript is that – as illustrated in Figure 6 - the craniofacial human ectoderm samples contain tissues related to the optic and otic placodes as well as to the neck (thyroid). These findings indicate a much greater region of the head and neck were surveyed than just the craniofacial prominences shown in Figure 1 and it is imperative to know the morphological boundaries that were used for the dissections.\nAs another example, on line 116, describing clusters expressing TUBB3 and MAP2: “We reasoned these clusters were derived from developing brain tissue not directly part of the craniofacial structures”. In removing these “brain clusters” their sample size decreased from 86,359 nuclei to 42,131, a drop of over half. These data indicate that their dissections of craniofacial tissue were not precise but included a majority of cells from other regions of the head. On the other hand, throwing away cells with expression of these two markers likely means they are losing cells from critical craniofacial structures, particularly the trigeminal ganglia. The authors indicate that the mouse data were obtained and processed in a similar manner, but it would be valuable to know if it required excluding Tubb3 and Map2 cells, and if so what percentage of mouse cells were removed from the analysis to help determine equivalence of dissections.\nThese considerations are especially important when the authors discuss the idea that there are multiple marker genes that were only expressed in a single species (Results: Section titled “Species-specific differences in marker gene expression during human and mouse craniofacial development”). Examination of the data provided in Table 20, for example, states that genes including ALX3, DLX2, DLX3, DLX5, and DLX6 are only detected as expressed in the human facial mesenchyme and not in the mouse “Only called in a single species”. This is certainly not supported by multiple other studies, raising concerns about the scope of the dissections, the data processing, and attention to detail in data analysis.\n3. The question of anatomy and the regions dissected also has a direct bearing on the parts of the paper discussing the identification and analysis of clusters 1 and 6, called as mesenchyme and cranial neural crest cells, respectively. The mesenchyme of the head is derived from two major tissue-types, the mesoderm and the neural crest. Fate mapping studies have indicated that the majority of the mesenchyme in the craniofacial prominences is of neural crest origin, whereas mesodermal derived mesenchyme occupies areas of the head that are above and behind the developing eye. There is also a core of mesoderm tissue in the arches, including arch 1, that will form muscles of the head and neck. Therefore, it would be expected that a significant portion of the mesenchyme in cluster 1 would have a neural crest origin. However, the ratios of neural crest versus mesoderm derived mesenchyme will vary depending on the amount of tissue included from more posterior regions of the head. This consideration again makes it imperative to understand the boundaries of the dissections performed.\nEven if a large amount of mesodermally-derived mesenchyme was present in cluster 1, there should still be a significant contribution of neural crest to this cluster. In this respect, when the authors subcluster the mesenchyme, they identify maxilla, mandibular and nasal prominence cells indicating the presence of neural crest derivatives. Therefore, cluster 6 should not be addressed as the only cranial neural crest cell population in their datasets. In fact, examination of cluster 6 suggests it may represent a few specialized neural crest derivatives – particularly Schwann cells as well as a few melanocytes. For example, as shown in Figure 3C, PLP1 - a prominent component of myelin - is expressed in all subclusters of cluster 6. Similarly, Mpz (myelin protein z) is broadly expressed in most subclusters. Together, these observations suggest that the mesenchyme has been mis-classified as not containing neural crest, while cluster 6 has been annotated as the major cranial neural crest population, whereas it represents highly specific populations of neural crest derivatives.\n4. The authors define many subclusters in the mesenchyme, ectoderm, and neural crest but do not adequately describe how they were assigned to specific regions or structures. For example, for the ectoderm, TBX18 was assigned as a fusion zone marker, but it is not clear why, nor is it clear what fusion zone they are discussing anatomically. These issues of mapping and assignment are compounded by an absence of appropriate verification analyses. Here the authors rely on published data of spatial transcriptomics of two sagittal sections through a CS13 human embryo to overlay their expression data which are then presented in several figures throughout the manuscript. There are two problems with this approach. The first is that they are relying on just two sections and so it is not possible to determine the 3-dimensional expression patterns of any gene or cluster, meaning that considerable information is lost regarding assignment and anatomy. The second issue is that the resolution of spatial mapping was a 55µm spot, whereas an approximately 2µm spot would be required to accurately map a transcript to a single cell layer. Therefore, the study fails to link subcluster identity with specific anatomy at a meaningful resolution and in the absence of any other verification studies this significantly reduces the power of the dataset.\n\nMinor issues\n1. Line 462-463. “with LNP1 being more posterior and LNP1 being more anterior”. One of those should not be LNP1.\n2. Line 625. “We found no enrichment in any cell types…”. This line is duplicated in the next sentence.\n3. The authors refer to Li et al (reference 12) throughout the manuscript concerning mouse craniofacial data. However, this reference describes genome analysis of human brain data.\n4. I have a number of suggestions for the images of the human sections. First, as the manuscript is focused on the head and face, it would be more useful to show only these structures rather than the whole embryo to reflect the focus of the paper. Second, orientation should be indicated on the panels: left/right (i.e. for Fig 4C), anterior/posterior, rostral/caudal. Third, it would be useful to have one or more annotated images describing the pertinent structures present in the images.\n5. The figure legends for the Supplemental figures are missing. There was also no explanation of the information provided in the Supplemental tables.\n6. Multiple instances: “this data” should be “these data”.", "reviewer_3": "The paper by Khouri-Farah et al., “Gene expression patterns of the developing human face at single cell resolution reveal cell type contributions to normal facial variation and disease risk” reports a single-nucleus RNA sequencing (snRNA-seq) atlas of human craniofacial development from craniofacial tissues of 24 embryos that comprise six key time points spanning 4–8 post-conception weeks. The Authors describe seven major cell types, including muscle progenitors and cranial neural crest cells (CNCCs), as well as multiple subtypes of ectoderm and mesenchyme. Comparative analyses reveal substantial conservation of major cell types, alongside human biased differences in gene expression programs. The Authors also claim that CNCCs exhibit the lowest marker gene conservation, underscoring their evolutionary plasticity. Spatial transcriptomics is also employed to further localize cell populations, which should provide a detailed view of their anatomical context. Lastly, the Authors attempt to link these developmental processes to genetic variation, identifying cell type-specific enrichments for common variants associated with facial morphology and rare variants linked to orofacial clefts.\n\nUndoubtedly, this atlas has the potential to offer critical insights into the cellular and genetic mechanisms shaping the human face, highlighting conserved and distinctly human aspects of craniofacial biology and linking them to the developmental origins of craniofacial disorders and the genetic basis of facial variation. This work could certainly provide an important resource for exploring craniofacial biology, with implications for developmental genetics, evolutionary biology, and congenital anomalies. Therefore, the study could be of interest and could fill a poorly explored aspect of human craniofacial biology by reporting studies conducted on human embryos as well as mouse embryos. However, there are multiple serious concerns regarding the current manuscript that need to be fully addressed. I list below some of the most salient issues that need to be resolved in order to increase the rigor and accuracy of the described findings, as well as supporting many of the conclusions.\n\nSingle-nucleus RNA Sequencing (snRNA-seq) Data Generation and Processing (Human Samples)\n• It would be important to show a cartoon schematizing the regions of the human embryonic heads that were dissected. The Authors state that they processed the “developing craniofacial region” (quote) of human embryos from 4 to 8 weeks post conception. Lacking more detailed information, the reader is left guessing which parts of the developing craniofacial region were processed.\n• How did the Authors generate the initial clustering of the major cell populations? Which Seurat version did they use? Which other packages were employed? All relevant information is lacking.\n• More details are required in the Methods sections describing how the single-cell data were embedded and clustered, e.g. how many highly variable genes, PC dimensions, and nearest neighbors used to generate the UMAPS appearing in the main figures. Was any batch correction used? Were distinct embeddings calculated for the subset of neural crest clusters?\n• It would also be important to provide separate clusters for each gestational time analyzed. This would highlight differences that are to be expected in cluster transcriptomes across the different gestational days, which are currently not clear.\n• Why do the Authors use only TUBB3 and MAP2 as markers for neuronal cells that they want to remove from the query set? While primarily known as a neuronal marker, TUBB3 expression is also observed in non-neuronal cells and tissues, including endothelial cells, fibroblasts, spermatogenic cells, and stem cells, suggesting a broader role beyond just neurons. The choice of neural markers appears to be ill-informed and extremely reductive.\n• Regrettably, there appears to be lack of rigor in tissue collection. It is of concern that the Authors lose half of their cells through the first filtering step. According to their statements, they started from 86,359 individual nuclei, but ended up with 42,131 remaining nuclei for their analyses. They concluded: “We reasoned these clusters were derived from developing brain tissue not directly part of the craniofacial structures. We therefore excluded these nuclei from downstream analyses resulting in a total of 42,131 remaining nuclei with an average of 4095 nuclei from each sample”. In summary, either the contamination of facial tissues by brain tissue was massive, or the Authors also removed cells that are not neural cells and thus could be integral part of head developing structures.\n• The Authors state: “Early samples had consistently higher mitochondrial reads (Fig. S2d), potentially reflecting their higher dependence on mitochondrial output or an artifact related to lower cell numbers in the processing of each sample.”\nThis assertion is superficial and must be revisited. It is well known that high mitochondrial signatures are a strong indicator of apoptosis. The Authors have not unequivocally ruled out the presence of cell death in their samples.\n\nSingle-nucleus RNA Sequencing (snRNA-seq) Data Generation and Processing (Mouse Samples)\n• Similar concerns are present regarding the isolation of mouse tissues. There is a lack of adequate information regarding the collection of the mouse tissues. Were mouse tissues also fresh frozen + thawed (like human tissues) or were they processed directly following tissue collection? Differences in tissue processing could heavily influence the snRNA-seq results obtained. There are serious concerns that the “species-specific” signatures the Author report might be due at least in part to differences in processing of the collected tissues, and they may not necessarily underpin species-specific identities.\nThe authors state that: “Samples were mechanically broken into single-cell suspensions” and that “Craniofacial prominences were collected in a very similar fashion to human samples”. What do the Authors mean by “similar”? Were neural tissue also removed in mouse craniofacial tissue dissections as it was done for human dissections? Or was there absence of neural contamination? All of these details, which have critical bearing on the reported findings of species-specific differences, are not provided.\n• The Authors compare their datasets to those reported by the group of Trevor Williams (Li, H., Jones, K.L., Hooper, J.E. & Williams, T. The molecular anatomy of mammalian upper lip and primary palate fusion at single cell resolution. Development 146, 2019). But the latter datasets were generated only from murine midfaces and only at E11.5. This substantial difference in the datasets should at least be mentioned.\n\nNeural Crest Cell Cluster (Cluster 6)\n• What is cluster 6? It is not at all clear whether this is a biologically defined group of cells or a somewhat “mixed bag of cells” cluster. The Authors seem not to be completely sure themselves: “We then turned to the not yet concretely identified cluster 6 […]. We did not detect any significant pathway enrichments for this particular cluster (201). “The wide variety of biological functions and specific disease enrichments all suggested that this cluster might be enriched for neural crest cells. (208-210)”.\nThese descriptions are rather confusing.\n• The Authors state that CNCC markers are significantly higher in cluster 6.\n“When we examined a module of genes from regulatory networks recently identified in cultured human and chimpanzee cranial neural crest cells (CNCCs)69, we found significantly higher expression in cluster 6 (Fig. 2D).”\nWhich test or package or statistics did they use to reach this conclusion? Again, no relevant information is provided.\n• Do all developmental stages contribute equally to the population in cluster 6? This information is lacking and is very important.\n• The Authors state: “SOX9 and COL20A1 were more specifically expressed across the (CNCC) clusters, but again in non-overlapping patterns.” (261-262).\nThis sentence is very confusing. What do the Authors mean?\n• There appears to be lack of a definite proof that Cluster 6 is truly a Neural Crest Cell cluster.\n\nSpatial Transcriptomics\n• There appears to be huge differences in between the two replicates shown for spatial transcriptomics (Fig. 3 E,F,G). For certain markers and modules the two replicates shown in the figure appear to be completely different from one another regarding the respective gene expression spatial features. Clearly\nthese human samples are precious and hard to collect, making it difficult to generate multiple biological replicates.. However, it is very hard to draw firm conclusions from these spatial validations.\n\nEnrichment of Disease Genes and Variants Linked to Facial Variation\nThese analyses are hard to follow and very cryptic. Disease associations are correlative. Certain statements and conclusions do not appear to be based on rigorous supporting data. One example is provided below. The Authors state:\n“Interestingly, while CNCCs certainly give rise to many of the downstream cell types and tissues, we found relatively few shape associations for CNCC subtypes” (531-532)\nThis statement is somewhat concerning, given that a rich body of literature directly links differences of facial shapes to variants in CNCC genes. The following review article comprehensively summarizes most of the relevant literature showing that variants in CNCC genes and enhancers underpin differences in face variation: Naqvi S, Hoskens H, Wilke F, Weinberg SM, Shaffer JR, Walsh S, Shriver MD, Wysocka J, Claes P. Decoding the Human Face: Progress and Challenges in Understanding the Genetics of Craniofacial Morphology. Annu Rev Genomics Hum Genet. 2022 Aug 31;23:383-412.\nInterestingly, the Auhtors cite this review article in their “Introduction”.\n\nMinor points\n• Inconsistency in size of figure labelling (A,B,C etc.), legends in multiple figures (partially) overlap, image quality in multiple figures is poor.\n• Lettering in some figures is exceedingly small and/or warped and thus very hard to read, even when zooming in onto the image or cartoon.\n• The Authors need to revise the nomenclature for genes. Mouse genes should not be capitalized.\n• Some of the references cited are not correct. For example, Ref 12 (line 146) is not correct.\n\nStyle and Grammar (only a few issues of many are reported below):\n• 105: “after filtering for doubles and quality of per nucleus data“. Confusing sentence.\n• 106/107: “Samples obtained from human embryos are more difficult to control what is collected”.\nSentence needs rewriting.\n• 256: “TFAP2A was expression was observed in all clusters”. Sentence needs rewriting.\n• 625-627: “We found no enrichment in any cell types for genes affected by de novo synonymous variants. We found no enrichment in any cell types for genes affected by de novo synonymous variants.”\nThe sentence is repeated twice.\n• 758: space and period missing: “embryo25By“\n• 815 f. Clip -> lip: “Failure of this region to fuse in humans has been suggested to cause cleft clip that could also involve the\nnostril region and primary palate”\n• 833: dominant or major: “Mesenchyme is by far the dominant major cell type”\n• 993: describeD: “as describe for human”\n• 1019: remove word duplication: “mouse genes were reduced and converted to those to those with one”\n• 1058: space and period missing: “MAGMA89Frontal”\nLastly, this reviewer suggests careful editing of the text – in its present form, the manuscript appears to have been rushed and, regrettably, it is overall sloppy.", "reviewer_4": "The authors present a comprehensive single-nucleus RNA sequencing (snRNA-seq) atlas of human craniofacial development, derived from 24 embryonic samples spanning six key time points during the embryonic period. Additionally, they integrate spatial transcriptomics data to further localize cell populations. Our findings shed light on the developmental origins of craniofacial disorders, the genetic basis of facial variation, and the evolutionary legacy of ancient hominins. Overall, this work provides a valuable new atlas for craniofacial biology, supported by comprehensive experiments.\nHowever, there are still areas that require improvement. Specifically, some paragraphs should be refined or removed to enhance readability. At present, several sections in the Results are overly redundant, making it difficult to focus on the key findings. I encourage the authors to revise the “Results” section for improved clarity and readability.\n\nMy comments are as follows:\n1. Beyond providing a more extensive atlas, have the authors identified any conclusions that Yankee, T.N. et al. (Ref 13) did not?\n2. The authors state, “Closer inspection of the genes strongly expressed in these clusters revealed many canonical neuronal genes, such as TUBB3 and MAP2.” How did they quantify \"strong\" expression? Furthermore, they excluded nuclei derived from developing brain tissue. How can they ensure that all unrelated clusters from craniofacial structures were successfully removed?\n3. Line 172: Why did the authors choose 100 marker genes per cluster? In contrast, in line 156, they selected only 10 marker genes to identify unknown clusters. What is the rationale behind these differing marker gene selection criteria?\n4. Lines 172–216: The authors should refine these paragraphs to highlight the main findings rather than listing excessive details.\n5. Line 973: It is unclear why the authors selected multiple resolution values ranging from 0 to 1 to identify clusters.\n6. Line 1034: The authors should clarify how they preprocessed the spatial transcriptomics data quantitatively and in detail. Did they filter out low-quality spots?\n7. Figures 2, 3, 4, 5, and 6 share a common issue: they contain excessive information, making them overwhelming. The authors should consider moving less critical subfigures to the supplementary materials. The main figures should emphasize distinct and important messages for readers.", "reviewer_5": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_6": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We thank the reviewers for their thoughtful commentary and feedback on our\nmanuscript and have diligently considered these suggestions to improve our work.\nSpecifically, we have considerably reduced the length of the text to refine the scope of\nthe manuscript, improve readability, and strengthen the rationale. We have also\nincluded more detailed explanations of the anatomy and region dissected in both the\nhuman and mouse samples that were used for snRNA-seq data collection. Furthermore,\nwe have reanalyzed the snRNA-seq data using previously presented as well as newly\nsuggested methods. To maximize the number of nuclei for downstream analysis, we\nperformed filtering on raw filtered barcodes from the single nucleus multiome\nexperiment, resulting in an increase of over 44,000 nuclei (130,427 in new analysis vs.\n86,359 in previous analysis). We also removed neuronal and radial glia populations\nusing a new systematic method of generating module scores, resulting in the resolution\nof roughly 54,000 additional nuclei (95,925 in new analysis vs. 42,131 in previous\nanalysis). When subclustering the putative CNCC’s, we labeled plausible SOX10+\nSchwann using a similar module scoring methodology as described for neuronal and\nradial glia cell types. We have also streamlined the comparisons with mouse\ncraniofacial samples and focused on conservation of cell subtype annotations on a large\natlas of mouse spatial transcriptomic data. Based on these new analyses, we have\nupdated all main and supplemental figures to reflect these changes in addition to\nsimplifying main figures to emphasize critical points and improve interpretability for the\nreader. We believe these substantial additions have improved the manuscript and\naddressed the reviewers’ concerns. We have also addressed each of the reviewers’\nspecific comments below (in blue). We also provide a version of the manuscript that\nhighlights all tracked changes in this revision compared to the original submission.\nLastly, we have substantially improved our online application that accompanies the\npaper to include the new clustering as well as spatial transcriptomics data. We hope the\nreviewers and editor find these revisions to be more than sufficient for publication."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have addressed my previous comments and present a substantially improved manuscript. I am supportive of publication. However, I encourage the authors to further tighten the text. The manuscript still appears to exceed the recommended word limit for this format by more than 50% (approximately 7,600 words for a 5,000-word limit). Condensing it to around 5,000 words would improve its focus, readability, and overall accessibility.", "reviewer_2": "The revised manuscript by Khouri-Farah et al describing a novel dataset derived from single-nucleus RNA-seq analysis of human embryonic face and head development across weeks 4-8 after conception is significantly improved from the previous version. However, several major concerns still remain that should be rectified.\n1. Although it has been shortened, the manuscript is still very long at about 7500 words through the Discussion section. The main value of the manuscript is that it describes the overall dataset that can be used as a resource for the community. However, the author’s attempt to cover multiple topics – the human data, the mouse data, the human versus mouse data, human disease associations, human facial shape variation, and Neanderthal anatomy – with the data they have obtained. This makes the analysis very broad and it reads as one list after another. Moreover, verification of the data can only be superficial given the limited access to human embryo tissue samples. Therefore, these limits in resolution and in mapping of human data onto human anatomy reduces the utility of using such datasets to study shape variation in a robust manner. The manuscript would be improved by moving some of these sections on anatomy to a different manuscript to shorten and focus the current version.\n2. The text of the abstract, introduction and results is still misleading in regards to the tissues used for the analysis, as for example, when the author’s state that they have “constructed …. an atlas of human craniofacial development by collecting craniofacial tissues from 24 individual human embryos” (Section Main, line 61). Although they provide caveats at the very end of the discussion, it would be more accurate to state throughout the manuscript that they have “collected head and neck tissue containing the craniofacial complex and used bioinformatic methods to remove presumptive neural tissue from the analysis.” This would provide readers with a more accurate picture of the data used throughout the manuscript.\n3. The mesenchyme, cluster 0, is the largest cluster they identify in the UMAP analysis. In my previous review I had highlighted that the mesenchyme of the head is derived from two major tissue-types, the mesoderm and the neural crest. Therefore, it would be expected that a very large fraction of cluster 0 actually represents neural crest derivatives. For example, one of the markers they use for cluster 0 is ALX1. Expression of this gene is known to be strongly associated with the craniofacial neural crest in human, mouse, and zebrafish (e.g. Pini et al. EMBO Mol Med 2020 “ALX1-related frontonasal dysplasia results from defective neural crest cell development and migration” and Iyyanar et al Front Cell Dev Bio 2022. “Alx1 Deficient Mice Recapitulate Craniofacial Phenotype and Reveal Developmental Basis of ALX1-Related Frontonasal Dysplasia”.\nMy interpretation of the author’s presentation is that only clusters 5 and 7 are neural crest. Perhaps the author’s postulate that clusters 5 and 7 are migratory neural crest cells and can be viewed differently than other neural crest derivatives. However, as cluster 5 is associated with cell types differentiating towards Schwann cells and melanocytes, they would appear to be no different that neural crest cells in cluster 0 that are differentiating towards skeletal elements. Therefore, it is imperative that the authors directly discuss the contribution of neural crest to cluster 0. I would suggest that altering the following statement on line 100 would be one of the simplest ways to accomplish this - although there are others:\n“ALX1 was most strongly expressed in cluster 0 (mesenchyme),”\nchanged to\n“ALX1 was most strongly expressed in cluster 0 (mesenchyme), a tissue -type which has a dual origin from the mesoderm and cranial neural crest during facial development”\nI would then also recommend changing:\n“We inferred RNA velocity and transition states to identify trajectories in the data across the newly identified main cell types. We found that the trajectories emanate from the SOX2+/SOX10- population, confirming that these cells likely represent an early progenitor population (Fig. 2c).”\nTo:\n“We inferred RNA velocity and transition states to identify trajectories in the data across the newly identified main cell types. We found that the trajectories emanate from the SOX2+/SOX10- population, confirming that these cells likely represent an early progenitor population for other neural crest populations in clusters 0 and 5 (Fig. 2c).”\nDefining the contribution of the neural crest to the craniofacial mesenchyme also has implications for how the data relating to critical tissues involved in craniofacial disease and shape variation are considered later in the results section.\n4. Spatial transcriptomics on two human embryo slide sections is still utilized for verification analysis in several main figures as well as in multiple Supplementary figures. As stated in the previous review the resolution of spatial mapping in the current study was a 55µm spot, whereas an approximately 2µm spot would be required to accurately map a transcript to a single cell layer. Given this issue, it is very difficult to interpret the data shown. I would recommend that arrows are used highlight and label the regions we should be focusing on, or else many of the panels should be dropped if they do not contribute to the analysis so that a reader can focus on the salient points of the study.", "reviewer_3": "The paper by Nagham Khouri-Farah and Alexandra Manche “Gene expression patterns of the developing human face at single-cell resolution reveal cell-type contributions to normal facial variation and disease risk” has been substantially improved after revision and is now stronger, shorter, and clearer.\nThe following issues still need to be addressed:\n1) The authors state that: “Samples were mechanically broken into single-cell suspensions” and that “Craniofacial prominences were collected in a very similar fashion to human samples”.\nWhat do the Authors mean by “similar”? Was the neural tissue removed when conducting the mouse craniofacial tissue dissections, as it was done for human dissections? Was there absence of neural contamination? All of these details, which have critical bearing on the reported findings of species-specific differences, have not yet been clearly provided.\nThe Authors state that they have included additional information about the collection and processing of mouse craniofacial tissue in the Methods Section of the revised version of the manuscript:\n“less significant contamination of neuronal cell types was observed in mouse data, which was identified and filtered as described for human.”\nThe point has not been fully addressed; it would help to clearly explain how the dissections were carried out. In the Methods Section it should be precisely specified which mouse craniofacial prominences were dissected and analyzed for the transcriptomic analysis. The text states “craniofacial prominences”. Was the whole mouse embryonic head/face homogenized, or which facial prominences were dissected?\n2) The Authors compare their datasets to those reported by the group of Trevor Williams (Li, H., Jones, K.L., Hooper, J.E. & Williams, T. The molecular anatomy of mammalian upper lip and primary palate fusion at single cell resolution. Development 146, 2019).\nBut the latter datasets were generated only from murine midfacial prominences and only at E11.5.\nThis substantial difference in the datasets should at least be mentioned.\nThe Authors state: “We have now included this additional detail in the Results section of the revised version of the manuscript.”\nIt does not appear that this information was included or cited in the revised version of the document. It should be made clear where this information was added within the Results section, as this reviewer still cannot find it.\n3) The revised manuscript still needs additional clean-up. For instance, capital letters are used to label the figure panels, but then lowercase letters are used in the text to refer to the figure panels. Also, format should be consistent throughout the text.\nFor instance, in the Figure Legends either the Authors add a dot or do not add a dot after the parenthesis when they list the figure panels. To give an example: some figure legends appear like this:\nA). B). C) D). E)\nIf the Authors cannot spend the time, or do not care to fix these details, they should ask for the help of a capable editor. Otherwise, the paper will still give the overall impression of sloppiness.", "reviewer_4": "The authors have resolved all of my questions. I recommend accepting this paper.", "reviewer_5": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_6": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We thank the reviewers for their time and thoughtful comments in evaluating our\nrevision. We appreciate they recognized our significant efforts to streamline the text and\nprovide additional insightful analyses of this substantial craniofacial resource. Indeed 4\nof the six reviewers had no additional comments or suggestions and recommended\npublication. However, there were still a few areas of contention related to source of\ntissue and identification of specific cell-types/subtypes based on single genes. On this\nfirst point we have not tried to oversell what tissues we have obtained and make claims\nabout precision of tissue collection that we cannot defend. This is in line with many\nprevious papers’ description of the source of human and mouse craniofacial samples,\nincluding a recent paper in this journal (Rajderkar et al 2024). We instead leveraged\nunbiased, systematic filtering and identification of clusters to present to the reader what\ncell types we obtained. That being said we have added a simple statement to further\ndescribe what portions of the head we tried to exclude in the Methods section.\nSpecifically we have added the statement “Primary human craniofacial prominences\nexcluding the eyes and more anterior tissue from CS12, CS13, CS14, CS16, CS17 and\nCS20, each stage represented by a minimum of 3 replicates, were obtained from\nHDBR.” Furthermore while we have make our best efforts to collect comparable\nsamples from both human and mouse there is no way to make them identical. We again\nerr on the side of caution and contend that samples are very similar but not identical.\nOn the second point, we apologize if there was confusion about the use of published\nmarker genes described by Li et al 2019. Since this paper was one of the first to use\nsingle cell gene expression on craniofacial samples and has been frequently cited by\nother papers in the field, we used those markers in a qualitative way in our initial\nexploratory analysis of the craniofacial clusters. In the example of ALX1, our analysis\nfound it had biased expression of this gene in the large cluster initially labelled as\nmesenchyme but was far from being a top marker (number 523 on our list). Instead our\nresults highlight many other mesenchymal associated genes which we provided in\nFigure 2B and Supplemental Table 1 of the previous revision. To reduce the possibility\nof confusion we have altered Figure 2A heading to “Li et al. 2019 Cell-Type Markers” to\nhighlight how they were identified and used."}, {"round": 3, "reviewer_comments": {"reviewer_2": "The second revision of the manuscript by Khouri-Farah et al has made a few modest changes compared to the previous version. Most of the other comments from the three remaining viewers who had recommendations for alterations have instead only been considered in the rebuttal letter. Because of this, I believe that a review of the new manuscript is not warranted. Instead, I would suggest 3 things moving forward.\n1.\nThe paper is too long. I find it hard to read and pull out the really important details.\nThis is also a point raised by reviewer 1 again:\n“The manuscript still appears to exceed the recommended word limit for this format by more than 50% (approximately 7,600 words for a 5,000- word limit). Condensing it to around 5,000 words would improve its focus, readability, and overall accessibility.”", "reviewer_3": "The Authors have tried to improve the manuscript to the best of their abilities.\nThis Reviewer understands that working with human embryos is challenging and that certain issues (e.g. use of perfectly reproducible dissections of intact tissues) cannot be completely resolved.\nThis Reviewer also counts on the Editors of Nature Communications to fix some of the typos, awkward sentences, and formatting issues that are still present in the revised text and figures - and to polish the manuscript."}, "author_response": ""}]

Decision letter