Profiling and Bioinformatic Analyses Indicate Differential circRNA and miRNA/isomiR Expression and Interactions

As a novel class of noncoding RNAs, circular RNAs (circRNAs) have been reported to play a role in various biological processes. Some circRNAs may serve as microRNA (miRNA) sponges, regulating transcription or splicing. Herein, we investigated the expression profiles and interactions of miRNAs/isomiRs and circRNAs in male patients with esophageal cancer. We found that some miRNA genes generated two deregulated miRNA products (miR-#-5p and miR-#-3p), and these products were consistently abnormally expressed. Some circRNAs were predicted to be miRNA sponges for specific miRNAs. Some of these typically showed opposing expression patterns in cancer tissues: one upregulated and the other downregulated. Although fewer miRNAs were predicted to interact with circRNAs, the number of predicted interactions would be substantially increased if detailed isomiRs were involved. High sequence similarity across multiple isomiRs suggested that they might interact with circRNAs, similar to the interaction of homologous miRNAs with circRNAs. At the isomiR level, due to the characteristics of the sequences and expression patterns involved, the cross-talk between different ncRNAs is complicated despite simplification of the isomiRs involved through clustering. We expect that our results may provide methods for further study of the cross-talk among ncRNAs and elucidate their biological roles in human diseases.


Introduction
Esophageal cancer is a common cancer with high cancer mortality [1]. The impact of this disease is most severe in southern Africa, eastern Africa, and eastern Asia [2]. In China, the incidence of esophageal cancer in some areas is high [3]. Noncoding RNAs (ncRNAs), including microRNAs (miRNAs), long noncoding RNAs (lncRNAs), and circular RNAs (circRNAs), may play crucial roles in the regulation of the occurrence and development of many cancers. Indeed, some ncRNAs have been proposed as biomarkers for the diagnosis of many cancers [4,5].
Among the ncRNAs, miRNAs have been shown to play crucial regulatory roles in a variety of biological processes, particularly those related to pathological and physiological pathways [6,7]. An increasing number of studies have shown that miRNA loci can generate miRNA variants (isomiRs) [8][9][10][11][12][13][14][15]. It has further been shown that these isomiRs are functional molecules that have expression and sequence heterogeneities [12,[16][17][18][19][20][21][22][23][24][25]. Some small RNAs have been shown to have important regulatory roles in occurrence and development of esophageal cancer. For example, miR-330-3p regulates PDCD4 and might have an oncogenic role in esophageal squamous cell carcinoma [26], circuit of miR-31, SOX4, and EZH2 increases tumor progression in invasive esophageal carcinomas [27], and miR-10b may target KLF4 and thus promote the migration and invasion of tumorous cells [28]. These findings suggest that the small RNAs play important roles in esophageal cancer. However, few studies have investigated isomiRs, despite the fact that the canonical miRNA sequence is only one specific variant of multiple possible isomiRs [14,29].
Another class of recently discovered animal RNAs, cir-cRNAs, forms covalently closed continuous loops lacking 5 2 BioMed Research International and 3 ends; these are also potential regulatory molecules [30]. Similar to miRNAs, circRNAs are thought to be crucial to cancer development and have thus been widely studied: a series of dysregulated circRNAs have been identified as important to the development of radiation resistance [31], circ 002059 is significantly downregulated in gastric carcinoma and may be a stable and novel biomarker [32], and circRNA 001569 may target miR-145 and thus contribute to colorectal cancer cell proliferation and invasion [33].
An understanding of miRNA-circRNA interactions is important for understanding the roles of ncRNAs in general, especially with respect to cross-talk among ncRNAs. Some circRNAs have been shown to have important regulatory roles in miRNA activity via sponging miRNAs. For example, the circRNA CDR1 acts as a miR-7 super-sponge [30,[34][35][36]. Similar to miRNAs, the expression profiles of circRNAs are specific to different cell types and tissues [30,37], indicating potential regulatory functions and roles in cellular processes. Some studies have suggested that circRNAs may be potential biomarkers for disease diagnosis [32,38], based on their potential biological roles and on their stable features, such as closed RNA loops.
Although miRNAs and circRNAs have been studied in many cancers, few have been investigated at the isomiR level, despite the fact that multiple isomiRs are the detailed products from miRNA locus. Herein, we aim to understand the expression profiles of miRNA/isomiR and circRNA in esophageal cancer using deep sequencing data, focusing on their interactions and potential roles in cancer development. We also propose an approach for the analysis of miRNA/isomiR-circRNA cross-talk.

Patient Samples.
A total of 5 paired tissue samples (tumor tissue and matched adjacent nontumor normal tissue) were taken from 5 male patients (age: 64.2 ± 2.95) with squamous esophageal cancer (in all cases, the distance between the tumor and nontumor normal tissue was >5 cm, Table S1). The 5 patients underwent surgical resection without prior radiotherapy or chemotherapy, in the Department of Cardiothoracic Surgery at the Third Affiliated Hospital of Soochow University, Changzhou, China. All samples from patients were HPV-negative and HPV has been suggested to play an important role in the development of esophageal cancer [39]. Informed consent was obtained from each subject. This study was approved by the ethics committee of Soochow University.

High-Throughput Sequencing and Analysis.
Total RNA was extracted from each sample using Trizol (Invitrogen, USA). Purified RNA was quantitatively and qualitatively analyzed with NanoDrop (Thermo Company, USA). CircRNA and small RNA sequencing libraries were constructed with HiSeq 2000 sequencing platform following the manufacturer's instructions.
CIRCexplorer software (version 1.1.7) [40] was used to identify circRNAs based on the fusion gene. After trimming adaptors, the small sequence fragments were mapped onto the human genome using Bowtie [41]. miRNAs/isomiRs were identified based on the canonical miRNA and pre-miRNA sequences annotated in miRBase version 21.0 (http://www.mirbase.org/) [42]. IsomiRs were identified according to the location of the canonical miRNA in pre-miRNA with ±4 nt [43], and further analyses were performed from the miRNA and isomiR levels, respectively. Differentially expressed circRNA and miRNA/isomiR profiles were obtained with edgeR, particularly those that were highly expressed. In the small RNA analysis, the canonical miRNA sequence and isomiRs from the specific miRNA locus were analyzed, respectively.

miRNA/isomiR-circRNA Interaction
Analysis. The potential binding sites of abnormally expressed miRNAs were predicted using MiRanda [44] based on deregulated cir-cRNAs. The potential relationships of different RNAs were constructed using Cytoscape version 3.4.0 [45]. It is possible to predict specific interactions at the isomiR levels, as well as based on the canonical miRNA sequence, or via clustering analysis of multiple isomiRs based on their close seed sequences (close sequence relationships can be obtained based on miRNA locus or independent miRNA locus) ( Figure 1). Potential miRNA-circRNA interactions were also predicted and obtained from starBase v2.0 [46].

Functional Analysis of Differentially Expressed miRNAs and isomiRs.
For the top differentially expressed miRNAs (the canonical miRNA and isomiRs), we performed functional analysis using CapitalBio Molecule Annotation System V4.0 (MAS, http://bioinfo.capitalbio.com/mas3/). Target mRNAs of the miRNAs were either identified from experimentally validated annotated in miRTarBase [47] or predicted with TargetScan [48]. Based on the same predictive methods, the target mRNAs of the isomiRs were predicted after clustering multiple isomiRs using seed sequences. We also further discussed the potential relationships among miRNAs/isomiRs, mRNAs, and circRNAs based on their functional relationships.

Statistical Analysis.
Differentially expressed circRNAs and miRNA/isomiRs were obtained using edgeR based on paired t-test results with FDR correction: that is, when adj value < 0.05 and the |log 2 (FC)| ≥ 2.0, we considered the RNA a deregulated species. Relative expressions were described using mean ± SD ( ± SD) based on normalized sequence counts. Specifically, the relative expression of the isomiR was further estimated based on percentage of multiple isomiRs in the specific miRNA locus, and those abundantly expressed isomiR species were specifically concerned in relevant analysis. Distributions of different small RNAs were estimated using Venny [49].  on the detected circRNAs from sequencing data (Figures S1 and S2). Among the small ncRNAs with abundant detection, many miRNA/isomiRs showed significant abnormal expression ( Figure S3), and the isoform profiles were consistent with the relative expression levels of the miRNA locus.

Deregulated Small RNAs Contributed to Basic Biological
Progress Based on Functional Analysis. Further functional analysis suggested that deregulated small RNAs might contribute to some basic biological progress via the negative regulation of target mRNAs ( Figure 2). For example, these small RNAs might be involved in processes such as singlemulticellular organism process, system development, cell development, and cell differentiation. Although the role of miRNAs in these pathways was indirect, their flexible regulatory contributions might be important to many relevant biological processes.

Predicted miRNA-circRNA Interaction Based on the
Canonical miRNA. Based on known mature miRNAs and predicted circRNAs, we identified miRNA-circRNA pairs that potentially interact. Of these, some miRNAs were homologous with high sequence similarity (some nucleotides might vary, but these varied nucleotides were not located in the functional seed sequence). Homologous miRNAs might have been generated by homologous miRNA genes or by multiple pre-miRNAs. A total of 93 predicted circRNAs and 1,399 miRNAs might interact. Of these, 51% of them were oppositely expressed: one was upregulated as the other was downregulated. Due to their many-to-many interactions, many RNAs were repeatedly estimated. We found that 511 miRNAs and 59 circRNAs were downregulated. The miRNA-circRNA network based on potential relationships indicated that some circRNAs and miRNAs had a star-like distribution, suggesting potential crucial network positions ( Figure 3).
We also found that target gene of miRNA may generate circRNA as ceRNA of the miRNA. For example, circ-GMPS   and miR-217 might interact, and GMPS was also a potential target of miR-217 according to TargetScan [50].

isomiR-circRNA Interaction at the isomiR Levels.
Based on the miRNA-circRNA pairs identified, we further investigated potential isomiR-circRNA interaction at the isomiR levels. For the dominantly expressed miRNA gene, both miR-#-5p and miR-#-3p might generate multiple isomiRs with various expressions (Figure 4(a)). Due to the expression characteristics of isomiRs, most isomiRs sequences are similar to the canonical miRNA sequence (Figure 4(b)), despite differences in expression and involvement with additional nontemplate nucleotides. Based on the potential binding characteristics of miRNA and circRNA, all of the relevant isomiRs from a given miRNA locus can also bind to the circRNA (Figure 4(c)). Thus, at the isomiR level, it is possible to detect many isomiR:circRNA pairs based on close sequence relationships. In order to survey these interactions, multiple isomiRs could first be clustered based on miRNA loci or different miRNA loci (some homologous miRNA loci generate similar isomiR sequences) (Figure 1(b)). Clustered specific miRNA hsa-miR-204-5p CDK14_hsa_circ_001834 sequences with enriched expressions and sequences could be used to predict miRNA-circRNA interactions. Thus, many types of small RNAs may interact with a specific circRNA, and these relationships are complex.

Discussion
Bioinformatic analysis indicates that 276 miRNAs bind to circRNAs according to starBase [46], and 2,603 mature miR-NAs have been reported and validated according to miRBase (version 21.0) [42]. Many miRNA gene families are involved in these interactions; this may be because homologous miRNAs always have higher similarity. Only 10.60% of all identified miRNAs potentially interact with circRNAs. Most miRNAs are not involved in miRNA-circRNA interactions, while miRNA-mRNA interactions are quite popular. miRNA-circRNA interactions are not specific, and many small RNAs may interact with different molecules. These interactions may be important to regulatory processes, particularly those regulated by the small RNA. Some circRNAs that have been identified as miRNA sponges binding miRNA and reducing miRNA activity include ciRS-7/CDR1 [30,34], circ-Sry [34], circ-ITCH [51,52], circ-Foxo3 [51], circ-HRCR [53], circ-HIPK3 [54], and hsa circ 001569 [33]. The ability of circRNA to act as miRNA sponges may also be important for the regulatory function of small RNAs. Herein, we aim to understand the two noncoding RNAs expression profiles and their interactions in male patients with esophagus cancer, focusing on the interactions at the isomiR level.
Based on sequence data, we collected several miRNA-circRNA pairs. Based on widespread miRNA gene families and gene clusters, many miRNA sequences are highly similar any may thus have similar expression patterns and functional relationships. Even though few miRNA-circRNA interactions were predicted when compared to the large numbers of RNA species, the details of the corresponding relationships vary substantially. However, similar to miRNA-mRNA and miRNA-miRNA interactions, the actual functional binding of the two ncRNAs in vitro remains unclear. Indeed, competitive relationships may also exist between multiple isomiRs. Of these isomiRs, many are repeatedly estimated due to high sequence similarity. miRNA-circRNA pairs are not as common as miRNA-mRNA pairs, but the role of potential miRNA sponges should be not ignored. For example, ciRS-7 abrogates miR-7 [55], circRNA 010567 suppress miR-141 [56], and circRNA 100290 may act as a sponge for the miR-29 gene family [57].
The miRNA-circRNA network that we constructed demonstrates the complex interaction patterns between the two ncRNAs ( Figure 2). This network would increase in complexity complex if we included the detailed isomiRs. A given miRNA locus may generate multiple isomiRs with difference sequence and expression patterns. These isomiRs mainly differ at the 3 ends. Thus, these isomiRs may potentially interact with circRNAs as well as their canonical miRNA sequence and isomiRs with high sequence similarity and similar expression relationships that are similar to homologous miRNAs [29]. Therefore, at the isomiR level, more complex interactions with circRNAs are possible, and miRNA sponges may absorb many types of small RNAs with similar sequences. The interesting phenomenon seems to enlarge interaction profiles as well as current pattern: a circRNA may be miRNA sponge for homologous miRNAs in a miRNA gene family, while multiple isomiRs can be considered a special larger gene family with high sequence similarity. At the isomiR level, the isomiR-circRNA interaction may be quite common and the network may be quite complex. Indeed, except for the interaction between isomiRs and circRNAs, interactions between miRNAs or isomiRs are also quite common, such as coordinated or restricted interactions, and these horizontal and vertical interactions increase RNA cross-talk. Moreover, based on the deregulated miRNAs and circRNAs we identified, the paired circRNAs may be derived from target gene of miRNA, although this seems rare compared to the number of miRNA-cirRNA interactions. CircRNA may suppress miRNA to increase the expression level of the miRNA target [34], simultaneously generating ceRNAs to control miRNAs to ensure gene expression.
In conclusion, based on the small RNA and circRNA expression profiles in male patients with esophagus cancer, some deregulated RNA species may show potential functional relationship, and their expression patterns are always opposed. Although fewer miRNAs are predicted to interact with circRNAs, the interaction between two ncRNAs is more common if detailed analysis is performed at the isomiR level. Similar to homologous miRNAs, isomiRs increase analysis and cross-talk complexity.
of Jiangsu (no. BK20171443), and Top-Notch Academic Programs Project of Jiangsu Higher Education Institutions (TAPP: PPZY2015A067). The work was also supported by the Qing-lan Project of Jiangsu Province and the Excellent Young Teacher Project of Nanjing Medical University, Changzhou High-Level Medical Talents Training Project (no. 2016CZBJ043), and the Applied Basic Research Program of Changzhou (CJ20159030) and Sponsored by NUPTSF (nos. NY215068 and NY217100) and the Priority Academic Program Development of Jiangsu Higher Education Institution (PAPD). Figure S1: clustering analysis showed expression distribution of circRNAs across different samples. Figure S2: Venn's distributions of expressed circRNAs between paired normal and tumor samples. Figure S3: clustering analysis and volcano distribution of miRNAs. Table S1: patients characteristics of 5 male patients with esophageal squamous cancer.