Identification of metabolomic changes and potential therapeutic targets during ovarian aging

Study population and sample collection

A total of 45 women were included in this study. The participants had undergone the first ART cycles. All the participants included in the study have provided written informed consent, and their cause of infertility was all tubal obstruction, patients with PCOS, endometriosis, premature ovarian insufficiency, and chromosome abnormality were excluded.

Follicular phase long-acting protocol: patients were injected 3.75 mg GnRH antagonist (GnRH-a Tryptorelin, Ferring, Germany) on the 2nd day of menstruation if the ultrasound did not find cysts and follicles >10 mm. The patients will visit the hospital 28 days after the injection to be examined on ultrasound+serum FSH, LH, E2, and progesterone (P). Administration of a starting dose of 150–300 IU of human recombinant FSH (rFSH, Gonal F, Serono Ltd., Switzerland), with subsequent adjustment of gonadotropin (Gn) use according to follicular growth. The starting dose is according to age, AFC, and basal hormone levels. Follicular development was detected by using vaginal ultrasound. Ovulation is performed using human chorionic gonadotropin (HCG) (2,000 IU) (Lizhu Ltd., China) and recombinant HCG (250 µg) (Serono Ltd., Switzerland) when at least one follicle was 18 mm in size, and FF and GCs were obtained 36–38 hours after HCG administration.

Depending on their age, patients were divided into three groups (Group A: 20–30 years old; Group B: 30–35 years old; Group C: 35–45 years old). Briefly, follicular aspirates were centrifuged at 1600 RMP for 15 minutes, the supernatant was FF and the sediment was the cell mixture. The purified GCs were obtained by density gradient centrifugation of the sediment as previously described [13]. Both FF and GCs were stored at −80°C until assayed.

Metabolite extraction method

GCs: add 100 µl of the whole target isotope internal standard to a 2 ml EP tube containing 1 × 10⁷ cells and vortex for 30 s to mix the cells. Then add 400 µl of acetonitrile: (methanol: ddH₂O mixed solution (2:2:1, V:v:v)) and vortex oscillation for 30 s. Place the centrifuge tube containing the sample, invade liquid nitrogen for rapid freezing for 5 minutes, take out the centrifuge tube, and freeze-thaw at room temperature on the double panel. Place the centrifuge tube in the 2 ml adapter again, and install it until grinding in the instrument. Oscillate at 60 Hz for 2 min. Repeat the last step twice, take out the centrifuge tube, centrifuge at 12000 rpm 4°C for 10 min, take 450 µl of the supernatant into a 2 ml centrifuge tube, and concentrate in vacuo to dryness. Accurately add 100 µl of acetonitrile (0.1% FA (1:9, v:v) solution (without 2-chlorophenylalanine)) (stored at −20°C) to reconstitute the sample, filter through a 0.22 µm membrane, and add the filtrate to detection bottle.

FF: thaw all samples at 4°C. Transfer 100 µL of each sample into 2 mL centrifuge tubes. Add 100 µL of mixed internal standard solution and 400 µL of methanol (−20°C), then vortex for 60 s. Centrifuge at 4°C for 10 min at 12000 rpm, and then transfer 500 µL of the supernatant from each sample into another 2 mL centrifuge tube. Samples were concentrated to dry in vacuum. Dissolve samples with 150 µL of 80% methanol solution, and centrifuge at 4°C for 10 min at 12000 rpm again to obtain the supernatant for LC-MS.

After the above procedure, 20 µL of each sample was used for quality control (QC), while the rest of the samples were then used for LC-MS analysis.

LC-MS/MS analysis

Chromatographic separations were performed on a Thermo Ultimate 3000 system equipped with an ACQUITY UPLC^®HSS T3 (150 × 2.1 mm, 1.8 µm, Waters, Milford, MA, USA) column maintained at 40°C. The autosampler temperature was 8°C. The analytes were eluted at 0.25 ml/min with a gradient of 0.1% formic acid solution (C) and 0.1% acetonitrile formatted solution (D) or 5 mm ammonium formatted solution (A) and acetonitrile solution (B). After equilibration, 2 µl was injected into each sample. The linear-gradient (v/v) of solvent B increased: 0–1 min, 2% B/D; 1–9 min, 2–50% B/D; 9–12 min, 50–98% B/D; 12–13.5 min, 98% B/D; 13.5–14 min, 98–2% B/D; 14–20 min, 2% d positive model (14–17 min, 2% b negative model).

ESI-MSn experiments were performed on a Thermo Q Exactive HF-X mass spectrometer using spray voltages of 3.5 kV and −2.5 kV in positive and negative modes, respectively. The casing and auxiliary gases were arbitrary units of 30 and 10, respectively. The capillary temperature was 325°C. The analyzer was fully scanned in the mass range m/z 81–1,000 with a mass resolution of 60,000. Data dependent acquisition (DDA) MS/MS experiments were performed using HCD scans. The normalized collision energy was 30 eV. A dynamic exclusion method was used to remove unwanted information from the MS/MS spectrum.

Statistical methods and bioinformatic analysis

The baseline and ovarian stimulation characteristics of patients were analyzed. Mann–Whitney tests were used for continuous variables and categorical variables were examined by Chi-squared tests. The results are shown as the mean ± standard deviation (SD) or % (n/N). All data were analyzed by SPSS 23.0 (IBM Corp., NY, USA), and the threshold was set as a two-tailed P < 0.05.

Comprehensive metabolomic data analysis was performed by R-software (Metabo Analyst 3.0). Firstly, sample normalization of the metabolomics data was performed. Principal component analysis (PCA) of metabolomics data was performed on a normalized dataset. Analysis of variance (ANOVA) was used to analyze the metabolites among the three groups, a P-value < 0.05 was considered as the differential metabolites. Next, Partial Least Squares Discriminant Analysis (PLSDA) was performed in Group A and Group C to decipher differences between them, and Variable Importance Projection (VIP) scores were generated in Metabo Analyst. To identify different metabolites between groups, the threshold was set as VIP ≥ 1, |Log2FC (fold change) |> 1 or <−1 and P-value < 0.05. Metabolic pathway enrichment analysis of the identified metabolites was carried out in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database.

Availability of data and materials

All data generated or analyzed during this study are included in this published article (and its Supplementary Tables).

Baseline characteristics and ovarian induction characteristics of patients

As shown in Table 1, no significant differences were found among the three groups in body mass index (BMI) (22.41 ± 0.48 vs. 21.51 ± 0.45 vs. 23.87 ± 0.55), basal serum FSH (5.87 ± 0.54 vs. 6.66 ± 0.42 vs. 6.40 ± 0.34), basal serum LH (5.21 ± 0.70 vs. 3.51 ± 0.76 vs. 4.38 ± 0.50), AMH (4.35 ± 0.98 vs. 3.16 ± 0.60 vs. 2.75 ± 0.48), AFC (16.53 ± 0.68 vs. 12.93 ± 1.37 vs. 11.80 ± 1.54) and basal endometrial thickness (5.05 ± 1.03 vs. 5.71 ± 1.07 vs. 5.94 ± 1.25).

Compared with group B, group A was associated with a lower total amount of gonadotropin (Gn), more MII oocytes, and more 2PN embryos. Compared with group C, patients in group A used fewer starting dosages of Gn (151.67 ± 16.27 vs. 212.50 ± 13.42), fewer total amount of Gn (2116.67 ± 234.60 vs. 3525.00 ± 286.63), more retrieved oocytes (15.53 ± 1.62 vs. 10.33 ± 1.29), more MII oocytes (13.40 ± 1.28 vs. 8.93 ± 1.18), more 2PN embryos (9.80 ± 0.79 vs. 6.33 ± 0.99) and more transferable embryos (5.78 ± 0.74 vs. 4.31 ± 0.61).

Profile of GC and FF metabolomics

After the normalization of the intensity of metabolites, PCA was performed. The PCA scatter plot of metabolomics data from the FF in the three groups showed that there was a cluster between Group A and Group C, but did not differ significantly between Group A and Group B (Figure 1A). Then, we performed a hierarchical clustering analysis of the screened differential metabolites in all the groups, distinguishing them by different colors. Regions with similar colors represent that metabolites in that region have similar functions or participate in a biological process together, while larger color differences indicate that the functions of metabolites in these regions are far apart. Here, 93 metabolites showed gradual variations in different groups, with 61 gradually elevating and 32 gradually declining (Figure 1B).

Figure 1. Profile of GC and FF metabolomics. (A) PCA scatter plot of metabolomics data from FF in three groups. Each point represents one sample, different groups are marked with different colors, and the area marked by the ellipse is the 95% confidence region of the sample points; (B) heatmap of hierarchical cluster analysis of FF metabolites, the former ten metabolites decreasing with aging, while the latter eight metabolites increasing with aging. The magnitude of relative content in the plot is shown by the difference in color, where the columns represent samples and the rows represent metabolites; (C) PCA scatter plot of metabolomics data from GC in three groups; (D) heatmap of hierarchical cluster analysis of FF metabolites, the former 61 metabolites increasing with aging, while the latter 32 metabolites decreasing with aging.

The PCA plot reflected the original state of GC metabolites, it showed that the three groups cannot be completely separated (Figure 1C). To investigate the differences in GC metabolites at different ages, ANOVA analysis was performed on 3 sets of metabolite data to screen for metabolites with significant differences, and the clustering of metabolites was visualized by heatmap, results indicated that there were ten metabolites declined with aging, while eight metabolites increased with aging (Figure 1D).

Metabolic differences in GC between group A and group C

PCA cannot ignore within-group errors, eliminate random errors unrelated to the purpose of the study, and is not conducive to finding differences between groups, therefore, PLSDA was utilized to obtain metabolites that differed significantly between Group A and Group C (Figure 2A). Then, differential metabolites in the two groups were screened, and we got 21 metabolites that differed greatly between Group A and Group C (Figure 2B): 11 were up-regulated and 10 were down-regulated (Supplementary Table 1). Heatmap showed the top 25 differential metabolites (Figure 2C), among them, five metabolites (24%) belong to amino acids, peptides, and analogs, three metabolites (14%) belong to carbohydrates and carbohydrate conjugates, and five metabolites (24%) cannot be classified to any sub-class, the rest of metabolites belonged to amines (4%), benzoic acids and derivatives (4%), carboxylic acids (5%), fatty acids and conjugates (5%), indoles (5%), methoxyphenyl (5%), pyrimidine nucleotide sugars (5%) and pyrimidine ribonucleotides (5%), respectively (Figure 2D).

Figure 2. Metabolic differences in GC between group A and group C. (A) PLSDA score plot; (B) volcano plot of differential metabolites, each point represents a differential metabolite, the horizontal coordinate represents the log2 fold change of the metabolite, and the vertical coordinate represents the P-value (log10 transformed); (C) heatmap of top 25 differential metabolites; (D) sub-classification of differential metabolites.

Metabolic differences in FF between group A and group C

Similarly, we used PLSDA to screen the differences between Group A and Group C (Figure 3A). Finally, a total of 72 metabolites were selected based on the threshold that we set, volcano plot was used for visualization. Heatmap showed the top 25 differential metabolites (Figure 3B, 3C), among them, 37 were down-regulated, and 35 were up-regulated in Group C, and palmitic acid was up-regulated (Supplementary Table 2). Of the 72 differential metabolites, 13 metabolites (18%) belong to fatty acyls and 11 metabolites (15%) belong to carboxylic acids and derivatives. Moreover, keto acids and derivatives and organooxygen compounds each made up 6% of the total, benzene and substituted derivatives, phenols and steroids, and steroid derivatives accounted for 4%, respectively (Figure 3D).

Figure 3. Metabolic differences in FF between group A and group C. (A) PLSDA score plot; (B) volcano plot of differential metabolites, each point represents a differential metabolite, the horizontal coordinate represents the log2 fold change of the metabolite, and the vertical coordinate represents the P-value (log10 transformed); (C) heatmap of top 25 differential metabolites; (D) sub-classification of differential metabolites.

KEGG metabolic pathway of differential metabolites in FF and GC

To further analyze the function of these differential metabolites, we performed a KEGG enrichment analysis. Through metabolic pathway analysis (MetPA), the associated metabolic pathways of the two groups of differential metabolites were analyzed, the hypergeometric test was used for data analysis and the pathway topology is the relative-betweenness centrality. The pathway impact value represents the importance of the pathway by summing the importance measures of the differential metabolites enriched to the pathway and is divided by the sum of the importance measures of all metabolites on the pathway, which is the pathway impact value. In brief, the larger the value, the more important the pathway.

For FF, there were 126 pathways enriched for the differential metabolites, among them, the P-value for 35 metabolic pathways was less than 0.1, and 22 pathways had P-value < 0.05. Interestingly, we found that the PPAR signaling pathway, linoleic acid metabolism, arginine biosynthesis, cysteine and methionine metabolism, cAMP signaling pathway, pantothenate and CoA biosynthesis, alanine, and aspartate and glutamate metabolism were also significantly enriched (Figure 4A).

Figure 4. KEGG metabolic pathway of differential metabolites in FF and GC. (A) KEGG pathway enrichment analysis of FF metabolites; (B) KEGG pathway enrichment analysis of GC metabolites. The vertical coordinate represents the corresponding pathway and the horizontal coordinate represents the impact value of the pathway. The size of the dot represents the number of differential metabolites in the target metabolic pathway, the larger the dot, the greater the number of metabolites enriched into that pathway. P-values are represented by color, the darker the color, the smaller the p-value.

For GC, 43 metabolic pathways were enriched, 18 pathways had a P-value < 0.1, and 8 pathways had a P-value < 0.05 (Figure 4B). Notably, the pentose phosphate pathway, nicotinate and nicotinamide metabolism, GnRH secretion, and estrogen signaling pathway were also involved in the enrichment of differential metabolites. We then took the intersection of the pathways notably enriched by the two sets of samples and found that 6 pathways were co-enriched by the differential metabolites of FF and GC: protein digestion and absorption, ABC transporters, GABAergic synapse, morphine addiction, synaptic vesicle cycle, cysteine, and methionine metabolism. Among them, we noticed that the pathway, GABAergic synapse, had a high impact value in both FF and GC enrichment. After screening the metabolites involved in this pathway, we found that GABA was down-regulated in GC, while its downstream metabolite succinic acid was down-regulated in FF (Figure 5A). Further ROC curve analysis was performed on these two metabolites and results showed that they all had a favorable predictive value (Figure 5B, 5C).

Figure 5. GABAergic synapse pathway model and ROC curve analysis. (A) KEGG pathway for GABAergic synapse, blue dots represent down-regulated metabolites (https://www.genome.jp/kegg/); (B) ROC curve analysis for GABA; (C) ROC curve analysis for succinic acid.