Identification of Intercellular Crosstalk between Decidual Cells and Niche Cells in Mice

Decidualization is a crucial step for human reproduction, which is a prerequisite for embryo implantation, placentation and pregnancy maintenance. Despite rapid advances over recent years, the molecular mechanism underlying decidualization remains poorly understood. Here, we used the mouse as an animal model and generated a single-cell transcriptomic atlas of a mouse uterus during decidualization. By analyzing the undecidualized inter-implantation site of the uterus as a control, we were able to identify global gene expression changes associated with decidualization in each cell type. Additionally, we identified intercellular crosstalk between decidual cells and niche cells, including immune cells, endothelial cells and trophoblast cells. Our data provide a valuable resource for deciphering the molecular mechanism underlying decidualization.


Introduction
Decidualization is a crucial step for human reproduction. It is a prerequisite for embryo implantation, placentation and pregnancy maintenance [1]. Decidualization initiates spontaneously in the secretory phase of the menstrual cycle, which is tightly controlled by ovarian estrogen and progesterone [2]. Upon decidualization, uterine stromal cells change into large epithelioid cells and secrete two protein markers, decidual prolactin (dPRL) and insulin-like growth factor-binding protein 1 (IGFBP1) [3]. Impaired decidualization may lead to intrauterine growth restriction (IUGR), repeated pregnancy loss (RPL) and severe pre-eclampsia (sPE) [4]. Therefore, it is imperative to understand the molecular mechanism of decidualization.
Due to ethical restrictions, studies on human decidualization are limited to the in vitro level. The in vivo investigation of decidualization heavily relies on mice. Slightly different from humans, the decidual reaction in mice is an embryo-dependent process. Prolactin family 8 subfamily a member 2 (Prl8a2), a paralog of human dPRL, is a well-established marker for mouse decidualization [5]. By using uterus-specific gene knockout mice, a number of genes have been proved to be required for decidualization [6]. Alternatively, several studies have analyzed global gene expression changes associated with decidualization by using high-throughput transcriptomic approaches [5,[7][8][9][10]. The limitation of these studies is that the bulk tissue was used. The decidua is a complex tissue consisting of many cell types, including decidual cells, stromal cells, endothelial cells and various immune cells. Thus, bulk-tissue methods were unable to accurately capture cell-type-specific gene expression changes. Moreover, the contribution of crosstalk between different cell types was ignored.
With advances in the single-cell RNA-seq techniques, it is now possible to analyze the global gene expression profile within highly heterogeneous tissues at a single-cell level [11]. In the present study, by using the state-of-the-art single-cell RNA-seq approach, we resolved all the cell types at the implantation site (decidualized) and the inter-implantation we resolved all the cell types at the implantation site (decidualized) and the inter-implantation site (undecidualized, served as control) of the mouse uterus on gestational day eight. Consequently, we were able to identify differential expression changes associated with decidualization in all the cell types. Additionally, we predicted intercellular crosstalk between the decidual cells and the niche cells. Our study provides a valuable resource for understanding the molecular mechanisms underlying decidualization.

A Single-Cell Atlas of Mouse Uterus on Gestational Day Eight
To create a cell-type resolved map of a mouse uterus upon decidualization, we performed single-cell RNA-seq analysis ( Figure 1A). The implantation site (IS, decidualized uterus) and the inter-implantation site (IIS, undecidualized uterus, served as a control) were collected from gestational day eight (GD8) ( Figure 1B). The whole uterus, which consists of the endometrium, the myometrium and the perimetrium, was subjected to single-cell dissociation ( Figure 1C). The embryo at the IS was also kept. Single-cell RNA-seq data were generated by using the 10× Genomics platform. After quality control, a total of 15087 cells (6793 for the IIS and 8294 for the IS) were obtained ( Figure 1D,E). In order to validate this single-cell RNA-seq dataset, we also generated a bulk RNA-seq dataset using the same samples ( Figure 1F). The cell-averaged single-cell RNA-seq data were highly accordant with the conventional bulk RNA-seq data (r = 0.7465 for the IIS and r = 0.7653 for the IS), indicative of the high quality of our single-cell RNA-seq data. Cells in green were live and cells in yellow were dead. A representative photo was provided, and a bar plot was shown on the right. n = 3; * p < 0.05. (D,E) Single-cell RNA-seq data pre-processing and quality control. Cells with detected genes of fewer than 200 or more than 6000 were removed. Only cells with total mitochondrial gene expression below 25% were kept. (F) Scatter plots showing the correlation between single-cell RNA-seq and bulk RNA-seq. For single-cell RNA-seq data, gene expression levels were averaged and normalized as transcript per million (TPM). For bulk RNA-seq data, gene expression levels were measured as fragments per kilobase per million (FPKM). Cells were stained with the AO/PI solution. Cells in green were live and cells in yellow were dead. A representative photo was provided, and a bar plot was shown on the right. n = 3; * p < 0.05. (D,E) Single-cell RNA-seq data pre-processing and quality control. Cells with detected genes of fewer than 200 or more than 6000 were removed. Only cells with total mitochondrial gene expression below 25% were kept. (F) Scatter plots showing the correlation between single-cell RNA-seq and bulk RNA-seq. For single-cell RNA-seq data, gene expression levels were averaged and normalized as transcript per million (TPM). For bulk RNA-seq data, gene expression levels were measured as fragments per kilobase per million (FPKM).  Hormone-responsive cells included epithelial cells expressing Epcam and Krt19 [16] ( Figure 2C), stromal cells expressing Hoxa11 [17] (Figure 2D), smooth muscle cells expressing Acta2 [18] and pericytes expressing Rgs5 [19] ( Figure 2E). We found two epithelial cell clusters, LE and GE. LE was luminal epithelial cells expressing Tacstd2 and GE was glandular epithelial cells expressing Foxa2 [20]. We identified five stromal cell clusters, S1, S1p, S2, S3 and S3p. The cells in S2 but not S1 expressed high levels of Hand2, implying that S2 was superficial stromal cells and S1 was deep stromal cells [21]. S1p was a subset of proliferating S1 with a high level of Mki67 ( Figure 2G). Notably, there was no proliferating subset for S2. S3 was decidual cells expressing decidualization marker genes Wnt4, Bmp2 and Prl8a2 [22][23][24]. S3p was intermediate decidual cells expressing proliferation marker gene Mki67 ( Figure 2G). Only one smooth muscle cell cluster and one pericyte cluster were found.
Finally, we aimed to discover novel markers for each cell type. We selected the genes that were expressed significantly higher in the cell type of interest than the other cell types using the Wilcoxon rank sum test. A complete list of marker genes is presented in Table S1.
We are particularly interested in GeneSet#4, the S3-enriched genes. On GD8, the entire IS, from the anti-mesometrial region to the mesometrial region, is fully decidualized. Meanwhile, the IIS shows no sign of decidualization. Thus, the S3-enriched genes could be validated using qTR-PCR with bulk tissues. The top 12 S3-enriched genes (A2m, Dmkn, Mt4, Tdo2, Mt3, Rrm2, Dio3, Serpinb6b, Ass1, Cystm1, Psca and Sfrp5) were selected ( Figure S1A). Although there was variation in the fold changes between the qRT-qPCR and the bulk RNA-seq, the up-regulation expression trends for all the genes tested were coincident between these two techniques ( Figure S1B). These data provide the validity of our single-cell RNA-seq data.

Cell-Cell Communication between Decidual Cells and Immune Cells
We investigated the abundance of each immune cell type at the IS compared to the IIS. For each cell type, proliferating and non-proliferating cells were summed. The χ 2 test was employed to assess the significance of difference between the two groups. By using the criteria of p < 0.05 and fold change > 2, the proportions of all the cell types are unchanged, except DC and B cells ( Figure 4A).

Cell-Cell Communication between Decidual Cells and Immune Cells
We investigated the abundance of each immune cell type at the IS compared to the IIS. For each cell type, proliferating and non-proliferating cells were summed. The χ 2 test was employed to assess the significance of difference between the two groups. By using the criteria of p < 0.05 and fold change > 2, the proportions of all the cell types are unchanged, except DC and B cells ( Figure 4A).
We investigated the breadth of transcriptional changes in each immune cell type by performing differential gene expression analysis. The proliferating cells were excluded. Using a logFC cutoff of 0.25 and a p value cutoff of 0.05, we identified 278, 380, 123, 181, 432 and 270 differentially expressed genes for T, B, NK, M DC and pDC, respectively (Figure 4B and Table S3). We then explored the biological implications of differentially expressed genes using gene ontology (GO) analysis. A complete list of the enriched GO terms is provided in Figure 4C. These data indicated that each immune cell type invokes distinct biological processes upon decidualization. We then used the CellChat software [28] to predict the ligand-receptor interactions between the decidual cells and the immune cells. We found a total of 68 ligand-receptor interaction pairs ( Figure 5A). Pathway analysis revealed that these ligand-receptor interactions were enriched among the ECM-receptor interaction (FDR = 1.00 × 10 −39 ), the PI3K-Akt signaling pathway (FDR = 1.00 × 10 −34 ), the Cytokine-cytokine receptor interaction (FDR = 1.00 × 10 −16 ), the regulation of actin cytoskeleton (FDR = 1.00 × 10 −14 ), the MAPK We investigated the breadth of transcriptional changes in each immune cell type by performing differential gene expression analysis. The proliferating cells were excluded. Using a logFC cutoff of 0.25 and a p value cutoff of 0.05, we identified 278, 380, 123, 181, 432 and 270 differentially expressed genes for T, B, NK, M DC and pDC, respectively ( Figure 4B and Table S3). We then explored the biological implications of differentially expressed genes using gene ontology (GO) analysis. A complete list of the enriched GO terms is provided in Figure 4C. These data indicated that each immune cell type invokes distinct biological processes upon decidualization.

Cell-Cell Communication between Decidual Cells and Endothelial Cells
By using the criteria of p < 0.05 and fold change > 2, we found that the proportion of LEC was unchanged, whereas the proportions of VEC were significantly increased (Figure

Cell-Cell Communication between Decidual Cells and Trophoblast Cells
Trophoblast cells may interact with decidual cells and pay an important role during the decidualization process. However, due to the relatively small number of embryonic cells at the implantation site, we did not find any fetal cell clusters in our single-cell RNAseq data. Alternatively, we re-analyzed a published single-cell RNA-seq dataset on mouse E7.5 embryos (GSM3457437) [29]. E7.5 is equivalent to gestational day eight in our study. By using marker gene Plac1 [30], we discovered a trophoblast cell cluster ( Figure 8A). These cells were further divided into ectoplacental cone (EPC) expressing Dlx3 [31] and trophoblast giant cells (TGC) expressing Prl3d1 [32] (Figure 8B).

Discussion
We have recently established a single-cell atlas of mouse uterus during embryo implantation using the 10× Genomics approach [33]. In this study, we systematically characterized the mouse uterus on GD8 at a single-cell resolution. We identified 21 distinct cell clusters in mouse uterus during decidualization, including 5 stromal cell clusters, 2 epithelial cell clusters, 1 smooth muscle cell cluster, 1 pericyte cluster, 4 endothelial cell clusters and 8 immune cell clusters. To our knowledge, our study is the first comprehensive single-cell transcriptomic atlas for mouse decidualization.
We used 2 mg/mL Collagenase II and 10 mg/mL Dispase II for a single-cell suspension preparation. Collagenase II is a crude collagenase preparation with weak trypsin-like activity and, thus, trypsin was not used in our method. It turned out that the cell viability was >80% and the percentage of cell clumps was <10%, indicative of the high efficiency of our method. Nevertheless, we noticed that the percentage of decidualized stromal cells was much lower than expected: S3 (decidual cells) was 4.4% and S3p (intermediate proliferating decidual cells) was 1.7%. As we observed lower cell viability at the implantation site compared to the inter-implantation site, we suspected that the collagenase treatment might selectively destroy the decidualized cells. This was in line with the finding that the highly expressed decidual-specific marker gene Prl8a2 contaminated all the other cells, indicative of the rupture of the decidual cells during the single-cell suspension preparation.
Moreover, we would like to note that the collagenase treatment might cause transcriptional disturbances, leading to artifacts [34]. We compared single-cell RNA-seq data against parallel bulk RNA-seq data. We found that a large number of genes were changed in both the IS and the IIS ( Figure S2A). The complete list of differentially expressed genes is listed in Tables S4 and S5. There were 1590 down-regulated genes and 1297 up-regulated genes shared by the IS and the IIS ( Figure S2B). We focused on up-regulated genes because we suspected that these genes might exhibit de novo expression upon collagenase treatment. We examined the expression pattern of 12 selected up-regulated genes in single-cell RNAseq data. We found that Egr1, Egr2, Egr3, Fos, Jun and Atf3 were ubiquitously expressed in nearly all cell types, whereas Ccl6, Ccl9, Ccl5, Ccl8, Ar and Tagln were expressed in a cell-type-restricted manner ( Figure S2C). We believe that these genes could be used as marker genes without any problem, although they are similar to "artificial markers". However, it will be problematic if they are chosen as differentially expressed genes, because the de novo expression levels might vary across each collagenase treatment. Our data suggest that single-cell RNA-seq data should be interpreted with caution.
In this study, the embryo at the IS was kept. We estimated that there were around 2000 cells per embryo. We usually obtained two million cells per sample in the single-cell dissociation procedure. Approximately 5000 cells were sequenced per sample by the 10× platform. By calculating the probability, we found that only~5 cell per embryo could be captured in our single-cell RNA-seq data. This was in line with our findings that there were no embryo-derived cell clusters in our single-cell RNA-seq data. The embryo-derived cells, if they existed at all, might be scattered among uterine cells in our single-cell RNA-seq data. Due to their small number, these cells might not have an obvious effect on the assigning cell type label for each cell cluster. Moreover, these cells might not interfere with the process for identifying differentially expressed genes either, because we set min.pct to 0.20, which means that we required the differentially expressed genes to be expressed in a minimum of 20% of cells within the cluster of interest.
In mice, embryo implantation initiates at midnight on GD4. Upon embryo implantation, the superficial stromal cells proximally surrounding the implantation chamber exhibit proliferation and spreading. These cells differentiate and form the primary decidual zone (PDZ) on the afternoon of gestational day five. The PDZ subsequently undergoes apoptosis and disappears by GD8. Meanwhile, the stromal cells next to the PDZ continue to proliferate and differentiate into decidual cells, forming the secondary decidual zone (SDZ). The SDZ is full developed on GD8. Here, we found five stromal cell clusters for the IS (S2, S1, S1p, S3 and S3p) and three stromal cell clusters for the IIS (S2, S1 and S1p). S2 was superficial stromal cells and S1/S1p was deep stromal cells, which is in line with our previous study [33]. S3 and S3p were IS-specific: S3 was decidual cells expressing the decidualization marker genes Wnt4, Bmp2 and Prl8a2, while S3p was intermediate decidual cells expressing the proliferation marker gene Mki67. Through pseudotime trajectory analysis, we found that S3p/S3 decidual cells originated from S1/S1p, but not from S2 cells. By using the Wilcoxon rank sum test, we identified two gene sets (GeneSet#3 and GeneSet#4) associated with decidualization. GeneSet#3 was gradually decreased during the decidualization process, whereas GeneSet#4 was gradually increased during the decidualization process. Based on GO analysis, the enriched terms for GeneSet#3 were cell cycle and proliferation, protein metabolism, cell organization and biogenesis, DNA metabolism and transport. The enriched GO terms for GeneSet#4 were transport, protein metabolism and cell cycle and proliferation. Our data provided clues for the development trajectory for decidualized stromal cells in SDZ.
We are particularly interested in GeneSet#4, which may represent the key molecular mechanism of decidualization. The best strategy to validate our findings is to obtain a pure cell population by using fluorescence activated cell sorting (FACS). Previously, CD13-specific antibodies were used to sort stromal cells by FACS from human endometrial biopsies [35,36]. However, the same method to sort decidual cells is yet to be established. GeneSet#4 was S3-enriched genes. On GD8, the IS but not the IIS is fully decidualized. Thus, the genes in GeneSet#4 should be up-regulated in the IS compared to the IIS at the bulk-tissue level. In this way, we were able to validate GeneSet#4, at least in part, by using qRT-PCR. The top 12 genes were tested. It turned out that all these genes were up-regulated in the IS compared to the IIS. This result might provide the validity of our single-cell RNA-seq data.
Immune cells play an important role in embryo implantation and decidualization in mice [37]. We investigated the abundance of each immune cell type at the IS compared to the IIS. By using the criteria of p < 0.05 and fold change > 2, the proportions of all the cell types are unchanged, except for DC and B cells. As B cells were a minor cell type (< 5%) in the uterus, we thus focused on DC cells. It has been demonstrated that the depletion of DC cells resulted in a severe impairment of the decidualization process, leading to embryo resorption [38]. We found that S3 decidual cells expressed secreted factors such as Tgfb2, Tgfa, Ptn, Mif, Mdk, Lgals9, Igf1, Fgf2 and C3, while their receptors were expressed on DC cells. Thus, these secreted factors might attract DC cells, leading to the accumulation of DC cells in the decidualized implantation site of the uterus. On the other hand, DC cells expressed secreted factors Tnf, Tgfb1, Nampt and Mif. As their receptors were expressed on S3 decidual cells, these secreted factors might contribute to the decidualization reaction of stromal cells.
Endothelial cells are in the niche of decidual cells. We hypothesized that endothelial cells might contribute to decidualization. We found that VEC cells were significantly increased, while LEC cells were slightly decreased upon decidualization. This was in line with the fact that the decidualization process is accompanied by angiogenesis [39]. Though cell-cell communication analysis, we found that angiogenesis was likely mediated by Vegfa and Vegfb expressed on S3 decidual cells and their receptors Vegfr1 and Vegfr2 expressed on VEC cells. VEC cells might promote decidualization by expressing soluble factors such as Hbegf. The contribution of VEC to decidualization deserves further investigation.
Previously, E7.5 embryos obtained from the GD8 uterus have been subjected to singlecell RNA-seq [29]. Through data integration, we inferred cell-cell communication between E7.5 embryos and the implantation site of the GD8 uterus by their expression of ligandreceptor pairs. We generated an interaction network for decidual cells from the uterus and EPC/TGC cells from the embryo. These interactions might provide clues for the physical function of decidual cells in placentation and fetal development.
In summary, we revealed intercellular crosstalk between decidual cells and niche cells, including immune cells, endothelial cells and trophoblast cells. This cell atlas of mouse uterus on GD8 provides an essential resource for understanding the molecular mechanism underlying decidualization.

Sample Collection
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 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 GD8, the implantation sites (decidualized uterus) and non-implantation sites (undecidualized uterus, served as a control) were collected separately.

Bulk RNA-Seq Analysis
The total RNA from uterine tissues were extracted with the TRIzol reagent (Invitrogen). RNA-seq libraries were generated by using the TruSeq RNA sample preparation kit (Illumina, San Diego, CA, USA) and sequenced on an Illumina HiSeq 2500 system with the paired-end 150-bp protocol. Raw data were analyzed using TopHat v2.0.4 [40] and Cufflinks v2.2.1 [41]. The mouse genome UCSC mm10 was used for reference. Gene expression levels were measured as fragments per kilobase per million (FPKM).

Single-Cell Dissociation of Mouse Uterus
Uterine samples from 3 mice for each group were pooled and minced with a blade. Samples were then incubated in dissociation buffer containing 2 mg/mL Collagenase II (#C6885, Sigma-Aldrich, St. Louis, MO, USA), 10 mg/mL Dispase II (#354235, Corning, Corning, NY, USA) and 50,000 U/mL DNase I (#DN25, Sigma-Aldrich, St. Louis, MO, USA) 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-micrometer 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, Lawrence, MA, USA) and counted using a Cellometer Auto 2000 instrument (#SD-100, Nexcelom Bioscience, Lawrence, MA, USA). 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 Chromium TM Single Cell B Chip (#1000073, 10× Genomics, Pleasanton, CA, USA), aiming at recovering 8000-10,000 cells. The Chromium Single Cell 3 Library and Gel Bead Kit v3 (#1000075, 10× Genomics, Pleasanton, CA, USA) was used for single-cell bar-coding, cDNA synthesis and library preparation, following the manufacturer's instructions. 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 (v2.19.0.316, Illumina, San Diego, CA, USA). These fastq files were aligned to the UCSC mm10 mouse reference genome by using the CellRanger software (v3.0.1, 10× Genomics, Pleasanton, CA, USA). The resulting gene counts matrix was analyzed with the R package Seurat (v3.1.3) [42]. 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 the NormalizeData and ScaleData functions, respectively. The top 2000 highest variable genes were used for the principal component analysis (PCA) and the optimal number of PCA components was determined using the JackStraw procedure. Single cells were clustered using 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 min.logfc being set to 0.25 and min.pct being set to 0.20. By using the same parameters, the FindMarkers function was used to find differential expressed genes in the IS compared to the IIS for each cell type.

Pseudotime Analysis
The Monocle2 package v2.18.0 was used for pseudotime analysis [27]. The count data and meta data were exported from the Seurat object and then imported into the CellDataSet object in Monocle2. Feature genes were selected by using the differentialGeneTest function. After dimension reduction by using the DDRTree algorithm, the orderCells function was used to infer the trajectory with default parameters. The reconstructed trajectory was visualized using the plot_cell_trajectory function.

Gene Ontology Analysis
Gene ontology (GO) analysis was performed as described previously [43]. GO terms were grouped according to the biological process category defined in the Mouse Genome Informatics (MGI) GOslim database [44]. To test for enrichment, a hypergeometric test was conducted and p < 0.05 was used as the significance threshold to identify enriched GO terms.

Pathway Analysis
Pathway enrichment analysis was conducted by using the Metascape v7.4 online tools [45]. The significance threshold for FDR was set at 0.05.

Cell-Cell Communication Analysis
The CellChat v1.1.0 R package [28] 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 using a permutation test. p < 0.05 was considered significant.

Validation by Quantitative RT-PCR
The TRIzol reagent (Invitrogen, Carlsbad, CA, USA) was used to extract total RNAs. The cDNAs were synthesized with the PrimeScript reverse transcriptase reagent kit (TaKaRa, Dalian, China). A quantitative RT-PCR was performed on Applied Biosys-tems 7500 (Life Technologies, Carlsbad, CA, USA) with the THUNDERBIRD SYBR qPCR Mix (Toyobo, Osaka, Japan). The Rpl7 gene was used as the reference gene for data normalization. A complete list of primer sequences is provided in Table S6.