A single-cell gene expression atlas of human follicular aspirates: Identification of leukocyte subpopulations and their paracrine factors.

A schematic ( Fig. 1a ) depicts an experimental procedure from sample preparation to identification of cell subpopulations. After integrating four follicular aspirate datasets from individual IVF patient samples, 7,609 cells filtered were dimensionally reduced and visualized on a uniform manifold approximation and projection (UMAP). Unsupervised clustering analysis using single-cell gene expression profiles revealed 19 cell clusters (clusters 0-18, Fig. 1b ). Feature plots for CYP11A1 and PTPRC represent follicular cells and leukocytes, respectively, in the UMAP ( Fig. 1c ). Each dot represents an individual cell, and the color gradient reflects gene expression level. The cells in the 19 clusters ( Fig. 1b ) were then re-grouped using the expression of cell type-specific marker genes noted in Supplementary Table 3 and described below. To identify follicular cells, we used CYP11A and HSD3B2 as these genes are involved in the first two steps of steroidogenesis, converting cholesterol to pregnenolone and pregnenolone to progesterone, respectively. Clusters 0, 2, 9, 11, 12, and 18 showed high levels of these genes ( Fig. 1c ), indicating steroidogenic follicular cells. Follicular cells consist of granulosa and theca cells. Thus, the follicular cell clusters were then re-grouped into the two cell types, granulosa (clusters 0, 2, 9, 11, and 12 ) and theca cells (cluster 18), using known granulosa cell marker genes [e.g., CYP19A1 ( 28 ) and INHA ( 29 )] and theca cell marker genes [e.g., COL1A1, COL1A2 , and COL3A1 ( 30 )] ( Fig. 1c and Supplementary Fig. 1 ). Next, leukocytes, red blood cells (RBCs), and endothelial cells were searched among the remaining clusters. Leukocyte-specific PTPRC was detected in clusters 1, 3, 4, 5, 6, 7, 8, 10, 13, 14, 15, 16, and 17 ( Fig. 1c ). The RBC-specific HBA1 expression was shown only in the cluster 16 ( Fig. 1c ). Meanwhile, none of the 19 clusters of the cell pool showed the significant expression of endothelial cell specific VWF and TIE1 ( Supplementary Fig. 1 ), suggesting few endothelial cells in follicular aspirates collected for our experiments. PTPRC-positive cells were then sub-clustered into ten subsets of leukocytes based on the expression of marker genes specific to various types of leukocytes ( Supplementary Table 3 and Fig. 1d ). These included M1- and M2-macrophages, NK and NKT cells, helper and cytotoxic T cells, neutrophils, baso/eosinophils, B cells, and dendritic cells ( Supplementary Fig. 2 ). In the end, the initial unsupervised 19 clusters were re-grouped into 13 groups as noted in Table 1 . These groups were re-labeled in the UMAP ( Fig. 1d ). These include granulosa cell (GC) and theca cell (TC) groups, the RBC/platelet group, and ten leukocyte subgroups. Table 1 also includes the number and proportion of cells in each group. Follicular cells account for 38.2% (37.9% and 0.3% of GC and TC, respectively), whereas RBC/platelets and leukocytes account for 61.8% of the total cell population. The percentages of leukocytes in four individual IVF patients are 86%, 47%, 56%, and 63%, respectively ( Supplementary Fig. 3a ). Differences in cell numbers and percentages among the four patients are shown in a dimensional reduction plot ( Supplementary Fig. 3b ) and two barplots ( Supplementary Fig. 3c ). To determine whether leukocyte populations identified from IVF patient samples were also present in the ovulatory follicle of normally cycling women, immunohistochemical analyses were performed to detect CD68 protein in dominant follicle tissues collected before and at defined hours after hCG administration. CD68 was chosen as this gene is known to be highly expressed in specific leukocyte populations, including macrophages and monocytes ( 31 ). Our scRNA-seq data analysis also showed the high expression of CD68 mRNA mainly in macrophages in IVF patient samples ( Fig. 2a ). Immunohistochemical analysis using dominant follicle tissues collected before and at defined times after hCG administration not only confirmed the presence of CD68-positive leukocytes but also indicate the increase in the number of CD68-positive leukocytes after hCG stimulation. For instance, before the LH surge, immuno-positive staining for CD68 was localized to only a few cells in the theca layer ( Fig. 2b - i ). However, after hCG administration, the distribution of CD68-positive cells was increased in the thecal layer ( Fig.2b - ii , iii ,& iv ). In the late ovulatory follicle, CD68-positive cells were also observed in the granulosa cell layer, although the expression was limited, predominantly on the basal side of the granulosa cell layer (arrows in Fig. 2b - iii & iv ). After ovulation, CD68-positive cells were detected throughout the granulosa-lutein and thecal cell layers (arrows in Fig. 2b - v ). Five clusters identified as GCs in scRNA-seq data indicate the heterogeneity of granulosa cell populations in human follicular aspirates. To delineate the developmental state of these cells, scRNA-seq data of the 5 clusters in the GC group were subjected to RNA velocity analysis. According to developmental projections among the 5 clusters, they were labeled as GC 1, GC 2, GC 3, GC 4, and GC 5. Results predicted that the differentiation of granulosa cells in response to hCG stimulation was initiated from the GC 1 group, directed toward GC 2 and then GC 5 and in the other branch, directed toward GC 3, GC 4, and GC 5 in a sequence ( Fig. 3a ). As shown in a histogram ( Fig. 3b ), the expression of genes known to be expressed in granulosa cells of early ovulatory follicles, such as PGR ( 23 ), INHBA/INHBB ( 31 ), and BMP3 ( 32 ), was found in GC 1 and GC 2 (with higher expression) groups. Whereas remaining GC 3, GC 4, and GC 5 showed the relatively higher expression for late ovulatory genes, including RUNX2 ( 33 , 34 ), ADAMTS1 ( 35 , 36 ), HSD11B1 ( 36 , 37 ), EDN2 ( 38 ), SLCO2A1 ( 23 ), and CD24 ( 39 ). These genes are known to be involved in various aspects of the ovulatory process including transcription regulation ( RUNX2 ), follicular constriction ( EDN2 ), prostaglandin synthesis and production ( CD24 and SLCO2A1 ), and cortisol production ( HSD11B1 ). Since cytokines are essential secretory proteins from leukocytes that can act on resident cells in targeted tissues as paracrine factors, we searched for cytokines exclusively expressed in leukocytes and their corresponding receptors expressed in granulosa cells. This is to determine how leukocytes directly modulate granulosa cell function, thereby promoting ovulation and luteinization in the human ovary. The expression of genes encoding cytokines, including interleukins (ILs), interferons (IFNs), and tumor necrosis factors (TNFs), and their respective receptors were screened in leukocytes and granulosa cells, respectively. As shown in dot plots in Fig. 4a & b , the intensive or higher percent expression of cytokines, IL1B, IL16, IFNG, TNF, LTA, LTB , and TNFSF10 was shown predominantly in several groups of leukocytes but barely in granulosa cell groups. In contrast, the expression of their receptors, IL1R1, IL1RAP, CD9, IFNGR1, IFNGR2, TNFRSF1A, LTBR , and TNFRSF10B was detected in granulosa and theca cell groups. The leukocyte group-specific expression of the ligands and the expression of their receptors in granulosa cell groups were also denoted in violin and feature plots ( Supplementary Fig. 3 ). To further verify that leukocyte-derived cytokines from IVF patient samples were also present in the ovulatory follicle of normally cycling women, the expression of IL1B was assessed using the dominant follicles collected at the late ovulatory phase. Immuno-positive staining for IL1B was localized to the cells inside the blood vessel or in the thecal layer of late ovulatory follicles, indicating that IL1B is secreted from leukocytes present in the ovulatory follicle (arrowheads in Fig. 4c ). Meanwhile, mRNA for IL1R1 , a receptor for IL1B, was detected predominantly in granulosa and theca cells in scRNA-seq data. Because the antibody for IL1R1 used in the present study was only working for Western blot application, the presence of IL1R1 protein in granulosa cells was determined using primary granulosa/lutein cells collected from IVF patient samples. Western blot analyses confirmed the presence of IL1R1 protein in granulosa/lutein cells. Moreover, hCG treatment increased the level of IL1R1 protein when cultured for 36h, suggesting that IL1R1 expression is upregulated by ovulatory hCG stimulation in granulosa cells of preovulatory follicles ( Fig.4d ). In addition to cytokines, other paracrine factors of leukocyte origin and their corresponding receptors in granulosa cells were discovered in the datasets of the four IVF patient samples. EGF-like peptides, AREG and EREG , were abundantly expressed in not only granulosa cells but also leukocytes (mainly in neutrophils and macrophages, respectively), whereas their shared receptor EGFR expression was specific to granulosa cells ( Fig. 4e ; Supplementary Fig. 4 ). Likewise, the expression of PTGS2 , the rate-limiting enzyme in the PG synthesis pathway, was detected in granulosa cells and leukocytes, mainly macrophages and neutrophils, while the expression of PG receptors, except PTGER2 , was specific to granulosa cells ( Fig. 4e ; Supplementary Fig. 4 ). The expression of urokinase-type plasminogen activator ( PLAU ) was detected explicitly in M2-macrophages, neutrophils, and dendritic cells, while its receptor PLAUR expression was detected in granulosa cells ( Fig. 4f ; Supplementary Fig. 4 ). Similarly, neuregulin-1 ( NRG1 ) expression was predominantly detected in M1- and M2-macrophages and neutrophils, whereas its receptors, EBRR2, ERBB3 , and ERBB4 , were detected in granulosa cells ( Fig. 4f ; Supplementary Fig. 4 ). Tissue proteases play an instrumental role in the periovulatory process by breaking down the basement membrane between granulosa cells/thecal layer and under the germinal epithelium to aid rupture of the exterior follicular wall at ovulation ( 40 ). Accumulating evidence showed that leukocytes secrete various tissue proteases ( 13 - 15 ). Therefore, we searched the expression of tissue proteases, including matrix metalloproteinases (MMPs), a disintegrin and metalloproteinases (ADAMs), ADAM with thrombospondin motif (ADAMTSs), and tissue inhibitor of MMPs (TIMPs) in leukocytes in the human follicular aspirate obtained at 36 h after hCG administration, and then visualized them in dot ( Fig. 5 ) and feature plots ( Supplementary Fig. 5 ). MMP2, MMP9, MMP10, MMP14, MMP17, MMP19, MMP23B , and MMP25 were detected in follicular aspirates ( Fig. 5a ). MMP9 and MMP25 mRNA expression was exclusively in leukocytes, while MMP10 and MMP2 expression was predominant in granulosa cells and theca cells, respectively. Meanwhile, MMP14, MMP17, MMP19 , and MMP23B expression was detected in both follicular cells and leukocytes, albeit with differential distribution among different populations of cells. To further support this finding, immunohistochemical analysis was performed for MMP9 in dominant follicle tissues collected during the late ovulatory period from normally cycling women. Positive immunostaining for MMP9 was localized to leukocytes distributed in/around the blood vessel and mainly in the thecal layer. MMP9 positive staining was also detected in leukocytes in the granulosa cell layer (arrowhead, Fig. 5b ). Of ADAM isoforms detected in the datasets, ADAM8, ADAM19 , and ADAM28 expression was predominant or exclusive in leukocyte subsets: macrophages, T, B, NK, and NKT cells. In contrast, ADAM9, ADAM10, and ADAM15 expression was detected in both follicular cells and leukocytes ( Fig. 5c ; Supplementary Fig. 5 ). Among the ADAMTS family, ADAMTS1 expression was detected in both follicular cells and various types of leukocytes, though overall expression levels were higher in follicular cells, particularly granulosa cells, than leukocytes. Other isoforms of ADAMTS detected in follicular aspirates were ADAMTS4 , ADAMTS5, ADAMTS6, ADAMTS10 , and ADAMTS17 , displaying specific distribution patterns and expression levels among follicular cells and leukocytes ( Fig. 5d ; Supplementary Fig. 5 ). Among four TIMPs, TIMP1, TIMP2 , and TIMP3 mRNA was detected in cells in follicular aspirates. Their expression levels were relatively higher than those of proteases and detected both in granulosa and theca cells as well as various types of leukocyte groups. TIMP2 expression is abundantly found in M1- and M2-macrophage groups, while TIMP3 expression is mainly distributed in follicular cells. The expression of TIMP4 is undetectable in human follicular aspirates ( Fig. 5e ; Supplementary Fig. 5 ). Previous studies using the monkey ovary reported increases in the expression of VEGFA , placental growth factor ( PGF ), insulin-like growth factor-1 ( IGF1 ), fibroblast growth factor-2 (FGF2), and Thrombospondin-1 ( THBS1 ) in granulosa and theca cells of ovulatory follicles ( 4 , 41 - 43 ). In those studies, these factors are suggested to be involved in promoting angiogenesis during ovulation and luteinization in the monkey ovary ( 4 , 41 - 43 ). As shown in dot, violin, and feature plots ( Fig. 6 ; Supplementary Fig. 6 ), follicular (granulosa and theca) cells expressed all the VEGF isoforms and PGF , but VEGFD mRNA was not detectable in follicular aspirates. Especially, the expression of VEGFC, PGF, FGF2 , and IGF1 was confined to follicular cells. Meanwhile, leukocytes also expressed VEGFA (macrophages and granulocytes), VEGFB (all types of cells in the aspirate), and THBS1 (macrophages, granulocytes, and dendritic cells).

Ovarian hyper-stimulation was induced by administering recombinant human FSH in individualized doses to patients. IVF patients were then given hCG ( 10 , 11 ) on Days 9 to 11. Dominant follicles 18-36 mm in diameter were aspirated 36h later from four women aged 26-38 with non-ovarian etiologies (male factor or egg donor), as routinely performed during the in vitro fertilization (IVF) cycle at the Bluegrass Fertility Clinic (Lexington, Kentucky) ( 11 , 23 ). The collection protocol for follicular aspirate samples from in vitro fertilization (IVF) patients was approved by the Institutional Review Board of the University of Kentucky Office of Research Integrity. The patient characteristics are shown in Supplementary Table 1 . Immediately after retrieval of cumulus-oocyte complexes, the cells in the aspirates were subjected to red blood cell lysis in an ACK RBC lysis buffer to remove erythrocytes. The remaining cells were resuspended and filtered through a 40-μm cell strainer to dissociate cell clumps into single cells. Whole follicles were collected from patients across the periovulatory period as previously described ( 23 ). Women (age 30-38 years) exhibiting regular menstrual cycles and had not taken hormonal contraceptives for at least three months before their enrollment in the study underwent laparoscopic sterilization. Women were monitored by transvaginal ultrasound for 2 to 3 menstrual cycles before surgery to ascertain cycle regularity and monitor the dominant follicle’s growth during the follicular phase. These patients were divided into three groups: Pre-(PO), early (EO), and late ovulatory (LO) phases. In the PO group, patients were subjected to surgery to collect follicles that reached >14 mm and ≤17.5 mm in diameter before the endogenous LH surge. The remaining patients were given recombinant hCG (Ovitrelle, 250 μg) and were divided into two groups: EO (surgery between 12h and 18h post-hCG) and LO (surgery between 18h and 34h post-hCG). To confirm whether these patients followed a normal hormonal pattern before the LH surge or after hCG administration, blood samples were taken at surgery and measured for serum progesterone and estradiol (the patient characteristics in Supplementary Table 2 ). The whole dominant follicles were excised by surgery (laparoscopic scissors) and processed for either immunohistochemical or gene expression analysis. The follicle was bisected, and mural granulosa cells were gently scraped off from the interior of the follicle by small tissue forceps. The follicular fluid and cell suspension were combined and centrifuged at 500 xg to pellet and collect granulosa cells. The Human Ethics Committee of the Sahlgrenska Academy at the University of Gothenburg approved the protocol using human tissues. All patients had given their informed written consent before participating. Human granulosa/lutein cells were obtained from aspirates of IVF patients. The experiments with hGLC were carried out as described previously ( 11 , 23 ). Briefly, immediately after retrieval of cumulus-oocyte complexes, the remaining cells in aspirates were subjected to Percoll gradient centrifugations to remove red blood cells. The isolated cells were first examined for morphology and counted under the microscope. The cells were then resuspended with OptiMEM media supplemented with 10% fetal bovine serum and antibiotic-antimycotic and then seeded onto culture plates (2.5 × 10 5 cells/ml). The cells were acclimatized for six days, changing media every 24h. At the end of acclimation, the cells were treated with or without hCG (1 IU/ml) in OptiMEM media supplemented with antibiotic-antimycotic for 0h, 12h, and 36h. The cells were collected for Western blot analysis. Only samples showing acceptable viability (above 85%) were subjected to the following processes. 0.5 x 10 6 cells were used to generate single-cell cDNA libraries using the Chromium Single Cell 3’ v3 Reagent kit at the Art & Science Imaging Center of the University of Kentucky. Briefly, individual cells were separated into droplets and barcoded in gel beads, including 10x cell barcodes, unique molecular identifier (UMI), and poly(dT) sequences. Reverse transcription reactions were engaged to generate barcoded full-length cDNA. Single-cell cDNA libraries were sequenced on one lane of an Illumina NovaSeq 6000 S4 flowcell as paired-reads with 28 cycles for read 1, 8 cycles for the index read, 150 cycles for read 2. Basecalling and demultiplexing of raw data were performed with the mkfastq command of the Cell Ranger Single Cell Software Suite (v3.1.0, 10x Genomics) at the Roy J. Carver Biotechnology Center at the University of Illinois. Using default parameters, the scRNA-seq output was gone through mkfastq and count pipelines in the Cell Ranger Single Cell Software Suite (v3.1.0). Reads were quantified using the human genome assembly GRCH38 (hg38). The cellranger count pipeline using the default setting called 2,067, 3,456, 2,452, and 1,173 cells in 4 IVF patient samples, respectively. Mean raw reads were 202,001, 112,324, 160,483, and 370,893 in each sample, and median genes per cells were 1,850, 1,444, 1,935, and 1,527 in each sample. The following analyses were implemented in R (v4.2.0) using the Seurat package (v.4.0.6)( 24 ). Using the Seurat package, feature barcode matrices from 4 IVF patient samples were imported using Read10x() and CreateSeuratObject(). Cells showing a minimum expression of 500 genes were filtered, and genes showing the expression in a minimum of 5 cells were applied. After log-normalization using NormalizeData(), four datasets were integrated following FindIntegrationAnchors() and IntegrateData() procedures with the first 20 principal components for weighting. Filtered cells were dimensionally reduced and visualized on a Uniform Manifold Approximation and Projection (UMAP)( 25 ). Unsupervised clustering analyses on the UMAP embedding were conducted using FindClusters() with ‘resolution = 0.05’, and cells were divided by four individual IVF patient samples in the same UMAP. The cells on the UMAP of 4 patient samples were re-clustered using cell-specific markers genes ( Supplementary Table 3 ). Cell-specific differential gene expression among cell groups determined was identified by FindMarker() or FindAllMarker() and visualized by FeaturePlot(), VlnPlot(), ad DoHeatmap() of the Seurat package and dittoDotPlot() and dittoBarPlot() from the dittoSeq package. Differential gene expression between granulosa cells and leukocyte subpopulations was implemented through a DEseq2 test in FindMarker(). The BAM files from 4 IVF patient datasets were individually transformed and merged into a loom format file by velocyto (v0.17.17) ( 26 ). The Seurat object including only clusters designated as granulosa cells in the samples was converted into a loom format file. These two loom files were combined and used for scVelo (v0.2.3) with default parameters ( 27 ). The velocities are projected onto the UMAP embedding and colored the cells by the cluster assignments exported from the Seurat object. Follicles were fixed in 4% formaldehyde, embedded in paraffin, sectioned (7 μm), and processed for immunostaining. Briefly, heat-induced epitope retrieval was performed in a Biocare Medical Decloaking chamber utilizing Dako’s low pH Target Retrieval Solution. Primary antibody incubation was carried out at 4°C overnight for CD68, IL1B, and MMP9 ( Supplementary Table 4 ). The antibody was detected using an appropriate Immpress alkaline phosphatase kit and Vector Red AP chromogen (Vector Laboratories) according to the manufacturer’s instructions. Slides were counterstained with hematoxylin. The negative control slides were prepared in an identical manner and processed without a primary antibody. Whole-cell lysates were isolated from cultured human granulosa/lutein cells, denatured, run on a 10% polyacrylamide gel, and then transferred onto a nitrocellulose membrane as described previously ( 23 ). The membrane was incubated overnight at 4°C in 5% skim milk/Tris-buffered saline including 0.1% Tween-20 solution containing primary antibodies against IL1R1 or ACTB ( Supplementary Table 4 ). The blots were incubated with the respective secondary HRP-conjugated antibody for 1h. Peroxidase activity was visualized using the Amersham ECL Prime Western Blotting Detection Reagent (Cytiva).

The LH surge triggers the migration of a myriad of leukocytes into the ovulatory follicle ( 4 , 9 , 10 ). The leukocytes recruited in periovulatory follicles have been considered critical mediators of the LH surge-induced changes, facilitating successful ovulation and CL formation in the ovary ( 44 ). Despite the vital roles of leukocytes in these processes, the identity and exact functions of these leukocytes in promoting follicular rupture and luteinization have remained largely elusive. The present study provides comprehensive information on subpopulations and their transcriptomic profiles of leukocytes present in follicular aspirates obtained from women undergoing the standardized IVF procedure due to male factor infertility. In addition to leukocytes, our results also identified multiple clusters of granulosa cells exhibiting varying differentiation status and theca cells in follicular aspirates. More importantly, the single-cell transcriptomic profiling of these heterogeneous cell populations suggests that these cells (e.g., leukocytes, granulosa cells, and theca cells) of ovulatory follicles work together to promote ovulation and luteal transformation in humans. Our study in 1994 first documented the presence of different subsets of leukocytes (macrophages, neutrophils, and T cells) in the follicle wall obtained from normally cycling women ( 45 ). This study also reported the increases in the density of these leukocytes at ovulation ( 45 ). Expanding this report, we revealed the presence of additional subsets of immune cells, including NK cells, NKT cells, B cells, and dendritic cells in follicular aspirates, suggesting that these leukocyte subpopulations are also recruited/activated and play critical roles in human ovulatory follicles. Similar to the original study ( 45 ), the current data identified macrophages as a major leukocyte population (34% of all leukocytes) that were divided into two subsets, M1-macrophage and M2-macrophage, which are typically considered as pro and anti-inflammatory cytokine-producing cells, respectively ( 46 ). Neutrophils and Baso/Eosinophil were also detected. A little over half of the leukocytes were lymphocytes, consisting of helper T, cytotoxic T, NK, NKT, and B cells. Although the smallest number, dendritic cells were also present. Leukocytes are involved in all facets of the inflammatory response from initiation to termination through the coordinated action of different subtypes of leukocytes throughout each phase of inflammation ( 47 ). The current data showed that at least ten different subsets of immune cells are present in the follicular aspirate at 36h after hCG administration, a time point around 2 h before follicle rupture. Therefore, our findings strongly support the generally accepted hypothesis that a restricted inflammatory response is a part of the ovulatory process in which a diverse array of immune cells are recruited. In support of this concept, our data also revealed that these leukocytes express a variety of inflammatory factors that can either directly regulate follicular cell functions by activating their receptors in granulosa and theca cells or stimulate tissue remodeling and angiogenesis required for follicle wall rupture and CL formation as discussed below. Previous studies showed that IL1B and TNF were present in the follicular fluid from IVF patients but not in spent media from granulosa-lutein cell cultures ( 45 , 48 ), suggesting that these cytokines are leukocyte-derived. Indeed, our data confirmed that IL1B and TNF expression was exclusive in specific subsets of leukocytes, mainly in macrophages. Meanwhile, their respective receptors, IL1R1 and TNFRSF1A, were detected in follicular cells and leukocytes. In support of our findings, previous microarray data listed IL1R1 and TNFRSF1A as differentially up-regulated genes in granulosa cells obtained 36 h after hCG administration compared to those obtained before the endogenous LH surge/hCG administration ( 36 , 49 ). In cultured human granulosa-lutein cells, IL1B increased progesterone (P4) production ( 50 , 51 ). TNF treatment also increased granulosa-lutein cell proliferation, estradiol (E2), and P4 production in vitro ( 52 ). Importantly, in animal models, both IL1B and TNF were shown to increase ovulation rate ( 3 , 21 , 53 , 54 ), indicating that these leukocyte-derived factors promote the ovulatory process. In addition to the cytokines mentioned above, we found that leukocytes express various paracrine factors that can act directly on follicular cells. Among those included are various EGF-like peptides (AREG, EREG, BTC, and NRG1) and PLAU. Previous studies have shown that EGF-like peptides are highly up-regulated in granulosa cells of ovulatory follicles and play vital roles in ovulation, oocyte maturation, and/or cumulus cell expansion ( 55 - 57 ). Besides these secretory proteins, the present study unveiled the expression of prostaglandin synthetase 2 (PTGS2), a rate-limiting enzyme in PG synthesis, in specific subsets of leukocytes. Previous studies from mice and monkeys have demonstrated that PTGS2 expression in ovulatory follicles is critical for successful ovulation ( 58 , 59 ). Subsequent studies showed that PGE2, through its receptors, functions as an essential mediator of ovulation ( 60 ). The current findings showed the expression of these genes in leukocytes and their respective receptors in follicular cells (e.g., EGFR, PTGER3, ERBB3, ERBB4, and PLAUR). It is also worth noting that PTGER2, PLAUR, and ERBB2 are expressed in leukocytes' subpopulations, indicating that these paracrine factors can act back on leukocytes. Together, these data indicated leukocytes as additional sources of these paracrine factors that can directly modulate the function of follicular cells as well as leukocytes, implicating their involvement in various aspects of the periovulatory process. The extracellular matrices (ECM) around ovulatory follicles undergo rapid and dramatic remodeling during ovulation and CL formation. The ECM remodeling requires the coordinated action of proteases and protease inhibitors ( 13 ). Previous studies have mainly focused on identifying the expression of various proteases and their inhibitors in follicular cells ( 4 ). However, it has also been speculated that leukocytes contribute to ECM remodeling by secreting these proteins. Indeed, not only did the current data confirm the expression of several members of MMP, ADAM, and ADAMTS families as well as TIMPs in follicular cells, but it also showed the expression of these proteins in various types of leukocytes. Many of these proteins were present in both follicular cells and leukocytes, but some were expressed predominantly either in follicular cells (e.g., MMP2, MMP10, MMP19 , and ADAMTS1 ) or leukocytes (e.g., MMP9, MMP17, MMP25, ADAM19 , and ADAM28 ). The current approach also allowed us to pinpoint exactly where each protease is expressed in a cell type-specific manner. For instance, our data showed that MMP2 is expressed in theca cells and M2-macrophages, while MMP10 expression is confined to granulosa cells. In contrast, TIMP1 expression was detected in all cell types present in follicular aspirates, while TIMP4 expression was not detectable in any cell types. Together, these data suggest that ovulatory tissue remodeling results from the collaborative actions of specific proteases and their inhibitors from both follicular cells and leukocytes. Along with massive tissue remodeling, rapid blood vessel generation (angiogenesis) also occurs in ovulatory follicles and newly forming CL. Leukocytes were shown to be involved in angiogenesis by secreting VEGFs in other systems (( 4 )( 4 )4). In the ovary, the ovulatory hCG-induced increases in the expression of VEGFA and PGF in follicular cells of monkey preovulatory follicles ( 37 ). Moreover, the follicular administration of antibodies for VEGFA and PGF resulted in compromised ovulation and CL formation in monkeys ( 41 , 42 ). Our data showed the expression of mRNA for VEGFA, VEGFB, VEGFC , and PGF in both granulosa cells and theca cells. In addition, VEGFA and VEGFB mRNA were detected in various types of leukocytes, identifying additional sources of VEGFs that can play a role in periovulatory angiogenesis in humans. Besides leukocytes, this study showed the presence of five different clusters of granulosa cells with differential gene expression profiles. RNA velocity and ovulatory marker gene expression analyses indicated the differentiation projection of these granulosa cell clusters from early (GC 1 and 2) to late ovulatory phases (GC 3, 4, and 5). These data suggested the temporal heterogeneity of the granulosa cell population. One possible explanation for this temporal heterogeneity might be that granulosa cells respond to ovulatory hCG stimulation with different rates depending on their location within the preovulatory follicle. This notion is based on the findings that LH/hCG receptors (LHCGR) are mainly expressed in the outer mural granulosa cell layer lining the basement membrane, but absent in the inner mural granulosa cells facing the antrum site ( 61 ). Therefore, it is conceivable that the LH surge/hCG administration first activates its receptor (LHCGR) expressed on the basement membrane side. Then, downstream signaling molecules (e.g., cAMP) generated by LHCGR-activated cells propagate from the outer mural granulosa cell layer to inner mural granulosa cells through gap junctions. This transition likely generates the temporal response gradient in ovulatory gene expression among granulosa cells within the follicle. In support of this notion, scRNA-seq analysis by Wu and colleagues ( 62 ) also reported nine different clusters of granulosa cells in follicular aspirates collected from a single follicle from each patient (n = 6 patients). However, at present, it is difficult to predict exactly where granulosa cells with the different expression profile are distributed within the ovulatory follicle. Therefore, further studies are needed to understand the granulosa cell heterogeneity in human ovulatory follicles. Besides granulosa cells, theca cells were also detected as one of the clusters, although they constituted the smallest number (0.3% of total cells). In summary, the current study provided invaluable information on leukocyte subpopulations and their gene expression profiles in IVF patient samples. From these data, we hypothesized that several subtypes of leukocytes recruited to ovulatory follicles secrete a diverse array of paracrine factors that can directly impact granulosa cell function by activating their receptors (e.g., IL1B, IL16, TNF, TNFSF10) or control the tissue remodeling (e.g., MMPs, ADAMs, ADAMTSs, TIMPs) and angiogenesis (e.g., VEGFA, VEGFB, THBS1), thereby promoting ovulation and luteinization ( Fig. 7 ). Given that changes in numbers, actions, and functions of immune cells in tissues are directly related to chronic or abnormal inflammation ( 63 , 64 ), the current datasets can also serve as a foundation for understanding chronic or abnormal inflammation-associated ovarian etiologies, including polycystic ovary syndrome, ovarian endometriosis, and ovarian cancer ( 65 - 67 ). Therefore, this study may contribute to designing novel strategies or treatments to improve female fertility.

Leukocytes in the ovary are essential in situ regulators involved in various aspects of ovarian function ( 1 ). They exert their actions through local secretion of paracrine factors, such as various cytokines that have been implicated to be important in steroidogenesis and oocyte maturation ( 2 - 5 ). A recent study showed significant changes in peripheral blood leukocyte profiles across the menstrual cycle ( 6 ). The total leukocyte count and the count of neutrophils were increased around ovulation compared to the day of menstruation ( 6 ). The mid-cycle luteinizing hormone (LH) surge induces dramatic morphological and physiological changes in and around preovulatory follicles, culminating in the rupture of the follicle wall, extrusion of a cumulus-oocyte complex, and rapid transformation of the ruptured follicle into the corpus luteum (CL). Since Espey proposed the hypothesis ( 7 ) that ovulation is an inflammatory reaction, evidence has accumulated to show the similarities between the ovulatory process and an acute inflammatory response ( 4 , 8 ). One of the critical features of the inflammatory response is the infiltration of leukocytes into the site of infection or injury. In the ovary, the LH surge increased the ovarian blood volume ( 9 ) and a rapid influx of millions of leukocytes into the ovulatory follicles ( 8 , 10 ). It has been postulated that the LH surge stimulates follicular cells of preovulatory follicles to secrete chemokines, which in turn induce the migration of leukocytes from the circulation into the ovary. In the human ovary, our previous studies reported that granulosa cells of preovulatory follicles express chemokines [e.g., CCL20 and CXCL12]. These chemokines potentially recruit leukocytes expressing CCR6 (receptor for CCL20; T and dendritic cells) and CXCR4 (receptor for CXCL12; T cells and other leukocytes) to migrate into the follicle ( 10 - 12 ). Leukocytes recruited to the ovulatory follicle are suggested to be involved in follicle wall rupture and angiogenesis during the periovulatory period ( 10 ) by producing and secreting tissue proteases, such as metalloproteinases (MMPs)( 13 - 15 ), and vascular endothelial growth factor (VEGF)( 16 , 17 ), respectively. In addition, a paracrine interaction between leukocytes and granulosa cells stimulated granulosa cell secretion of granulocyte-colony stimulating factor (G-CSF), known to induce neutrophil differentiation and activation ( 18 ). Previous studies using animal models provided experimental evidence supporting the functional significance of leukocytes in the ovulatory process ( 19 - 22 ). The perfusion of leukocytes into rat ovaries showed a 3-fold increase in the ovulation rate ( 19 ). On the contrary, the depletion of macrophages by intrabursal injection of clodronate liposome, which induces macrophages’ apoptosis, blocked ovulation ( 20 ). Neutralizing antibodies for either neutrophils ( 21 ) or their activating factor, interleukin-8 [CXCL8 ( 22 )] showed reduced ovulation rates by 27% and 50.5%, respectively. Despite their vital roles in the ovulatory process, the types of leukocytes recruited to the ovary and their specific functions in ovulation and CL formation remain largely elusive, especially in humans. Conducting such studies has been challenging due to; 1) the lack of technical approaches to simultaneously characterize and identify different subtypes of leukocytes recruited into the periovulatory follicles and 2) extremely limited access to periovulatory follicles in the human ovary from normally cycling women. Recent advances in single-cell RNA sequencing (scRNA-seq) technology allow us to overcome the first limitation: it identifies cell subpopulations in a heterogeneous cell pool and individual cells’ gene expression profiles. By performing scRNA-seq analyses on follicular aspirate samples obtained from women with no-ovarian infertility etiologies undergoing the IVF procedure, the present study determined the subpopulation of leukocytes and identified the transcriptomic profile in each subset of leukocytes and follicular cells collected at 36 h after hCG administration. Next, the transcriptomic profile data were further analyzed to identify leukocyte-secreted factors that 1) can directly communicate with follicular cells and 2) are involved in tissue remodeling and angiogenesis. Importantly, to relate these findings to the normal ovulatory process, we examined the localization of leukocytes and leukocyte-specific expression of CD68, IL1B , and MMP9 in dominant follicle tissues obtained throughout the periovulatory period from normally cycling women.