All experimental procedures were approved by the Animal Ethical and Welfare Committee of Zhejiang Chinese Medical University (IACUC-20210503-05). A total of 6 female Sprague Dawley (SD) rats (3-week-old, and 50 $\pm$ 5 g in weight) had free access to adequate food and water and were adaptively fed for 1 week. All rats were randomly separated into POF and control groups. To establish POF models, rats were administered intraperitoneally with 80 mg/kg 4-vinylcyclohexene diepoxide (VCD) solution, which had been dissolved in 2.5 mL/kg/day sesame oil, for 15 consecutive days. Control rats did not receive any treatment. The ovaries were collected under anesthesia 45 days after the initial treatment, and stored at –80 °C.

RNA-seq analysis of the ovaries was conducted by Aksomics Inc. (Shanghai, China). Briefly, total RNA was extracted using TRIzol reagent (Invitrogen, catalog number15596026, Thermo Fisher, New York, NY, USA), and twice treated with DNase I (Ambion, catalog number EN0521, Thermo Fisher, New York, NY, USA) for 30 min at 37 °C. Next, 3 µg extracted RNA was subjected to RiboMinus Eukaryote Kit (catalog number A15020, life technologies, New York, NY, USA) to remove any ribosomal RNA, and subsequently treated with RNase R (catalog number RNR07250, epicentre, Epicentre Biotechnologies, Madison, WI, USA). Purified RNA was treated with RNase R and purified with Trizol. Next, a RNA sequencing library was constructed, and raw sequencing data were achieved on the Illumina NovaSeq 6000. Quantile normalization and data preprocessing were performed using the limma computational approach [18]. Our circRNA RNA-seq data were uploaded into the Gene Expression Omnibus repository with the accession number GSE221728.

Transcriptome microarrays from POF granulosa cells (n = 3) and controls (n = 3) were acquired from the GSE33423 dataset in the Gene Expression Omnibus repository [19]. This dataset was based on the Affymetrix platform.

Expression levels of circRNAs and mRNAs were compared between POF and control groups using limma [18]. For differentially expressed circRNAs, the criteria were set as $|$fold change (FC)$|$ $>$1.5 and p-value 0.05. Furthermore, mRNAs with differential expressions were screened based upon $|$FC$|$ $>$1.2 and p-value 0.05.

The CircAtlas 2.0 database, which is an integrated resource of 1,007,087 circRNAs from 1070 RNA-seq profiles acquired from 19 normal tissue specimens among 6 vertebrate species [20], was adopted to infer the interactions between circRNAs and miRNAs. Through utilizing three online databases comprising miRanda [21], miRDB [22], and miTarBase [23], miRNAs that targeted down- or upregulated mRNAs were estimated and intersected through any two databases. After integrating the circRNA–miRNA and miRNA–miRNA relationships, Cytoscape software (version3.7.2, The Cytoscape Consortium, San Diego, CA, USA) was used to establish the competitive endogenous RNA (ceRNA) networks [24].

To evaluate the interactions between proteins derived from mRNAs with differential expression, they were imported into the STRING online tool using the criteria of required confidence (combined score) $>$0.2 [25]. Protein–protein interaction (PPI) network and subnetwork were constructed using Cytoscape software [24].

Gene Ontology [26] and Kyoto Encyclopedia of Genes and Genomes (KEGG) [27] analyses were used to analyze the functional annotations of mRNAs with differential expression or those from the ceRNA network. A p-value was computed by Fisher’s exact test. KEGG pathways were integrated and visualized using the Pathview web tool [28].

To gain insights into the biological processes as well as to predict the signaling pathways underlying POF-specific necroptosis genes, gene set enrichment analysis (GSEA) was performed to determine the differences in predefined gene sets between low and high expression groups using GSEA software (GSEA4.3.2, Broad Institute, Massachusetts, USA), which was acquired from the Broad Institute website, to divide the median expression value of each gene [29]. Terms with a false discovery rate (FDR) 0.05 were regarded as being statistically different.

The abundance levels of 28 immune cell types were estimated by conducting single-cell GSEA (ssGSEA) [30]. Marker genes of immune signatures were acquired from previously reported literature [31, 32]. Then, the transcriptional levels of immune checkpoints, chemokines, HLA molecules, and immune receptors were measured.

Student’s t test was performed to compare the two groups. Pearson’s test was executed for correlation analysis. All statistical analyses were implemented utilizing appropriate R packages, and a p-value 0.05 was regarded as being statistically different.

To identify the circRNAs linked to POF, the current study implemented circRNA sequencing analysis of POF (n = 3) and control (n = 3) rats. CircRNAs with differential expression in POF versus controls were screened based upon the criteria of $|$FC$|$ $>$1.5 and p-value 0.05. Consequently, 1266 circRNAs displayed notable upregulation in POF, with 1283 displaying downregulation (Fig. 1A–C). Especially, according to $|$FC$|$, we listed the top twenty up- and downregulated circRNAs in POF versus the controls (Fig. 1D; Table 1).

Fig. 1.

Characterization of circular RNAs (circRNAs) with differential expression in premature ovarian failure (POF) rats (n = 3) versus controls (n = 3). (A,B) Scatter and volcano plots illustrate the circRNAs with differential expression in POF compared to controls based upon the criteria of $|$FC$|$ $>$1.5 together with p-value 0.05. Blue denotes downregulated circRNA, while red indicates upregulated circRNA. (C) Heatmap illustrates the expression levels of circRNAs with differential expression across POF and control specimens. Colors from blue to red represent downregulated to upregulated circRNA levels. (D) Heatmap depicts the top twenty up- or downregulated circRNAs in POF compared to controls.

Using the GSE33423 dataset, we acquired transcriptome microarrays of granulosa cells from POF rats (n = 3) and controls (n = 3). In accordance with the criteria of $|$FC$|$ $>$1.2 and p-value 0.05, 1101 upregulated and 1168 downregulated mRNAs were determined in POF versus the controls (Fig. 2A–C). Fig. 2D and Table 2 exhibit the top twenty up- and downregulated mRNAs in POF.

Fig. 2.

Characterization of mRNAs with differential expression in POF (n = 3) versus controls (n = 3) in the GSE33423 dataset. (A,B) Scatter and volcano plots displaying the mRNAs with differential expression in POF compared to control specimens after applying the criteria of $|$FC$|$ $>$1.2 together with p-value 0.05. Blue represents downregulated mRNA and red indicates upregulated mRNA. (C) Heatmap exhibits the expression levels of mRNAs with differential expression across POF and controls. Colors from blue to red denote downregulated to upregulated mRNA levels. (D) Heatmap illustrates the top twenty up- or downregulated mRNAs in POF versus control specimens.

Next, the current study further probed the biological significance of mRNAs with abnormal expressions. As illustrated in Fig. 3A, diverse metabolic (carboxylic acid, lipid, organic acid, cholesterol, sterol, etc.) and biosynthetic (steroid, isoprenoid, lipid, cholesterol, etc.) processes were notably enriched. In addition, mRNAs with differential expression were especially linked to intracellular membrane-bounded organelles (Fig. 3B) and transaminase activity (Fig. 3C). For KEGG pathway enrichment analysis, metabolic (cysteine and methionine, 2-Oxocarboxylic acid, phenylalanine, etc.), biosynthesis (steroid, terpenoid backbone, valine leucine and isoleucine, phenylalanine tyrosine, tryptophan, etc.), and necroptosis pathways were prominently correlated to mRNAs with differential expression (Fig. 3D). The accumulated evidence unveiled that necroptosis participates in gynecological endocrine diseases [33]. Nonetheless, no current studies have uncovered the function of necroptosis in POF. Notably, the necroptosis pathway was visualized, and among mRNAs with differential expression, IFNA1, TRADD, H2AFY2, STAT5A, SPATA2, PYGL, PPID, PLA2G4A, and HMGB1 were enriched in this pathway (Fig. 3E), thereby demonstrating the potential implication of necroptosis in POF.

Fig. 3.

Exploration of the biological significance of differentially expressed mRNAs in POF. (A) Biological process, (B) cellular component, (C) molecular function, and (D) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway terms for mRNAs with differential expression. (E) Visualization of necroptosis pathway enriched by mRNAs with differential expression.

Next, the circRNA-targeted miRNAs with differential expressions were inferred from the CircAtlas 2.0 database. Totally, there were 82 circRNA–miRNA relationships in the POF dataset (Supplementary Table 1). By integrating any two databases (miRDB, miRanda, and miTarBase), a total of 271 miRNAs were estimated that potentially targeted the downregulated mRNAs (Fig. 4A) and 3 miRNAs that potentially targeted the upregulated mRNAs (Fig. 4B). Thereafter, circRNA–miRNA–mRNA regulatory networks were comprised (Supplementary Table 2). Fig. 4C depicts the circRNA–miRNA–mRNA network of downregulated necroptosis genes, comprising rno_circRNA_004995-rno-miR-148b-5p-H2afy2, rno_circRNA_016998-rno-miR-29a-5p-Hmgb1, and rno_circRNA_017593-rno-miR-29a-5p-Hmgb1. In addition, Fig. 4D illustrates the circRNA–miRNA–mRNA network of upregulated necroptosis genes, which is composed of rno_circRNA_015900-rno-miR-935-Stat1, rno_circRNA_007946-rno-miR-328a-3p-Stat5a, rno_circRNA_007947-rno-miR-328a-3p-Stat5a, rno_circRNA_005064-rno-miR-18a-5p-Stat1, rno_circRNA_005064-rno-miR-18a-5p-Stat5a, rno_circRNA_005115-rno-miR-22-3p-Stat1, rno_circRNA_009028-rno-miR-342-5p-Stat1, rno_circRNA_011240-rno-miR-1224-Stat5a, and rno_circRNA_016078-rno-miR-711-Stat5a.

Fig. 4.

Establishment of necroptosis-based circRNA–miRNA–mRNA networks in POF. (A) Venn diagram depicting the microRNAs (miRNAs) of the downregulated mRNAs by integrating any two databases from miRDB, miRanda, and miTarBase. (B) Venn diagram exhibiting the miRNAs of the upregulated mRNAs by integrating any two databases from miRDB, miRanda, and miTarBase. (C) The circRNA–miRNA–mRNA network of downregulated necroptosis genes. Hexagon: circRNA; diamond: miRNA; circle: mRNA. (D) The circRNA–miRNA–mRNA network of upregulated necroptosis genes.

Then, we investigated the biological implications of mRNAs derived from the circRNA–miRNA–mRNA network. Consequently, they exhibited prominent associations with biosynthetic (sterol, cholesterol, isoprenoid, steroid, lipid, etc.) and metabolic (cholesterol, sterol, alpha-amino acid, etc.) processes (Fig. 5A). Moreover, carboxylic acid binding, protein phosphatase binding, etc. (Fig. 5B) and intracellular membrane-bounded organelles (Fig. 5C) were notably enriched. Various biosynthetic (steroid, terpenoid backbone, phenylalanine tyrosine, tryptophan, arginine, etc.) and metabolic (2-Oxocarboxylic acid, phenylalanine, cysteine, methionine, tyrosine, etc.) pathways were notably linked to mRNAs that formed the circRNA–miRNA–mRNA network (Fig. 5D). Specifically, we visualized the tyrosine metabolism pathway that was enriched by Aldh3b1 and Got1 (Fig. 5E).

Fig. 5.

Biological implications of mRNAs derived from the circRNA–miRNA–mRNA network. (A) Biological process, (B) molecular function, (C) cellular component, and (D) KEGG pathway terms for mRNAs derived from the circRNA–miRNA–mRNA network. (E) Visualization of the tyrosine metabolism pathway enriched by mRNAs derived from the circRNA–miRNA–mRNA network.

Supplementary Fig. 1 illustrates the complex interactions between the differentially expressed mRNAs. Then, a PPI subnetwork of POF-specific-necroptosis genes was conducted, which comprised STAT1, STAT5A, PLA2G4A, HMG1L1, HMGB1, AGER, EGFR, HDAC7, IFNA1, IL10RB, IL27RA, PYGL, SOCS1, TRADD, CXCL10, DDX5, EZH2, FADS2, FER, H2AFY2, HIST1H2AF, IFI44L, IL27, IRGM, MX1, NFKB2, PAFAH2, PEMT, PGM2L1, PGR, PHKA2, and PLB1 (Fig. 6A). Their differential expression was demonstrated in POF versus the controls (Fig. 6B). Afterward, the biological processes and KEGG pathways underlying POF-specific necroptosis genes were probed. Here, H2AFY2, HMGB1, STAT1, and STAT5A exhibited positive associations with the negative regulation of mitotic nuclear division and female meiotic nuclear division, while it also regulated the centrosome duplication processes and Fanconi anemia, homologous recombination, and the DNA replication pathways. Furthermore, negative correlations were found relating to the receptor signaling pathway, via JAK–STAT, in skeletal muscle fiber development, and cellular-modified amino acid catabolic processes alongside steroid biosynthesis, glycosaminoglycan degradation, and 2-Oxocarboxylic acid metabolism (Fig. 6C–J).

Fig. 6.

Identification of POF-specific necroptosis genes and any underlying molecular mechanisms. (A) Protein–protein interaction (PPI) subnetwork depicts POF-specific necroptosis genes. (B) Heatmap exhibits the transcriptional levels of POF-specific necroptosis genes across POF and control specimens. Colors from blue to red indicate downregulated to upregulated levels of POF-specific necroptosis genes. (C,D) Gene set enrichment analysis (GSEA unveils the biological processes and KEGG pathways significantly correlated to H2AFY2. (E,F) GSEA unveils the biological processes and KEGG pathways that exhibit significant correlations to HMGB1. (G,H) GSEA predicts the biological processes and KEGG pathways significantly associated with STAT1. (I,J) GSEA reveals the biological processes and KEGG pathways significantly linked to STAT5A.

Utilizing the ssGSEA algorithm, we estimated the abundance of 28 immune cell types across POF and control specimens (Fig. 7A). Significant relationships between immune cells were observed (Fig. 7B). A higher abundance of most immune cells was exhibited in POF samples compared to the control, such as activated B cell, activated CD8${}^{+}$ T cell, activated dendritic cell, CD56 bright and CD56 dim natural killer cells, central memory CD4${}^{+}$ and CD8${}^{+}$ T cell, effector memory CD8${}^{+}$ T cell, macrophage, mast cell, T follicular helper cell, and Type 1 T helper cells (Fig. 7C). In addition, the current study observed differences between the POF and control samples in the transcriptional levels of immune checkpoints, chemokines, HLA molecules, and immune receptors (Fig. 7D–G).

Fig. 7.

Immunological features in POF. (A) Heatmap illustrates the abundance of diverse immune cell types across POF and control specimens. Colors from blue to red denote low to high infiltration of immune cells. (B) Heatmap exhibits the relationships between diverse immune cell types across POF and controls. Blue denotes the negative correlation coefficient, while red is the positive correlation coefficient. The shade of the colors is proportional to the correlation. p-value is marked in boxes. (C) Box plot illustrates the abundance score of diverse immune cell types in POF versus control specimens. (D) Heatmap visualizes the transcriptional levels of immune checkpoints across POF and controls. Colors from blue to red denote downregulated to upregulated transcriptional levels. (E–G) Heatmaps depict the transcriptional levels of (E) chemokines, (F) Human leukocyte antigen (HLA) molecules, and (G) immune receptors across POF and control specimens.

Further analysis exhibited that POF-specific necroptosis genes were significantly linked to most immune cells (Fig. 8A), immune checkpoints (Fig. 8B), chemokines (Fig. 8C), HLA molecules (Fig. 8D), and immune receptors (Fig. 8E).

Fig. 8.

Associations between POF-specific necroptosis genes and POF immunological features. (A) Heatmap illustrates the relationships between POF-specific necroptosis genes and diverse immune cell types in POF and control specimens. Blue denotes the negative correlation coefficient, while red is the positive correlation coefficient. The shade of colors is proportional to the correlation. (B–E) Heatmaps exhibit the associations between POF-specific necroptosis genes and (B) immune checkpoints, (C) chemokines, (D) HLA molecules, and (E) immune receptors in POF and control samples.