Two healthy pregnant women for termination of pregnancy (TOP) surgery at the eighth week due to psychosocial reasons were recruited in this study. The TOP surgery was performed at the Department of Obstetrics and Gynaecology, The University of Hong Kong-Shenzhen Hospital. Before the surgery, the participants received genetic counselling, and signed a consent form for the study, which was approved by the Institutional Review Board of The University of Hong Kong-Shenzhen Hospital ([2021]208) and the Institutional Review Board of BGI (BGI-IRB 22118).

The placentas were collected by vacuum aspiration. The villi samples were sorted immediately after the surgery based on the villous morphology. Fetal membranes and tissues with blood clots were discarded. About 1 cm × 1 cm × 1 cm of the fresh villi samples were washed in a 1.5 ml tube containing 1 × pre-cooled phosphate buffered saline (PBS) (GIBCO, 70011044, from 10 × PBS diluted with diethylpyrocarbonate treated water) on ice to remove the residual blood. The samples were then mixed with pre-cooled Tissue-Tek optimum cutting temperature (OCT) compound (Sakura, 4583), snapped frozen in liquid nitrogen and stocked at −80°C. All these procedures were completed within 30 min after sample collection.

For paraformaldehyde (PFA)-fixed samples, the villi were trimmed to less than 1 cm × 1 cm × 2 cm block, washed with PBS and immersed in 4% PFA (BOSTER, AR1069) at a 10:1 ratio of volume for 16 h at 4°C (Figure 1A). The fixed tissues were then dehydrated in 10% sucrose for 4–6 h, 20% sucrose overnight and 30% sucrose for 24 h at 4°C. The dehydrated samples were embedded in OCT as mentioned above. The fresh villi and the PFA-fixed samples were transferred in dry-ice. Cryosections were conducted sagitally at a thickness of 10 μm in a cryostat (DAKEWE Cryostat Microtome, 6250). Three successive cryosections of the largest intersecting surface were obtained from each sample. One of the cryosection was used for Hematoxylin and Eosin (H&E) staining and the other two sections were used for Stereo-seq as replicates.

For the OCT-embedded fresh samples, Stereo-seq libraries were prepared and sequenced using the STOmics Gene Expression kit S1 (BGI, 1000028493) following the standard protocol V1.1 as previously described (Chen et al., 2022a). Briefly, tissue sections were immediately adhered to a Stereo-seq chip at −20°C, and warmed by fingertip to melt the OCT. Then the chip was fixed in methanol for 30 min at −20°C. After fixation, single-strand DNA dye (Thermo fisher, Q10212) was added to image the nucleus under a Ti-7 Nikon Eclipse microscope. Then the permeabilization enzymes were added to release the transcripts under a condition with pH = 2. After reverse transcription, the residual tissues were digested, and the cDNA was released from the chips and purified using the AMPure XP beads (Vazyme, N411-03). The indexed Stereo-seq libraries were quantified by the Qubit ssDNA Assay Kit (Thermo Fisher Scientific, Q10212). The PFA-fixed samples were prepared with an extra de-crosslinking procedure by immersing the chips into the Tris(hydroxymethyl)aminomethane-ethylendiaminetetraacetic acid (Tris-EDTA) buffer (pH = 10) containing 5% RNase inhibitor at 70°C for 1 h before permeabilization. After de-crosslinking and permeabilization, the chips were re-fixed in pre-cooled methanol at −20°C for 15 min. The captured RNA were reversely transcribed to cDNA. After amplification, fragmentation and cyclization, the DNA nanoballs (DNB) were sequenced on a MGI DNBSEQ-T10 sequencer (MGI, Shenzhen, China) according to the manufacturer’s protocol (50 bp for read 1, 100 bp for read 2) (Chen et al., 2022a).

The fastq files containing sequenced reads generated by the MGI DNBSEQ-T10 sequencer were processed using the STOmics Analysis Workflow (SAW) v4.1.0 (https://github.com/BGIResearch/SAW) (Chen et al., 2022a) to generate spatial gene expression matrices (Figure 1B). The coordinate identity (CID) and the unique molecular identifiers (UMIs) sequences are in read 1 while the cDNA sequences are in read 2. Firstly, CID were mapped to the designed coordinates of the in situ captured chip, and low-quality sequences were filtered out based on the UMIs quality score. Then read 2 was aligned to the reference genome (hg38) using STAR (Dobin et al., 2013). Aligned reads with MAPQ < 10 were filtered out, and retained reads were further mapped to corresponding genes. After merging barcode reads count, a CID-containing expression profile matrix was generated for further quality control. The Stereo-seq raw data were analyzed at the resolution of bin50 (50 × 50 DNB, or 37.5 μm × 37.5 μm) in this study, with the transcripts from the same genes combined into an aggregated profile within each bin as described earlier (Chen et al., 2022a). First, the distributions of both log-transformed total counts and log-transformed n genes by counts were calculated, and the presence of bimodal distributions suggested strong background signals. After that, bin50s with the percent counts in top 100 genes reaching 100% were removed. The remaining bins were then clustered by using the Spatially Constrained Clustering (SCC). The clusters generated from background signals were filtered out. The Wilcoxon test was conducted to identify the differentially expressed genes (DEG). The dissociative counts from dead cells were filtered, since most of the top genes with high expression levels were mitochondrial and ribosomal genes. After filtering out the low-quality clusters, the final expression matrix was used for the downstream analysis.

Principal component analysis (PCA) of the expression matrix were conducted using the Scanpy package (Wolf et al., 2018). SCC was performed on the dimensional-reduced data with consideration of both gene expression similarities and neighborhood relations in spatial distribution. In brief, the spatial k-nearest neighbor (KNN) graph, produced by the Squidpy (Holgersen et al., 2021), was integrated with the KNN graph which was built from the transcriptomic data by the Scanpy. The graph including gene expression and spatial information was then passed to the Leiden algorithm for further clustering, and the markers of each output clusters were identified by the rank_genes_groups function of Scanpy with the Wilcoxon method. All clusters were annotated by combining the information from the marker genes and canonical markers of villi cell types.

The gene regulatory network (GRN) was constructed following the standard procedures of the pySCENIC pipeline (Aibar et al., 2017). Co-expressing genes from the bin50 expression matrix with the human HUGO Gene Nomenclature Committee (hgnc) (Bruford et al., 2020) transcription factors (TFs) were identified using the GENIE3 algorithm (Huynh-Thu et al., 2010). To exclude false-positive results from the co-expression networks, the cisTarget database was then adopted to confirm the putative direct-binding targets of each TF, and filtered out the indirect targets lacking motif support. The verified co-expression modules were regulons with significant motif enrichment of the corresponding regulator. Finally, the AUCell algorithm was conducted to score the regulon activities of all bins and identify cell subgroups with high subnetwork activities (Aibar et al., 2017).

To identify cell type abundances in each spatial spot, we conducted deconvolution analysis by the cell2location software (Kleshchevnikov et al., 2022). A high-quality scRNA-seq reference panel is necessary for the accurate estimation of cell types in the spatial transcriptomic data, and here the scRNA-seq data from a previous study (Vento-Tormo et al., 2018) was adopted for the spatial data deconvolution. We selected the data sampled from placenta tissues and filtered out cells not belonging to SCT, VCT, EVT, FB, HB, and Endo for further analysis. First, the cell2location identified reference cell type signatures from the scRNA-seq profile using the negative binomial regression. Then, the cell-type gene-signature model was used to decompose the spatially resolved gene-count matrices into the reference cell-type signatures, thereby estimating the absolute cell abundance of individual cell types across all spatial spots.

To study crosstalk between cell clusters mediated by ligand-receptor complexes, cell-cell communication analysis was performed using the CellPhoneDB software (https://github.com/Teichlab/cellphonedb). The CellPhoneDB datasets of ligands, receptors and their interactions were adopted, and the statistically significant enriched ligand-receptor interactions were identified based on the expression levels of ligands and coordinate receptors (Efremova et al., 2020).

The Stereo-seq data of fresh and PFA-fixed treated villi were dramatically different in terms of tissue morphology and data quality. The cryosections of the fresh villi embedded in OCT showed poor morphology and changed tissue shapes, making it difficult to determine the structure. The H&E staining results of the sections showed that the outer SCT layer were ruptured. Many cells, especially the inner core cells beneath the SCT were lost (Supplementary Figure S1A). Approximate 1.8 G reads were generated in each section of the fresh sample with a Q30 bases of 89.5%. The number of gene types (mean gene type) detected in each bin50 area (37.5 μm × 37.5 μm) was 399 and the number of transcripts (mean molecular identifier, MID) detected in each bin50 was 820 (Supplementary Table S1). We observed strong diffusion of the transcriptome signals of the fresh villi from the sequencing results (Figure 1C), possibly due to RNAs spillover during the permeabilization of villi. In contrast, H&E staining of the PFA-fixed villi showed intact villi structures and good cell morphology, especially the integrity of the cells in the inner core of villi (Supplementary Figure S1B). Approximate 5.5 G reads were generated from the PFA fixed villi with a Q30 bases of 90.9%. The mean gene type of bin50 was 617 and the MID of bin50 was 1.59 K (Supplementary Table S1). Compared to the fresh villi, the sequencing results of the cryosections of the PFA-fixed villi showed dramatic reduction of the RNA signal diffusion (Figure 1D). Due to the severe RNA diffusion, trophoblast-specific markers could not be identified in the fresh villi (Supplementary Figure S2A). In contrast, Specific RNA signals of trophoblast-specific markers were identified in the PFA-fixed villi (Supplementary Figure S2B).

The PFA fixed villi were analyzed in detail for functional regions since they had better data quality and tissue morphology. A series of quality control (QC) steps were performed to remove background noise and low-quality bins that were mainly composed of mitochondrial and ribosomal RNAs. After data QC, the cell clusters corresponded well to the ssDNA image (Supplementary Figure S3). After cell bin clustering by SCC and cell type annotation, four major cell clusters were identified using the classic markers of placenta cells. To demonstrate the fine structure of a villus, a chip region with a representative integral villi structure was selected for detailed spatial transcriptome analysis (Figures 1D, 2A). Markers of SCT namely, CGA (Vento-Tormo et al., 2018), PSG1 and CSH1 (Sheridan et al., 2021) were identified in the outermost layer markers; VCT were identified in the inner layers of SCT using EGFR (Fock et al., 2015), ITGA6, and LRP5 (Shannon et al., 2022) as markers; markers of EVT (MMP2, FN1, and HLA-C) (Papuchova et al., 2019) were identified in an anchoring villus-like cell column structure (Figure 2B). Besides trophoblasts, markers of FB including COL1A1, COL1A2 and PITX2 (Vento-Tormo et al., 2018) were identified in the core of the villi (Figure 2B). We also identified the markers of HB cells (MTF, MAF and POU2F2) and Endo cells (ZEB1, ZBTB20 and PECAM1) (Papuchova et al., 2019) in the core of villi (Figure 2C). Since the FB, HB and Endo were mixed in the core of villi, we denoted them as the inner core cells. All the cell clusters were then deconvoluted by the Cell2location using the RNA-seq data from a previous study (Vento-Tormo et al., 2018). The results suggested that the HB and the Endo cells were dispersed among the FBs (Figure 2D). A localization map of different cell types in villi is shown in Figure 2E. Applying the same methodology, similar cell types and spatial localizations were found in the PFA-fixed villi of another chip (Supplementary Figure S4).

Gene ontology (GO) annotation was used to find the DEG and the involved biological processes in different trophoblast (Supplementary Figure S5; Supplementary Table S2). The most significantly upregulated GO biological processes in EVT versus VCT included cell migration, cell adhesion and cytoskeleton, while the pathways of cell proliferation such as protein translation, folding and negative regulation of degradation were downregulated in EVT. This result was in accordance with the previous data (Apps et al., 2011). Compared to VCT, the most significantly upregulated GO pathways in SCT included regulation of protein transport, hormone secretion, negative regulation of immune system process, viral entry into host cell, response to ER stress and apical plasma membrane, whereas the positive regulation of growth pathways was downregulated.

To investigate cell communications, we analyzed the expression level of ligand-receptor pairs in different cell types of the villi. SCT and VCT exhibited the highest number of ligand-receptor pairs, while the inner core cells have the least number of ligand-receptor pairs (Figure 3A; Supplementary Table S3). As expected, the ligand DLK localized mainly to FB, while its receptors, NOTCH2 and NOTCH3, were expressed in all trophoblasts, implying a role of the interaction in maintenance of FB through paracrine signaling from trophoblast (Figure 3B). The interaction between PGF and FLT1 was abundant in villi, both expressed in nearly all trophoblasts. The imbalance of this interaction is an indicator for PE (Maynard et al., 2003). The imprinted Igf2-Igf2r pairs have been well studied. In mice, Igf2 is expressed highly in feto-placental endothelial cells, and is important in maintaining trophoblast development, particularly the SynT-II layer (Sandovici et al., 2022). Our study found that IGF2 was abundant in FB and VCT, while IGF2R was in all cell types. Wnt ligands such as Wnt5A and WNT7A were highly expressed in SCT, while the receptors of Wnt, including FRZB and LDLR, were expressed in SCT and inner core cells (Figure 3B). Wnt induces SynT-II cells fusion in mice through induction of Gcm1 and SynB expression (Zhu et al., 2017). Here, Wnt may function on SCT through an autocrine action. EGFR was also highly expressed in SCT and VCT, while the ligands were expressed in nearly all the villous cells (Figure 3B). The expression of EGFR increases with differentiation of the SCT (Rab et al., 2013). With the Stereo-seq data, we directly mapped the location of various ligand receptors in the villi, and showed their co-localization with corresponding ligands inferring the physical contact between specific cells with ligand-receptor pairs (Figure 3C).

Using the SCENIC analysis, we identified a total of 316 regulons in EVT, VCT, SCT and inner core cells. The top 20 regulons specific in each cell type were identified based on the regulon specificity score (RSS). The heatmap generated from the Z-scores of these cell-type-specific regulons shows that the cell clusters based on the regulons AUC scores (Figure 4A; Supplementary Table S4). A group of TFs reported to be specifically expressed in the placenta was used to confirm our analysis of spatial regulatory activity (Figure 4B). For example, TFAP2C, MSX2, FOSL1 and HMGA1 are expressed in VCT (Renaud et al., 2014; Hornbachner et al., 2021). We showed consistent findings on the expression and spatial localization of these TFs (Figures 4B,C). However, we did not find the expression of FOSL1 and HMGA1 in the EVT, which is against the report that FOSL1 is also expressed in the proliferative EVT cell columns, while HMGA1 in the EVT invading the decidua (Renaud et al., 2014). TCF7L2, TEAD1, ASCL2, and SMAD3 have been reported to be expressed specifically in EVT (Meinhardt et al., 2014; Dong et al., 2022; Haider et al., 2022). These TFs were also highly expressed in EVT in our results (Figure 4B). We also confirmed the dominant expressions and spatial distribution of AHR, CEBPA, TBX3 and NFE2L3 in the nuclei of SCT (Figures 4B,C), consistent with previous studies (Simmons et al., 2008; Iqbal et al., 2021; Chen et al., 2022b; Dong et al., 2022). GATA6, ZEB1, RUNX1, and ERG are fibroblast or endothelial cells-related TFs (Molkentin, 2000; Ponder et al., 2016; Yuan et al., 2020; Zhang et al., 2021b; Caporarello et al., 2022). We confirmed that these TFs were mostly expressed in the inner core cells (Figures 4B,C). To further validate our results, we acquired the protein expression data of the above TFs from the Human Protein Atlas (https://www.proteinatlas.org/), which showed the concordant expression in cell types and localization (Supplementary Figure S6).

By GO analysis, the regulons upregulated in SCT relative to VCT were related to hormones transport and secretion, and response to hormones (Figure 4D; Supplementary Table S5). This is consistent with the fact that the SCT synthesize and secrete a variety of hormones, including hCGA and hCGB, CSH1, PSGs, PGF, CRH and others (Jaremek et al., 2021). Meanwhile, the regulons associating with the regulation of mitotic cell cycle were downregulated, which is concordant with the exit from the mitosis before differentiation of SCT. Pathways of nutrient transport and mammary gland development were also identified (Figure 4D). In EVT, pathways such as the epithelial to mesenchyme transition pathway, epithelial development and differentiation, extracellular matrix organization were identified as previously described (Illsley et al., 2020; Papuchova and Latos, 2022), which were associated with the invasive activities of EVT (Figure 4D; Supplementary Table S5).

As the outermost cell type, SCT has direct contact with the maternal interface. We found that the DEG of SCT contained the genes involving in viral entrance to cells and viral life cycles, implying possible vertical transmission of virus from mother to fetus through SCT (Supplementary Figure S7; Supplementary Table S6). To evaluate the placental receptors of virus transmission, we mapped the expression of a panel of 20 viral receptors and co-receptors for TORCH (Toxoplasma gondii, Other agents, Rubella, Cytomegalovirus & Herpes simplex virus), severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), ZIKA, and other infectious diseases in the villi (Figure 5A). Generally, EGFR, SDC1, HSPA5, MYH9, CD46 were strongly expressed in the SCT but also widely expressed in other cell types. In contrast, ITGAV, CD55, TFRC, SCARB1, HSPG2, and BSG were only mainly expressed in SCT.

Human cytomegalovirus (HCMV) may infect SCT to induce monocyte activation through the receptor EGFR and the co-receptor integrin proteins (Wang et al., 2003; Chan and Guilbert, 2005). However, this theory is still controversial, as the virion may across the syncytium by transcytosis instead of direct infection (Fisher et al., 2000; Maidji et al., 2006). We found that SCT exhibited little expression of the neonatal Fc receptor (FCGRT), which helps to transport IgG for passive immunity. Yet, the key integrin receptors, α1β1 and αVβ3, together with EGFR had moderate or strong expression in SCT (Figures 5A,B), which provides the molecular basis for viral infection. The protein expressions of α1β1 and αVβ3 in SCT were also verified by the data from the Human Protein Atlas (Supplementary Figure S8).

Till now, it is still controversial whether the SARS-CoV-2 can be vertically transmitted (Islam et al., 2020; Rodrigues et al., 2020; Yoon et al., 2020; Kotlyar et al., 2021). We assessed the expression of a series of SARS-CoV-2 receptors in villi. Our results indicated a low expression of ACE2 and TMPRSS2 in all cell types in the villi (Figures 5A,B), consistent with previous single-cell analysis studies (Lu et al., 2020; Pique-Regi et al., 2020). However, BSG and HSPA5 were expressed highly in the SCT as previously reported (Dong et al., 2021; Lee et al., 2022) (Figures 5A,B), inferring the possibility of vertical transmission of SARS-CoV-2.

In SCT, the transporter genes were strongly upregulated, inferring their roles in transportation of nutrition, gas, and hormones between mother and fetus (Supplementary Figure S7; Supplementary Table S6). Specifically, we analyzed the transporter genes relevant to drug transportation, including 52 members in ATP-binding cassette (ABC) family (Sarkadi et al., 2006) and 395 members in solute carrier protein (SLC) families (Schlessinger et al., 2013). A total of 47 transporters were enriched in SCT, most of which have been previously reported (Figure 5C; Supplementary Figure S9). ABCB1, for example, has been reported to localize at the apical membrane of SCT to back-transport xenobiotics and drugs to maternal circulation (Sun et al., 2006). ABCG2 is also localized predominantly to the apical surface of SCT, and limits the transfer of nitrofurantoin (Merino et al., 2005), cimetidine (Pavek et al., 2005) and glyburide (Gedeon et al., 2008) to fetal circulation. Moreover, SLC6A4 (Wang et al., 2007), SLC22A11 (Shimizu et al., 2005), SLC2A1, SLC3A2, SLC4A2 (Nishimura and Naito, 2008), SLC7A5 (Jia et al., 2022) protein or RNA have been shown to localize to SCT. These were all detected in our data (Figures 5C,D; Supplementary Figure S9). We also found transporters that had not been elucidated in the placenta, such as SLC40A1, SLC19A2 and SLC31A2 (Figure 5D; Supplementary Figure S10), which could be candidates for future study on blood-placental barrier.