Single-cell transcriptomic analysis of immune cell dynamics in the healthy human endometrium.

Kaixing Chen: Writing – original draft, Visualization, Validation, Software, Resources, Project administration, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. Qiaoni Yu: Writing – original draft, Visualization, Resources, Formal analysis, Data curation, Conceptualization. Qing Sha: Writing – review & editing, Resources, Investigation, Conceptualization. Junyu Wang: Writing – review & editing, Visualization, Supervision, Investigation, Conceptualization. Jingwen Fang: Writing – review & editing, Supervision, Investigation, Conceptualization. Xin Li: Supervision, Conceptualization. Xiaokun Shen: Supervision, Conceptualization. Binqing Fu: Supervision, Conceptualization. Chuang Guo: Writing – review & editing, Validation, Supervision, Resources, Project administration, Funding acquisition.

All authors consent to the publication of this article.

All the data used in our manuscript are sourced from publicly available datasets and all figures are all generated from data analysis. The datasets have been explicitly acknowledged in the Availability of data and materials section of our manuscript. Therefore, an additional Ethical Approval may not be required.

This work was supported by the 10.13039/501100001809 National Natural Science Foundation of China grants ( 32270978 to C.G.), We thank the USTC supercomputing center and the School of Life Science Bioinformatics Center for providing computing resources for this project.

To investigate the changes in the composition and transcriptome of endometrial immune cell populations throughout the reproductive cycle, we collected seven publicly available scRNA-seq datasets obtained from the endometrial samples of 52 total healthy women [ 1 , 5 , 9 , 11 , 12 , 14 , 15 ]. Specifically, samples from the proliferative (n = 4) and secretory (n = 11) stages were obtained on days 1–14 and 15–28 of the menstrual cycle, respectively. The average pregnancy duration of early pregnancy samples (n = 28) varied from 6 to 14 gestational weeks. Likewise, the pregnancy duration of late pregnancy samples (n = 9) ranged from 33 to 40 gestational weeks ( Fig. 1 A). After removing cell doublets and filtering out low-quality cells, we applied Seurat to integrate the scRNA-seq datasets of endometrial samples from different stages of the human female reproductive cycle. This process resulted in a combined total set of 230,049 single-cell transcriptomes, including 22,115 cells obtained in the proliferative stage, 67,098 cells from the secretory stage, 83,111 cells from early pregnancy, and 57,725 cells from late pregnancy. The resulting gene expression matrix was normalized, and a subsequent hierarchical clustering analysis revealed the presence of 47 distinct clusters, which were visualized using uniform manifold approximation and projection (UMAP) plots. Cell lineages were identified based on predominant markers ( Table S3 ) [ 9 ], including three NK cell subsets (NK1, NK2, and NK3 cells), two macrophage subsets (Mac1 and Mac2 cells), three T-cell subsets (CD4 + T, CD8 + T, and Treg cells), dendritic cells (DCs), plasma cells, granulocytes, and nonimmune cells (including stromal cells, endothelial cells, fibroblasts, epithelial cells, perivascular cells, and trophoblasts) ( Fig. 1 B–D, and Figs. S1A and B ). In addition, this analysis also detected proliferative cell subsets, such as proliferative NK cells (NKp), proliferative macrophages (Mpro), and proliferative T cells (Tpro). Comparison with cell types previously annotated in the literature showed remarkable consistency with the cell identities of populations in these datasets, thus confirming the validity of our cluster annotations [ 1 , 5 , 9 ] ( Figs. S2A–C ). In addition, stage-specific subsets of cells in our atlas were preserved through the cell-typing process. For example, endometrial stromal (eS) cells largely appeared in the proliferative stage, while decidual stromal (dS) cells only appeared in the secretory stage, early pregnancy, and late pregnancy ( Fig. 1 B, C and S3A ), which aligns well with the known time of differentiation from eS cells to dS cells [ 39 ]. Trophoblasts were also detected only in early and late pregnancy since they develop to initiate the invasion process after embryo implantation ( Fig. 1 B, C and Fig. S3A ). These collective results indicated that the single-cell transcriptome atlas of human endometrial cells from different stages of the reproductive cycle was reliable for further analysis of changes in immune cell composition, transcriptome, and intercellular interactions. The immune cells found in the endometrium play crucial roles in facilitating a successful pregnancy, including maintaining immune tolerance, regulating trophoblast invasion, promoting fetal growth, and fighting infections [ 7 , [40] , [41] , [42] ]. Thus, we next investigated the distribution of major immune cell subsets (NK cells, macrophages, and T cells) throughout the reproductive cycle. Notably, during early pregnancy, these cell subsets all increased significantly ( Fig. S3B ), indicating their importance in facilitating normal pregnancy, particularly implantation. We further investigated the compositional changes in specific cell subsets within each major immune cell subset. We found that the proportion of NK1 cells significantly increase during early pregnancy, while NK3 cells, which secrete high levels of cytokines, significantly increase in late pregnancy and maintain a high proportion with NK2 cells ( Fig. 1 E– S3C ), both of which were corroborated in Whettlock's study using flow cytometry [ 43 ]. Additionally, CD4 + T and CD8 + T cells were observed to increase significantly during late pregnancy ( Fig. S3D ). Regarding the macrophage population, Mac1 with pro-inflammatory characteristics predominates during the secretory stage, while the anti-inflammatory Mac2 accounted for a higher proportion in early pregnancy ( Figs. S4A and B ). Therefore, these results reflect that the anti-inflammatory environment which aids in the placenta formation in early pregnancy [ 3 , 44 ] comes from the contribution of multiple cells, such as an increase in anti-inflammatory macrophages and a decrease in pro-inflammatory lymphocytes such as NK3 and CD8 + T cells. As the most abundant immune cells in pregnancy, NK cells are the primary subject of our investigation into functional changes throughout the reproductive cycle [ 3 ]. To determine the gene expression differences between NK cell subsets across stages, we utilized consensus nonnegative matrix factorization (cNMF) and hierarchical clustering strategies, resulting in the identification of 8 distinct stage-specific gene modules (GMs) ( Fig. 2 A) (see Methods). Functional enrichment analysis revealed that GM1, associated with ATP metabolism, and GM2, linked to interferon response, were highly expressed in early pregnancy ( Fig. 2 C). GM3, involved in cell proliferation, was notably expressed in NKp cells. Beyond this, GM4, related to immunity and actin activity, was enriched in secretory stage, which GM5, associated with protein synthesis, was enriched in late pregnancy. GM6, linked to immune migration, show enrichment in non-pregnant stages. Therefore, GM4 and GM6 in NK cells may be involved in cell migration during tissue reconstruction in non-pregnant stages, particularly in NKp and NK3 cells ( Fig. 2 B, and S5 ). Moreover, the significant increase in protein synthesis during late pregnancy suggests enhanced cytokine secretion by NK cells. Within GM1, genes involved in Glycolysis ( TPI1 , GAPDH ), ATP synthesis and catabolism ( ATP5PF , ATP5MF , ATP5MC2 , ATP5MG ), and Aerobic respiration ( UQCRQ , COX7A2 ) were highly expressed in NK cells during early pregnancy ( Fig. S7A ). Our findings thus support previous work which characterized the active glycolytic metabolism amongst NK1 cells within early pregnancy and furthermore, reveal that oxidative phosphorylation activity of NK cells is more active in early pregnancy than in other stages. This was further confirmed by single-cell metabolism analysis ( Fig. S7B ). The origin of NK cells remains unclear [ 7 ], leading us to focus on NKp cell-related GM3 of cell proliferation ( MKI67 , TOP2A ). We found that NKp cells show the strongest proliferative capacity in the non-pregnant stages, which was further confirmed with the highest G2M phase ratio of cell cycle phases ( Fig. 2 D). Our previous work demonstrated that NKp cells have the potential to differentiate into other NK cell subsets in early pregnancy [ 5 ]. To further elucidate this phenomenon, we applied widely used trajectory tools CytoTRACE (Differentiated activities), PAGA (Connectivities) and DPT (Pseudotime) to NK cells of four stages. We found that NK cells of all four stages showed similar differentiation programs, from NKp cells to NK1 to NK2 to NK3 cells ( Fig. 2 E and F), which was also found in each of these four stages ( Fig. S8 ). Furthermore, we found that NKp of non-pregnant stages had a higher differentiation potential compared to NKp of pregnancy ( Fig. 2 G). Meanwhile, the proportion of NKp decreased gradually as the pregnancy progresses ( Fig. 1 E). These results suggest that NKp cells of non-pregnant stages may have a greater proliferative and differentiation capacity than that of pregnancy. Given the possibility that peripheral blood NK (pbNK) cells may act as a source of endometrial NK cells [ 45 ], we assessed the expression levels of GM6 (Leukocyte migration: CCL5 , CXCR4 ) in endometrial NK cells. We discovered that the NK3 subset exhibited high gene expression related to chemotactic properties, such as CXCR4 . For annotated pbNK cells obtained from Vento-Tormo R. et al., GM6 is highly expressed by CD16 + pbNK cells ( Figs. S9A and B ). Moreover, CXCL12 (one ligand of CXCR4 ) was found to have high expression levels in eS cells and fibroblasts during non-pregnant stages ( Fig. S9C ), thereby suggesting that NK3 cells could be recruited by these cells. We also conducted GM analysis on both macrophages and T cells. We discovered that similar to GM2 of NK cells, both macrophages and T cells exhibit the strongest type II interferon (IFN-γ) response in early pregnancy (macrophage: GM5, T cell: GM1) ( Figs. S10 and 11 ). By using CytoSig to predict cytokine signaling activity, we further substantiated this observation ( Fig. 3 A). Further investigation on the genes contained in these interferon-associated GMs, we discovered that during early pregnancy, Mac1 and Mac2 cells express chemokines such as CCL3 and CCL4 ( Fig. 3 B). Additionally, NK cells and T cells display elevated levels of MHC class I molecules, while macrophages display increased levels of MHC class II molecules. These findings suggest that immune cells exhibited an activated state and potential antiviral functions during early pregnancy. Fig. 3 Type II interferon response of major immune cell subsets. A. Box plots predicting the type II interferon signaling activity of major immune cell subsets in four stages. B. Dot plots of the expression of type II interferon response genes of major immune cell subsets in four stages. Statistical significance between different stages was evaluated with the two-sided Wilcoxon rank sum test. In A), the median and interquartile range (IQR), with whiskers extending to 1.5 × IQR, are shown in box plots. **** P  < 0.0001, *** P  < 0.001, ** P  < 0.01, * P  < 0.05, ns ≥ 0.05. Fig. 3 Type II interferon response of major immune cell subsets. A. Box plots predicting the type II interferon signaling activity of major immune cell subsets in four stages. B. Dot plots of the expression of type II interferon response genes of major immune cell subsets in four stages. Statistical significance between different stages was evaluated with the two-sided Wilcoxon rank sum test. In A), the median and interquartile range (IQR), with whiskers extending to 1.5 × IQR, are shown in box plots. **** P  < 0.0001, *** P  < 0.001, ** P  < 0.01, * P  < 0.05, ns ≥ 0.05. In late pregnancy, we observed that NK cells exhibit robust protein synthesis capacity (GM5) and primarily consist of NK2 and NK3 cells, which demonstrate heightened cytokine secretion capabilities ( Fig. 1 , Fig. 2 C, and Fig. S5 ). Mac1 cells demonstrate an increase during late pregnancy and exhibit a strong immune chemotaxis capacity (GM4) ( Figs. S4C and S10 ). T cells, serving as the predominant immune cell type in late pregnancy, express genes associated with interleukin response (GM3) and cell matrix adhesion (GM4) ( Figs. S3B and S11 ). These findings suggest that late pregnancy fosters an actively immunoregulatory microenvironment. Cell-cell interactions between immune cells and non-immune cells at the maternal-fetal interface play a vital role in pregnancy success [ 9 , 46 , 47 ], we next explored the correlation between cell-cell interactions and the maternal events during different stages of reproductive cycle. Initially, we generated a set of stage-specific cytokines, with approximately 93 % of them encoding secreted proteins. Potential stage-specific receptor/ligand pairs were then predicted in different cell subsets during different stages with CellTalkDB (see Methods). Among the four stages, the late pregnancy exhibited the highest diversity of ligand-receptor pairs (n = 20 receptor/ligand pairs) as predicted, followed by early pregnancy (n = 19) and the proliferative stage (n = 12), while the secretory stage had the fewest predicted pairs (n = 5) ( Fig. 4 A, and Figs. S13A and B ). Fig. 4 Cell-cell interactions of the endometrium in four stages. A. Circos plots of the receptor/ligand pairs between different cell subsets in four stages. B. The number of top three signaling pathway-related LR pairs in each stage. C. Dot plots of the expression of ligand-receptor pairs between different cell type pairs. D. Spatial correlation of ligand and receptor expression during four stages. Mean correlation value of each ligand-receptor pair during each stage was calculated with pearson correlation coefficient (PCC) in ST data (2 slides for proliferative stage, 2 slides for secretory stage, 8 slides for early pregnancy, and 3 slides for late pregnancy). Fig. 4 Cell-cell interactions of the endometrium in four stages. A. Circos plots of the receptor/ligand pairs between different cell subsets in four stages. B. The number of top three signaling pathway-related LR pairs in each stage. C. Dot plots of the expression of ligand-receptor pairs between different cell type pairs. D. Spatial correlation of ligand and receptor expression during four stages. Mean correlation value of each ligand-receptor pair during each stage was calculated with pearson correlation coefficient (PCC) in ST data (2 slides for proliferative stage, 2 slides for secretory stage, 8 slides for early pregnancy, and 3 slides for late pregnancy). Next, we utilized CellPhoneDB [ 9 ], a tool developed on the context of endometrial microenvironment, to explore the specific cell-cell interactions between non-immune and immune cell during each stage. These specific cell-cell interactions were then categorized into functionally related signaling pathways, highlighting enriched pathways for cellular interaction pairs at each stage ( Table S8 ). For instance, the HMGB1 and CXCL pathway were prominent in the proliferative stage, the LIFR pathway in the secretory stage, the CCL pathway in early pregnancy, and the TGFb pathway in late pregnancy ( Fig. 4 B). Notably, various chemokine-related interactions during early pregnancy specifically involved immune cells and EVTs, such as CCL3-CCR1 and CCL4-CCR1 ( Fig. 4 C). Previous studies have indicated that immune cells regulate EVT migration through the chemokine-CCR1 pathway in early pregnancy [ 48 , 49 ], thereby validating the reliability of our cell-cell interaction analysis. In the proliferative stage, we found that the ligand HMGB1 , associated with cell proliferation and migration, is primarily expressed in NK2, NK3, and CD4 + T cells, while its receptor SDC1 is expressed in fibroblasts and eS cells ( Fig. 4 C). To validate this stage-specific interaction, we utilized spatial transcriptomic (ST) data across four stages and found that the HMGB1-SDC1 pair exhibits strongest interaction in the proliferative stage ( Fig. 4 D), with co-localization observed between NK cells and eS cells ( Fig. S14A ). Additionally, we observed that regions with HMGB1-SDC1 interaction exhibit stronger proliferation signals, characterized by high expression of MKI67 and TOP2A ( Figs. S14B and C ). In the secretory stage, we identified the LIF - IL6ST pair specifically between NK2 cells and stromal cells ( Fig. 4 C and D). LIF is believed to promote decidualization [ 50 ], and we found that regions enriched for the LIF-IL6ST pair were associated with high expression of decidualization-related genes, such as IGFBP1 ( Figs. S14D–G ). During late pregnancy, we identified stage-specific ligand-receptor pairs including TGFB1 - TGFBR1 , MIF - EGFR , and TIMP3 - CD44 in both single-cell and ST data ( Fig. 4 C, D, and S15 ), which are associated with inhibiting cell migration and growth [ [51] , [52] , [53] ]. Given that various reproductive diseases occur at specific stages, such as recurrent pregnancy loss (RPL) in early pregnancy [ 5 ] and preeclampsia in late pregnancy [ 20 ], we assessed the clinical relevance of stage-specific gene expression profiles in immune cells. We collected risk gene sets of five reproductive diseases from the Harmonizome database [ 36 ], including endometriosis, implantation failure, gestational diabetes (GD), RPL, and preeclampsia ( Table S9 ). With stage-specific genes of each cell type, we obtained stage-specific risk genes and then calculated enrichment scores of each cell type during each stage with Fisher's exact test, which measured the association between the risk genes and the stage-specific genes. Notably, Mac1 in early pregnancy was associated with the high expression of multiple risk genes for various diseases, such as GD, RPL, and preeclampsia ( Fig. 5 A). Fig. 5 Clinical relevance between reproductive diseases and immune cells across reproductive stages. A. The enrichment score of disease risk genes within stage-specific genes of immune cell subsets. Statistical significance calculated by Fisher's exact test shows clinical relevance with stage-specific genes of each cell type. -log10 ( P value) was defined as the enrichment score. * means P  < 0.05. B. The Enrichment score of immune cell subsets in different reproductive cycle stages. C. The number of stage-specific risk genes in immune cell subsets during different reproductive cycle stages. Fig. 5 Clinical relevance between reproductive diseases and immune cells across reproductive stages. A. The enrichment score of disease risk genes within stage-specific genes of immune cell subsets. Statistical significance calculated by Fisher's exact test shows clinical relevance with stage-specific genes of each cell type. -log10 ( P value) was defined as the enrichment score. * means P  < 0.05. B. The Enrichment score of immune cell subsets in different reproductive cycle stages. C. The number of stage-specific risk genes in immune cell subsets during different reproductive cycle stages. We found that in early pregnancy, all presented immune cell types showed a high correlation with the RPL risk gene set ( Fig. 5 B), which is consistent with the occurrence of RPL typically in early pregnancy. To investigate the relationship between stage-specific risk genes and disease onset, we utilized additional single-cell data from RPL patients collected during early pregnancy from a previous study [ 5 ]. Differential expression analysis revealed that majority of the risk DEGs were stage-specific ( Fig. S16A ). Compared to non-stage-specific risk genes, all stage-specific risk genes were significantly upregulated in RPL patients ( Fig. S16B ), indicating a potential association between the upregulation of stage-specific risk genes and RPL. Additionally, we identified stage-specific risk genes associated with endometriosis and implantation failure during the menstrual cycle, as well as with GD and preeclampsia during pregnancy ( Fig. 5 A and B). Further investigation into stage-specific risk genes revealed genes related to interferon response and stage-specific interactions ( Fig. 5 C). We also assessed the clinical relevance of stage-specific gene expression profiles in non-immune cells. Notably, genes that are considered pivotal for the appearance of endometriosis, for instance, CDKN1B , HOXA10 , and WNT4 [ [54] , [55] , [56] ], are found to be enriched in the secretory stage ( Fig. S17A ). To further facilitate researchers in exploring potential disease predictive biomarkers and therapeutic targets, we have compiled lists of disease-associated stage-specific risk genes for all cell types across stages ( Table S10 ).

For the scRNA-seq data we collected, 4 stages of the reproductive cycle were examined: proliferative (n = 4), secretory (n = 11), early pregnancy (n = 28), and late pregnancy (n = 9). For sample sources, 1 of proliferative donors and 1 of secretory donors are deceased organ donors (within 1 h of circulatory arrest) from Garcia-Alonso, L. et al. ( Nat. Genet. , 2021), and the other 13 of non-pregnant donors are live donors; Decidual tissue of Vento-Tormo et al. ( Nature , 2018) was obtained from elective terminations of normal pregnancies between 6 and 14 gestational weeks. In the other 3 datasets of early pregnancy, Guo et al. ( Cell Discov. , 2021), Wang et al. ( Genomics Proteomics Bioinformatics , 2021), and Chen et al. ( Front. Immunol. , 2021), decidual tissues were obtained from elective terminations of apparently normal pregnancies and CD45 + cells were sorted for further sequencing analysis. Tissues of late pregnancy were obtained immediately after delivery. scRNA-seq data used in our study was obtained with 10x Genomics libraries. More information of cohorts could be found in Table S1 . For the raw sequencing SRA files provided by the articles, fastq-dump (v2.8.0) was used here to download and convert to FASTQ.GZ files, and then Cell Ranger (v4.0.0) and the GRCh38 human reference genome provided by it were used for sequence comparison. For data quality control, in this study, six scRNA-seq datasets were first merged and retained 36,601 standard genes using the Merge function of Seurat (v4.0.5) [ 21 ], and the data were subsequently partitioned into smaller datasets using the difference in the samples to which they belonged. After preprocessing with Seurat's NormalizeData and ScaleData functions, this study used DoubletFinder (v2.0.3) [ 22 ] to screen out double cells with a default setting of 7.5 %. Subsequently, cells with detected gene counts between 500 and 6000 and with less than 25 % mitochondrial gene expression were retained in this study. In addition, we retained genes expressed in at least 10 cells and simultaneously removed mitochondrial genes. Due to severe batch effects in direct follow-up analysis, different integration approaches were applied here, including Seurat, Harmony (v0.1.0) [ 23 ] and scVI (v0.13.0) [ 24 ]. We first run these packages based on their default parameters. We observed that scVI exhibited poor batch effect removal in our data, while Harmony showed relatively better performance in batch effect correction ( Fig. S1C ). However, we noticed an abnormal connection between immune and non-immune cells in the UMAP space, such as uterine smooth muscle cells (uSMCs) and NK cells ( Fig. S1C ), which finally prompted us to choose Seurat. The full-step integration process included splitting the dataset (SplitObject), normalizing the data (NormalizeData), obtaining highly expressed genes (FindVariableFeatures), obtaining highly expressed genes with a high number of replicates based on the splitting results (SelectIntegrationFeatures), determining integration anchors (FindIntegrationAnchors), and integrating (IntegrateData), while adjusting the split.by parameter of the SplitObject() function in the Seurat standard integration process, the nfeatures parameter in the FindVariableFeatures and SelectIntegrationFeatures functions, the anchor.features parameter in the FindIntegrationAnchors function, and the IntegrateData function's sample.tree parameter in the FindIntegrationAnchors function and sample.tree parameter in the IntegrateData function to further optimize the integration results. The integration also included Seurat's standard process for scRNA-seq downstream data processing and analysis, and this study used the integrated data matrix for normalization (ScaleData), PCA principal component analysis (RunPCA), UMAP downscaling (RunUMAP), calculation of neighborhoods (FindNeighbors), and cell clustering (FindClusters, resolution = 1.0). In this study, the cell types were broadly classified mainly based on the cell type characteristic genes in the Extended Data provided in the article [ 9 ]. Here we did not regress out with Cell Cycle genes in integration step, which leaded to identify our proliferative cell types. Additionally, using the Jaccard index (Ratio of the intersection size of two sets to the size of the concatenated set), our study calculated the similarity between the natural clusters and the original annotations of the datasets to confirm the cluster-cell type correspondence. When batch effect appeared in the further detailed annotation, we still used Seurat to integrate before dimensionality reduction. For clusters hard to be annotated with existing markers, we used the FindMarkers of Seurat to obtain the DEGs of the cluster. After cell annotation was completed, we removed some cells that included (a) NK cells with high expression of HSP genes, which may be due to experimental manipulation causing cellular stress; (b) clusters of cells expressing both T cells and Mac signature genes; and (c) clusters of T cells with high expression of antibody-related proteins. Pro- and anti-inflammatory gene signatures were obtained from Leader A.M. et al. [ 25 ] ( Table S4 ). Pro- and anti-inflammatory score were then calculated for all macrophages using the AddModuleScore function in Seurat. The M1-like and M2-like signatures consistently aligned with phenotypes of endometrial macrophage subsets across all four stages ( Fig. S4A ). For further comparison and visualization, we normalized the pro- and anti-inflammatory scores of these cells to a 0–1 scale. Gene module (GM) analysis is a common strategy for studying transcriptomic heterogeneity from single-cell data [ 26 , 27 ]. We generated gene modules based on cNMF and hierarchical clustering, referring to these studies. The analysis steps are listed as follows. (1) Input of NMF analysis. To obtain NMF programs associated with celltypes and stages, we identified differentially expressed genes (DEGs) of each cell subset across stages and DEGs of each stage across cell subsets, excluding mitochondrial (MT) and heat shock protein (HSP) genes. (2) Gene program analysis. We then input the expression matrix containing the aforementioned DEGs and all cells into cNMF (v1.3.2) analysis, setting the K parameter to obtain NMF gene programs according to the number of DEGs (10–20 for NK cells, 7–12 for macrophages and T cells). To minimize overlap between NMF gene programs, we selected the top n genes based on NMF coefficients (top 50 for NK cells, top 20 for macrophages and T cells). We identified 20 gene programs in NK cells, 12 gene programs in macrophages and 10 gene programs in T cells. (3) Gene module analysis. To group gene programs with similar expression patterns, we calculated the average expression value of each gene program in individual cells (NK cells or Macrophages or T cells) and computed the Pearson correlation coefficients (PCC) between different gene programs. Based on PCC, we performed hierarchical clustering to obtain GMs. Initially, we identified 10 GMs in NK cells, 7 GMs in macrophages, and 6 GMs in T cells. After filtering out GMs that did not show significantly differences between cell types or stages (Fold change ≥1.1, P  < 0.05), we finally obtained 8 GMs in NK cells, 5 GMs in macrophages and 4 GMs in T cells ( Tables S5, S6, and S7 ). Each GM contains all genes of each NMF gene program cluster. Further functional enrichment analysis of these GMs was performed using ClusterProfiler (v3.14.3) [ 28 ]. Input of NMF analysis. To obtain NMF programs associated with celltypes and stages, we identified differentially expressed genes (DEGs) of each cell subset across stages and DEGs of each stage across cell subsets, excluding mitochondrial (MT) and heat shock protein (HSP) genes. Gene program analysis. We then input the expression matrix containing the aforementioned DEGs and all cells into cNMF (v1.3.2) analysis, setting the K parameter to obtain NMF gene programs according to the number of DEGs (10–20 for NK cells, 7–12 for macrophages and T cells). To minimize overlap between NMF gene programs, we selected the top n genes based on NMF coefficients (top 50 for NK cells, top 20 for macrophages and T cells). We identified 20 gene programs in NK cells, 12 gene programs in macrophages and 10 gene programs in T cells. Gene module analysis. To group gene programs with similar expression patterns, we calculated the average expression value of each gene program in individual cells (NK cells or Macrophages or T cells) and computed the Pearson correlation coefficients (PCC) between different gene programs. Based on PCC, we performed hierarchical clustering to obtain GMs. Initially, we identified 10 GMs in NK cells, 7 GMs in macrophages, and 6 GMs in T cells. After filtering out GMs that did not show significantly differences between cell types or stages (Fold change ≥1.1, P  < 0.05), we finally obtained 8 GMs in NK cells, 5 GMs in macrophages and 4 GMs in T cells ( Tables S5, S6, and S7 ). Each GM contains all genes of each NMF gene program cluster. Further functional enrichment analysis of these GMs was performed using ClusterProfiler (v3.14.3) [ 28 ]. Respiratory metabolic activity of proliferating NK cells was calculated using scMetabolism [ 29 ]. Activity of Glycolysis/Gluconeogenesis, Citrate cycle (TCA cycle), and Oxidative phosphorylation could be calculated with sc.metabolism.Seurat() function. We use recommended parameters, method for ‘AUCell’, metabolism.type for ‘KEGG’. The differentiation potential of NK cells at each of the four stages was calculated using CytoTRACE (v0.3.3) [ 30 ], respectively. CytoTRACE outperformed previous methods and nearly 19,000 annotated gene sets for resolving 52 experimentally determined developmental trajectories. We run CytoTRACE with raw count matrix. The graphs of predicted order were created by plotCytoTRACE function with UMAP coordinates. In contrast to many other pseudotime analysis tools, CytoTRACE does not require an artificial definition of the differentiation starting point. Specially, we calculated CytoTRACE score of NK cells in all stages in Fig. 2 E and F, each stage in Fig. S8 , and only NKp cells of all stages in Fig. 2 G. Then we used PAGA [ 31 ] and DPT [ 32 ] as another trajectory evidence with Scanpy (v1.9.1) under Python (v3.7.11). PAGA graph was calculated based on the seurat integrated pca matrix produced in Fig. 1 . Fig. 1 Single-cell transcriptional landscape of four stages of the reproductive cycle. A. A schematic outline depicting the workflow for data collection from published literature and subsequent integrated analysis. The number of samples and the number of single-cell transcriptomes collected in each stage (proliferative, secretory, early pregnancy, and late pregnancy stages) are indicated. B. Uniform manifold approximation and projection (UMAP) embeddings of integrated single-cell transcriptomes of four stages of samples (dataset batches shown in right bottom). Cells are colored by cell subsets, and dashed circles indicate the major cell types. NK, natural killer cells; ILC, innate lymphocyte cells; p/pro, proliferative; Mac, macrophages; Mac3, maternal macrophages; HB, Hofbauer cells; DC, dendritic cells; T, T cells; Teff, effector T cells; Tex, exhausted T cells; MAIT, Mucosal-associated invariant T cells; dS, decidual stromal cells; eS, endometrial stromal cells; PV, perivascular cells; Epi, epithelial cells; Endo, endothelial cells; Fibro, fibroblasts; uSMC, uterine smooth muscle cells; SCT, syncytiotrophoblast; VCT, villous cytotrophoblast; EVT, extravillous trophoblast. C. UMAP embeddings of different stages illustrating no obvious batch effect in this integrated atlas. D. Violin plots of canonical markers (columns) for major cell types (rows). E. Bar plots of the cell subset composition of NK cells in four stages. Fig. 1 Fig. 2 Transcriptome changes in NK cells in four stages. A. Pearson correlation coefficients (PCC) of the average gene expression of NMF gene programs in all NK cells. Each row or column is a NMF gene program. Example genes for each gene module are shown on the right side of the heatmap. B. The average gene expression of GM1, GM3, and GM6 in different NK cell subsets in four stages. C. Gene Ontology (GO) analysis of gene modules in NK cells. D. The G2M phase ratio of cell cycle in NK cells during four stages. E. Trajectory analysis of NK cells from all four stages based on CytoTRACE. 'Less diff.' indicates greater differentiation potential. F. Box plots of NK cell differentiation potential calculated by CytoTRACE of all four stages. G. Box plots of NKp differentiation potential calculated by CytoTRACE in four stages. Statistical significance between different stages was evaluated with the two-sided Wilcoxon rank sum test. In B), E) and G), the median and interquartile range (IQR), with whiskers extending to 1.5 × IQR, are shown in the plots. **** P  < 0.0001, *** P  < 0.001, ** P  < 0.01, * P  < 0.05, ns ≥ 0.05. Fig. 2 Single-cell transcriptional landscape of four stages of the reproductive cycle. A. A schematic outline depicting the workflow for data collection from published literature and subsequent integrated analysis. The number of samples and the number of single-cell transcriptomes collected in each stage (proliferative, secretory, early pregnancy, and late pregnancy stages) are indicated. B. Uniform manifold approximation and projection (UMAP) embeddings of integrated single-cell transcriptomes of four stages of samples (dataset batches shown in right bottom). Cells are colored by cell subsets, and dashed circles indicate the major cell types. NK, natural killer cells; ILC, innate lymphocyte cells; p/pro, proliferative; Mac, macrophages; Mac3, maternal macrophages; HB, Hofbauer cells; DC, dendritic cells; T, T cells; Teff, effector T cells; Tex, exhausted T cells; MAIT, Mucosal-associated invariant T cells; dS, decidual stromal cells; eS, endometrial stromal cells; PV, perivascular cells; Epi, epithelial cells; Endo, endothelial cells; Fibro, fibroblasts; uSMC, uterine smooth muscle cells; SCT, syncytiotrophoblast; VCT, villous cytotrophoblast; EVT, extravillous trophoblast. C. UMAP embeddings of different stages illustrating no obvious batch effect in this integrated atlas. D. Violin plots of canonical markers (columns) for major cell types (rows). E. Bar plots of the cell subset composition of NK cells in four stages. Transcriptome changes in NK cells in four stages. A. Pearson correlation coefficients (PCC) of the average gene expression of NMF gene programs in all NK cells. Each row or column is a NMF gene program. Example genes for each gene module are shown on the right side of the heatmap. B. The average gene expression of GM1, GM3, and GM6 in different NK cell subsets in four stages. C. Gene Ontology (GO) analysis of gene modules in NK cells. D. The G2M phase ratio of cell cycle in NK cells during four stages. E. Trajectory analysis of NK cells from all four stages based on CytoTRACE. 'Less diff.' indicates greater differentiation potential. F. Box plots of NK cell differentiation potential calculated by CytoTRACE of all four stages. G. Box plots of NKp differentiation potential calculated by CytoTRACE in four stages. Statistical significance between different stages was evaluated with the two-sided Wilcoxon rank sum test. In B), E) and G), the median and interquartile range (IQR), with whiskers extending to 1.5 × IQR, are shown in the plots. **** P  < 0.0001, *** P  < 0.001, ** P  < 0.01, * P  < 0.05, ns ≥ 0.05. In our study, CytoSig [ 33 ] was used for the prediction of multiple cytokine signaling activities. The single-cell raw matrices of NK cells, macrophages and T cells were first normalized by the required log2(TPM/10 + 1), and then the processed matrices were output as.txt files, which were converted to.txt.gz files using the gzip command at the command line. In this study, we used.txt.gz files as the standard inputs to perform the signal activity prediction of 51 cytokines in the CytoSig database in the command line, and finally, we selected the.Zscore (regression coefficient/standard error) files as the result for subsequent usage. We first obtained a well-known gene list of cytokines [ 34 ] and then performed a preliminary screening and retained them if their average expression after normalization in single or multiple stages was greater than 0.5. We obtained stage-specific cytokines with significantly higher expression in one stage than in other stages and a fold change more than or equal to 1.05. We then performed hierarchical cluster on expression of cytokine ligand genes using ComplexHeatmap (v2.11.1) and obtained 4 clusters. We used the CellTalkDB [ 35 ] to obtain the receptors corresponding to these cytokine ligand genes. One gene was believed to be expressed in one cell type when it was expressed in at least 50 % of the cell type. Then, the potential receptor/ligand and cell type pairs were defined when cell type A expressed ligand gene X and cell type B expressing the corresponding receptor gene Y in the same stage. Finally, we showed all potential cell-cell interactions we predicted here with igraph (v1.3.1). We further predicted interaction pairs between immune and non-immune cells in another way by using CellPhoneDB (v3.0.0) 9 with a combined interaction pairs database which includes CellPhoneDB database and CellTalkDB database [ 35 ]. To better investigate cell-cell interactions at different stages, we collected spatial transcriptomics data corresponding to the stages of single-cell data. Cohort information of the spatial transcriptomic data can be found in Table S2 . Pathological annotation information for the slides was obtained from the original research. The visualization of spatial gene expression was based on log-normalized data. The spatial correlation of ligand and receptor genes was calculated in a smoothed gene expression matrix using the PCC. In the smoothed gene expression matrix, each gene expression value for each index cell was the mean gene expression taken from the cells in the respective index cell's neighborhood (a circle of cells adjacent to the index cell). The decidualization score was calculated using the decidualization gene set (GO: 0046697) with the AddModuleScore function. In analysis part of disease risk genes correlation with stage-specific genes, we used Fisher's exact test. Risk gene sets were obtained from Harmonizome database [ 36 ]. Compared to other tests, Fisher's exact test is a more appropriate exact test method for our small sample size (number of genes) [ 37 ]. In addition, Fisher's exact test is usually used to calculate the relative enrichment score [ 38 ]. The function for Fisher's exact test is fisher.test in R. Statistical significance in this study was primarily determined using the two-sided Wilcoxon rank sum test. The Benjamini-Hochberg method was applied to adjust for multiple comparisons. The Chi-Squared test and Fisher's exact test were utilized for disease association analysis.

To elucidate the changes that endometrial immune cells undergo throughout the reproductive cycle, we constructed an integrated single-cell transcriptomic atlas spanning various stages. This revealed the composition of immune cell types in the endometrium at each stage, such as a significant increase in NK1 cells during early pregnancy and a marked rise in T cells during late pregnancy. Utilizing NMF analysis, we identified stage-specific immune cell gene modules, including proliferation and chemotaxis modules during the menstrual cycle, glucose metabolism and interferon response modules in early pregnancy, and protein synthesis modules in late pregnancy. Additionally, we uncovered a potential NK cell differentiation trajectory that remains stable across the four stages. When exploring the relationship between stage-specific changes and immune-non-immune cell-cell interactions, we discovered interactions linked to the promotion and inhibition of cell proliferation, migration, and decidualization, which we further validated using spatial transcriptomic data. Finally, we collected risk genes associated with reproductive diseases and identified stage-specific risk genes ( Fig. 6 ). These findings may provide a good reference and potential biomarkers for future research on endometrial diseases.

Immune cell infiltration of the endometrium has been proposed as an essential step for promoting a successful pregnancy and immune homeostasis [ 4 ]. The large majority of published studies have used scRNA-seq to examine endometrial cell heterogeneity exclusively in pregnancy or nonpregnancy stages. In this study, we compiled scRNA-seq datasets from several cohorts obtained at four distinct stages of the reproductive cycle, including the proliferative, secretory, early pregnancy, and late pregnancy stages. Through advanced bioinformatic analysis, we characterized the immune landscape at each stage to identify changes in the composition, function, and cell-cell interactions of endometrial cell populations ( Fig. 6 ). Analysis of proportions revealed changes in the proportions of NK1 and NK2 subsets during four stages of reproductive cycle differ from those observed in NK3 cells. Identification of stage-specific GMs revealed that these immune cells have a strong type II interferon response in early pregnancy. Additionally, we found cell-cell interactions related to proliferation promotion in the proliferative stage and tissue remodeling in early pregnancy. Our integrative analysis of immune cell subsets during the reproductive cycle provided insights into the advancement of reproductive processes. Fig. 6 Graph abstract of integrated single-cell and spatial transcriptomics data analysis in healthy endometrium across reproductive cycle. ▲ represents a significant increase in the proportion of this cell type at the indicated stage, while * indicates a significant association of this cell type with reproductive disease risk at the indicated stage. Created with BioGDP.com . Fig. 6 Graph abstract of integrated single-cell and spatial transcriptomics data analysis in healthy endometrium across reproductive cycle. ▲ represents a significant increase in the proportion of this cell type at the indicated stage, while * indicates a significant association of this cell type with reproductive disease risk at the indicated stage. Created with BioGDP.com . Decidual NK cells have been thought to be derived from CD34 + hematopoietic stem cells, mature from immature endometrial NK cells, or to migrate from peripheral blood NK cells [ 7 , 45 , 46 , 57 , 58 ]. We found that NKp cells have a strong proliferative capacity and differentiation potential during non-pregnant stages, and this signal diminishes during pregnancy. The uterus maintains a low oxygen concentration in the early stages of implantation, while appearing to rise in early and mid-pregnancy [ 59 ]. We speculate that low oxygen concentrations during non-pregnant stages favor glycolytic metabolism and proliferation of NKp [ 60 ], whereas oxygen concentrations during early pregnancy favor a strong oxidative phosphorylation metabolism and functions of NK cells. We also observed that the relative proportions of NK1 and NK2 cells showed an increase in early pregnancy and a decrease in late pregnancy. Both existing studies [ 5 , 9 , 14 ] and our results show that NK1 and NK2 cells have similar transcriptome profiles compared with NKp cells. In contrast, NK3 cells have the least transcriptome similarity to NKp cells, while they exhibit a chemokine pattern similar to that of pbNK cells. Several studies of decidua in mid- and late pregnancy failed to separate uterine leukocytes from maternal or fetal blood leukocytes [ 61 ], indicating that uterine leukocytes are similar to peripheral blood leukocytes at this time, which may be associated with an increase in the proportion of NK3 in late pregnancy. Thus, we speculate that NK cells may derive from NKp and pbNK cells, while NK1/2 and NK3 may arise from different precursor cells. During the process of placentation, EVT migrates into the maternal uterus to remodel spiral arteries [ 9 ]. We identified the cytokines CCL3 , CCL4 , CCL5 , ICAM1 , and SPP1 during early pregnancy, which may help to promote EVT invasion. SPP1 was previously reported to be produced by NK cells and macrophages [ 9 , 61 ], and its expression was higher in macrophages than in NK cells. Excessive invasion of EVT is harmful, and the rate of EVT invasion decreases during pregnancy [ 62 ]. TGFB1 and TIMP3 were identified during late pregnancy in our results, which are both reported to regulate EVT invasion [ 62 ]. Immune cells (NK cells, macrophages, T cells, and DCs) and non-immune cells (fibroblasts, endothelial cells, stromal cells, and other trophoblast cells) appear to be involved in this process. Mast cells and neutrophils may also regulate EVT invasion [ 62 ]. By integrating our data with risk gene sets associated with common reproductive diseases, we elucidated the clinical relevance of stage-enriched gene expression profiles. For instance, risk genes for RPL preferentially expressed in early pregnancy were significantly upregulated in RPL patients. In a recent study [ 20 ], maternal blood samples were utilized to identify and verify cfRNA transcriptomic changes that are linked to preeclampsia. Genes that can differentiate between patients with preeclampsia and healthy individuals in early pregnancy were identified. By examining intersections of stage-specific risk genes and cfRNA genes, two candidate risk genes ( CLU and TGFBR1 ) were identified, especially TGFBR1, which was highly expressed in both immune cells and non-immune cells during the secretory stage ( Figs. S17B and C ). Thus, our integrated analysis underscores the potential clinical significance of profiling stage-specific gene sets. Our work has some limitations and room for improvement. Our study included a relatively small size in certain stages, such as the proliferative stage and late pregnancy, necessitating further validation with larger cohorts to enhance the credibility of our findings. Since the late pregnancy samples were obtained immediately after delivery, we cannot rule out the possibility that delivery may have introduced unknown changes to the transcriptomic profiles of endometrial cell types. Although we visualized the co-localization of stage-specific pairs of ligands and receptors using spatial transcriptomic data, further validation at the protein level (e.g., multiplex immunofluorescence) is still required. In summary, our new findings in this study provide a reference for further investigation of the multifaceted changes in immune cells in different stages of the healthy female reproductive cycle and provide clues for the prediction and treatment of pregnancy-related diseases.

The reproductive cycle can be broadly classified into the menstrual cycle, which encompasses the proliferative and secretory stages [ 1 ], and the stages of pregnancy, including the early, mid-, and late pregnancy [ 2 ]. Many significant occurrences in the reproductive process are closely linked with the endometrium, including embryo implantation, pregnancy, and labor. Such instances are accompanied by a dynamic process involving shedding, regeneration, and differentiation of the endometrial tissue [ 3 ]. A comprehensive characterization of the endometrium in healthy individuals throughout the reproductive cycle can facilitate the understanding of the normal transitions and variations in the endometrial microenvironment during different stages of the reproductive cycle, and serve as a foundation for exploring pathological changes in reproductive processes. The immune microenvironment of the endometrium significantly influences the success of pregnancy [ 4 , 5 ]. Previous research indicates a notable increase in endometrial immune cell proportion from 8.2 % during the proliferative stage to 31.7 % during early pregnancy [ 6 ], among which the most predominant are NK cells, macrophages, and T cells [ 7 ]. Single-cell RNA sequencing (scRNA-seq) has been increasingly used to characterize the endometrial microenvironment at various stages of the reproductive cycle [ 8 ], including identifying new subsets of cell types and studying cell development [ 1 , [9] , [10] , [11] , [12] ]. During early pregnancy, decidual NK cells were observed to comprise three distinct subsets, including dNK1-3 [ 9 ]. dNK1 cells are noted for expressing high levels of killer immunoglobulin-like receptors (KIRs) capable of binding to HLA-C. During early pregnancy, both dNK1 and dNK2 cells are believed to interact with EVTs, with the presence of LILRB1 and NKG2A which encode proteins binding to HLA-G and HLA-E. Meanwhile, dNK3 and dNK2 cells are observed to produce more cytokines than dNK1 cells [ 13 ]. However, a comprehensive understanding of the dynamic changes in the endometrial immune environment spanning the reproductive cycle is still lacking. Recent studies on single-cell transcriptome and peripheral blood transcriptome have provided valuable insights into uterine endometrial diseases such as recurrent pregnancy loss (RPL) [ 5 , [14] , [15] , [16] ], endometriosis [ 17 , 18 ], and preeclampsia [ 19 , 20 ], proposing potential new therapeutic targets and diagnostic approaches. These studies have uncovered intricate cellular and molecular landscapes underlying reproductive diseases and physiological processes, elucidating dysregulated immune responses, altered cell-cell interactions, and disrupted tissue homeostasis. Furthermore, they have identified potential therapeutic targets for modulating immune responses, apoptosis, inflammation, and tissue regeneration. For instance, Li Qian and colleagues found that various genes regulated by NK cells in EVTs are dysregulated in preeclampsia. They also discovered that these genes may indicate increased risk in other pregnancy diseases [ 19 ]. Preeclampsia typically occurs after 20 weeks of gestation. However, Mira N. Moufarrej and colleagues identified peripheral blood transcriptome genes that can predict the risk of preeclampsia in early pregnancy (5–12 weeks) [ 20 ]. Endometrial diseases are often associated with dysregulation of gene expression. Concurrently, transcriptomic characteristics within the endometrium undergo significant changes across different stages of the reproductive cycle. However, the relationship between these alterations and the risk of endometrial diseases remains unclear. In this study, we conducted an advanced bioinformatic analysis of large-scale public scRNA-seq datasets obtained from endometrial cells of healthy women across four stages of the reproductive cycle. We identified the stage-specific transcriptomic characteristics of major endometrial immune cells (including NK cells, macrophages, and T cells) during different stages. Subsequently, we explored the relationship between cell-cell interactions based on stage-specific differentially expressed genes and maternal events across different stages. We found the strongest immune-stromal interactions in late pregnancy which may be associated with the shaping of the microenvironment leading to labor. Finally, we identified some stage-specific risk genes by using the above stage-specific genes and the risk genes of common reproductive diseases. This study provides a relatively comprehensive perspective of the dynamic landscape in the endometrial immune microenvironment throughout the reproductive cycle, which can facilitate studies investigating pathological changes to improve the diagnosis and treatment of endometrium-related ailments.

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:Chuang Guo reports financial support was provided by 10.13039/501100001809 National Natural Science Foundation of China grants. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.