Comprehensive circular RNA expression profiling with associated ceRNA network in orbital venous malformation

Patients and samples

Nine samples were taken from OVM patients diagnosed with pathological verification at the Ninth People’s Hospital, Shanghai Jiao Tong University School of Medicine, between January 2017 and June 2019 (Appendix 1), as previously reported [9]. Nine normal orbital vascular specimens were collected from patients with orbital trauma or cosmetic surgery (Appendix 2), as previously reported [9]. Specimens were only used with informed consent from the patients, and this study was approved by the research ethics committee of the Ninth People’s Hospital, Shanghai Jiao Tong University School of Medicine.

RNA extraction, library preparation, and sequencing

Total RNA was isolated from the OVM samples of four patients and four normal orbital vascular tissues using a mirVana miRNA Isolation Kit (Ambion, Austin, TX) following the manufacturer’s instructions. Briefly, the tissue samples were lysed in a denaturing lysis solution, and then the RNA was extracted with Acid-Phenol: Chloroform to remove most of the other cellular components. Further purification was performed on a glass-fiber filter to yield either total RNA or miRNA-enriched fractions. Finally, solutions formulated specifically for miRNA retention was used to avoid the loss of small RNA. Then, cDNA libraries were constructed using a TruSeq Stranded Total RNA Kit with Ribo‐Zero Gold (Illumina, San Diego, CA) following the manufacturer’s protocol. Quality validation of the libraries was performed using an Agilent 2100 Bioanalyzer (Agilent Technologies), and high-throughput RNA-seq of these libraries was conducted on the Illumina sequencing platform (HiSeqTM 2500).

For miRNA sequencing, a Small RNA Sample Prep Kit (Illumina) was employed to construct libraries following the manufacturer’s instructions. The validity of the library quality control was verified by an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA), and sequencing was conducted on a HiSeq X Ten.

Differential expression RNA analysis

We used DESeq software packages [16] to analyze the differentially expressed mRNAs, circRNAs, and miRNAs. The cutoff criterion of differentially expressed RNAs was |log2 FC (fold change) |≥1 and p<0.05. Then, hierarchical cluster analysis was conducted using the pheatmap package (Version 1.0.8) [17] with these differentially expressed RNAs.

qRT-PCR

Total RNA was isolated from OVM samples and normal orbital vascular tissues using the mirVana miRNA Isolation Kit (Ambion), following the manufacturer’s protocol. Thereafter, the RNA was reverse transcribed into cDNA using the PrimeScript RT reagent kit (Takara Bio, Shiga, Japan). Then, qRT-PCR was performed using SYBR Green PCR master mix (Life Technologies, Carlsbad, CA) and an RT-PCR system (Applied Biosystems, Waltham, MA). The housekeeping gene ACTB was the control for PCR product quantification and normalization. The primers for circRNAs and mRNAs are listed in Table 1.

Construction of circRNA‐mRNA coexpression network

The coexpression analysis of the circRNA-mRNA network was constructed according to the differential expression levels of mRNAs and circRNAs in OVM and normal control samples. We used the Pearson’s correlation coefficient of mRNAs and circRNAs >0.80 and p<0.05 to construct the coexpression network and Cytoscape software (Version 3.2.1) to visualize it.

Construction of a circRNA-associated ceRNA network

The exploration of potential circRNA-associated ceRNA networks was analyzed according to the obtained differentially expressed mRNAs, circRNAs, and miRNAs. CircRNA–miRNA pairs and miRNA–mRNA pairs were predicted using the miRanda tools. These two pairs sharing the same miRNAs were identified as candidate circRNA–miRNA–mRNA competing interactions. Expression correlations between circRNAs, miRNAs, and mRNAs with a threshold of p value <0.05 and Pearson’s correlation coefficient >0.85 or < −0.85 were filtered out. Thereafter, the circRNA-associated ceRNA network was constructed using Cytoscape software (Version 3.2.1).

GO and KEGG analysis

The functions of the differentially expressed mRNAs identified from the circRNA-associated ceRNA network were analyzed using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. GO analysis used the biologic process (BP), cellular component (CC), and molecular function (MF) to mark genes. KEGG analysis provided annotation information for gene signal transduction and disease pathways. The GO terms and KEGG pathways with p<0.05 were considered to be enriched with statistical significance.

Representative clinical characteristics of OVM patients

In our study, nine OVM tissue samples and nine normal control samples were included; of these, four OVM tissue samples and four normal control samples were sequenced and the other five OVM tissue and five normal control samples were used to validate the sequencing by qRT-PCR. Representative pictures of the patient’s appearance are shown in Figure 1A,B. In magnetic resonance imaging (MRI), hypointensity is apparent in T1-weighted images (T1WIs) and hyperintensity in T2-weighted images (T2Wis; Figure 1C,D). Normal control tissues were collected from orbital trauma or cosmetic surgery, as shown with OVM samples in Figure 1E,F. In addition, hematoxylin and eosin (H&E) staining (Appendix 3) and immunohistochemical staining (Appendix 3) of a vascular endothelial cell marker, CD34, showed that both OVM and control samples were vascular tissue.

Identification and characterization of circRNAs

The expression profile of circRNAs in OVM and normal control tissues was detected by RNA-seq, and 9,375 circRNAs were identified. Among these circRNAs, 5,106 were reported in circBase, and 4,269 identified circRNAs were newly discovered (Figure 2A,B). Sequence length analysis revealed that these circRNAs ranged in length from less than 200 nucleotides (nt) to more than 2,000 nt, and most were 200–500 nt or over 2,000 nt (Figure 2C). The circRNA transcripts were broadly distributed throughout the genome, with more on chromosomes 1, 2, and 3 (Figure 2D). In addition, exonic circRNAs accounted for the largest type (Figure 2E), and the number of exons per circRNA mainly ranged from 1 to 6 (Figure 2F). We used reads per million (RPM) values to estimate the expression level of the circRNAs, and there were no significant differences in each sample (Figure 3A,B).

Identification of differentially expressed RNAs

To identify the potential biological role of circRNAs in OVM, we analyzed differentially expressed mRNAs, circRNAs, and miRNAs with thresholds of |log2 FC| ≥1 and p<0.05. RNA-seq analysis identified 3449 aberrantly expressed mRNAs in OVM, including 2,175 upregulated and 1,274 downregulated mRNAs (Figure 4A, Appendix 4). Meanwhile, 45 upregulated and 144 downregulated circRNAs were found in OVM tissues compared with normal controls (Figure 4B, Appendix 5). In addition, 156 miRNAs were upregulated and 168 miRNAs were downregulated in the OVM samples (Figure 4C, Appendix 6). Hierarchical clustering analysis showed that differentially expressed mRNA, circRNA, and miRNA expression patterns among samples were distinguishable (Figure 4D,E).

Validation of differentially expressed mRNAs and circRNAs by qRT-PCR

We performed qRT-PCR to verify the reliability of the differentially expressed circRNAs and mRNAs identified in the OVM tissues and normal controls by RNA-seq. Five mRNAs (CALD1, RHOA, AES, CTSD, and LAPTM5), and three circRNAs—circRNA_2189 (hsa_circ_0000437 officially), circRNA_3654 (hsa_circ_0000745 officially), and circRNA_3192 (hsa_circ_0003856 officially), were randomly selected. The results revealed that CALD1, RHOA, and AES were downregulated in OVM samples and that CTSD and LAPTM5 were upregulated in OVM samples, which is in accordance with the RNA-seq results (Figure 5A,C). The qRT-PCR results of circRNA_2189, circRNA_3654, and circRNA_3192 were in line with the RNA-seq data (Figure 5B,D), suggesting that the differentially expressed circRNAs identified via RNA-seq were reliable.

Construction of circRNA-associated ceRNA networks

To investigate the role of circRNAs in OVM, the expression correlations between circRNA and mRNA and ceRNAs (circRNA–miRNA–mRNA) were identified according to the differentially expressed RNAs, including 54 circRNAs, 1,566 mRNAs, 107 miRNAs, and bioinformatic prediction results. We first constructed a coexpression network of the top 100 circRNA–mRNA interactions by Cytoscape 3.2.1 (Figure 6). Then, two circRNA-associated ceRNA networks—circRNA_8224 (upregulated in OVM)–miRNAs (downregulated in OVM)–mRNAs (upregulated in OVM) and circRNA_1678 (downregulated in OVM)–miRNAs (upregulated in OVM)–mRNA (downregulated in OVM)—were constructed (Figure 7). These ceRNA network results may provide a new underlying mechanism for OVM.

Functional analyses of dysregulated mRNAs in the ceRNA network

GO and KEGG pathway analyses further explored the potential function of mRNAs involved in the ceRNA network. Based on the GO analysis results, mRNAs with upregulated expression levels were significantly enriched in BPs, such as extracellular matrix organization and positive regulation of actin nucleation; CCs, such as the plasma membrane and adherens junction; and MFs, such as lysophosphatidic acid receptor activity and cell adhesion molecule binding. The downregulated mRNAs were mainly associated with BPs, such as mesenchymal migration and muscle contraction; CCs, such as the Z disc and stress fiber; and MFs, such as actin binding and cadherin binding (Figure 8A-B, Appendix 7 and Appendix 8). According to the results of the KEGG pathway analyses, we identified upregulated mRNAs that participate in pathways, such as mucin-type O-glycan biosynthesis, glycosaminoglycan degradation, and the PI3K-Akt signaling pathway. The downregulated mRNAs were involved in pathways, such as regulation of the actin cytoskeleton, the calcium signaling pathway, and the cAMP signaling pathway (Figure 8C-D, Appendix 9 and Appendix 10). These functional analyses suggested that the circRNA-associated ceRNA network may have multiple regulatory functions in OVM pathogenesis. In addition, we found that most of the genes in the mucin-type O-glycan biosynthesis, glycosaminoglycan degradation and PI3K-Akt signaling pathway were upregulated (Appendix 11), and representative genes in these signaling pathways were verified by qRT-PCR (Figure 8E). We further applied immunohistochemical (IHC) staining to p-PI3K and p-Akt, and the results demonstrated that p-PI3K and p-Akt presented stronger signals in OVM samples (Figure 8F). These data indicate that these pathways may account for OVM.