Transcriptome analysis during 4-vinylcyclohexene diepoxide exposure-induced premature ovarian insufficiency in mice

Chemical

4-Vinylcyclohexene diepoxide (V820469-5 ml) and sesame oil (S905724-50 ml) were purchased from Macklin, Inc. (Macklin, Shanghai, China). VCD was mixed with sesame oil and administered to mice via intraperitoneal injection.

Feeding conditions for experimental animals

A total of 30 specific pathogen free (SPF) seven-week-old female C57BL/6 mice were purchased from the Jiangsu Huachuang Sino Pharmaceutical Technology Co., Ltd. (Jiangsu, China). The mice were acclimatized for one week before the formal experiment, kept in separate cages (six cages, with five mice per cage), and had free access to food and water. The environment provided a normal circadian rhythm and a constant temperature and humidity (24 ± 2 °C/40%). Ethics committee permission was obtained from Lanzhou University’s First Hospital for all animal experiments (LDYYLL2023-105).

Animal model establishment and sample collection

Complete randomization was used to divide the mice into three groups: a control group (n = 10), a V15 experimental group (n = 10), and a V30 experimental group (n = 10). The control group was injected with sesame oil, and the two experimental groups were given VCD dissolved in sesame oil at a concentration of 160 mg/kg via intraperitoneal injection. The injection required for 15 consecutive days. The dose of VCD were based on published reports (Cao et al., 2020) and our pre-experiments. Simultaneous vaginal cytology cell smears were made during the 15-day modeling period for observation of the estrous cycle.

After modeling, the mice were fed normally for different time durations. For the V15 group, serum and ovaries were harvested on day 15 after injection. For the CON and V30 groups, serum and ovaries were collected on day 30. All mice were euthanized via transperitoneal injection of pentobarbital in accordance with the American Veterinary Medical Association (AVMA) Guidelines for Euthanasia of Animals before serum collection and organ (ovary) harvesting. After euthanasia, the eyes of the mice were immediately removed to collect blood, and the blood was allowed to sit overnight in a 4 °C refrigerator. The serum was collected by centrifugation at 3,000 rpm for 10 min. The left ovaries of the mice were stored in a −80 °C freezer, and the right ovaries were completely immersed in tissue fixative for 24 h before being used for subsequent experiments.

Hormone measurement

Serum samples from 10 mice in each group (a total of 30 mice) were used to measure hormone levels. The serum anti-mullerian hormone (AMH), follicle stimulating hormone (FSH) and estrogen (E₂) levels of the mice were determined using matching enzyme-linked immunosorbent assay (ELISA) kits purchased from the Jingmei Biotech Company (Jiangsu, China). The standard curve range for AMH was 0–180pg/ml, and the lower limit of quantification was 11.25 pg/ml; the standard curve range for FSH was 16 U/L, and the lower limit of quantification was 1U/L; the standard curve range for E₂ was 0–200 pmol/L, and the lower limit of quantification was 12.5 pmol/L. The test was performed in strict accordance with the kit instructions.

Follicle count

The right side of the mouse ovary was sectioned (n = 10 from each group). Ovaries were immersed in paraformaldehyde and fixed by paraffin embedding. The embedded ovaries were serially sectioned (5 µm thick sections), and the sections were stained with hematoxylin-eosin (HE). Observation of ovarian morphology was carried out under an orthoptic microscope (OLYMPUS, BX51; Olympus, Tokyo, Japan). The number of follicles at each stage was counted in one of every five sections, and only follicles with visible oocyte nuclei, including primordial follicles, primary follicles (PFs), secondary follicles (SFs), antral follicles (AFs) and the corpus luteum (CLs), were recorded. Ovarian tissue sections in this article were photographed using E3ISPM06300KPA (Hangzhou ToupTek Photonics Co., Ltd., Hangzhou, China).

RNA-seq

Ovarian tissue RNA was extracted using a TRIzol kit (Invitrogen; n = 5 from each group). The RNA Nano 6000 Assay Kit on the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA) was used to evaluate RNA integrity. Poly(A)-tailed mRNA was enriched using oligo (dT) magnetic beads. The obtained mRNA was then fragmented using a fragmentation buffer with divalent cations. Fragmented mRNA served as a template for the synthesis of the first cDNA strand using random oligonucleotides as primers in the presence of M-MuLV reverse transcriptase. Following RNaseH-mediated RNA degradation, the second cDNA strand was synthesized using DNA polymerase I and dNTPs. Purified double-stranded cDNA was subjected to end repair, A-tailing, and adapter ligation, followed by size selection using AMPure XP beads targeting cDNA fragments ranging from 370 to 420 bp. PCR amplification was performed on the selected fragments, followed by purification using AMPure XP beads to obtain the final library. After library construction, preliminary quantification was conducted using a Qubit 2.0 Fluorometer, and libraries were diluted to 1.5 ng/µl. The insert size of the libraries was assessed using an Agilent 2100 bioanalyzer. Libraries underwent quantitative real-time PCR (qRT-PCR) for accurate quantification of effective library concentration (library concentration > 1.5 nM) to ensure library quality. Then, the cDNA libraries were sequenced using the Illumina Nova Seq platform. The obtained raw data were filtered to eliminate reads with adapters, poly-N (which denotes undetermined base information), or low-quality reads (reads with Qphred ≤ 5 bases that accounted for more than 50% of the entire read length). HISAT2 v2.0.5 was used to compare the reads with the reference genome. The differential expression analysis was conducted using DESeq2 software (version 1.20.0). The Benjamin-Hochberg procedure (B-H FDR) was used with a screening criterion of | log2 (Fold Change) | ≥ 1 and p-adjusted (padj) ≤ 0.05 was used to detect DEGs between the treatment and control groups. The statistical power of these experimental data was as follows: RNASeqPower of CON vs V15 is 0.98984, and CON vs V30 is 0.99862. Correlation coefficients within and between group samples were also calculated based on FPKM values for all genes in each sample and produced as heatmaps (Fig. S1).

GO and KEGG analyses of DEGs

The ClusterProfiler software was used to analyze the differential gene sets for GO function enrichment and KEGG pathway enrichment. Gene Ontology (GO) enrichment analysis was used to determine the protein functions of the genes, including biological process (BP), cellular component (CC) and molecular function (MF), with padj ≤ 0.05 as the threshold of significance for enrichment. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was used to determine the gene function, genomic information, and functional information of the DEGs, with padj ≤ 0.05 serving as the threshold for significant enrichment.

PPI network construction

The STRING protein interaction database (http://string-db.org) was used to analyze the protein–protein interaction (PPI) networks, which were then visually edited with Cytoscape software (version 3.10) and analyzed via the degree method to identify hub nodes in the networks. Hub genes may include genes that play relatively important roles in biological processes.

RT-PCR

RT-PCR was used to analyze several key genes screened in the RNA-seq results (n = 5 from each group). Total RNA was extracted using the TRIzol method (Invitrogen, USA) and all steps were performed at room temperature (20 °C–25 °C). The RNA was reverse transcribed to obtain cDNA using the PrimeScript RT Kit and the genomic DNA scavenger (Takara, Japan), and the experimental procedure was carried out in full accordance with the instructions of the kit. The cDNA was extracted using the SYBR Green PCR kit (Takara, Japan) on a StepOnePlus Real-Time PCR system (Applied Biosystems, Waltham, MA, USA) for RT-PCR. Table 1 lists the primer sequences that were designed.

Statistical analysis

All the statistical analyses were performed using GraphPad Prism 8 (GraphPad Software, La Jolla, CA, USA). One-way ANOVA was used to determine significant differences among the three groups. P-values less than 0.05 were considered to indicate statistically significant. All the data were expressed as the mean ±SD. Experimental manipulation and data analysis were carried out by two separates people.

VCD affects sex hormone secretion and follicular growth in mice

Previous research has shown that the difference in the ovaries of VCD-model mice occurs 30 days after the end of drug administration (Cao et al., 2020). The other commonly used mouse models of POI, such as cyclophosphamide (CTX), usually show difference in the ovaries on the 15th day (Dai et al., 2023). Based on these previous results, the present study used the 15th and 30th days as the observation time points and collected the ovaries and serum at those time points for subsequent experiments (Fig. 1A). The results showed AMH and E₂ levels were significantly lower in both the V15 and V30 group than in the CON (Fig. 1B). FSH tended to increase in the experimental groups, but the difference was not statistically significant (Fig. 1B). Next, the numbers of primordial follicles, primary follicles (PFs), secondary follicles (SFs), antral follicles (AFs), and corpora luteum (CL) was counted in the three groups (Figs. 1C and 1D). The results showed the number of follicles of all five species was significantly reduced in both V15 and V30. During the modeling period, the estrous cycles of each group of mice were monitored, revealing that the estrous cycle rhythm was disrupted in the mice injected with VCD (Fig. S2 and Table S1). These results are consistent with the clinical features of POI. Based on the differences in the apparent morphological characteristics of the mice in the two experimental groups, these results suggest that the damage caused by VCD to the ovaries of the mice intensified over time.

Transcriptomics-based identification of DEGs

A transcriptome analysis was performed to reveal the differences in transcript profiles between the control group and experimental groups. A principal component analysis (PCA) revealed notable differences between samples, with tightly clustered samples from each group being relatively separated from one another, indicating significant distinctions in gene expression profiles between the groups (Fig. 2A). The control group was used as a reference point to identify disparities in gene expression. Then, a heatmap (Fig. 2B) and volcano plots (Figs. 2C and 2D) were generated to visualize the differences between the groups. The results showed that nine genes were upregulated and 458 genes were downregulated in the V15 group in contrast to the control group (Fig. 2C). Similarly, 488 genes were upregulated and 1,861 genes were downregulated in the V30 group (Fig. 2D). This finding suggests that VCD affected gene expression, with a significant difference in the expressed gene profiles of the transcriptomes observed between the groups. The complete results of the difference analysis are shown in Tables S2 and S3.

Enrichment analysis of DEGs

Enrichment analyses of the differentially expressed genes were then performed. The top 10 significant enrichment results of the GO terms and the KEGG pathways in the V15 and V30 groups are shown in Figs. 3 and 4 (if there were fewer than 10 significant results, the first 10 results are displayed). In group V15, too few genes were upregulated, so no significant enrichment results were obtained (Figs. 3A and 3C). The GO and KEGG pathway enrichment analyses results showed the downregulated DEGs were related mainly to the motor protein family and mitosis-related molecular activities, such as cilium movement and axoneme assembly (Figs. 3B and 3D). Richer results were obtained from the GO and KEGG pathway analyses for group V30: the downregulated DEGs were associated with the cell cycle; DNA damage repair-related pathways such as Fanconi anemia; homologous recombination, and steroid anabolism, including meiosis, ovarian steroidogenesis, and oocyte homologous recombination (Figs. 4B and 4D). The upregulated DEGs were significantly associated with steroid biosynthesis (Figs. 4A and 4C). There were common pathways between the two groups, but they also had unique differences; more results were obtained at V30, which may indicate that the effects of VCD on the mouse ovary increase over time.

PPI Network Analyses of DEGs

The differential gene protein interaction networks were analyzed using the interactions in the STRING database and two PPI networks based on the differential genes of V15 vs CON and V30 vs CON were produced using Cytoscape_v3.10.0. The hub nodes were identified from this PPI network based on node degree, and the results were arranged in descending order of node degree (Table S6). For V15 vs CON, the top three key genes were Oog3, Oog4, Oog2 (Fig. 5A). And according to the analysis results, C87414, Pramef1, C87499, Gm10436, Gm13084, Gm13023, Gm13103, Oog1, C87977, Pramef20, C87414-2, Gm3106, and D5Ertd577e share the same node degree. Oog3, Oog4, Oog2 and Oog1 belong to a family of genes, all of which have been found to be located on chromosome 4 and are expressed throughout oogenesis and preimplantation embryonic development (Dadé et al., 2003). No studies related to female reproduction are available for the remaining genes. For V30 vs CON, the five genes with the most commonalities were KIF11, BUB1B, CDC20, BUB1, and CENPE (Fig. 5B), which all play important roles in meiosis in oocytes (Duesbery et al., 1997; Chen et al., 2010; Liu et al., 2010; Kovacovicova et al., 2016). This PPI network construction and analysis provides insights into the possible causative agent of VCD.

Expression levels of target genes in the ovaries

To validate the RNA sequencing results and explore the possible genes and pathways related to altered ovarian function in mice during VCD modeling, 10 key genes were selected: BUB1, BUB1B, CDC45, CYP11A1, CYP17A1, HSD17B1, HSD3B6, FSHR, BRCA2, and FANCA. These genes differed significantly between the groups according to the results of the difference analysis (Tables S2 and S3). Based on the GO and KEGG results, these 10 genes were clustered in pathways with strong physiological relevance to POI, such as the cell cycle, DNA damage and steroid anabolism. The results of quantitative real-time PCR were consistent with the results of the RNA-seq analysis. BUB1, BUB1B, and CDC45 were mapped to the oocyte meiotic cycle pathway. As shown in Figs. 6A, 6B and 6C, the expression levels of these genes showed a downward trend and they were statistically significant in the V30 group. BRCA2 and FANCA are key genes of the Fanconi anemia pathway, and their expression was also significantly lower in V30 than in the control mice (Figs. 6D and 6E). Notably, the expression trends of key genes involved in ovarian steroid hormone synthesis were different between V15 and V30 mice, with the expression levels of CYP17A1 and HSD3B6 increasing in the V30 group (Figs. 6F and 6G), and the expression level of HSD17B1 gradually decreasing, reaching significance at V15 (Fig. 6H). CYP11A1 and FSHR did not significantly differ between the groups (Figs. 6I and 6J). These findings underscore the dynamic alterations in gene expression patterns associated with VCD-induced ovarian dysfunction, shedding light on potential molecular mechanisms underlying the pathogenesis of POI.