Browse the corpus

Spatial multiomic landscape of the human placenta at molecular resolution. Successful pregnancy relies directly on the placenta's complex, dynamic, gene-regulatory networks. Disruption of this vast collection of intercellular and intracellular programs leads to pregnancy complications and developmental defects. In the present study, we generated a comprehensive, spatially resolved, multimodal cell census elucidating the molecular architecture of the first trimester human placenta. We utilized paired single-nucleus (sn)ATAC (assay for transposase accessible chromatin) sequencing and RNA sequencing (RNA-seq), spatial snATAC-seq and RNA-seq, and in situ sequencing and hybridization mapping of transcriptomes at molecular resolution to spatially reconstruct the joint epigenomic and transcriptomic regulatory landscape. Paired analyses unraveled intricate tumor-like gene expression and transcription factor motif programs potentially sustaining the placenta in a hostile uterine environment; further investigation of gene-linked cis-regulatory elements revealed heightened regulatory complexity that may govern trophoblast differentiation and placental disease risk. Complementary spatial mapping techniques decoded these programs within the placental villous core and extravillous trophoblast cell column architecture while simultaneously revealing niche-establishing transcriptional elements and cell-cell communication. Finally, we computationally imputed genome-wide, multiomic single-cell profiles and spatially characterized the placental chromatin accessibility landscape. This spatially resolved, single-cell multiomic framework of the first trimester human placenta serves as a blueprint for future studies on early placental development and pregnancy.

Successful gestation and the healthy development of the human embryo directly depend on the placenta, the first fetal organ1. Complications in placenta development disturb maternal and fetal health, leading to diseases such as gestational hypertension and diabetes, preterm birth and fetal growth restriction, pre-eclampsia and eclampsia, and spontaneous pregnancy loss2. Placental cells form the multifarious maternal–fetal interface (MFI), ensuring local and systemic immunological and vascular adaptation of the mother, vital for the adequate nutrition and well-being of the rapidly growing fetus3. Correct placentation and communication between maternal cells and trophoblasts (TBs) are fundamental to shaping a pregnancy-sustaining environment; it is interesting that common features are increasingly being discovered between tumor and placental cells4. Pioneering studies performed on human placenta have used genomics tools such as single-cell RNA-seq (scRNA-seq), scATAC-seq, low-resolution, non-single-cell spatial transcriptomics and targeted spatial antibody staining to characterize the heterogeneous landscape of the MFI and define cell type-specific patterns of gene expression and chromatin accessibility5–7. However, much work remains to better understand this complex interface.

TAC-seq, low-resolution, non-single-cell spatial transcriptomics and targeted spatial antibody staining to characterize the heterogeneous landscape of the MFI and define cell type-specific patterns of gene expression and chromatin accessibility5–7. However, much work remains to better understand this complex interface. In the present study, we applied joint snRNA-seq and snATAC-seq along with three complementary, single-cell spatial multiomics methods—Slide-tags8, STARmap in situ sequencing (ISS)9 and STARmap in situ hybridization (ISH)9—on the same first trimester placentas to reconstruct a comprehensive, spatially resolved, single-cell, multiomic cell census and molecular architecture of the early human placenta. Using our spatial multiomic data, we identified a new TB marker, FOXP1, and proposed yet unknown genes potentially regulating invasiveness and immunomodulation. We investigated and spatially explored transcriptional regulators potentially vital to placental function and decoded location-dependent gene expression and motif enrichment gradients throughout the placenta. We also described domains of regulatory chromatin (DORCs) comprising interwoven regulatory networks linking cis-regulatory elements (CREs) to gene expression. These CRE–gene links identified two origin roots for TB differentiation along with drivers of lineage determination and suggested a genome-wide association study (GWAS) hit for excessive vomiting in extravillous trophoblasts (EVTs). We next mapped cell types and their interactions, proposing an EVT-forming TP63-MYCNUT–MYCN axis and a human endogenous retrovirus (ERVH48-1) that paradoxically prevents syncytiotrophoblast (STB) formation. Finally, we computationally imputed a fully resolved, multiome-wide, spatial dataset mapping single-cell transcriptomes and chromatin accessibility, allowing us to spatially explore peak–gene linkages in the placenta. Taken together, we describe a comprehensive, spatially resolved interplay between intrinsic gene-regulatory networks and extrinsic intercellular communication while simultaneously providing a spatial, single-cell, multiomic architecture of the first trimester human placenta at molecular resolution.

To map gene expression and chromatin accessibility from the same cells in space, we first isolated human placenta nuclei from eight donors at early gestation (weeks 6–11) and cut tissue sections from these placentas (Fig. 1a and Supplementary Table 1). We next performed paired snRNA-seq and snATAC-seq on dissociated single nuclei to resolve cellular profiles across the placenta (Fig. 1b–d and Supplementary Table 1). Subsequently, we used three complementary spatial profiling technologies. We applied Slide-tags8, a spatial single-cell multiomic technique, to spatially contextualize our multiomic single-cell measurements via genome-wide chromatin accessibility and transcriptome landscapes from the same cells (Fig. 1e, Supplementary Table 1 and Methods). To address its sparse spatial sampling framework, we performed two additional complementary spatial mapping technologies: an image-based in situ RNA-seq method (STARmap-ISS)9 to map cell types in space using 1,001 highly variable genes selected from joint snRNA-seq and snATAC-seq analysis, followed by in situ RNA hybridization (STARmap-ISH)9 with 48 manually curated landmark genes to capture larger cross-sections of the placenta (Fig. 1f,g, Supplementary Table 1 and Methods).Fig. 1A single-cell transcriptomic and epigenomic reconstruction of the early human placenta.a, Placental structure and schematic of experimental design. First trimester placental villi are composed of an epithelial TB bilayer surrounding a villous core, the latter containing mesenchymal cells (FIBs), endo cells and fetal macrophages (HBCs). The vCTBs fuse to form an overlying hormone-producing STB layer in contact with maternal blood, representing both the placental barrier and the transport unit between the fetal and maternal systems. At villous tips, vCTBs proliferate and differentiate into EVTs that invade maternal layers to locally adapt the maternal immune system and transform uterine spiral arteries into large vessels of low resistance. Placental single cells were isolated for paired snRNA-seq–snATAC-seq analysis.

he fetal and maternal systems. At villous tips, vCTBs proliferate and differentiate into EVTs that invade maternal layers to locally adapt the maternal immune system and transform uterine spiral arteries into large vessels of low resistance. Placental single cells were isolated for paired snRNA-seq–snATAC-seq analysis. Sections of donor-matched frozen placental tissue were utilized for spatial snRNA-seq–snATAC-seq (Slide-tags), STARmap-ISS and STARmap-ISH analysis. Donor information, QC metrics and cell counts can be found in Supplementary Table 1. b, UMAP plots of snATAC-seq- and snRNA-seq-based clustering. Cells are colored according to cluster identification. Cluster composition can be found in Supplementary Table 1. c, UMAP representation of clusters of the placenta. Cells are colored according to cluster identification. d, UMAP plot showing all sequenced single-cell samples. Cells are colored according to sample. e, Slide-tags spatial visualization of identified clusters across samples W8-2, W9 and W11 with adjacent hematoxylin and eosin stains of the tissue section used in the experiment. Spatial multiomic profiles were obtained from 1,923 single cells, with 1,662 passing stricter joint ATAC and RNA filters. Clusters were manually annotated and revealed small subsets of cells not found in the main multiomic dataset. Cluster composition can be found in Supplementary Table 1. f, Multiomics-derived cell types identified in space with STARmap-ISS via Seurat integration on samples W8-2, W9, W11 and W7-1. Cluster composition can be found in Supplementary Table 1. g, Whole-section STARmap-ISH of 48 landmark genes, divided into 6 groups (n = 3). All groups were hybridized to each sample; here, we show representative images of each group for samples W7-2 (groups 1 and 3), W8-2 (groups 2 and 4) and W11 (groups 5 and 6). EBs, erythroblasts; EVT-prog, EVT-progenitor; myel., myeloid; STB-prog, STB-progenitor; unk., unknown.

ark genes, divided into 6 groups (n = 3). All groups were hybridized to each sample; here, we show representative images of each group for samples W7-2 (groups 1 and 3), W8-2 (groups 2 and 4) and W11 (groups 5 and 6). EBs, erythroblasts; EVT-prog, EVT-progenitor; myel., myeloid; STB-prog, STB-progenitor; unk., unknown. a, Placental structure and schematic of experimental design. First trimester placental villi are composed of an epithelial TB bilayer surrounding a villous core, the latter containing mesenchymal cells (FIBs), endo cells and fetal macrophages (HBCs). The vCTBs fuse to form an overlying hormone-producing STB layer in contact with maternal blood, representing both the placental barrier and the transport unit between the fetal and maternal systems. At villous tips, vCTBs proliferate and differentiate into EVTs that invade maternal layers to locally adapt the maternal immune system and transform uterine spiral arteries into large vessels of low resistance. Placental single cells were isolated for paired snRNA-seq–snATAC-seq analysis. Sections of donor-matched frozen placental tissue were utilized for spatial snRNA-seq–snATAC-seq (Slide-tags), STARmap-ISS and STARmap-ISH analysis. Donor information, QC metrics and cell counts can be found in Supplementary Table 1. b, UMAP plots of snATAC-seq- and snRNA-seq-based clustering. Cells are colored according to cluster identification. Cluster composition can be found in Supplementary Table 1. c, UMAP representation of clusters of the placenta. Cells are colored according to cluster identification. d, UMAP plot showing all sequenced single-cell samples. Cells are colored according to sample. e, Slide-tags spatial visualization of identified clusters across samples W8-2, W9 and W11 with adjacent hematoxylin and eosin stains of the tissue section used in the experiment. Spatial multiomic profiles were obtained from 1,923 single cells, with 1,662 passing stricter joint ATAC and RNA filters. Clusters were manually annotated and revealed small subsets of cells not found in the main multiomic dataset. Cluster composition can be found in Supplementary Table 1. f, Multiomics-derived cell types identified in space with STARmap-ISS via Seurat integration on samples W8-2, W9, W11 and W7-1. Cluster composition can be found in Supplementary Table 1. g, Whole-section STARmap-ISH of 48 landmark genes, divided into 6 groups (n = 3).

composition can be found in Supplementary Table 1. f, Multiomics-derived cell types identified in space with STARmap-ISS via Seurat integration on samples W8-2, W9, W11 and W7-1. Cluster composition can be found in Supplementary Table 1. g, Whole-section STARmap-ISH of 48 landmark genes, divided into 6 groups (n = 3). All groups were hybridized to each sample; here, we show representative images of each group for samples W7-2 (groups 1 and 3), W8-2 (groups 2 and 4) and W11 (groups 5 and 6). EBs, erythroblasts; EVT-prog, EVT-progenitor; myel., myeloid; STB-prog, STB-progenitor; unk., unknown.

composition can be found in Supplementary Table 1. f, Multiomics-derived cell types identified in space with STARmap-ISS via Seurat integration on samples W8-2, W9, W11 and W7-1. Cluster composition can be found in Supplementary Table 1. g, Whole-section STARmap-ISH of 48 landmark genes, divided into 6 groups (n = 3). All groups were hybridized to each sample; here, we show representative images of each group for samples W7-2 (groups 1 and 3), W8-2 (groups 2 and 4) and W11 (groups 5 and 6). EBs, erythroblasts; EVT-prog, EVT-progenitor; myel., myeloid; STB-prog, STB-progenitor; unk., unknown. Overall, we generated matched single-nucleus transcriptomes and single-nucleus epigenomes from 36,456 cells that passed stringent quality control (QC) criteria and filtering, with 98.6% fetal and 1.4% maternal cells (Extended Data Fig. 1a–c, Supplementary Fig. 1 and Methods). We performed unsupervised Louvain clustering and dimensionality reduction visualization on a concatenated cell-by-gene + peak matrix to identify 17 major clusters in total, which we manually annotated using canonical marker gene expression (Fig. 1c and Extended Data Fig. 1d–f). We identified all major cell types, including villous cytotrophoblasts (vCTB1–3, vCTBp), EVTs, STBs, endothelial cells (endos), myeloid cells, fetal fibroblasts (FIB1–2) and maternal decidual fibroblasts (mat. FIBs) (Fig. 1c). To further resolve clusters, we manually subclustered myeloid cells and EVTs (Extended Data Fig. 1g–j). Unknown 1 and myeloid unknown were almost exclusively from the same donor and exhibited FIB2 and Hofbauer cell (HBC) phenotypes, respectively; they probably correspond to these cell types but were differentiated by lower cell quality and viability (Supplementary Fig. 1). They may also comprise unknown cell populations. Unknown 2 expressed select marker genes associated with both FIBs and TBs and could represent either a doublet cluster or an unknown cell state (Extended Data Fig. 1d,e).

these cell types but were differentiated by lower cell quality and viability (Supplementary Fig. 1). They may also comprise unknown cell populations. Unknown 2 expressed select marker genes associated with both FIBs and TBs and could represent either a doublet cluster or an unknown cell state (Extended Data Fig. 1d,e). The vCTBs were separated into stem/progenitor (vCTB1,3), proliferative (vCTBp) and predominant (vCTB2) phenotypes. EVTs were divided into three phenotypes: progenitor (EVT1), cell column (EVT2) and mature secretory distal (EVT3)10. EVT2 expressed LY6E, a tumor-associated gene that induces embryonic lethality in mice owing to compromised fetal–maternal vascularization when knocked out11 (Extended Data Fig. 1k). PDGFRB+ FIB1 represented a perivascular subcluster, whereas the common PDGFRA+ FIB2 shared various gene expression patterns with cancer-associated fibroblasts including CXCL14, DCN and SPON1 (ref. 12) (Extended Data Fig. 1e). In EVTs, we identified multiple unique genes promoting cancer cell invasiveness such as QSOX1 and propose a distinct set of genes orchestrating the maternal immune cell landscape including FGFR1 (refs. 13,14) (Extended Data Fig. 1l). In STBs, BRAF and TBX3 may contribute to immune tolerance by recruiting myeloid-derived suppressor cells and forming an immunosuppressive environment15 (Extended Data Fig. 1l). Beyond these genes, we identified a new list of cluster-specific marker genes with yet unknown roles in placentation (Extended Data Fig. 1m and Supplementary Table 2).

ay contribute to immune tolerance by recruiting myeloid-derived suppressor cells and forming an immunosuppressive environment15 (Extended Data Fig. 1l). Beyond these genes, we identified a new list of cluster-specific marker genes with yet unknown roles in placentation (Extended Data Fig. 1m and Supplementary Table 2). Finally, we identified FOXP1 as a potential marker for vCTBs; it has previously been implicated in embryonic stem cell (ESC) lineage maintenance16. We used immunofluorescence (IF), STARmap-ISH and western blotting to confirm FOXP1 vCTB-specific expression in primary cells and three-dimensional (3D) TB organoids (TB-ORGs). Subsequent small interfering (si)RNA treatment of FOXP1 in TB stem cells (TSCs) reduced proliferation (CCNA2) and activated STB differentiation (CGB and ENDOU) of TSCs (Fig. 2a and Extended Data Fig. 2a).Fig. 2Dynamics of gene regulation in the human placenta.a, Representative IF (n = 4) image (FOXP1/ENDOU; FOXP1/DAPI) of FOXP1 in a W7-1 sample demonstrating nuclear expression restricted to vCTB layers: representative FOXP1 mRNA expression (STARmap-ISH, W7-2) in vCTBs; western blot of isolated placental cell subtypes (villous core, vCTBs, STBs, EVTs, n = 4 donors) detecting FOXP1 and selected specification markers for TB populations (TP63-vCTB; CGβ-STB; HLA-G-EVT) confirming FOXP1 expression in vCTBs. b, Comparison of gene expression (RNA) and gene activity score (ATAC) UMAP plots for the canonical TB marker genes PAGE4 (vCTB) and CYP19A1 (STB). c, RNA and ATAC UMAP plots of new EVT-specific genes uncovered by snATAC-seq, including DIO2 (ensuring local availability of active T3) and MGAT5 (reinforcing TGFβ pathway activity) along with STARmap-ISH spatial expression of DIO2 (n = 3). d, Heatmap showing variable TF motifs across all clusters. Differentially accessible peaks were identified using a two-sided Wilcoxon’s test (FDR ≤ 0.1, log2(FC) ≥ 0.5). Intercluster motif enrichments were calculated via a hypergeometric test to generate P values. e, Significantly enriched TF motifs via ChromVAR enrichment analysis. Specific motifs and enrichment values can be found in Supplementary Table 5. f, Plotted tracks showing regions of open chromatin in enriched motifs of interest identified by ChromVAR, namely TP73 (tumor suppressor preferentially methylated in TBs), SMARCC1 (implicated in STB differentiation and tumor proliferation, no known role in EVTs), SNAI1 (EMT/invasion) and MESP2 (Notch mesodermal specification factor).

ks showing regions of open chromatin in enriched motifs of interest identified by ChromVAR, namely TP73 (tumor suppressor preferentially methylated in TBs), SMARCC1 (implicated in STB differentiation and tumor proliferation, no known role in EVTs), SNAI1 (EMT/invasion) and MESP2 (Notch mesodermal specification factor). TP73 and SMARCC1 overlapped with prior intercluster motif enrichments, whereas SNAI1 and MESP2 were newly uncovered by ChromVAR. Accompanying IF (n = 4) staining in a W9 sample shows prominent SMARCC1 (magenta) expression in EVTs alongside villous core-expressed vimentin (orange). Nuclei are stained with DAPI. g, Positive TF regulators identified based on snRNA-seq gene expression. Specific motifs and enrichment values can be found in Supplementary Table 5. h, Spatial contextualization of cell type-specific enriched motifs for TP63, SNAI1, SMARCC1 and RUNX1 in samples W8-2, W9 and W11 using Slide-tags. Boxes depict magnified enrichment of TP63 motifs in vCTBs and enrichment of SNAI1, SMARCC1 and RUNX1 motifs in EVTs. Max., Maximum.Source data

und in Supplementary Table 5. h, Spatial contextualization of cell type-specific enriched motifs for TP63, SNAI1, SMARCC1 and RUNX1 in samples W8-2, W9 and W11 using Slide-tags. Boxes depict magnified enrichment of TP63 motifs in vCTBs and enrichment of SNAI1, SMARCC1 and RUNX1 motifs in EVTs. Max., Maximum.Source data a, Representative IF (n = 4) image (FOXP1/ENDOU; FOXP1/DAPI) of FOXP1 in a W7-1 sample demonstrating nuclear expression restricted to vCTB layers: representative FOXP1 mRNA expression (STARmap-ISH, W7-2) in vCTBs; western blot of isolated placental cell subtypes (villous core, vCTBs, STBs, EVTs, n = 4 donors) detecting FOXP1 and selected specification markers for TB populations (TP63-vCTB; CGβ-STB; HLA-G-EVT) confirming FOXP1 expression in vCTBs. b, Comparison of gene expression (RNA) and gene activity score (ATAC) UMAP plots for the canonical TB marker genes PAGE4 (vCTB) and CYP19A1 (STB). c, RNA and ATAC UMAP plots of new EVT-specific genes uncovered by snATAC-seq, including DIO2 (ensuring local availability of active T3) and MGAT5 (reinforcing TGFβ pathway activity) along with STARmap-ISH spatial expression of DIO2 (n = 3). d, Heatmap showing variable TF motifs across all clusters. Differentially accessible peaks were identified using a two-sided Wilcoxon’s test (FDR ≤ 0.1, log2(FC) ≥ 0.5). Intercluster motif enrichments were calculated via a hypergeometric test to generate P values. e, Significantly enriched TF motifs via ChromVAR enrichment analysis. Specific motifs and enrichment values can be found in Supplementary Table 5. f, Plotted tracks showing regions of open chromatin in enriched motifs of interest identified by ChromVAR, namely TP73 (tumor suppressor preferentially methylated in TBs), SMARCC1 (implicated in STB differentiation and tumor proliferation, no known role in EVTs), SNAI1 (EMT/invasion) and MESP2 (Notch mesodermal specification factor). TP73 and SMARCC1 overlapped with prior intercluster motif enrichments, whereas SNAI1 and MESP2 were newly uncovered by ChromVAR. Accompanying IF (n = 4) staining in a W9 sample shows prominent SMARCC1 (magenta) expression in EVTs alongside villous core-expressed vimentin (orange). Nuclei are stained with DAPI. g, Positive TF regulators identified based on snRNA-seq gene expression. Specific motifs and enrichment values can be found in Supplementary Table 5. h, Spatial contextualization of cell type-specific enriched motifs for TP63, SNAI1, SMARCC1 and RUNX1 in samples W8-2, W9 and W11 using Slide-tags.

are stained with DAPI. g, Positive TF regulators identified based on snRNA-seq gene expression. Specific motifs and enrichment values can be found in Supplementary Table 5. h, Spatial contextualization of cell type-specific enriched motifs for TP63, SNAI1, SMARCC1 and RUNX1 in samples W8-2, W9 and W11 using Slide-tags. Boxes depict magnified enrichment of TP63 motifs in vCTBs and enrichment of SNAI1, SMARCC1 and RUNX1 motifs in EVTs. Max., Maximum. Source data

are stained with DAPI. g, Positive TF regulators identified based on snRNA-seq gene expression. Specific motifs and enrichment values can be found in Supplementary Table 5. h, Spatial contextualization of cell type-specific enriched motifs for TP63, SNAI1, SMARCC1 and RUNX1 in samples W8-2, W9 and W11 using Slide-tags. Boxes depict magnified enrichment of TP63 motifs in vCTBs and enrichment of SNAI1, SMARCC1 and RUNX1 motifs in EVTs. Max., Maximum. Source data Beyond gene expression, there is limited research on the chromatin accessibility dynamics of the human placenta7. In our dataset, snATAC-seq successfully resolved all major cell types (Fig. 1b); to quantitatively assess these data, we calculated gene activity scores, a metric based on aggregate local chromatin accessibility of genes. Gene expression and gene activity scores of cluster-specific genes were well correlated, supporting the relevance of our snATAC-seq data (Fig. 2b, Extended Data Fig. 2b and Supplementary Fig. 2a). It is interesting that certain genes showed increased cluster-specific chromatin accessibility not mirrored in their expression, potentially because open chromatin precedes transcription (Supplementary Fig. 2b). Further examination of gene activity scores unraveled new cell type-specific genes probably involved in EVT lineage determination (MGAT5), T3-dependent EVT-decidua crosstalk (DIO2), anti-inflammation (ANXA1) and angiogenesis (AGTR1)10,17–19 (Fig. 2c and Supplementary Fig. 2c). EVT-expressed ATP11A, regulating a ‘don’t eat me’ cell membrane signal, may conceal exposure to maternal immune cells and support EVT fusion10,20 (Supplementary Fig. 2c). We found candidate microRNAs (miRNAs) with heightened accessibility in TBs including miR7973-1 (STBs) and miR23B (vCTBs), the expression of which we confirmed through quantitative (q)PCR (Extended Data Fig. 2c, Supplementary Table 3 and Methods). The miR7973-1 is embedded in the intron of STB marker CYP19A1 and targets TB markers ADAM19 (EVT) and FOXP1 (vCTB), potentially safeguarding STB differentiation21. Antioxidant miR23B could help alleviate intervillous maternal blood flow-associated oxidative stress22. Further functional validation of identified miRNAs is crucial for the advancement of this burgeoning field. Finally, we generated a list of gene activity score-derived marker genes (Extended Data Fig. 2d and Supplementary Table 3).

23B could help alleviate intervillous maternal blood flow-associated oxidative stress22. Further functional validation of identified miRNAs is crucial for the advancement of this burgeoning field. Finally, we generated a list of gene activity score-derived marker genes (Extended Data Fig. 2d and Supplementary Table 3). We next aimed to identify key lineage-determining regulatory elements. We first identified transcription factor (TF) motifs differentially enriched between clusters (Methods); to investigate their potential importance, we examined each motif’s corresponding TF gene expression and activity score (accessibility) (Fig. 2d, Extended Data Fig. 2e and Supplementary Table 4). TFs previously implicated in TB development and maintenance (TFAP2C, TP63 and TEAD4) were enriched in vCTBs relative to EVTs and STBs (Fig. 2d and Supplementary Table 4). EVTs were enriched for members of the invasion-related AP-1 signaling network relative to other TBs, including the synergistic FOS and JUN families, STB differentiation-implicated NFE2 and epithelial-to-mesenchymal transition (EMT) facilitator JDP2 (refs. 7,23,24) (Fig. 2d and Supplementary Table 4). Remarkably, distal EVTs (EVT3) and STBs showed lower accessibility of JDP2 than other TBs, indicating that transcription of JDP2 may halt on reaching terminal differentiation states (Extended Data Fig. 2e). The nuclear factor κ-light chain enhancer of activated B cells pathway members BACH1 and BACH2 were enriched in EVTs and may impact immune tolerance7,25 (Supplementary Table 4). EVTs showed striking enrichment for SMARCC1, implicated in STB differentiation and tumor proliferation, but with an unspecified role in EVT differentiation26 (Fig. 2d and Supplementary Table 4). Enriched STB motifs with unclear roles included ESRRG, ATF4 and the CEBP family (Fig. 2d and Supplementary Table 4).

able 4). EVTs showed striking enrichment for SMARCC1, implicated in STB differentiation and tumor proliferation, but with an unspecified role in EVT differentiation26 (Fig. 2d and Supplementary Table 4). Enriched STB motifs with unclear roles included ESRRG, ATF4 and the CEBP family (Fig. 2d and Supplementary Table 4). Within EVT subclusters, EVT3 was enriched for BACH1 and FOSB, which showed corresponding EVT3-specific gene expression indicating potential markers of distal differentiation (Extended Data Fig. 2e and Supplementary Table 4). Among the vCTB subclusters, vCTB2 was enriched for members of the KLF/SP family, which have obscure TB-related roles, and PBX3, a proto-oncogenic TF that caused decreased invasion when downregulated in pre-eclampsia (PE)27 (Supplementary Table 4). The vCTB3 was markedly enriched for CTCF, previously implicated in PE28 (Fig. 2d and Supplementary Table 4). The epigenomic landscape of the placental stromal core has not been well described; relative to FIB1, FIB2 was prominently enriched for many of the same motifs as TBs, as well as malignancy-related TFAP4 and EMT-related ASCL2, components of the Wnt pathway crucial for placental function29–31 (Fig. 2d and Supplementary Table 4). FIB1 enrichments comprised pericyte-associated EBF1 (downregulated in spontaneous preterm birth), newly identified immunosuppressive NR2F6 and malignancy-related invasion regulator ZNF148, all without clear roles in the placenta32–35 (Fig. 2d and Supplementary Table 4).

unction29–31 (Fig. 2d and Supplementary Table 4). FIB1 enrichments comprised pericyte-associated EBF1 (downregulated in spontaneous preterm birth), newly identified immunosuppressive NR2F6 and malignancy-related invasion regulator ZNF148, all without clear roles in the placenta32–35 (Fig. 2d and Supplementary Table 4). We subsequently used ChromVAR36 to calculate TF motif enrichment on a per-cell basis, accounting for Tn5 insertion bias, and then examined the correlation between enriched motifs and corresponding TF gene expression and activity score (Fig. 2e, Extended Data Fig. 2e and Supplementary Table 5). Highly variable motifs identified by ChromVAR corresponded with previously identified motifs (Fig. 2e,f and Supplementary Fig. 3a). New motifs identified by ChromVAR included ZEB1, implicated in EMT37, and TB-expressed ID3 and ID4, which may impact TB invasion38 (Supplementary Fig. 3a). Other TB motifs consisted of SNAI1 and SNAI2 (EMT/invasion) along with MESP1 and MESP2 (Wnt and Notch mesodermal specification factors, respectively)39–42 (Fig. 2f). We next identified 40 positive TF regulators, defined as variant ChromVAR-identified motifs averaged across clusters, with a positive correlation between motif accessibility changes and expression of the corresponding genes (Fig. 2g and Supplementary Table 5). These correlations suggest potential regulatory interactions. Using gene activity score in lieu of gene expression produced 55 positive TF regulators, including an approximately 50% overlap with the gene expression-derived regulators (Extended Data Fig. 2f and Supplementary Table 5). Several previously mentioned motifs were deemed positive TF regulators, along with TB regulator GCM1 (ref. 43), vascular regulator HOXA13 (ref. 44) and EVT invasion regulator HMGA1 (ref. 45) (supporting its recently discovered importance in PE) (Fig. 2g and Extended Data Fig. 2e,f). We also identified potentially key motifs involved in PE and malignancy, with roles in immunomodulation, invasiveness and stromal remodeling (Supplementary Fig. 3b).

r HOXA13 (ref. 44) and EVT invasion regulator HMGA1 (ref. 45) (supporting its recently discovered importance in PE) (Fig. 2g and Extended Data Fig. 2e,f). We also identified potentially key motifs involved in PE and malignancy, with roles in immunomodulation, invasiveness and stromal remodeling (Supplementary Fig. 3b). Finally, we used Slide-tags8 to better understand how variation in chromatin structure and gene expression manifests spatially across the placenta (Figs. 1e and 2h, Supplementary Figs. 4 and 5 and Supplementary Table 6). We observed finely resolved motif deviation score enrichments between closely located tissue microenvironments, such as the highly EVT-specific RUNX1 motif enrichment in close physical proximity to vCTBs (Fig. 2h). We investigated which genes and motifs showed location-dependent expression and enrichment across all cells and within cell types, elucidating spatial expression and enrichment gradients throughout the placenta (Supplementary Table 6).

EVT-specific RUNX1 motif enrichment in close physical proximity to vCTBs (Fig. 2h). We investigated which genes and motifs showed location-dependent expression and enrichment across all cells and within cell types, elucidating spatial expression and enrichment gradients throughout the placenta (Supplementary Table 6). Seeking to unravel the relationship between genes and distal CREs, we linked distal peaks to genes in cis and identified 43,622 significant peak–gene associations that represent potential enhancer–gene relationships, with an average of 54 peaks per gene (Fig. 3a, Supplementary Table 7 and Methods). Following ref. 46, we distinguished the 1,057 regions with >10 significant peak–gene associations as DORCs, areas that point to regulatory locus complexity for certain genes (Supplementary Table 7). Essential vascular regulators, FOXF1 and RTL1, had the most associated peaks with 183 and 179, respectively47,48 (Fig. 3b and Supplementary Table 7). Other highly associated DORCs included stromal markers DLK1, CXCL14 and COL6A3, as well as maternally expressed imprinted genes MEG3 and MEG8, the abnormal expression of which has been implicated, for both, in TB dysfunction49,50 (Fig. 3c, Extended Data Fig. 3a, Supplementary Fig. 6 and Supplementary Table 7). Many previously mentioned tumor invasion genes, immunomodulation-related genes and snATAC-seq-identified genes were found among DORCs (Fig. 3b,c, Extended Data Fig. 3a and Supplementary Table 7). We next calculated DORC scores across cell types, defined by the normalized sum of counts in all significantly correlated peaks per gene for all cells, and found that DORCs were highly cell type specific (Fig. 3c, Extended Data Fig. 3a,b and Supplementary Table 7). We identified TB-specific DORCs potentially vital to TB function such as KANK1 (tumor-suppressing YAP1 regulator), PARD6B (cell polarity regulator involved in tumor growth), FOXI3 (metastasis related) and MYCN/MYCNUT (growth promoter countered by TP63 (expressed in vCTBs), which might restrict MYCNUT/MYCN expression to EVTs)26,51–54 (Fig. 3c, Extended Data Fig. 3a, Supplementary Fig. 6 and Supplementary Table 7). In EVTs, we identified TEA domain proteins (TEAD) family binder VGLL3; EVTs are known to express high levels of its target IGFBP3, indicating that a VGLL3–TEAD1 axis might regulate insulin-like growth factor (IGF) at the maternal–fetal interface55,56 (Extended Data Fig. 3a,b, Supplementary Fig. 6 and Supplementary Table 7).

Ts, we identified TEA domain proteins (TEAD) family binder VGLL3; EVTs are known to express high levels of its target IGFBP3, indicating that a VGLL3–TEAD1 axis might regulate insulin-like growth factor (IGF) at the maternal–fetal interface55,56 (Extended Data Fig. 3a,b, Supplementary Fig. 6 and Supplementary Table 7). In stroma, FOXF1 was excluded from mat. FIB, indicating that FOXF1 might be important in early fetal mesenchymal lineage determination57. It is interesting that we identified cardiac lineage regulator TBX5 in FIB1 and FIB2, supporting the rising evidence for a link between cardiac and placental development58 (Extended Data Fig. 3a and Supplementary Table 7).Fig. 3Identifying interactions between CREs and genes.a, Heatmap showing chromatin accessibility (left) and gene expression (right) of 43,622 significantly linked CRE–gene pairs across all cell types (top row). Peak–gene interactions were clustered using k-means clustering into 25 bins (k = 25). Significant peak–gene links can be found in Supplementary Table 7. b, Track visualization of FOXF1 and QSOX1, which showed additional chromatin accessibility in: endo, FIB1, FIB2, Unknown 1 (FOXF1) and EVT1-3 (QSOX1). Each track is accompanied by all significant peak–gene links identified for that gene. c, UMAP plots of DORC score and gene expression (RNA) showing correlated expression patterns of ANXA1 (vCTB3, STB), KANK1 (vCTB1–3), HLA-G (EVT2/3), QSOX1 (EVT2/3), MYCN (EVT2) and DLK1 (FIB1/2). DORC scores are defined by the normalized sum of counts in all significantly correlated peaks per gene for all cells. DORCs and DORC scores can be found in Supplementary Table 7. d, Stream plot of chromatin potential revealing a dynamic TB differentiation process with multifocal roots. e, CellRank applied to chromatin potential (Methods) revealing initial (left) and terminal (right) TB differentiation states. Both vCTB2 and EVT2 were identified as initial states, whereas EVT3, STB and vCTB2 were classified as terminal states. f, Heritability enrichment score (Escore) of the number of miscarriages trait for different cell types based on (1) sc-linker (ABC + Roadmap), (2) Multiome and (3) sc-linker (Multiome). Enrichment scores were calculated using S-LDSC (FDR < 0.05). Numerical results are reported in Supplementary Table 9. g, Enrichment analysis of GWAS hits for nine placenta-related traits in SNPs in peaks linked to any gene in a cell type (Multiome). Enrichment scores were calculated using S-LDSC (FDR < 0.05).

ultiome). Enrichment scores were calculated using S-LDSC (FDR < 0.05). Numerical results are reported in Supplementary Table 9. g, Enrichment analysis of GWAS hits for nine placenta-related traits in SNPs in peaks linked to any gene in a cell type (Multiome). Enrichment scores were calculated using S-LDSC (FDR < 0.05). The numerical results are reported in Supplementary Table 10. h, Enrichment analysis of GWAS hits for nine placenta-related traits in SNPs in peaks linked to cell type-specific genes (sc-linker (Multiome)). Enrichment scores were calculated using S-LDSC (FDR < 0.05). The numerical results are reported in Supplementary Table 10.

ical results are reported in Supplementary Table 10. h, Enrichment analysis of GWAS hits for nine placenta-related traits in SNPs in peaks linked to cell type-specific genes (sc-linker (Multiome)). Enrichment scores were calculated using S-LDSC (FDR < 0.05). The numerical results are reported in Supplementary Table 10. a, Heatmap showing chromatin accessibility (left) and gene expression (right) of 43,622 significantly linked CRE–gene pairs across all cell types (top row). Peak–gene interactions were clustered using k-means clustering into 25 bins (k = 25). Significant peak–gene links can be found in Supplementary Table 7. b, Track visualization of FOXF1 and QSOX1, which showed additional chromatin accessibility in: endo, FIB1, FIB2, Unknown 1 (FOXF1) and EVT1-3 (QSOX1). Each track is accompanied by all significant peak–gene links identified for that gene. c, UMAP plots of DORC score and gene expression (RNA) showing correlated expression patterns of ANXA1 (vCTB3, STB), KANK1 (vCTB1–3), HLA-G (EVT2/3), QSOX1 (EVT2/3), MYCN (EVT2) and DLK1 (FIB1/2). DORC scores are defined by the normalized sum of counts in all significantly correlated peaks per gene for all cells. DORCs and DORC scores can be found in Supplementary Table 7. d, Stream plot of chromatin potential revealing a dynamic TB differentiation process with multifocal roots. e, CellRank applied to chromatin potential (Methods) revealing initial (left) and terminal (right) TB differentiation states. Both vCTB2 and EVT2 were identified as initial states, whereas EVT3, STB and vCTB2 were classified as terminal states. f, Heritability enrichment score (Escore) of the number of miscarriages trait for different cell types based on (1) sc-linker (ABC + Roadmap), (2) Multiome and (3) sc-linker (Multiome). Enrichment scores were calculated using S-LDSC (FDR < 0.05). Numerical results are reported in Supplementary Table 9. g, Enrichment analysis of GWAS hits for nine placenta-related traits in SNPs in peaks linked to any gene in a cell type (Multiome). Enrichment scores were calculated using S-LDSC (FDR < 0.05). The numerical results are reported in Supplementary Table 10. h, Enrichment analysis of GWAS hits for nine placenta-related traits in SNPs in peaks linked to cell type-specific genes (sc-linker (Multiome)). Enrichment scores were calculated using S-LDSC (FDR < 0.05). The numerical results are reported in Supplementary Table 10.

ical results are reported in Supplementary Table 10. h, Enrichment analysis of GWAS hits for nine placenta-related traits in SNPs in peaks linked to cell type-specific genes (sc-linker (Multiome)). Enrichment scores were calculated using S-LDSC (FDR < 0.05). The numerical results are reported in Supplementary Table 10. The vCTB differentiation trajectories have been explored through scRNA-seq5,7 but have lacked accompanying chromatin accessibility, which comprises select enhancer activities that foreshadow specific gene expression. We employed chromatin potential46, a technique that predicts future cell states on a greater timescale than RNA-based methods, in combination with CellRank59 to calculate transition probabilities, identifying initial (EVT2 and vCTB2) and terminal (EVT2,3 and STB) states (Fig. 3d,e and Extended Data Fig. 3c). The classification of vCTB2 as an initial and terminal state may reflect its role as the predominant vCTB subtype, constituting a source for self-renewal, proliferation and STB/EVT differentiation (Fig. 3e and Extended Data Fig. 3c). EVT2′s inclusion as an initial state, paired with the stream plot’s bidirectional flow of EVT2 into EVT2 and EVT3, supports EVT2 as a source for EVT3 differentiation and EVT2 self-renewal. Together, these findings suggest both EVTs and vCTBs as progenitors for EVT differentiation. Finally, we identified lineage drivers for each terminal state; highly correlated drivers included both canonical genes and previously mentioned new genes QSOX1 and ANXA1 (Supplementary Fig. 7 and Supplementary Table 8).

newal. Together, these findings suggest both EVTs and vCTBs as progenitors for EVT differentiation. Finally, we identified lineage drivers for each terminal state; highly correlated drivers included both canonical genes and previously mentioned new genes QSOX1 and ANXA1 (Supplementary Fig. 7 and Supplementary Table 8). Expanding on our multimodal analyses, we analyzed potential disease risk for nine UK Biobank (UKBB) pregnancy-related traits (Supplementary Table 9 and Methods). We assessed disease heritability informativeness of four SNP annotation types (LDSC-SEG, sc-linker (activity-by-contact (ABC) + Roadmap), Multiome, sc-linker (Multiome)), examined the enrichment of GWAS hits for pregnancy-related traits in these four annotation types and performed gene set enrichment analysis (GSEA) of placental cell-type gene programs using multi-marker analysis of genomic annotation (MAGMA) (Methods).

ivity-by-contact (ABC) + Roadmap), Multiome, sc-linker (Multiome)), examined the enrichment of GWAS hits for pregnancy-related traits in these four annotation types and performed gene set enrichment analysis (GSEA) of placental cell-type gene programs using multi-marker analysis of genomic annotation (MAGMA) (Methods). We observed, overall, higher heritability enrichments on average across pregnancy-related traits and cell types for the Multiome and sc-linker (Multiome) annotations (Extended Data Fig. 3d and Supplementary Table 9). For the trait of number of miscarriages, we observed significant enrichment signals for vCTBs, STBs, endos and macrophages (Fig. 3f and Supplementary Table 9). We also performed MAGMA GSEA of cell type-specific programs and observed the strongest enrichment signal for FIB1 and gestational hypertension (P = 0.0009) (Supplementary Table 10). As a secondary analysis, we performed an enrichment analysis of the top GWAS hits (P < 0.05) (Fig. 3g,h and Supplementary Table 10). Of the 4,827 total GWAS hits across these pregnancy-related traits, 324 hits were found in gene-linked peaks in some placental cell types. Excessive vomiting during pregnancy showed the highest average enrichment across cell types (Fig. 3g and Supplementary Table 10). One notable GWAS hit for excessive vomiting during pregnancy is rs117659937, located in the intronic region of TP53INP2 and residing in a peak linked to TP53INP2 in EVT3 (Extended Data Fig. 3e). A potential role for TP53 in EVT differentiation and/or function is unknown. In cell-type gene program-specific enrichment analyses, a variety of enrichments was observed across seven out of nine traits (Fig. 3h and Supplementary Table 10).

INP2 and residing in a peak linked to TP53INP2 in EVT3 (Extended Data Fig. 3e). A potential role for TP53 in EVT differentiation and/or function is unknown. In cell-type gene program-specific enrichment analyses, a variety of enrichments was observed across seven out of nine traits (Fig. 3h and Supplementary Table 10). Biological intricacies predominantly manifest in spatial cellular arrangements of niches, physiological gradients and cell–cell interactions, which can be elucidated through single-cell spatial technologies. To map the early placenta, we performed STARmap on snap-frozen human placenta tissue sections from four of the eight previously described donors: W7-1, W8-2, W9 and W11 (Methods). First, we utilized STARmap-ISS and measured the spatial expression of 1,001 genes, composed of both highly variable genes from paired snRNA-seq and snATAC-seq and manually selected markers, and imaged a representative region of sectioned tissue spanning across the entire placenta (Fig. 1f and Supplementary Table 11). Cell-type labels were projected on to spatially resolved cells and manually verified (Fig. 4a–e, Extended Data Fig. 4a–c and Methods). Out of 17 clusters, 12 were predicted with 16,179 cells (Fig. 4a, Extended Data Fig. 5 and Supplementary Table 1). In addition, we identified a new list of differentially expressed genes (DEGs) from each STARmap cluster to examine spatial expression patterns of discovered marker genes (Fig. 4f and Supplementary Table 12). To investigate additional manually selected genes across large whole-tissue sections, STARmap-ISH analyses on 48 landmark genes were performed on 3 samples (W7-2, W8-2 and W11) for a total of 868,920 cells (Fig. 1g, Supplementary Figs. 8 and 9 and Supplementary Tables 1 and 11). Classic placental morphology across all samples was visualized with marker genes (Fig. 4d,e and Extended Data Fig. 4b,c). STARmap-ISH enabled additional spatial verification of Multiome-identified cell clusters and identification of TB-expressed tumor and immunomodulation-associated genes (Fig. 4e and Extended Data Fig. 6). STARmap-ISS and STARmap-ISH exhibited similar gene expression patterns, reinforcing ISS results.Fig. 4Spatial transcriptomic landscape of the early human placenta.a, UMAP representation of STARmap-ISS-derived clusters of the placenta. QC metrics can be found in Supplementary Table 1. STARmap genes are detailed in Supplementary Table 11. b, UMAP showing all cells of STARmap-ISS (1,001 genes).

inforcing ISS results.Fig. 4Spatial transcriptomic landscape of the early human placenta.a, UMAP representation of STARmap-ISS-derived clusters of the placenta. QC metrics can be found in Supplementary Table 1. STARmap genes are detailed in Supplementary Table 11. b, UMAP showing all cells of STARmap-ISS (1,001 genes). Cells are colored according to the individual sample. Individual placental samples W7-1, W9, W11 and W8-2 contributed equally to identified cell clusters. c, Multiomics-derived cell types identified in space with STARmap-ISS via Seurat integration on sample W8-2. d, Canonical marker gene spatial expression for vCTB (PAGE4), STB (CGA), EVT (NOTUM) and the stromal core (VIM) identified by STARmap-ISS (n = 4) on sample W8-2 with insets 1–4 showing magnified villous core (1) and cell columns (2–4). e, Canonical marker gene spatial expression in W8-2 for FIB (COL3A1), STB (CYP19A1), EVT (HAPLN3, AOC1), HBC (CD163) and endo (KDR) identified by STARmap-ISS (n = 4), along with detection of multiple genes in STARmap-ISH allowing distinct spatial allocation of cell clusters. f, Heatmap showing DEGs identified by STARmap-ISS across clusters; each row represents the top ten DEGs for that cluster. DEGs for each cluster are listed in Supplementary Table 12.

ified by STARmap-ISS (n = 4), along with detection of multiple genes in STARmap-ISH allowing distinct spatial allocation of cell clusters. f, Heatmap showing DEGs identified by STARmap-ISS across clusters; each row represents the top ten DEGs for that cluster. DEGs for each cluster are listed in Supplementary Table 12. a, UMAP representation of STARmap-ISS-derived clusters of the placenta. QC metrics can be found in Supplementary Table 1. STARmap genes are detailed in Supplementary Table 11. b, UMAP showing all cells of STARmap-ISS (1,001 genes). Cells are colored according to the individual sample. Individual placental samples W7-1, W9, W11 and W8-2 contributed equally to identified cell clusters. c, Multiomics-derived cell types identified in space with STARmap-ISS via Seurat integration on sample W8-2. d, Canonical marker gene spatial expression for vCTB (PAGE4), STB (CGA), EVT (NOTUM) and the stromal core (VIM) identified by STARmap-ISS (n = 4) on sample W8-2 with insets 1–4 showing magnified villous core (1) and cell columns (2–4). e, Canonical marker gene spatial expression in W8-2 for FIB (COL3A1), STB (CYP19A1), EVT (HAPLN3, AOC1), HBC (CD163) and endo (KDR) identified by STARmap-ISS (n = 4), along with detection of multiple genes in STARmap-ISH allowing distinct spatial allocation of cell clusters. f, Heatmap showing DEGs identified by STARmap-ISS across clusters; each row represents the top ten DEGs for that cluster. DEGs for each cluster are listed in Supplementary Table 12.

ified by STARmap-ISS (n = 4), along with detection of multiple genes in STARmap-ISH allowing distinct spatial allocation of cell clusters. f, Heatmap showing DEGs identified by STARmap-ISS across clusters; each row represents the top ten DEGs for that cluster. DEGs for each cluster are listed in Supplementary Table 12. All samples contained comparable placental structures comprising villi and EVT cell columns (CCs); the latter was visualized by uniformly expressed canonical genes and genes identified by our multiomic analysis, including NOTUM, QSOX1 and BACH2 (Figs. 4d and 5a and Extended Data Fig. 4b). Notably, W8-2 exhibited CCs of different EVT maturation status (CC type 1, CC type 2 and CC type 3) composed of varying percentages of EVT1 (CC type 1), EVT2 (CC type 2) or EVT3 (CC type 3) (Extended Data Fig. 7). We detected multiple EVT-specific genes with decreased expression as lineage differentiation progressed, including CDK7, ERVH48-1, HAPLN3, ITGB4, LY6E, MYCNUT and ATP11A (Figs. 4e and 5a,b and Extended Data Figs. 1k and 7). This may indicate timely, spatially limited roles in EVT development. CC type 3 showed high levels of ADAM19, ITGB4, AOC1, SLCO2A1 and PAPPA2 associated with a mature EVT predominantly detected in the maternal decidua10. We detected DORC-identified MYCNUT in CC type 1/2 EVTs, supporting our hypothesis of TP63-MYCNUT/MYCN-regulated EVT formation (Fig. 5a). STARmap analyses identified yet unknown vCTB-enriched genes MSI2, ERVH48-1, EFEMP1 (previously identified vCTB2 lineage driver) and KANK1; the last was identified by DORC analyses (Fig. 5b). MSI2 is crucial to ESC self-renewal and might represent a new TF in human placentation60. ERVH48-1, a truncated variant of the placenta-specific human endogenous retrovirus (HERV) family localized to EVT and vCTB, may prevent premature TB fusion by competing with ERVW-1 for receptor-binding sites; this mechanism would be in direct contrast to STB-inducing HERV family members61. Surprisingly, prime epithelial and TB markers KRT7 and KRT19 were strongly increased in EVTs whereas KRT23 was restricted to vCTBs (Extended Data Fig. 8a). We also detected certain genes shared between EVTs and either STBs or vCTBs (Extended Data Fig. 8b). Notably, we detected high levels of tissue factor pathway inhibitor (TFPI) (EVTs) and its paralog TFPI2 (STBs) (Fig. 5c). Strong gene expression of these two major anticoagulants may maintain homeostasis of the fibrinolytic system at both the intervillous space and the MFI.

either STBs or vCTBs (Extended Data Fig. 8b). Notably, we detected high levels of tissue factor pathway inhibitor (TFPI) (EVTs) and its paralog TFPI2 (STBs) (Fig. 5c). Strong gene expression of these two major anticoagulants may maintain homeostasis of the fibrinolytic system at both the intervillous space and the MFI. STBs expressed typical markers GDF15, CGA and CYP19A1; however, we also detected BMP2, MSX1 and PXDN (Fig. 5d). BMP2 has predominantly been described in EVT functions62; consequently, STBs might represent a source for BMP2 or BMP2 may have a unique role in STB function. PXDN contributes to collagen IV-associated extracellular matrix (ECM) assembly and, like BMP2, might be secreted by STBs into the intervillous space to support ECM formation63. It is interesting that STBs strongly expressed the transcriptional repressor MSX1, a paralog of the STB lineage repressor MSX2, indicating that MSX1 might have an opposite function in TB differentiation.Fig. 5Spatial mapping of the early human placenta at molecular resolution.a, Spatial verification of EVT-expressed genes identified in previous analyses, including QSOX1 (ATAC), BACH2 (motif enrichment) and MYCNUT (DORC) via STARmap-ISS (n = 4) in samples W8-2 and W9. Insets depict magnified areas co-visualizing genes of interest with the canonical marker genes PAGE4 (vCTB), VIM (villous core) and NOTUM (EVT), depicting EVT-specific expression of QSOX1, BACH2 and MYCNUT. b, Spatial characterization of new, vCTB-specific genes identified by STARmap-ISS (n = 4) (MSI2, ERVH48-1 and EFEMP1) and DORC (KANK1) analyses in samples W8-2 and W9. Insets depict magnified areas co-visualizing genes of interest with the canonical marker genes CGA (STB), PAGE4 (vCTB) and NOTUM (EVT), verifying vCTB-specific (MSI2, EFEMP1 and KANK1) and vCTB-/EVT-specific (ERVH48-1) expression. Please note the absence of these genes in STB. Stippled lines demarcate vCTBs from the villous core. c, Spatial characterization of TFPI and TFPI2, two crucial hemodynamic modulators expressed predominantly in EVTs (TFPI) and STBs (TFPI2), via STARmap-ISS (n = 4) across samples W8-2, W9 and W7-1. The insets depict magnified areas with co-visualization of TFP1, TFP2, VIM (villous core) and NOTUM (EVT). Stippled lines demarcate vCTBs from the villous core. d, Spatial visualization of canonical (GDF15) and new (MSX1, PXDN and BMP2) STB-specific markers via STARmap-ISS (n = 4) in samples W8-2, W9 and W7-1.

nd W7-1. The insets depict magnified areas with co-visualization of TFP1, TFP2, VIM (villous core) and NOTUM (EVT). Stippled lines demarcate vCTBs from the villous core. d, Spatial visualization of canonical (GDF15) and new (MSX1, PXDN and BMP2) STB-specific markers via STARmap-ISS (n = 4) in samples W8-2, W9 and W7-1. The insets depict magnified areas with co-visualization of genes of interest and the canonical markers PAGE4 (vCTB) and VIM (villous core). Stippled lines demarcate vCTBs from the villous core.

nd W7-1. The insets depict magnified areas with co-visualization of TFP1, TFP2, VIM (villous core) and NOTUM (EVT). Stippled lines demarcate vCTBs from the villous core. d, Spatial visualization of canonical (GDF15) and new (MSX1, PXDN and BMP2) STB-specific markers via STARmap-ISS (n = 4) in samples W8-2, W9 and W7-1. The insets depict magnified areas with co-visualization of genes of interest and the canonical markers PAGE4 (vCTB) and VIM (villous core). Stippled lines demarcate vCTBs from the villous core. a, Spatial verification of EVT-expressed genes identified in previous analyses, including QSOX1 (ATAC), BACH2 (motif enrichment) and MYCNUT (DORC) via STARmap-ISS (n = 4) in samples W8-2 and W9. Insets depict magnified areas co-visualizing genes of interest with the canonical marker genes PAGE4 (vCTB), VIM (villous core) and NOTUM (EVT), depicting EVT-specific expression of QSOX1, BACH2 and MYCNUT. b, Spatial characterization of new, vCTB-specific genes identified by STARmap-ISS (n = 4) (MSI2, ERVH48-1 and EFEMP1) and DORC (KANK1) analyses in samples W8-2 and W9. Insets depict magnified areas co-visualizing genes of interest with the canonical marker genes CGA (STB), PAGE4 (vCTB) and NOTUM (EVT), verifying vCTB-specific (MSI2, EFEMP1 and KANK1) and vCTB-/EVT-specific (ERVH48-1) expression. Please note the absence of these genes in STB. Stippled lines demarcate vCTBs from the villous core. c, Spatial characterization of TFPI and TFPI2, two crucial hemodynamic modulators expressed predominantly in EVTs (TFPI) and STBs (TFPI2), via STARmap-ISS (n = 4) across samples W8-2, W9 and W7-1. The insets depict magnified areas with co-visualization of TFP1, TFP2, VIM (villous core) and NOTUM (EVT). Stippled lines demarcate vCTBs from the villous core. d, Spatial visualization of canonical (GDF15) and new (MSX1, PXDN and BMP2) STB-specific markers via STARmap-ISS (n = 4) in samples W8-2, W9 and W7-1. The insets depict magnified areas with co-visualization of genes of interest and the canonical markers PAGE4 (vCTB) and VIM (villous core). Stippled lines demarcate vCTBs from the villous core.

of canonical (GDF15) and new (MSX1, PXDN and BMP2) STB-specific markers via STARmap-ISS (n = 4) in samples W8-2, W9 and W7-1. The insets depict magnified areas with co-visualization of genes of interest and the canonical markers PAGE4 (vCTB) and VIM (villous core). Stippled lines demarcate vCTBs from the villous core. Various collagen genes were identified in our DORC studies in the top 10% of peak–gene links and several genes were highly expressed in the villous core (Fig. 4e, Extended Data Fig. 4c and Supplementary Table 7). Notably, some collagen genes were also shared between villous core cells and EVTs, probably contributing to the fibrin-matrix-type fibrinoid constituting the maternofetal border as part of the decidua basalis (Extended Data Fig. 8c).

ressed in the villous core (Fig. 4e, Extended Data Fig. 4c and Supplementary Table 7). Notably, some collagen genes were also shared between villous core cells and EVTs, probably contributing to the fibrin-matrix-type fibrinoid constituting the maternofetal border as part of the decidua basalis (Extended Data Fig. 8c). Next, we augmented our resolution to whole-transcriptome scale by imputing single-cell transcriptomic profiles using a mutual nearest neighbors imputation method to borrow information for STARmap data from paired snRNA-seq and snATAC-seq reference data in joint embeddings (Fig. 6a and Methods). We tested and determined the optimal nearest neighbors (Fig. 6a and Extended Data Fig. 9a). Finally, we imputed 33,357 gene expression profiles for 16,179 cells in our STARmap dataset. We evaluated transcriptomic imputation performance by using genes in both STARmap-ISH and STARmap-ISS analyses (Fig. 6b and Extended Data Fig. 9b). Using other STARmap-ISH genes that are not part of STARmap-ISS analyses, we observed a strong correlation between ISH expression and imputed gene expression (Fig. 6c, Extended Data Fig. 9c and Supplementary Fig. 10). We next leveraged the fact that each single-cell transcriptome was paired with a single-cell chromatin accessibility landscape, allowing us to impute epigenetic profiles from multiome snATAC-seq reference data. We evaluated this imputation by comparing gene activity scores of STARmap genes with their measured STARmap-ISS expression (Fig. 6d and Extended Data Fig. 9d). In addition, we found that we could identify DORC scores, spatially resolving complex peak–gene interactions across the placenta (Fig. 6d).Fig. 6A reconstructed spatial multiomic landscape of the early human placenta.a, Schematic of the imputation workflow for generating 33,357 gene expression profiles across 16,179 STARmap-ISS cells. b, Comparison of STARmap-ISS (n = 4), imputation-based gene expression (Imputed-RNA) and STARmap-ISH (n = 3) for the canonical marker genes CYP19A1 (STB) and HAPLN3 (EVT) in sample W8-2, validating imputation. Stippled lines demarcate the villous core from the vCTB layers. Scale bar, 50 µm. c, Visualization of Imputed-RNA compared with STARmap-ISH (n = 3) for HMGB2, TFRC and SLCO2A1 in sample W8-2. Please note a similar expression pattern between detected gene expression (STARmap-ISH) and imputed gene expression (Imputed-RNA). Stippled lines demarcate the villous core from vCTB layers. Scale bar, 50 µm.

. c, Visualization of Imputed-RNA compared with STARmap-ISH (n = 3) for HMGB2, TFRC and SLCO2A1 in sample W8-2. Please note a similar expression pattern between detected gene expression (STARmap-ISH) and imputed gene expression (Imputed-RNA). Stippled lines demarcate the villous core from vCTB layers. Scale bar, 50 µm. d, Spatial visualization of imputed DORC scores (Imputed-DORC), gene activity scores (Imputed-ATAC) and gene expression (STARmap-ISS) in sample W8-2 for EVT markers (QSOX1 and HAPLN3), pan-TB markers (KRT7 and KRT19), vCTB/EVT marker (ERVH48-1) and FIB marker (COL3A1). e, Pie chart showing ligand–receptor interactions in samples W7-1, W9, W11 and W8-2. More information about identified ligand–receptor interactions can be found in Supplementary Table 13. f, Spatial visualization of ligand–receptor pairs EGFR (vCTB)–AREG (FIB) and BMP6 (HBC)–BMPR1B (FIB) via STARmap-ISS (n = 4) in samples W8-2, W9 and W7-1 to support CellChat findings. The stippled lines demarcate the villous core from vCTB layers.

ied ligand–receptor interactions can be found in Supplementary Table 13. f, Spatial visualization of ligand–receptor pairs EGFR (vCTB)–AREG (FIB) and BMP6 (HBC)–BMPR1B (FIB) via STARmap-ISS (n = 4) in samples W8-2, W9 and W7-1 to support CellChat findings. The stippled lines demarcate the villous core from vCTB layers. a, Schematic of the imputation workflow for generating 33,357 gene expression profiles across 16,179 STARmap-ISS cells. b, Comparison of STARmap-ISS (n = 4), imputation-based gene expression (Imputed-RNA) and STARmap-ISH (n = 3) for the canonical marker genes CYP19A1 (STB) and HAPLN3 (EVT) in sample W8-2, validating imputation. Stippled lines demarcate the villous core from the vCTB layers. Scale bar, 50 µm. c, Visualization of Imputed-RNA compared with STARmap-ISH (n = 3) for HMGB2, TFRC and SLCO2A1 in sample W8-2. Please note a similar expression pattern between detected gene expression (STARmap-ISH) and imputed gene expression (Imputed-RNA). Stippled lines demarcate the villous core from vCTB layers. Scale bar, 50 µm. d, Spatial visualization of imputed DORC scores (Imputed-DORC), gene activity scores (Imputed-ATAC) and gene expression (STARmap-ISS) in sample W8-2 for EVT markers (QSOX1 and HAPLN3), pan-TB markers (KRT7 and KRT19), vCTB/EVT marker (ERVH48-1) and FIB marker (COL3A1). e, Pie chart showing ligand–receptor interactions in samples W7-1, W9, W11 and W8-2. More information about identified ligand–receptor interactions can be found in Supplementary Table 13. f, Spatial visualization of ligand–receptor pairs EGFR (vCTB)–AREG (FIB) and BMP6 (HBC)–BMPR1B (FIB) via STARmap-ISS (n = 4) in samples W8-2, W9 and W7-1 to support CellChat findings. The stippled lines demarcate the villous core from vCTB layers.

ied ligand–receptor interactions can be found in Supplementary Table 13. f, Spatial visualization of ligand–receptor pairs EGFR (vCTB)–AREG (FIB) and BMP6 (HBC)–BMPR1B (FIB) via STARmap-ISS (n = 4) in samples W8-2, W9 and W7-1 to support CellChat findings. The stippled lines demarcate the villous core from vCTB layers. Finally, we applied CellChat64 to infer 4,263 ligand–receptor interactions using both transcriptional profiles and spatial context (Fig. 6e, Extended Data Fig. 10a, Supplementary Table 13 and Methods). CellChat identified canonical extraplacental ligand–receptor interactions between villous core cells and TBs, including ligands AREG (FIB2), HGF (FIB2) and HBEGF (HBC) with the receptors EGFR and MET (vCTB2), vital to vCTB survival and STB formation65,66 (Fig. 6f, Extended Data Fig. 10a,b and Supplementary Table 13). We found ligand–receptor interactions between PTN (EVT2) and migration-related SDC4 (EVT3), indicating a potential signaling pathway contributing to EVT invasiveness67,68. Other new ligand–receptor interactions included ANXA (FIB2)–TPR1 (HBC), BMP6 (HBC)–BMPR1B (FIB2) and VEGFA (FIB2)–KDR (endo), revealing mutual crosstalk between villous core cells probably supporting angiogenesis and stromal cell function, which require further investigation (Fig. 6f, Extended Data Fig. 10a and Supplementary Table 13).

Seeking to unravel the relationship between genes and distal CREs, we linked distal peaks to genes in cis and identified 43,622 significant peak–gene associations that represent potential enhancer–gene relationships, with an average of 54 peaks per gene (Fig. 3a, Supplementary Table 7 and Methods). Following ref. 46, we distinguished the 1,057 regions with >10 significant peak–gene associations as DORCs, areas that point to regulatory locus complexity for certain genes (Supplementary Table 7). Essential vascular regulators, FOXF1 and RTL1, had the most associated peaks with 183 and 179, respectively47,48 (Fig. 3b and Supplementary Table 7). Other highly associated DORCs included stromal markers DLK1, CXCL14 and COL6A3, as well as maternally expressed imprinted genes MEG3 and MEG8, the abnormal expression of which has been implicated, for both, in TB dysfunction49,50 (Fig. 3c, Extended Data Fig. 3a, Supplementary Fig. 6 and Supplementary Table 7). Many previously mentioned tumor invasion genes, immunomodulation-related genes and snATAC-seq-identified genes were found among DORCs (Fig. 3b,c, Extended Data Fig. 3a and Supplementary Table 7). We next calculated DORC scores across cell types, defined by the normalized sum of counts in all significantly correlated peaks per gene for all cells, and found that DORCs were highly cell type specific (Fig. 3c, Extended Data Fig. 3a,b and Supplementary Table 7). We identified TB-specific DORCs potentially vital to TB function such as KANK1 (tumor-suppressing YAP1 regulator), PARD6B (cell polarity regulator involved in tumor growth), FOXI3 (metastasis related) and MYCN/MYCNUT (growth promoter countered by TP63 (expressed in vCTBs), which might restrict MYCNUT/MYCN expression to EVTs)26,51–54 (Fig. 3c, Extended Data Fig. 3a, Supplementary Fig. 6 and Supplementary Table 7). In EVTs, we identified TEA domain proteins (TEAD) family binder VGLL3; EVTs are known to express high levels of its target IGFBP3, indicating that a VGLL3–TEAD1 axis might regulate insulin-like growth factor (IGF) at the maternal–fetal interface55,56 (Extended Data Fig. 3a,b, Supplementary Fig. 6 and Supplementary Table 7).

Biological intricacies predominantly manifest in spatial cellular arrangements of niches, physiological gradients and cell–cell interactions, which can be elucidated through single-cell spatial technologies. To map the early placenta, we performed STARmap on snap-frozen human placenta tissue sections from four of the eight previously described donors: W7-1, W8-2, W9 and W11 (Methods). First, we utilized STARmap-ISS and measured the spatial expression of 1,001 genes, composed of both highly variable genes from paired snRNA-seq and snATAC-seq and manually selected markers, and imaged a representative region of sectioned tissue spanning across the entire placenta (Fig. 1f and Supplementary Table 11). Cell-type labels were projected on to spatially resolved cells and manually verified (Fig. 4a–e, Extended Data Fig. 4a–c and Methods). Out of 17 clusters, 12 were predicted with 16,179 cells (Fig. 4a, Extended Data Fig. 5 and Supplementary Table 1). In addition, we identified a new list of differentially expressed genes (DEGs) from each STARmap cluster to examine spatial expression patterns of discovered marker genes (Fig. 4f and Supplementary Table 12). To investigate additional manually selected genes across large whole-tissue sections, STARmap-ISH analyses on 48 landmark genes were performed on 3 samples (W7-2, W8-2 and W11) for a total of 868,920 cells (Fig. 1g, Supplementary Figs. 8 and 9 and Supplementary Tables 1 and 11). Classic placental morphology across all samples was visualized with marker genes (Fig. 4d,e and Extended Data Fig. 4b,c). STARmap-ISH enabled additional spatial verification of Multiome-identified cell clusters and identification of TB-expressed tumor and immunomodulation-associated genes (Fig. 4e and Extended Data Fig. 6). STARmap-ISS and STARmap-ISH exhibited similar gene expression patterns, reinforcing ISS results.Fig. 4Spatial transcriptomic landscape of the early human placenta.a, UMAP representation of STARmap-ISS-derived clusters of the placenta. QC metrics can be found in Supplementary Table 1. STARmap genes are detailed in Supplementary Table 11. b, UMAP showing all cells of STARmap-ISS (1,001 genes).

Placentation is a complex and tightly regulated process that remains poorly understood, especially during the critical first trimester of pregnancy. In the present study, we used snRNA-seq and snATAC-seq to conduct an expansive investigation into previously uncharacterized, spatially resolved gene expression and CRE-linked epigenetic programs probably involved in continuous placental evolution and adaptation, including tumor invasion and immune evasion mechanisms. To spatially verify our findings, we adopted a unique strategy through the use of three different spatial technologies (Slide-tags, STARmap-ISS and STARmap-ISH) and imputed 33,357 gene multiomic profiles. We found that our three-pronged approach allowed the complementary techniques to synergistically account for each of their respective limitations, enhancing the breadth of our mapping and the reliability of our dataset while striking a time- and cost-effective balance between throughput and resolution.

tiomic profiles. We found that our three-pronged approach allowed the complementary techniques to synergistically account for each of their respective limitations, enhancing the breadth of our mapping and the reliability of our dataset while striking a time- and cost-effective balance between throughput and resolution. Recently, invaluable efforts have helped spatially deconstruct the human placenta6,7. Greenbaum et al.6 investigated how intravascular and perivascular EVTs might promote an environment conducive to spiral artery remodeling through 37-plex antibody panels and whole spatial transcriptomics. Arutyunyan et al.7 performed single-cell multiomics on frozen implantation-site tissue blocks, providing a detailed look at endovascular, interstitial and placental bed EVT invasion deep into maternal uterine layers with 55-µm 10x Visium spots composed of 1–10 cells for spatial analysis. Distinct from these studies, we focused on the villous core of the placenta (FIBs, HBCs and endos) and nearby YB subtypes (earlier-stage EVTs, vCTBs, STBs). We added true single-cell and molecular spatial resolution and extensive characterization of TF motif accessibility, CRE–gene linkages and tumor and immune-related genes across placental cell types in both single-cell and spatial analyses, an aspect that has needed further investigation. Ultimately, computational methods for integrating different data modalities are still in their early stages of development, especially for spatially phenotyping cells. Comprehensively validating findings such as those in our study will require relevant, functional in vitro and in vivo models alongside perturbation assays, all of which necessitate further optimization.

different data modalities are still in their early stages of development, especially for spatially phenotyping cells. Comprehensively validating findings such as those in our study will require relevant, functional in vitro and in vivo models alongside perturbation assays, all of which necessitate further optimization. In sum, our study is particularly relevant for future studies seeking to utilize multimodal analyses integrating cell morphology, gene expression and epigenetic regulation. Our selection of various analytical packages was guided by the belief that they were each well equipped for their respective task; as the field evolves, there may emerge more effective combinations of packages for handling the complexities of multimodal single-cell data and further exploration in this direction is essential. Ultimately, our adaptable and extensively modular experimental and computational framework identifies numerous genes and pathways for future functional studies and can be easily utilized to chart comprehensive cell atlases for entire organs in various species and disease models. This will greatly facilitate research into areas such as development, evolution and disease.

First trimester placental tissue from weeks 6–7 (n = 4), weeks 8–9 (n = 5) and weeks 10–11 (n = 3) was obtained from legal pregnancy terminations. Utilization of tissues and experimental procedures were approved by the ethics boards of the Medical University of Vienna (no. 084/2009). Written informed consent was required from donating women and no compensation was offered. Placental tissue was processed within 2 h of collection. Placental trees were cut from the chorionic plate, placed into optimal cutting temperature (OCT)-containing plastic molds, frozen on dry ice and stored at −80 °C.

First trimester placental tissue from weeks 6–7 (n = 4), weeks 8–9 (n = 5) and weeks 10–11 (n = 3) was obtained from legal pregnancy terminations. Utilization of tissues and experimental procedures were approved by the ethics boards of the Medical University of Vienna (no. 084/2009). Written informed consent was required from donating women and no compensation was offered. Placental tissue was processed within 2 h of collection. Placental trees were cut from the chorionic plate, placed into optimal cutting temperature (OCT)-containing plastic molds, frozen on dry ice and stored at −80 °C. Placental cells were collected by four consecutive enzymatic digestion steps as described elsewhere with minor modifications10. Single placentae were washed with Mg2+/Ca2+-free Hanks balanced salt solution (HBSS, Gibco), cut from the chorionic membranes and minced further into small pieces (2–3 mm). Three digestions were performed in prewarmed HBSS containing 0.25% trypsin (Gibco) and 1.25 mg ml−1 of DNase I (Sigma-Aldrich) for 10, 15 and 15 min in a water bath at 37 °C. Digestions were stopped with 10% fetal bovine serum (FBS), filtered through 100-µm cell strainers, pelleted (1,500 rpm for 5 min at 4 °C), pooled and washed twice with HBSS. To remove cell debris, cells were loaded on top of Percoll gradients (10–70% (v:v)) and cells between 15% and 60% were collected. Meanwhile, the remaining tissue was further digested using 0.5 mg ml−1 of DNase I and 1 mg ml−1 of Collagenase IV (Sigma-Aldrich) in HBSS for 30 min in a shaking water bath, pelleted and washed twice with HBSS. All cells collected after Percoll purification and digestion step 4 were pooled, frozen in cell banker 2 (0.5–5 × 106 cells per ml; Zenoaq) and stored at −80 °C.

using 0.5 mg ml−1 of DNase I and 1 mg ml−1 of Collagenase IV (Sigma-Aldrich) in HBSS for 30 min in a shaking water bath, pelleted and washed twice with HBSS. All cells collected after Percoll purification and digestion step 4 were pooled, frozen in cell banker 2 (0.5–5 × 106 cells per ml; Zenoaq) and stored at −80 °C. For selecting vCTBs (digestion 2 and 3), STBs and EVTs (digestion 1) and villous core cells (digestion 4), cells were selected from the above-described digestion steps and, if required, further purified using phycoerythrin-labeled specific antibodies/anti-phycoerythrin MicroBeads (MACS Miltenyi Biotec): EVTs were separated from STBs by sorting with HLA-G-PE (Exbio, cat. no. PE1P-292-C100, 1:20). The vCTBs, isolated with digestion steps 2 and 3 were further purified using CD49f/ITGA6-PE (BioLegend, cat. no. 313612, 1:20). Cell populations were pelleted and stored at −80 °C.

coerythrin MicroBeads (MACS Miltenyi Biotec): EVTs were separated from STBs by sorting with HLA-G-PE (Exbio, cat. no. PE1P-292-C100, 1:20). The vCTBs, isolated with digestion steps 2 and 3 were further purified using CD49f/ITGA6-PE (BioLegend, cat. no. 313612, 1:20). Cell populations were pelleted and stored at −80 °C. TB organoids were established, cultivated and differentiated as published previously10. Briefly, isolated vCTBs (weeks 6–7) were resuspended with TB-ORG stemness medium containing advanced Dulbecco’s modified Eagle’s medium (DMEM)/F12 (Invitrogen), 10 mM Hepes, 1× B27 (Gibco), 1× insulin–transferrin–selenium–ethanolamine (ITS-X) (Gibco, cat. no. 51500056.ITS-X), 2 mM Glutamax (Gibco), 0.05 mg ml−1 of gentamicin (Gibco), 2 µM A8301 (Tocris), 50 ng ml−1 of recombinant human epidermal growth factor (rhEGF; R&D Systems) and 3 µM CHIR99021 (Tocris). For the first TB-ORG formation, 5 µM Rock inhibitor (Santa Cruz, cat. no. Y27632) was also added. Cell/medium suspension was mixed with Matrigel (growth factor reduced, Corning) and 40-µl drops were placed into the center of 24-well culture dishes. The medium was routinely changed after 2–4 d. Organoids were split every 6–8 d as described previously10,69. For EVT differentiation, TB-ORGs were split and incubated in EVT differentiation medium consisting of advanced DMEM/F12, 10 mM Hepes, 1× B27, 1× ITS-X, 2 mM Glutamax, 0.05 mg ml−1 of gentamicin, 5 µM A8301 and 50 ng ml−1 of rhEGF. After 5 d, TB-ORGs were incubated with advanced DMEM/F12 for 1 h and further incubated for another 5 d with either DIFF-1 or DIFF-3 differentiation conditions: DIFF-1 TB-ORGs were incubated with advanced DMEM/F12, 10 mM Hepes, 1× B27, 1× ITS-X, 2 mM Glutamax, 0.05 mg ml−1 of gentamicin and 5 µM A8301. DIFF-3 TB-ORGs were incubated with advanced DMEM/F12, 10 mM Hepes, 1× B27, 1× ITS-X, 2 mM Glutamax, 0.05 mg ml−1 of gentamicin and 5 ng ml−1 of recombinant human transforming growth factor (TGF)-β1 (Abcam, cat. no. ab50036). Afterwards, organoid-derived EVTs were collected. Briefly, after dissolving Matrigel with cell recovery solution (Corning) for 45 min at 4 °C, cells were further separated with TrypLE for 10 min at 37 °C, filtered and pelleted and EVTs were enriched using HLA-G-phycoerythrin/anti-phycoerythrin MicroBeads. Subsequently, EVTs were counted, and equal cell numbers (0.2 × 106 cells per pellet) were pelleted and stored at −80 °C.

ution (Corning) for 45 min at 4 °C, cells were further separated with TrypLE for 10 min at 37 °C, filtered and pelleted and EVTs were enriched using HLA-G-phycoerythrin/anti-phycoerythrin MicroBeads. Subsequently, EVTs were counted, and equal cell numbers (0.2 × 106 cells per pellet) were pelleted and stored at −80 °C. Isolated vCTBs (weeks 6–7 of gestation) were plated on to fibronectin-coated (20 µg ml−1, Millipore) 6-well dishes (0.8 × 106 cells per well) in TSC medium consisting of advanced DMEM/F12, 10 mM Hepes, 1× B27, 1× ITS-X (100×), 2 mM Glutamax, 0.05 mg ml−1 of gentamicin, 2 µM A8301, 50 ng ml−1 of rhEGF, 3 µM CHIR99021 and 5 µM Rock inhibitor. After the second passage, the cells were split (10–15% confluency) and transfected with ON-TARGETplus siRNAs against FOXP1 (L-004256-01-0005) or nontargeting siRNAs (D-001810-10-0020) using Lipofectamine RNAiMAX (Invitrogen) according to the manufacturer’s instructions. The medium containing fresh siRNA was changed after 2–3 d. After 8 d, the cells were harvested for qPCR analyses.

y) and transfected with ON-TARGETplus siRNAs against FOXP1 (L-004256-01-0005) or nontargeting siRNAs (D-001810-10-0020) using Lipofectamine RNAiMAX (Invitrogen) according to the manufacturer’s instructions. The medium containing fresh siRNA was changed after 2–3 d. After 8 d, the cells were harvested for qPCR analyses. RNA isolation (PeqGold Trifast; PegLab), reverse transcription (RevertAid H Minus Reverse Transcriptase, Thermo Fisher Scientific) and qPCR (7500 Fast Real-time PCR system, Applied Biosystems) were performed as instructed by the manufacturers. For detection of miRNAs, RNA was isolated from relevant placental cell-type populations derived from additional placental samples (n = 4 donors, weeks 7–9 of gestation) according to the All-Prep RNA/DNA/miRNA Universal Kit (QIAGEN, cat. no. 80224) instructions and transcribed into complementary DNA using the TaqMan Advanced miRNA cDNA Synthesis Kit (Thermo Fisher Scientific, cat. no. 28007). The following TaqMan Gene Expression Assays (ABO) were used: CGB (ABI, cat. no. Hs00361224_gH), ENDOU (ABI, cat. no. Hs_00195731_m1), CCNA2 (ABI, cat. no. Hs00996788_m1), TP63 (ABI, cat. no. Hs00978340), miRNA7973 (ABI, cat. no. 483150_mir) and miRNA23B (ABI, cat. no.480150_mir); signals were normalized to TATA-box-binding protein (TBP: ABI, cat. no. 4333769F) and 18S (ABI, cat. no. Hs03928985_g1), respectively.

ENDOU (ABI, cat. no. Hs_00195731_m1), CCNA2 (ABI, cat. no. Hs00996788_m1), TP63 (ABI, cat. no. Hs00978340), miRNA7973 (ABI, cat. no. 483150_mir) and miRNA23B (ABI, cat. no.480150_mir); signals were normalized to TATA-box-binding protein (TBP: ABI, cat. no. 4333769F) and 18S (ABI, cat. no. Hs03928985_g1), respectively. EVT pellets were directly lysed in boiling sodium dodecylsulfate (SDS) loading buffer (0.15 M Tris, pH 6.8, 5% SDS, 25% glycerol, 1.6 M β-mercaptoethanol, 1 mg ml−1 of Bromophenol Blue) and separated on SDS–polyacrylamide gel electrophoresis (PAGE) gels, transferred on to Hybond-P poly(vinylidene difluoride) (GE Healthcare) membranes, blocked with 5% nonfat dry milk in Tris-buffered saline (TBS) containing 0.1% Tween (TBS-T) for 1 h at room temperature (RT) and incubated overnight at 4 °C with anti-FOXP1 (R&D Systems, 1:1,000), anti-HLA-G (Santa Cruz, 1:1,000), anti-P63 (Cell Signaling, 1:1,000), anti-CGβ (DAKO, 1:500) and anti-glyceraldehyde 3-phosphate dehydrogenase (Cell Signaling, 1:1,000) diluted in 5% bovine serum albumin (BSA)/TBS-T. Afterwards, membranes were washed 3× with TBS-T and incubated with anti-mouse horseradish peroxidase (HRP; Cell Signaling, 1:10,000) and anti-rabbit HRP (Cell Signaling, 1:10,000) for 1 h. Subsequently, membranes were washed 3× with TBS-T. Signals were developed using WesternBright Chemiluminescence Substrate Quantum (Biozym) and visualized with a ChemiDoc Imaging System (BioRad) using Image Lab 6.0 software.

roxidase (HRP; Cell Signaling, 1:10,000) and anti-rabbit HRP (Cell Signaling, 1:10,000) for 1 h. Subsequently, membranes were washed 3× with TBS-T. Signals were developed using WesternBright Chemiluminescence Substrate Quantum (Biozym) and visualized with a ChemiDoc Imaging System (BioRad) using Image Lab 6.0 software. IF of paraffin sections was performed as described previously70. Sections (2.5 µm) of paraffin-embedded placental tissue were deparaffinized and rehydrated, followed by antigen retrieval with citrate buffer, pH 6 (Sigma-Aldrich), using a KOS microwave histostation (Milestone). Sections were clamped with cover plates into vertical staining stations (Shandon), incubated with blocking solution (5% normal goat serum/TBS-T) for 1 h and further incubated with the following primary antibodies overnight at 4 °C, diluted in 5% normal goat serum/TBS-T: anti-FOXP1 (R&D Systems, 1:200), Endou (Sigma-Aldrich, 1:250), anti-SMARCC1 (Santa Cruz, 1:500) and anti-vimentin (Abcam, 1:200). The next day, the sections were washed 3× with TBS-T, incubated with appropriate secondary antibodies (anti-mouse 488 (cat. no. A11011), anti-rabbit 568 (cat. no. A-21069) (Alexa, Molecular Probes, 1:1,000)) for 1 h and nuclei were stained with 1 µg ml−1 of DAPI (Roche). Sections were embedded using fluoromount G (Soubio), analyzed by fluorescence microscopy (Olympus BX50) and digitally photographed (CellP software, Olympus).

ouse 488 (cat. no. A11011), anti-rabbit 568 (cat. no. A-21069) (Alexa, Molecular Probes, 1:1,000)) for 1 h and nuclei were stained with 1 µg ml−1 of DAPI (Roche). Sections were embedded using fluoromount G (Soubio), analyzed by fluorescence microscopy (Olympus BX50) and digitally photographed (CellP software, Olympus). Sample single-cell suspensions were thawed in a 37 °C water bath before incremental addition of warm DMEM + 10% FBS. Cells were filtered through a 40-μm filter before resuspension in phosphate-buffered saline (PBS) + 0.04% BSA (Invitrogen). Next, nuclei were isolated by following the Nuclei Isolation for Single Cell Multiome ATAC + Gene Expression Sequencing protocol (10x Genomics). After nuclei isolation, nuclei were processed and single-cell libraries were prepared using the Chromium Next GEM Single Cell Multiome ATAC + Gene Expression protocol (10x Genomics), loading 6,000 nuclei per lane. Pooled libraries were then sequenced using NextSeq High Output Cartridge kits and a NextSeq 550 sequencer (Illumina). RNA libraries were sequenced using the configuration of: R1, 28 cycles; R2, 44 cycles; Index1, 10 cycles; and Index2, 10 cycles. ATAC libraries were sequenced using the configuration of: R1, 30 cycles; R2, 30 cycles; Index1, 8 cycles; Index2, 16 cycles (+8 dark cycles).

idge kits and a NextSeq 550 sequencer (Illumina). RNA libraries were sequenced using the configuration of: R1, 28 cycles; R2, 44 cycles; Index1, 10 cycles; and Index2, 10 cycles. ATAC libraries were sequenced using the configuration of: R1, 30 cycles; R2, 30 cycles; Index1, 8 cycles; Index2, 16 cycles (+8 dark cycles). Slide-tags was performed as described in ref. 8. Barcoded bead arrays (also known as pucks) were fabricated using beads with the following sequence: 5′-TTT-PC-GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTJJJJJJJJTCTTCAGCGTTCCCGAGAJJJJJJJNNNNNNNVVA30. Fresh frozen tissues were cryo-sectioned to 20 μm on a cryostat (Leica, cat. no. CM1950). A precooled 5.5-mm squircle, custom-made biopsy punch was used to isolate regions of interest from tissue sections. The punched tissue regions were then melted on to the barcoded bead array, placed on ice and 6–10 µl of dissociation buffer (82 mM Na2SO4, 30 mM K2SO4, 10 mM glucose, 10 mM Hepes, 5 mM MgCl2) was placed on top of the puck so that the buffer covered the whole puck. The puck was then placed under an ultraviolet (365 nm) light source (0.42 mW mm−2, Thorlabs, cat. nos. M365LP1-C5 and LEDD1B) for 30 s, to cleave the spatial barcodes. After photocleavage, the puck was incubated for 7.5 min and then placed into a 12-well plate (Corning, cat. no. 3512). Using a 200-µl pipette, 10× 200-µl aliquots of extraction buffer (Dissociation Buffer, 1% Kollidon VA64, 1% Triton X-100, 0.01% BSA, 666 units ml−1 of RNase inhibitor (Biosearch Technologies, cat. no. 30281-1)) were dispensed on to the puck for a total volume of 2 ml. Dispensed extraction buffer was triturated up and down on the puck 10–15× to release the tissue. This step was repeated until the tissue was completely removed from the puck. The puck was removed and mechanical dissociation of the supernatant was performed using a 1-ml pipette and 20–30 triturations. Dissociated nuclei were removed from the well and the well was rinsed twice with 1 ml of wash buffer (82 mM Na2SO4, 30 mM K2SO4, 10 mM glucose, 10 mM Hepes, 5 mM MgCl2, 50 µl of RNase inhibitor (Biosearch Technologies, cat. no. 30281-1)) which was added to the nuclei suspension. Wash buffer was added to the tube to a final volume of 20 ml. This 20 ml was mixed and divided equally into another 50-ml falcon tube. Nuclei were spun in a precooled swinging bucket centrifuge at 600g for 10 min at 4 °C. After centrifugation, 19.5 ml of supernatant was removed, leaving 500 µl in each tube. The pellet was resuspended and pooled.

e tube to a final volume of 20 ml. This 20 ml was mixed and divided equally into another 50-ml falcon tube. Nuclei were spun in a precooled swinging bucket centrifuge at 600g for 10 min at 4 °C. After centrifugation, 19.5 ml of supernatant was removed, leaving 500 µl in each tube. The pellet was resuspended and pooled. This pooled suspension was then filtered using a precooled 40-µm cell strainer (Corning, cat. no. 431750). DAPI (Thermo Fisher Scientific, cat. no. 62248) was added to the filtered solution at a 1:1,000 dilution and incubated for 5–7 min at 4 °C. This was then centrifuged at 200g for 10 min at 4 °C. The supernatant was removed, leaving ~10 µl of pellet. The pellet was resuspended and nuclei were counted manually using a C-Chip Fuchs-Rosenthal disposable hemocytometer (INCYTO, cat. no. DHC-F01-5).

solution at a 1:1,000 dilution and incubated for 5–7 min at 4 °C. This was then centrifuged at 200g for 10 min at 4 °C. The supernatant was removed, leaving ~10 µl of pellet. The pellet was resuspended and nuclei were counted manually using a C-Chip Fuchs-Rosenthal disposable hemocytometer (INCYTO, cat. no. DHC-F01-5). Counted nuclei were loaded into the 10x Genomics Chromium controller using the Chromium Next GEM Single Cell Multiome ATAC + Gene Expression Reagent Bundle (10x Genomics, cat. no. PN-1000283). The Chromium Next GEM Single Cell Multiome ATAC + Gene Expression CG000338 Rev F user guide was used according to the manufacturer’s recommendations with modifications. During step 4.1, 1 µl of 0.329 µM spike-in primer (5′-GTGACTGGAGTTCAGACGT-3′) was added. For spatial barcode libraries, a customized PCR protocol was used: 5 µl of cleaned supernatant from step 4.3, 50 µl of NEBNext High-Fidelity 2× PCR Master Mix (NEB, cat. no. M0541S), 2.5 µl of 10 μM STAG_iP7_x oligo (5′-CAAGCAGAAGACGGCATACGAGATNNNNNNNNNNGTGACTGGAGTTCAGACGT*G*T-3′), 2.5 µl of 10 μM P5-STAG_ip5_x oligo (5′-AATGATACGGCGACCACCGAGATCTACACNNNNNNNNNNACACTCTTTCCCTACACGACGC*T*C-3′), 40 µl of UltraPure DNase/RNase-Free Distilled Water (Invitrogen, cat. no. 10977015). In this sample, 15 PCR cycles were performed according to the protocol used in the Chromium Next GEM Single Cell 3′ Reagent Kits v.3.1 (Dual Index) with Feature Barcode technology for Cell Surface Protein CG000317 Rev C user guide step 4.1.

of UltraPure DNase/RNase-Free Distilled Water (Invitrogen, cat. no. 10977015). In this sample, 15 PCR cycles were performed according to the protocol used in the Chromium Next GEM Single Cell 3′ Reagent Kits v.3.1 (Dual Index) with Feature Barcode technology for Cell Surface Protein CG000317 Rev C user guide step 4.1. Gene expression and spatial barcode libraries were sequenced together on an Illumina Nextseq 1000/2000 instrument using a p2 100-cycle kit (Illumina, cat. no. 20046811). Spatial barcodes were also further sequenced on their own lane of an Illumina Nextseq 1000/2000 instrument using a p2 100-cycle kit. ATAC libraries were also sequenced together using a p2 100-cycle kit. A glass-bottomed 12-well culture plate was pretreated with γ-methacryloxypropyltrimethoxysilane. The plate was further coated with poly(d-lysine) solution. OCT-embedded placenta tissues were cut into 10-µm, sections then fixed with 4% paraformaldehyde in PBS for 15 min and permeabilized by prechilled methanol at −20 °C for 20 min before hybridization.

plate was pretreated with γ-methacryloxypropyltrimethoxysilane. The plate was further coated with poly(d-lysine) solution. OCT-embedded placenta tissues were cut into 10-µm, sections then fixed with 4% paraformaldehyde in PBS for 15 min and permeabilized by prechilled methanol at −20 °C for 20 min before hybridization. STARmap-ISS 1,001 genes were chosen from a combination of highly variable genes elicited from our multiomic dataset and canonical markers, whereas STARmap-ISH genes were chosen from canonical markers and tumor/DNA damage-repair genes (Supplementary Table 9). SNAIL probes were designed and the probe library was constructed according to ref. 9. The probes were dissolved in ultrapure RNase-free water and pooled to the final concentration of 5 nM per probe. The probe mixture was heated at 40 °C for 15 min and then equilibrated to 37 °C. Tissue samples were removed from −20 °C tissues, equilibrated to RT and treated with 10 mM Tris, pH 7.5 for 10 min. The samples were then washed by PBS-TR (0.1% Tween-20, 0.1 U µl−1 of SUPERase•In in PBS) and incubated in hybridization buffer (2× saline–sodium citrate (SSC), 10% formamide, 1% Tween-20, 20 mM ribonucleoside vanadyl complex, 0.1 mg ml−1 of yeast transfer RNA, 0.1 U ul−1 of SUPERase•In, SNAIL probes with 5 nM per probe) at 40 °C with gentle shaking for 48 h. The samples were washed by PBS-TR twice for 20 min at 37 °C and then by 4× SSC in PBS-TR for 20 min at 37 °C, following a rinse by PBS-TR at RT. The SNAIL padlock probes annealed to the samples were ligated by incubating with T4 DNA ligation mixture (1:10 dilution of T4 DNA ligase, 0.2 mg ml−1 of BSA, 0.5 U ul−1 of SUPERase•In) for 2 h at RT with gentle agitation, followed by a 5-min wash of PBS-TR twice. The samples were then incubated in the RCA mixture (1:50 dilution of Phi29 DNA polymerase, 250 µM dNTP, 20 µM 5-(3-aminoallyl)-dUTP, 0.2 mg ml−1 of BSA, 0.2 U µl−1 of SUPERase•In) for 2 h at 30 °C with gentle agitation, followed by a 5-min wash of PBS-TR twice. Subsequently, the samples were treated with 25 mM acrylic acid N-hydroxysuccinimide ester for 2 h at RT with agitation, rinsed by PBS-T (0.1% Tween-20 in PBS) once and incubated with monomer buffer (4% acrylamide, 0.2% bis-acrylamide, 2× SSC in H2O) for 15 min for polymerization pretreatment.

-min wash of PBS-TR twice. Subsequently, the samples were treated with 25 mM acrylic acid N-hydroxysuccinimide ester for 2 h at RT with agitation, rinsed by PBS-T (0.1% Tween-20 in PBS) once and incubated with monomer buffer (4% acrylamide, 0.2% bis-acrylamide, 2× SSC in H2O) for 15 min for polymerization pretreatment. The buffer was removed and 30 µl of monomer mixture (0.1% ammonium persulfate, 0.1% tetramethylethylenediamine in monomer buffer) was directly added to the center of each sample, which was immediately covered by a coverslip (no. 2 coverslip was coated with Gel-Slick Solution according to the manufacturer’s instructions) and allowed to polymerize in ambiance for 1 h. The tissue–gel hybrid was washed with PBS-T twice and cleared by proteinase K digestion mixture (50 mM Tris, pH 7.5, 100 mM NaCl, 1% SDS, 0.2 mg ml−1 of proteinase K in H2O) at 37 °C overnight. On the next day, the samples were treated with a dephosphorylation mixture (1:100 dilution of shrimp alkaline phosphatase, 0.2 mg ml−1 of BSA, 1:10 dilution of CutSmart buffer in H2O) and rinsed with PBS-T. For ISH detection, the 19-nt fluorescent oligo complementary to DNA amplicon was diluted at 100 nM in 1× SSC dissolved in PBS-T and samples incubated at RT for 30 min, then washed by PBS-T 3× for 5 min each before imaging.

The buffer was removed and 30 µl of monomer mixture (0.1% ammonium persulfate, 0.1% tetramethylethylenediamine in monomer buffer) was directly added to the center of each sample, which was immediately covered by a coverslip (no. 2 coverslip was coated with Gel-Slick Solution according to the manufacturer’s instructions) and allowed to polymerize in ambiance for 1 h. The tissue–gel hybrid was washed with PBS-T twice and cleared by proteinase K digestion mixture (50 mM Tris, pH 7.5, 100 mM NaCl, 1% SDS, 0.2 mg ml−1 of proteinase K in H2O) at 37 °C overnight. On the next day, the samples were treated with a dephosphorylation mixture (1:100 dilution of shrimp alkaline phosphatase, 0.2 mg ml−1 of BSA, 1:10 dilution of CutSmart buffer in H2O) and rinsed with PBS-T. For ISH detection, the 19-nt fluorescent oligo complementary to DNA amplicon was diluted at 100 nM in 1× SSC dissolved in PBS-T and samples incubated at RT for 30 min, then washed by PBS-T 3× for 5 min each before imaging. For ISS detection, it was conducted according to ref. 9. Each cycle of sequencing started with the treatment of stripping buffer (60% formamide, 0.1% Triton X-100) for 5 min and triple washing with PBS-T for 5 min. The samples were then incubated with a sequencing mixture (0.2 mg ml−1 of BSA, 10 µM reading probe, 5 µM fluorescent decoding probe, 1:25 dilution of T4 DNA ligase) for 3 h at RT. Subsequently, the samples were triple washed with washing and imaging buffers (2× SSC, 10% formamide) for 10 min before proceeding to imaging. DAPI staining was performed after cycle 6 of imaging for cell segmentation. Images were acquired by Leica Stellaris 5 confocal microscope with a 405 diode, white light laser and ×40 oil-immersed objective (numerical aperture 1.3).

hing and imaging buffers (2× SSC, 10% formamide) for 10 min before proceeding to imaging. DAPI staining was performed after cycle 6 of imaging for cell segmentation. Images were acquired by Leica Stellaris 5 confocal microscope with a 405 diode, white light laser and ×40 oil-immersed objective (numerical aperture 1.3). We sought to obtain a collection of placental samples representing the first trimester of pregnancy after major structures have been formed (weeks 6–11), as well as across multiple different individuals (n = 9 donors); we felt that this number of donors would be sufficient given the scarce and precious nature of first trimester tissue. For follow-up experiments, we analyzed at least three independent biological sample sets. As the purpose of the present study was to uncover broad genomic underpinnings of the human first trimester placenta, randomization and blinding were not relevant to the present study. Sequencing reads were processed with Cellranger Arc in Cumulus on the Terra platform71 (https://app.terra.bio) using the Human GRCh38 sequences (GENCODE v.32/Ensembl 98), Cellranger Arc reference 2.0.0. For the integrated analysis of snATAC-seq and RNA-seq data, the ArchR package was employed72. The human reference genome hg38 was used for aligning ATAC-seq reads. Arrow files were generated with the addGeneScoreMat parameter set to TRUE to calculate gene scores for each cell.

sembl 98), Cellranger Arc reference 2.0.0. For the integrated analysis of snATAC-seq and RNA-seq data, the ArchR package was employed72. The human reference genome hg38 was used for aligning ATAC-seq reads. Arrow files were generated with the addGeneScoreMat parameter set to TRUE to calculate gene scores for each cell. Initial QC was applied to RNA-seq data, filtering out cells with <500 or >6,000 detected genes. The filtered RNA-seq data were integrated into the ArchR project as a gene expression matrix while excluding mitochondrial chromosome reads. Various QC metrics, including transcription start site (TSS) enrichment and total number of fragments (nFrags), were evaluated and visualized using ridges and violin plots. We utilized the ArchR function addDoubletScores to compute doublet scores and subsequently eliminated them. To address batch effects, the Harmony algorithm was applied to latent semantic indexing (LSI), reduced-dimensionality datasets for both RNA- and ATAC-seq73. Descriptive statistics, including mean and median values for the number of fragments, number of unique molecular identifiers (UMIs) and number of detected genes, were calculated. Subsequently, cells were clustered based on a combined reduced-dimensionality dataset of the harmonized RNA-seq and ATAC-seq data, employing a clustering resolution of 0.2. Uniform Manifold Approximation and Projection (UMAP) was used for visualization, with cluster labels applied based on known marker genes.

lculated. Subsequently, cells were clustered based on a combined reduced-dimensionality dataset of the harmonized RNA-seq and ATAC-seq data, employing a clustering resolution of 0.2. Uniform Manifold Approximation and Projection (UMAP) was used for visualization, with cluster labels applied based on known marker genes. In addition to the primary analysis, we conducted a comparative study using principal component analysis (PCA) on the RNA modality instead of LSI to support the robustness of our clustering. First, we applied batch correction with Harmony and then extracted the top 30 principal components (PCs). Both the PCA results and the ATAC–LSI data were scaled. After this, we integrated the modalities, computed the nearest neighbors and identified Louvain clusters with a resolution of 0.6, aiming to match the number of clusters observed in the LSI analysis. The results of this clustering approach, in comparison to the LSI–LSI clusters, are presented in Supplementary Fig. 11. Myeloid cells and EVTs were annotated based on specific marker transcripts. These cells were then filtered and included in new ArchR objects using the BiocGenerics package74. To identify subclusters within these cell populations, we used Seurat’s FindClusters function with the original Louvain clustering algorithm on the existing low-dimensional representation with resolution parameter set to 0.1. UMAP embeddings were then calculated for this subproject using the same low-dimensional data, applying a minimum distance of 0.8.

e cell populations, we used Seurat’s FindClusters function with the original Louvain clustering algorithm on the existing low-dimensional representation with resolution parameter set to 0.1. UMAP embeddings were then calculated for this subproject using the same low-dimensional data, applying a minimum distance of 0.8. To find marker genes for each subcluster, we performed Wilcoxon’s rank-sum tests, adjusting for confounding variables such as TSS enrichment and log(transformed fragment counts) (log10(nFrags)). We kept markers with a false discovery rate (FDR) of 0.1 or lower and a log2(fold-change) (log2(FC)) of 1.25 or higher. For manual annotation, we examined the top 20 genes for each subcluster based on their FDR values. This analysis identified three unique subclusters within the EVT group and four within the myeloid cells. These newly annotated subclusters were then integrated back into the original ArchR project. To distinguish cells with maternal or fetal origin, we used Freemuxlet, a genotype-free demultiplexing pipeline75. Using the 1000 Genome project variant call sites, we identified the cell origins based on the SNPs detected76. We calculated DEGs using the getMarkerFeatures function with bias = c(‘TSSEnrichment’, ‘log10(nFrags)’) and testMethod = ‘wilcoxon’ as the parameters. Gene scores were calculated with the addGeneScoreMatrix function with default parameters.

To distinguish cells with maternal or fetal origin, we used Freemuxlet, a genotype-free demultiplexing pipeline75. Using the 1000 Genome project variant call sites, we identified the cell origins based on the SNPs detected76. We calculated DEGs using the getMarkerFeatures function with bias = c(‘TSSEnrichment’, ‘log10(nFrags)’) and testMethod = ‘wilcoxon’ as the parameters. Gene scores were calculated with the addGeneScoreMatrix function with default parameters. To investigate intercluster regulatory differences, we leveraged the ArchR package and specifically employed its getMarkerFeatures function to identify marker features. In this analysis, one targeted cellular cluster was compared against a set of related clusters, which were used as the background. Differential accessibility was computed using Wilcoxon’s rank-sum test while accounting for biases such as sequencing depth and the number of ATAC-seq fragments. For annotating these marker features with TF motifs, the addMotifAnnotations function was utilized, drawing from the cis-BP database. Motif enrichment in these annotated features was subsequently performed using the peakAnnoEnrichment function. We selected motifs showing significant enrichment by applying a FDR cutoff of 0.1 and an absolute log2(FC) threshold of 0.5.

fs, the addMotifAnnotations function was utilized, drawing from the cis-BP database. Motif enrichment in these annotated features was subsequently performed using the peakAnnoEnrichment function. We selected motifs showing significant enrichment by applying a FDR cutoff of 0.1 and an absolute log2(FC) threshold of 0.5. We employed ArchR’s implementation of ChromVAR to assess per-cell TF activity. Initially, motif annotations from the cis-BP database were added to the ArchR project. We then generated a background peak set for bias correction (addBgdPeaks). A motif-focused deviations matrix was computed (addDeviationsMatrix), followed by extraction and visualization of variable motif deviations (getVarDeviations). The top 25 motifs were further analyzed and mapped to cell clusters and UMAP embeddings. To investigate the motifs’ relationship with gene expression, we visualized them in the gene score and gene expression matrices on the UMAP embeddings. After ChromVAR analysis, we employed ArchR to further pinpoint positive transcription factor regulators. We used the Motif Matrix within our existing ArchR project, grouped it by cell clusters and calculated a maxDelta metric for each motif to measure activity variations. Next, we correlated two key matrices: the Gene Score Matrix and the Gene Expression Matrix, each against the Motif Matrix. The correlation was executed using ArchR’s correlateMatrices function with reduced ‘LSI_Combined’ dimensions.

After ChromVAR analysis, we employed ArchR to further pinpoint positive transcription factor regulators. We used the Motif Matrix within our existing ArchR project, grouped it by cell clusters and calculated a maxDelta metric for each motif to measure activity variations. Next, we correlated two key matrices: the Gene Score Matrix and the Gene Expression Matrix, each against the Motif Matrix. The correlation was executed using ArchR’s correlateMatrices function with reduced ‘LSI_Combined’ dimensions. Motifs were then ranked and filtered based on criteria including correlation values, adjusted P values and the maxDelta metric. Those meeting the criteria were labeled as potential positive regulators. Expanding on co-accessibility analyses, we implemented ArchR’s ‘peak-to-gene links’ feature to correlate peak accessibility with gene expression using integrated snRNA-seq data. We set a maximum distance of 1 million base-pairs (1 M bp) for linkage and used the Gene Expression Matrix for the correlation. We ran addPeak2GeneLinks with the ‘LSI_Combined’ reduced dimensions to add these peak–gene links to the existing ArchR project. The links were then retrieved using getPeak2GeneLinks, applying a correlation cutoff of 0.45 and a resolution of 10. To visualize these peak–gene associations, we selected marker genes such as NOTCH1, CDX2 and ELF5, among others, and generated browser tracks using plotBrowserTrack.

links to the existing ArchR project. The links were then retrieved using getPeak2GeneLinks, applying a correlation cutoff of 0.45 and a resolution of 10. To visualize these peak–gene associations, we selected marker genes such as NOTCH1, CDX2 and ELF5, among others, and generated browser tracks using plotBrowserTrack. To define DORCs, we followed Ma et al.’s approach46 and ranked genes by the number of significantly associated peaks (±50 kb around TSSs), where we used 10 peaks per gene as cutoffs. Then, we re-calculated peak–gene association by expanding the window to ±500 kb around the TSSs. To calculate DORC scores, we first normalized peak counts by the total number of unique fragments in peaks per cell. Then, we defined the DORC scores for a given gene as the sum of counts in all significantly correlated peaks per gene to obtain the cell × DORC score matrix. To calculate chromatin potential, we calculated the distance (Di,j) between the chromatin profile of a given cell (Catac,i) and the gene expression profile of each cell (Crna,i,j). The arrow length was defined by normalizing Di,j.

To define DORCs, we followed Ma et al.’s approach46 and ranked genes by the number of significantly associated peaks (±50 kb around TSSs), where we used 10 peaks per gene as cutoffs. Then, we re-calculated peak–gene association by expanding the window to ±500 kb around the TSSs. To calculate DORC scores, we first normalized peak counts by the total number of unique fragments in peaks per cell. Then, we defined the DORC scores for a given gene as the sum of counts in all significantly correlated peaks per gene to obtain the cell × DORC score matrix. To calculate chromatin potential, we calculated the distance (Di,j) between the chromatin profile of a given cell (Catac,i) and the gene expression profile of each cell (Crna,i,j). The arrow length was defined by normalizing Di,j. To identify initial states and terminal states, we used CellRank59 where chromatin potential was used instead of the RNA velocity. After computing the transition matrix with compute_transition_matrix, we carried out the Schur decomposition with the compute_schur function. With compute_macrostates we identified 21 macrostates. Functions predict_initial_states and predict_terminal_states identified two and three states, respectively. To identify lineage driver genes, we found genes with high correlation values with fate probabilities of the respective lineages using the compute_lineage_drivers function.

crostates we identified 21 macrostates. Functions predict_initial_states and predict_terminal_states identified two and three states, respectively. To identify lineage driver genes, we found genes with high correlation values with fate probabilities of the respective lineages using the compute_lineage_drivers function. We sought to use our multimodal dataset to analyze potential disease risks for nine UKBB pregnancy-related traits with full summary statistics and average N = 169,000 (refs. 77–80). For each placental cell type, we generated four different types of SNP annotations by combining the information on peak–gene links and cell type-specific gene expression from the 10× Multiome data. These SNP annotations include: (1) LDSC-SEG81 annotation comprising SNPs in a 100-kb window around genes specifically enriched in expression across cell types; (2) sc-linker (ABC + Roadmap)82 annotation comprising SNPs linked to cell type-specific genes using enhancer–gene links from Roadmap83 and ABC84 in placenta biosamples; (3) Multiome annotation comprising SNPs in peaks linked to any gene in a cell type using the ArchR72 multiome peak–gene linking approach; and (4) modified sc-linker (Multiome) annotation comprising SNPs linked to cell type-specific genes using the ArchR peak–gene linking method. We assessed disease heritability informativeness of these annotations using a stratified linkage disequilibrium (LD) score regression (S-LDSC85) framework, conditional on a set of 86 baseline (baseline-LD v.2.1) annotations comprising coding, conserved, broad epigenomic annotations from ENCODE86 and Roadmap Epigenomics83 and LD-related annotations. As secondary analyses, we also examined the enrichment of GWAS hits for pregnancy-related traits in the above four types of placental cell-type annotations and the MAGMA GSEA87 of these placental cell-type gene programs.

d, broad epigenomic annotations from ENCODE86 and Roadmap Epigenomics83 and LD-related annotations. As secondary analyses, we also examined the enrichment of GWAS hits for pregnancy-related traits in the above four types of placental cell-type annotations and the MAGMA GSEA87 of these placental cell-type gene programs. Sc-linker is a method for computing a disease heritability enrichment score for a set of genes representing a gene program82. Using sc-linker, gene programs were defined from differential expression analysis of snRNA-seq data and SNPs were first linked to genes using a combination of DNase, histone mark and 3D contact-driven, driver enhancer–gene links from Roadmap Epigenomics83 and ABC84. In the present study, we employed a modified sc-linker strategy, sc-linker (Multiome), that links variants to genes using the ArchR peak–gene linking strategy. This was tested for gene programs determined from cell-type gene programs73 computed from the RNA-seq component of the 10× Multiome data, as well as for the gene program consisting of all genes. In addition to the Multiome peak–gene linking strategy and Roadmap and ABC enhancer–gene linking strategy, we used a standard 100-kb window-based strategy as suggested in ref. 81. Then, an enrichment score was computed for the SNPs based on the heritability enrichment of the SNPs obtained from S-LDSC85. More specifically, for each gene set G, a set of probabilistic weights between 0 and 1 was constructed for each SNP based on the confidence of their influencing any gene in G. The sc-linker method computes heritability enrichment estimates \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E}_{G}= \% ({h}^{2}(G))/ \% {{\mathrm{SNP}}}(G)$$\end{document}EG=%(h2(G))/%SNP(G) for a gene set G and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E}_{{{\mathrm{ALL}}}}= \% ({h}^{2}({{\mathrm{ALL}}}))/ \% {{\mathrm{SNP}}}({{\mathrm{ALL}}})$$\end{document}EALL=%(h2(ALL))/%SNP(ALL) for a gene program representing all genes.

fonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E}_{{{\mathrm{ALL}}}}= \% ({h}^{2}({{\mathrm{ALL}}}))/ \% {{\mathrm{SNP}}}({{\mathrm{ALL}}})$$\end{document}EALL=%(h2(ALL))/%SNP(ALL) for a gene program representing all genes. Here %(h2(A)) corresponds to the fraction of heritability captured by variants linked to genes in program A and %SNP(A) corresponds to the fraction of variants out of all common and low-frequency 1000 Genomes76 variants that are linked to the genes in program A. Finally, the enrichment score for A was computed as the difference \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E}_{G}$$\end{document}EG − \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E}_{{{\mathrm{ALL}}}}$$\end{document}EALL, where subtracting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E}_{{{\mathrm{ALL}}}}$$\end{document}EALL controls for the baseline level of heritability enrichment for SNPs that influence any gene (as most SNPs do not influence any genes). P values were obtained for the null hypothesis using a block jackknife procedure85.

greek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E}_{{{\mathrm{ALL}}}}$$\end{document}EALL controls for the baseline level of heritability enrichment for SNPs that influence any gene (as most SNPs do not influence any genes). P values were obtained for the null hypothesis using a block jackknife procedure85. In addition to the sc-linker enrichment analysis, we performed disease-related enrichment of the gene programs, using the enrichment of GWAS hits, defined by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\% {{\mathrm{overlap}}}({{\mathrm{GWAS}}\; {\mathrm{hits}}})/ \% {{\mathrm{overlap}}}({{\mathrm{control}}\; {\mathrm{variants}}})$$\end{document}%overlap(GWAShits)/%overlap(controlvariants) where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\% {{\mathrm{overlap}}}(A)$$\end{document}%overlap(A) is defined by the fraction of variants linked to genes by Multiome peak–gene or other enhancer–gene or window-based linking strategies with respect to the set of variants in A; the set of control variants, in this case, was defined by the common and low-frequency variants (minor allele count ≥5) from the 1000 Genomes Project. We additionally performed MAGMA gene set level analysis for the gene programs identified from the RNA-seq component of the Multiome data, following the pipeline suggested in ref. 87.

ntrol variants, in this case, was defined by the common and low-frequency variants (minor allele count ≥5) from the 1000 Genomes Project. We additionally performed MAGMA gene set level analysis for the gene programs identified from the RNA-seq component of the Multiome data, following the pipeline suggested in ref. 87. We used Cell Ranger-arc v.2.0.2 mkfastq (10x Genomics) to generate demultiplexed FASTQ files from the raw sequencing reads. We aligned these reads to the human GRCh38 genome and quantified gene counts as UMIs using Cell Ranger-arc count (ATAC and RNA-seq) and Cell Ranger count (v.6.1.2) (10x Genomics). The union of the list of cell barcodes called as cells by each of these counts was then used to generate filtered gene expression matrices. Filtered gene expression matrices from each sample were further processed using Seurat (v.4.3.0)88 together with R (v.4.1.1). Each Cell Ranger-called cell was mapped to a spatial coordinate as previously described8 by using DBSCAN for spatial barcode noise removal, followed by allocation of spatial coordinates using the spatial barcode, UMI-weighted centroid for each cell. SCTransform (v.0.3.5) was used to facilitate normalization, identifying the top 3,000 highly variable genes for each sample. Dimensionality was condensed using PCA down to 50 dimensions for each sample. Subsequent batch effect removal in the PCA space was accomplished with Harmony (v.0.1.1)73. For visualization, we deployed UMAP89, employing the top 20 Harmony-adjusted PCs. Shared nearest neighbors were identified using the same PCs. Cluster detection leveraged the Louvain method via ‘FindClusters‘ with a set resolution of 1.5. Initial cluster assignment insights were obtained through Seurat’s Label Transfer, referencing our labeled snRNA data. Furthermore, clusters with mitochondrial gene percentages >50% were filtered out. Definitive clusters were corroborated via canonical marker gene expression as seen in Supplementary Fig. 5. The top 50 differentially expressed, cluster-specific marker genes were generated using the FindAllMarkers function and ordered by log2(FC).

e, clusters with mitochondrial gene percentages >50% were filtered out. Definitive clusters were corroborated via canonical marker gene expression as seen in Supplementary Fig. 5. The top 50 differentially expressed, cluster-specific marker genes were generated using the FindAllMarkers function and ordered by log2(FC). ATAC data analysis was conducted using ArchR72 (v.1.0.3), drawing on the sorted combined fragments matrix and utilizing EnsDb.Hsapiens.v.86 for annotations. Dimensionality reduction was achieved using iterative LSI, with Harmony rectifying any batch effects. Chromatin accessibility peaks were initially discerned using MACS2 (v.2.2.7.1)90, leading to the generation of a peak set consistent across the identified cell clusters. Unraveling the potential transcriptional regulatory networks, the cis-BP motif database was employed to annotate TF-binding motifs within these peaks. Distinct chromatin accessibility markers for cellular clusters were discerned, with biases like TSS enrichment and log(transformed fragment counts) in mind. A Wilcoxon’s test highlighted significant variations in peak accessibility across distinct cellular states. TF motifs and their activities were visualized and assessed at the cellular level using ChromVAR36 (v.1.16.0). Assessment of genes and motif deviation scores that display spatial autocorrelation was conducted using Moran’s I with the package spdep91 (v.1.2-7). P values were adjusted for multiple comparisons using the Benjamini–Hochberg method.

ATAC data analysis was conducted using ArchR72 (v.1.0.3), drawing on the sorted combined fragments matrix and utilizing EnsDb.Hsapiens.v.86 for annotations. Dimensionality reduction was achieved using iterative LSI, with Harmony rectifying any batch effects. Chromatin accessibility peaks were initially discerned using MACS2 (v.2.2.7.1)90, leading to the generation of a peak set consistent across the identified cell clusters. Unraveling the potential transcriptional regulatory networks, the cis-BP motif database was employed to annotate TF-binding motifs within these peaks. Distinct chromatin accessibility markers for cellular clusters were discerned, with biases like TSS enrichment and log(transformed fragment counts) in mind. A Wilcoxon’s test highlighted significant variations in peak accessibility across distinct cellular states. TF motifs and their activities were visualized and assessed at the cellular level using ChromVAR36 (v.1.16.0). Assessment of genes and motif deviation scores that display spatial autocorrelation was conducted using Moran’s I with the package spdep91 (v.1.2-7). P values were adjusted for multiple comparisons using the Benjamini–Hochberg method. To start, images for individual tiles were deconvoluted using Huygens (v.23.04) to enhance image signals and suppress background noises. Next, we proceeded with image registration for different sequencing rounds, employing the 3D Fourier transform implemented through the functions available in numpy.fft and Scipy92. Crosscorrelation between pairs of images at all translational offsets was computed and the position associated with the highest correlation coefficient was identified. This position was then utilized to translate image volumes to compensate for the offset. During this procedure, the first sequencing round was used as the reference, and subsequent rounds were registered to be aligned with it.

computed and the position associated with the highest correlation coefficient was identified. This position was then utilized to translate image volumes to compensate for the offset. During this procedure, the first sequencing round was used as the reference, and subsequent rounds were registered to be aligned with it. After registration, individual dots were identified in each color channel on the first round of sequencing. Then, starfish (v.0.2.2)93 was used to process sequencing signals and registered images were converted into starfish format for downstream analysis. Amplicon dots were identified by BlobDetector in starfish, where the signal of each dot was treated as a Gaussian kernel. BlobDetector convolved kernels of multiple sizes and picked the best fit for each spot. The minimum s.d. and the maximum s.d. of Gaussian kernels were set as 0.5 and 10.0 voxel volume, respectively, owing to the predominant spot sizes falling within this specific range. The number of kernels was set as 10. Spots exhibiting local maxima in the lowest 10% of intensity values were eliminated to remove low-quality signals. Then, identified spots were decoded using the PerRoundMaxChannel function to pick the channel with maximum signal intensity within 3 × 3 × 3 voxel3 volume search radius for each round. A channel vector was generated for each spot by selecting channels from all six rounds. This vector was subsequently translated into a gene barcode using a barcode codebook, wherein each channel is linked with a two-base sequence encoding. Spots with undetectable signals in certain rounds were excluded from the analysis, as were spots with gene barcodes not present in the barcode codebook.

all six rounds. This vector was subsequently translated into a gene barcode using a barcode codebook, wherein each channel is linked with a two-base sequence encoding. Spots with undetectable signals in certain rounds were excluded from the analysis, as were spots with gene barcodes not present in the barcode codebook. ClusterMap (v.0.0.1)94 was employed for 2D cell segmentation of identified transcripts. Given the extensive scale of the data, ClusterMap was executed tile by tile. Image-free segmentation mode was utilized to comprehensively capture transcripts in regions of both nuclei and cytoplasm, with ‘xy_radius’ set as 55 and ‘cell_num_threshold’ set as 10−4. Transcripts that were not assigned to any cells were removed for downstream analysis. Segmented cells from individual tiles were integrated for subsequent single-cell analysis using Scanpy (v.1.9.4)95. Cells with <80 transcripts were labeled as low quality and excluded. After cell filtering, gene expression raw counts were normalized to 10,000, log(transformed) and scaled. PCA was conducted to obtain the reduced dimensions. Harmony73 was employed to integrate cells from different tiles and the corrected PCs were used to construct a neighborhood graph and generate UMAPs. Leiden clustering was performed using the PCs with a resolution of 1.0.

malized to 10,000, log(transformed) and scaled. PCA was conducted to obtain the reduced dimensions. Harmony73 was employed to integrate cells from different tiles and the corrected PCs were used to construct a neighborhood graph and generate UMAPs. Leiden clustering was performed using the PCs with a resolution of 1.0. Single cells of STARmap were annotated by leveraging labels transferred from single-nucleus multiomics reference data using Seurat (v.4). First, integration anchors of reference data were identified using reciprocal PCA. Then, transfer anchors that connected reference single-nucleus data and STARmap data were identified using the ‘TransferData’ function with the top 30 dimensions. Last, labels from reference data were transferred to STARmap single-cell data using MapQuery. After label transfer, differential expression analysis was conducted utilizing Wilcoxon’s test of Scanpy function rank_genes_groups. P values were corrected using the Benjamini–Hochberg method. We conducted the imputation of unmeasured genes by learning from snRNA-seq data following the integration of STARmap and single-cell multiomics reference data, employing a strategy similar to that described in ref. 96.

After label transfer, differential expression analysis was conducted utilizing Wilcoxon’s test of Scanpy function rank_genes_groups. P values were corrected using the Benjamini–Hochberg method. We conducted the imputation of unmeasured genes by learning from snRNA-seq data following the integration of STARmap and single-cell multiomics reference data, employing a strategy similar to that described in ref. 96. Initially, we performed intermediate mapping to determine the optimal parameters for imputation. Specifically, we aligned the STARmap and single-cell multiomics reference data using Seurat (v.4), as detailed in the previous analysis. For each cell in the STARmap data, we computed its nearest neighbors within the single-cell multiomics dataset, representing a set of the most similar cells in the reference dataset. For each of the 1,001 genes that overlapped between STARmap and single-cell multiomics reference data, we imputed the expression level as the average expression within the neighborhoods of individual cells. We explored various neighborhood sizes, ranging from 5 to 400, by calculating Pearson’s correlation between the imputed expression levels and observed expression levels in STARmap data for individual genes. We then assessed the imputation performance using accumulated Pearson’s correlation values and selected the 100 nearest neighbors based on this score.

s, ranging from 5 to 400, by calculating Pearson’s correlation between the imputed expression levels and observed expression levels in STARmap data for individual genes. We then assessed the imputation performance using accumulated Pearson’s correlation values and selected the 100 nearest neighbors based on this score. Subsequently, we performed a final imputation. We calculated the nearest neighbors and their distances within the integration data, as previously described. Then, we determined each gene’s imputed expression levels as the weighted average of gene expression levels in the neighborhood, where the weights were proportional to the reciprocal of distances. The output is a matrix measuring 16,179 cells by 33,357 genes. Apart from imputing the expression levels, we applied the same strategy to impute gene activities and DORC scores. We leveraged the previously mentioned nearest neighbors and distances, calculated from the integrated data, to determine the weighted average of gene activities and DORC scores for neighboring cells with matching epigenetic profiles in the single-cell multiome reference dataset. The imputed DORC scores were subsequently normalized to 10,000 per cell. These imputed gene activities and DORC scores were integrated with the imputed expression profiles for a comprehensive analysis.

and DORC scores for neighboring cells with matching epigenetic profiles in the single-cell multiome reference dataset. The imputed DORC scores were subsequently normalized to 10,000 per cell. These imputed gene activities and DORC scores were integrated with the imputed expression profiles for a comprehensive analysis. We employed CellChat64 to infer cell–cell interactions. For this purpose, we used the expression profiles of 33,357 imputed genes as input for CellChat and determined overexpressed genes and interactions using the human interaction database within CellChat. After this, we calculated communication probabilities while applying a distance constraint with a maximum interaction length of ligands set to 10 μm. We set the trim fraction to 0.1 and established a minimum threshold of 10 cells in each group for qualified communication. Only significant interactions were considered for subsequent analysis. For ISH images, we initiated the process by detecting spots within the raw images and generating multiplexed images, encompassing multiple genes within a single tissue. The initial step involved aligning images from two rounds of acquisition using the 2D Fourier transform, a technique consistent with the image registration approach utilized for STARmap data, as previously described.

n the raw images and generating multiplexed images, encompassing multiple genes within a single tissue. The initial step involved aligning images from two rounds of acquisition using the 2D Fourier transform, a technique consistent with the image registration approach utilized for STARmap data, as previously described. Subsequently, we employed a method to identify regions with local maxima, which serve as signal spots for visualization. Specifically, we applied maximum and minimum filters from the scipy.ndimage.filters library to locate potential local maxima and minimal values within neighborhoods sized at 15 pixels. Regions where the difference between local maxima and minimal values exceeded a defined threshold of 100 were identified and marked as local maxima. Connected regions were delineated as a result of this identification process. All tissue samples used for the present study were obtained with written informed consent from all participants and approval from the Medical University of Vienna Ethics Committee. All human tissue samples and experiments were approved by the Institutional Review Board at Broad Institute of MIT and Harvard and Massachusetts General Hospital. All experiments followed all relevant guidelines and regulations. Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

First trimester placental tissue from weeks 6–7 (n = 4), weeks 8–9 (n = 5) and weeks 10–11 (n = 3) was obtained from legal pregnancy terminations. Utilization of tissues and experimental procedures were approved by the ethics boards of the Medical University of Vienna (no. 084/2009). Written informed consent was required from donating women and no compensation was offered. Placental tissue was processed within 2 h of collection.

Placental cells were collected by four consecutive enzymatic digestion steps as described elsewhere with minor modifications10. Single placentae were washed with Mg2+/Ca2+-free Hanks balanced salt solution (HBSS, Gibco), cut from the chorionic membranes and minced further into small pieces (2–3 mm). Three digestions were performed in prewarmed HBSS containing 0.25% trypsin (Gibco) and 1.25 mg ml−1 of DNase I (Sigma-Aldrich) for 10, 15 and 15 min in a water bath at 37 °C. Digestions were stopped with 10% fetal bovine serum (FBS), filtered through 100-µm cell strainers, pelleted (1,500 rpm for 5 min at 4 °C), pooled and washed twice with HBSS. To remove cell debris, cells were loaded on top of Percoll gradients (10–70% (v:v)) and cells between 15% and 60% were collected. Meanwhile, the remaining tissue was further digested using 0.5 mg ml−1 of DNase I and 1 mg ml−1 of Collagenase IV (Sigma-Aldrich) in HBSS for 30 min in a shaking water bath, pelleted and washed twice with HBSS. All cells collected after Percoll purification and digestion step 4 were pooled, frozen in cell banker 2 (0.5–5 × 106 cells per ml; Zenoaq) and stored at −80 °C.

For selecting vCTBs (digestion 2 and 3), STBs and EVTs (digestion 1) and villous core cells (digestion 4), cells were selected from the above-described digestion steps and, if required, further purified using phycoerythrin-labeled specific antibodies/anti-phycoerythrin MicroBeads (MACS Miltenyi Biotec): EVTs were separated from STBs by sorting with HLA-G-PE (Exbio, cat. no. PE1P-292-C100, 1:20). The vCTBs, isolated with digestion steps 2 and 3 were further purified using CD49f/ITGA6-PE (BioLegend, cat. no. 313612, 1:20). Cell populations were pelleted and stored at −80 °C.

TB organoids were established, cultivated and differentiated as published previously10. Briefly, isolated vCTBs (weeks 6–7) were resuspended with TB-ORG stemness medium containing advanced Dulbecco’s modified Eagle’s medium (DMEM)/F12 (Invitrogen), 10 mM Hepes, 1× B27 (Gibco), 1× insulin–transferrin–selenium–ethanolamine (ITS-X) (Gibco, cat. no. 51500056.ITS-X), 2 mM Glutamax (Gibco), 0.05 mg ml−1 of gentamicin (Gibco), 2 µM A8301 (Tocris), 50 ng ml−1 of recombinant human epidermal growth factor (rhEGF; R&D Systems) and 3 µM CHIR99021 (Tocris). For the first TB-ORG formation, 5 µM Rock inhibitor (Santa Cruz, cat. no. Y27632) was also added. Cell/medium suspension was mixed with Matrigel (growth factor reduced, Corning) and 40-µl drops were placed into the center of 24-well culture dishes. The medium was routinely changed after 2–4 d. Organoids were split every 6–8 d as described previously10,69. For EVT differentiation, TB-ORGs were split and incubated in EVT differentiation medium consisting of advanced DMEM/F12, 10 mM Hepes, 1× B27, 1× ITS-X, 2 mM Glutamax, 0.05 mg ml−1 of gentamicin, 5 µM A8301 and 50 ng ml−1 of rhEGF. After 5 d, TB-ORGs were incubated with advanced DMEM/F12 for 1 h and further incubated for another 5 d with either DIFF-1 or DIFF-3 differentiation conditions: DIFF-1 TB-ORGs were incubated with advanced DMEM/F12, 10 mM Hepes, 1× B27, 1× ITS-X, 2 mM Glutamax, 0.05 mg ml−1 of gentamicin and 5 µM A8301. DIFF-3 TB-ORGs were incubated with advanced DMEM/F12, 10 mM Hepes, 1× B27, 1× ITS-X, 2 mM Glutamax, 0.05 mg ml−1 of gentamicin and 5 ng ml−1 of recombinant human transforming growth factor (TGF)-β1 (Abcam, cat. no. ab50036). Afterwards, organoid-derived EVTs were collected. Briefly, after dissolving Matrigel with cell recovery solution (Corning) for 45 min at 4 °C, cells were further separated with TrypLE for 10 min at 37 °C, filtered and pelleted and EVTs were enriched using HLA-G-phycoerythrin/anti-phycoerythrin MicroBeads. Subsequently, EVTs were counted, and equal cell numbers (0.2 × 106 cells per pellet) were pelleted and stored at −80 °C.

Isolated vCTBs (weeks 6–7 of gestation) were plated on to fibronectin-coated (20 µg ml−1, Millipore) 6-well dishes (0.8 × 106 cells per well) in TSC medium consisting of advanced DMEM/F12, 10 mM Hepes, 1× B27, 1× ITS-X (100×), 2 mM Glutamax, 0.05 mg ml−1 of gentamicin, 2 µM A8301, 50 ng ml−1 of rhEGF, 3 µM CHIR99021 and 5 µM Rock inhibitor. After the second passage, the cells were split (10–15% confluency) and transfected with ON-TARGETplus siRNAs against FOXP1 (L-004256-01-0005) or nontargeting siRNAs (D-001810-10-0020) using Lipofectamine RNAiMAX (Invitrogen) according to the manufacturer’s instructions. The medium containing fresh siRNA was changed after 2–3 d. After 8 d, the cells were harvested for qPCR analyses.

RNA isolation (PeqGold Trifast; PegLab), reverse transcription (RevertAid H Minus Reverse Transcriptase, Thermo Fisher Scientific) and qPCR (7500 Fast Real-time PCR system, Applied Biosystems) were performed as instructed by the manufacturers. For detection of miRNAs, RNA was isolated from relevant placental cell-type populations derived from additional placental samples (n = 4 donors, weeks 7–9 of gestation) according to the All-Prep RNA/DNA/miRNA Universal Kit (QIAGEN, cat. no. 80224) instructions and transcribed into complementary DNA using the TaqMan Advanced miRNA cDNA Synthesis Kit (Thermo Fisher Scientific, cat. no. 28007). The following TaqMan Gene Expression Assays (ABO) were used: CGB (ABI, cat. no. Hs00361224_gH), ENDOU (ABI, cat. no. Hs_00195731_m1), CCNA2 (ABI, cat. no. Hs00996788_m1), TP63 (ABI, cat. no. Hs00978340), miRNA7973 (ABI, cat. no. 483150_mir) and miRNA23B (ABI, cat. no.480150_mir); signals were normalized to TATA-box-binding protein (TBP: ABI, cat. no. 4333769F) and 18S (ABI, cat. no. Hs03928985_g1), respectively.

EVT pellets were directly lysed in boiling sodium dodecylsulfate (SDS) loading buffer (0.15 M Tris, pH 6.8, 5% SDS, 25% glycerol, 1.6 M β-mercaptoethanol, 1 mg ml−1 of Bromophenol Blue) and separated on SDS–polyacrylamide gel electrophoresis (PAGE) gels, transferred on to Hybond-P poly(vinylidene difluoride) (GE Healthcare) membranes, blocked with 5% nonfat dry milk in Tris-buffered saline (TBS) containing 0.1% Tween (TBS-T) for 1 h at room temperature (RT) and incubated overnight at 4 °C with anti-FOXP1 (R&D Systems, 1:1,000), anti-HLA-G (Santa Cruz, 1:1,000), anti-P63 (Cell Signaling, 1:1,000), anti-CGβ (DAKO, 1:500) and anti-glyceraldehyde 3-phosphate dehydrogenase (Cell Signaling, 1:1,000) diluted in 5% bovine serum albumin (BSA)/TBS-T. Afterwards, membranes were washed 3× with TBS-T and incubated with anti-mouse horseradish peroxidase (HRP; Cell Signaling, 1:10,000) and anti-rabbit HRP (Cell Signaling, 1:10,000) for 1 h. Subsequently, membranes were washed 3× with TBS-T. Signals were developed using WesternBright Chemiluminescence Substrate Quantum (Biozym) and visualized with a ChemiDoc Imaging System (BioRad) using Image Lab 6.0 software.

IF of paraffin sections was performed as described previously70. Sections (2.5 µm) of paraffin-embedded placental tissue were deparaffinized and rehydrated, followed by antigen retrieval with citrate buffer, pH 6 (Sigma-Aldrich), using a KOS microwave histostation (Milestone). Sections were clamped with cover plates into vertical staining stations (Shandon), incubated with blocking solution (5% normal goat serum/TBS-T) for 1 h and further incubated with the following primary antibodies overnight at 4 °C, diluted in 5% normal goat serum/TBS-T: anti-FOXP1 (R&D Systems, 1:200), Endou (Sigma-Aldrich, 1:250), anti-SMARCC1 (Santa Cruz, 1:500) and anti-vimentin (Abcam, 1:200). The next day, the sections were washed 3× with TBS-T, incubated with appropriate secondary antibodies (anti-mouse 488 (cat. no. A11011), anti-rabbit 568 (cat. no. A-21069) (Alexa, Molecular Probes, 1:1,000)) for 1 h and nuclei were stained with 1 µg ml−1 of DAPI (Roche). Sections were embedded using fluoromount G (Soubio), analyzed by fluorescence microscopy (Olympus BX50) and digitally photographed (CellP software, Olympus).

Sample single-cell suspensions were thawed in a 37 °C water bath before incremental addition of warm DMEM + 10% FBS. Cells were filtered through a 40-μm filter before resuspension in phosphate-buffered saline (PBS) + 0.04% BSA (Invitrogen). Next, nuclei were isolated by following the Nuclei Isolation for Single Cell Multiome ATAC + Gene Expression Sequencing protocol (10x Genomics). After nuclei isolation, nuclei were processed and single-cell libraries were prepared using the Chromium Next GEM Single Cell Multiome ATAC + Gene Expression protocol (10x Genomics), loading 6,000 nuclei per lane. Pooled libraries were then sequenced using NextSeq High Output Cartridge kits and a NextSeq 550 sequencer (Illumina). RNA libraries were sequenced using the configuration of: R1, 28 cycles; R2, 44 cycles; Index1, 10 cycles; and Index2, 10 cycles. ATAC libraries were sequenced using the configuration of: R1, 30 cycles; R2, 30 cycles; Index1, 8 cycles; Index2, 16 cycles (+8 dark cycles).

Slide-tags was performed as described in ref. 8. Barcoded bead arrays (also known as pucks) were fabricated using beads with the following sequence: 5′-TTT-PC-GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTJJJJJJJJTCTTCAGCGTTCCCGAGAJJJJJJJNNNNNNNVVA30. Fresh frozen tissues were cryo-sectioned to 20 μm on a cryostat (Leica, cat. no. CM1950). A precooled 5.5-mm squircle, custom-made biopsy punch was used to isolate regions of interest from tissue sections. The punched tissue regions were then melted on to the barcoded bead array, placed on ice and 6–10 µl of dissociation buffer (82 mM Na2SO4, 30 mM K2SO4, 10 mM glucose, 10 mM Hepes, 5 mM MgCl2) was placed on top of the puck so that the buffer covered the whole puck. The puck was then placed under an ultraviolet (365 nm) light source (0.42 mW mm−2, Thorlabs, cat. nos. M365LP1-C5 and LEDD1B) for 30 s, to cleave the spatial barcodes. After photocleavage, the puck was incubated for 7.5 min and then placed into a 12-well plate (Corning, cat. no. 3512). Using a 200-µl pipette, 10× 200-µl aliquots of extraction buffer (Dissociation Buffer, 1% Kollidon VA64, 1% Triton X-100, 0.01% BSA, 666 units ml−1 of RNase inhibitor (Biosearch Technologies, cat. no. 30281-1)) were dispensed on to the puck for a total volume of 2 ml. Dispensed extraction buffer was triturated up and down on the puck 10–15× to release the tissue. This step was repeated until the tissue was completely removed from the puck. The puck was removed and mechanical dissociation of the supernatant was performed using a 1-ml pipette and 20–30 triturations. Dissociated nuclei were removed from the well and the well was rinsed twice with 1 ml of wash buffer (82 mM Na2SO4, 30 mM K2SO4, 10 mM glucose, 10 mM Hepes, 5 mM MgCl2, 50 µl of RNase inhibitor (Biosearch Technologies, cat. no. 30281-1)) which was added to the nuclei suspension. Wash buffer was added to the tube to a final volume of 20 ml. This 20 ml was mixed and divided equally into another 50-ml falcon tube. Nuclei were spun in a precooled swinging bucket centrifuge at 600g for 10 min at 4 °C. After centrifugation, 19.5 ml of supernatant was removed, leaving 500 µl in each tube. The pellet was resuspended and pooled.

of UltraPure DNase/RNase-Free Distilled Water (Invitrogen, cat. no. 10977015). In this sample, 15 PCR cycles were performed according to the protocol used in the Chromium Next GEM Single Cell 3′ Reagent Kits v.3.1 (Dual Index) with Feature Barcode technology for Cell Surface Protein CG000317 Rev C user guide step 4.1. Gene expression and spatial barcode libraries were sequenced together on an Illumina Nextseq 1000/2000 instrument using a p2 100-cycle kit (Illumina, cat. no. 20046811). Spatial barcodes were also further sequenced on their own lane of an Illumina Nextseq 1000/2000 instrument using a p2 100-cycle kit. ATAC libraries were also sequenced together using a p2 100-cycle kit.

A glass-bottomed 12-well culture plate was pretreated with γ-methacryloxypropyltrimethoxysilane. The plate was further coated with poly(d-lysine) solution. OCT-embedded placenta tissues were cut into 10-µm, sections then fixed with 4% paraformaldehyde in PBS for 15 min and permeabilized by prechilled methanol at −20 °C for 20 min before hybridization.

STARmap-ISS 1,001 genes were chosen from a combination of highly variable genes elicited from our multiomic dataset and canonical markers, whereas STARmap-ISH genes were chosen from canonical markers and tumor/DNA damage-repair genes (Supplementary Table 9). SNAIL probes were designed and the probe library was constructed according to ref. 9. The probes were dissolved in ultrapure RNase-free water and pooled to the final concentration of 5 nM per probe. The probe mixture was heated at 40 °C for 15 min and then equilibrated to 37 °C. Tissue samples were removed from −20 °C tissues, equilibrated to RT and treated with 10 mM Tris, pH 7.5 for 10 min. The samples were then washed by PBS-TR (0.1% Tween-20, 0.1 U µl−1 of SUPERase•In in PBS) and incubated in hybridization buffer (2× saline–sodium citrate (SSC), 10% formamide, 1% Tween-20, 20 mM ribonucleoside vanadyl complex, 0.1 mg ml−1 of yeast transfer RNA, 0.1 U ul−1 of SUPERase•In, SNAIL probes with 5 nM per probe) at 40 °C with gentle shaking for 48 h. The samples were washed by PBS-TR twice for 20 min at 37 °C and then by 4× SSC in PBS-TR for 20 min at 37 °C, following a rinse by PBS-TR at RT. The SNAIL padlock probes annealed to the samples were ligated by incubating with T4 DNA ligation mixture (1:10 dilution of T4 DNA ligase, 0.2 mg ml−1 of BSA, 0.5 U ul−1 of SUPERase•In) for 2 h at RT with gentle agitation, followed by a 5-min wash of PBS-TR twice. The samples were then incubated in the RCA mixture (1:50 dilution of Phi29 DNA polymerase, 250 µM dNTP, 20 µM 5-(3-aminoallyl)-dUTP, 0.2 mg ml−1 of BSA, 0.2 U µl−1 of SUPERase•In) for 2 h at 30 °C with gentle agitation, followed by a 5-min wash of PBS-TR twice. Subsequently, the samples were treated with 25 mM acrylic acid N-hydroxysuccinimide ester for 2 h at RT with agitation, rinsed by PBS-T (0.1% Tween-20 in PBS) once and incubated with monomer buffer (4% acrylamide, 0.2% bis-acrylamide, 2× SSC in H2O) for 15 min for polymerization pretreatment.

-min wash of PBS-TR twice. Subsequently, the samples were treated with 25 mM acrylic acid N-hydroxysuccinimide ester for 2 h at RT with agitation, rinsed by PBS-T (0.1% Tween-20 in PBS) once and incubated with monomer buffer (4% acrylamide, 0.2% bis-acrylamide, 2× SSC in H2O) for 15 min for polymerization pretreatment. The buffer was removed and 30 µl of monomer mixture (0.1% ammonium persulfate, 0.1% tetramethylethylenediamine in monomer buffer) was directly added to the center of each sample, which was immediately covered by a coverslip (no. 2 coverslip was coated with Gel-Slick Solution according to the manufacturer’s instructions) and allowed to polymerize in ambiance for 1 h. The tissue–gel hybrid was washed with PBS-T twice and cleared by proteinase K digestion mixture (50 mM Tris, pH 7.5, 100 mM NaCl, 1% SDS, 0.2 mg ml−1 of proteinase K in H2O) at 37 °C overnight. On the next day, the samples were treated with a dephosphorylation mixture (1:100 dilution of shrimp alkaline phosphatase, 0.2 mg ml−1 of BSA, 1:10 dilution of CutSmart buffer in H2O) and rinsed with PBS-T.

For ISH detection, the 19-nt fluorescent oligo complementary to DNA amplicon was diluted at 100 nM in 1× SSC dissolved in PBS-T and samples incubated at RT for 30 min, then washed by PBS-T 3× for 5 min each before imaging. For ISS detection, it was conducted according to ref. 9. Each cycle of sequencing started with the treatment of stripping buffer (60% formamide, 0.1% Triton X-100) for 5 min and triple washing with PBS-T for 5 min. The samples were then incubated with a sequencing mixture (0.2 mg ml−1 of BSA, 10 µM reading probe, 5 µM fluorescent decoding probe, 1:25 dilution of T4 DNA ligase) for 3 h at RT. Subsequently, the samples were triple washed with washing and imaging buffers (2× SSC, 10% formamide) for 10 min before proceeding to imaging. DAPI staining was performed after cycle 6 of imaging for cell segmentation. Images were acquired by Leica Stellaris 5 confocal microscope with a 405 diode, white light laser and ×40 oil-immersed objective (numerical aperture 1.3).

We sought to obtain a collection of placental samples representing the first trimester of pregnancy after major structures have been formed (weeks 6–11), as well as across multiple different individuals (n = 9 donors); we felt that this number of donors would be sufficient given the scarce and precious nature of first trimester tissue. For follow-up experiments, we analyzed at least three independent biological sample sets. As the purpose of the present study was to uncover broad genomic underpinnings of the human first trimester placenta, randomization and blinding were not relevant to the present study. Sequencing reads were processed with Cellranger Arc in Cumulus on the Terra platform71 (https://app.terra.bio) using the Human GRCh38 sequences (GENCODE v.32/Ensembl 98), Cellranger Arc reference 2.0.0. For the integrated analysis of snATAC-seq and RNA-seq data, the ArchR package was employed72. The human reference genome hg38 was used for aligning ATAC-seq reads. Arrow files were generated with the addGeneScoreMat parameter set to TRUE to calculate gene scores for each cell. Initial QC was applied to RNA-seq data, filtering out cells with <500 or >6,000 detected genes. The filtered RNA-seq data were integrated into the ArchR project as a gene expression matrix while excluding mitochondrial chromosome reads. Various QC metrics, including transcription start site (TSS) enrichment and total number of fragments (nFrags), were evaluated and visualized using ridges and violin plots.

genes. The filtered RNA-seq data were integrated into the ArchR project as a gene expression matrix while excluding mitochondrial chromosome reads. Various QC metrics, including transcription start site (TSS) enrichment and total number of fragments (nFrags), were evaluated and visualized using ridges and violin plots. We utilized the ArchR function addDoubletScores to compute doublet scores and subsequently eliminated them. To address batch effects, the Harmony algorithm was applied to latent semantic indexing (LSI), reduced-dimensionality datasets for both RNA- and ATAC-seq73. Descriptive statistics, including mean and median values for the number of fragments, number of unique molecular identifiers (UMIs) and number of detected genes, were calculated. Subsequently, cells were clustered based on a combined reduced-dimensionality dataset of the harmonized RNA-seq and ATAC-seq data, employing a clustering resolution of 0.2. Uniform Manifold Approximation and Projection (UMAP) was used for visualization, with cluster labels applied based on known marker genes.

n the raw images and generating multiplexed images, encompassing multiple genes within a single tissue. The initial step involved aligning images from two rounds of acquisition using the 2D Fourier transform, a technique consistent with the image registration approach utilized for STARmap data, as previously described. Subsequently, we employed a method to identify regions with local maxima, which serve as signal spots for visualization. Specifically, we applied maximum and minimum filters from the scipy.ndimage.filters library to locate potential local maxima and minimal values within neighborhoods sized at 15 pixels. Regions where the difference between local maxima and minimal values exceeded a defined threshold of 100 were identified and marked as local maxima. Connected regions were delineated as a result of this identification process.

We sought to obtain a collection of placental samples representing the first trimester of pregnancy after major structures have been formed (weeks 6–11), as well as across multiple different individuals (n = 9 donors); we felt that this number of donors would be sufficient given the scarce and precious nature of first trimester tissue. For follow-up experiments, we analyzed at least three independent biological sample sets. As the purpose of the present study was to uncover broad genomic underpinnings of the human first trimester placenta, randomization and blinding were not relevant to the present study.

Sequencing reads were processed with Cellranger Arc in Cumulus on the Terra platform71 (https://app.terra.bio) using the Human GRCh38 sequences (GENCODE v.32/Ensembl 98), Cellranger Arc reference 2.0.0. For the integrated analysis of snATAC-seq and RNA-seq data, the ArchR package was employed72. The human reference genome hg38 was used for aligning ATAC-seq reads. Arrow files were generated with the addGeneScoreMat parameter set to TRUE to calculate gene scores for each cell. Initial QC was applied to RNA-seq data, filtering out cells with <500 or >6,000 detected genes. The filtered RNA-seq data were integrated into the ArchR project as a gene expression matrix while excluding mitochondrial chromosome reads. Various QC metrics, including transcription start site (TSS) enrichment and total number of fragments (nFrags), were evaluated and visualized using ridges and violin plots. We utilized the ArchR function addDoubletScores to compute doublet scores and subsequently eliminated them.

Initial QC was applied to RNA-seq data, filtering out cells with <500 or >6,000 detected genes. The filtered RNA-seq data were integrated into the ArchR project as a gene expression matrix while excluding mitochondrial chromosome reads. Various QC metrics, including transcription start site (TSS) enrichment and total number of fragments (nFrags), were evaluated and visualized using ridges and violin plots. We utilized the ArchR function addDoubletScores to compute doublet scores and subsequently eliminated them. To address batch effects, the Harmony algorithm was applied to latent semantic indexing (LSI), reduced-dimensionality datasets for both RNA- and ATAC-seq73. Descriptive statistics, including mean and median values for the number of fragments, number of unique molecular identifiers (UMIs) and number of detected genes, were calculated. Subsequently, cells were clustered based on a combined reduced-dimensionality dataset of the harmonized RNA-seq and ATAC-seq data, employing a clustering resolution of 0.2. Uniform Manifold Approximation and Projection (UMAP) was used for visualization, with cluster labels applied based on known marker genes.

lculated. Subsequently, cells were clustered based on a combined reduced-dimensionality dataset of the harmonized RNA-seq and ATAC-seq data, employing a clustering resolution of 0.2. Uniform Manifold Approximation and Projection (UMAP) was used for visualization, with cluster labels applied based on known marker genes. In addition to the primary analysis, we conducted a comparative study using principal component analysis (PCA) on the RNA modality instead of LSI to support the robustness of our clustering. First, we applied batch correction with Harmony and then extracted the top 30 principal components (PCs). Both the PCA results and the ATAC–LSI data were scaled. After this, we integrated the modalities, computed the nearest neighbors and identified Louvain clusters with a resolution of 0.6, aiming to match the number of clusters observed in the LSI analysis. The results of this clustering approach, in comparison to the LSI–LSI clusters, are presented in Supplementary Fig. 11.

Myeloid cells and EVTs were annotated based on specific marker transcripts. These cells were then filtered and included in new ArchR objects using the BiocGenerics package74. To identify subclusters within these cell populations, we used Seurat’s FindClusters function with the original Louvain clustering algorithm on the existing low-dimensional representation with resolution parameter set to 0.1. UMAP embeddings were then calculated for this subproject using the same low-dimensional data, applying a minimum distance of 0.8. To find marker genes for each subcluster, we performed Wilcoxon’s rank-sum tests, adjusting for confounding variables such as TSS enrichment and log(transformed fragment counts) (log10(nFrags)). We kept markers with a false discovery rate (FDR) of 0.1 or lower and a log2(fold-change) (log2(FC)) of 1.25 or higher. For manual annotation, we examined the top 20 genes for each subcluster based on their FDR values. This analysis identified three unique subclusters within the EVT group and four within the myeloid cells. These newly annotated subclusters were then integrated back into the original ArchR project.

To distinguish cells with maternal or fetal origin, we used Freemuxlet, a genotype-free demultiplexing pipeline75. Using the 1000 Genome project variant call sites, we identified the cell origins based on the SNPs detected76.

We calculated DEGs using the getMarkerFeatures function with bias = c(‘TSSEnrichment’, ‘log10(nFrags)’) and testMethod = ‘wilcoxon’ as the parameters. Gene scores were calculated with the addGeneScoreMatrix function with default parameters.

To investigate intercluster regulatory differences, we leveraged the ArchR package and specifically employed its getMarkerFeatures function to identify marker features. In this analysis, one targeted cellular cluster was compared against a set of related clusters, which were used as the background. Differential accessibility was computed using Wilcoxon’s rank-sum test while accounting for biases such as sequencing depth and the number of ATAC-seq fragments. For annotating these marker features with TF motifs, the addMotifAnnotations function was utilized, drawing from the cis-BP database. Motif enrichment in these annotated features was subsequently performed using the peakAnnoEnrichment function. We selected motifs showing significant enrichment by applying a FDR cutoff of 0.1 and an absolute log2(FC) threshold of 0.5.

We employed ArchR’s implementation of ChromVAR to assess per-cell TF activity. Initially, motif annotations from the cis-BP database were added to the ArchR project. We then generated a background peak set for bias correction (addBgdPeaks). A motif-focused deviations matrix was computed (addDeviationsMatrix), followed by extraction and visualization of variable motif deviations (getVarDeviations). The top 25 motifs were further analyzed and mapped to cell clusters and UMAP embeddings. To investigate the motifs’ relationship with gene expression, we visualized them in the gene score and gene expression matrices on the UMAP embeddings. After ChromVAR analysis, we employed ArchR to further pinpoint positive transcription factor regulators. We used the Motif Matrix within our existing ArchR project, grouped it by cell clusters and calculated a maxDelta metric for each motif to measure activity variations. Next, we correlated two key matrices: the Gene Score Matrix and the Gene Expression Matrix, each against the Motif Matrix. The correlation was executed using ArchR’s correlateMatrices function with reduced ‘LSI_Combined’ dimensions. Motifs were then ranked and filtered based on criteria including correlation values, adjusted P values and the maxDelta metric. Those meeting the criteria were labeled as potential positive regulators.

Expanding on co-accessibility analyses, we implemented ArchR’s ‘peak-to-gene links’ feature to correlate peak accessibility with gene expression using integrated snRNA-seq data. We set a maximum distance of 1 million base-pairs (1 M bp) for linkage and used the Gene Expression Matrix for the correlation. We ran addPeak2GeneLinks with the ‘LSI_Combined’ reduced dimensions to add these peak–gene links to the existing ArchR project. The links were then retrieved using getPeak2GeneLinks, applying a correlation cutoff of 0.45 and a resolution of 10. To visualize these peak–gene associations, we selected marker genes such as NOTCH1, CDX2 and ELF5, among others, and generated browser tracks using plotBrowserTrack. To define DORCs, we followed Ma et al.’s approach46 and ranked genes by the number of significantly associated peaks (±50 kb around TSSs), where we used 10 peaks per gene as cutoffs. Then, we re-calculated peak–gene association by expanding the window to ±500 kb around the TSSs. To calculate DORC scores, we first normalized peak counts by the total number of unique fragments in peaks per cell. Then, we defined the DORC scores for a given gene as the sum of counts in all significantly correlated peaks per gene to obtain the cell × DORC score matrix. To calculate chromatin potential, we calculated the distance (Di,j) between the chromatin profile of a given cell (Catac,i) and the gene expression profile of each cell (Crna,i,j). The arrow length was defined by normalizing Di,j.

We sought to use our multimodal dataset to analyze potential disease risks for nine UKBB pregnancy-related traits with full summary statistics and average N = 169,000 (refs. 77–80). For each placental cell type, we generated four different types of SNP annotations by combining the information on peak–gene links and cell type-specific gene expression from the 10× Multiome data. These SNP annotations include: (1) LDSC-SEG81 annotation comprising SNPs in a 100-kb window around genes specifically enriched in expression across cell types; (2) sc-linker (ABC + Roadmap)82 annotation comprising SNPs linked to cell type-specific genes using enhancer–gene links from Roadmap83 and ABC84 in placenta biosamples; (3) Multiome annotation comprising SNPs in peaks linked to any gene in a cell type using the ArchR72 multiome peak–gene linking approach; and (4) modified sc-linker (Multiome) annotation comprising SNPs linked to cell type-specific genes using the ArchR peak–gene linking method. We assessed disease heritability informativeness of these annotations using a stratified linkage disequilibrium (LD) score regression (S-LDSC85) framework, conditional on a set of 86 baseline (baseline-LD v.2.1) annotations comprising coding, conserved, broad epigenomic annotations from ENCODE86 and Roadmap Epigenomics83 and LD-related annotations. As secondary analyses, we also examined the enrichment of GWAS hits for pregnancy-related traits in the above four types of placental cell-type annotations and the MAGMA GSEA87 of these placental cell-type gene programs.

We used Cell Ranger-arc v.2.0.2 mkfastq (10x Genomics) to generate demultiplexed FASTQ files from the raw sequencing reads. We aligned these reads to the human GRCh38 genome and quantified gene counts as UMIs using Cell Ranger-arc count (ATAC and RNA-seq) and Cell Ranger count (v.6.1.2) (10x Genomics). The union of the list of cell barcodes called as cells by each of these counts was then used to generate filtered gene expression matrices. Filtered gene expression matrices from each sample were further processed using Seurat (v.4.3.0)88 together with R (v.4.1.1). Each Cell Ranger-called cell was mapped to a spatial coordinate as previously described8 by using DBSCAN for spatial barcode noise removal, followed by allocation of spatial coordinates using the spatial barcode, UMI-weighted centroid for each cell. SCTransform (v.0.3.5) was used to facilitate normalization, identifying the top 3,000 highly variable genes for each sample. Dimensionality was condensed using PCA down to 50 dimensions for each sample. Subsequent batch effect removal in the PCA space was accomplished with Harmony (v.0.1.1)73. For visualization, we deployed UMAP89, employing the top 20 Harmony-adjusted PCs. Shared nearest neighbors were identified using the same PCs. Cluster detection leveraged the Louvain method via ‘FindClusters‘ with a set resolution of 1.5. Initial cluster assignment insights were obtained through Seurat’s Label Transfer, referencing our labeled snRNA data. Furthermore, clusters with mitochondrial gene percentages >50% were filtered out. Definitive clusters were corroborated via canonical marker gene expression as seen in Supplementary Fig. 5. The top 50 differentially expressed, cluster-specific marker genes were generated using the FindAllMarkers function and ordered by log2(FC).

ATAC data analysis was conducted using ArchR72 (v.1.0.3), drawing on the sorted combined fragments matrix and utilizing EnsDb.Hsapiens.v.86 for annotations. Dimensionality reduction was achieved using iterative LSI, with Harmony rectifying any batch effects. Chromatin accessibility peaks were initially discerned using MACS2 (v.2.2.7.1)90, leading to the generation of a peak set consistent across the identified cell clusters. Unraveling the potential transcriptional regulatory networks, the cis-BP motif database was employed to annotate TF-binding motifs within these peaks. Distinct chromatin accessibility markers for cellular clusters were discerned, with biases like TSS enrichment and log(transformed fragment counts) in mind. A Wilcoxon’s test highlighted significant variations in peak accessibility across distinct cellular states. TF motifs and their activities were visualized and assessed at the cellular level using ChromVAR36 (v.1.16.0).

Assessment of genes and motif deviation scores that display spatial autocorrelation was conducted using Moran’s I with the package spdep91 (v.1.2-7). P values were adjusted for multiple comparisons using the Benjamini–Hochberg method.

To start, images for individual tiles were deconvoluted using Huygens (v.23.04) to enhance image signals and suppress background noises. Next, we proceeded with image registration for different sequencing rounds, employing the 3D Fourier transform implemented through the functions available in numpy.fft and Scipy92. Crosscorrelation between pairs of images at all translational offsets was computed and the position associated with the highest correlation coefficient was identified. This position was then utilized to translate image volumes to compensate for the offset. During this procedure, the first sequencing round was used as the reference, and subsequent rounds were registered to be aligned with it.

all six rounds. This vector was subsequently translated into a gene barcode using a barcode codebook, wherein each channel is linked with a two-base sequence encoding. Spots with undetectable signals in certain rounds were excluded from the analysis, as were spots with gene barcodes not present in the barcode codebook. ClusterMap (v.0.0.1)94 was employed for 2D cell segmentation of identified transcripts. Given the extensive scale of the data, ClusterMap was executed tile by tile. Image-free segmentation mode was utilized to comprehensively capture transcripts in regions of both nuclei and cytoplasm, with ‘xy_radius’ set as 55 and ‘cell_num_threshold’ set as 10−4. Transcripts that were not assigned to any cells were removed for downstream analysis.

Segmented cells from individual tiles were integrated for subsequent single-cell analysis using Scanpy (v.1.9.4)95. Cells with <80 transcripts were labeled as low quality and excluded. After cell filtering, gene expression raw counts were normalized to 10,000, log(transformed) and scaled. PCA was conducted to obtain the reduced dimensions. Harmony73 was employed to integrate cells from different tiles and the corrected PCs were used to construct a neighborhood graph and generate UMAPs. Leiden clustering was performed using the PCs with a resolution of 1.0. Single cells of STARmap were annotated by leveraging labels transferred from single-nucleus multiomics reference data using Seurat (v.4). First, integration anchors of reference data were identified using reciprocal PCA. Then, transfer anchors that connected reference single-nucleus data and STARmap data were identified using the ‘TransferData’ function with the top 30 dimensions. Last, labels from reference data were transferred to STARmap single-cell data using MapQuery. After label transfer, differential expression analysis was conducted utilizing Wilcoxon’s test of Scanpy function rank_genes_groups. P values were corrected using the Benjamini–Hochberg method.

We conducted the imputation of unmeasured genes by learning from snRNA-seq data following the integration of STARmap and single-cell multiomics reference data, employing a strategy similar to that described in ref. 96. Initially, we performed intermediate mapping to determine the optimal parameters for imputation. Specifically, we aligned the STARmap and single-cell multiomics reference data using Seurat (v.4), as detailed in the previous analysis. For each cell in the STARmap data, we computed its nearest neighbors within the single-cell multiomics dataset, representing a set of the most similar cells in the reference dataset. For each of the 1,001 genes that overlapped between STARmap and single-cell multiomics reference data, we imputed the expression level as the average expression within the neighborhoods of individual cells. We explored various neighborhood sizes, ranging from 5 to 400, by calculating Pearson’s correlation between the imputed expression levels and observed expression levels in STARmap data for individual genes. We then assessed the imputation performance using accumulated Pearson’s correlation values and selected the 100 nearest neighbors based on this score.

We employed CellChat64 to infer cell–cell interactions. For this purpose, we used the expression profiles of 33,357 imputed genes as input for CellChat and determined overexpressed genes and interactions using the human interaction database within CellChat. After this, we calculated communication probabilities while applying a distance constraint with a maximum interaction length of ligands set to 10 μm. We set the trim fraction to 0.1 and established a minimum threshold of 10 cells in each group for qualified communication. Only significant interactions were considered for subsequent analysis.

For ISH images, we initiated the process by detecting spots within the raw images and generating multiplexed images, encompassing multiple genes within a single tissue. The initial step involved aligning images from two rounds of acquisition using the 2D Fourier transform, a technique consistent with the image registration approach utilized for STARmap data, as previously described. Subsequently, we employed a method to identify regions with local maxima, which serve as signal spots for visualization. Specifically, we applied maximum and minimum filters from the scipy.ndimage.filters library to locate potential local maxima and minimal values within neighborhoods sized at 15 pixels. Regions where the difference between local maxima and minimal values exceeded a defined threshold of 100 were identified and marked as local maxima. Connected regions were delineated as a result of this identification process.

All tissue samples used for the present study were obtained with written informed consent from all participants and approval from the Medical University of Vienna Ethics Committee. All human tissue samples and experiments were approved by the Institutional Review Board at Broad Institute of MIT and Harvard and Massachusetts General Hospital. All experiments followed all relevant guidelines and regulations.

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at 10.1038/s41591-024-03073-9.

Supplementary InformationSupplementary Figs. 1–11, Tables 1–13 and References. Reporting Summary Supplementary Table 1Information about the human primary samples isolated from first trimester placentas and used in the present study along with accompanying assays: Multiome (combined snRNA-seq and snATAC-seq), Slide-tags, STARmap-ISS and STARmap-ISH. Each assay is accompanied by relevant QC metrics. Supplementary Table 2DEGs discovered by snRNA-seq. Supplementary Table 3Differentially accessible genes discovered by snATAC-seq. Supplementary Table 4Intercluster motif enrichment comparisons. Differentially accessible peaks were identified using a two-sided Wilcoxon’s test (FDR ≤ 0.1, log2(FC) ≥ 0.5). Intercluster motif enrichments were calculated via a hypergeometric test to generate P values. For instance, ‘vCTBvsEVT’ lists motifs enriched across peaks that are more accessible in vCTB clusters compared with EVT clusters. Supplementary Table 5ChromVAR motif enrichment analysis and accompanying positive TF regulators. Correlation between variables (motif enrichment, gene expression, gene accessibility) was calculated using Pearson’s correlation coefficient and statistical significance was assessed using P values adjusted by the Benjamini–Hochberg method (Padj < 0.01). Supplementary Table 6Cluster-specific differentially expressed marker genes across all clusters identified by Slide-tags, along with spatially autocorrelated genes and TF motifs with location-dependent expression and enrichment as identified by Slide-tags analysis using Moran’s I (includes correlations across all cells as well as within individual cell types). P values were adjusted by the Benjamini–Hochberg method (Padj < 0.05). Cluster numbering can be found in Supplementary Fig. 5. Supplementary Table 7Inferred peak–gene links, DORCs and per-cluster DORC scores. SnRNA-seq-identified tumor invasion and immunomodulation-related genes overlapped with DORCs, namely RUNX1, C12orf75, QSOX1, RASGRF2, PLXNB2, JAK1 and MYCN. Overlap with snATAC-seq-identified genes included ATP11A, DIO2, ANXA1, KLF6, ASCL2, NR2F6, SNAI2, TCF21, FLI1, PITX1, GRHL1 and NFIX.

and per-cluster DORC scores. SnRNA-seq-identified tumor invasion and immunomodulation-related genes overlapped with DORCs, namely RUNX1, C12orf75, QSOX1, RASGRF2, PLXNB2, JAK1 and MYCN. Overlap with snATAC-seq-identified genes included ATP11A, DIO2, ANXA1, KLF6, ASCL2, NR2F6, SNAI2, TCF21, FLI1, PITX1, GRHL1 and NFIX. Supplementary Table 8, 12Supplementary Table 8 Lineage drivers for chromatin potential/CellRank-derived terminal states (EVT3, STB, vCTB2) calculated by CellRank analyses. P values of the two-sided Fisher transformation tests were calculated and adjusted using the Storey–Tibshirani procedure for multiple hypothesis comparisons. Suppplementary Table 12 DEGs identified by STARmap-ISS. Statistics were derived using Wilcoxon’s test (two sided). P values were adjusted using the Benjamini–Hochberg method. Supplementary Table 9Description of UKBB pregnancy-related traits and accompanying heritability enrichments across cell types. Supplementary Table 10MAGMA GSEA of cell type-specific programs for each cell type as well as enrichment analysis of top GWAS hits for each pregnancy-related trait across cell types. Supplementary Table 11Genes, 1,001, for STARmap-ISS and 48 genes for STARmap-ISH and accompanying probe sequences used for ISS and ISH. Supplementary Table 13Ligand–receptor interactions across clusters and samples. Interaction (communication) probabilities and significance were computed by permutation test (P < 0.05). Supplementary Figs. 1–11, Tables 1–13 and References. Reporting Summary

Supplementary Table 8, 12Supplementary Table 8 Lineage drivers for chromatin potential/CellRank-derived terminal states (EVT3, STB, vCTB2) calculated by CellRank analyses. P values of the two-sided Fisher transformation tests were calculated and adjusted using the Storey–Tibshirani procedure for multiple hypothesis comparisons. Suppplementary Table 12 DEGs identified by STARmap-ISS. Statistics were derived using Wilcoxon’s test (two sided). P values were adjusted using the Benjamini–Hochberg method. Supplementary Table 9Description of UKBB pregnancy-related traits and accompanying heritability enrichments across cell types. Supplementary Table 10MAGMA GSEA of cell type-specific programs for each cell type as well as enrichment analysis of top GWAS hits for each pregnancy-related trait across cell types. Supplementary Table 11Genes, 1,001, for STARmap-ISS and 48 genes for STARmap-ISH and accompanying probe sequences used for ISS and ISH. Supplementary Table 13Ligand–receptor interactions across clusters and samples. Interaction (communication) probabilities and significance were computed by permutation test (P < 0.05). Supplementary Figs. 1–11, Tables 1–13 and References. Reporting Summary Information about the human primary samples isolated from first trimester placentas and used in the present study along with accompanying assays: Multiome (combined snRNA-seq and snATAC-seq), Slide-tags, STARmap-ISS and STARmap-ISH. Each assay is accompanied by relevant QC metrics. DEGs discovered by snRNA-seq. Differentially accessible genes discovered by snATAC-seq.

Information about the human primary samples isolated from first trimester placentas and used in the present study along with accompanying assays: Multiome (combined snRNA-seq and snATAC-seq), Slide-tags, STARmap-ISS and STARmap-ISH. Each assay is accompanied by relevant QC metrics. DEGs discovered by snRNA-seq. Differentially accessible genes discovered by snATAC-seq. Intercluster motif enrichment comparisons. Differentially accessible peaks were identified using a two-sided Wilcoxon’s test (FDR ≤ 0.1, log2(FC) ≥ 0.5). Intercluster motif enrichments were calculated via a hypergeometric test to generate P values. For instance, ‘vCTBvsEVT’ lists motifs enriched across peaks that are more accessible in vCTB clusters compared with EVT clusters. ChromVAR motif enrichment analysis and accompanying positive TF regulators. Correlation between variables (motif enrichment, gene expression, gene accessibility) was calculated using Pearson’s correlation coefficient and statistical significance was assessed using P values adjusted by the Benjamini–Hochberg method (Padj < 0.01). Cluster-specific differentially expressed marker genes across all clusters identified by Slide-tags, along with spatially autocorrelated genes and TF motifs with location-dependent expression and enrichment as identified by Slide-tags analysis using Moran’s I (includes correlations across all cells as well as within individual cell types). P values were adjusted by the Benjamini–Hochberg method (Padj < 0.05). Cluster numbering can be found in Supplementary Fig. 5.

F motifs with location-dependent expression and enrichment as identified by Slide-tags analysis using Moran’s I (includes correlations across all cells as well as within individual cell types). P values were adjusted by the Benjamini–Hochberg method (Padj < 0.05). Cluster numbering can be found in Supplementary Fig. 5. Inferred peak–gene links, DORCs and per-cluster DORC scores. SnRNA-seq-identified tumor invasion and immunomodulation-related genes overlapped with DORCs, namely RUNX1, C12orf75, QSOX1, RASGRF2, PLXNB2, JAK1 and MYCN. Overlap with snATAC-seq-identified genes included ATP11A, DIO2, ANXA1, KLF6, ASCL2, NR2F6, SNAI2, TCF21, FLI1, PITX1, GRHL1 and NFIX. Supplementary Table 8 Lineage drivers for chromatin potential/CellRank-derived terminal states (EVT3, STB, vCTB2) calculated by CellRank analyses. P values of the two-sided Fisher transformation tests were calculated and adjusted using the Storey–Tibshirani procedure for multiple hypothesis comparisons. Suppplementary Table 12 DEGs identified by STARmap-ISS. Statistics were derived using Wilcoxon’s test (two sided). P values were adjusted using the Benjamini–Hochberg method. Description of UKBB pregnancy-related traits and accompanying heritability enrichments across cell types. MAGMA GSEA of cell type-specific programs for each cell type as well as enrichment analysis of top GWAS hits for each pregnancy-related trait across cell types. Genes, 1,001, for STARmap-ISS and 48 genes for STARmap-ISH and accompanying probe sequences used for ISS and ISH.

Description of UKBB pregnancy-related traits and accompanying heritability enrichments across cell types. MAGMA GSEA of cell type-specific programs for each cell type as well as enrichment analysis of top GWAS hits for each pregnancy-related trait across cell types. Genes, 1,001, for STARmap-ISS and 48 genes for STARmap-ISH and accompanying probe sequences used for ISS and ISH. Ligand–receptor interactions across clusters and samples. Interaction (communication) probabilities and significance were computed by permutation test (P < 0.05).

Source Data Fig. 2aUncropped and unprocessed blots corresponding to Fig. 2a. Please note that. for the lower plot detecting p63, HLA-G and CGβ, an additional internal loading control has been added (TOPOIIbeta). The red stippled lines demarcate the protein bands shown in Fig. 2a. Source Data Extended Data Fig. 2aUncropped and unprocessed blots corresponding to Extended Data Fig. 2a. The red stippled lines demarcate the protein bands shown in Extended Data Fig. 2a. Uncropped and unprocessed blots corresponding to Fig. 2a. Please note that. for the lower plot detecting p63, HLA-G and CGβ, an additional internal loading control has been added (TOPOIIbeta). The red stippled lines demarcate the protein bands shown in Fig. 2a. Uncropped and unprocessed blots corresponding to Extended Data Fig. 2a. The red stippled lines demarcate the protein bands shown in Extended Data Fig. 2a.