Circular RNAs in endometriosis analyzed through RNA sequencing and bioinformatics for expression profile

Participants and sample collection

In this study, a total of three individuals with endometriosis (case group) aged between 25 to 35 years in the proliferative phase, alongside three individuals with normal endometrium devoid of endometriosis (control group) aged between 26 to 38 years in the proliferative phase, were included. Tissue specimens from snap-frozen cyst walls of ovarian endometriosis and corresponding normal endometrial tissues were collected for RNA sequencing analysis and qRT-PCR, samples are quickly frozen in liquid nitrogen after collection and stored in liquid nitrogen tanks. All participants had consistent menstrual cycles and had abstained from hormonal therapy for a period of six months prior to the procurement of the samples. Additionally, individuals in the control group showed no signs of endometrial tumors and lacked a histological diagnosis of adenomyosis. Specimen diagnoses were confirmed through histopathology, and menstrual cycle phases were verified using both the last menstrual period and histological assessments. Authorization for the conduct of this research was provided by the institutional ethical review committee of the first People’s Hospital of Wuhu City (approval number: YYLL20220004). Full and voluntary written consent was obtained from participants before all tissue samples were obtained, and the study strictly adhered to the ethical guidelines set out in the Declaration of Helsinki and its subsequent revisions or similar ethical codes. The statistical power of this experimental design, calculated in RNASeqPower is 0.90.

RNA isolation and quality control

The frozen tissue samples underwent total RNA extraction employing the TRIzol reagent (Invitrogen, Waltham, MA, USA) in accordance with the manufacturer’s prescribed protocol. Subsequently, the RNA concentration in each sample was quantified utilizing the NanoDrop 2000 spectrophotometer (NanoDrop, Thermo Fisher Scientific, Waltham, MA, USA). Furthermore, the integrity of the extracted RNA was validated through 1% formaldehyde-based denaturing gel electrophoresis.

circRNA sequencing

The responsibility of circular RNA (circRNA) sequencing was assigned to the Gene Denovo Company, situated in Guangzhou, China. The initial phase involved preparing circRNA sequencing libraries, which kicked off with isolating total RNA from individual endometriosis tissue samples. Prior to library construction, the total RNA mixture underwent depletion of both ribosomal RNAs (rRNAs) and linear RNAs. This was accomplished using the Ribo-Zero Gold rRNA Removal Kit sourced from Illumina (San Diego, CA, USA) and the Ribonuclease R Kit procured from Epicentre (Madison, WI, USA), respectively. Following this, double-stranded cDNA was synthesized employing the polymerase chain reaction (PCR) technique. Subsequently, this cDNA underwent a series of processing steps, including adapter ligation, A-tailing, and end-repair, utilizing the NEBNext poly(A) mRNA Magnetic Extraction Module from New England Biolabs (Ipswich, MA, USA). Then, a 12-cycle PCR amplification was carried out, and the resulting products were purified with AMPure XP beads manufactured by Beckman (Brea, CA, USA). The quality assessment of the sequencing libraries was performed using the DNA 1000 Assay Kit from Agilent Technologies. The circRNA sequencing process was executed on the Illumina NovaSeq 6000 platform.

The raw sequencing data underwent processing to exclude those adapter reads that were of low quality, resulting in clean reads. These cleaned reads were then aligned to the rRNA database using Bowtie2. The reads that successfully mapped to rRNA were subsequently discarded, and the remaining reads were mapped to the reference genome through the utilization of TopHat2 version 2.0.3.1214. Finally, the unmapped reads were selected for the purpose of identifying circRNAs.

Differentially expressed circRNA analysis

Following QC procedures, the raw sequencing data underwent filtration to remove concatenated sequences and excessively short fragments. Following the trimming process, the refined dataset was aligned to the reference genome (GRCh37/hg19) by employing the STAR software tool. This alignment step involved mapping the cleaned and processed data to the specified reference genome, utilizing the capabilities of STAR. The subsequent step involved differential expression analysis of the circRNA-seq data, which was executed within the R statistical environment (version 4.2.0) utilizing the edgeR package. CircRNAs that demonstrated a significant fold change (either ≥ 2.0 or ≤ 0.5), a p-value of less than 0.05, and an FDR q-value of less than 0.05, when compared between the case and control groups, were designated as differentially expressed circRNAs. Annotation of these circRNAs was achieved by cross-referencing with circBase, with those that remained unannotated being categorized as novel circRNAs. Unmapped sequencing data to the reference genome were further scrutinized for circRNAs using CIRCexplorer2 software to detect reverse splicing events. Visualization of these differentially expressed circRNAs was achieved using CIRCOS visualization software in circular format.

qRT-PCR

To confirm the expression levels of the circRNAs identified through sequencing, qRT-PCR validation was undertaken. A total of six circRNAs, consisting of the three most upregulated and three most downregulated circRNAs annotated in circBase, were chosen for validation. The primers required for this validation process were procured from Sangon Biotech, located in Shanghai, China. For the purpose of qRT-PCR analysis, The sample is ground in liquid nitrogen and placed in Trizol for RNA extraction, the entire process is carried out in an environment containing DNase and RNase, including subsequent experiments. RNA concentration is measured using a Nanodrop (Thermo Fisher Scientific, Waltham, MA, USA). cDNA synthesis was carried out using the All-in-One First-Strand cDNA Synthesis Kit offered by GeneCopoeia Inc (Santa Cruz, CA, USA), with one µg of total RNA as the starting material, the total reaction system is 20 µl, which includes one µg RNA, four µl RT mix, one µl Priming oligonucleotide, and then filled up to 20 µl with DEPC water. The machine program is set as follows: 25 °C for 10 min, 55 °C for 40 min, and 85 °C for 5 min. Subsequently, the qRT-PCR assays were conducted utilizing the All-in-One qPCR Mix from GeneCopoeia Inc (Rockville, MD, USA) on the ABI 7500HT qRT-PCR System manufactured by Applied Biosystems (Foster City, CA, USA). The thermal cycling protocol for the qRT-PCR encompassed an initial denaturation step at 95 °C for 30 s, followed by 40 cycles of denaturation at 95 °C for 5 s and annealing/extension at 60 °C for 30 s, the reaction system is as follows, with a total volume of 20 µl, including 10 µl of SYBR Green mix, one µl of 10 µM primers each, two µl of cDNA, and finally topped up with DEPC water to 20 µl. The threshold cycle (Ct) value represented the point at which a positive amplification signal was detected. GAPDH served as the internal reference gene for normalization. Expression levels were quantified employing the 2-ΔΔCt method, all checks repeated three times, and the primer sequences are presented in Table 1.

GO and KEGG pathway analysis

Utilizing the DAVID tool, we delved into the potential functions and roles of the differentially expressed circRNAs and their corresponding parental genes. We executed a Gene Ontology (GO) enrichment analysis, encompassing biological processes (BP), cellular components (CC), and molecular functions (MF), drawing upon data from the Gene Ontology database. Additionally, to identify the significantly enriched pathways associated with the parental genes of these circRNAs, we performed a KEGG pathway enrichment analysis, leveraging the extensive information contained within the KEGG database.

CircRNA-miRNA-mRNA network analysis

To delve into the modus operandi of circular RNAs (circRNAs) serving as decoys for microRNAs (miRNAs), thereby exerting an indirect influence on the synthesis of their respective messenger RNAs (mRNAs), we utilized miRanda software (accessible at http://mirtoolsgallery.tech/mirtoolsgallery/node/1055) and RNAhybrid software (located at https://bibiserv.cebitec.uni-bielefeld.de/rnahybrid/) to anticipate the miRNAs and mRNAs that are the targets of the most perturbed circRNAs. Subsequent to this anticipation, we employed Cytoscape software (version 3.9.0, hosted at https://cytoscape.org/) to visually represent the intricate interplay within the intricate circRNA/miRNA/mRNA networks, reshuffling word order and incorporating synonyms for clarity.

Statistical analysis

Data analysis was conducted utilizing SPSS Statistics (Version 20, Chicago, IL, USA). The continuous data were averaged, and both addition and subtraction operations were applied. Additionally, the standard deviation was calculated. To examine the disparities between two groups, a t-test was executed. The outcomes unveiled a noteworthy distinction between the two groups (p < 0.05).

Quality control of circRNA sequencing

Using high-throughput deep sequencing, we first examined the circRNA expression profiles in six paired tissues collected from patients with ovarian endometriosis and control subjects. Among the host genes of these candidate circRNAs, 97.5% originated from exonic regions, while the remainder were located in intronic and unannotated regions (Fig. 1A). Next, we conducted Pearson correlation analysis to assess the correlation among all samples. As illustrated in Fig. 1B, we observed highly correlated gene expression patterns of circRNAs in both the control group (C1, C2, C3) and the case group (M4, M5, M6), respectively. The correlation heat map of circRNA expression levels across different samples also revealed a strong correlation between the internal samples of the control group and the case group (Fig. 1C). Subsequently, PCA analysis was performed, revealing a clear separation between the case group and the control group (Fig. 1D).

CircRNA expression profile in endometriosis

A scatter plot (depicted in Fig. 2A) was employed to compare the circRNA expression profiles between the ovarian endometriosis group and a cohort of three corresponding normal tissues. Statistically significant differentially expressed circRNAs were identified based on stringent criteria: fold changes ≥ 2.0 or ≤ 0.5, p-values < 0.05, and FDR q-values < 0.05. This analysis revealed a total of 679 circRNAs exhibiting marked differential expression, with 371 circRNAs being significantly upregulated and 308 circRNAs being noticeably downregulated. These findings were further illustrated in volcano plots (Fig. 2B) and a clustered heatmap (Fig. 2C), showcasing the upregulation and downregulation patterns. Additionally, the distribution pattern of these differentially expressed circRNAs across the human chromosomes was visualized through a heatmap (Fig. 2D).

Validation of differentially expressed circRNAs

To verify the outcomes derived from RNA-sequencing (RNA-seq), we selected the three most prominently increased and decreased circRNAs for quantitative reverse transcription-polymerase chain reaction (qRT-PCR) assessment on an independent cohort with 50 endometriosis patients. Our findings revealed a statistically noteworthy upregulation of hsa_circ_000005, hsa_circ_000002, and hsa_circ_000016 within ovarian endometriosis tissues (Figs. 3A–3C). Conversely, the expression patterns of hsa_circ_000011, hsa_circ_000010, and hsa_circ_000021 exhibited significant downregulation in the identical ovarian endometriosis samples (Figs. 3D–3F), thereby reinforcing the initial RNA-seq discoveries.

GO and KEGG pathway analyses of changed circRNAs

Using GO and KEGG analyses, the predicted potential function and association of the circRNAs that showed differential expression and their corresponding parental genes were examined. The enriched GO terms for the differentially expressed circRNAs were ranked, highlighting the top eight terms. Notably, the enriched GO terms pertaining to biological processes primarily focused on muscle attachment and the regulation of cAMP-mediated signaling (Fig. 4A). Regarding cellular components, the enriched GO terms were predominantly linked to the cytoskeleton and cytoplasmic vesicles (Fig. 4B). In terms of molecular functions (MF), the enriched GO terms were primarily associated with microtubule binding and the structural constituents of the cytoskeleton (Fig. 4C). Additionally, the enriched scatter diagram illustrated the top eight pathways identified through KEGG pathway enrichment analysis (Fig. 4D), indicating potential associations of the differentially expressed genes with Vitamin B6 metabolism and the synthesis, secretion, and function of Parathyroid hormone.

CircRNA-miRNA-mRNA network

The circRNA–miRNA–mRNA regulatory network was analyzed using miRanda and RNAhybrid software. The most upregulated circRNA, hsa_circ_000005, and the most downregulated circRNA, hsa_circ_000011, were prioritized in the network. Both circRNAs were predicted to interact with multiple miRNAs, with at least four miRNAs binding to each circRNA. The potential functional network involving hsa_circ_000005 and hsa_circ_000011 is illustrated in Fig. 5.