CXCL10-based gene cluster model serves as a potential diagnostic biomarker for premature ovarian failure

Data collection

Transcriptome data associated with POF was obtained from the Gene Expression Omnibus (GEO) database (dataset: GSE39501). This dataset comprised gene expression data from three POF mouse samples (cases, including GSM970254, GSM970255, and GSM970256) and three wild-type mouse samples (controls, including GSM970248, GSM970249, and GSM970250). The POF mouse samples were ovaries obtained from C57BL/6 mice with induced POF, while the control mouse samples were ovaries obtained from wild-type C57BL/6 mice. The three POF mouse samples were further divided into high and low CXCL10 expression groups using as a threshold the median expression level of CXCL10 (4.059386919). The data underwent standardization using RunNorm, a tool developed based on the Limma tool (Ritchie et al., 2015).

Principal component analysis (PCA)

PCA is a standard technique for reducing dimensionality in multivariate statistical analysis, emphasizing the most variance-contributing dataset features. By using PCA, we reduced dataset dimensionality and created a two-dimensional scatter map for an initial sample distribution overview (Li et al., 2022). PCA was performed using the PCA expression plot feature in RStudio. The transcriptional levels of CXCL10 in POF and control groups were visualized using ViolinPlot, which generates violin plots based on the result file and grouping information output with the RunNorm application.

Receiver operating characteristic (ROC) curve

The ROC curve is a common clinical practice tool that can be used to evaluate the performance of classifiers. An area under the curve (AUC) value exceeding 0.5 and approaching 1, indicates better diagnostic accuracy. AUC values ranging from 0.5 to 0.7 indicate low accuracy, values ranging from 0.7 to 0.9 indicate moderate accuracy, and values exceeding 0.9 indicate high accuracy. ROC analysis was performed to assess the optimum cutoff value, sensitivity, specificity, and AUC with a 95% confidence interval (CI). The ROC curve of CXCL10 was plotted and analyzed using the PlotROC application, based on the R language package pROC (Sachs, 2017).

The Least Absolute Shrinkage and Selection Operator (LASSO) regression is a commonly used data mining method in machine learning known for variable selection and complexity adjustment within generalized linear models to mitigate multicollinearity in regression analysis. We utilized Baiyin Cloud to implement RunLASSO for selecting feature genes impacted by CXCL10 in POF with diagnostic potential. This application, based on the R language LASSO package, establishes a LASSO model (Friedman, Hastie & Tibshirani, 2010) while generating Lambda plots, LASSO model diagrams, and ROC curves. Transcriptional expression of model genes was depicted using PlotBox showcasing them in a box plot.

Gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Set Enrichment Analysis (GSEA) pathway analysis

The “clusterProfiler” package was used to perform GO, KEGG, and GSEA pathway analyses of the differentially expressed genes (DEGs) to annotate biological functions and enriched pathways. Enrichment results with false discovery rate (FDR) <0.05 were considered significant. POF dysregulation genes influenced by CXCL10 underwent enrichment analysis using RunMutiGroupclusterProfiler, an application developed based on the clusterProfiler package in R language (Yu et al., 2012), for GO and KEGG. The significance of the enrichment results was based on a p-value < 0.05. The enrichment results for biological processes (BP) or KEGG pathways were visualized using the PlotClusterBubble application on Baiyin Cloud.

Cell proliferation assay

The cell counting kit-8 (CCK-8) method was used for assessing cell proliferation. KGN cells were cultured in a 96-well plate for 24 or 48 h. Then, 10 μL of CCK-8 medium (HY-K0301; MCE, Concord, CA, USA) was added to each well. After a 1-h incubation at 37 °C with 5% CO₂, the optical density (OD) at 450 nm was measured. Each measurement was conducted in triplicate.

Cell culture and transfection

Human ovarian granulosa cells (KGN cells) were obtained from the American Type Culture Collection (ATCC) and cultured in RPMI-1640/F12 medium (Gibco, Waltham, MA, USA) supplemented with 10% fetal bovine serum (FBS; Gibco, Waltham, MA, USA), 100 U/mL penicillin, and 100 mg/mL streptomycin. The cells were grown at 37 °C with 5% CO₂. CXCL10 short interfering RNAs (siRNAs) were synthesized by Shanghai Gene Pharma Co., Ltd. (Shanghai, China). Lipfectmine3000 (Thermo Fisher, Waltham, MA, USA) was used to transfect the siRNAs at a final concentration of 50 nM, following the manufacturer’s instructions.

For cisplatin (DDP) treatment (Wang et al., 2022), the KGN cells were treated with different concentrations of substance (0, 10, 20, 40, and 80 μg/mL) for 24 or 48 h. The utilized siRNA sequences were as follows: Si-CXCL10, 5′-GCCTTATCTTTCTGACTCTAA-3′; Si-NC, 5′-CCTAAGGTTAAGTCGCCCTCG-3′.

RNA extraction and quantitative reverse transcription (RT-qPCR)

Total RNA from KGN cells was extracted using TRIzol^® Reagent (Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer’s instructions. DNA was digested using Dnase I (18068015; Thermo Fisher Scientific, Waltham, MA, USA) at 25 °C for 1 h. Contamination assessment was performed using a 20-minute agarose gel electrophoresis at 130 V. The quality of the RNA was evaluated using NanoDrop One (Thermo Fisher Scientific, Waltham, MA, USA), targeting an A260/280 ratio of 1.8–2.0. RNA integrity was evaluated using an Agilent 2100 bioanalyzer, with an RNA integrity number (RIN) ranging from 7 to 10. The SPUD assay was used to test for inhibition, according to the SPUD Assay for Detection of Assay Inhibitors Protocol (Merck, Rahway, NJ, USA). For each microgram of total RNA, 0.5 µg of primer or primer-adaptor was added to 2 µg of total RNA in a sterile Rnase-free microcentrifuge tube.

The mixture, in a total volume of ≤15 µL of water, was heated to 70 °C for 5 min.

Reverse transcription of cDNA was performed using 200 units of M-MLV reverse transcriptase (M170A; Promega, Madison, WI, USA) in a 25 μL volume, incubating for 60 min at 37 °C with random primers or 42 °C with other primers or primer-adaptors. RT-qPCR was carried out on ABI7500 Real-Time PCR Detection System with a 20 μL reaction system, including 2 μL of cDNA and AceQ Universal SYBR qPCR Master Mix (Q511; Vazyme, Nanjing, China) containing 50 μmol/L dNTP, 1.5 mmol/L Mg²⁺, 20 U AceTaq DNA polymerase according to the protocol. The PCR process comprised an initial pre-denaturing step at 95 °C for 2 min, followed by 42 cycles of 95 °C for 30 s, 58 °C for 30 s, and 72 °C for 30 s.

The dissolution curve must be unimodal. Standard curves with slope and y-intercept were developed, aiming for an r² value of 0.99. The dynamic range of the PCR reaction should exhibit linearity, with a Cq variation at the lower limit of 10. The limit of detection (LOD) was determined to be 2.5 target molecules with a 95% confidence interval. The 2^−ΔΔCt method was used to determine the relative expression of CXCL10 (NM_001565.4), normalized to GAPDH (NM_002046.7), using the Applied Biosystems ViiA™ 7 Real-Time PCR System, with GAPDH serving as control. The primer sequences used for GAPDH and CXCL10 were as follows: 5′-GTGGCATTCAAGGAGTACCTC-3′, 5′-TGATGGCCTTCGATTCTGGATT-3′, 5′-AAGTATGACAACAGCCTCAAG-3′ and 5′-TCCACGATACCAAAGTTGTC-3′.

The amplicon length was approximately 100 bp, corresponding to the gene’s exon, which was verified in NCBI. The comparative 2^−ΔΔCt method assessed the stability of each gene by determining the standard deviation of Cq differences. No Cq value was detected for the No Template Control (NTC). Outliers, defined as values significantly deviating from the other two values within three technical replicates were excluded from statistical analysis. Each experiment involved three biological and technical replicates.

Western blot

Cells were suspended in CellScope’s RIPA lysis buffer, and protein concentration was determined using the bicinchoninic acid assay (Thermo Fisher Scientific, Waltham, MA, USA). A 40 μg protein sample was loaded on a 10% gel for sodium dodecyl polyacrylamide gel electrophoresis (SDS-PAGE). The proteins were transferred to a polyvinylidene fluoride (PVDF) membrane which was blocked using 5% skim milk and subsequently incubated with the following antibodies: anti-GAPDH (ab9484; Abcam, Cambridge, UK), anti-CXCL10 (ab306587; Abcam, Cambridge, UK), anti-peroxisome proliferator-activated receptor (PPAR) β (ab310323; Abcam, Cambridge, UK), anti-long-chain acyl-coenzyme A synthetases 1 (ACSL1, ab177958; Abcam, Cambridge, UK), anti-caspase (ab32351; Abcam, Cambridge, UK), anti-B-cell lymphoma 2 (Bcl-2, ab182858; Abcam, Cambridge, UK), and anti-Bcl-2-associated X protein (Bax, ab32503; Abcam, Cambridge, UK). The blots were then incubated with a secondary goat anti-rabbit IgG (H+L) antibody (A16104; Thermo Fisher, Waltham, MA, USA) for 1 h at room temperature. Immunoreactive bands were visualized using an enhanced chemiluminescence reaction (Pierce, Appleton, WI, USA) following standard protocols.

Annexin V-fluorescein isothiocyanate (FITC)/propidium iodide (PI) assay

KGN cells were placed in a six-well culture dish and cultured in an incubator at 37 °C with 5% CO₂. After 48 h of culture, cells from each group were then treated with 5 μL of annexin V-FITC and 5 μL of PI solutions (E-CK-A211; Elabscience, Houston, TX, USA) at 37 °C for 20 min in the absence of light. Apoptotic cells were counted using a flow cytometer (BD Biosciences Co. Ltd, Franklin Lakes, NJ, USA).

Enzyme-linked immunosorbent assay (ELISA)

Human estradiol (E2) levels were detected using human estradiol ELISA kit (MM-0777H1; Jiangsu Meimian Industrial Co., Ltd., Jiangsu, China) according to the instructions of the manufacturer. The absorbance was measured at 450 nm using a microplate reader.

Statistical analysis

We performed statistical analysis using SPSS 21.0 software (SPSS Inc., Chicago, IL, USA). Data are expressed as mean ± standard deviation (SD). To assess differences between two or more groups, we employed an unpaired Student’s t-test or ${\chi }^2$-test, and a one-way analysis of variance, respectively. Tukey’s multiple comparison test was applied when necessary. Each experiment was conducted with a minimum of three biological replicates, and statistical significance was defined as p-levels < 0.05.

Expression of CXCL10 in POF

The PCA analysis revealed that the expression of CXCL10 was higher in the POF group than in the control group (Fig. 1A). Additionally, the violin plot demonstrated a significant increase in CXCL10 expression levels in the POF group compared to the control group (Fig. 1B). ROC curve analysis using PlotROC revealed that CXCL10 had a relatively accurate diagnostic potential for POF with an AUC value of 1 (Fig. 1C). Overall, these results indicate that the expression of CXCL10 is dysregulated in POF.

Global regulation of CXCL10 in POF

To understand the global regulation of CXCL10 in POF, we conducted a differential expression analysis to identify DEGs in the control group of POF using Debylimma.

This analysis revealed 3,362 up-regulated and 3,969 down-regulated DEGs (Fig. 2A, Table S1). Furthermore, we found 1,304 up-regulated and 1,315 down-regulated DEGs between the high and low CXCL10 expression groups of POF (Fig. 2B, Table S2). The overlap of DEGs between the POF-control group and the high- and low-CXCL10 expression groups was determined using PlotVenn, as shown in Fig. 2C. A total of 786 DEGs, uniformly up-regulated or down-regulated, were identified as POF-related genes dysregulated by CXCL10. These genes were further analyzed using QuadrantPlot (Fig. 2D).

Enrichment analysis

To explore the potential BP and pathways affected by the CXCL10 gene in POF progression, we conducted GO and KEGG analyses on DEGs using the RunMulti Group cluster Profiler. GO enrichment analysis revealed significant involvement of these DEGs in cell chemotaxis, negative regulation of response to DNA damage stimulus, and diadenosine polyphosphate metabolic processes (Fig. 3A, Table S3). KEGG analysis showed significant enrichment of these genes in the PPAR signaling pathway, adiponectin signaling pathway, drug metabolism-cytochrome P450, and metabolism of xenobiotics by cytochrome P450 (Fig. 3B, Table S3). Gene Set Enrichment Analysis (GSEA) indicated strong activation of the PPAR signaling pathway in POF (Fig. 3C, FDR = 0.007 and enrichment score = 0.813). The map generated using RunPathview displayed significant changes in genes related to the PPAR signaling pathway, including LXRa, ACS, CYP4A1, ACO, CPT-1, MCAD, POAR, and ADIPO. Among these, ACS was the most upregulated in POF (Fig. 4, Table S2). These results indicated that CXCL10 activates the PPAR signaling pathway in POF.

Diagnostic efficacy of CXCL10-based clinical model

ROC curves for CXCL10-dysregulated genes in POF were analyzed using PlotROC. RunLASSO performed characteristic gene selection for these genes, identifying three feature genes, namely Cxc110, Raf1, and ITGA2, with non-zero regression coefficients (lambda.min = 0.0006). The Lambda diagram (Fig. 5A) illustrates the diagnostic performance of the CXCL10-based model across different Lambda values, with the best Lambda at 0.001 and the minimum at 0.0006. The LASSO model diagram (Fig. 5B) indicates the model’s confidence in predicting POF diagnosis. Finally, the ROC curve (Fig. 5C) demonstrates that the CXCL10-based model possesses excellent diagnostic efficiency (AUC = 1) for POF in the training set, serving as a potential diagnostic biomarker for POF.

Moreover, the expression of the model genes significantly differed between the POF and control groups. Additional ROC curves (Fig. S1) further support the diagnostic potential of CXCL10-based model genes. The CXCL10 expression was higher in the case group than in the control group, indicating its association with POF. Similarly, lower Itga2 and Raf1 expression in the case group suggests their potential as diagnostic markers for POF (Fig. 5D).

CXCL10 affects the progression of POF by influencing the PPAR signaling pathway

To verify the role of CXCL10 in POF, we constructed a DDP-induced KGN cell model of injury, silencing the expression of CXCL10 in these cells. Firstly, we determined that DDP at concentrations exceeding 10 μg/mL significantly reduced cell proliferation at both 24 and 48 h (Fig. 6A), indicating that this concentration of DDP can induce KGN cell damage.

Si-CXCL10#2 was chosen for further experiments due to its efficacy in silencing CXCL10 in KGN cells, as confirmed by RT-qPCR analysis (Fig. 6B). RT-qPCR and WB analyses revealed increased expression of CXCL10 in the POF cell model group compared to the control group. However, CXCL10 expression was significantly reduced in the POF+Si-CXCL10 group compared to the si-NC group (Figs. 6C and 6D). Estradiol secretion was significantly inhibited in the KGN cell model of injury after 48 h of culture, but this inhibition was alleviated upon interference with CXCL10 expression (Fig. 6E). These findings demonstrated that CXCL10 interference had a protective effect against DDP-induced granulosa cell damage, preserving granulosa cell function.

Subsequently, we examined the effect of CXCL10 on apoptosis. The results indicated a substantial increase in the apoptosis rate of KGN cells in the POF group compared to the control group. However, silencing CXCL10 resulted in significant inhibition of apoptosis when compared with the si-NC group (Fig. 7A). This was further corroborated by WB analysis, which revealed that the expressions of pro-apoptotic genes Bax and caspase increased while the anti-apoptotic gene Bcl-2 decreased in the POF group compared to the control group. Compared to the Si-NC group, the expressions of Bax and caspase were decreased in the Si-CXCL10 group, and the expression of Bcl-2 was significantly increased (Fig. 7B).

In addition, we examined the effect of CXCL10 on the PPAR signaling pathway, revealing a significant increase in the expressions of PPARβ and ACSL1 in the POF group compared to the control group. However, CXCL10 silencing led to decreased levels of PPARβ and ACSL1 compared to the Si-NC group (Fig. 7C), indicating that CXCL10 impacts the PPAR signaling pathway.