All participants, aged between 30 and 40 years, were recruited from the Endocrinology Clinic of Changhai Hospital from June 2021 to October 2021. Patients with POI (n=3) and healthy control (HC, n=3) women with normal ovarian reserve were recruited. According to the ESHRE and Chinese guidelines, the eligibility criteria for POI included secondary amenorrhea for at least 4 months and serum basal FSH > 25 IU/L (on two occasions with more than 1 month interval) before 40 years of age. Women with chromosomal abnormalities, known gene mutations, history of ovarian surgery, radiotherapy or chemotherapy, recurrent spontaneous abortion, endometriosis, autoimmune diseases, and infection within three months before enrollment were excluded. This study was approved by the Institutional Review Board of Changhai Hospital. All participants signed the written informed consent forms.

The preparation and cell suspension of PBMC were performed as previously described (24). Cell capture, cDNA amplification, and library preparation were performed using Chromium Single Cell 3’ Library and Gel Bead Kit v2 (10x Genomics). Then, the libraries were sequenced on the Illumina Nova Seq 6000 platform with a 150-bp paired-end mode.

The Cell Ranger software pipeline (version 5.0.0) provided by 10×Genomics was used to preprocess the raw sequencing data across samples, producing a matrix of gene counts versus cells. The gene-barcode matrix of unique molecular identifier (UMI) counts was then analyzed with Seurat (version 4.0.0) (25) for quality control, normalization, dimensional reduction, batch effect removal, clustering, and visualization. To remove low-quality cells and multiple captures, the following criteria were applied: 1- Cells with UMI/gene numbers more than the limit of median value +/- 2 fold of median absolute deviations were filtered out, assuming a Guassian distribution of each cells’ UMI/gene numbers. 2- We discarded low-quality cells where >10% of the counts belonged to mitochondrial genes. 3- We removed the potential duplicates identified by the Doublet Finder package (version 2.3.0) (25). The count matrix was log-normalized, and the top 3000 most variable genes were identified for dimensional reduction. Most of the variable genes were selected using the “Find variable Genes” function. Principal component analysis (PCA) was performed to reduce the dimensionality with “Run PCA” function in Seurat. Graph-based clustering was performed to cluster cells according to their gene expression profile using the “Find Clusters” function. The “Find All Markers” function was used to identify marker genes of each cluster. The R package Single R (26), a novel computational method for unbiased cell type recognition of scRNA-seq was used to independently infer the cell of origin for each single cell and identify cell types based on the reference transcriptomic dataset ‘Human Primary Cell Atlas’ (27).

A pseudo-bulk expression profile was generated by summing the UMI counts of all cells with the same combination of cell type and sample. We then performed differential expression analysis for each cell type on the sample level using the R package DESeq2. The R package cluster Profiler (v.3.16.0) was used for function over-representation analyses of DEGs with FDR < 0.05.

DEGs were identified using the “Find Markers” function in Seurat. P value < 0.05 and |log₂ fold change| > 2 was set as the threshold for significantly differential expression. GO enrichment and KEGG pathway enrichment analyses of DEGs were performed using the R package cluster Profiler (28) and based on the hypergeometric distribution.

Cell Chat was utilized to infer and analyze global cell-cell communications from scRNA-seq data. For cells in HC or POF, gene expression data for corresponding cells were utilized as input, and cell-cell communication probability was calculated between interactions of ligands.

The number of captured cells in the 6 samples was 13367 for P1, 8601 for P2, 9351 for P3, 9372 for N1, 7347 for N2, and 8648 for N3, which respectively decreased to 12615, 8081, 8904, 8682, 6442, and 7966 after quality control ( Figure 1A ). Data from multiple sequencing were merged using a cell ranger pipeline, and all cells were sorted according to the number of detected genes. Approximately, 1710 genes were identified in these 6 samples, and more genes were identified in N3 (a normal control) than in other samples ( Figure 1B ).

Pairwise correlation values of inter-groups ranged from 0.94 to 0.99 ( Figure 1C ), revealing high consistency in cell type distribution among samples and inter-group heterogeneity.

We then performed unsupervised clustering to separate the cells into 22 distinct clusters, which were visualized by Uniform Manifold Approximation and Projection for Dimension Reduction (UMAP) ( Figure 2A ). Based on the previously validated canonical marker genes, these clusters were assigned to 9 cell lineages, including T cell subsets (naive T cells, CD4+ T cells, naive CD8+ T cells, effector CD4+ T cells, and effector CD8+ T cells), B cell subsets (naive B cells and plasma B cells), NK cells, monocyte cells, and DCs. The distribution of these cells in HC and POI samples is shown in Figure 2A . T cells constituted the majority of all cells. B cells, NK cells, monocytes, and other cell types constituted small percentages of cells, It is consistent with the characteristics of the composition of human peripheral blood cells. The heatmap confirmed the distinct gene expression patterns among the clusters ( Figure 2B ).

As shown in Figure 2C , a significant decrease in the abundance of CD14+monocytes (p-value=0.0013, two-sided t-test) and proliferating T cells (p-value=0.047, two-sided t-test), and NK cells (p-value=0.035, two-sided t-test) was identified in POI compared with HC. In contrast, the abundance of plasma B cells showed a significant increase (p-value=0.0057, two-sided t-test) in POI. No significant difference was observed in the abundance of CD4+ T cells (p-value=0.066, two-sided t-test), while a significant reduction was observed in the abundance of CD8+ T cells in patients with POI (p-value=0.02, two-sided t-test). The CD4/CD8 ratio was significantly higher in patients with POI. Moreover, three atypical cell types, including CD8+ NKT-like cells, MAIT cells, and γδ T cells, were decreased in POI patients; however, results were not statistically significant. Notably, the abundance of Treg cells was increased in patients with POI.

In total, 24 downregulated genes and 7 upregulated genes between HC and POI were found in PBMCs ( Figure 3A ). IGHM and LYZ were respectively the most significantly upregulated and downregulated genes among all cell clusters of POI patients. Among the upregulated genes, IGHM, IGKC, IGHD, and IGLC2 were predominantly expressed in B cells and DCs of POI patients corresponding to clusters 5, 8, 17, 18, and 22 ( Figure 3B ). IGKC, IFITM1, CD69, and JUND were upregulated in almost all cell types of POI patients compared with HC ( Figure 3B ). In contrast, downregulated genes, such as LYZ, GNLY, VCAN, and S100A9, were mostly related to NK cells, CD8+ T cells, and monocytes ( Figure S1 ).

To further investigate the transcriptomic differences between HC and POI, pathway enrichment analysis was performed using the DEGs. Upregulated genes were enriched in B cell receptor signaling pathway, IL-17 signaling pathway, parathyroid hormone synthesis, secretion and action, and osteoclast differentiation pathways. Downregulated genes were mainly enriched in NK cell-mediated cytotoxicity, cytolysis, and cell-killing functions ( Figure 3C ).

The results of DEG analysis for all clusters showed a decrease in monocytes, especially CD14+ monocytes, suggesting that monocytes may be affected in POI. Therefore, we re-analyzed the cluster of monocytes and divided cells into 5 subsets ( Figure 4A ), including CMono0, CMono1, CMono2, CMono3, and CMono4 ( Figure 4B ). Surprisingly, we found that the proportion of Mono0-subset showed a significant decrease while the proportion of other subsets remained stable.

Each subset exhibited a distinct gene expression pattern. CMono0, expressed high levels of classical monocyte marker genes, including CTSS, FCN1, VCAN, KTD12, CYP1B1, and IFI30, indicating that this subset was predominantly classical monocytes. CMono2 expressed high levels of neutrophil-related genes, including CXCL8, IFIT3, IFIT2, IFIT1, NAMPT, IFITM2, G0S2, FCGR3B, CXCR2, and CMTM2. CMono3 expressed high levels of DC marker genes, while CMono4 expressed high levels of basophil marker genes ( Figure 4C ). All cells were phagocytes and there was a group of innate myeloid immune cells.

Next, 47 DEGs from all subsets of monocyte were identified between HC and POI. The highly up-regulated DEGs, such as IGLC2, G0S2, FCGR3B, IFITM1, LTB, IGKC, and IFITM2, were predominantly found in classical monocytes in POI ( Figure 4D ). These DEGs were significantly enriched in the infection and T cell differentiation pathways ( Figures 4E, F ). Interestingly, ISG15, a gene from these pathways, was highly upregulated in the CMono1-subset.

Interestingly, 8 DEGs from the classical monocytes, including 2 upregulated genes (IGKC and IFITM1) and 6 downregulated genes (LYZ, AL627171.2, VCAN, CTSS, ITGB2, and PSAP) also appeared in the DEGs from all cell clusters, indicating that classical monocytes were the main driving component for these DEGs ( Figure 4G ). Our results were consistent with those from previous studies. For example, a recent study suggested that the classical monocytes overexpress interferon-induced transmembrane protein 1 in perinatal hepatic inflammation (29), which is a highly specific immune marker of normal endometrial stroma (30).

NK cells are an important group of immune cells, which are widely distributed in human peripheral lymphoid organs and in the circulatory system.

CD56Dim,CD56bright is the main subtype of NK cells in the blood, and CD16+CD56Dim is the main subtype of NK cells in peripheral blood. In our dataset, two main groups of NK cells have been identified, namely c3 and c16 ( Figure 2A ). They are also named CD16^bright CD56^Dim and CD16^Dim CD56^bright according to their markers ( Figures 5C, E ). The average proportions of the two clusters in HC were 10.5% and 0.65%, while in POI, their proportions decreased to 6.6% and 0.64%.

Next, we investigated the DEGs in NK cells between POI and HC. In total, 19 upregulated genes and 8 downregulated genes were identified ( Figure 5A ). KLRC1, which encodes the killer cell lectin-like receptor C1 and is preferentially expressed in NK cells, was the most upregulated DEG in the NK cells of POI patients. This protein forms a complex with another family member, KLRD1/CD94, and has been implicated in the recognition of the MHC class I HLA-E molecules in NK cells. CD3E was another highly upregulated gene. It is a cell surface protein involved in signal transduction to T cell antigen receptor complex; therefore, it plays a crucial role in the adaptive immune response. Enrichment analysis showed that the differentially expressed genes were enriched primarily in the IL-17 signaling pathway, cellular adhesion molecules, and in pathways related to NK cell-mediated cytotoxicity. ( Figure 5B ). Differential gene expression is responsible for the regulation of biological functions such as cell cycle progression, apoptosis, and immune response ( Figure 5D ).

In addition, the upregulation of chemokine CCL5, together with the upregulation of HLA-DPA1, HLA-DRB1, and HLA-DRA suggested that compared with HCs, NK cell-mediated response and cytotoxicity were notable in patients with POI ( Figure 5A ) (31).

B cells are the second-largest type of lymphocytes in humans and mainly participate in humoral immunity. In human peripheral blood, B cells differentiate mainly into the memory B cells and plasma cells. Plasma cells, also known as effector B cells, mainly secrete antibodies and cytokines.

A significantly higher abundance of plasma cells was observed in POI compared with HC ( Figures 6A, B ). Plasma cells were divided into two clusters, namely c17 and c18. In humans, B cells can be divided into naive mature B cells (CD19+CD27-) and memory B cells (CD19+CD27+) based on CD19 and CD27 expression. After antigen stimulation, mature B cells differentiate into either memory B cells or plasma cells. Mature B cells lose some of their surface antigens upon differentiation into plasma cells. We found that B cell markers, including CD19 and MS4A1 (an encoding CD20), were markedly downregulated, while the plasma cell marker, CD38, was highly expressed ( Figure 6E ).

We identified DEGs related to plasma cells between patients with POI and HCs ( Figure 6D ). In patients with POI, S100A9, and CCL5 were upregulated significantly ( Figure 6D ). Enrichment analysis showed that the upregulated genes were mainly enriched in the ribosome, primary immunodeficiency, and T cell receptor signaling pathways, while downregulated genes were mainly enriched in protein processing in the endoplasmic reticulum and protein export ( Figure 6C ).

Treg cells play a key role in maintaining self-tolerance and preventing autoimmune diseases. Previous data showed that Treg abundance is decreased in patients with POI (32).

In this study POI samples were labeled with CD4- CD25+ FOXP3+ and a stable expression of FOXP3 gene. Then, we reanalyzed Treg cells and divided them into 3 subsets ( Figures 7A, B ), including Treg-0, Treg-1, and CD4^- FOXP3⁺ Treg cells ( Figure 7C ). All Treg cell subsets displayed a higher proportion in POI ( Figure 7D ).

DEGs from Treg cells between patients with POI and HCs were identified ( Figure 7E ). Only three genes were upregulated in POI, including IFITM1, JUN, and PPBP. Enrichment analysis showed that the upregulated genes were mainly enriched in B cell receptor signaling pathway, while the downregulated genes were mainly enriched in immune pathways including allograft rejection, autoimmune thyroid disease, and RIG−I−like receptor signaling pathway ( Figure 7F ).

We compared the outgoing and incoming interactions strength to identify cell populations with significant alterations in sending or receiving signals between HC and POI. We observed that naive B cells and classical monocytes were the major affected sources and targets in patients with POI compared with HC ( Figure 8A ). Naive B cells had higher levels of outgoing and incoming interaction strength, whereas classical monocytes showed reduced levels of cell communication, suggesting an enhanced innate immune response and functional impairment of classical monocyte in POI. Besides, CD8+ NKT-like cells and NK cells showed markedly reduced levels of outgoing and incoming interaction strength. These findings are consistent with autoimmune ovarian damage caused by an increased abundance of autoantibody-producing B cells (33). Similarly, our results are consistent with defective monocyte polarization and impaired NK cell activity in POI reported by Hoek et al. (23). We subsequently compared the overall relative strength of the identified signaling pathways in each cell cluster by its aggregated outgoing and incoming signals ( Figure 8B ). Significant signaling pathways were ranked based on differences in the overall communication probability in the indicated networks between HC and POI. Four pathways (RESISTIN, OX40, NMU, and WNT) were unique in HC, while 2 pathways (EPO and TNF) were unique in POI ( Figures 8C, D ). In addition, there were 29 pathways significantly enriched in HC, and 5 pathways significantly enriched in POI ( Figure S2 ). Among them, some pathways that are normally originating from classical monocytes, such as COMPLEMENT and BMP, were originating from NK cells in POI. In contrast, the receivers of some pathways, such as APRIL, LIGHT, and IL4, changed from monocytes to NK cells and naive B cells. Taken together, these results strongly indicate that the intercellular communication networks among the cell clusters were largely disrupted in POI.