Introduction

Female infertility and premature ovarian failure (POF) have become a global concern due to environmental pollution, tremendous life and work stress, and other factors. Women with POF have almost zero ability to conceive naturally, and they respond poorly to various ovulation induction regimens; in addition, the low estrogen levels resulting from the decline in ovarian function not only lead to the early onset of aging symptoms but also increase the risk of osteoporosis and coronary heart disease, and the mortality rate of this population is up to twice that of the general population.1 The incidence of POF has been reported to be 1% among women under the age of forty years and 0.1% among women under the age of thirty years, and the incidence has increased annually.2,3 For a more accurate and comprehensive description of the POF condition, in 2016, the European Society of Human Reproduction and Embryology proposed replacing the term “premature ovarian failure (POF)” POF with “premature ovarian insufficiency (POI)”. POI is defined as a decline in gonadal function due to the lack of functional follicles in the ovaries of women under forty years of age, with elevated gonadotropin levels, a decrease in estrogen levels, and perimenopausal symptoms including scant menstruation or amenorrhea. Early detection and intervention in infertile women with POI with the goal of delaying POI is of great clinical importance.

Transcriptional characterization of human oocytes and granulosa cells (GCs) in early growing and ovulating follicles has been previously described.4 The follicle microenvironment, which consists of cells and extracellular matrix components, is important for oocyte development and ovulation.5,6 For example, granulosa cells (GCs), which secrete cytokines, regulate meiosis in mammalian follicles before ovulation.7,8 However, there have been few reports on the differences in the follicular microenvironments of patients with different ovarian reserve function on the basis of the grouping of patients according to the level of ovarian function and with the inclusion of naturally aged and advanced-aged populations.

Human folliculogenesis is one of the most complex biological processes in vivo, and there is a strong clinical need to reproduce folliculogenesis protocols in vitro for fertility preservation.9 The complexity is reflected not only in the balance between follicular growth and remodeling, but also in the high degree of synchronization of different types of follicular cells within a single follicle, with the common goal of generating competent mature oocytes.10 Here, we included populations with early-onset reduced ovarian reserve function, advanced age and normal ovarian reserve and studied these three groups of women as the research object by obtaining a single follicle after superovulation induced by the hCG trigger in an in vitro fertilization (IVF) cycle. Then, we used single-cell RNA sequencing (scRNA-seq) technology to analyze the cell populations and related factors that play a key role in the occurrence and development of POI in these three groups and to explore the mechanism of POI occurrence.

Compared with traditional bulk RNA sequencing methods, scRNA-seq can comprehensively reveal gene expression among individual cells, thereby identifying different cell types and their cell type-specific gene signatures, and can reveal previously unknown cell types and subtypes, which may contribute to our understanding of rare cells and the functional roles of cells, especially in organs and tissues containing multiple cell types such as ovaries, tumors, lungs, and livers.11–16 Some articles studying single-cell sequencing and the follicular microenvironment, Wu et al4 firstly studied the follicular microenvironment before ovulation and mapped the follicular microenvironment, which provided an important theoretical basis for the subsequent study. In addition, Kristine et al collected somatic cells from preovulatory follicular fluid from patients undergoing in vitro fertilization and used scRNA-seq techniques to study the association between the ovarian sensitivity index and gene expression. This study revealed new information regarding the differences in the follicular gene expression between hypo- and normo-responders.17 However few studies have been conducted on the follicular microenvironment of patients with different ovarian reserve functions. This study was conducted with patients with different ovarian functions as the research object, aiming to reveal the mechanism of the occurrence and development of POI to provide important theoretical support for its intervention.

Materials and MethodsStudy Population

Three discarded follicular fluid samples were collected from three POI patients who underwent IVF/intracytoplasmic sperm injection (ICSI) at the Tianjin Central Hospital of Obstetrics and Gynecology/Nankai University Affiliated Maternity Hospital (Tianjin, China). The diagnostic criteria for POI were (1) age <40 years, (2) oligomenorrhea or amenorrhea for at least four months, and (3) at least two test results showing a basal FSH level of >25 IU/L (testing interval >4 weeks).18 Patients with ovarian dysfunction caused by factors such as ovarian surgery or medication were excluded from our study. Additionally, six discarded follicular fluid samples were obtained from three normal individuals and three older individuals. The normal group included individuals of reproductive age with normal ovarian reserve function who underwent IVF or ICSI due to fallopian tube damage or male factors. The older group included individuals with advanced maternal age (≥40 years).18 All patients underwent controlled ovarian stimulation (COS) using an antagonist protocol. Individuals with type I or type II diabetes, impaired thyroid, renal or liver function, congenital adrenal hyperplasia, endometriosis, hypothalamic amenorrhea, or chromosomal abnormalities were excluded from the study. Experiments involving human samples were carried out according to the Declaration of Helsinki. Ethical approval for the study was granted by the institutional review board of the Reproductive Center of Tianjin Central Hospital of Obstetrics and Gynecology/Nankai University Affiliated Maternity Hospital (approval no. 2022KY020). Informed written consent was obtained from all participants.

Sample Collection and scRNA-Seq Library Preparation

For the patients, recombinant follicle-stimulating hormone (Gonal-F, Merck Serono, Darmstadt, Germany or Puregon, NV Organon, Oss, The Netherlands) treatments were started on days 2–3 of the menstrual cycle or 30–35 days after gonadotropin-releasing hormone (GnRH) agonist (leuprorelin acetate microsphere for injection, Lizhu, China) injection, with the option to adjust the dose according to the response after a 4-day stimulation treatment. In the GnRH antagonist protocol, when the leading follicle reached 13–14 mm in diameter, 0.25 mg of the GnRH antagonist (Orgalutran, NV Organon, Oss, The Netherlands) daily was administered until the day of hCG administration. hCG (Ovidrel, Merck Serono, Darmstadt, Germany) or hCG (chorionic gonadotropin for injection, Lizhu, China) with triptorelin acetate injection (Decapeptyl, Ferring, Kiel, Germany) was administered when one or two leading follicles reached 18 mm in diameter. Transvaginal ultrasound-guided oocyte retrieval was performed 34–36 h later.

The first perforated follicle (diameter: 18–20 mm) containing stage MII cumulus-oocyte complexes (COCs) was used for each subject. The follicular fluid was clear, with no visible blood contamination, and the follicular fluid was collected separately. The cumulus cells, which are the main component of the COC, were gently removed with 1% hyaluronidase and then washed with phosphate-buffered saline (PBS). The cells were then resuspended in GIVF medium (Vitrolife, Sweden). Follicular fluid cells and mural granulosa cells (MGCs) isolated from the same follicle were analyzed by scRNA-seq. The collected cells were centrifuged and the cell pellet was resuspended in PBS (HyClone, Marlborough, MA, USA). Afterward, the GEXSCOPE® red blood cell lysis buffer (RCLB, Singleron Biotechnologies, Nanjing, China) was added and the mixture [Cell: RCLB=1:2 (volume ratio)] was incubated at room temperature for 5–8 min to remove red blood cells. The trypan assay was used to assess cell viability. We counted the cells using a TC20 automatic cell counter (Bio-Rad, Hercules, CA, USA).4 Single-cell preparations were suspended in PBS at a concentration of 1×105 cells/mL and transferred to a microfluidic chip (GEXSCOPE Single Cell RNA-seq Kit, Singleron Biotechnologies). The scRNA-seq libraries were constructed (Singleron Biotechnologies, Jiangsu, China) and sequenced on an Illumina HiSeq X10 instrument using 150-bp paired-end reads. The remaining follicular fluid from the same patient was used to extract cells for subsequent experiments.

Quality Control, Dimensionality Reduction and Clustering (Scanpy)

Scanpy v1.8.1 is used for quality control, dimensionality reduction and clustering in Python 3.7. For each sample dataset, we filtered the expression matrix according to the following criteria: 1) cells with fewer 200 genes or the top 2% of genes were excluded; 2) cells with the top 2% of UMIs were excluded; 3) cells with 20% of mitochondrial content were excluded; and 4) genes expressed in fewer than 5 cells were excluded.

Generation of Single-Cell Gene Expression Matrices

The raw reads were processed to generate gene expression matrices with CeleScope (https://github.com/singleron-RD/CeleScope). First, low-quality reads, poly-A tails, and adaptor sequences were removed from the raw data. The reads were then aligned to the GRCh38 human reference genome (Ensembl version 92) with STAR (version 020201). Using featureCounts (version 1.6.2), reads from the same gene and cell barcodes were pooled together, and the number of unique molecular identifiers (UMIs) per gene per cell was counted. Cell selection was performed with the UMI-tools algorithm. The numbers and quality of the cells were evaluated to identify effective cells and evaluate their gene expression.

Dimensionality Reduction, Unsupervised Clustering, Differential Gene Expression Screening, and Functional Analyses

Cell normalization and filtering were performed with the Seurat package. Principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) were used to reduce data dimensionality and to describe relationships between single cells. For unsupervised clustering, K-nearest neighbors (KNNs) were calculated, and a shared nearest neighbor (SNN) graph was constructed with FindNeighbor. FindClusters was then used to identify cell clusters through modularity optimization methods with the Louvain algorithm. The FindAllMarkers function of the Seurat package and the Wilcoxon rank-sum test (with P values adjusted for multiple testing with Bonferroni correction) were used to identify the differentially expressed genes (DEGs) in each cell type compared with other cell types. Significant DEGs were identified using a |fold change| of >0.25 and a false discovery rate (FDR) of <0.05 as cutoff thresholds. The package clusterProfiler (3.14.3) was used to perform Gene Ontology (GO) analysis on the DEGs, and significant biological processes were identified with the thresholds “p value cutoff=0.05” and “q value cutoff=0.05”.

scRNA-Seq Signature Score

For gene scoring analysis, we used Seurat’s AddModuleScore function to compare different gene profiles in the subgroups.19 The immunosuppressive Characterization score was defined as the average expression of a series of immune checkpoint inhibitors and immunosuppressive molecules, including CD244, CD160, CTLA4, PDCD1, TIGIT, LAYN, LAG3, HAVCR2, CD274, CD47, CD96, ENTPD1, VSIR, BTLA, EBI3, IL2RB, IL2RA, and IL2RG.20,21 The dysfunctional signature score was defined as the average expression of LAG3, HAVCR2, PDCD1, PTMS, FAM3C, IFNG, AKAP5, CD7, PHLDA1, ENTPD1, SNAP47, TNS3, CXCL13, RDH10, DGKH, KIR2DL4, LYST, MIR155HG, RAB27A, CSF1, TNFRSF9, CTLA4, CD27, CCL3, ITGAE, PAG1, TNFRSF1B, GALNT1, GBP2, MYO7A, and TIGIT.22 A list of transcription factors involved in T-cell exhaustion, including BATF, BCL6, BHLHE40, BTLA, CD200, EOMES, ETV1, FOXP3, HIF4A, HOPX, ID2, ID3, IFI16, IKZF, IKZF3, NR4A1, NR4A2, NR4A3, PRDM1, RBPJ, SOX4, STAT3, TBX21, TCF7, TOX, TOX2, VDR, ZBED2, ZNF683, IFI16, DRAP1, and ETS1, was also included.23

Cell–Cell Interactions

Crosstalk analysis between various cell types was performed with the CellChat package (v.0.5.5) with CellChatDB.huma as the ligand–receptor interaction reference database. The function computeCommunProbPathway was used to infer cell–cell communication at the signaling pathway level. Interactions between signaling pathways in various cell subtypes were investigated with the identifyCommunicationPatterns function.

Determination of the Developmental Trajectory

To determine potential lineage differentiation between different GC subtypes, we performed trajectory analysis with the Monocle 2 (version 2.18.0) algorithm.24 To this end, a CellDataSet object was constructed with the newCellDataSet function with the expressionFamily set as the negbinomial.size function. Dimensionality reduction was performed with the DDRTree algorithm (max_components parameters=4) on the basis of the expression of the top 3000 highly variable genes. Next, the cell trajectory was described by means of the orderCells function. The inferred cell trajectories were then visualized by means of the plot_cell_trajectory function. To visualize genes whose expression levels changed along a pseudotime trajectory, we used the plot_pseudotime_heatmap function and transcription factor genes with an adjusted P value of <0.01 and a fold change of >1.5.

Western Blotting

The cells for Western blotting were isolated from the remaining follicular fluid the of same patients using density gradient centrifugation at 2000 rpm for 30 min with human peripheral blood lymphocyte separation medium (Beyotime Biotechnology, China). Following centrifugation, the middle layer was collected and washed with PBS. Western blot experiments were conducted as previously described.25 The membranes were incubated with the following primary antibodies overnight at 4 °C: anti-C5AR1 (1:1000 dilution, cat no. 21316-1-AP, Proteintech) and anti-β-actin (1:5000 dilution, cat no. KM9001, Sungene Biotech, Tianjin, China). The membranes were subsequently incubated with secondary antibodies (1:5000 dilution, cat nos. SE134, Solarbio) for 1 h at room temperature. The bands were visualized with chemiluminescence reagents (Millipore). The protein levels were quantified via ImageJ software (version 1.48; National Institutes of Health, Washington, DC, USA). Three independent Western blots were performed to determine the expression of the different proteins relative to that of β-actin.

Statistical Methods

All raw letter analyses were performed in R, and P<0.05 was considered statistically significant. Data were analyzed with SPSS 26.0 software. ANOVA was used for between-group comparisons of continuous variables. Nonparametric tests were used for between-group comparisons of continuous variables that did not conform to a normal distribution or did not pass the chi-square test. The LSD test was used for two-by-two comparisons among multiple groups, and P<0.05 was considered to indicate statistical significance. Graphs were prepared with GraphPad Prism 9.0.

ResultsBasic Patient Information

To determine the pathogenesis of POI, we performed scRNA-seq on nine follicular fluid samples, three each from POI patients, older individuals, and normal individuals. The age ranges were 28–32, 41–45 and 31–35 years for the normal, older and POI individuals, respectively. Other clinical data, including BMI (body mass index), bFSH (basal follicle-stimulating hormone), AMH (anti-Müllerian hormone) and the number of oocytes are shown in Table 1.

Table 1 The Clinical Information of the 9 Individuals

Data Quality Control

To analyze the differences in single-cell mapping across the groups according to ovarian reserve function and advanced age, we analyzed the cell populations and related factors that play key roles in POI occurrence. To evaluate the mechanism of POI occurrence, we collected 3 samples from each group, and a total of 9 samples were subjected to snRNA-seq. After strict quality control and filtering, the data were subjected to the next step of data analysis (Table 2).

Table 2 Cell Proportion Data for Each Individual

Mapping Preovulatory Follicular Fluid SequencingCell Type Annotation

A total of 87,323 cells were isolated and classified as T cells, B cells, neutrophils, basophils, mononuclear phagocytes (MPs), and GC clusters (Figure 1A). The cell types we found are consistent with published studies on molecular mapping and cellular composition.4 UMAP was used to reveal the distribution of these six clusters in each of the 9 individuals (Figure 1B). The corresponding marker genes for each cell type are shown in Figure 1C. Comparison of the proportions of these six groups in samples of patients with POI, older and normal individual revealed that the fractional ratio group of MP in the POI sample was smaller than the other groups (Figure 1D and Table 2). We constructed heatmaps of the communication networks among the cell subpopulations in the three groups (Figure 1E), and they revealed that the relationship between MPs and granulosa cells was similar in the POI and NOR groups, and stronger in both groups than in the old group. The relationship between MPs and MPs was similar in the POI and old groups, and stronger in both groups than in the NOR group. These findings could suggest that the correlation between MPs and granulosa cells or and MPs cell may plays an important role in the development of POI. Therefore, we conducted further studies on MPS cell functions.

Figure 1 Single-cell atlas of follicular fluid samples from normal individuals, older individuals, and POI patients. (A) UMAP plot of the 87,323 cells from the nine samples based on the five cell clusters, and (B) nine individuals. (C) Bubble plot of the canonical cell type markers of the five main cell types. (D) The distribution of the five main cell clusters in samples from normal individuals, older individuals, and POI patients. (E) The heat maps of the communication networks among cell subpopulations in the three groups.

ScRNA-Seq Reveals MP Heterogeneity Across Normal Individuals, Older Individuals, and POI Patients

To further characterize the differences between patients with POI and normal or older individuals, MP cells were classified into macrophage, mature dendritic cell (DC), monocyte, cDC1, and cDC2 subclusters (Figure 2A and B). The data showed that the number of monocytes in the samples of the POI patients was higher than in the other groups (Figure 2C). Pathway enrichment analysis of the DEGs that found in samples from POI patients compared with the other groups revealed that the response to interferon-gamma was enriched in the downregulated genes (Figure 2D). KEGG pathway analysis showed that in samples from POI patients, antigen processing and presentation pathways were enriched in the downregulated genes (Figure 2E).

Figure 2 MP subclusters in samples from normal individuals, older individuals, and POI patients. (A) UMAP plot of the five MP subtypes in the nine samples. (B) Bubble plot of canonical cell type markers of the five main MP cell types. (C) Line graph of the proportions of the five MP cell subclusters in samples from normal individuals, older individuals, and POI patients. (D and E) GO and KEGG analyses of the differentially expressed genes identified in monocytes from POI patients vs older individuals and the normal group.

Differential Gene Enrichment Analysis of Monocyte Cell Subpopulations Among the Three Groups

Next, we investigated the DEGs in monocytes between POI and the remaining two groups to summarize the top five upregulated and downregulated DEGs among these genes. We found that the gene that was simultaneously downregulated in the POI group compared with the normal and older groups was HLA-DRB5, suggesting that the downregulation of the HLA-DRB5 gene in POI may contribute to the occurrence of POI (Figure 3A). The expression of HLA-DRB5 in both the POI and older groups was lower than that in the normal group (Figure 3B). In addition the results obtained from different samples are generally consistent with the results from different groups (Figure 3C). Moreover, we observed that HLA-DRB5 was expressed mainly in monocytes (Figure 3D).

Figure 3 Differential gene enrichment analysis of monocyte cell subpopulations among the three groups. (A) DEGs identified in the monocyte cell subpopulations between samples from POI patients vs normal individuals and older individuals. (B) The expression of HLA-DRB5 in the POI, older and the normal group. (C) The expression of HLA-DR15 in each sample. (D) HLA-DRB5 was expressed in MPs.

Simulated Time Series Analysis

On the basis of the MP expression profiles among the three groups of samples, the cells were subjected to time series analysis of the transformation of the cell states of different subpopulations in the process of disease development and treatment, and an algorithm was used to simulate the process of the transformation of the cell states to determine the overall trajectory of the changes in the expression of the genes and to identify the genes with significant differences to explore the mechanism of POI occurrence. We used the algorithm to simulate the transition of the cell state during the treatment process. The temporal differentiation trajectory revealed that different monocytes play important roles in the beginning and end stages of differentiation (Figure 4A). The C5AR1 gene is highly expressed in monocytes and C5AR1 expression may play an important role in disease development (Figure 4B). C5AR1 expression gradually decreased with pseudotime (Figure 4B1–B3).

Figure 4 Pseudotime analysis of MPs from normal individuals, older individuals, and POI patients. (A) Pseudotime analysis reveals the progression of MPs. (B1) Heatmap of AUC matrix clustering of top 30 regulon in each cell. (B2) Top 8 genes that change with pseudotime of the state 1–3. (B3) Top 8 genes that change with pseudotime in different clusters.

Expression of C5AR1 in Follicle Microenvironments

In MP cells, the expression of C5AR1 was greater in the POI group than in the old and NOR groups, and the old group had presented greater than the NOR group did. The expression ratios of the POI and old groups were similar, and the two groups were greater than those of the NOR group (Figure 5A). In addition the results obtained from different samples are generally consistent with the results from different groups (Figure 5B). C5AR1 was highly expressed in POI follicular fluid cells, and these results were further confirmed by Western blotting (Figure 5C).

Figure 5 Expression of C5AR1 in follicle microenvironments from three group. (A) The expression of C5AR1 in the three groups. (B) The expression of C5AR1 in the nine samples. (C) The expression of C5AR1 in follicle fluid cells from the three groups was detected using Western blotting. (*P<0.05, **P<0.01).

Discussion

Single-cell transcriptome sequencing allows simultaneous sequencing of the transcriptome of each individual cell in a population of nearly one million cells, greatly improving our understanding of intercellular heterogeneity. Single-cell sequencing technology allows better identification of specific subpopulations of cells and exploration of their potential roles in physiopathology.26

We constructed preovulatory follicular fluid cellular profiles of infertile women with different ovarian reserve functions by single-cell sequencing. Granulosa cells (GCs), which are follicular somatic cells that provide essential nutrients and growth factors, are also known to secrete the steroids needed for the progression of folliculogenesis. GC dysfunction can initiate follicle atresia and apoptosis and can eventually lead to POI.27,28 Previous studies have shown that dysregulation of Epg5 in GCs can trigger POI pathogenesis.29 However, owing to the small number of other cells, the mechanism of intercellular action has rarely been studied, limiting further research to some extent. We found that the interaction between MPs and granulosa cells was stronger in the POI group than in the other two groups. These findings suggest that the correlation between MPs and granulosa cells may lead to the development of POI.

Our research revealed that, compared with that in each of the other groups, the response to interferon-gamma in the POI group was enriched among the downregulated genes. Our preliminary findings suggest that excess inflammation plays a crucial role in both ovarian aging and POI.30 We hypothesized that an attenuated interferon response might lead to a downstream inflammatory response, contributing to the development of POI. Interestingly, a high percentage of the cells in the MPs were monocytes. Comparison revealed that the gene HLA-DRB5, which was simultaneously downregulated in the POI group compared with the normal and older groups, was highly expressed in monocytes. Therefore, it is hypothesized that HLA-DRB5 plays an important role in delaying the process of POI. Although autoimmune abnormalities appear to be involved in the development of POI, only a few studies with respect to human leukocyte antigen (HLA) exist. Human leukocyte antigen-DR15 is a haplotype that results in the expression of DRB1 and DRB5, most frequently HLA-DRB1*15:01 and HLADRB5*01:01. Previous studies have demonstrated the association of the HLA-DR15 haplotype with multiple sclerosis, a chronic T-cell-mediated autoimmune disease.31,32 HLA-DR15 status has also been reported to be an important marker in autoimmune bone marrow disease.33 Research on the role of HLA-DR15 in multiple sclerosis has revealed that patients display T-cell proliferative responses to myelin-associated autoantigens, which are related to the presence and level of HLA-DR15 expression.34 In addition, a recent study revealed a lack of HLA-DRB5 expression and a lower antigen-presenting ability in ST memory B-cell subsets, as well as the macrophage and DC subsets of ACPA-RA.35 Inconsistent with the findings of previous studies, these results imply that a specific HLA haplotype (A * 24:02-C * 03:03-B * 35:01) constitutes a susceptibility factor for apparently isolated POI in Japanese women.36 Therefore, we speculate that autoimmune B and T-cell subsets are not major players in the pathogenesis of POI. The presence of the HLA-DRB5 genotype in POI patients might also indicate autoproliferation of Th1 cells to present autoantigens, contributing to disease activity.

In addition, our study revealed that different monocytes play important roles in the beginning and end stages of differentiation. C5AR1 could be used to predict the development of this disease to a certain extent. Specifically, the immunomodulatory effects of Complement C5a occur through the activation of C5a receptor 1 (C5aR1) in myeloid cells, such as macrophages, monocytes, DCs, and neutrophils.37 C5a-C5aR1 signaling frequently reshapes the immunosuppressive tumor microenvironment and increases the potential for cancer progression via molecular and cellular events.38 Recent studies have revealed the contribution of C5a-C5aR1 signaling to tumor-associated immune responses through the induction of bioactive molecules, including IL-1β, IL-6, and TNF-α.39,40 Furthermore, crosstalk between C5aR1 and toll-like receptors (TLRs) enhances the production of TNF-α and IFN-γ by natural killer and natural killer T cells.40 C5aR1 blockade also impairs tumor growth and diminishes the recruitment of myeloid-derived suppressor cells (MDSCs) into tumors in lung cancer-bearing mice.

Studies of the ovaries of women have revealed that the complement C5a/C5aR1 axis contributes to the programming of the immunosuppressive phenotype of tumor-associated macrophages in solid tumors and represents a promising immunomodulatory target for treating high-grade serous ovarian cancer.41 The activation of C5a-C5aR1 signaling has also been shown to have prognostic effects. Overexpression of C5aR1 in non-small cell lung cancer predicts poor prognosis.42 In breast cancer, C5aR1 expression is associated with increased tumor size, increased proliferation rates, the presence of lymph node metastasis, and advanced clinical stages.43 Similarly, patients with C5aR1-positive breast cancer have lower survival rates than those with C5aR1-negative breast cancer. C5aR1 expression has also been associated with the disease prognosis in patients with gastric cancer, hepatocellular carcinoma, and urothelial and renal cell carcinoma, with lower overall survival rates among patients with high tumor C5aR1 expression.44–46 Therefore, the abundance of C5aR1 could predict an inferior prognosis in HGSCs, and incorporating PD-L1 may serve as a novel predictive biomarker to guide therapeutic options. However, the mechanism of its effect on POI and its predictive role have not been well studied. Interestingly, in this study we found that C5aR1 was overexpressed in POI patients. To a certain extent, C5AR1 could be used to predict the development of a POI. However, their correlation and specific mechanism of action require in vivo and in vitro verification with larger sample sizes.

This study has several limitations. First, the sample size was small; therefore, a large-scale study is needed to verify our findings. Second, to ensure the consistency of follicular fluid, we did not include COC granulosa cells; that these cells play an important role in POI cannot be ignored, and subsequent studies will need to include these cells to further analyze and validate the mechanism of POI occurrence.

In conclusion, our study provides the first comprehensive single-cell transcriptomic atlas of the microenvironment in human preovulatory follicles according to ovarian reserve function. Our findings revealed that the interaction between monocytes and granulocytes may contribute to the development of POI. We demonstrated a lack of HLA-DRB5 expression as well as impaired “antigen processing and presentation activity” in POI. To a certain extent, C5AR1 could be used to predict the development of a POI. However the correlation between C5AR1 and POI requires verification with more samples.

Data Sharing Statement

Publicly available datasets were analyzed in this study. The raw data have been uploaded to the Genome Sequence Archive for Human (https://ngdc.cncb.ac.cn/gsa-human/submit/hra/subHRA007647/finishedOverview). Other data will be made available to the editors of the journal for review or upon request.

Informed Consent Statement

Informed consent was obtained from all of the subjects involved in the study.

Acknowledgments

The authors thank all the doctors, nurses, and embryologists at the Reproductive Medicine Center of Tianjin Central Hospital of Gynecology Obstetrics for their help in collecting the data.

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Funding

This study was funded by the Tianjin Health Research Project (Grant No. TJWJ2022QN086), the Cellular Ecology Haihe Laboratory “Unveils the Commander-in-Chief” Program (Grant No. HH22KYZX0024) and the Tianjin Science and Technology Planning Project (Grant Number: 21JCZDJC00330).

Disclosure

All authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References

1. Qin C, Chen Y, Lin Q, Yao J, Wu W, Xie J. The significance of polymorphism and expression of oestrogen metabolism-related genes in Chinese women with premature ovarian insufficiency. Reprod Biomed Online. 2017;35(5):609–615. doi:10.1016/j.rbmo.2017.07.007

2. Laml T, Schulz-Lobmeyr I, Obruca A, Huber JC, Hartmann BW. Premature ovarian failure: etiology and prospects. Gynecol Endocrinol. 2000;14(4):292–302. doi:10.3109/09513590009167696

3. Webber L, Davies M, Anderson R, et al. ESHRE Guideline: management of women with premature ovarian insufficiency. Hum Reprod. 2016;31(5):926–937. doi:10.1093/humrep/dew027

4. Wu H, Zhu R, Zheng B, et al. Single-Cell Sequencing Reveals an Intrinsic Heterogeneity of the Preovulatory Follicular Microenvironment. Biomolecules. 2022;12(2):231. doi:10.3390/biom12020231

5. Kinnear HM, Tomaszewski CE, Chang AL, et al. The ovarian stroma as a new frontier. Reproduction. 2020;160(3):R25–r39. doi:10.1530/rep-19-0501

6. Zhou Y, Ding X, Wei H. Reproductive immune microenvironment. J Reprod Immunol. 2022;152:103654. doi:10.1016/j.jri.2022.103654

7. Jaffe LA, Egbert JR. Regulation of Mammalian Oocyte Meiosis by Intercellular Communication Within the Ovarian Follicle. Annu Rev Physiol. 2017;79(1):237–260. doi:10.1146/annurev-physiol-022516-034102

8. Li R, Albertini DF. The road to maturation: somatic cell interaction and self-organization of the mammalian oocyte. Nat Rev Mol Cell Biol. 2013;14(3):141–152. doi:10.1038/nrm3531

9. Fiorentino G, Cimadomo D, Innocenti F, et al. Biomechanical forces and signals operating in the ovary during folliculogenesis and their dysregulation: implications for fertility. Hum Reprod Update. 2023;29(1):1–23. doi:10.1093/humupd/dmac031

10. Telfer EE, Andersen CY. In vitro growth and maturation of primordial follicles and immature oocytes. Fertil Steril. 2021;115(5):1116–1125. doi:10.1016/j.fertnstert.2021.03.004

11. Hwang B, Lee JH, Bang D. Single-cell RNA sequencing technologies and bioinformatics pipelines. Exp Mol Med. 2018;50(8):1–14. doi:10.1038/s12276-018-0071-8

12. MacParland SA, Liu JC, Ma XZ, et al. Single cell RNA sequencing of human liver reveals distinct intrahepatic macrophage populations. Nat Commun. 2018;9(1):4383. doi:10.1038/s41467-018-06318-7

13. Travaglini KJ, Nabhan AN, Penland L, et al. A molecular cell atlas of the human lung from single-cell RNA sequencing. Nature. 2020;587(7835):619–625. doi:10.1038/s41586-020-2922-4

14. Astapova O, Minor BMN, Hammes SR. Physiological and Pathological Androgen Actions in the Ovary. Endocrinology. 2019;160(5):1166–1174. doi:10.1210/en.2019-00101

15. Richani D, Dunning KR, Thompson JG, Gilchrist RB. Metabolic co-dependence of the oocyte and cumulus cells: essential role in determining oocyte developmental competence. Hum Reprod Update. 2021;27(1):27–47. doi:10.1093/humupd/dmaa043

16. Gong X, Zhang Y, Ai J, Li K. Application of Single-Cell RNA Sequencing in Ovarian Development. Biomolecules. 2022;13(1):47. doi:10.3390/biom13010047

17. Roos K, Rooda I, Keif RS, et al. Single-cell RNA-seq analysis and cell-cluster deconvolution of the human preovulatory follicular fluid cells provide insights into the pathophysiology of ovarian hyporesponse. Front Endocrinol. 2022;13:945347. doi:10.3389/fendo.2022.945347

18. Ferraretti AP, La Marca A, Fauser BC, Tarlatzis B, Nargund G, Gianaroli L. ESHRE consensus on the definition of ‘poor response’ to ovarian stimulation for in vitro fertilization: the Bologna criteria. Hum Reprod. 2011;26(7):1616–1624. doi:10.1093/humrep/der092

19. Satija R, Farrell JA, Gennert D, Schier AF, Regev A. Spatial reconstruction of single-cell gene expression data. Nat Biotechnol. 2015;33(5):495–502. doi:10.1038/nbt.3192

20. Wykes MN, Lewin SR. Immune checkpoint blockade in infectious diseases. Nat Rev Immunol. 2018;18(2):91–104. doi:10.1038/nri.2017.112

21. Hu Q, Hong Y, Qi P, et al. Atlas of breast cancer infiltrated B-lymphocytes revealed by paired single-cell RNA-sequencing and antigen receptor profiling. Nat Commun. 2021;12(1):2186. doi:10.1038/s41467-021-22300-2

22. Li H, van der Leun AM, Yofe I, et al. Dysfunctional CD8 T Cells Form a Proliferative, Dynamically Regulated Compartment within Human Melanoma. Cell. 2019;176(4):775–789.e18. doi:10.1016/j.cell.2018.11.043

23. Hu H, Miao YR, Jia LH, Yu QY, Zhang Q, Guo AY. AnimalTFDB 3.0: a comprehensive resource for annotation and prediction of animal transcription factors. Nucleic Acids Res. 2019;47(D1):D33–d38. doi:10.1093/nar/gky822

24. Trapnell C, Cacchiarelli D, Grimsby J, et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol Apr. 2014;32(4):381–386. doi:10.1038/nbt.2859

25. Han Y, Luo H, Wang H, Cai J, Zhang Y. SIRT1 induces resistance to apoptosis in human granulosa cells by activating the ERK pathway and inhibiting NF-κB signaling with anti-inflammatory functions. Apoptosis. 2017;22(10):1260–1272. doi:10.1007/s10495-017-1386-y

26. van Beek JJP, Rescigno M, Lugli E. A fresh look at the T helper subset dogma. Nat Immunol. 2021;22(2):104–105. doi:10.1038/s41590-020-00858-1

27. Sun Z, Zhang H, Wang X, et al. TMCO1 is essential for ovarian follicle development by regulating ER Ca(2+) store of granulosa cells. Cell Death Differ. 2018;25(9):1686–1701. doi:10.1038/s41418-018-0067-x

28. Yeung CK, Wang G, Yao Y, et al. BRE modulates granulosa cell death to affect ovarian follicle development and atresia in the mouse. Cell Death Dis. 2017;8(3):e2697. doi:10.1038/cddis.2017.91

29. Liu W, Chen M, Liu C, et al. Epg5 deficiency leads to primary ovarian insufficiency due to WT1 accumulation in mouse granulosa cells. Autophagy. 2023;19(2):644–659. doi:10.1080/15548627.2022.2094671

30. Han Y, Yao R, Yang Z, et al. Interleukin-4 activates the PI3K/AKT signaling to promote apoptosis and inhibit the proliferation of granulosa cells. Exp Cell Res. 2022;412(1):113002. doi:10.1016/j.yexcr.2021.113002

31. Rasmussen HB, Kelly MA, Clausen J. Additive effect of the HLA-DR15 haplotype on susceptibility to multiple sclerosis. Mult Scler. 2001;7(2):91–93. doi:10.1177/135245850100700203

32. Prat E, Tomaru U, Sabater L, et al. HLA-DRB5*0101 and -DRB1*1501 expression in the multiple sclerosis-associated HLA-DR15 haplotype. J Neuroimmunol. 2005;167(1–2):108–119. doi:10.1016/j.jneuroim.2005.04.027

33. Saunthararajah Y, Nakamura R, Nam JM, et al. HLA-DR15 (DR2) is overrepresented in myelodysplastic syndrome and aplastic anemia and predicts a response to immunosuppression in myelodysplastic syndrome. Blood. 2002;100(5):1570–1574. doi:10.1182/blood.V100.5.1570.h81702001570_1570_1574

34. Mohme M, Hotz C, Stevanovic S, et al. HLA-DR15-derived self-peptides are involved in increased autologous T cell proliferation in multiple sclerosis. Brain. 2013;136(Pt 6):1783–1798. doi:10.1093/brain/awt108

35. Wu X, Liu Y, Jin S, et al. Single-cell sequencing of immune cells from anticitrullinated peptide antibody positive and negative rheumatoid arthritis. Nat Commun. 2021;12(1):4977. doi:10.1038/s41467-021-25246-7

36. Ayabe T, Ishizuka B, Maruyama T, et al. Association of primary ovarian insufficiency with a specific human leukocyte antigen haplotype (A*24:02-C*03:03-B*35:01) in Japanese women. Sex Dev. 2011;5(5):235–240. doi:10.1159/000330122

37. Klos A, Tenner AJ, Johswich KO, Ager RR, Reis ES, Köhl J. The role of the anaphylatoxins in health and disease. Mol Immunol. 2009;46(14):2753–2766. doi:10.1016/j.molimm.2009.04.027

38. Markiewski MM, DeAngelis RA, Benencia F, et al. Modulation of the antitumor immune response by complement. Nat Immunol. 2008;9(11):1225–1235. doi:10.1038/ni.1655

39. Cavaillon JM, Fitting C, Haeffner-Cavaillon N. Recombinant C5a enhances interleukin 1 and tumor necrosis factor release by lipopolysaccharide-stimulated monocytes and macrophages. Eur J Immunol. 1990;20(2):253–257. doi:10.1002/eji.1830200204

40. Corrales L, Ajona D, Rafail S, et al. Anaphylatoxin C5a creates a favorable microenvironment for lung cancer progression. J Immunol. 2012;189(9):4674–4683. doi:10.4049/jimmunol.1201654

41. Zhang C, Cao K, Yang M, et al. C5aR1 blockade reshapes immunosuppressive tumor microenvironment and synergizes with immune checkpoint blockade therapy in high-grade serous ovarian cancer. Oncoimmunology. 2023;12(1):2261242. doi:10.1080/2162402x.2023.2261242

42. Gu J, Ding JY, Lu CL, et al. Overexpression of CD88 predicts poor prognosis in non-small-cell lung cancer. Lung Cancer. 2013;81(2):259–265. doi:10.1016/j.lungcan.2013.04.020

43. Imamura T, Yamamoto-Ibusuki M, Sueta A, et al. Influence of the C5a-C5a receptor system on breast cancer progression and patient prognosis. Breast Cancer. 2016;23(6):876–885. doi:10.1007/s12282-015-0654-3

44. Hu WH, Hu Z, Shen X, Dong LY, Zhou WZ, Yu XX. C5a receptor enhances hepatocellular carcinoma cell invasiveness via activating ERK1/2-mediated epithelial-mesenchymal transition. Exp Mol Pathol. 2016;100(1):101–108. doi:10.1016/j.yexmp.2015.10.001

45. Wada Y, Maeda Y, Kubo T, Kikuchi K, Eto M, Imamura T. C5a receptor expression is associated with poor prognosis in urothelial cell carcinoma patients treated with radical cystectomy or nephroureterectomy. Oncol Lett. 2016;12(5):3995–4000. doi:10.3892/ol.2016.5137

46. Maeda Y, Kawano Y, Wada Y, et al. C5aR is frequently expressed in metastatic renal cell carcinoma and plays a crucial role in cell invasion via the ERK and PI3 kinase pathways. Oncol Rep Apr. 2015;33(4):1844–1850. doi:10.3892/or.2015.3800