Deciphering mouse uterine receptivity for embryo implantation at single‐cell resolution

2.1. Sample collection

Adult CD‐1 mice of the specific‐pathogen free (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 the mating was confirmed the next morning by the presence of a vaginal plug. The day of the vaginal plug was denoted as gestation day (GD) 1. The whole uterus was obtained at 0900 h on GD3 (pre‐receptive stage) and GD4 (receptive stage) respectively. Success of breeding and early embryo development was further confirmed by recovering morula from the oviduct on GD3. On GD4, one horn of the uterus was flushed with saline and early embryo development was further confirmed by recovering blastocysts from the uterus. The other intact uterine horn was used for sample collection. All animal procedures were approved by the Institutional Animal Care and Use Committee of South China Agricultural University.

2.2. Haematoxylin and eosin staining

Uterine tissues were fixed in 4% paraformaldehyde solution for 24 h. After paraffin processing, tissues were cut into 6‐μm sections and stained with haematoxylin and eosin.

2.3. Immunohistochemistry

Paraformaldehyde‐fixed paraffin‐embedded uterine samples were cut into 6 μm sections. Antigen retrieval was performed by heating the slides in 10 mM citrate buffer for 10 min. Endogenous peroxidase activity was removed by using 3% H₂O₂. After blocking with 10% horse serum in PBS, sections were incubated with anti‐FOXA2 (1:200 dilution, #ab108422, Abcam) or normal rabbit IgG in 10% horse serum overnight at 4°C. The signal was developed by the 3,3'‐diaminobenzidine (DAB)‐HRP reaction system. Sections were counterstained with haematoxylin. The positive signal was visualized as a dark brown colour.

2.4. Bulk RNA‐seq analysis

The total RNA from uterine tissues was extracted with the TRIzol reagent (Invitrogen). RNA‐seq libraries were generated by using the TruSeq RNA sample preparation kit (Illumina) and sequenced on a HiSeq 2500 system (Illumina). Raw data were trimmed by the fastp program ¹⁷ with the parameters ‘‐w 4 ‐q 20 ‐u 50’ to obtain clean reads. Clean reads were mapped to mouse genome UCSC mm10 by using Hisat2 software v2.0.5 ¹⁸ with default parameters. HTSeq v0.13.5 ¹⁹ was subsequently employed to convert aligned short reads into read counts for each gene. Gene expression levels were normalized as transcripts per kilobase million (TPM) using an in‐house PERL script.

2.5. Single‐cell dissociation of mouse uterus

The uterine tissues from 3 mice for each group were pooled and minced with a blade. Tissues were then incubated in the dissociation buffer containing 2 mg/ml Collagenase II (#C6885, Sigma‐Aldrich), 10 mg/ml Dispase II (#354235, 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 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 stained 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%.

2.6. 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™ Single Cell B Chip (#1000073, 10x Genomics), aiming at recovering 8000–10,000 cells. The Chromium Single Cell 3’ Library & Gel Bead Kit v3 (#1000075, 10x 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 a NovaSeq 6000 system (Illumina) configured with the paired‐end 150 bp protocol at a sequencing depth of approximately 400 million reads.

2.7. Single‐cell RNA‐seq data analysis

Raw data of bcl files from the NovaSeq 6000 system were converted to fastq files using the bcl2fastq2 tool v2.19.0.316 (Illumina). These fastq files were aligned to the mm10 mouse reference genome by using the CellRanger software v3.0.1 (10x Genomics). The resulting gene counts matrix was processed with the R package Seurat v3.1.3. ²⁰ Cells 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 clean 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 graph‐based algorithm in PCA space and visualized using the t‐distributed stochastic neighbour 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 20% 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.pct being set to 0.20 and min.logfc being set to 0.25.

2.8. 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.

2.9. Gene network reconstruction

Gene network reconstruction was performed by using the STRING online tool v11.0. ²³  The threshold score for gene‐gene interaction was set to 0.9. The Cytoscape software v2.8.1 was used to display the network. ²⁴  The Network Analyzer plugin for Cytoscape was used to compute degree distribution. ²⁵ Hub genes in the network were selected by using a defined degree cut‐off value of mean + 2 × SD.

2.10. Cell‐cell communication analysis

The single‐cell RNA‐seq data for mouse E3.5 blastocysts were downloaded from the GEO database (GSM4026212) ²⁶ and re‐analysed with the same pipeline as described above. E3.5 is equivalent to GD4. In order to analyse the cell‐cell communication between blastocysts and uterus, we merged E3.5 blastocyst data and GD4 uterus data into a single Seurat object, from which a meta file as well as a count file was generated. These 2 files were used as input for the CellphoneDB software v2.1.4 to infer cell‐cell communications based on ligand‐receptor interactions with default parameters. ²⁷ p < 0.05 were considered significant.

3.1. Identification of cell types in mouse uterus

To create a cell‐type resolved map of uterine receptivity in mice, we performed single‐cell RNA‐seq analysis (Figure 1A). Pre‐receptive uterus and receptive uterus were collected from gestational days (GD) 3 and 4 respectively (Figure 1B). Success of breeding and early embryo development was confirmed by recovering morula from the oviduct on GD3 and blastocysts from the uterus on GD4 (Figure 1C). The whole uterus, which is consist of endometrium, myometrium and perimetrium, was used for single‐cell suspension preparation (Figure 1D). Single‐cell RNA‐seq data were generated by using the 10x Genomics platform. After quality control, a total of 9337 cells (4108 cells for GD3 uterus and 5229 cells for GD4 uterus) were obtained (Figure 1E). In order to validate this single‐cell RNA‐seq data set, we also generated a bulk RNA‐seq data set using the same samples. The cell‐averaged single‐cell RNA‐seq data were highly accordant with the conventional bulk RNA‐seq data (r = 0.8062 for GD3 and r = 0.7988 for GD4), indicative of high quality of our single‐cell RNA‐seq data (Figure 1F).

Unsupervised clustering analysis revealed 19 distinct cell clusters for all cells from GD3 and GD4 combined (Figure 2A). Major cell types were defined using the expression of known cell‐type‐specific genes, with stromal cells expressing Hoxa11, ¹¹ epithelial cells expressing Epcam and Krt19, ²⁸ smooth muscle cells expressing Acta2, ²⁹ pericytes expressing Rgs5, ³⁰ endothelial cells expressing Pecam1 ³¹ and immune cells expressing Ptprc. ³² We found 3 stromal cell clusters, S1, S2 and S1p. Cells in S1 but not in S2 expressed high levels of Hand2 (Figure 2B), implying that S1 was superficial stromal cells and S2 was deep stromal cells. ³³ S1p was linked to S1 and expressed high level of Mki67, suggesting that S1p was a subset of proliferating S1 cells. There were 2 epithelial cell clusters, E and Ep (Figure 2C). Ep was a subset of proliferating E cells with high expression of Mki67. Only one smooth muscle cell cluster (SMC) and one pericyte cluster (PC) were identified (Figure 2D). Endothelial cells have 4 clusters (Figure 2E): VEC and its proliferating subset VECp were vascular endothelial cells expressing Vwf, ³¹ while LEC and its proliferating subset LECp were lymphatic endothelial cells expressing Prox1. ³¹ The 8 immune cell clusters are macrophages (M, Ptprc⁺Adgre1^{+ 34} ), dendritic cells (DC, Ptprc⁺Adgre1⁻Itgax^{+ 34} ), plasmacytoid dendritic cells (pDC, Ptprc⁺Siglech^{+ 35} ), natural killer cells (NK, Ptprc⁺Nkg7⁺Cd3e^{− 36} ), natural killer T cells (NKT, Ptprc⁺Nkg7⁺Cd3e^{+ 36} ), T cells (T, Ptprc⁺Nkg7⁻Cd3e^{+ 36} ), B cells (B, Ptprc⁺Cd79a^{+ 36} ), a proliferating subset of mixed macrophages and dendritic cells (M/DCp, Ptprc⁺Adgre1⁺Mki67⁺ or Ptprc⁺Itgax⁺Mki67⁺), and a proliferating subset of mixed natural killer cells, natural killer T cells and T cells (NK/NKT/Tp, Ptprc⁺Nkg7⁺Mki67⁺ or Ptprc⁺Cd3e⁺Mki67⁺) (Figure 2F–H).

We then re‐analysed a single‐cell RNA‐seq data set for non‐pregnant adult mouse uterus from the mouse cell atlas project. ³⁷ Using the same computational pipeline, we identify LE, GE, S1, S2, SMC, PC, VEC, LEC, M, DC, pDC and NK. However, NKT, T or B were not found (Fig. S1). This analysis indicated that the markers used for receptive uterus can also be used for non‐pregnant adult mouse uterus.

We next 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 heat map depicting the top 10 marker genes for each cell type is shown in Fig. S2A. A complete list of marker genes was presented in Table S1. We picked up a single‐marker gene for all cell types except DC and SMC (Fig. S2B). We found that these novel marker genes were as good as, if not better than, canonical marker genes.

Finally, we searched for endometrial epithelial stem/progenitor cells and endometrial mesenchymal stem cells. The marker genes for endometrial epithelial stem/progenitor cells were Cd44, ³⁸ Lrg5, ³⁹ Fut4 ⁴⁰ and Aldh1a1, ⁴¹ while the marker genes for endometrial mesenchymal stem cells were Susd2 ⁴² and Ncam. ⁴³ We found that there was little overlap in expression pattern between these marker genes. We also found that these potential stem/progenitor cells were not clustered but rather scattered among other cells in the TSNE plots (Fig. S3). Further validations are needed to confirm our findings.

3.2. Cell type‐specific transcriptional changes associated with uterine receptivity

We investigated the abundance of each cell type (Figure 3A). The proportions of S2, SMC and PC were unchanged in GD4 uterus compared to GD3 uterus. The proportion of S1 was significantly increased, whereas the proportions of E and immune cells (M, DC, pDC, NK, NKT and T) were significantly decreased. Notably, although the proportion of VEC was unchanged, the proportion of LEC was significantly reduced in GD4 uterus compared to GD3 uterus.

We investigated the breadth of transcriptional changes in each cell type by performing differential gene expression analysis. Using a logFC cut‐off of 0.25 and a p value cut‐off of 0.05, we identified 739, 399, 173, 414, 367, 396, 412, 464, 130, 298, 135, 303, 332 and 306 differentially expressed genes for E, S1, S2, SMC, PC, VEC, LEC, M, DC, pDC, NK, NKT, T and B respectively (Figure 3B and Table S2). We then explored the biological implications of differentially expressed genes using gene ontology (GO) analysis. A complete list of enriched GO terms was provided in Figure 3C. These data indicated that each cell type invokes distinct biological processes in order to participate in the establishment of uterine receptivity.

3.3. Dissecting luminal and glandular epithelium from single‐cell RNA‐seq data

We focussed on uterine epithelium (Figure 4A), because it is the first maternal contact for embryo implantation. ⁴⁴ Using Tacstd2 as luminal epithelium (LE) marker and Foxa2 as glandular epithelium (LE) marker, ⁴⁵ we divided uterine epithelium into 4 cell clusters: LE, proliferating LE (LEp), GE and proliferating GE (GEp) (Figure 4B and Figure 4C). LEp and GEp were found only in GD3 uterus. We observed that Foxa2 was only expressed in a small portion of GE cells in our single‐cell RNA‐seq data. By using Immunohistochemistry, we found that Foxa2 was actually expressed in almost all GE cells (Figure 4D). This was likely a ‘dropout’ phenomenon commonly seen in single‐cell RNA‐seq experiments: some mRNA molecules might be lost due to their tiny initial amounts in the library preparation step, leading to false zero values. ⁴⁶

In LE, we found that a total of 749 genes were differentially expressed, of which 296 genes were down‐regulated and 453 genes were up‐regulated on GD4 in comparison to GD3 (Figure 5A and Table S3). Functional clustering analysis categorized the differentially expressed genes into 12 biological processes (Figure 5B): cell adhesion (2.3%), stress response (5.2%), cell death (4.9%), transport (9.8%), cell cycle & proliferation (7.2%), signal transduction (7.9%), cell‐cell signalling (1.0%), cell organization & biogenesis (10.2%), developmental processes (12.5%), protein metabolism (12.7%), DNA metabolism (1.4%), RNA metabolism (11.2%) and other biological processes (13.7%). Based on the enrichment test, 5 out of these GO terms, namely cell death (p = 0.00018), cell cycle & proliferation (p = 0.000015), cell organization & biogenesis (p = 0.0000049), developmental processes (p = 0.00011) and protein metabolism (p = 0.00028), were significantly enriched among differentially expressed genes (p < 0.01). Additionally, we reconstructed a network for differentially expressed genes by using the STRING online tool. This network had 242 genes with 518 edges (Figure 5C). Within this network, some nodes, known as hub genes, are highly connected compared to others. Using a defined cut‐off value, we identified 7 hub genes: Pabpc1, Ep300, Rhoa, Fos, Ppp2ca, Fau and Stat3.

In GE, we identified a total of 595 differentially expressed genes, of which 149 genes were down‐regulated and 446 genes were up‐regulated on GD4 compared to GD3 (Figure 6A and Table S4). According to gene ontology (GO), differentially expressed genes can be categorized into 12 biological processes (Figure 6B): cell adhesion (1.0%), stress response (4.2%), cell death (3.7%), transport (11.5%), cell cycle & proliferation (5.5%), signal transduction (6.4%), cell‐cell signalling (0.6%), cell organization & biogenesis (7.6%), developmental processes (8.9%), protein metabolism (14.6%), DNA metabolism (1.0%), RNA metabolism (8.6%) and other biological processes (26.4%). Based on hypergeometric test, transport (p = 0.00038) and protein metabolism (p = 0.00000015) were significantly enriched. We constructed a gene network for differently expressed genes, consisting of 139 nodes connected via 252 edges (Figure 6C). We found that Uba52, Rbx1, Akt1, Fau, Rhoa and Ywhae were hub genes of this network.

3.4. Validating single‐cell RNA‐seq data by laser capture microdissection‐coupled microarray data

Previously, microarray analysis was conducted on LE and GE isolated by laser capture microdissection (LCM) from uterus at 1600 h on GD3 and GD4 of pseudo‐pregnancy. ⁴⁷ In order to further validate our single‐cell RNA‐seq data set, a comparison analysis was performed. LCM microarray data were downloaded from GEO database under the accession number GSE48239. Fold change >2 was used to select differently expressed genes. In LE, we identified 33 overlap down‐regulated genes and 26 overlap up‐regulated genes respectively (Figure 7A). In GE, we found 12 overlap down‐regulated genes and 118 overlap up‐regulated genes respectively (Figure 7B). Notably, p < 0.05 was reached for all comparison, providing validity of our single‐cell RNA‐seq data.

3.5. Inferring cell‐cell communication between blastocysts and the receptive uterus

On GD4, the embryo, which is denoted as E3.5 blastocysts, enters the uterus from the oviduct. We re‐analysed a published 10x single‐cell RNA‐seq data set on mouse E3.5 blastocysts ²⁶ (Figure 8A). By using canonical marker genes, 3 major cell types were identified: trophectoderm (TE), inner cell mass/epiblast (ICM/EPI) and primitive endoderm (PE) (Figure 8B). To investigate the crosstalk between E3.5 blastocysts and GD4 uterus, we used CellPhoneDB to predict the ligand‐receptor interactions between distinct cell types (Figure 8C). Of special interest, 70 TE‐LE/TE‐GE ligand‐receptor interaction pairs were identified (Figure 8D). Considering spatial relationships between cell types, these ligand‐receptor interactions might play a crucial role in the establishment of uterine receptivity.

3.6. Comparison of uterine receptivity between mice and humans at single‐cell resolution

Finally, we compared our single‐cell RNA‐seq data with the a recently published single‐cell RNA‐seq data set which was generated from the pre‐receptive phase (day 17 of the menstrual cycle, D17) and the receptive phase (day 22 of the menstrual cycle, D22) of human endometrium by using the 10x Genomic method. ⁴⁸  The same computational pipeline was employed to re‐analyse this data set (Figure 9A). We identified a total of 12 cell clusters, namely luminal epithelial cells (LE), glandular epithelial cells (GE), ciliated epithelial cells (cE), proliferating mixed luminal and glandular epithelial cells (LE/GEp), stromal cells (S1), proliferating stromal cells (S1p), smooth muscle cells (SMC), proliferating smooth muscle cells (SMCp), mixed macrophages and dendritic cells (M/DC), mixed NK, NKT and T cells (NK/NKT/T), mixed vascular and lymphatic endothelial cells (VEC/LEC) and proliferating mixed vascular and lymphatic endothelial cells (VEC/LECp) (Figure 9B–G). We then attempted to compare differentially expressed genes associated with uterine receptivity between mice and humans in 4 stringently comparable cell types, namely LE, GE, S1 and SMC. For down‐regulated genes, a significant overlap was found in LE, GE and SMC, but not in S1 (Figure 9H). For up‐regulated genes, a significant overlap was found in all these 4 cell types (Figure 9I). These data highlighted the similarity in uterine receptivity between mice and humans.