A spatially resolved timeline of the human maternal–fetal interface - Nature.com

kembaliui.blogspot.com

Abstract

Beginning in the first trimester, fetally derived extravillous trophoblasts (EVTs) invade the uterus and remodel its spiral arteries, transforming them into large, dilated blood vessels. Several mechanisms have been proposed to explain how EVTs coordinate with the maternal decidua to promote a tissue microenvironment conducive to spiral artery remodelling (SAR)^1,2,3. However, it remains a matter of debate regarding which immune and stromal cells participate in these interactions and how this evolves with respect to gestational age. Here we used a multiomics approach, combining the strengths of spatial proteomics and transcriptomics, to construct a spatiotemporal atlas of the human maternal–fetal interface in the first half of pregnancy. We used multiplexed ion beam imaging by time-of-flight and a 37-plex antibody panel to analyse around 500,000 cells and 588 arteries within intact decidua from 66 individuals between 6 and 20 weeks of gestation, integrating this dataset with co-registered transcriptomics profiles. Gestational age substantially influenced the frequency of maternal immune and stromal cells, with tolerogenic subsets expressing CD206, CD163, TIM-3, galectin-9 and IDO-1 becoming increasingly enriched and colocalized at later time points. By contrast, SAR progression preferentially correlated with EVT invasion and was transcriptionally defined by 78 gene ontology pathways exhibiting distinct monotonic and biphasic trends. Last, we developed an integrated model of SAR whereby invasion is accompanied by the upregulation of pro-angiogenic, immunoregulatory EVT programmes that promote interactions with the vascular endothelium while avoiding the activation of maternal immune cells.

Main

Normal development during healthy pregnancy depends on a complex interplay between maternal cells and placental trophoblasts that ultimately transforms the womb into a specialized niche capable of meeting the metabolic demands of a growing semi-allogeneic fetus while maintaining maternal tolerance^1,2. After implantation, the decidua is invaded by EVTs. EVTs and maternal cells remodel uterine spiral arteries into highly dilated vessels with minimal smooth muscle where EVTs have partially replaced the maternal endothelium³. In healthy pregnancies, SAR results in low-resistance vessels that deliver blood to the intervillous space at low flow velocities that prevent damage to the placental architecture⁴. Conversely, impaired SAR, low numbers of tolerogenic maternal cells and abnormal decidual invasion of EVTs have each been implicated in placenta-related obstetric complications, including preeclampsia, intrauterine growth restriction and preterm birth⁵. Therefore, a detailed investigation of the cell population dynamics at the maternal–fetal interface is key to understanding the biology of normal pregnancy and obstetric complications.

Owing to the poor feasibility of controlled studies in pregnant humans, much of what is known about maternal–fetal tolerance and SAR is based on small-animal models⁶. Although some similarities exist, key facets of haemochorial placentation in humans are primate-specific⁷. For example, EVT giant cells in mice do not replace the vascular endothelium and are thought to play a minor part in SAR compared to maternal uterine natural killer (NK) cells⁸. The extensive degree of EVT invasion in humans is thought to be an evolutionary adaptation that permitted upright, bipedal locomotion while maintaining adequate blood flow in the third trimester when brain development accounts for 60% of fetal metabolic needs⁹.

The study of human decidual remodelling is further complicated by additional inherent challenges. First, cell composition and structure are temporally dynamic. Therefore, aggregating data across different gestational ages (GAs) or observing a single time point may be misleading. Second, these dynamics are spatially coordinated in the local tissue microenvironment. For example, periarterial decidual NK cells are thought to contribute to SAR by initiating smooth muscle breakdown and by secreting chemokines that attract invading EVTs, whereas phagocytic macrophages are thought to facilitate clearance of the resultant apoptotic debris¹⁰. Overall, the formation of the human maternal–fetal interface involves sophisticated spatiotemporal coordination such that tissue composition, structure and function are inextricably coupled.

With this in mind, we constructed a multimodal spatiotemporal atlas of the human maternal–fetal interface. We leveraged archival tissue banks to assemble a cohort of maternal decidua from 66 women who underwent elective terminations of otherwise healthy pregnancies at 6–20 weeks of gestation, constituting a large single-cell study of the maternal–fetal interface. We performed subcellular imaging with multiplexed ion beam imaging by time-of-flight (MIBI-TOF)¹¹ using a 37-plex antibody panel designed to identify the location, lineage and function of all major maternal and fetal cells.

We also profiled the transcriptome of arteries, decidua and EVTs. To understand how SAR relates to local decidual composition, we developed new algorithms to quantify vascular morphology that enabled us to assign a remodelling score to each individual artery. We discerned which changes in decidual composition, transcriptome and structure were preferentially driven by GA, SAR or both. Overall, the frequencies and spatial distribution of maternal immune cells exhibited a strong temporal dependence that enabled us to predict GA exclusively on the basis of these features.

By contrast, EVT invasion and perivascular localization were the dominant drivers of SAR in the tissue microenvironment, and these processes correlated with extensive shifts in arterial transcription. Given these findings, we used our atlas to characterize the temporal nature of intravascular EVT invasion in the decidua basalis. The accumulation of perivascular EVTs around arteries preceded smooth muscle loss and the appearance of intravascular EVTs, which is consistent with a model whereby perivascular EVTs intravasate into the artery lumen. Taken together, these investigations support a cooperative interplay in the first half of pregnancy in which temporally dependent changes in decidual function permit placental EVTs to extensively alter the maternal uterine vasculature.

Multiplexed imaging of the decidua

As part of the Human BioMolecular Atlas Program initiative, we created a spatiotemporal tissue atlas of the human maternal–fetal interface in the first 20 weeks of pregnancy (Fig. 1a). The goal of this study was to comprehensively define the structure and composition of decidua and to understand how it evolves in the first two trimesters with respect to two axes: GA and maternal SAR. We first assembled a large retrospective cohort of archival formalin-fixed, paraffin-embedded placenta and decidua tissue from 66 individuals who underwent elective termination of pregnancies with no known fetal abnormalities. Archival tissue blocks were manually screened by a perinatal pathologist. Tissue sections stained with haematoxylin and eosin (H&E) were analysed to determine which samples contained decidua, and mostly regions that contained anchoring villi were selected. Then, regions of decidua that contained spiral arteries were demarcated, cored and assembled into two tissue microarrays (TMAs) of 1 mm and 1.5 mm cores. The final dataset included samples for 6–20 weeks of gestation (13.72 ± 4.8 weeks, mean ± SD) from 66 women of varying parity (1.45 ± 1.72), age (28.17 ± 5.9 years), body–mass index (28.19 ± 7.3 kg m^–2) and ethnicity (Fig. 1b–e and Supplementary Table 1). Owing to inherent limitations in how the tissue was procured, precise anatomical locations could not be determined. However, 61 out of 66 tissue blocks contained placental villi, which suggested that the majority of this cohort was sampled from the decidua basalis (Supplementary Table 1 and Methods).

Previous studies of intact tissue that examined only one or a few cell populations at a time reported shifts in maternal immune cells towards tolerogenic states that are permissive to invasion by fetal EVTs¹². To gain a more complete picture of the complex cell–cell interactions that establish maternal tolerance in the first half of pregnancy, we combined the strengths of targeted subcellular imaging with antibodies and spatial transcriptomics on serial co-registered sections to construct a comprehensive composite model of SAR and decidual remodelling (Fig. 1f).

For MIBI-TOF, we designed and validated a 37-plex antibody panel to map the functional state and location of all major maternal and fetal cell populations (Fig. 1f and Methods). This panel included canonical lineage-defining markers and ten functional markers previously implicated in maternal immune tolerance^13,14,15 (Fig. 1f).

For spatial transcriptomics, we used the NanoString GeoMx Digital Spatial Profiler (DSP) for whole transcriptome analysis of arteries, EVTs and decidua. Immunofluorescence imaging of TMAs stained with antibodies for HLA-G, vimentin (VIM) and smooth muscle actin (SMA) were used to define regions of interest (ROIs) specific for each of these histological features (Methods). In total, we collected whole transcriptome data from 13 individual arteries, their adjacent decidua, 5 samples of interstitial and 3 samples of intravascular EVTs (19 cores from 17 individuals; Methods).

For cell segmentation, we used an optimized version of our previously validated deep-learning pipeline that was refined for decidua-specific cell types using 93,000 manual annotations (Methods). In total, we identified 495,349 segmented cells across 211 images and classified them into 25 cell populations (Fig. 2a,b, Methods and Extended Data Fig. 1). Functional marker expression in these populations was determined using per-marker thresholds (Methods). Noteworthy histological features—such as arteries, vessels, glands, the cell columns and decidual tissue boundaries—were manually annotated in collaboration with a perinatal pathologist.

**Fig. 2: Multiplexed imaging of human decidua reveals the immune-tolerance-conducive composition of the maternal–fetal interface.**

Non-immune maternal (structural) cells accounted for the majority (56.3%) of all segmented events in the decidua and were predominantly composed of decidual fibroblasts (60.5%) and myofibroblasts (24.8%), with smaller contributions from vascular endothelial cells (7.6%) and glandular epithelial cells (7.1%; Fig. 2b). Notably, we observed a new, rare subset of TIGIT⁺ glandular cells (0.34% of glandular cells; Supplementary Information). Consistent with previous work¹⁶ that quantified maternal populations in situ, maternal immune cells (31% of all cells) were dominated by macrophages (47.6% of immune cells) and NK cells (42.6% of immune cells), with minor contributions from T cells (8% of immune cells), dendritic cells (1.3% of immune cells) and mast cells (0.5% of immune cells).

Decidual macrophages ubiquitously expressed both CD163 and CD206, which is consistent with a M2-polarized, tolerogenic phenotype (Fig. 2g). In line with previous work¹⁷ that showed pregnancy-specific recruitment, 77% of these cells expressed DC-SIGN (Fig. 2a). We further classified DC-SIGN⁺ macrophages into three subsets on the basis of CD11c and HLA-DR expression: Mac2a, Mac2b and Mac2c. DC-SIGN^– macrophages (Mac1a and Mac1b) were subclustered on the basis of CD68 expression (Fig. 2a,b).

Four subsets of NK cells (NK1, NK2, NK3 and NK4) were classified on the basis of the combinatorial expression of CD57, CD11c and CD8. NK1 cells (CD57^–CD16^low) constituted the largest NK cell population present, making up 59.7% of NK cells (Fig. 2a,b). A new CD57⁺ population of decidual NK cells (NK2, 25.8% of NK cells) that had only been previously identified in peripheral blood during pregnancy¹⁸ was also identified. Moreover, most of this population expressed the tissue-residency marker CD49a (94.2%; Extended Data Fig. 2a,b). As described below, the frequency and spatial distribution of NK2 cells suggests that they play a distinct role in SAR (Extended Data Fig. 2c).

T cells consisted of CD8⁺ cells (53.2% of T cells), NK T (NKT) cells (28.8% of T cells), CD4⁺ cells (17.1% of T cells) and sparse numbers of regulatory T (T_reg) cells (Fig. 2a,b). We identified a PD-1⁺-activated population of T_reg cells with a TIM-3⁺LCK⁺ subset that accounted for >50% of this population (Fig. 2a,b). Notably, both T_reg cells and NKT cells were the most proliferative cell populations. Together with CD8⁺ NK cells, T_reg cells and NKT cells expressed granzyme B more frequently than any other cell type. TIGIT was most frequently expressed by T_reg cells—a rare subset that has been suggested to bind PVR (also known as CD155) on EVTs². This interaction has been observed in the tumour microenvironment¹⁹ and may serve a similar role in driving maternal–fetal tolerance.

Fetal cells (12.6% of all cells) were primarily composed of four subsets of EVTs that were delineated on the basis of the combinatorial expression of HLA-G, CK7, CD57 and CD56 (Fig. 2a). HLA-G⁺ EVTs were CK7⁺ (EVT1a), CK7^– (EVT1b) or CD56⁺ (EVT1c) (Fig. 2c–f). EVT2 lacked HLA-G and were CD57^–CK7^low and were located predominantly at the base of attaching cell columns. Taken together, these data provide spatial context to previous work that used dissociated samples^15,20; that is, an ensemble of functional states in fetal and maternal cells are collectively aligned to maintain a tolerogenic niche.

SAR correlates with local composition

Perfusion of the intervillous space by uterine spiral arteries is the sole source of oxygen and nutrients to the growing fetus after the establishment of arterial flow. During the first half of pregnancy, these vessels undergo an extensive remodelling process that culminates in dilated, non-contractile vessels depleted of smooth muscle where the maternal endothelium has been partially replaced by EVTs. Abnormal SAR is associated with obstetric complications such as intrauterine growth restriction and preeclampsia⁵. However, it is still not fully understood which cell populations directly participate in SAR, how this process is locally regulated and to what extent these changes are synchronized with GA.

We therefore used our spatiotemporal atlas to construct a SAR trajectory to understand how this relates to temporal changes in decidua composition and structure²¹. Using artery size, smooth muscle disruption, endothelial continuity and EVT infiltration, we manually assigned each artery to one of five sequential remodelling stages on the basis of previously published criteria (Methods and Fig. 3a). To ensure that scoring was not biased by demographics of the individuals or the composition of neighbouring arteries and stroma, scoring was performed by blinded experts on cropped images in which only the artery of interest was visible. Out of 588 arteries, 186 were unremodelled and assigned to stage 1 (Fig. 3b,c). Stage 2 arteries (300 arteries) were characterized by moderate smooth muscle disruption and endothelial swelling (Fig. 3d,e). Stage 3 arteries (43 arteries) exhibited more dilation, smooth muscle loss and early endothelial disruption (Fig. 3f,g). Progression to stage 4 (34 arteries) was marked by the presence of EVTs within the arterial lumen (Fig. 3h,i). Fully remodelled stage 5 arteries (25 arteries) were identified on the basis of their very large size, near-complete smooth muscle loss and EVT endothelization (Fig. 3j,k, Extended Data Fig. 3a and Supplementary Table 2).

**Fig. 3: SAR progression significantly influences maternal–fetal interface composition.**

Although SAR correlated to some extent with GA (Spearman’s ρ = 0.28, P = 1.5 × 10⁻¹²), in many cases these were discordant. For example, at least one late-stage artery (stage 4–5) was present in 40% of week 8 samples, whereas minimally remodelled arteries were present throughout (Fig. 3l). Moreover, SAR staging of arteries from the same individual often varied significantly between tissue cores (32% of individuals had arteries that differed by at least two stages), which suggested that aspects of SAR are locally regulated by the tissue microenvironment (Fig. 3l and Extended Data Fig. 3b).

We exploited this discordance between SAR and GA to identify changes in decidual composition that were predominantly driven by one or the other. We first developed a quantitative scheme to assign a continuous remodelling score. For each artery, we extracted 35 parameters that described the same aspects of arterial morphology used for manual scoring (Fig. 3m, Methods and Extended Data Fig. 3c,d). Together with manual staging, we used this profile to construct a highly resolved pseudotime trajectory of SAR using linear discriminant analysis (LDA). This trajectory assigned a continuous remodelling score to each artery (δ) (Fig. 3m and Methods). Thus, each artery in our dataset could be mapped along a temporal or remodelling trajectory using GA or δ, respectively (Extended Data Fig. 4a–c).

Applying linear regression to these values per image, we determined which aspects of decidual remodelling were preferentially correlated with GA or SAR (Methods and Extended Data Fig. 4d,e). The frequency of decidual EVTs was better correlated with SAR, whereas changes in the proportion of maternal immune cells were mostly driven by GA (Fig. 3n). A notable exception to the latter correlation was observed within the NK cell compartment, in which the ratio of NK2 cells (CD57⁺) to NK1 cells (CD57^–) decreased with SAR progression (Fig. 3n).

To further investigate this finding, we examined how the spatial distribution of NK cells near arteries changed as SAR progressed (Methods). Notably, NK2 cells were the only subset of maternal immune cells to preferentially localize around arteries (Supplementary Table 3). NK2 cell accumulation around arteries spiked specifically at stage 2 of SAR, when smooth muscle swelling and disruption are maximal (P =2 × 10⁻³; Extended Data Fig. 2c). Notably, CD57 expression in human NK cells is associated with a cytotoxic phenotype in tumours¹⁸, which suggests that this subset could serve a similar role in mediating early smooth muscle disruption during SAR.

To create a transcriptional trajectory that integrated with our spatial proteomics data, we used NanoString DSP on serial sections of the TMA imaged by MIBI-TOF. We collected whole transcriptome profiles of 13 individual arteries at various stages of remodelling and their adjacent decidua (a total of 26 ROIs; Methods). Matching these samples with their respective MIBI-TOF images enabled us to assign a remodelling score and GA to each transcriptome profile (Methods). We then used a methodology similar to that presented in Fig. 3n to categorize genes that displayed temporal expression trends as correlated with either SAR, GA or both (Methods).

For arteries, 2,935 out of 18,695 genes exhibited significant trends, with most genes preferentially correlating with SAR (Fig. 3o). Within this group of genes, we identified 78 temporally synchronized gene ontology pathways, including modules related to vessel remodelling and translation (Methods and Supplementary Table 4). These pathways exhibited both monotonic and biphasic trends (Fig. 3o), which showed that SAR is a composite of interrelated processes that occur continuously and episodically. We identified 185 genes that peaked at stage 2 of remodelling before subsequently declining. This expression pattern correlated with perivascular enrichment of NK2 cells as indicated by our MIBI data (Extended Data Fig. 2c). In addition, this group of genes was enriched for genes related to collagen fibril organization and responses to bone morphogenic protein (Fig. 3o and Extended Data Fig. 5). Consistent with cell growth and subsequent apoptosis of arterial smooth muscle, translation-related genes followed a biphasic trend, peaking at around stage 3 of remodelling (Fig. 3o and Extended Data Fig. 5). We also observed continual downregulation of genes involved in Notch signalling as SAR progressed (Fig. 3o). Taken together, these multimodal data provide a fully integrated atlas of decidual remodelling that describes tissue structure, single-cell function and changes in transcriptional programmes.

Immune composition correlates with GA

We next interrogated these data to identify GA-dependent, temporal changes in decidual composition. This analysis revealed a substantial shift from a lymphoid-dominant to myeloid-dominant landscape. Images at weeks 6–8 (Fig. 4c,e) showed NK cells and T cells exhibiting cytotoxic (Fig. 4d) and immunosuppressive (Fig. 4e) phenotypes and greatly outnumbering macrophages (Fig. 4b,c). By contrast, images from weeks 16 to 20 were dominated by interstitial EVTs (Fig. 4a,f,g) and tolerogenic macrophages (Fig. 4h). To further evaluate this relationship, we asked whether immune cell composition in the decidua alone could be used to predict GA. Using immune features that were preferentially associated with GA (Fig. 3n), we trained and validated a ridge regression model on a per-image basis using a random 70/30 test–train split (Extended Data Fig. 6a). Notably, the model predicted GA in the withheld test set within 19 days of the true value (R² = 0.7; Fig. 4i). On inspecting the model weights, the relative contribution of decidual immune cells was consistent with the observed shift in the proportion of myeloid and lymphoid cells. Relative frequencies of T cells and NK cells were negatively correlated with GA, whereas total macrophage frequency was positively correlated with GA (Fig. 4j). Notably, a regression model for predicting SAR (δ) based on the same immune cell frequencies performed poorly (R² = 0.05; Extended Data Fig. 6b), which reinforced our hypothesis that these changes are driven by GA and not SAR.

**Fig. 4: A lymphoid-to-myeloid shift in immune-compartment composition is tightly correlated with GA.**

Using computational approaches validated in previous work to identify significant spatial enrichment of two cell types^22,23,24 (Methods), we observed that the majority of significant pairwise enrichments involved EVT, NK cells and macrophages (Extended Data Fig. 6c and Supplementary Table 3). Again, by examining these relationships on a per-image basis, we were able to distinguish spatial relationships that evolved dynamically with respect to GA (Methods). Of these relationships, the pregnancy-specific Mac2a population was involved in the largest number of pairwise enrichments, becoming more enriched around several NK cell and EVT subsets, even though NK cells were in decline.

Upregulation of tolerogenic markers with GA

Having examined the influence of GA and SAR in driving changes in the frequency of cell populations in the decidua, we next used a similar approach to understand how these two time axes correlate with shifts in functional marker expression. These data revealed three overarching trends. First, both SAR and GA were associated with dynamic changes in IDO-1 expression. We identified a GA-driven decline in IDO-1⁺ glandular cells, in line with previous observations of IDO-1⁺ glandular cells in the first trimester but not at term²⁵. We also observed a SAR-driven decline in IDO-1⁺ dendritic cells and an increase in IDO-1⁺ vascular endothelium that was comparably correlated with both GA and SAR (Fig. 5b,d). Second, consistent with the cell frequency analysis (Fig. 3n) in which NK1 cells exhibited a preferential increase with SAR, NK1 cells also exhibited a similar increase in Ki67⁺ frequency, becoming more proliferative as SAR progressed (Fig. 5a). Third, functional shifts in innate immunity were preferentially correlated with GA. All five macrophage populations upregulated either TIM-3 and/or its cognate ligand galectin-9 (GAL-9) with GA (Fig. 5a,b). This trend was most prominent in the Mac2a and Mac2b populations, in which a tightly correlated upregulation of both TIM-3 and GAL-9 was observed (Fig. 5e,g,h and Supplementary Information).

**Fig. 5: Coordinated upregulation of tolerogenic functional markers with GA.**

Notably, GAL-9 upregulation was also detected in fibroblasts at 12–20 weeks GA (Fig. 5b,g–i). In previous work, interactions between maternal immune and stromal cell populations have been implicated in the promotion of fetal tolerance²⁶. With this in mind, we next sought to determine whether the GAL-9⁺ fibroblasts subset was biased to colocalize within specific spatial niches. To answer this question, we quantified their frequency within ten tissue microenvironments that were identified by clustering the cell-type compositions of the closest neighbours of each cell (Methods). GAL-9⁺ fibroblasts were strongly biased to colocalize with CD57⁺ NK cells (NK2, microenvironment G; Fig. 5f). Notably, this trend was accompanied by a GA-dependent increase in the expression of inducible nitric oxide synthase (iNOS) in NK2 cells (Fig. 5a). Both TIM-3 and GAL-9 have been implicated in the suppression of antitumour surveillance by impairing the activity of cytotoxic NK cells and T cells in various human cancers²⁷. Together with the transient perivascular enrichment of NK2 cells observed in early SAR, these findings suggest that expression of these proteins by macrophages and fibroblasts could have a concerted tolerizing role with fetal EVTs to attenuate immune cytotoxicity subsequent to NK-cell-dependent disruption of arterial smooth muscle.

Spatiotemporal EVT phenotype

Although it is accepted that EVT plugs formed during early gestation near the intervillous space are derived from the cell columns, it is unclear whether intravascular EVTs found deeper in the decidua basalis share this origin. These cells could also arrive at this location through intravasation—a migratory route by which interstitial EVTs within the decidua invade the vascular lumen by transiting across the arterial wall³ (Extended Data Fig. 7a).

To examine this possibility, we used our spatiotemporal atlas to quantify how the phenotype and spatial distribution of EVTs evolve with respect to SAR. First, we manually defined feature masks in our images to demarcate cell columns and three decidual compartments: interstitial, perivascular and intravascular (Fig. 6a–c and Extended Data Fig. 7b). We then quantified EVT frequency in each. Together with our SAR temporal trajectory, we first used these data to answer a question that has been qualitatively explored in previous work²⁸. That is, whether the initial accumulation of EVTs is in the perivascular compartment (adjacent to arteries) or within the intravascular compartment. Perivascular EVTs were consistently present earlier in SAR than intravascular EVTs (Fig. 6d). Furthermore, out of all arteries with intravascular EVTs, 75% also had perivascular EVTs, which is a higher percentage than would be expected if retrograde migration was the primary source of intravascular EVTs (Extended Data Fig. 7c). In arteries where both types of EVTs were present, the ratio of intravascular to perivascular EVTs followed a smooth trend with respect to SAR, such that intravascular EVTs increased at the expense of perivascular EVTs (R² = 0.5, P = 9 × 10⁻¹²; Extended Data Fig. 7d).

**Fig. 6: Spatiotemporal EVT distributions suggest that intravasation is the predominant route of EVT invasion in superficial decidua.**

Loss of smooth muscle and endothelium have defining roles in determining the extent of SAR. Using morphometrics to quantify these cell types, we examined how the integrity of these concentric layers relates to EVT enrichment (Methods). Perivascular EVTs were consistently present at an earlier stage, whereas intravascular EVTs appeared only after 80% of smooth muscle was lost (Extended Data Fig. 7e). Notably, intravascular EVTs increased proportionally with endothelial loss, whereas perivascular EVTs were present irrespective of this process (Extended Data Fig. 7f). Together, these findings align with a process whereby near-complete loss of smooth muscle permits perivascular EVTs to invade the artery lumen through discontinuities in the vascular endothelium²⁹.

To further evaluate this model, we next examined whether the EVT phenotype shifts progressively in a manner consistent with a migratory route that passes through the decidua into the artery lumen. The proportion of EVT subsets in cell column, interstitial, perivascular and intravascular compartments shifted systematically in a manner consistent with this proposed path of migration (Extended Data Fig. 8a,b). Examining functional marker expression within each compartment again revealed a progressive shift in EVT phenotype. Cell columns were distinctly enriched for proliferative (Ki67⁺) and CD57⁺ EVTs. A progressive decrease in these markers across interstitial, perivascular and intravascular EVTs was accompanied by a concomitant increase in PD-L1 that peaked in the intravascular compartment (Fig. 6e, Extended Data Fig. 8c and Methods).

These analyses collectively align with a spatial trajectory whereby decidual invasion of cell column EVTs is synchronized with a downregulation of CD57 and Ki67 and an upregulation of HLA-G. Perivascular accumulation of EVTs occurs early in SAR before the appearance of intravascular EVTs and any loss in endothelium. In this model, as the endothelial barrier is lost, perivascular EVTs invading the artery lumen upregulate CD56 and PD-L1 (Extended Data Fig. 8d–h and Methods). Notably, owing to its observational nature, this study cannot rule out alternative mechanisms by which detaching EVTs arising from EVT plugs at the intervillous space undergo retrograde migration.

Irrespective of the migration route, the distinct changes in phenotypic markers measured by MIBI-TOF suggest that arterial invasion is accompanied by a shift in EVT transcriptional programmes. With this in mind, we used NanoString DSP to measure the transcriptomes of interstitial and intravascular EVTs. We found 274 differentially expressed genes (DEGs) (Extended Data Figs. 9 and 10a and Methods). In addition to confirming previous work noting an upregulation of NCAM1, JAG1 and LAIR2 specifically in intravascular EVTs^30,31,32, we identified transcriptional changes in genes important for extracellular matrix remodelling and angiogenesis (Fig. 6f). For example, MMP12, MMP15 and ADAMTS4 were specifically upregulated in intravascular EVTs (mean log(fold change) of 11.67 for MMP12, 9.74 for MMP15, and 7.88 for ADAMTS4), which suggested that these proteins play a significant role in late-stage SAR. In addition, arterial invasion was accompanied by a shift from VEGFR1 (also known as FLT1) to VEGFR3 (also known as FLT4) (Fig. 6f).

Intravascular EVTs also upregulated DKK3, C5orf30 (also known as MACIR) and CD24 (Fig. 6f), which have each been shown in previous work to have roles in fetal viability, tumour invasion or immune tolerance^33,34,35. Similarly, we observed an accompanying downregulation of genes associated with invasion in intravascular EVTs, such as MGAT5. With respect to immune modulation, C5orf30 is a potent immunometabolic regulator that has been shown to inhibit macrophage-mediated tissue damage in rheumatoid arthritis³³. Similarly, CD24 binding to Siglec-10 was recently found in many cancers to promote immune evasion by serving as an anti-phagocytic, ‘don’t eat me’ signal³⁶. Taken together, our multimodal approach paints a picture of a highly regulated and controlled process. We observe a transcriptional shift away from a more invasive phenotype (SERPINE1 and CORO6) in interstitial EVTs towards genes implicated in vascular remodelling in intravascular EVTs. Notably, this was accompanied by an increase in immunoregulatory modules that allow EVTs to be in continuous contact with maternal blood while avoiding immune activation (Fig. 6e,f).

To understand how these changes promote SAR, we investigated potential cell–cell interactions between intravascular EVTs and arterial cells using NicheNet³⁷ (Methods). We identified ten protein–protein interactions between EVTs and arteries that were predicted to affect 121 downstream targets (Fig. 6g,h, Extended Data Fig. 10b and Supplementary Table 5). For example, interactions between EVT JAG1 and arterial Notch were predicted to drive downstream changes in arterial MEOX1 and MT2A, which have been implicated in endothelial dysfunction and apoptosis^38,39. Similarly, CGA–FSHR and LEP–LEPR interactions correlated with changes in arterial hormone receptors (LHCGR) and several cell adhesion targets, respectively (Fig. 6h and Extended Data Fig. 10b). Notably, among the most prominent downstream targets were the olfactory receptors OR51E2 and the human-specific OR7A5, the expression of which outside the olfactory bulb has been thought to regulate blood pressure and angiogenesis^40,41.

CD24–SELP was the second most significant interaction and had several targets related to blood vessel function and formation (Fig. 6h and Extended Data Fig. 10b). Notably, reduced CD24 mRNA levels in bulk placental samples has been associated with higher preterm preeclampsia risk³⁵; however, EVT-specific CD24 expression has not been previously reported. Given that abnormal decidual and SAR are thought to play a major part in preeclampsia⁴², we sought to determine whether other genes involved in EVT invasion and vascular remodelling had previously been implicated. To do this, we first compared our list of EVT DEGs with genes found previously to be differentially expressed in decidua samples from women diagnosed with preeclampsia⁴³. We found that 31% of EVT DEGs (12 genes) are differentially expressed in preeclamptic decidua (Fig. 6f and Supplementary Table 6). Notably, FN1 and FLT1, which have been proposed as biomarkers for early prediction of preeclampsia^44,45, were markedly downregulated in intravascular EVTs.

Half of the NichNet interactions and 19 downstream targets overlapped with this list of preeclampsia DEGs⁴³ (Fig. 6h and Extended Data Fig. 10b). These included WNT10B, a newly identified accelerator of EVT invasion⁴⁶, and OR51E2, a target of CD24–SELP signalling that also exhibited the highest regulatory potential. With respect to the latter, SELP is notable for being differentially expressed in peripheral blood-cell-free RNA from patients with preeclampsia⁴⁷. Taken together, our transcriptomics approach validated and complemented the stepwise changes in EVT phenotype seen in our spatial atlas while revealing pathways that are perturbed in pregnancy-related disorders.

Discussion

Decidualization is a fascinating process with no other normative precedent in human biology. In this process, the structure and function of the maternal endometrium transforms to promote the regulated invasion of genetically dissimilar fetal cells. The decidua plays a dual role by permitting EVT invasion in the first trimester and later limiting it by inducing EVT apoptosis⁴⁸. EVT invasion can also be limited by morphological changes such as EVT fusion, which leads to polyploidization that limits invasion owing to nuclear size⁴⁹. Given the lack of tractable and relevant animal models and the inability to study decidualization prospectively, our understanding of it is immature relative to other areas of human physiology. Therefore, our study aimed to understand how global, temporally dependent changes in decidual composition are coupled to local regulation of vascular remodelling in pregnancy. Initial invasion of placental EVTs is prompted by a shift towards a permissive milieu, whereas progression of SAR depends on the subsequent migration and perivascular accumulation of EVTs, where they are thought to participate in cooperative cell–cell interactions with maternal fibroblasts, NK cells and macrophages². Thus, the formation of the maternal–fetal interface is mediated by global, temporally dependent cues that serve as a gating function for remodelling processes that are regulated in the local tissue microenvironment.

With this paradigm in mind, we set out to delineate which aspects of the first half of pregnancy are driven globally by GA and how this relates to SAR. In the study of placentation and SAR, an ideal sampling strategy might use elective caesarean hysterectomies from normal pregnancies performed across GA in an ethnically diverse patient population. As ethical considerations prohibit this approach, previous work has used a range of sample types that each have their own strengths and weaknesses. Here we utilized archival tissue from elective terminations with no known pregnancy complications. This enabled us to examine these questions in a large, ethnically diverse cohort that is well-distributed with respect to GA. As tissue procured during terminations is fragmented, anatomical registration for determining whether these tissue blocks were sampled from central or peripheral regions of the decidua basalis was not feasible.

Using LDA, image morphometrics and expert annotations, we assigned quantitative remodelling scores to every spiral artery in these images. These targeted multiplexed imaging data were complemented by spatially co-registered tissue transcriptomics. This multimodal dataset enabled us to reveal how cell frequency and function, tissue organization and transcriptional programmes in maternal decidua, arteries and EVTs change with SAR and GA.

Our analyses of these changes determined that GA is the predominant driver of maternal immune cell recruitment. Progressive decreases in the numbers of NK cells and T cells drive a transition at 12–14 weeks GA from a lymphoid-dominant to myeloid-dominant decidua enriched for iNOS⁺ NK cells, IDO-1⁺ vascular endothelium and DC-SIGN⁺ macrophages that express both TIM-3 and GAL-9. Notably, this relationship between immune composition and GA was strong enough to allow us to predict GA within 19 days exclusively on the basis of immune population frequencies.

By contrast, all EVT subsets and only two maternal cell populations (NK1 and NK2) preferentially correlated with SAR. Higher remodelling scores were correlated with more EVTs, more NK1 cells and fewer NK2 cells. A sharp accumulation in NK2 cells around arteries was observed early in the remodelling process around the time smooth muscle disruption had been initiated. NK1 and NK2 cells primarily differed in that the latter express CD57, a marker associated with a cytotoxic phenotype. Higher proportions of presumptively more reactive NK2 cells early in SAR around arteries aligns well with previous results⁵⁰ that have suggested that decidual NKs initiate early disruption of arterial smooth muscle through the secretion of GrB, MMP2 and MMP9. Likewise, the proportional gains seen here as SAR progresses of less reactive NK1 cells and invasive EVTs are consistent with the tolerizing effects of HLA-G, which has previously been shown^51,52 to decrease NK cell cytotoxicity and induce the production of interleukin-6 and interleukin-8 through the binding of HLA-G to KIR2DL4, LILRB1 and LILRB2. Taken together, these data suggest that maternal and fetal cells have cooperative, interdependent roles with SAR, transitioning through NK-dependent and EVT-dependent phases.

We also examined a lingering question³ in the field: the path of migration taken by intravascular EVTs in the decidua basalis. In line with early work based on studies of 8–18 week hysterectomy specimens processed in toto³, we found that perivascular EVTs accumulated before intravascular EVTs. By comparing the cellular composition within cytotrophoblast cell columns of anchoring villi, decidua and arteries, we observed a sequential and coordinated shift in EVT frequency and phenotype, which provided support for a model in which EVTs enter spiral arteries from nearby decidua by crossing the arterial wall.

Notably, previous studies of samples from caesarean hysterectomies identified morphological evidence of arterial extravasation³. Given the observational nature of this study and the limited preservation of tissue structure at the intravillous space, an extravasation model in which EVTs migrate in a retrograde manner after entering spiral arteries directly at the basal plate cannot be excluded. We also note the possibility that following intravasating arteries in the decidua, EVTs could migrate upstream to reach the upper third of the myometrium. This idea would be consistent with previous studies in which perivascular trophoblasts become increasingly scarce as a function of myometrial depth⁵³. Therefore, it is conceivable that both processes may be at play in different regions of the decidua as pregnancy progresses.

Previous single-cell and bulk sequencing studies of decidua have characterized the transcriptome of decidual cells; however, they were performed using dissociated tissue, agnostic to spatial context and the local extent of SAR^2,12. Correlating spatial morphology and tissue composition with targeted tissue transcriptomics enabled us to observe how the transcriptome evolves with respect to SAR. In arteries, our analysis revealed a downregulation of Notch signalling, tissue organization and cohesion with SAR, which was accompanied by a burst of translation-related activity around stage 2 of remodelling. By comparing interstitial EVTs with intravascular EVTs, our analyses revealed genes upregulated in the interstitial populations that shed light on how EVTs facilitate immune tolerance. Almost one-third of DEGs between interstitial and intravascular EVTs overlapped with DEGs in preeclamptic decidua samples. Given the significant contribution that abnormal vascular remodelling and EVT invasion are thought to play in preeclampsia, this work serves as a valuable resource for contextualizing preeclampsia-related changes in future studies.

Notably, many of these pathways are also associated with cancer progression. Formation of the maternal–fetal interface is an organized and controlled invasive process that is sometimes viewed as a template for understanding invasive and immunosuppressive properties of tumours⁵⁴. Both processes involve a genetically dissimilar invasive cell type (haploidentical EVTs versus clonal, mutated cancer cells), extracellular matrix remodelling and recruitment of a wide variety of tolerogenic immune cells, including M2-polarized macrophages and proliferating T_reg cells. The intersection of anchoring placental villi and maternal decidua morphologically resembles the invasive margin of carcinomas and contains trophoblast cells that express high levels of immunomodulatory proteins and growth factors implicated in tumour severity, including PD-L1, IDO-1, TIM3, HER2 and EGFR^27,55,56. In addition to these phenotypic and structural similarities, recent work that revealed mosaicism and clonal mutations in normal-term placentas demonstrated that this phenotypic overlap is even manifest at a genomic level⁵⁷.

Overall, we anticipate that this spatiotemporal atlas of the early human maternal–fetal interface will provide a normative framework for elucidating aetiological perturbations in maternal–fetal tolerance and SAR in pregnancy complications. Likewise, this work may also serve as a template for understanding how immune tolerance, tissue remodelling and angiogenesis are aberrantly recruited and synergized during tumour progression. With this in mind, we plan in future studies to extend this comparative approach to archival tissue in the context of obstetric complications to further elucidate cellular interactions involved in the regulation of SAR and EVT invasion.

Methods

Retrospective cohort design

The study cohort comprised decidua tissue from archival formalin-fixed, paraffin embedded (FFPE) blocks, sampled after elective pregnancy terminations from an outpatient clinic located within a large public hospital affiliated with an academic medical centre. Patients at this clinic reflect a diverse population. Although the patient population is predominantly low-income, women of all economic backgrounds are cared for at the clinic.

In the clinic, an ultrasound examination is performed to estimate GA, and a medical history is taken and logged as an electronic medical record (electronic clinical works) or handwritten forms. A board-certified gynaecologist reviewed medical records and specifically extracted the following details: age, ethnicity, body–mass index, gravidity, parity, previous terminations, smoking status, medications, HIV status, history of preeclampsia, chronic hypertension, diabetes mellitus, renal disease, autoimmune disease, multifetal pregnancy, and congenital anomalies (Supplementary Table 1). For procedures occurring at <14 weeks GA, suction aspiration is routinely used. For procedures at >14 weeks GA, a combination of suction aspiration and grasping forceps is used. After the procedure, tissue samples are routinely sent to pathology.

TMA construction

Whole tissue sections from individuals who underwent elective termination at 6–20 weeks of gestation were first reviewed by H&E staining to identify samples containing decidual tissue and spiral arteries. These regions were manually demarcated and assessed for suitability. Blocks containing decidua with vessels were selected, cored with a bore needle and assembled into the TMA used in this study. Archival tissue blocks from 74 individuals were initially selected by a board-certified perinatal pathologist (G.R.) to be included in the TMAs. The first TMA consisted of 205 cores (including 3 tonsil cores, 1 endometrium core and 1 myometrium core) of 1 mm in diameter and the second contained 86 cores of 1.5 mm in diameter). Unfortunately, cores from eight individuals did not end up containing decidua, and there was not sufficient tissue in the block for additional re-coring. We therefore had to exclude these samples from the analysis. The final cohort included 66 individuals, an exhaustive list of which is provided in Supplementary Table 1. Images from samples from six individuals did not have arteries and therefore were not included in analyses related to spiral arteries. Information on the histological characteristics of the blocks retrieved, including the presence of cell column anchoring villi, is in Supplementary Table 1. High-resolution scans of each core were uploaded to the Stanford Tissue Microarray Database (http://tma.im/cgi-bin/home.pl), a collaborative internal platform for designing, viewing, scoring and analysing TMAs. Sequential recuts of the main experiment were stained with H&E to aid in choosing the imaging ROIs and analysing data.

Antibody preparation

Antibody staining was validated as previously described^11,58. In brief, each reagent was first tested using single-plex chromogenic immunohistochemistry (IHC) using multiple positive and negative FFPE tissue controls before metal conjugation. Antibodies were then conjugated to isotopic metal reporters as previously described^{11,22,23,24,58} with the exception of biotin-conjugated anti-PD-L1, for which a metal-conjugated secondary antibody was used. The performance of metal-conjugated antibody reagents were then tested within the complete MIBI-TOF staining panel under conditions identical to those in the main study and compared with representative single-plex chromogenic IHC to confirm equivalent performance. Representative stains and information for each marker is provided in the Supplementary Information and in Supplementary Table 7, respectively. After conjugation, antibodies were diluted in Candor PBS Antibody Stabilization solution (Candor Bioscience). Antibodies were either stored at 4 °C or lyophilized in 100 mM d-(+)-trehalose dehydrate (Sigma Aldrich) with ultrapure distilled H₂O for storage at −20 °C. Before staining, lyophilized antibodies were reconstituted in a buffer of Tris (Thermo Fisher Scientific), sodium azide (Sigma Aldrich), ultrapure water (Thermo Fisher Scientific) and antibody stabilizer (Candor Bioscience) to a concentration of 0.05 mg ml^–1. Information on the antibodies, metal reporters and staining concentrations is in Supplementary Table 7.

Tissue staining

Tissues were sectioned (4 μm in thickness) from tissue blocks on gold and tantalum-sputtered microscope slides. Slides were baked at 70 °C for 20 min followed by deparaffinization and rehydration with washes in xylene (3 times), 100% ethanol (2 times), 95% ethanol (2 times), 80% ethanol (once), 70% ethanol (once) and ddH₂O with a Leica ST4020 Linear Stainer (Leica Biosystems). Tissues next underwent antigen retrieval, which was carried out by submerging sides in 3-in-1 Target Retrieval solution (pH 9, Dako Agilent) and incubating them at 97 °C for 40 min in a Lab Vision PT Module (Thermo Fisher Scientific). After cooling to room temperature, slides were washed in 1× PBS IHC washer buffer with Tween 20 (Cell Marque) with 0.1% (w/v) BSA (Thermo Fisher). Next, all tissue samples underwent two rounds of blocking, the first to block endogenous biotin and avidin with an Avidin/Biotin Blocking kit (BioLegend). Tissue samples were then washed with wash buffer and blocked for 1 h at room temperature with 1× TBS IHC wash buffer with Tween 20 and 3% (v/v) normal donkey serum (Sigma-Aldrich), 0.1% (v/v) cold fish skin gelatin (Sigma Aldrich), 0.1% (v/v) Triton X-100 and 0.05% (v/v) sodium azide. The first antibody cocktail was prepared in 1× TBS IHC wash buffer with Tween 20 and 3% (v/v) normal donkey serum (Sigma-Aldrich) and filtered through a 0.1 μm centrifugal filter (Millipore) before incubation with tissue overnight at 4 °C in a humidity chamber. After overnight incubation, slides were washed for 2 min in wash buffer. The next day, the antibody cocktail was prepared as described (Supplementary Table 7) and incubated with the tissues for 1 h at 4 °C in a humidity chamber. After staining, slides were washed twice for 5 min in wash buffer and fixed in a solution of 2% glutaraldehyde (Electron Microscopy Sciences) solution in low-barium PBS for 5 min. Slides were washed in low-barium PBS for 20 s then, using a linear stainer, through 0.1 M Tris at pH 8.5 (3 times), ddH₂O (2 times) and then dehydrated by washing in 70% ethanol (once), 80% ethanol (once), 95% ethanol (2 times) and 100% ethanol (2 times). Slides were dried under vacuum before imaging.

MIBI-TOF imaging

Imaging was performed using a custom MIBI-TOF instrument with a Xe⁺ primary ion source, as previously described^22,58. In total, 222 808 × 808 µm FOVs were acquired at approximately 600 nm resolution using an ion dose of 7 nA × h mm^–2. After excluding 11 FOVs that contained necrotic or non-decidual tissue, or consisted of duplicate tissue regions, the final dataset consisted of 211 FOVs from 66 individuals.

Low-level image processing

Multiplexed image sets were extracted, slide background-subtracted, denoised and aggregate filtered as previously described^{22,23,24,58,59}. For several markers, a background channel consisting of signal from the mass 128 channel was used. All parameters used as inputs for low-level processing are listed in Supplementary Table 7.

Feature annotation

Large tissue features were manually annotated in collaboration with a perinatal pathologist. Pseudo-coloured MIBI images stained with H3 to identify cell nuclei, VIM for decidual stromal cells, SMA and CD31 for vessels, cytokeratin 7 (CK7) for glands and the fetal cell columns, and HLA-G for EVTs were used to guide annotation. Serial H&E sections, and a H&E recut of the entire block, if necessary, were additionally used to supplement annotation. Labelling was performed in ImageJ and the annotated features were exported as binary TIF masks.

Single-cell segmentation

The Mesmer segmentation algorithm⁶⁰ was adapted specifically to segment the cells in our dataset. First, training data were generated using a subset of 15 images out of 211 in our cohort, in addition to 10 decidua MIBI-TOF images from titration data. In total, 1,024 × 1,024 pixel crops were selected to encompass the range of different cell morphologies present. The markers H3, VIM, HLA-G, CD3, CD14 and CD56 were used to capture the major cell lineages present. Subsequently, a team of annotators parsed these images to identify the location of each unique cell using DeepCell Label, custom annotation software specifically developed for this task⁶⁰ (https://github.com/vanvalenlab/deepcell-label). The manually annotated images were used to generate partially overlapping crops of 256 × 256 pixels from each image. In total, training data included 1,600 distinct crops with 93,000 cells. This dataset was used to retrain the Mesmer segmentation model, modifying the architecture to accept six distinct channels of input. The output from the network was then post-processed using the default model settings (Supplementary Information).

Segmentation post-processing

Examination of the images revealed that glandular cells and chorionic villus trophoblasts did not express any markers included in the training data; namely these cells were predominantly CK7⁺. This resulted in effectively nuclear-only segmentation being predicted by the convolutional neural network within these features. To account for this, segmented cells that overlapped with the gland mask were expanded radially by 5 pixels, and those in the cell column mask by 2 pixels. The number of pixels used for expansion was optimized to approximate the observed cell size, which was based on a systematic inspection of three images per GA. Objects <100 pixels in area were deemed noncellular and excluded from subsequent analyses. The final number of segmented events per FOV is provided in Supplementary Table 8.

Single-cell phenotyping and composition

Single-cell expression data were extracted for all cell objects and area-normalized. Single-cell data were linearly scaled with a scaling factor of 100 and ArcSinh-transformed with a co-factor of 5. All mass channels were normalized to the 99th percentile. To assign decidual cell populations (≥70% cell area in decidua) to a lineage, the clustering algorithm FlowSOM (Bioconductor FlowSOM package in R)²⁹ was used, which separated cells into 100 clusters based on the expression of 19 canonical lineage-defining markers (Supplementary Information). Clusters were further classified into 21 cell populations, with proper lineage assignments ensured by manual examination of overlayed FlowSOM cluster identity with lineage-specific markers. Clusters containing non-biologically meaningful or distinct signals were assigned the label ‘other’. T_reg cells were identified by thresholding T cells (FlowSOM clusters 43, 53 and 63) with the CD3 signal ≥ the mean CD3 expression of CD4⁺ T cells and >0.5 the normalized expression of FOXP3. Mast cells were identified as cells for which normalized expression of tryptase was >0.9. Mac2b (CD11c⁺) cells were identified as macrophages with >0.5 normalized expression of CD11c. Placental macrophages (Hofbauer cells) were defined as CD14⁺ >0.5 cells located within the cell column. Cells from FlowSOM clusters 4, 5 and 15 ubiquitously and predominantly expressed CK7 and were reassigned to the EVT2 subset if located within the cell column feature mask or as glandular cells otherwise (Supplementary Information). These thresholds were selected based on the distribution of lineage marker expression (Supplementary Information) and on systematic examination of the images by eye as expression patterns varied significantly between markers. For a comprehensive list of all single cells, their morphological features, markers expression, lineage classification, among others, see the Data availability section.

Definition of thresholds for functional marker positivity

Cells were considered positive for a functional marker if their scaled expression level was greater than or equal to a set threshold, as previously described²². Thresholds for individual functional markers were determined on the basis of examining the images by eye, as expression patterns varied significantly between markers (Supplementary Table 9 and Supplementary Information). To set the per-marker thresholds, five images for each functional marker were reviewed, and increasing threshold values were examined using custom software. Subsequently, cells defined as negative for a marker based on the determined threshold value were re-examined to ensure that the thresholds were representative. For Ki67 positivity, only cells that had a nucleus in the image were considered. Ki67 values were not normalized to the cell size because the Ki67 signal is exclusive to nuclei.

Two-colour IHC

Before staining, FFPE sections were incubated at 70 °C for 1 h. After deparaffination and antigen retrieval (Dako, S2367) was performed, endogenous horseradish peroxidase and alkaline phosphatase were blocked using BLOXALL (Vector Laboratories, SP-6000-100) for 30 min, followed by blocking buffer solution (95% 1× TBS IHC wash buffer with Tween 20, 1% Triton 10%, 1% gelatin 10%, 2% horse serum and 1% sodium azide 20 mg ml^–1) for 1 h at room temperature. Double staining was performed using CD57 (mouse IgG) paired with CD49a (rabbit IgG). Sections were incubated at 4 °C overnight with the antibodies CD57 (clone NK/804, Abcam, ab269771; titre, 0.5 µg ml^–1) and CD49a (clone E9K2J, CST, 15574T; titre, 1:1,500). The following day, secondary antibody (ImmPRESS Duet reagent; HRP anti-rabbit IgG and AP anti-mouse IgG; Vector Laboratories) was applied for 10 min at room temperature. Antibodies were revealed with Vector Blue AP substrate (Vector Laboratories, SK-5300) for 10 min in the dark followed by DAB HRP substrate (Vector Laboratories, SK-4105) for 40 s. For subsequent analyses and colour deconvolution, single-plex staining for CD57 and CD49a were performed on one slide each. For details on the method, buffers and solutions, refer to ref. ⁶¹.

The IHC slides were scanned using a NanoZoomer Digital Pathology Scanner 2.0RS (Hamamatsu) and analysed using QuPath (v.0.4.0). To score CD57⁺ NK cells (NK2) for expression of the tissue-residency marker CD49a, the two colours in the IHC slides were deconvolved with QuPath using single-plex staining as colour references. CD57⁺ NK cells were then manually annotated, in decidual regions only, by a board-certified pathologist. These cells were than manually scored for CD49a expression and counted.

Blinded manual artery staging

Arteries were categorized into five remodelling stages based on criteria adapted from a previously proposed four-stage model²¹. These criteria were used to describe spiral arteries observed in H&E and single-channel IHC images and were adapted to suit multiplexed MIBI data (Fig. 3a, details in Extended Data Fig. 3a). In total, 600 arteries were categorized according to these criteria by a single reviewer using only crops of MIBI pseudocolour overlays (SMA, VIM, CD31, H3 and HLA-G), including only the artery (as defined by a feature mask) and any EVTs in the lumen. The reviewer was blinded to the rest of the image, serial H&E sections, GA and any clinical data. Twelve partially captured arteries were excluded from the final dataset of 588 arteries.

Automated digitization of artery morphological features

The same format of cropped artery MIBI images that were manually scored by the reviewer were used to calculate a set of geometric parameters for several selected features. These features described the organization and structure of the vessel wall, the continuity of the endothelium and its thickness, and the presence and structure of intravascular EVTs. To capture these features, a structure of concentric circles we termed the ‘onion’ structure was defined, with the outer circle of this structure enclosing the artery and the inner circles dividing it into layers. This structure is described below using the two-dimensional cylindrical coordinate system, with the radial axis r, azimuthal (angular) axis ø, and origin of the axis at point (x,y). Point (x,y) is the user-defined artery centre. For an artery in the binary mask M, the following algorithm was used to create the onion structure (Extended Data Fig. 3c).First, define a circle enclosing the artery, centred at point (x,y) with radius a as follows: (x,y) was taken as the user-defined artery centre point; a, the radius is defined as the maximum distance between (x,y) and the edge of M, rounded up to the nearest integer multiple of n, such that a = I × n for an integer I. n is a user-defined thickness parameter for the onion layers

Second, define the inner circles comprising the onion layers by dividing the radius a of the outer circle into I equal sections of length n, creating layers along the radial r axis. The radii of the inner circles are then defined as 0,1 × n,2 × n,…(I – 1) × n.

Third, divide the onion into k equal sectors along the ø axis. k is a user-defined integer.

Fourth, subdivide each sector into segments. The sectors are internally divided by the circles, creating parts with four corners and four sides, with the two sides being straight (sector dividers), and the two sides being arcs (parts of circle circumferences). The arcs are replaced with secants (straight line connecting the ends of the arc), turning the segment into a trapezoid. The parameters n = 10 pixels and k = 100 were used to allow for segments large enough to contain a sufficient number of pixels to average the expression over.

Geometrical and protein morphology features were then extracted for each artery onion. For geometrical features, the following parameters were defined: (1) radius, the maximum distance between any pixel within the mask and the closest pixel on the edge of the mask; (2) perimeter, the Euclidean distance between all adjacent pixels on the edge of the artery mask; and (3) area, the total number of pixels within the artery mask.

For the protein morphology features, for the markers CD31, CK7, H3, HLA-G, SMA and VIM, the following parameters were defined. (1) Average signal: the weighted average over segments of marker expression, in which the weight of a segment corresponds to the number of pixels it contains. The weighted average was used to avoid smaller inner segments having a disproportionate effect on the average.(2) Thickness: for each sector, we calculated the distance d between the inner-most segment positive for the marker and the outer-most positive segment. Positivity was measured by comparing the mean signal over pixels the segment to a user-defined threshold.The mean and standard deviation of thickness were calculated as the mean and standard deviation of d over all sectors. (3) Radial coverage: the percentage of sectors positive for marker signal. A sector was considered positive if the mean signal over sector pixels acceded a user-defined threshold. (4) Jaggedness: this feature measures the extent jaggedness of an artery outline. To do so, first, a previously described skeletonization function⁶² is applied to the artery mask, and this function returns a ‘skeleton’ of the artery outline. This skeleton also assigns values to the outline pixels based on their distance from the core shape. Then, two different binarization thresholds are chosen: a non-branch threshold (a high value = 60 pixels, which indicates a greater topological distance) and a ‘branch’ threshold (a low value = 5 pixels, which indicates a smaller topological distance). The ratio between the total number of non-branch and branch pixels is the jaggedness.

Calculation of continuous SAR remodelling score δ

A supervised dimensionality reduction technique based on LDA⁶³ (https://github.com/davidrglass) was applied using the per-artery digitized morphological features and manually assigned remodelling stage labels as inputs. All artery morphology feature values were standardized (mean subtracted and divided by the standard deviation) and all arteries were used as the training data. The LDA output was as follows (Supplementary Table 3): the optimal linear combination of a subset of features that maximized the separation by manual stage between arteries in LDA space; and the coordinates of each artery in LDA space.

To define the SAR trajectory, a fourth-degree polynomial was fitted to the artery coordinates in LDA space. To determine the optimal degree of the polynomial, polynomials with degrees 1–6 were fitted, and the degree that minimized the P value for separating δ distributions between arteries grouped using the manual remodelling stage (Extended Data Fig. 4c) was selected. The polynomial fit was implemented using the MATLAB function fit and resulted in the following polynomial: f(x) = 0.0005 × x⁴ – 0.01227 × x³ + 0.1363 × x²– 0.4354 × x – 0.7425. The polynomial was then numerically interpolated on a dense 10⁴-point grid, and the distance from each artery point in LDA space to the polynomial was calculated using this grid and the MATLAB exchange function distance2curve⁶⁴. δ per artery was then calculated as the line integral from the curve origin to closest point to the artery on the curve (Extended Data Fig. 4a, inset). This integral was numerically calculated using a custom MATLAB script. δ values were linearly rescaled to the range 1–5 using the MATLAB function rescale.

Cell-type frequency as a function of GA and SAR

To examine cell-type frequencies within the decidua as a function of GA and SAR (Figs. 3 and 4), per-image cell frequency tables were constructed in which cell-type frequencies were calculated as the proportion of cells in the decidua feature mask of that image. Cells located in other feature masks (artery, gland, vessel or cell column masks) were not counted, nor were cells of an unassigned type (‘other’). To focus these analyses on cell populations strictly found in the decidua, muscle and glandular cells were also excluded; these cell types occasionally extended outside their artery and gland feature masks, respectively. Cell frequency as a function of GA for a cell type was defined as the per-image proportion values for that cell type, as a function of the GAs associated with the images. Similarly, cell frequency as a function of SAR for a cell type was defined as the per-image proportions of that cell type, as a function of the mean δ values per image. For the volcano plot in Fig. 3n, we fitted a linear regression model to the two above-described functions. All linear regression models were implemented using the MATLAB function fitlm and the volcano plot only shows points for which regression R² ≥ 0.05. R² and P values for all δ-based and GA-based regressions are provided in Supplementary Table 13. The ratio between R² in the two regression models was used to classify trends as GA-driven, SAR-driven or synchronized. For example, the increase in EVTs out of all cells, R_EVT, was classified as GA-driven because R² for R_EVT as a function of δ was 0.3 but only 0.1 for R_EVT as a function of GA (Extended Data Fig. 4d and Supplementary Table 10). Another example is the increase in macrophages out of immune cells, I_sumMac: it was classified as GA-driven because R² for I_sumMac as a function of GA was 0.6 but only 0.1 for I_sumMac as a function of δ (Extended Data Fig. 4e and Supplementary Table 10). To determine the trend sizes depicted in Fig. 3n, the following calculation was used: denote the per-image frequencies of a cell type as V, and the corresponding per image temporal stamps (either GA or mean image δ) as X. Trend size is then calculated as the difference between the first and last time point in units of the mean: \(\frac{V(\max (X))-V(\min (X))}(V)}\).

NanoString GeoMx DSP

The experiment was performed using NanoString Technologies according to company manuals, details are below.

Slide preparation

Serial sections of the TMAs were cut into 5 µm FFPE sections and were mounted on SuperFrost Plus slides (Fisher Scientific, 12-550-15), air dried and baked overnight at 60 °C. Slides were then processed as specified by the NanoString GeoMx DSP Slide Preparation User Manual (NanoString Technologies, MAN-100 7). In brief, slides were dewaxed, underwent antigen retrieval and treated with proteinase K (Ambion, 2546) at 1 µg ml^–1 concentration. Slides were then post-fixed. For RNA probe hybridization, slides were placed in a slide rack with Kimwipes damped with 2× SSC lining the bottom. Each slide was treated with 200 µl of NanoString Technologies whole transcriptome RNA probe mix at a concentration of 4 nM per probe in 1× buffer R (NanoString Technologies). A Hybridslip (Grace Biolabs, 714022) was applied over each slide. Slides were incubated at 37 °C overnight. After hybridization, slides were dipped in a 2× SSC with 0.1% Tween 20 (Teknova, T0710) to remove the coverslips. They were then washed twice in 2× SSC and 50% formamide (ThermoFisher AM9342) at 37 °C for 25 min followed by two washes in 2× SSC for 5 min each at room temperature. Slides were blocked in buffer W (NanoString Technologies) at room temperature for 30 min, followed by the application of 200 µl morphology marker mix for 1 h. Details of the morphology markers are provided in Supplementary Table 7.

Sample collection

Sample collection was performed as indicated in the GeoMx DSP instrument user manual (MAN-10088-03). Slides were loaded into the GeoMx DSP instrument and scanned. For each tissue sample, we selected ROIs corresponding to one of the following categories: artery (13), decidua (13), interstitial EVT (5), intravascular EVT (3); in total, 34 ROIs were selected (Supplementary Table 11). Morphology markers for SMA and VIM were used in conjunction with a serial H&E section to provide tissue context and to locate arteries and decidua on the platform. Artery, decidua and intravascular EVT ROIs were selected using the geometric selection tool, and interstitial EVTs were selected using a HLA-G⁺ mask. Intravascular EVTs were identified as HLA-G⁺ cells located within arteries. Each ROI was collected into a single well in a 96-well plate.

GeoMx DSP NGS library preparation and sequencing

Each GeoMx sample or well was uniquely indexed using an i5 × i7 dual-indexing system from Illumina. In total, 4 µl of a GeoMx DSP sample was used in a PCR reaction with 1 µM of i5 primer, 1 µM i7 primer and 1× NSTG PCR master mix. For the PCR amplification reaction, each 96-well plate was placed in a thermocycler programmed with the following protocol: 37 °C for 30 min, 50 °C for 10 min, 95 °C for 3 min, 18 cycles of 95 °C for 15 s, 65 °C for 60 s, 68 °C for 30 s, and final extension of 68 °C for 5 min. PCR assays were purified with two rounds of AMPure XP beads (Beckman Coulter) at 1.2× bead-to-sample ratio. Libraries were paired-end sequenced (2 × 75) on a NextSeq550 with up to 400 million total aligned reads.

Normalization and scaling of GeoMx counts data

Raw counts from each gene in each sample were extracted from the NanoString GeoMx NGS processing pipeline (Supplementary Table 11). Quality control was done according to the NanoString data analysis manual (MAN-10154-01) with default parameters as indicated in the manual. For each EVT sample, the counts were normalized using one of the manufacturer’s recommended approaches for normalizing GeoMx data: dividing all genes in each sample by the 75th percentile of expression in that sample, followed by multiplication by an identical scaling factor for all samples: the geometric mean of all 75th percentiles. This approach eliminates differences in counts between samples due to ROI-specific properties such as size and RNA-binding efficiency. The background due to nonspecific binding per sample was approximated with the geometric mean of the 100 negative control probes included in the probe mix, as recommended by NanoString Technologies. The above-described normalization step eliminated the correlation between background and ROI size for EVT samples. For artery and decidua samples, normalization was complicated by the fact that the ROI size was tightly correlated with SAR stage and therefore biologically meaningful trends in the data. This led to the correlation between ROI size and background not being entirely eliminated by normalization. We therefore used a background subtraction correction technique before normalization as recommended in the NanoString Technologies manual for such cases. The correction was performed by subtracting the geometric mean of negative probes from gene counts on a per-sample basis and proceeding with normalization as previously described.

Gene expression in artery as a function of GA and SAR

In brief, for each gene, we performed polynomial regressions of gene expression with δ and GA as the independent variables and used regression P values to determine which genes were trending and the ratio of regression R² values to classify the trends as detailed below.

The NanoString Technologies RNA probes panel contains probes for 18,696 transcripts. For this analysis on artery samples, only genes with background-subtracted, normalized counts ≥10 in at least two arteries were considered. This resulted in 14,471 expressed genes. Each artery sample was assigned a remodelling score δ based on the δ of the sampled artery in the MIBI data. If several arteries were sampled, the assigned δ was the average δ values of the sampled arteries. Endothelial loss and SMA loss per sample were calculated similarly based on the corresponding MIBI values (Supplementary Table 11). The following steps were then performed on artery samples.

For all expressed genes, gene expression as a function of GA was defined as the background-subtracted and normalized counts for that gene, as a function of the GAs associated with the samples. Similarly, expression as a function of SAR for a gene was defined as the per-sample background-subtracted and normalized counts of that gene, as a function of the δ values per sample. A second-degree polynomial regression model was then fitted to the two above-described functions. The reason for using a second-degree polynomial instead of linear regression was to allow the regression models to capture non-monotonic trends in gene expression. All regression models were implemented using the MATLAB function fitnlm. Expression fold change was defined as the ratio between the maximum and the minimum of expression values. The centre of mass (COM) of the expression trajectory of a gene as a function of t (t being either GA or δ) was defined as the weighted mean of t values, where the weights are the expression values at the respective t.

Genes with a P value ≤ 0.05 and fold change ≥2 for either GA or δ regression were classified as trending genes. The ratio between R² in the two regression models was used to classify trending genes as GA-driven, SAR-driven or synchronized. Trending genes with \({\log }_{2}({R}_{\delta }^{2}/{R}_}^{2})\ge 1\) and \({R}_{\delta }^{2}\ge 0.05\) were classified as SAR-driven, whereas genes with \({\log }_{2}({R}_{\delta }^{2}/{R}_}^{2})\le 1\) and \({R}_}^{2}\ge 0.05\) were classified as GA-driven. Other trending genes were classified as synchronized (Supplementary Table 4).

For visualization only, two fitted expression trajectories (one as a function of GA and another as a function of δ) were calculated per gene. These fitted expression trajectories were calculated as the values of the fitted second-degree polynomial model at five evenly spaced values of GA and δ, respectively. To compare fitted expression trajectories between genes, they were normalized by Z-scoring their value per gene (Fig. 3o).

See Supplementary Information for further details about analysis of NanoString data in decidua ROIs.

Coordinated gene expression by pathways in the artery

We set out to find gene pathways with coordinated expression trends among our genes of interest: genes trending with δ in arteries. To find these coordinated pathways, we first defined the pathways and then defined temporal coordination.

To define pathways, we used the R package msigdbr to obtain the lists of genes per pathway for the Gene Ontology by Biological Process database (7,481 pathways). We then cross-referenced the list of genes for each pathway with the genes of interest and discarded pathways with an intersection of fewer than ten genes. For the remaining pathways, we examined whether the pathway genes that appeared in the gene set of interest exhibited coordination in their expression as function of δ.

A group of coordinated genes was defined as a group of genes for which the COMs were significantly closer to each other than one would expect at random (see previous section for the definition of COM). Using the spread of COMs as a measure for coordination allowed us to leverage the raw data rather than fitted gene expression trajectories while still maintaining robustness against noise.

To calculate the extent of coordination between a group of N genes, we first calculated their median COM, denoted COM_med. Then, their COM dispersal was defined as the median of the absolute deviations from COM_med for the N genes, denoted CD. To determine whether the CD for the gene group, CD_group, is significantly smaller than expected at random, we calculated the randomly expected CD, denoted CD_rand. This was done by selecting N random genes without replacement and calculating their CD, 10⁵ times to estimate the null distribution. The random CD_rand was then calculated as the median over the CDs for randomized gene sets. The coordination score for our N genes group was then defined as log₂(CD_rand/CD_group). The P value for the coordination score was defined as the number of times a randomized CD was smaller than CD_group, divided by the number of randomizations (10⁵). (1/number of randomizations) was then added to all P values to account for the finite number of randomizations. q values were calculated using the Benjamini and Hochberg method on P values, implemented using the MATLAB function mafdr.

The CD, coordination scores, P values and q values were calculated as described above for all 7,481 pathways. Pathways with coordination score ≥ 1.5 and P value ≤ 0.05 were considered to be coordinated (Supplementary Table 4).

Ridge regression for predicting GA from immune composition

Ridge regression was implemented using the sklearn Python package (sklearn.linear_model.Ridge, RidgeCV). Per-image immune frequencies were rescaled to the range 0–1 before model fitting using the sklearn scaling function. Images with fewer than ten immune cells were excluded (n = 8). A randomly derived test–train split of 30/70 was used, and GA distribution was verified to be equally represented in the test and train sets (Extended Data Fig. 6a). Ridge regression adds a regularization penalty to the loss function to prevent over or under representation of correlated variables, such as immune cell populations. The penalty used for the test set (0.81) was selected using leave-one-out cross-validation on the training set.

Cell–cell and cell–artery spatial enrichment analysis

To identify preferential colocalization of maternal immune cells in decidua, we measured the spatial proximity enrichment for all cell-type pairs, which evaluates the spatial organization of cell types relative to each other, as previously described^22,23,24. Cells located in non-decidual feature masks (artery, gland, vessel or cell column masks) were not included in this analysis. The distances in pixels between all pairs of cells were calculated in each image. The resulting per-image distance matrices were binarized with a distance threshold (100 pixels or 39 μm in our case), and pairs of cells closer than 100 pixels from each other were considered a close interaction. To evaluate the number of close interactions between two cell types, this proximity matrix was subset column-wise by cell type A and subset row-wise by cell type B. The sum of the resulting submatrix quantified the number of close interactions between the cells of types A and B. To evaluate the significance of the number of close interactions, given the total number of cells in the image, tissue architecture and composition across the cohort, and total number of cells of types A and B in the image, a bootstrapping approach was used. For each of 100 bootstrapping iterations, the location of cells of type A was randomized across all cell locations (of any type) in the image while their total number was preserved. The number of close interactions with cells of type B was calculated for each randomized iteration. Repetitions of this process approached a null distribution for the number of close interactions between cells A and B. The enrichment score for cells A around cells B in the image was then calculated as the Z-score of the measured number of close interactions between A and B when Z-scored together with the random bootstraps. This analysis was extended to incorporate enrichment of cell types around spiral arteries. For each cell, the distance to the nearest spiral artery was considered. An additional column was added to the proximity matrix described above, which thresholded distances between cells and arteries with the same 100 pixel threshold. The above-described bootstrapping approach also provided a null distribution for artery proximity. Tools for this analysis were written in Python, with the bootstrapping accelerated using Cython. An intuitive, easy-to-use Jupyter Notebook interface was created to allow for easy implementation of this algorithm. For per-image spatial enrichment scores, see the Data availability statement. The code for this analysis is available at GitHub (https://github.com/angelolab/ark-analysis).

Cell–cell and cell–artery enrichment temporal trends and trending with GA or SAR or constant

For examining cell–cell and cell–artery enrichment within the decidua as a function of GA and SAR (Extended Data Figs. 6c and 2c), per-image enrichment score matrices E were calculated as described in the previous section, in which E_i,j is the enrichment score of cell type i around cell type j in the image. Enrichment as a function of GA was defined as the per-image enrichment, as a function of the GAs associated with the images. Similarly, enrichment as a function of SAR was defined as the per-image enrichment, as a function of the mean δ values per image. We fitted a linear regression model to the two above-described functions. All linear regression models were implemented using the MATLAB function fitlm. R² and P values for all δ-based and GA-based regressions are provided in Supplementary Table 3. The ratio between R² in the two regression models was used to classify trends as GA-driven, SAR-driven or synchronized like in Fig. 3n. Extended Data Fig. 6c only shows points for which regression R² ≥ 0.05, P value ≤ 0.05, maximal absolute value of linear fit ≥ 2. Trends including muscle, fibroblast, myofibroblast, glandular, other and endothelial cells were not considered in this analysis. For determining trend sizes, the following calculation was used: denote the linear fit to per-image enrichment scores as V, and the corresponding per-image temporal stamps (either GA or mean image δ) as X. Trend size is then calculated as (V(max(X)) – V(min(X))).

To determine whether two cell types were significantly enriched around each other throughout the cohort, we averaged their pairwise enrichment over all images. The pair was considered enriched if the absolute value of mean enrichment was ≥2 (Supplementary Table 3).

In Extended Data Fig. 6c, the following cell–cell enrichments were not plotted for clarity: enrichment of a cell type around itself (for example, T_reg cells around T_reg cells); enrichments including muscle, fibroblast, myofibroblast, glandular, other and endothelial cell types; enrichment trends that are SAR-driven.

Functional marker positivity rate per cell type as a function of GA and SAR

For examining cell-type-specific temporal trends in the expression of functional markers (Fig. 5a), 48 combinations of cell-type functional marker were selected. The selected combinations were those for which the positivity frequency Z-score exceeded 0.5 (Fig. 2a, right). For each of these combinations, the frequency of cells positive for the functional marker was calculated as the number of cells positive for the marker (see the methods section ‘Definition of thresholds for functional marker positivity’) out of the total number of cells of the same cell type in the image. All cells except those located within the cell column mask were included to focus the analysis on functional marker trends of maternal cells and EVTs that had infiltrated the decidua. For glandular cells, the location was further restricted to the glands mask. The frequency of cells positive for a functional marker as a function of GA, for a cell type, was defined as the per-image positivity proportion values as a function of the GAs associated with the images. Similarly, marker positivity frequency as a function of SAR for a cell type was defined as the per-image proportions of that cell type positive for the marker, as a function of the mean δ values per image. For the volcano plot in Fig. 5a, we fitted a linear regression model to the two above-described functions. All linear regression models were implemented using the MATLAB function fitlm, and the volcano plot only shows points for which regression R² > 0.05. R² and P values for all δ-based and GA-based regressions are provided in Supplementary Table 10. For determining trend sizes depicted in Fig. 5a, the following calculation was used: denote the linear fit to the per-image marker positivity proportion of a cell type as V, and the corresponding per-image temporal stamps (either GA or mean image δ) as X. Trend size is then calculated as the difference between the first and last time point in units of the mean: \(\frac{V({\rm{\max }}({\rm{X}}))-V({\rm{\min }}(X))}(V)}\).

Cellular microenvironments

For each cell in the dataset, we defined a ‘neighbourhood’ consisting of its 25 closest neighbouring cells, as measured by Euclidean distance between x/y centroids, excluding cells that were not in the decidua (that is, cells that overlapped with any artery, gland, anchoring villous or vessel feature masks). We clustered these cellular neighbourhoods on the basis of their composition of the 26 cell populations as identified previously using FlowSOM. For clustering, we used the scikit-learn implementation of k-means algorithm with k = 20 to identify neighbourhoods characterized by the presence of rare cell populations. Selected clusters were merged on the basis of similarity when hierarchically clustered, a threshold of 0.5 when comparing Euclidean distances between k-means cluster centroids and manual inspection of the cluster assignment when overlaid on the images. Based on these approaches, we defined 12 distinct decidual cellular microenvironments, 10 of which are shown in Fig. 5f (microenvironments characterized by predominantly stromal cell populations, fibroblasts (3 in Supplementary Table 12) and myofibroblasts (6 in Supplementary Table 12) are not shown in the heatmap).

Definition of anatomical EVT location and associated arteries

Cell column EVTs were defined as EVTs located within cell column masks, intravascular EVTs were located within artery masks, and interstitial EVTs were located in the decidua. Perivascular EVTs were defined as interstitial EVTs located within 50 pixels of the edge of an artery, as defined by the radial expansion of the artery masks (Extended Data Fig. 7b). Arteries were said to have perivascular or intravascular EVTs (Fig. 6d and Extended Data Fig. 7e,f) if the number of EVTs in the appropriate artery compartment was >5.

SMA and endothelium-loss scores

The loss scores presented in Extended Data Fig. 7e,f were based on digitized morphological features. For SMA, the average feature was used, whereas for endothelium, the radial coverage of CD31 was used (see methods section ‘Automated digitization of artery morphological features’). The values for each of the two features were then divided by their maximum across arteries and subtracted from 1 to obtain a loss score. The resulting values were then linearly rescaled to the range 0–1 using the MATLAB function rescale.

Characterization of EVTs by compartment

To further characterize EVT composition per spatial compartment (cell column, interstitial, perivascular or intravascular), we first quantified EVT-subtype frequency per compartment. One image was excluded (16_31762_20_8) owing to abnormal tissue morphology (Extended Data Fig. 8a). We then compared the distance from the nearest artery between EVT subtypes (Extended Data Fig. 8b). For this analysis, only images that contained all four EVT types were considered, and the cell-to-artery distance was measured from the cell centroid as detected by segmentation to the closest point on the border of the artery mask.

We then set out to assess the extent of similarity between EVTs by compartment in terms of expression of functional markers. For this analysis we used a LDA based method, similarly to our calculation of the continuous SAR remodelling score δ for compartment-wise analysis of EVT types. The input table for LDA consisted of MIBI-measured marker expression values per EVT. The following lineage and functional markers expressed by EVTs were included: CD56, CD57, HLA-G, CK7, PD-L1 and Ki67. EVTs were labelled by spatial compartment as cell column, interstitial, perivascular or intravascular. Marker expression values were standardized (mean subtracted and divided by the standard deviation), and cell column, interstitial and intravascular location labels per EVT were used for training the LDA model. Perivascular EVTs were withheld as a test set. Owing to the small number or features (markers), a one-dimensional LDA was calculated to produce a single coordinate, LD1. LD1 was the optimal linear combination of a subset of markers to maximize the separation by compartment between EVTs (Supplementary Table 14). LD1 values were subsequently calculated for the withheld test set of perivascular EVTs (Supplementary Table 14). The distributions of LD1 values per compartment indicated that perivascular EVTs are similar to interstitial and intravascular EVTs, with a median value between the two (Extended Data Fig. 8c). This result implies that perivascular EVTs could be an intermediate state between interstitial and intravascular, in line with the intravasation model whereby interstitial EVTs invade the artery lumen.

Origin of CD56⁺ EVTs in the intravascular compartment

The frequency of CD56⁺ EVTs was highest in the intravascular compartment (Extended Data Fig. 8a). Furthermore, the frequency of CD56⁺ EVTs increased with SAR both in the perivascular and intravascular compartment (Extended Data Fig. 8d,e). However, the increase in the intravascular compartment was steeper (Extended Data Fig. 8f). We set out to test whether this was compatible with intravasation whereby the source for intravascular EVTs is perivascular EVTs or whether additional sources of intravascular CD56⁺ EVTs were needed to account for their increase. These alternative sources could be cell proliferation or extravasation, whereby EVTs migrate in a retrograde manner after entering spiral arteries directly at the basal plate. Under the intravasation model, the steep increase in CD56⁺ EVTs between the perivascular and intravascular compartments would be explained by CD56 upregulation upon arterial invasion. We proposed that this should involve an intermediate state of perivascular EVTs defined by moderate levels of CD56. To test this hypothesis, we compared the average CD56 intensity of perivascular and intravascular EVT1a and EVT1b cells for each artery (for arteries that initiated remodelling: δ ≥ 2). This analysis detected a significant increase in CD56 expression between the perivascular and intravascular compartment by EVT1a and EVT1b cells (sided Wilcoxon signed rank test P value = 5 × 10⁻³; Extended Data Fig. 8g). An alternative explanation for the increasing frequency of CD56⁺ EVT1c cells within arteries could be proliferation. However, only 0.5% of intravascular EVT1c cells were Ki67⁺, a lower frequency than 9.6% and 1.8% for intravascular EVT1a and EVT1b cells, respectively, which suggested that proliferation is not a primary contributor (Extended Data Fig. 8h).

DEGs in EVTs

DEGs between intravascular and interstitial EVTs were identified using the Bioconductor package limma⁶⁵ (linear models for microarray data) after consulting with the NanoString statistics team. Using the default parameters in limma on 75th percentile normalized counts, 131 upregulated genes and 143 downregulated genes were found (false-discovery rate cutoff = 0.1, log fold change cutoff = 2). Genes with log fold change ≥ 2.3 or ≤ 2.3 and adjusted P value ≤ 0.05 are shown in the heatmap in Fig. 6f. A complete list of DEGs is shown in Extended Data Fig. 9 and Supplementary Table 11. For IHC validation of differential expression for selected genes, see Extended Data Fig. 10a.

NicheNet analysis

We used the NicheNet R package to predict ligand–receptor interactions between intravascular EVTs and arteries. The analysis was performed by following the vignette available at GitHub (https://github.com/saeyslab/nichenetr/blob/master/vignettes/ligand_activity_geneset.md).

NicheNet requires three input gene lists to predict ligands in sender cells that are likely to interact with receptors in receiver cells and by doing so affect the expression of a set of genes of interest. These three gene lists are: genes of interest, genes expressed in sender cells and genes expressed in receiver cells.

For our analysis, we wanted to check which ligands expressed in intravascular EVTs are likely to be causing temporal gene expression trends with remodelling in arteries. To do so, we defined the genes of interest as all genes trending with remodelling in arteries (Fig. 3o). The genes expressed in receiver cells were defined as all genes expressed in arteries (see previous sections for the definition of ‘expressed’), and genes expressed in sender cells were defined as genes differentially expressed between interstitial and intravascular EVTs and higher in intravascular EVTs (Supplementary Table 11). NicheNet analysis was performed as described in the vignette to prioritize ligands and to infer corresponding receptors and downstream targets (Extended Data Fig. 10b). The inferred targets were manually classified according to their known function using the Gene Cards database (https://www.genecards.org) and survey of the literature. A list of references for all classifications are provided in Supplementary Table 5.

Statistical analyses

Throughout the paper, unless indicated otherwise, the Kruskal–Wallis test was used. It was implemented using the MATLAB function KruskalWallis. All linear regression models were implemented using the MATLAB function fitlm unless stated otherwise. The sided Wilcoxon signed-rank test for paired analysis was implemented using the MATLAB function signrank. MATLAB v.2020b was used throughout the article for statistical analysis.

Ethics statement

All human samples were acquired and all experiments were approved by Institutional Review Board protocol number 46646 “Assessing normal expression patterns of immune and non-immune markers across tissue types with multiplexed ion beam imaging” at Stanford University. Per this protocol from the Institutional Review Board at Stanford University, the consent to use archival deidentified tissue was not required. All experiments followed all relevant guidelines and regulations.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability

MIBI data are available at the Human BioMolecular Atlas Program (https://doi.org/10.35079/hbm585.qpdv.454). The same MIBI data in a browsable format, along with segmentation masks, extracted features, cell phenotype masks (CPMs), cell–cell and cell–artery spatial-enrichment scores per image, a table enumerating all single cells in this study and provides their location, morphological characteristics (such as size and shape), marker expression, FlowSOM cluster assignment and cell-type assignment, are available at FigShare (https://doi.org/10.6084/m9.figshare.16663465). H&E images of tissue blocks from the cohort with annotations are available at Dryad (https://doi.org/10.5061/dryad.v15dv41zp). The following public databases were used in this study (see Methods for details): the Gene Ontology Biological Process database (http://geneontology.org/) and the Gene Cards database (https://www.genecards.org/). Source data are provided with this paper.

Code availability

The following software were used for analysis in the article: ImageJ, MAUI (https://github.com/angelolab/MAUI) for low-level image processing; DeepCell v.0.5.0 (https://deepcell.readthedocs.io/en/master/index.html) for cell segmentation; Ark Analysis for cell–cell spatial enrichment (https://github.com/angelolab/ark-analysis); and QuPath (v.0.4.0) for two-colour IHC analysis. Custom code for this study is available at FigShare (https://doi.org/10.6084/m9.figshare.16663465).

References

Menkhorst, E., Winship, A., van Sinderen, M. & Dimitriadis, E. Human extravillous trophoblast invasion: intrinsic and extrinsic regulation. Reprod. Fertil. Dev. 28, 406–415 (2016).

Article CAS PubMed Google Scholar
Vento-Tormo, R. et al. Single-cell reconstruction of the early maternal–fetal interface in humans. Nature 563, 347–353 (2018).

Article ADS CAS PubMed PubMed Central Google Scholar
Pijnenborg, R., Vercruysse, L. & Vercruysse, M. The uterine spiral arteries in human pregnancy: facts and controversies. Placenta 27, 939–958 (2006).

Article CAS PubMed Google Scholar
Burton, G. J., Woods, A. W., Jauniaux, E. & Kingdom, J. C. P. Rheological and physiological consequences of conversion of the maternal spiral arteries for uteroplacental blood flow during human pregnancy. Placenta 30, 473–482 (2009).

Article CAS PubMed PubMed Central Google Scholar
Brosens, I., Pijnenborg, R., Vercruysse, L. & Romero, R. The “Great Obstetrical Syndromes” are associated with disorders of deep placentation. Am. J. Obstet. Gynecol. 204, 193 (2011).

Article PubMed Google Scholar
Carter, A. M. Animal models of human placentation—a review. Placenta 28, S41–S47 (2007).

Article PubMed Google Scholar
Schmidt, A., Morales-Prieto, D. M., Pastuschek, J., Fröhlich, K. & Markert, U. R. Only humans have human placentas: molecular differences between mice and humans. J. Reprod. Immunol. 108, 65–71 (2015).

Article CAS PubMed Google Scholar
Gaynor, L. M. & Colucci, F. Uterine natural killer cells: functional distinctions and influence on pregnancy in humans and mice. Front. Immunol. 0, 467 (2017).

Article Google Scholar
Rockwell, C., Vargas, E. & Moore, L. G. Human physiological adaptation to pregnancy: inter- and intraspecific perspectives. Am. J. Hum. Biol. 15, 330–341 (2003).

Article PubMed Google Scholar
Moffet, A. & Colucci, F. Uterine NK cells: active regulators at the maternal–fetal interface. J. Clin. Invest. 124, 1872–1879 (2014).

Article Google Scholar
Angelo, M. et al. Multiplexed ion beam imaging of human breast tumors. Nat. Med. 20, 436–442 (2014).

Article CAS PubMed PubMed Central Google Scholar
Yang, F., Zheng, Q. & Jin, L. Dynamic function and composition changes of immune cells during normal and pathological pregnancy at the maternal–fetal interface. Front. Immunol. 10, 2317 (2019).

Article CAS PubMed PubMed Central Google Scholar
Chang, R.-Q., Li, D.-J. & Li, M.-Q. The role of indoleamine-2,3-dioxygenase in normal and pathological pregnancies. Am. J. Reprod. Immunol. 79, e12786 (2018).

Article PubMed Google Scholar
Li, M. et al. Tim-3/CTLA-4 pathways regulate decidual immune cells–extravillous trophoblasts interaction by IL-4 and IL-10. FASEB J. 35, e21754 (2021).

Article CAS PubMed Google Scholar
Hu, X. et al. Newly characterized decidual Tim-3⁺ Treg cells are abundant during early pregnancy and driven by IL-27 coordinately with Gal-9 from trophoblasts. Hum. Reprod. 35, 2454–2466 (2020).

Article CAS PubMed PubMed Central Google Scholar
Krop, J. et al. Imaging mass cytometry reveals the prominent role of myeloid cells at the maternal–fetal interface. iScience 25, 104648 (2022).

Article ADS CAS PubMed PubMed Central Google Scholar
Kämmerer, U. et al. Unique appearance of proliferating antigen-presenting cells expressing DC-SIGN (CD209) in the decidua of early human pregnancy. Am. J. Pathol. 162, 887–896 (2003).

Article PubMed PubMed Central Google Scholar
Nielsen, C. M., White, M. J., Goodier, M. R. & Riley, E. M. Functional significance of CD57 expression on human NK cells and relevance to disease. Front. Immunol. https://doi.org/10.3389/fimmu.2013.00422 (2013).
Harjunpää, H. & Guillerey, C. TIGIT as an emerging immune checkpoint. Clin. Exp. Immunol. 200, 108–119 (2020).

Article PubMed Google Scholar
Huhn, O. et al. Distinctive phenotypes and functions of innate lymphoid cells in human decidua during early pregnancy. Nat. Commun. 11, 381 (2020).

Article ADS CAS PubMed PubMed Central Google Scholar
Smith, S. D., Dunk, C. E., Aplin, J. D., Harris, L. K. & Jones, R. L. Evidence for immune cell involvement in decidual spiral arteriole remodeling in early human pregnancy. Am. J. Pathol. 174, 1959–1971 (2009).

Article CAS PubMed PubMed Central Google Scholar
Keren, L. et al. A structured tumor–immune microenvironment in triple negative breast cancer revealed by multiplexed ion beam imaging. Cell 174, 1373–1387.e19 (2018).

Article CAS PubMed PubMed Central Google Scholar
McCaffrey, E. F. et al. The immunoregulatory landscape of human tuberculosis granulomas. Nat. Immunol. 23, 318–329 (2022).

Article CAS PubMed PubMed Central Google Scholar
Risom, T. et al. Transition to invasive breast cancer is associated with progressive changes in the structure and composition of tumor stroma. Cell 185, 299–310.e18 (2022).

Article CAS PubMed PubMed Central Google Scholar
Sedlmayr, P. et al. Localization of indoleamine 2,3-dioxygenase in human female reproductive organs and the placenta. Mol. Hum. Reprod. 8, 385–391 (2002).
Sharma, S., Godbole, G. & Modi, D. Decidual control of trophoblast invasion. Am. J. Reprod. Immunol. 75, 341–350 (2016).

Article CAS PubMed Google Scholar
Silva, I. G. et al. The Tim-3–galectin-9 secretory pathway is involved in the immune escape of human acute myeloid leukemia cells. eBioMedicine 22, 44–57 (2017).

Article Google Scholar
Pijnenborg, R., Bland, J. M., Robertson, W. B. & Brosens, I. Uteroplacental arterial changes related to interstitial trophoblast migration in early human pregnancy. Placenta 4, 397–413 (1983).

Article CAS PubMed Google Scholar
van Gassen, S. et al. FlowSOM: using self-organizing maps for visualization and interpretation of cytometry data. Cytometry A 87, 636–645 (2015).

Article PubMed Google Scholar
Founds, S. A., Fallert-Junecko, B., Reinhart, T. A. & Parks, W. T. LAIR2-expressing extravillous trophoblasts associate with maternal spiral arterioles undergoing physiologic conversion. Placenta 34, 248–255 (2013).

Article CAS PubMed Google Scholar
Hunkapiller, N. M. et al. A role for Notch signaling in trophoblast endovascular invasion and in the pathogenesis of pre-eclampsia. Development 138, 2987–2998 (2011).

Article CAS PubMed PubMed Central Google Scholar
Zhang, P. Phenotypic switch of endovascular trophoblasts in decidual vasculopathy with implication for preeclampsia and other pregnancy complications. Fetal Pediatr. Pathol. 40, 543–562 (2021).

Article PubMed Google Scholar
Dorris, E. R. et al. The autoimmune susceptibility gene C5orf30 regulates macrophage-mediated resolution of inflammation. J. Immunol. 202, 1069–1078 (2019).

Article CAS PubMed Google Scholar
Kano, J., Wang, H., Zhang, H. & Noguchi, M. Roles of DKK3 in cellular adhesion, motility, and invasion through extracellular interaction with TGFBI. FEBS J. https://doi.org/10.1111/FEBS.16529 (2022).
Sammar, M. et al. Expression of CD24 and Siglec-10 in first trimester placenta: implications for immune tolerance at the fetal–maternal interface. Histochem. Cell Biol. 147, 565–574 (2017).

Article CAS PubMed Google Scholar
Barkal, A. A. et al. CD24 signalling through macrophage Siglec-10 is a target for cancer immunotherapy. Nature 572, 392–396 (2019).
Browaeys, R., Saelens, W. & Saeys, Y. NicheNet: modeling intercellular communication by linking ligands to target genes. Nat. Methods 17, 159–162 (2020).

Article CAS PubMed Google Scholar
Northcott, J. M., Czubryt, M. P. & Wigle, J. T. Vascular senescence and ageing: a role for the MEOX proteins in promoting endothelial dysfunction. Can. J. Physiol. Pharmacol. 95, 1067–1077 (2017).

Article CAS PubMed Google Scholar
Ma, H. et al. HMBOX1 interacts with MT2A to regulate autophagy and apoptosis in vascular endothelial cells. Sci. Rep. 5, 15121 (2015).

Article ADS CAS PubMed PubMed Central Google Scholar
Pluznick, J. L. et al. Olfactory receptor responding to gut microbiotaderived signals plays a role in renin secretion and blood pressure regulation. Proc. Natl Acad. Sci. USA 110, 4410–4415 (2013).

Article ADS CAS PubMed PubMed Central Google Scholar
Kim, S. H. et al. Expression of human olfactory receptor 10J5 in heart aorta, coronary artery, and endothelial cells and its functional role in angiogenesis. Biochem. Biophys. Res. Commun. 460, 404–408 (2015).

Article CAS PubMed Google Scholar
Meekins, J. W., Pijnenborg, R., Hanssens, M., MCFadyen, I. R. & van Asshe, A. A study of placental bed spiral arteries and trophoblast invasion in normal and severe pre-eclamptic pregnancies. Br. J. Obstet. Gynaecol. 101, 669–674 (1994).

Article CAS PubMed Google Scholar
Tong, J. et al. Transcriptomic profiling in human decidua of severe preeclampsia detected by RNA sequencing. J. Cell. Biochem. 119, 607–615 (2018).

Article CAS PubMed Google Scholar
Tarca, A. L. et al. Maternal whole blood mRNA signatures identify women at risk of early preeclampsia: a longitudinal study. J. Matern. Fetal Neonatal Med. 34, 3463–3474 (2021).

Article CAS PubMed Google Scholar
Luttun, A. & Carmeliet, P. Soluble VEGF receptor FLT1: the elusive preeclampsia factor discovered? J. Clin. Invest. 111, 600–602 (2003).

Article CAS PubMed PubMed Central Google Scholar
Takahashi, H. et al. Extravillous trophoblast invasion accelerated by WNT3A, 5A, and 10B via CD44. J. Matern. Fetal Neonatal Med. 34, 3377–3385 (2021).

Article CAS PubMed Google Scholar
Moufarrej, M. N. et al. Early prediction of preeclampsia in pregnancy with cell-free RNA. Nature 602, 689–694 (2022).

Article ADS CAS PubMed PubMed Central Google Scholar
von Rango, U. et al. Apoptosis of extravillous trophoblast cells limits the trophoblast invasion in uterine but not in tubal pregnancy during first trimester. Placenta 24, 929–940 (2003).

Article Google Scholar
Wolf, K. et al. Physical limits of cell migration: control by ECM space and nuclear deformation and tuning by proteolysis and traction force. J. Cell Biol. 201, 1069–1084 (2013).

Article CAS PubMed PubMed Central Google Scholar
Naruse, K. et al. Localization of matrix metalloproteinase (MMP)-2, MMP-9 and tissue inhibitors for MMPs (TIMPs) in uterine natural killer cells in early human pregnancy. Hum. Reprod. 24, 553–561 (2009).

Article CAS PubMed Google Scholar
Hunt, J., Petroff, M., McIntire, R. & Ober, C. HLA-G and immune tolerance in pregnancy. FASEB J. 19, 681–693 (2005).

Article CAS PubMed Google Scholar
Tilburgs, T., Evans, J. H., Crespo, Â. & Strominger, J. The HLA-G cycle provides for both NK tolerance and immunity at the maternal–fetal interface. Proc. Natl Acad. Sci. USA 112, 13312–13317 (2015).

Article ADS CAS PubMed PubMed Central Google Scholar
Ramsey, E. M. & Donner, M. W. Placental vasculature and circulation in primates. Placent. Vasculariz. Blood Flow 3, 217–233 (1988).

Article Google Scholar
Kshitiz et al. Evolution of placental invasion and cancer metastasis are causally linked. Nat. Ecol. Evol. 3, 1743–1753 (2019).

Article CAS PubMed Google Scholar
Bulmer, J. N., Thrower, S. & Wells, M. Expression of epidermal growth factor receptor and transferrin receptor by human trophoblast populations. Am. J. Reprod. Immunol. 21, 87–93 (1989).

Article CAS PubMed Google Scholar
Wright, J. K. et al. HER1 signaling mediates extravillous trophoblast differentiation in humans. Biol. Reprod. 83, 1036–1045 (2010).

Article CAS PubMed Google Scholar
Coorens, T. H. H. et al. Inherent mosaicism and extensive mutation of human placentas. Nature 592, 80–85 (2021).

Article ADS CAS PubMed PubMed Central Google Scholar
Keren, L. et al. MIBI-TOF: a multiplexed imaging platform relates cellular phenotypes and tissue structure. Sci. Adv. 5, eaax5851 (2019).

Article ADS CAS PubMed PubMed Central Google Scholar
Baranski, A. et al. MAUI (MBI Analysis User Interface)—an image processing pipeline for multiplexed mass based imaging. PLoS Comput. Biol. 17, e1008887 (2021).

Article CAS PubMed PubMed Central Google Scholar
Greenwald, N. F. et al. Whole-cell segmentation of tissue images with human-level performance using large-scale data annotation and deep learning. Nat. Biotechnol. 40, 555–565 (2022).

Article CAS PubMed Google Scholar
Spence, A. et al. Dual antibody immunohistochemistry staining V.1. protocols.io https://doi.org/10.17504/protocols.io.81wgbyoryvpk/v1 (2023).
Howe, N. Better Skeletonization version 1.0.0.0. MatLab Central File Exchange https://www.mathworks.com/matlabcentral/fileexchange/11123-better-skeletonization (2023).
Tsai, A. G. et al. Multiplexed single-cell morphometry for hematopathology diagnostics. Nat. Med. 26, 408–417 (2020).

Article CAS PubMed PubMed Central Google Scholar
D’Errico, J. distance2curve version 1.1.0.0. MatLab Central File Exchange https://www.mathworks.com/matlabcentral/fileexchange/34869-distance2curve (2023).
Ritchie, M. E. et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43, e47 (2015).

Article PubMed PubMed Central Google Scholar

Download references

Acknowledgements

We thank M. Amouzgar, C. Liu, N. Vivanco, A. Moore, E. McCaffrey, D. Glass, T. Risom, J. P. Oliveria, K. O'Neill and C. Coutifaris for comments. M.A. is supported by 5U54CA20997105, 5DP5OD01982205, 1R01CA24063801A1, 5R01AG06827902, 5UH3CA24663303, 5R01CA22952904, 1U24CA22430901, 5R01AG05791504 and 5R01AG05628705 from the NIH, W81XWH2110143 from the DOD, and other funding from the Bill and Malinda Gates Foundation, the Cancer Research Institute, the Parker Center for Cancer Immunotherapy, and the Breast Cancer Research Foundation. S.G. is supported by the Bill and Melinda Gates Foundation (OPP1113682). I.A. is an awardee of the Weizmann Institute of Science–Israel National Postdoctoral Award Program for Advancing Women in Science. E.S. is supported by National Science Scholarship, Agency for Science, Technology, and Research (A*STAR), Singapore. Illustrations in Figs. 1a and 3a, parts of Figs. 3m and 6g and Extended Data Fig. 8a,b were commissioned by SciStories.

Author information

These authors contributed equally: Shirley Greenbaum, Inna Averbukh, Erin Soon

Authors and Affiliations

Department of Pathology, Stanford University, Stanford, CA, USA

Shirley Greenbaum, Inna Averbukh, Erin Soon, Alex Baranski, Noah F. Greenwald, Adam Kagel, Marc Bosse, Zumana Khair, Shirley Kwok, Shiri Warshawsky, Hadeesha Piyadasa, Mako Goldston, Angie Spence, Matt van de Rijn & Michael Angelo
Department of Obstetrics and Gynecology, Hadassah-Hebrew University Medical Center, Jerusalem, Israel

Shirley Greenbaum
Immunology Program, Stanford University, Stanford, CA, USA

Erin Soon
Department of Pathology, University of Californica San Francisco, San Francisco, CA, USA

Gabrielle Rizzuto
Cancer Biology Program, Stanford University, Stanford, CA, USA

Noah F. Greenwald
Department of Obstetrics Gynecology and Reproductive Sciences, University of California San Francisco, San Francisco, CA, USA

Eleni G. Jaswa
Division of Biology and Bioengineering, California Institute of Technology, Pasadena, CA, USA

Geneva Miller, Morgan Schwartz, Will Graf & David Van Valen
Department of Obstetrics and Gynecology, Stanford University, Stanford, CA, USA

Virginia D. Winn
Department of Pathology, Memorial Sloan Kettering Cancer Center, New York, NY, USA

Travis Hollmann
Department of Molecular Cell Biology, Weizmann Institute of Science, Rehovot, Israel

Leeat Keren

Authors

Shirley Greenbaum

View author publications

You can also search for this author in PubMed Google Scholar
Inna Averbukh

View author publications

You can also search for this author in PubMed Google Scholar
Erin Soon

View author publications

You can also search for this author in PubMed Google Scholar
Gabrielle Rizzuto

View author publications

You can also search for this author in PubMed Google Scholar
Alex Baranski

View author publications

You can also search for this author in PubMed Google Scholar
Noah F. Greenwald

View author publications

You can also search for this author in PubMed Google Scholar
Adam Kagel

View author publications

You can also search for this author in PubMed Google Scholar
Marc Bosse

View author publications

You can also search for this author in PubMed Google Scholar
Eleni G. Jaswa

View author publications

You can also search for this author in PubMed Google Scholar
Zumana Khair

View author publications

You can also search for this author in PubMed Google Scholar
Shirley Kwok

View author publications

You can also search for this author in PubMed Google Scholar
Shiri Warshawsky

View author publications

You can also search for this author in PubMed Google Scholar
Hadeesha Piyadasa

View author publications

You can also search for this author in PubMed Google Scholar
Mako Goldston

View author publications

You can also search for this author in PubMed Google Scholar
Angie Spence

View author publications

You can also search for this author in PubMed Google Scholar
Geneva Miller

View author publications

You can also search for this author in PubMed Google Scholar
Morgan Schwartz

View author publications

You can also search for this author in PubMed Google Scholar
Will Graf

View author publications

You can also search for this author in PubMed Google Scholar
David Van Valen

View author publications

You can also search for this author in PubMed Google Scholar
Virginia D. Winn

View author publications

You can also search for this author in PubMed Google Scholar
Travis Hollmann

View author publications

You can also search for this author in PubMed Google Scholar
Leeat Keren

View author publications

You can also search for this author in PubMed Google Scholar
Matt van de Rijn

View author publications

You can also search for this author in PubMed Google Scholar
Michael Angelo

View author publications

You can also search for this author in PubMed Google Scholar

Contributions

S.G. assembled the tissue cohort, performed tissue staining and MIBI experiments, annotated images, analysed and interpreted MIBI data and wrote the manuscript. I.A. analysed and interpreted MIBI data and NanoString data and wrote the manuscript. E.S. assisted with tissue staining and MIBI experiments, annotated images, performed NanoString experiments, analysed and interpreted MIBI and NanoString data and wrote the manuscript. G.R. advised on the cohort design, assembled the tissue cohort and annotated images. A.B. wrote software for artery digital morphometrics. N.F.G., G.M., M.S., W.G. and D.V.V. developed and wrote software for cell segmentation and clustering. A.K. wrote software for cell–cell spatial-enrichment measurement. M.B. performed two-colour IHC experiments, advised on experimental design and reagent validation. M.G. and A.S. performed IHC experiments. H.P. and V.D.W. advised on experimental design and reagent validation. E.G.J. assembled cohort patient metadata. L.K. advised on computational analyses. Z.K. prepared and validated reagents. S.K. constructed the TMAs. S.W. annotated images. T.H. validated reagents and advised on experimental design. M.v.d.R. oversaw TMA construction. M.A. conceived the study, advised on experimental design and data analyses, interpreted data and wrote the manuscript.

Corresponding authors

Correspondence to Shirley Greenbaum or Michael Angelo.

Ethics declarations

Competing interests

M.A. and is a named inventor on patent US20150287578A1, which covers the mass spectrometry approach utilized by MIBI-TOF to detect elemental reporters in tissue using secondary ion mass spectrometry. M.A. is a board member and shareholder in IonPath, which develops and manufactures the commercial MIBI-TOF platform. All other authors declare no competing interests.

Peer review

Peer review information

Nature thanks Anniina Farkkila and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

Additional information

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Extended data figures and tables

Extended Data Fig. 1 Representative H&E and MIBI images for cell phenotyping.

a. Serial sections with cell phenotype maps. Each cell colored by its phenotype assignment for GAs: 6, 12, 20 weeks. VIM, vimentin; SMA, smooth muscle actin; Placental Mac, Placental Macrophage; MyoFib, Myofibroblast. Representative images of n = 55 FOVs (6 weeks), n = 12 FOVs (12 weeks), n = 33 FOVs (20 weeks). Scale bar, 100 µm.

Extended Data Fig. 2 Validation of CD57+ NK2 population.

a. Bar plot of double positive cells (CD57+CD49a+) percent out of CD57+ cells per patient. b. Representative examples of double positive (CD57+CD49a+, left) and single positive (CD57+CD49a-, right) staining using dual color IHC for CD57 and CD49a on decidua samples, 595 cells in total. Representative images of n = 574 double positive cells and n = 21 single positive cells. Scale bar, 20 µm. c. Pairwise enrichment between NK2 cells and arteries. Left: mean enrichment per image by δ, n = 209 images. Red line – fitted linear regression, regression p-value. Right: Violin plot of NK2-artery enrichment scores by early (1 ≤ δ < 3) and late (3 ≤ δ < 5) remodeling stage. p-value Kruskal-Wallis. Early: n = 108, min = −3.6, max= 7.07, center = 0.13; Late: n = 15, min = −1.77, max=0.41, center = −0.45.

Source Data

Extended Data Fig. 3 SAR manual and digitized staging.

a. Flowchart for manual staging of arteries. b. Bar plot of manually-determined artery stage frequencies per patient. c. Schematic of digitization of artery morphology features. d. Heatmap of quantified digitized morphological features, mean across SAR stages. Columns Z-scored and hierarchically clustered.

Source Data

Extended Data Fig. 4 SAR trajectory.

a. Scatter of arteries in LDA space colored by manually-assigned stage. Polynomial fit: remodeling trajectory. Inset: matching each artery point a_i to the SAR trajectory by finding the nearest point along trajectory b_i. The continuous SAR score δ was then defined as the distance from origin x₀ to b_i along the trajectory curve. b. Scatter plot of arteries in LDA space colored by δ. c. Violins of artery δs by manual stage d,e. Scatter plots by δ (left) and GA (right) of EVT frequency out of all cells (d), macrophage frequency out of immune cells (e) Red lines, fitted linear regression, regression p-values on top.

Source Data

Extended Data Fig. 5 SAR-trending GO pathways.

a. 3 SAR (δ)-trending GO pathways, showing normalized expression of genes in the pathway by SAR (δ): mean over genes +/− SEM.

Source Data

Extended Data Fig. 6 Ridge regression for predicting GA.

a. Distribution of GA (in days) across the whole dataset, training (70%) set, and test (30%) set. b. Predicted versus actual δ for ridge regression trained on GA-associated immune features, for the withheld test set (30%). Line; best fit, shaded region; 1 SD. c. Pairwise spatial enrichment relationships, including: temporally coordinated with GA and GA/SAR synchronized. Trend size color coded: red indicates increasing spatial pairwise enrichment, blue indicates decreasing enrichment. Arrow length represents mean spatial enrichment (Z-score) across the entire dataset. Node size represents its number of pairwise relationships.

Source Data

Extended Data Fig. 7 EVT distribution.

a. Two hypotheses for intravascular EVT invasion. (Left) Intravasation: orange arrows indicate movements of EVTs from the cell column of the anchoring villi into the decidua and through the wall of the artery into the lumen. (Right) Extravasation: red arrows indicate movement of EVTs from the fetal villi through the intervillous space into the artery. b. Anatomical locations of EVT populations in the decidua. c. Distribution of artery δs, based on the presence of perivascular and/or intravascular EVTs. [Violins left to right] min = 1, max = 2.65, center = 1.7; min = 1.19, max = 3, center = 1.84; min = 1.77, max = 5, center = 3.2; min = 1.6, max = 4.5, center = 2.4. d. Scatter plot of log₂(Intravascular/Perivascular) ratio by δ, for arteries with both perivascular and intravascular EVTs present. Black line, fitted linear regression. Regression p-value on top. e. Percentage of arteries with ≤ a given SMA loss (s) threshold, by perivascular or intravascular EVTs present. Arteries were considered to have perivascular or intravascular EVT if the number of EVT in the appropriate compartment was ≥ 5. f. Percentage of arteries with scores ≤ a given endothelial loss (e) threshold, by perivascular or intravascular EVTs present. Arteries were considered to have perivascular or intravascular EVT if the number of EVT in the appropriate artery compartment was ≥ 5.

Source Data

Extended Data Fig. 8 EVT Phenotype by Location.

a. Frequency of EVT populations by anatomical location. b. Violin plot of distance from artery (in pixels) by EVT type. EVT1c: N = 209, min = 12.8, max = 1536.7, center (median) = 173.9; EVT1a: N = 7123, min = 10, max = 2028.7, center = 363.6; EVT1b: N = 5908, min = 7.2, max = 1989.5, center = 376.45; EVT1c: N = 185, min= 9.4, max. = 1818.1, center = 475.6. c. Violin plot of the distribution of LD1 for EVTs, by anatomical location. Center lines inside violins indicate mean. Anchoring: N = 8906, min = −3.6, max = 3.9, center = −1; Interstitial: N = 38395, min = −3.4, max = 4.4, center = 0; Perivascular: N = 1097, min = −3.14, max = 4, center = 0.08; Intravascular: N=4040, min = −3.15, max = 4.6, center = 2. d-h: d,e Scatter plots of perivascular (d) and intravascular (e) EVT1c (CD56+) frequency by δ. Red lines, fitted linear regression, regression p-values on top. f. Bar plot of the EVT1c frequency increase rates (regression slopes from d,e). Error bars, 95% confidence interval for regression slopes. g. Paired-by-artery CD56 expression in EVT1a&b, between perivascular and intravascular compartments. Arteries with δ ≥ 2 included. p = 5e-03, one sided paired Wilcoxon signed rank test, z-statistic = 2.5 h. Proportion of Ki67+ intravascular EVT by type.

Source Data

Extended Data Fig. 9 Differential gene in EVT.

a. Full heatmap for differentially expressed genes between intravascular and interstitial EVTs showing gene expression (Z-score), (logFC > 2, adj p-value <0.05 using Limma).

Source Data

Extended Data Fig. 10 EVT- artery interactions.

a. IHC validation on decidua serial sections for protein counterparts for 3 differentially expressed genes (JAG1, C5ORF30, and EBI3) with controls: HLA-G, CD56, and H&E. 1 section per marker. Scale bar, 100 µm. b. Full NicheNet output: 10 EVT-ligands best predicting receivers expressed in arteries, ranked by Pearson correlation coefficient or the EVT ligand activity ranking metric. Genes differentially expressed in Preeclamptic decidua samples indicated in red.

Source Data

Supplementary information

Supplementary Information

Full legends for Supplementary Tables 1–14, and Supplementary Methods.

Reporting Summary

Supplementary Table 1

Patients and block histology. See Supplementary Information PDF for full legend.

Supplementary Table 2

Artery properties and staging.

Supplementary Table 3

Cell–artery and cell–cell spatial enrichment.

Supplementary Table 4

NanoString DSP analysis of artery ROIs.

Supplementary Table 5

NicheNet analysis.

Supplementary Table 6

Preeclampsia indicated genes.

Supplementary Table 7

Antibodies.

Supplementary Table 8

Number of segmented events per image.

Supplementary Table 9

Positivity binary threshold for functional markers.

Supplementary Table 10

Volcano plots for Figs. 3 and 5.

Supplementary Table 11

NanoString DSP counts and meta-data.

Supplementary Table 12

Cell microenvironment lineage composition.

Supplementary Table 13

NanoString DSP analysis of decidua ROIs.

Supplementary Table 14

EVT by compartment LDA analysis.

Source data

Source Data Fig. 1

Source Data Fig. 2

Source Data Fig. 3

Source Data Fig. 4

Source Data Fig. 5

Source Data Fig. 6

Source Data Extended Data Fig. 2

Source Data Extended Data Fig. 3

Source Data Extended Data Fig. 4

Source Data Extended Data Fig. 5

Source Data Extended Data Fig. 6

Source Data Extended Data Fig. 7

Source Data Extended Data Fig. 8

Source Data Extended Data Fig. 9

Source Data Extended Data Fig. 10

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and Permissions

About this article

Cite this article

Greenbaum, S., Averbukh, I., Soon, E. et al. A spatially resolved timeline of the human maternal–fetal interface. Nature 619, 595–605 (2023). https://ift.tt/70Wugdz

Download citation

Received: 03 September 2021
Accepted: 08 June 2023
Published: 19 July 2023
Issue Date: 20 July 2023
DOI: https://ift.tt/70Wugdz

Comments

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.

Adblock test (Why?)

"interface" - Google News
July 19, 2023 at 10:18PM
https://ift.tt/6mvGsjA

A spatially resolved timeline of the human maternal–fetal interface - Nature.com
"interface" - Google News
https://ift.tt/TGBKI7x
https://ift.tt/XdVS6vJ

Bagikan Berita Ini

A spatially resolved timeline of the human maternal–fetal interface - Nature.com

Abstract

Main

Multiplexed imaging of the decidua

SAR correlates with local composition

Immune composition correlates with GA

Upregulation of tolerogenic markers with GA

Spatiotemporal EVT phenotype

Discussion

Methods

Retrospective cohort design

TMA construction

Antibody preparation

Tissue staining

MIBI-TOF imaging

Low-level image processing

Feature annotation

Single-cell segmentation

Segmentation post-processing

Single-cell phenotyping and composition

Definition of thresholds for functional marker positivity

Two-colour IHC

Blinded manual artery staging

Automated digitization of artery morphological features

Calculation of continuous SAR remodelling score δ

Cell-type frequency as a function of GA and SAR

NanoString GeoMx DSP

Slide preparation

Sample collection

GeoMx DSP NGS library preparation and sequencing

Normalization and scaling of GeoMx counts data

Gene expression in artery as a function of GA and SAR

Coordinated gene expression by pathways in the artery

Ridge regression for predicting GA from immune composition

Cell–cell and cell–artery spatial enrichment analysis

Cell–cell and cell–artery enrichment temporal trends and trending with GA or SAR or constant

Functional marker positivity rate per cell type as a function of GA and SAR

Cellular microenvironments

Definition of anatomical EVT location and associated arteries

SMA and endothelium-loss scores

Characterization of EVTs by compartment

Origin of CD56+ EVTs in the intravascular compartment

DEGs in EVTs

NicheNet analysis

Statistical analyses

Ethics statement

Reporting summary

Data availability

Code availability

References

Acknowledgements

Author information

Authors and Affiliations

Contributions

Corresponding authors

Ethics declarations

Competing interests

Peer review

Peer review information

Additional information

Extended data figures and tables

Supplementary information

Source data

Rights and permissions

About this article

Cite this article

Share this article

Comments

0 Response to "A spatially resolved timeline of the human maternal–fetal interface - Nature.com"

Post a Comment

Origin of CD56⁺ EVTs in the intravascular compartment