Whole exome sequencing reveals novel variants associated with diminished ovarian reserve in young women

Background: Diminished ovarian reserve is one of the most important causes of female infertility. In the etiology study of DOR, besides age, it is known that chromosomal abnormality, radiotherapy, chemotherapy and ovarian surgery can result in DOR. For young women without obvious risk factors, gene mutation should be considered as a possible cause. However, the specific molecular mechanism of DOR has not been fully elucidated.

Methods: In order to explore the pathogenic variants related to DOR, twenty young women under 35 years old affected by DOR without definite factors damaging ovarian reserve were recruited as the research subjects, and five women with normal ovarian reserve were recruited as the control group. Whole exome sequencing was applied as the genomics research tool.

Results: As a result, we obtained a set of mutated genes that may be related to DOR, where the missense variant on GPR84 was selected for further study. It is found that GPR84^Y370H variant promotes the expression of proinflammatory cytokines (TNF-α, IL12B, IL-1β) and chemokines (CCL2, CCL5), as well as the activation of NF-κB signaling pathway.

Conclusion: In conclusion, GPR84^Y370H variant was identified though analysis for WES results of 20 DOR patients. The deleterious variant of GPR84 could be the potential molecular mechanism of non-age-related pathological DOR through its role in promoting inflammation. The findings of this study can be used as a preliminary research basis for the development of early molecular diagnosis and treatment target selection of DOR.

GRAPHICAL ABSTRACT

Introduction

Ovarian reserve, refers to the remaining follicles or oocytes in the ovaries, which is not infinite for a woman from birth. Diminished ovarian reserve, commonly known as DOR, which is characterized by reduction of follicle pool, sometimes along with menstrual abnormalities, ovulation dysfunction and reduced organ, and systemic function owing to low estrogen levels in severe cases. The incidence of DOR ranges from 6% to 64% among infertile women of different ages, which is increasing in recent years (Hu et al., 2020). The assessment of one’s ovarian reserve includes antral follicle count (AFC), serum level of anti-Müllerian hormone (AMH), follicle-stimulating hormone (FSH), inhibin B and so on. In 2020, the Practice Committee of the American Society for Reproductive Medicine (ASRM) suggested that currently AMH and AFC would be the most sensitive and reliable markers of ovarian reserve (Practice Committee of the American Society for Reproductive Medicine. Electronic address and Practice Committee of the American Society for Reproductive Medicine, 2020). The diagnostic criteria of DOR are not uniform among different clinics across the world. According to the Bologna criteria suggested by ESHRE in 2011, DOR is diagnosed when serum level of AMH <0.5–1.1 ng/mL and AFC < 5–7, which has been the reference of all criteria since then.

DOR is an important cause of female infertility since the depletion of follicle is reported to be associated with ovulation disorder. In assisted reproductive technology for infertile women, DOR indicates poor prognosis because it may result in poor ovarian response (POR) to ovarian stimulation during IVF/ICSI-ET. The pathological mechanism of DOR can be divided into two types: a smaller initial follicle pool and accelerated depletion of follicles (Jaswa et al., 2021). It is well established that ovarian reserve decreases with age. Over the total lifetime, a woman’s ovarian reserve will certainly meet the criteria for DOR when approaching menopause (Greene et al., 2014). It is called “Age-related DOR” when age is the sole factor contributing to DOR (Liu et al., 2021). Other than this physiological condition, there are various etiology factors that accelerate the depletion of follicles, affect ovarian reserve, and thus bring about pathologic DOR, including genetic abnormity, diseases causing ovarian injury such as endometriosis and pelvic infection, iatrogenic factors such as ovarian surgery, chemotherapy and pelvic irradiation, environmental and lifestyle factors such as smoking (Collaborators, 2021).

Among all the causes, genetic factors appear to be important and have not been fully elucidated. So far it has been reported that chromosome abnormality (e.g., 45, X mosaicism), gene mutation [e.g., FMR1 premutation (Tang and Yu, 2020)], polymorphisms [GDF9 (Wang et al., 2013), FSHR (Livshyts et al., 2009), and ESR1 (Livshyts et al., 2013)], differential gene expression [e.g., AMH(Skiadas et al., 2012) and epigenetic effects (Nilsson et al., 2012)] were related to the disease. In many cases, DOR is a multifactorial disease, but some idiopathic cases occurring in young women may have an underlying genetic component (Greene et al., 2014). It is of great significance to explore the genetic etiology of DOR for the early identification and intervention. Novel sequencing technologies have been substantial for expediting our understanding of the genetics of diseases. Over the past decade, the emerging next-generation sequencing technology (NGS) has shown its promising and efficient application in genetic study of human reproduction diseases. In order to explore the underlying genetic factors involved in DOR, whole exome sequencing (WES) is an ideal approach since it only focuses on the exons or protein coding regions in the genome and is more cost-effective than whole genome sequencing (WGS).

In the sample selection, most of the previous studies used peripheral blood as the sample source (Wang et al., 2019; Liu et al., 2020; Rouen et al., 2022), which could only detect genetic variations, while somatic mutations will not be detected. The most ideal experimental sample in our study would be oocytes, but the oocytes of patients are so precious that they must be reserved for subsequent in vitro fertilization and cannot be used as sequencing samples. In the ovaries, every oocyte and the surrounding granulosa and theca cells form the functional unit-the follicle, which provides a protective and suitable environment for the oocyte. Ovarian granulosa cells (GCs) are situated closely adjacent to the oocyte and there are plentiful and various material exchange and signal transduction between them. Furthermore, GCs are the synthesis sites of vital hormones including estrogen and AMH, and thus play a key role in folliculogenesis and ovarian function (Strauss and Williams, 2019). Many studies have suggested that the apoptosis or dysfunction of GCs has a negative effect on the quantity or quality of oocytes, thus resulting in ovarian insufficiency. Therefore, mural GCs in follicular fluid which are readily available in IVF/ICSI procedure may be an appropriate choice for studying the genetic abnormalities especially somatic mutation associated with DOR. Somatic mutations occurring in GCs will have a great impact on ovarian reserve, so we proposed to use ovarian GCs as the sequencing sample in order to find out more significant variants.

In order to identify deleterious variants occurring in young DOR patients without other obvious pathogenic factors, so as to provide clues and basis for the early identification and risk prediction of DOR, we enrolled 20 DOR patients who underwent IVF/ICSI in the reproductive center of Tianjin Central Hospital of Gynecology Obstetrics, and performed whole exome sequencing of the ovarian granulosa cells collected from the patients to determine causal genetic variances of DOR.

Materials and methods

Ethical approval

The research was approved by the reproductive ethics committee of Tianjin Central Hospital of Gynecology Obstetrics (No. ZY2021001, 06/01/2021), and the study was conducted in accordance with the declaration of Helsinki. All patients had given their informed consent for the use of their samples for future researches.

Human subjects

From May 2018 to September 2022, of all age groups there were 382 DOR patients without definite factors impacting ovarian reserve, out of 8,715 patients undergoing IVF/ICSI-ET in the reproductive center of Tianjin Central Hospital of Gynecology Obstetrics. From these females, twenty women affected by DOR and five women with normal ovarian reserve were recruited in the study. The clinical assessment and definition of DOR in our study were based on the Bologna criteria. DOR was diagnosed when the serum level of AMH ≤1.1 ng/mL and AFC < 5-7, while matched women with normal ovarian reserve were included into the control group. All the included cases were younger than 35 years old and had a normal karyotype of 46, XX. Females with FMR1 premutation were excluded, as were patients with PCOS, history of ovarian surgery, pelvic infections, chemotherapy, pelvic radiation, autoimmune disease, and other severe medical conditions.

Sample preparation

Primary ovarian GCs were isolated refer to red blood cell lysis protocol for RBC lysis as described below. Patients’ follicular fluid collected on the day of oocyte pick-up during the IVF/ICSI procedure was centrifuged at 450 × g for 10 min at room temperature. Decant the supernatant and resuspend the cell pellet in 3 mL of PBS. Add 30 mL of red blood cell lysis buffer (Solarbio Science & Technology, Beijing, China), and then invert to mix. Incubate on ice for 15 min, during which invert gently every 5 min. Centrifuge at 450 × g for 6 min at 4°C and decant the supernatant. Resuspend the pellet gently in 20 mL of red blood cell lysis buffer and pellet GCs by centrifugation at 450 × g for 6 min at 4°C.

Genomic DNA extraction, isolation and whole exome sequencing

Total genomic DNA was extracted and isolated by conventional phenol-chloroform extraction. Before sequencing, the concentration of DNA was measured by Qubit Fluorometer and the integrity and purification were detected by Agarose Gel Electrophoresis. Library preparation and sequencing were carried out following certified protocols of the Beijing Genomics Institute (BGI). Briefly, the genomic DNA was randomly fragmented and adapters were added to both ends of the fragments followed by PCR amplification. The qualified library from each sample was respectively loaded on BGISEQ-500 sequencing platforms to perform high-throughput sequencing. The raw image files produced by the sequencer were processed using BGISEQ-500 base caller with default parameters, and the sequencing data were stored in FASTQ files.

Filtering, mapping, and local realignment

The raw reads were preprocessed as follows: 1) The reads containing adapters were removed; 2) The reads containing more than 50% low-quality bases (base quality ≤5) were filtered; 3) The reads containing more than 10% N bases were removed for the downstream analysis. The high-quality reads were mapped to the human reference genome (GRCh38) using Burrows-Wheeler Aligner (BWA v0.7.15) (Li and Durbin, 2010) and duplicates were marked and discarded via Picard-tools (v2.5.0). Then, the reads were subjected to local realignment around insertions/deletions (indels) and base quality score recalibration (BQSR) using the Genome Analysis Toolkit (GATK v3.7) (McKenna et al., 2010) following the best practice pipeline (DePristo et al., 2011).

Variant calling, annotation, and filtering

Single-nucleotide variants (SNVs) and indels were called using GATK3.7 (HaplotypeCaller) and annotated with SnpEff software (Cingolani et al., 2012). To filter irrelevant variants, we selected only variants whose predicted impact on the protein were “high” or “moderate”. The variants with minor allele frequency (MAF) ≥ 5% in 1,000 Genomes database were removed. The pathogenicity scores of the variants were calculated using various packages, including SIFT (Ng and Henikoff, 2003), Mutation assessor (Schwarz et al., 2010), PolyPhen2 (Adzhubei et al., 2013), and MetaSVM(Zaucha et al., 2020), and the variants with SIFT score ≤0.05 or PolyPhen2 ≥ 0.909 or Mutation Assessor score ≥1.9 or MetaSVM = “deleterious” were selected for downstream analysis.

Sanger verification

Approximately 400–800 bp DNA fragments containing the selected variant site were amplified by PCR using high fidelity DNA polymerase KOD-Plus (Toyobo, Osaka, Japan) and the primers in Supplementary Table S3 according to the manual. PCR products were purified by 1% agarose gel and gel Midi Purification Kit (Tiangen, Beijing, China) and sequenced by Sanger sequencing on ABI 3730XL sequencers at Sangon Biotech (Shanghai, China). DNA sequences were aligned with NCBI by BLASTN programs and Chromas software.

Cell culture

KGN cells were purchased from Procell Life Science&Technology (Wuhan, China). Cells were cultured in DMEM/F12 medium (Biological Industries, USA) supplemented with 10% FBS (Biological Industries) and 100 IU/mL penicillin/streptomycin (Gibco, Grand Island, NY) in a humidified atmosphere of 5% CO₂ at 37°C.

Construction of stable cell lines

The full-length coding sequence and mutated sequence of human GPR84 was amplified by PCR using GCs’ cDNA as template from normal females. The sequences were respectively cloned into plasmid pLV-EF1α-MCS-IRES-Bsd (Biosettia). All plasmids were validated by Sanger sequencing. The procedures of lentivirus vectors package in HEK293 cells and infection were described previously (Chang et al., 2015). Cells infected with lentivirus were screened with 2.5 μg/mL of blasticidin (Thermo Fisher Scientific, Waltham, MA).

RNA extraction and real-time quantitative PCR (qRT-PCR)

Total RNA was extracted with Trizol reagent (Thermo Fisher Scientific) according to the protocol, and reversely transcribed with Hifair^® II 1st Strand cDNA Synthesis Kit (Yeasen Biotech, Shanghai, China). Real-time quantitative polymerase chain reaction (qRT-PCR) analysis was performed using SYBR Green SuperMix (Yeasen Biotech) on LightCycler96 system (Roche, Basel, Switzerland). The primer sequences are listed in Supplementary Table S4. The relative gene expression was calculated by the 2^−ΔΔCt method with β-actin as a loading control.

Protein extraction and Western blot

Cells were lysed in RIPA lysis buffer with protease inhibitor cocktail (Roche) and phosphatase inhibitor cocktail (Sigma-Aldrich, St Louis, MO). The protein concentration was measured by BCA Protein Assay Kit (Thermo Fisher Scientific). The same amount of protein was loaded into SDS-PAGE gels for electrophoresis, and transferred onto a PVDF membrane (Merck, Darmstadt, Germany). After blocking with 5% defatted milk for 1 h at room temperature, the membrane was incubated with primary antibodies: GPR84 (1:2500, Invitrogen, Waltham, MA, USA), NF-κB p65 (1:1,000, CST, Danvers, MA, USA), β-actin (1:5000, Affinity Biosciences), α-tubulin (1:1,000, CST), Lamin A/C (1:2500, Proteintech, Wuhan, China) at 4°C overnight. The HRP-conjugated secondary antibodies (goat anti-rabbit IgG, goat anti-mouse IgG (1:2500, ZSGB-BIO, Beijing, China) were using for incubation of the membranes for 1 h at room temperature. The signals were visualized with an ECL detection system and photographed by Tanon-5200 (Tanon, Shanghai, China). The relative protein level was quantified with Image J 1.48v (NIH, Bethesda, MD).

Statistical analysis

Statistical analysis was performed in Prism 8.0.1 (GraphPad Software, San Diego, CA, USA). Quantitative data were presented as means ± SEM, and the differences among the groups were analyzed by Student’s two-tailed t-test. Differences are considered statistically significant at *p < 0.05; **p < 0.01; ***p < 0.001 and ****p < 0.0001; ns means no significance.

Results

Clinical information of included subjects

The baseline characteristics of included women were summarized in Supplementary Tables S1, S2. There were no significant differences in baseline characteristics such as age, BMI, and infertility period between the DOR and control group.

Landscape of variations in DOR

To explore the landscape of variations related to DOR, we performed WES of 20 DOR cases and 5 controls who underwent IVF/ICSI in the reproductive center of Tianjin Central Hospital of Gynecology Obstetrics. On average, 231 million reads passed quality control for each sample, and the mean sequencing depth on target regions was 235.89×. On average per sequencing individual, 98.86% of the targeted bases had at least 10× coverage. To filter irrelevant variations, we designed a systemic variation selection framework as shown in Figure 1. Briefly, an average number of 370,759 variants passed the GATK quality control in each individual, including 313,744 SNVs and 57,015 indels. Firstly, the variants overlapped with those identified in the normal samples were removed, and 127,310 variants left, after which the 87,634 variants with MAF <5% in the 1,000 Genomes database were chosen for the downstream analysis. Next, the non-synonymous SNVs and frameshift indels were selected followed by further choosing variants identified as deleterious using various pathogenicity prediction software, including SIFT ≤0.05, Polyphen2 ≥ 0.909, Mutation assessor ≥1.9 and MetaSVM = “Deleterious”.

FIGURE 1

FIGURE 1. Framework of WES data analysis. Of the qualified variants, the variants 1) identified in normal samples; 2) MAF >5% in the 1,000 Genome database; 3) in the non-coding regions; 4) predicted as no impact were removed, and top 20 mutated variants were finally reported.

As a result, 822 variants involving 731 genes were identified related to DOR, including 481 (58.52%) SNVs and 341 (41.48%) indels (Figure 2A). Of 481 SNVs (Figure 2B), 258 were C>T (49.33%), 89 were T>C (17.02%), 61 were C>G (11.66%), 52 were C>A (9.94%), 36 were T>G (6.88%), 27 were T>A (5.16%). In total, the variants included 466 missense variants, 15 SNVs at splice site, 207 frameshift deletion, 102 frameshift insertion, and 32 indels at splice site (Figure 2C). On average for each individual, there were 25 missense variants, 15 frameshift deletions, 8 frameshift insertions and 3 variants at splice site (Figure 2D). The distribution of different types of variants in each patient was shown in Figure 2E. High-throughput technologies, such as WES, can generate a large list of genes of interest. However, it is difficult to extract biological meaning from the list due to the noisy nature of biological processes. To explore the underlying relevance among these mutated genes, we conducted pathway enrichment analysis and functional annotation through KEGG database (Kanehisa et al., 2019), Reactome database (Gillespie et al., 2022) and Gene Ontology (GO) database (Ashburner et al., 2000; Gene Ontology, 2021) (Figure 2F). It was found that the mutated genes were enriched in KEGG pathways including protein digestion and absorption, metabolic pathways, beta-Alanine metabolism, ECM-receptor interaction, and cholesterol metabolism; the genes were enriched in Reactome pathways including metabolism, collagen biosynthesis and modifying enzymes, collagen formation, collagen chain trimerization, and ECM organization. The GO analysis revealed the mutated genes may be involved in calcium ion binding, collagen fibril organization, ECM structural constituent conferring tensile strength, calmodulin binding, and Z disc.

FIGURE 2

FIGURE 2. Summary of whole exome sequencing results of 20 DOR patients. (A) The percentage of SNVs and indels in filtered variants; (B) The classes of SNVs and their frequencies; (C) The frequencies of different types of variants in total samples; (D) The violin plot of the number of different types of variants in each patient; (E) The distribution of different types of variants in each patient; (F) Functional analysis of 731 genes using KEGG, Reactome, and GO database.

Top mutated variants in DOR

To identify reliable variants related to DOR, the 731 mutated genes were sorted according to the variation frequency in 20 samples, and the top 20 mutated genes were finally reported, which were GIGYF2, ODF1, DENND4B, ANKRD36, EFCAB2, LRP8, ZAN, ANKLE1, KAZALD1, MUC19, C16orf52, CNTN5, FAM228B, GPR84, KIAA0040, OR5H15, RETSAT, ZNF527, GOLGA6L4, and NPIPB11 (Figure 3A). GIGYF2 has the highest prevalence of amino acids altering variants identified in our DOR cohort with frameshift insertions of GC in the same genomic position of mutated samples. KAZALD1 and GRP84 only contained missense variants; ODF1, DENND4B, LRP8, ZNF527 and NPIPB11 had multiple variants in the genes; ANKRD36, MUC19, CNTN5, KIAA0040, RETSAT and GOLGA6L4 had mixed types of variants in different samples. The details of the identified variants in top 20 genes were listed in Table 1. To explore the underlying relevance of top 20 mutated genes, we conducted pathway enrichment analysis and functional annotation (Figure 3B). It was found that the mutated genes were enriched in Reactome pathways including retinoid metabolism and transport, metabolism of fat-soluble vitamins, and sensory perception. Then, we investigated the interactions among these mutated genes using GENEMANIA database (Warde-Farley et al., 2010), as shown in Figure 3C. There were connections such as co-expression or protein-protein interaction among these genes which were reported before, although not clearly elucidated yet.

FIGURE 3

FIGURE 3. Top 20 mutated genes in DOR. (A) The variant characteristics of top 20 mutated genes in each patient. The twenty DOR patients included for WES were represented by “P01”-“P20”. (B) Enrichment analysis of top 20 mutated genes using Reactome pathways. (C) Protein-protein interactions of mutated genes identified using GENEMANIA database.

TABLE 1

TABLE 1. The variants of top 20 mutated genes.

GPR84^Y370H variant promotes the inflammatory response in human ovarian granulosa cells

According to the Sanger sequencing results and biological functions of the top mutated genes, the missense variant (NM_020370.2:p.Tyr370His/c.1108T>C) in GPR84 was selected as the target for further study. The characteristics of GPR84 variant were summarized in Figure 4. The missense variants detected in GPR84 by Sanger sequencing were exactly the same as the WES results (Figure 4A). To rule out polymorphism or genetic instability at the variant site, we performed amino acid sequence alignment of GPR84 among different species. As shown in Figure 4B, tyrosine at position 370 is 100% conserved across all species, which suggested that the variant occurred at this site may be pathogenic. In order to see the expression of GPR84 in ovaries which has not been reported in databases or researches, we detected the mRNA and protein level of GPR84 in different tissue and organs in BALB/c mouse. The results showed that GPR84 was expressed in the ovary to some extent (Figure 4C), so we believe that it can be further studied as a research target for the etiology of DOR. To determine the exact location and feature of the variant site, we predicted the protein structure of GPR84 by RoseTTAFold (Baek et al., 2021). The results showed that the variant site was located within the cell membrane, close to the intracellular site. The electrically neutral, hydrophobic amino acid tyrosine of a six-membered ring mutates to the positively charged, hydrophilic histidine of a five-membered ring (Figure 4D).

FIGURE 4

FIGURE 4. Characteristics of the GPR84 variant. (A) Sanger sequencing validation of the heterozygous c.1108T > C variant in the GPR84 gene. The red box indicates the variant site. (B) Amino acid sequence alignment of GPR84 among different species. The blue box indicates the mutated amino acid. Tyrosine at position 370 is 100% conserved across all species. (C) The expression of GPR84 in different tissue and organs in BALB/c mouse were measured by qRT-PCR and Western blot. (D) Location of the mutated amino acid in the GPR84 protein predicted by RoseTTAFold. The extracellular segments are shown in cyan; the cytoplasmic segments are shown in violet; the transmembrane segments are shown in yellow. The wildtype amino acid tyrosine is shown in red, and the mutated amino acid histidine is shown in blue.

To further explore the effect of GPR84^Y370H variant on GCs, we constructed the GPR84^Y370H overexpressed (shown as GPR84^WT/mut) and GPR84 overexpressed (shown as GPR84^WT/WT) stable cell lines in KGNs (human granulosa-like tumor cell line), while the GPR84^WT/WT group was regarded as the control (Figures 5A, B). GRP84 was reported as an enhancer of inflammatory signaling (Recio et al., 2018b). To study the influence of GPR84^Y370H variant on its proinflammatory effect, KGNs were pretreated with LPS (0.1 μg/mL) for 2 h to set up an inflammation environment before stimulation with GPR84’s biological ligand lauric acid (200 μM) for 24 h. The mRNA level of inflammatory cytokines (TNF-α, IL-6, IL12B, IL-1β) and chemokines (CCL2, CCL5, CXCL1) were measured by qRT-PCR. As shown in Figure 5C, GPR84^Y370H variant promotes the expression of TNF-α, IL12B, IL-1β, CCL2 and CCL5. In order to explore the effect of GPR84^Y370H variant on the activation of inflammation-associated NF-κB signaling pathway, we conducted nuclear and cytoplasmic protein extraction on the cells that have been treated as described above, and tested NF-κB subunit p65 in nucleus and cytosol by Western blot. The results showed that the protein level of p65 in nucleus were increased in GPR84^WT/mut group compared with GPR84^WT/WT group (Figure 5D), which indicated that GPR84^Y370H variant promoted the activation of NF-κB.

FIGURE 5

FIGURE 5. GPR84^Y370H variant promotes the inflammatory response in human ovarian granulosa cells. (A, B) The construction of GPR84 (shown as GPR84^WT/WT) and GPR84^Y370H (shown as GPR84 ^WT/mut) overexpressed KGN cells (human granulosa-like tumor cell line) were verified by qRT-PCR and Western blot respectively. (C, D) The GPR84^WT/WT and GPR84^WT/mut overexpressed KGN cells were treated with LPS (0.1 μg/mL) for 2 h before stimulation with lauric acid (200 μM) for 24 h. (C) The mRNA level of inflammatory cytokines (TNF-α, IL-6, IL12B, IL-1β) and chemokines (CCL2, CCL5, CXCL1) were measured by qRT-PCR. (D) The protein level of p65 in nucleus were increased in GPR84^WT/mut group compared with GPR84^WT/WT group. Lamin A/C and α-Tubulin were used as loading controls in nucleus and cytosol respectively. Data are represented as mean ± SEM. Statistical significance was assessed by two-tailed unpaired t-test. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.

Discussion

In addition to age, factors such as history of ovarian surgery (Zhao et al., 2020), radiotherapy, chemotherapy (Fenton, 2015), and endometriosis (Muzii et al., 2018; Ashrafi et al., 2019) can lead to the decline of female ovarian reserve. However, there still are many young females who have no definite factors damaging ovarian function also suffer from diminished ovarian reserve, resulting in infertility, poor ovarian response in IVF, adverse pregnancy outcomes and so on. DOR in young women is easily neglected, which can bring about serious and irreversible impacts on female fertility. Therefore, it is of great significance to screen out susceptible population of DOR at the early stage so that they can get guidance on fertility timing and necessary treatment as soon as possible. Premature ovarian insufficiency (POI) is the loss of normal ovarian function before 40 years old (Fenton, 2015), which has a similar pathological basis with DOR. It has been speculated that the two diseases share some etiological mechanism in common (Zhao et al., 2019). With the rapid development of next-generation sequencing technology, the association of POI etiology with genetic mutations has been supported by many researches (Jiao et al., 2018), and numerous potential pathogenic mutations have been found (Qin et al., 2015; Jaillard et al., 2020; Chon et al., 2021). However, there are still few studies on the genetic molecular etiology of DOR.

Here, we performed whole exome sequencing in 20 young DOR patients and 5 women with normal ovarian reserve to explore the potential genetic mechanism of DOR. In the sample selection, most of the previous studies used peripheral blood as the sample source (Wang et al., 2019; Liu et al., 2020; Rouen et al., 2022), which could detect the genomic variation, but the somatic mutation occurred locally in the ovary will not be detected. The most ideal experimental sample should be oocytes, but the oocytes of patients are so precious that they must be reserved for subsequent in vitro fertilization and cannot be used as sequencing samples. It is confirmed that GCs play a key role in oocyte development, follicular growth and maturation (Strauss and Williams, 2019). Somatic mutations occur in GCs will have a great impact on ovarian reserve, so we proposed to use ovarian GCs as the sequencing sample in order to find out more significant variants.

Though a series of filtering steps for WES data, we identified 822 functional deleterious variants involving 731 genes. Pathway enrichment analysis and GO functional annotation showed that most of the genes were enriched in pathway such as the biosynthesis and metabolism of ECM, collagen, calcium, and lipid, which are associated with ovarian function. The variants in the top 20 genes with highest variant frequencies among the 20 samples were verified by Sanger sequencing. It is confirmed that the reported variants do occur in ODF1, EFCAB2, ZAN, KAZALD1 and GPR84. Through literatures reviewing, it is found that most of these genes had not been reported as related to ovarian reserve previously. By studying their basic biological functions and verifying the variants though Sanger sequencing, the missense variant in the gene GPR84 was selected as the target for further study.

G protein-coupled receptor 84 (GPR84), also known as inflammation-related G protein-coupled receptor EX33, is a member of metabolite-sensing G protein-coupled receptors (GPCRs) family with medium-chain fatty acid (MCFA) as its natural ligand (Recio et al., 2018a; Luscombe et al., 2020). GPR84 functions as an enhancer of inflammatory signaling in macrophages once inflammation is established, characterized by enhanced expression of phosphorylated Akt, p-ERK, and p65 nuclear translocation and elevated levels of the inflammatory mediators TNFα, IL-6, IL-12B, CCL2, CCL5, and CXCL1 (Recio et al., 2018b; Wojciechowicz and Ma’ayan, 2020). In our study, it is found that GPR84^Y370H variant enhanced inflammation response in ovarian granulosa cells under proinflammatory stimulation. It has been reported that chronic low-grade inflammation is a major contributor to the pathogenesis of premature ovarian insufficiency (Agacayak et al., 2016; Huang et al., 2019). Therefore, it can be speculated that the variants of GPR84 in ovarian GCs would affect ovarian reserve through its role in promoting inflammation, and thus become a potential etiological molecular mechanism of DOR.

There are still some limitations of the study. Firstly, many cases of POI are familial and hereditary, but the heredity of DOR, which is more common in clinics, has not been clearly reported. In this study, only sporadic DOR patients’ samples were collected for sequencing, while no parental DNA samples were obtained, so it is difficult to determine whether the variant is hereditary or not. Secondly, the definite impacts and mechanisms of GPR84^Y370H variant on GCs, oocytes and follicles need thorough study in vitro and in vivo to determine whether they can really impair ovarian reserve. Thirdly, whole exome sequencing only focuses on the protein coding region and the junction between exons and introns, whereas the role of non-coding regulatory regions needs other researches to study.

In conclusion, we conducted WES on 20 young females with DOR and obtained a set of mutated genes related to DOR, where the GPR84^Y370H variant was identified. The deleterious variants of GPR84 could be the potential molecular mechanism of non-age-related pathological DOR through its role in promoting inflammation. The findings of this study can be used as a preliminary research basis for the development of early molecular diagnosis and treatment target selection of DOR.

Data availability statement

The data presented in the study are deposited into the Sequence Read Archive database, accession number PRJNA943550.

Ethics statement

The studies involving human participants were reviewed and approved by The reproductive ethics committee of Tianjin Central Hospital of Gynecology Obstetrics, Tianjin, China. The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

Author contributions

The study was designed and directed by YP and YZ. The clinical information and samples were collected by NL and ZY. The data were analyzed by YP and NL. The experiments were performed by NL. The original manuscript was drafted by NL with contributions of review by YP and final version by all authors. All authors contributed to interpretation of findings.

Funding

This work was supported by the grants from the Bilateral Inter-Governmental S&T Cooperation Project from Ministry of Science and Technology of China [grant number 2018YFE0114300] and the National Natural Science Foundation of China [grant numbers 81972882, 32070752].

Acknowledgments

We thank the patients and the normal controls for their support during this research. We appreciate Prof. Rong Xiang, Prof. Yi Shi, Longlong Wang, Yan Fan (Nankai University) for giving valuable suggestions for the experiment and manuscript. We are grateful to all the medical staff in the Reproductive Medicine Center of Tianjin Central Hospital of Obstetrics and Gynecology for their help in obtaining the samples. Thanks to the Robetta server for help in protein structure prediction.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2023.1154067/full#supplementary-material

References

Adzhubei, I., Jordan, D. M., and Sunyaev, S. R. (2013). Predicting functional effect of human missense mutations using PolyPhen-2. Curr. Protoc. Hum. Genet. Chapter 7, Unit7.20. doi:10.1002/0471142905.hg0720s76

PubMed Abstract | CrossRef Full Text | Google Scholar

Agacayak, E., Yaman Goruk, N., Kusen, H., Yaman Tunc, S., Basaranoglu, S., Icen, M. S., et al. (2016). Role of inflammation and oxidative stress in the etiology of primary ovarian insufficiency. Turk J. Obstet. Gynecol. 13, 109–115. doi:10.4274/tjod.00334

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashburner, M., Ball, C. A., Blake, J. A., Botstein, D., Butler, H., Cherry, J. M., et al. (2000). Gene ontology: Tool for the unification of biology. The gene Ontology consortium. Nat. Genet. 25, 25–29. doi:10.1038/75556

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashrafi, M., Arabipoor, A., Hemat, M., and Salman-Yazdi, R. (2019). The impact of the localisation of endometriosis lesions on ovarian reserve and assisted reproduction techniques outcomes. J. Obstet. Gynaecol. 39, 91–97. doi:10.1080/01443615.2018.1465898

PubMed Abstract | CrossRef Full Text | Google Scholar

Baek, M., Dimaio, F., Anishchenko, I., Dauparas, J., Ovchinnikov, S., Lee, G. R., et al. (2021). Accurate prediction of protein structures and interactions using a three-track neural network. Science 373, 871–876. doi:10.1126/science.abj8754

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, A., Chen, Y., Shen, W., Gao, R., Zhou, W., Yang, S., et al. (2015). Ifit1 protects against lipopolysaccharide and D-galactosamine-induced fatal hepatitis by inhibiting activation of the JNK pathway. J. Infect. Dis. 212, 1509–1520. doi:10.1093/infdis/jiv221

PubMed Abstract | CrossRef Full Text | Google Scholar

Chon, S. J., Umair, Z., and Yoon, M. S. (2021). Premature ovarian insufficiency: Past, present, and future. Front. Cell. Dev. Biol. 9, 672890. doi:10.3389/fcell.2021.672890

PubMed Abstract | CrossRef Full Text | Google Scholar

Cingolani, P., Platts, A., Wang Le, L., Coon, M., Nguyen, T., Wang, L., et al. (2012). A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly. (Austin) 6, 80–92. doi:10.4161/fly.19695

PubMed Abstract | CrossRef Full Text | Google Scholar

Collaborators, G. B. D. T. (2021). Spatial, temporal, and demographic patterns in prevalence of smoking tobacco use and attributable disease burden in 204 countries and territories, 1990-2019: A systematic analysis from the global burden of disease study 2019. Lancet 397, 2337–2360. doi:10.1016/S0140-6736(21)01169-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Depristo, M. A., Banks, E., Poplin, R., Garimella, K. V., Maguire, J. R., Hartl, C., et al. (2011). A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat. Genet. 43, 491–498. doi:10.1038/ng.806

PubMed Abstract | CrossRef Full Text | Google Scholar

Fenton, A. J. (2015). Premature ovarian insufficiency: Pathogenesis and management. J. Midlife Health 6, 147–153. doi:10.4103/0976-7800.172292

PubMed Abstract | CrossRef Full Text | Google Scholar

Gene Ontology, C. (2021). The gene ontology resource: Enriching a GOld mine. Nucleic Acids Res. 49, D325–D334. doi:10.1093/nar/gkaa1113

PubMed Abstract | CrossRef Full Text | Google Scholar

Gillespie, M., Jassal, B., Stephan, R., Milacic, M., Rothfels, K., Senff-Ribeiro, A., et al. (2022). The reactome pathway knowledgebase 2022. Nucleic Acids Res. 50, D687–D692. doi:10.1093/nar/gkab1028

PubMed Abstract | CrossRef Full Text | Google Scholar

Greene, A. D., Patounakis, G., and Segars, J. H. (2014). Genetic associations with diminished ovarian reserve: A systematic review of the literature. J. Assist. Reprod. Genet. 31, 935–946. doi:10.1007/s10815-014-0257-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, S., Xu, B., and Jin, L. (2020). Perinatal outcome in young patients with diminished ovarian reserve undergoing assisted reproductive technology. Fertil. Steril. 114, 118–124 e1. doi:10.1016/j.fertnstert.2020.02.112

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, Y., Hu, C., Ye, H., Luo, R., Fu, X., Li, X., et al. (2019). Inflamm-aging: A new mechanism affecting premature ovarian insufficiency. J. Immunol. Res. 2019, 8069898. doi:10.1155/2019/8069898

PubMed Abstract | CrossRef Full Text | Google Scholar

Jaillard, S., Bell, K., Akloul, L., Walton, K., Mcelreavy, K., Stocker, W. A., et al. (2020). New insights into the genetic basis of premature ovarian insufficiency: Novel causative variants and candidate genes revealed by genomic sequencing. Maturitas 141, 9–19. doi:10.1016/j.maturitas.2020.06.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Jaswa, E. G., Mcculloch, C. E., Simbulan, R., Cedars, M. I., and Rosen, M. P. (2021). Diminished ovarian reserve is associated with reduced euploid rates via preimplantation genetic testing for aneuploidy independently from age: Evidence for concomitant reduction in oocyte quality with quantity. Fertil. Steril. 115, 966–973. doi:10.1016/j.fertnstert.2020.10.051

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiao, X., Ke, H., Qin, Y., and Chen, Z. J. (2018). Molecular genetics of premature ovarian insufficiency. Trends Endocrinol. Metab. 29, 795–807. doi:10.1016/j.tem.2018.07.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanehisa, M., Sato, Y., Furumichi, M., Morishima, K., and Tanabe, M. (2019). New approach for understanding genome variations in KEGG. Nucleic Acids Res. 47, D590–D595. doi:10.1093/nar/gky962

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H., and Durbin, R. (2010). Fast and accurate long-read alignment with Burrows-Wheeler transform. Bioinformatics 26, 589–595. doi:10.1093/bioinformatics/btp698

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, H., Wei, X., Sha, Y., Liu, W., Gao, H., Lin, J., et al. (2020). Whole-exome sequencing in patients with premature ovarian insufficiency: Early detection and early intervention. J. Ovarian Res. 13, 114. doi:10.1186/s13048-020-00716-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, H., Jiang, C., La, B., Cao, M., Ning, S., Zhou, J., et al. (2021). Human amnion-derived mesenchymal stem cells improved the reproductive function of age-related diminished ovarian reserve in mice through Ampk/FoxO3a signaling pathway. Stem Cell. Res. Ther. 12, 317. doi:10.1186/s13287-021-02382-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Livshyts, G., Podlesnaja, S., Kravchenko, S., Sudoma, I., and Livshits, L. (2009). A distribution of two SNPs in exon 10 of the FSHR gene among the women with a diminished ovarian reserve in Ukraine. J. Assist. Reprod. Genet. 26, 29–34. doi:10.1007/s10815-008-9279-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Livshyts, G., Podlesnaja, S., Kravchenko, S., and Livshits, L. (2013). Association of PvuII polymorphism in ESR1 gene with impaired ovarian reserve in patients from Ukraine. Reprod. Biol. 13, 96–99. doi:10.1016/j.repbio.2013.01.178

PubMed Abstract | CrossRef Full Text | Google Scholar

Luscombe, V. B., Lucy, D., Bataille, C. J. R., Russell, A. J., and Greaves, D. R. (2020). 20 Years an orphan: Is GPR84 a plausible medium-chain fatty acid-sensing receptor? DNA Cell. Biol. 39, 1926–1937. doi:10.1089/dna.2020.5846

PubMed Abstract | CrossRef Full Text | Google Scholar

Mckenna, A., Hanna, M., Banks, E., Sivachenko, A., Cibulskis, K., Kernytsky, A., et al. (2010). The genome analysis Toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20, 1297–1303. doi:10.1101/gr.107524.110

PubMed Abstract | CrossRef Full Text | Google Scholar

Muzii, L., Di Tucci, C., Di Feliciantonio, M., Galati, G., Di Donato, V., Musella, A., et al. (2018). Antimullerian hormone is reduced in the presence of ovarian endometriomas: A systematic review and meta-analysis. Fertil. Steril. 110, 932–940 e1. doi:10.1016/j.fertnstert.2018.06.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Ng, P. C., and Henikoff, S. (2003). Sift: Predicting amino acid changes that affect protein function. Nucleic Acids Res. 31, 3812–3814. doi:10.1093/nar/gkg509

PubMed Abstract | CrossRef Full Text | Google Scholar

Nilsson, E., Larsen, G., Manikkam, M., Guerrero-Bosagna, C., Savenkova, M. I., and Skinner, M. K. (2012). Environmentally induced epigenetic transgenerational inheritance of ovarian disease. PLoS One 7, e36129. doi:10.1371/journal.pone.0036129

PubMed Abstract | CrossRef Full Text | Google Scholar

PRACTICE COMMITTEE OF THE AMERICAN SOCIETY FOR REPRODUCTIVE MEDICINE. ELECTRONIC ADDRESS AND PRACTICE COMMITTEE OF THE AMERICAN SOCIETY FOR REPRODUCTIVE MEDICINE (2020). Testing and interpreting measures of ovarian reserve: A committee opinion. Fertil. Steril. 114, 1151–1157. doi:10.1016/j.fertnstert.2020.09.134

PubMed Abstract | CrossRef Full Text | Google Scholar

Qin, Y., Jiao, X., Simpson, J. L., and Chen, Z. J. (2015). Genetics of primary ovarian insufficiency: New developments and opportunities. Hum. Reprod. Update 21, 787–808. doi:10.1093/humupd/dmv036

PubMed Abstract | CrossRef Full Text | Google Scholar

Recio, C., Lucy, D., Iveson, P., Iqbal, A. J., Valaris, S., Wynne, G., et al. (2018a). The role of metabolite-sensing G protein-coupled receptors in inflammation and metabolic disease. Antioxid. Redox Signal 29, 237–256. doi:10.1089/ars.2017.7168

PubMed Abstract | CrossRef Full Text | Google Scholar

Recio, C., Lucy, D., Purvis, G. S. D., Iveson, P., Zeboudj, L., Iqbal, A. J., et al. (2018b). Activation of the immune-metabolic receptor GPR84 enhances inflammation and phagocytosis in macrophages. Front. Immunol. 9, 1419. doi:10.3389/fimmu.2018.01419

PubMed Abstract | CrossRef Full Text | Google Scholar

Rouen, A., Rogers, E., Kerlan, V., Delemer, B., Catteau-Jonard, S., Reznik, Y., et al. (2022). Whole exome sequencing in a cohort of familial premature ovarian insufficiency cases reveals a broad array of pathogenic or likely pathogenic variants in 50% of families. Fertil. Steril. 117, 843–853. doi:10.1016/j.fertnstert.2021.12.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwarz, J. M., Rodelsperger, C., Schuelke, M., and Seelow, D. (2010). MutationTaster evaluates disease-causing potential of sequence alterations. Nat. Methods 7, 575–576. doi:10.1038/nmeth0810-575

PubMed Abstract | CrossRef Full Text | Google Scholar

Skiadas, C. C., Duan, S., Correll, M., Rubio, R., Karaca, N., Ginsburg, E. S., et al. (2012). Ovarian reserve status in young women is associated with altered gene expression in membrana granulosa cells. Mol. Hum. Reprod 18, 362–371. doi:10.1093/molehr/gas008

PubMed Abstract | CrossRef Full Text | Google Scholar

Strauss, J. F., and Williams, C. J. (2019). “Chapter 8 - ovarian Life cycle,” in Yen and jaffe's reproductive endocrinology. Editors J. F. STRAUSS, and R. L. BARBIERI Eighth Edition (Philadelphia: Elsevier).

Google Scholar

Tang, R., and Yu, Q. (2020). The significance of FMR1 CGG repeats in Chinese women with premature ovarian insufficiency and diminished ovarian reserve. Reprod. Biol. Endocrinol. 18, 82. doi:10.1186/s12958-020-00645-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, T. T., Ke, Z. H., Song, Y., Chen, L. T., Chen, X. J., Feng, C., et al. (2013). Identification of a mutation in GDF9 as a novel cause of diminished ovarian reserve in young women. Hum. Reprod. 28, 2473–2481. doi:10.1093/humrep/det291

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Q., Li, D., Cai, B., Chen, Q., Li, C., Wu, Y., et al. (2019). Whole-exome sequencing reveals SALL4 variants in premature ovarian insufficiency: An update on genotype-phenotype correlations. Hum. Genet. 138, 83–92. doi:10.1007/s00439-018-1962-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Warde-Farley, D., Donaldson, S. L., Comes, O., Zuberi, K., Badrawi, R., Chao, P., et al. (2010). The GeneMANIA prediction server: Biological network integration for gene prioritization and predicting gene function. Nucleic Acids Res. 38, W214–W220. doi:10.1093/nar/gkq537

PubMed Abstract | CrossRef Full Text | Google Scholar

Wojciechowicz, M. L., and Ma'Ayan, A. (2020). GPR84: An immune response dial? Nat. Rev. Drug Discov. 19, 374. doi:10.1038/d41573-020-00029-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Zaucha, J., Heinzinger, M., Tarnovskaya, S., Rost, B., and Frishman, D. (2020). Family-specific analysis of variant pathogenicity prediction tools. Nar. Genom Bioinform 2, lqaa014. doi:10.1093/nargab/lqaa014

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, M., Feng, F., Chu, C., Yue, W., and Li, L. (2019). A novel EIF4ENIF1 mutation associated with a diminished ovarian reserve and premature ovarian insufficiency identified by whole-exome sequencing. J. Ovarian Res. 12, 119. doi:10.1186/s13048-019-0595-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, F., Lan, Y., Chen, T., Xin, Z., Liang, Y., Li, Y., et al. (2020). Live birth rate comparison of three controlled ovarian stimulation protocols for in vitro fertilization-embryo transfer in patients with diminished ovarian reserve after endometrioma cystectomy: A retrospective study. J. Ovarian Res. 13, 23. doi:10.1186/s13048-020-00622-x

PubMed Abstract | CrossRef Full Text | Google Scholar