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 [26] and immune cells expressing Ptprc [27] (Fig. 2B).
Hormone-responsive cells included epithelial cells expressing Epcam and Krt19 [28] (Fig. 2C), stromal cells expressing Hoxa11 [29] (Fig. 2D), smooth muscle cells expressing Acta2 [11] and pericytes expressing Rgs5 [30] (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 [31]. 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 [32]. S3 was primary decidual zone stromal cells expressing decidualization marker genes Wnt4 [33] (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 [26] (Fig. 2H). There were 9 immune cell clusters (Fig. 2I, J). Included were natural killer cells (NK, Ptprc+Nkg7+Cd3e−) [27], proliferating natural killer cells (NKp, Ptprc+Nkg7−Cd3e+Mki67+), T cells (T, Ptprc+Nkg7−Cd3e+) [27], B cells (B, Ptprc+Cd79a+Ms4a1+) [27], macrophages (M, Ptprc+Adgre1+) [34], proliferating macrophages (M, Ptprc+Adgre1+Mki67+), dendritic cells (DC, Ptprc+Itgax+) [34], proliferating dendritic cells (DC, Ptprc+Itgax+Mki67+) and plasmacytoid dendritic cells (pDC, Ptprc+Siglech+) [35].
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 [36]. 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 [37]. 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 [38], epiblast (EPI) expressing Pou5f1 and Nanog [39], extraembryonic ectoderm or ectoplacental cone (EXE/EPC) expressing Elf5 and Cdx2 [39], and trophoblast giant cells (TGC) expressing Gata2 [37] and Prl3d1 [40] (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 [41].
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).
We predicted the ligand-receptor interactions between PDZ (S3) and VEC/LEC. We found a total of 34 ligand-receptor interaction pairs (Fig. 10A). Pathway analysis revealed that these ligand-receptor interactions were enriched among PI3K-Akt signaling pathway (FDR = 1.00 × 10–23), Rap1 signaling pathway (FDR = 1.00 × 10–16), Ras signaling pathway (FDR = 1.00 × 10–14), Cytokine-cytokine receptor interaction (FDR = 1.00 × 10–14), HIF-1 signaling pathway (FDR = 1.26 × 10–7), Hippo signaling pathway (FDR = 1.00 × 10–6), Endocytosis (FDR = 3.16 × 10–6), Regulation of actin cytoskeleton (FDR = 6.31 × 10–6), mTOR signaling pathway (FDR = 1.58 × 10–5), MAPK signaling pathway (FDR = 1.58 × 10–5), ECM-receptor interaction (FDR = 3.98 × 10–4), TGF-beta signaling pathway (FDR = 3.98 × 10–4), NOD-like receptor signaling pathway (FDR = 5.01 × 10–3), Toll-like receptor signaling pathway (FDR = 1.26 × 10–2), AMPK signaling pathway (FDR = 2.00 × 10–2), Wnt signaling pathway (FDR = 3.16 × 10–2) and Jak-STAT signaling pathway (FDR = 3.98 × 10–2) (Fig. 10B).