Embryo implantation into the uterus is a crucial step for human reproduction. A hypothesis has been proposed that the molecular circuit invented by trophoblasts for invasive embryo implantation during evolution might be misused by cancer cells to promote malignancy. Unfortunately, our current understanding of the molecular mechanism underlying embryo implantation is far from complete.
Here we used the mouse as an animal model and generated a single-cell transcriptomic atlas of the embryo implantation site of mouse uterus at the invasion phase of embryo implantation on gestational day 6. We revealed 23 distinct cell clusters, including 5 stromal cell clusters, 2 epithelial cell clusters, 1 smooth muscle cell cluster, 2 pericyte clusters, 4 endothelial cell clusters, and 9 immune cell clusters. Through data analysis, we identified differentially expression changes in all uterine cell types upon embryo implantation. By integrated with single-cell RNA-seq data from E5.5 embryos, we predicted cell–cell crosstalk between trophoblasts and uterine cell types.
Our study provides a valuable resource for understanding of the molecular mechanism of embryo implantation.
Embryo implantation in humans is interstitial. It consists of the following three phases: embryo apposition, attachment, and invasion. Immediately after apposition and attachment to endometrial epithelium, the blastocyst penetrates through the epithelium, followed by the basal lamina, and invades into the stroma . It has been well established that excessive trophoblast invasion may lead to the pathogenesis of placenta creta , while shallow trophoblast invasion may result in pre-eclampsia and intrauterine growth restriction . Interestingly, it has long been recognized that there are striking similarities between trophoblast cells and invasive cancer cells, both of which share proliferative, migratory and invasive properties, as well as the capacity to confer immune privilege . A hypothesis has been proposed that the molecular circuit invented by trophoblasts for invasive embryo implantation during evolution might be misused by cancer cells to promote malignancy . Unfortunately, our current understanding of the mechanisms underlying embryo implantation is far from complete.
Due to ethical restrictions and experimental difficulties, in vivo analysis of embryo implantation heavily relies on mice . Slightly different from humans, embryo implantation in mice is eccentric. The apposition phase of mouse embryo implantation is thought to occur from the morning to the midnight of gestational day 4 (GD4). During this phase, the embryo seeks its position on the luminal epithelium of uterus. After that, the embryo attaches to the receptive uterine epithelium. A firm connection between the blastocyst and uterine luminal epithelium is established on the morning of GD5. The invasion phase of embryo implantation occurs on GD6, when the blastocyst penetrates through the epithelium and invades into the stroma [7, 8]. Previously, two studies have analyzed global gene expression changes associated with the invasion phase of embryo implantation in mice by using high-throughput transcriptomic approaches [9, 10]. The limitation of both studies was that the whole uterus was used. The uterus is a complex structure consisting of many cell types, including luminal and glandular epithelial cells, stromal cells, smooth muscle cells, endothelial cells, and various immune cells . Thus, whole uterus transcriptomic studies were unable to accurately capture cell-type-specific gene expression changes.
In the present study, by using the-state-of-the-art single-cell RNA-seq approach , we resolved all cell types at the embryo implantation site of the mouse uterus on GD6. Through data analysis, we identified differentially expression changes in all cell types at the invasion phase of embryo implantation. We also predicted cell–cell crosstalk between trophoblasts and uterine cell types. Our study provides a valuable resource for understanding of the molecular mechanism of embryo implantation.
Materials and methods
Adult CD-1 mice of the SPF grade were used in this study. All mice were caged under light-controlled conditions (14 h/10 h light/dark cycles) with free access to regular food and water. Female mice were mated with fertile males and success of mating was confirmed the next morning by the presence of a vaginal plug. The day of the vaginal plug was denoted as gestation day 1 (GD1). On GD6, the implantation sites and inter-implantation sites (served as a control) were collected separately. All animal procedures were approved by the Institutional Animal Care and Use Committee of South China Agricultural University (No. 2020B078, approved on 29/09/2020).
Bulk RNA-seq analysis
The total RNAs from uterine tissues were extracted with the TRIzol reagent (Invitrogen). RNA-seq libraries were generated by using the TruSeq RNA sample preparation kit (Illumina) and sequenced on an Illumina HiSeq 2500 system. Using UCSC mm10 mouse genome as reference, raw data were analyzed using TopHat v2.0.4  and Cufflinks v2.2.1  as described previously . Gene expression levels were measured as fragments per kilobase per million (FPKM).
Single-cell dissociation of mouse uterus
Single-cell suspension was prepared as described previously [16, 17]. Briefly, the uterine tissues from 3 mice for each group were pooled and minced with a blade. Tissues were then incubated in dissociation buffer containing 2 mg/ml Collagenase II (#C6885, Sigma-Aldrich), 10 mg/ml Dispase II (#354,235, Corning) and 50,000 U/ml DNase I (#DN25, Sigma-Aldrich) for up to 30 min at 37 °C in a shaking incubator. The digestion progress was checked every 5 min with a microscope until a single cell suspension was achieved. The single-cell suspension was then passed through a 40-μm cell strainer to remove undigested tissues. Cells were spun down at 250 g at 4 °C for 4 min and the pelleted cells were washed using centrifugation. In order to measure cell viability, cells were strained with AO/PI solution (#CS2-0106, Nexcelom Bioscience) and counted using a Cellometer Auto 2000 instrument (#SD-100, Nexcelom Bioscience). The single-cell suspension was carried forward to single-cell RNA-seq only if the cell viability was > 80% and the percentage of cell clumps was < 10%.
Single-cell RNA-seq library preparation and sequencing
The final concentration of single-cell suspension was adjusted to 1000 cells/μl and a volume of 15 µl was loaded into one channel of the ChromiumTM Single Cell B Chip (#1,000,073, 10 × Genomics), aiming at recovering 8000–10,000 cells. The Chromium Single Cell 3' Library & Gel Bead Kit v3 (#1,000,075, 10 × Genomics) was used for single-cell bar-coding, cDNA synthesis and library preparation, following the manufacturer's instructions provided as the Single Cell 3' Reagent Kits User Guide Version 3. Library sequencing was performed on an Illumina novaseq 6000 system configured with the paired-end 150-bp protocol at a sequencing depth of approximately 400 million reads.
Single cell RNA-seq data analysis
Raw data bcl files from the Illumina NovaSeq 6000 platform were converted to fastq files using the bcl2fastq tool (v22.214.171.1246). These fastq files were aligned to the mm10 mouse reference genome by using the CellRanger software (v3.0.1, 10 × Genomics). The resulting gene counts matrix was analyzed with the R package Seurat (v3.1.3) . Cell with fewer than 200 or greater than 6000 unique genes, as well as cells with greater than 25% of mitochondrial counts, were excluded. Meanwhile, genes expressed in fewer than 3 cells were removed. Following data filtering, the gene counts matrix was normalized and scaled by using NormalizeData and ScaleData, respectively. The top 2000 highest variable genes were used for the principal component analysis (PCA) and the optimal number of PCA components was determined by the JackStraw procedure. Single cells were clustered by the K-nearest neighbor (KNN) graph algorithm in PCA space and visualized using the t-distributed stochastic neighbor embedding (tSNE) dimensional reduction technique. The cell type label for each cell cluster was manually assigned based on canonical cell markers. The FindAllMarkers function was used to identify novel marker genes for each cluster with a minimum of 10% of cells expressing the gene within the cluster and a minimum logFC threshold of 0.25. In order to find differential expressed genes in the same cell type between pre-receptive uterus and receptive uterus, the FindMarkers function in Seurat was used with min.logfc being set to 0.25 and min.pct being set to 0.20.
Monocle2 package v2.18.0  was used for pseudotime analysis. The count data and meta data were export from the Seurat object and then import into the CellDataSet object in Monocle2. Feature genes were selected by using the differentialGeneTest function. After The dimension reduction by using the DDRTree algorithm, the orderCells function was used to infer the trajectory with default parameters. The reconstructed trajectory was visualized by the plot_cell_trajectory function.
Gene ontology analysis
Gene ontology (GO) analysis was performed as described previously . GO terms were grouped according to the biological process category in the Mouse Genome Informatics (MGI) GOslim database . To test for enrichment, a hypergeometric test was conducted and P < 0.05 was used as significance threshold to identify enriched GO terms.
Pathway enrichment analysis was conducted by using the Metascape v7.4 online tools . The significance threshold for FDR was set at 0.05.
Cell–cell communication analysis
The CellChat v1.1.0 software  was used to infer cell–cell communication based on ligand-receptor interaction with default parameters. For each ligand-receptor pair, CellChat assigned a communication probability value by the law of mass action based on the average expression values of a ligand by one cell group and that of a receptor by another cell group. The statistical significance of communication probability values was assessed by a permutation test. P < 0.05 were considered statistically significant.
A single-cell atlas of mouse uterus on gestational day 6
To create a cell-type resolved map of mouse uterus at the invasion phase of embryo implantation, we performed single-cell RNA-seq analysis (Fig. 1A). The implantation site (IS) and the inter-implantation site (IIS) of mouse uterus were collected on gestational day 6 (GD6). The whole uterus, which is consist of endometrium, myometrium and perimetrium, was used. The embryo at IS was also kept. Single-cell RNA-seq data were generated by using the 10× Genomics platform. After quality control, a total of 16,257 cells (7065 for IIS and 9192 for IS) were obtained (Fig. 1B). In order to validate this single-cell RNA-seq dataset, we also generate a bulk RNA-seq dataset using the same samples. It turned out that the cell-averaged single-cell RNA-seq data were highly accordant with the conventional bulk RNA-seq data (r = 0.7523 for IIS and r = 0.7759 for IS), indicative of high quality of our single-cell RNA-seq data (Fig. 1C).
Unsupervised clustering analysis revealed 23 distinct cell clusters for all cells from IIS and IS combined (Fig. 2A). Major cell types were defined using the expression of known cell type-specific genes, with hormone-responsive cells expressing Pgr and Esr1 [24, 25], endothelial cells expressing Pecam1  and immune cells expressing Ptprc  (Fig. 2B).
Hormone-responsive cells included epithelial cells expressing Epcam and Krt19  (Fig. 2C), stromal cells expressing Hoxa11  (Fig. 2D), smooth muscle cells expressing Acta2  and pericytes expressing Rgs5  (Fig. 2E). We found 2 epithelial cell clusters, LE and GE. LE was luminal epithelial cells expressing Tacstd2 and GE was glandular epithelial cells expressing Foxa2 . We identified 5 stromal cell clusters, S1, S1p, S2, S2p and S3. Cells in S2 but not S1 expressed high levels of Hand2, implying that S2 was superficial stromal cells and S1 was deep stromal cells . S3 was primary decidual zone stromal cells expressing decidualization marker genes Wnt4  (Fig. 2F). S1p and S2p were a subset of proliferating S1 and S2 with high level of Mki67, respectively (Fig. 2G). Only 1 smooth muscle cell cluster, SMC, was found. Meanwhile, 2 pericyte clusters, PC and its proliferating subset PCp, were identified. Endothelial cells had 4 clusters: VEC and its proliferating subset VECp are vascular endothelial cells expressing Sox17, while LEC and its proliferating subset LECp are lymphatic endothelial cells expressing Prox1  (Fig. 2H). There were 9 immune cell clusters (Fig. 2I, J). Included were natural killer cells (NK, Ptprc+Nkg7+Cd3e−) , proliferating natural killer cells (NKp, Ptprc+Nkg7−Cd3e+Mki67+), T cells (T, Ptprc+Nkg7−Cd3e+) , B cells (B, Ptprc+Cd79a+Ms4a1+) , macrophages (M, Ptprc+Adgre1+) , proliferating macrophages (M, Ptprc+Adgre1+Mki67+), dendritic cells (DC, Ptprc+Itgax+) , proliferating dendritic cells (DC, Ptprc+Itgax+Mki67+) and plasmacytoid dendritic cells (pDC, Ptprc+Siglech+) .
Finally, we aimed to discover novel markers for each cell type. We selected genes that were expressed significantly higher in the cell type of interest than the other cell types by Wilcoxon rank sum test. A complete list of these marker genes was presented in Additional file 1: Table S1.
Reconstruction of developmental trajectory for primary decidual zone
In our single-cell RNA-seq data, we identified 5 clusters of stromal cells: S1 (deep stromal cells), S1p (proliferating deep stromal cells), S2 (superficial stromal cells), S2p (proliferating superficial stromal cells) and S3 (primary decidual zone stromal cells, PDZ). We selected signature genes for each cell cluster by using Wilcoxon rank sum test. After the removal of redundancy, we identified a total of 1784 signature genes (Additional file 2: Table S2). Through heatmap, we grouped all these signature genes into 4 gene sets (Fig. 3A). Gene set 1 with 403 genes were S1-specific. Gene ontology analysis showed that these genes were enriched in cell adhesion (P = 4.20 × 10–11), developmental processes (P = 1.10 × 10–4), cell organization and biogenesis (P = 3.27 × 10–4), stress response (P = 4.88 × 10–4) and protein metabolism (P = 3.87 × 10–2). Gene set 2 with 352 genes were decreased in S3 compared to its intermediate S2p. These genes were enriched in DNA metabolism (P = 1.00 × 10–11), cell cycle and proliferation (P = 1.00 × 10–11), cell organization and biogenesis (P = 2.59 × 10–11) and protein metabolism (P = 3.97 × 10–3). Gene set 3 of 300 genes were S2-specific. Based on GO, enriched terms were protein metabolism (P = 3.39 × 10–11), developmental processes (P = 6.05 × 10–8) and cell cycle & proliferation (P = 8.09 × 10–3). Gene set 4 of 729 genes were unchanged or increased in S3 compared to its intermediate S2p. Enriched GO terms were protein metabolism (P = 2.49 × 10–8), RNA metabolism (P = 3.23 × 10–7), DNA metabolism (P = 4.42 × 10–5), transport (P = 1.96 × 10–3) cell organization and biogenesis (P = 2.90 × 10–5), and cell cycle and proliferation (P = 4.91 × 10–3).
The expression of known marker genes for PDZ (S3) were examined. We found that PDZ expressed pan-stromal cell markers Pgr, Esr1 and Hoxa11, as well as the superficial stromal cell marker Hand2. Additionally, PDZ expressed decidualization marker Wnt4, but not Prl8a2 . PDZ ceased proliferation from its intermediate S2p, showing low expression of Mki67. Interestingly, although PDZ expressed mesenchymal marker Vim, they were also positive for epithelial marker Krt19 (Fig. 3B).
To further reveal the relationship between these 5 stromal cell clusters, pseudotime trajectory analysis was conducted. Cells were arranged in a pseudotime manner with a pedigree reconstruction algorithm for biological processes based on transcriptional similarity. We found 2 paths of interest: (1) primary decidual zone formation, i.e. S2- > S2p/S3; and (2) secondary decidual zone formation, i.e. S1- > S1p (Fig. 3C).
Cell–cell communication between primary decidual zone and trophoblast giant cells
The cell–cell communication between primary decidual zone and trophoblast giant cells represents the key mechanism of embryo implantation. However, due to the relatively small number of embryonic cells at the implantation site, we did not find any embryo-related cell clusters in our single-cell RNA-seq data. Alternatively, we re-analyzed a published single-cell RNA-seq dataset on mouse E5.5 blastocysts . E5.5 is equivalent to gestational day 6 in our study. By using canonical marker genes, 4 major cell types were identified (Fig. 4A). The 4 cell clusters were visceral endoderm (VE) expressing Apob and Amn , epiblast (EPI) expressing Pou5f1 and Nanog , extraembryonic ectoderm or ectoplacental cone (EXE/EPC) expressing Elf5 and Cdx2 , and trophoblast giant cells (TGC) expressing Gata2  and Prl3d1  (Fig. 4B–E).
By using the CellChat software, we predicted the ligand-receptor interactions between PDZ (S3) and TGC. We found a total of 20 ligand-receptor interaction pairs (Fig. 4F). Based on pathway analysis, these ligand-receptor interactions were enriched among ECM-receptor interaction (FDR = 1.00 × 10–18), PI3K-Akt signaling pathway (FDR = 1.00 × 10–18), Regulation of actin cytoskeleton (FDR = 1.26 × 10–9), Rap1 signaling pathway (FDR = 1.58 × 10–6), Ras signaling pathway (FDR = 7.94 × 10–5), MAPK signaling pathway (FDR = 1.26 × 10–4), Phospholipase D signaling pathway (FDR = 7.94 × 10–3) and Cytokine-cytokine receptor interaction (FDR = 3.98 × 10–2) (Fig. 4G).
The impact of PDZ stromal cells on other stromal cells
We investigated the abundance of major stromal cell types at IS compared to IIS. The χ2 test was employed to assess the significance of difference between two groups. By using the criteria of P < 0.05 and fold change > 2, we found that the proportion of S1 was unchanged, whereas the proportion of S2 significantly decreased in IS compared to IIS. Meanwhile, S1p and S2p were almost exclusively detected in IS (Fig. 5A), and so was S3 (0.0% in IIS vs 3.1% in IS, data not shown). According to our trajectory analysis, S2p was originated from S2 and S3 was originated from S2p. Surprisingly, the proportion of all these S2-lineage cells (S2 + S2p + S3) together in IS was still less than that of S2 in IIS (24.6% vs 36.8%). This phenomenon was likely due to cell loss at PDZ from GD6. In fact, PDZ disappears by GD8 .
We investigated the breadth of transcriptional changes in each stromal cell type by performing differential gene expression analysis. Using a logFC cutoff of 0.25 and a pvalue cutoff of 0.05, we identified 328 and 581 differentially expressed genes for S1 and S2, respectively (Fig. 5B and Additional file 3: Table S3). We then explored the biological implications of differentially expressed genes using gene ontology (GO) analysis. Our results indicated that similar functional changes occurred during embryo implantation in both S1 and S2 (Fig. 5C).
In order to determine the impact of PDZ (S3) on S1 and S2, we used the CellChat software to predict the ligand-receptor interactions. Only secreted factors from PDZ were considered. We found a total of 35 ligand-receptor interaction pairs (Fig. 6A). Pathway analysis revealed that these ligand-receptor interactions were enriched among Cytokine-cytokine receptor interaction (FDR = 1.26 × 10–10), PI3K-Akt signaling pathway (FDR = 2.00 × 10–9), Regulation of actin cytoskeleton (FDR = 3.16 × 10–9), Hippo signaling pathway (FDR = 3.98 × 10–9), ECM-receptor interaction (FDR = 3.98 × 10–8), Rap1 signaling pathway (FDR = 1.00 × 10–6), MAPK (FDR = 3.16 × 10–6), Wnt signaling pathway (FDR = 2.00 × 10–5), Ras signaling pathway (FDR = 2.00 × 10–5), mTOR signaling pathway (FDR = 2.51 × 10–5), TGF-beta signaling pathway (FDR = 2.51 × 10–5), Endocytosis (FDR = 6.31 × 10–5) and Phagosome (FDR = 7.94 × 10–3) (Fig. 6B).
The impact of primary decidual zone stromal cells on immune cells
By using the criteria of P < 0.05 and fold change > 2, we found that the proportion of NKp, Mp and DCp were significantly increased at IS compared IIS (Fig. 7A). Using a logFC cutoff of 0.25 and a pvalue cutoff of 0.05, we identified 286, 294, 302, 495, 263 and 380 differentially expressed genes for T, B, NK, M, DC and pDC respectively (Fig. 7B and Additional file 3: Table S3). Differentially expressed genes were further characterized by GO analysis (Fig. 7C). Our data indicated that each immune cell type invoked distinct biological processes in order to accommodate embryo implantation.
We predicted the ligand-receptor interactions between PDZ (S3) and immune cells (T, B, NK, M, DC and pDC) using the CellChat software. We found a total of 26 ligand-receptor interaction pairs (Fig. 8A). Pathway analysis revealed that these ligand-receptor interactions were enriched among Cytokine-cytokine receptor interaction (FDR = 1.58 × 10–7), PI3K-Akt signaling pathway (FDR = 1.26 × 10–6), MAPK signaling pathway (FDR = 3.16 × 10–4), Endocytosis (FDR = 3.98 × 10–4), ECM-receptor interaction (FDR = 3.98 × 10–4), TGF-beta signaling pathway (FDR = 3.98 × 10–4), Rap1 signaling pathway (FDR = 7.94 × 10–4), Hippo signaling pathway (FDR = 2.51 × 10–3), Regulation of actin cytoskeleton (FDR = 7.94 × 10–3), Toll-like receptor signaling pathway (FDR = 1.00 × 10–2), Ras signaling pathway (FDR = 1.26 × 10–2), Natural killer cell mediated cytotoxicity (FDR = 1.58 × 10–2), Jak-STAT signaling pathway (FDR = 3.16 × 10–2) and NOD-like receptor signaling pathway (FDR = 3.98 × 10–2) (Fig. 8B).
The impact of primary decidual zone on endothelial cells
By using the criteria of P < 0.05 and fold change > 2, we found that the proportions of both VEC and LEC were significantly decreased in IS compared to IIS, which is in line with the fact that PDZ is avascular [42,43,44]. Interestingly, the proportion of VECp was significantly increased in IS compared to IIS, indicating that the PDZ is inducer of angiogenesis in the uterus (Fig. 9A). Using a logFC cutoff of 0.25 and a P value cutoff of 0.05, we identified 263 and 404 differentially expressed genes for VEC and LEC, respectively (Fig. 9B and Additional file 3: Table S3). These differentially expressed genes were further characterized by GO analysis (Fig. 9C).
Embryo implantation is a crucial step for human embryo implantation. Previously, we performed single-cell RNA analysis of the mouse uterus during the apposition phase and attachment phase of embryo implantation on gestational days (GD) 4–5 [16, 17]. Here, by using the mouse as an animal model, we profiled the single-cell transcriptome for 16,257 cells from mouse uterus at the invasion phase on GD6. We revealed 23 distinct cell clusters, including 5 stromal cell clusters, 2 epithelial cell clusters, 1 smooth muscle cell cluster, 2 pericyte clusters, 4 endothelial cell clusters, and 9 immune cell clusters. To the best of our knowledge, the present study is the first to highlight the transcriptome landscape associated with the invasion phase of embryo implantation at single-cell resolution.
In order to accommodate embryo implantation, an important change in mouse uterus is the formation of primary decidual zone (PDZ). The first sign of PDZ formation occurs on the afternoon of GD5. PDZ is fully established by GD6 [45,46,47]. In this study, we found that the stromal cells at IIS can be divided into 2 cell types, superficial stromal cells expressing Hand2 and deep stromal cells which are negative for Hand2. Through pseudotime trajectory analysis, we confirmed that PDZ at IS was derived from superficial stromal cells. Superficial stromal cells start proliferate to become intermediate PDZ cells with high expression of Wnt4, Bmp2 and Mki67. Notably, Prl8a2, the marker for secondary decidual zone , was not expressed. These cells ceased proliferation and became PDZ cells . The expression of another decidual marker Bmp2  was reduced during this process. Interestingly, although PDZ cells expressed mesenchymal marker Vim, they were also positive for epithelial marker Krt19. Our result is in line with previous findings showing that PDZ is avascular and epithelioid in nature [42,43,44]. In fact, PDZ is thought to function as a partial permeability barrier to safeguard the implanting embryo, as only molecules smaller than 45 kDa are freely permeable though PDZ .
The trophectoderm (TE) of mouse embryo is created by the end of the pre-implantation period. After implantation, TE develops into two different types, polar TE and mural TE. Polar TE differentiates into extraembryonic ectoderm (ExE) and ectoplacental cone (EPC), while mural TE differentiates into trophoblast giant cells (TGC). TGC cells are in closest contact to the uterus during embryo implantation . Thus, the interaction between TGC and PDZ represents the key mechanism of embryo implantation. Notably, PDZ cells expressed Ptn and Mdk, while the corresponding receptor Ncl was expressed in TGC. Currently, little is known about the role of Ptn and Mdk in regulating embryo implantation. Additionally, we found that Fgf2 and Angptl4 expressed in PDZ might function via their receptors Fgfr2 and Cdh5 expressed in TGC. In mouse uterus, Fgf2 in the superficial stromal cells maintains proliferation of epithelial cells . Because Fgfr2 is also expressed in TGC, Fgf2 in PDZ might contribute to the proliferation of TGC. Angptl4 acts as a secretory protein to regulate angiogenesis in various tissues including the uterus . Vascular endothelial cadherin Cdh5 is also expressed trophoblast cells . Our findings supported the findings showing that predecidual stromal cells have distinctive characteristics of pericytes  and trophoblasts mimic the endothelial cells . On the other hand, TGC cells expressed Pdgfa, while the corresponding receptors Pdgfra and Pdgfrb were expressed in PDZ. Indeed, proteome profiling confirmed the presence of Pdgfa in trophoblast supernatant . Pdgfra and Pdgfrb are regarded as canonical cell markers for uterine stromal cells . Therefore, TGC might regulate proliferation and differentiation of PDZ via Pdgfa. Apart from secreted signaling, we found that TGC might also crosstalk with PDZ through ECM-receptor interaction. Our data provided clues for the mechanisms underlying embryo implantation from the aspect of cell–cell communication.
Besides the emergence of PDZ, by comparing IS with IIS, we found apparent cell type abundance changes upon embryo implantation. In particular, the proliferating subsets of S1, S2, NK, M, DC and VEC were significantly increased. There were also massive gene expression changes in these cell types. Considering spatial relationships between cell types, uterine cell types other than PDZ are unable to directly communicate with embryonic cells during embryo implantation. Therefore, these changes were likely caused by PDZ indirectly via secreted signaling. Indeed, we found that many soluble factors from PDZ, such as Wnt4, Wnt5a and Wnt6, might have an influence on S1 and S2 cells. Wnt4 is most abundant in the decidual cells. In uterus-specific Wnt4 knockout mice, the embryos were able to attach to the uterine luminal epithelium, but they failed to successfully invade into the uterine stroma . During embryo implantation, Wnt5a is highly localized in stromal cells at the mesometrial pole. Mice with uterine inactivation of Wnt5a show impaired embryo implantation and decidualization . Using Wnt6-null mice, it has been shown that Wnt6 is critical for normal stromal cell proliferation during decidualization . The idea that Wnt signaling might act in a paracrine way between PDZ and other stromal cells deserves further investigation. Additionally, PDZ-secreted Tgfb2 and Tgfb3 might target immune cells. Immune cells accumulate at embryo implantation sites and some of them such as DC cells play an important role in embryo implantation . Increasing evidence supports the involvement of TGFβ signaling in uterine function including embryo implantation . However, whether the TGFβ-mediated PDZ-DC interaction is functional during embryo implantation is still unknown. Moreover, we found that PDZ-secreted Vegfa, Vegfb and Vegfd might target endothelial cells. PDZ is avascular, but angiogenesis is accompanied with the emergence of SDZ. It has been well established that uterine immune cells are responsible to angiogenesis in SDZ . Thus, we suggested that PDZ, by synthesizing VEGFs, might be another inducer of angiogenesis besides immune cells. Altogether, our data highlighted the role of PDZ as a key mediator for uterine response in different cell types during embryo implantation.
In conclusion, this study provided a comprehensive single-cell transcriptome atlas for mouse uterus at the invasion phase of embryo implantation. Our data present a valuable resource for deciphering the molecular mechanism underlying embryo implantation.
Availability of data and materials
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.
Lee KY, DeMayo FJ. Animal models of implantation. Reproduction. 2004;128:679–95.
Ferretti C, Bruni L, Dangles-Marie V, Pecking AP, Bellet D. Molecular circuits shared by placental and cancer cells, and their implications in the proliferative, invasive and migratory capacities of trophoblasts. Hum Reprod Update. 2007;13:121–41.
Moreno-Moya JM, Franchi NA, Martinez-Escribano S, Martinez-Conejero JA, Bocca S, Oehninger S, et al. Transcriptome of early embryonic invasion at implantation sites in a murine model. Reprod Fert Develop. 2015;9:78.
Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28:511–5.
Hilton HG, Rubinstein ND, Janki P, Ireland AT, Bernstein N, Fong NL, et al. Single-cell transcriptomics of the naked mole-rat reveals unexpected features of mammalian immunity. PLoS Biol. 2019;17:e3000528.
Gendron RL, Paradis H, Hsieh-Li HM, Lee DW, Potter SS, Markoff E. Abnormal uterine stromal and glandular function associated with maternal reproductive defects in Hoxa-11 null mice. Biol Reprod. 1997;56:1097–105.
Kirkwood PM, Gibson DA, Smith JR, Wilson-Kanamori JR, Kelepouri O, Esnal-Zufiaurre A, et al. Single-cell RNA sequencing redefines the mesenchymal cell landscape of mouse endometrium. FASEB J. 2021;35:e21285.
Fu DJ, De Micheli AJ, Bidarimath M, Ellenson LH, Cosgrove BD, Flesken-Nikitin A, et al. Cells expressing PAX8 are the main source of homeostatic regeneration of adult mouse endometrial epithelium and give rise to serous endometrial carcinoma. Dis Models Mech. 2020;13:8.
Franco HL, Dai D, Lee KY, Rubel CA, Roop D, Boerboom D, et al. WNT4 is a key regulator of normal postnatal uterine development and progesterone signaling during embryo implantation and decidualization in the mouse. FASEB J. 2011;25:1176–87.
Cheng S, Pei Y, He L, Peng G, Reinius B, Tam PPL, et al. Single-Cell RNA-Seq Reveals Cellular Heterogeneity of Pluripotency Transition and X Chromosome Dynamics during Early Mouse Development. Cell Rep. 2019;26:2593–607.
Tan J, Raja S, Davis MK, Tawfik O, Dey SK, Das SK. Evidence for coordinated interaction of cyclin D3 with p21 and cdk6 in directing the development of uterine stromal cell decidualization and polyploidy during implantation. Mech Dev. 2002;111:99–113.
Paria BC, Zhao X, Das SK, Dey SK, Yoshinaga K. Zonula occludens-1 and E-cadherin are coordinately expressed in the mouse uterus with the initiation of implantation and decidualization. Dev Biol. 1999;208:488–501.
Alva JA, Zovein AC, Monvoisin A, Murphy T, Salazar A, Harvey NL, et al. VE-Cadherin-Cre-recombinase transgenic mouse: a tool for lineage analysis and gene deletion in endothelial cells. Dev Dyn. 2006;235:759–67.
Munoz-Fernandez R, de la Mata C, Prados A, Perea A, Ruiz-Magana MJ, Llorca T, et al. Human predecidual stromal cells have distinctive characteristics of pericytes: Cell contractility, chemotactic activity, and expression of pericyte markers and angiogenic factors. Placenta. 2018;61:39–47.
Lyall F, Bulmer JN, Duffie E, Cousins F, Theriault A, Robson SC. Human trophoblast invasion and spiral artery transformation: the role of PECAM-1 in normal pregnancy, preeclampsia, and fetal growth restriction. Am J Pathol. 2001;158:1713–21.
Gellersen B, Wolf A, Kruse M, Schwenke M, Bamberger AM. Human endometrial stromal cell-trophoblast interactions: mutual stimulation of chemotactic migration and promigratory roles of cell surface molecules CD82 and CEACAM1. Biol Reprod. 2013;88:80.
This research was funded by National Natural Science Foundation of China (32070845 and 31771665), Guangdong Natural Science Funds for Distinguished Young Scholars (2021B1515020079), Innovation Team Project of Guangdong University (2019KCXTD001), and Guangdong Special Support Program (2019BT02Y276).
Authors and Affiliations
Guangdong Laboratory for Lingnan Modern Agriculture, College of Veterinary Medicine, South China Agricultural University, Tianhe District, No. 483 Wushan Road, Guangzhou, 510642, China
Jia-Peng He, Qing Tian, Qiu-Yang Zhu & Ji-Long Liu
J-LL supervised the study. J-LL designed the experiments. J-PH performed the experiments. J-PH, QT, Q-YZ and J-LL analyzed the data. J-PH and J-LL wrote the paper. All authors read and approved the final manuscript.
The complete list of differentially expressed genes in IS compared to IIS for all cell types (logFC > 0.25 and P < 0.05).
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/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.