Circular RNA profile of infantile hemangioma by microarray analysis

Background Circular RNAs (circRNAs) are a recently identified class of noncoding RNAs that participate in several physiological processes. However, the expression of circRNAs in infantile hemangioma (IH) remains unknown. Methods The profile of circRNAs was assessed by microarray in four pairs of IH and adjacent skin tissues. The expression of circRNAs was validated by quantitative reverse transcription polymerase chain reaction (qRT-PCR). Furthermore, circRNA-microRNAs (miRNA)-mRNA networks were constructed using bioinformatics tools. Results 234 up- and 374 down- regulated circRNAs were identified in IH by microarray. Among them, the expression of two up-regulated circRNAs (hsa_circRNA_100933 and hsa_circRNA_100709) and one down-regulated circRNA (hsa_circRNA_104310) was confirmed by qRT-PCR. In addition, 3,019 miRNA response elements (MREs) of circRNAs were predicted, and two circRNA-miRNA-mRNA networks were constructed, including 100 and 94 target genes of hsa_circRNA_100933 and hsa_circRNA_104310, respectively. GO and pathway analysis showed that both networks participated in angiogenesis and vascular development-related biological processes. Conclusions This is the first study to reveal the profiling of circRNAs in IH and pave the way for further characterization of the role of circRNAs in the pathogenesis of IH.


Introduction
Infantile hemangioma (IH) is a common vascular tumor in infancy, especially in females and premature with an increasing incidence of 3-10% [1]. Most IHs possess a unique life cycle composed of proliferation and involution stage; the proliferation is usually rapidly in the first year, followed by spontaneous resolution until 3-7 year old [2]. Notably, the rapid growth during proliferation stage may lead to ulceration, scarring, disfiguring, or even life-threatening complications [3].
Angiogenesis and vasculogenesis have been recognized as underlying neovascularization in IH [4,5]. The histopathological feature of proliferative IH includes masses of endothelial cells and closely packed capillaries. A high level of circulating endothelial progenitor cells (EPCs) support the hypothesis that IHs originate from intrinsic EPC [6,7]. Moreover, the theory of placenta origin suggests that the impairment of placenta during gestation may contribute to endothelial embolization and proliferation [8]. Accumulating evidence suggests that the placenta antigen glucose transporter 1 (GLUT 1) expresses in endothelial cells of IH and is an immunochemical marker for IH [9]. The extrinsic factors, such as hypoxia and acidosis and intrinsic factors, such as hypoxia-inducible factor 1 alpha (HIF1-α), vascular endothelial growth factors (VEGFs), and other angiogenic factors have been suggested to contribute to the development of IH [10,11]. Despite multiple hypotheses, the pathogenesis of IH is largely unknown.
Based on RNA sequencing data, a majority of human transcripts are non-coding RNAs (ncRNAs), with < 2% of the genome encoding the proteins [12]. An increasing volume of evidence suggests that ncRNAs are vital regulators in the vascular development of various disorders. Our previous study comparing the profile of long noncoding RNAs (lncRNAs) in IH and adjacent normal tissues revealed that lncRNAs are involved in angiogenesis and vascular biology [13]. Moreover, microRNAs (miRNAs), known as the common group of small ncRNAs, have been proved to play a pivotal role in IH pathogenesis. A recent study showed that chromosome 19 miRNA cluster (C19MC) was overexpressed in both IH tissues and blood samples, identifying C19MC as a candidate biomarker of IH [14]. MiR-382 is overexpressed in IH endothelial cells and contributes to IH development [15]. The inhibition of miR-130a suppresses the cell growth and angiogenesis of IH [16]. Nonetheless, the mechanism underlying the regulation of miRNA in IH necessitates further studies. Circular RNAs (circRNAs) are a recently identified class of ncRNAs that participate in several physiological processes [17]. Although circRNAs were discovered in 1972, they were mistaken for derivatives of mis-splicing [18,19]. The biology and function of circRNAs have aroused attention until 2012 thanks to high throughput technology [20]. Importantly, cir-cRNAs have been demonstrated to serve as miRNA sponges for regulating gene expression [21,22]. Some circRNAs compete with mRNAs for the same miRNA binding sites to form an integral network of posttranscriptional regulation, termed as the competing endogenous RNAs (ceRNAs) network [23,24]. Recently, circRNAs profiling have been identified in several kinds of disorders, especially in tumors, such as cutaneous squamous cell carcinoma [25], epithelial ovarian carcinoma [26], basal cell carcinoma [27]. However, the circRNAs involved in IH are yet to be explored. Interestingly, a recent study identified a circRNA profile in human umbilical endothelial cells under hypoxic or normoxic conditions and indicated that circRNA cZNF292 was up-regulated by hypoxia and controlled the formation of tubers and spheroid sprouting of endothelial cells [28]. These characteristics indicate that circRNAs might play an important role in the etiology of IH.
Therefore, in this study we identified the profile of circRNAs in four pairs of proliferative IH and adjacent skin tissues by high-throughput microarray. Subsequently, a total of 8 circRNAs were examined using quantitative reverse transcription polymerase chain reaction (qRT-PCR) in 10 paired IH and adjacent skin tissues. The most significant up-(hsa_circRNA_ 100933) and down-regulated (hsa_circRNA_104310) circRNA were chosen for further study. Next, we predicted the interaction between circRNAs and miRNAs and constructed hsa_ circRNA_100933-miRNA-mRNA and hsa_circRNA_104310-miRNA-mRNA networks using the bioinformatics approach.

Ethics statement
This study was approved by the Institutional Ethics Review Board of Shandong Provincial Hospital Affiliated to Shandong University (No. 2016-206). Written informed consents were obtained from the guardians on behalf of the children enrolled in this study.

Patients and samples
The proliferative IH was diagnosed by two independent doctors according to clinical and pathological features. IH samples (Table 1)  RNA isolation and quality control RNA from each sample was isolated by TRIzol reagent (Invitrogen life technologies, Carlsbad, USA) based on the manufacturer's instructions. The purity and concentrations of total RNA samples were determined with NanoDrop ND-1000. RNA integrity was tested by denaturing agarose gel electrophoresis.

Microarray and data analysis
The complete microarray data were deposited in Gene Expression Omnibus (GEO) database (accession number GSE93522 [NCBI tracking system #18429863]). Sample labeling and array hybridization were performed based on the Arraystar's protocol. Total RNA from each sample was treated with RNase R (Epicenter, Madison, USA) to enrich circular RNAs and remove linear RNAs. Subsequently, fluorescent cRNA were obtained by amplification and transcription of circular RNAs through random priming (Arraystar Super RNA Labeling Kit). The labeled cRNAs (pmol Cy3/μg cRNA) were purified by RNeasy Mini Kit (Qiagen) and assessed with NanoDrop ND-1000. Then blocking agent (10× 5uL) and fragmentation buffer (25× 1uL) were added to make each labeled cRNA fragmented, followed by heating at 60˚C for 30 min. Labeled cRNA was diluted with hybridization buffer (2× 25uL). Hybridization solution (50 uL) was added into the gasket slide. Finally, microarray slide was incubated for 17 h at 65˚C in an Agilent Hybridization Oven. The hybridized arrays were washed, fixed, and scanned using the Agilent Scanner G2505C. R software limma package was used for data normalization and calculation. The foldchange between the IH and adjacent skin tissues was calculated. The statistical significance was calculated by t-test. circRNAs with fold-change !2and P-values 0.05 were regarded as a significant differential expression. qRT-PCR validation Total RNAs of 10 paired IH and adjacent skin tissues were reverse transcribed to synthesize cDNA by Gene Amp PCR System 9700 (Applied Biosystems, Carlsbad, USA). The selected cir-cRNAs should meet the criteria that fold change was more than 4 and raw intensity of each sample was greater than 100. In addition, circRNAs with miRNA response elements (MREs) related to tumor progression and angiogenesis reported in the literatures were selected preferentially. Finally, 8 circRNAs were selected, including 5 up-regulated circRNAs (circRNA_ 100709, circRNA_102116, circRNA_051239, cricRNA_100933 and circRNA_102039) and three down-regulated circRNAs (circRNA_023016, circRNA_001654 and circRNA_104310), for amplification by specific divergent primers (S1 Table). The housekeeping gene Human β-actin was selected as an internal reference. All cDNAs were assembled in the ViiA 7 realtime PCR System (Applied Biosystems) in triplicates each. The reaction was based on the following procedures: 95˚C, 10 min; 40 PCR cycles (95˚C, 10 s; 60˚C, 1 min). The single-peak appearance of the melt curve demonstrated the specificity of PCR primers. The relative level of circRNAs was computed by the 2 -ΔΔCt method.

Construction of circRNA-miRNA-mRNA networks
The miRNA response elements (MREs) of circRNAs was predicted using Arraystar's miRNA target prediction software based on TargetScan [29] and miRanda [30]. The interaction of circRNAs and miRNAs were displayed by Cytoscape 3.01. The target genes of circRNAs were obtained by the integration of miRWalk 2.0 (miRWalk, miRanda, PITA, Targetscan) [31] and mRNA profile (GSE78811) previously [13]. GO analysis was used to predict the potential functions of the target genes by the Database for Annotation, Visualization, and Integrated Discovery (DAVID) online tool [32,33], including the categories of biological process (BP), molecular function (MF) and cellular component (CC). KEGG Ontology-Based Annotation System (KOBAS) 2.0 was used to proceed the pathway analysis of target genes, which incorporates five pathway databases [KEGG Pathway, PID, BioCyc, Reactome, and Panther] [34].

Statistical analysis
All data were presented as the mean ±standard error of mean (SEM). Student's t-test (twotailed) was used for data analysis. GraphPad Prism 5 and Microsoft Office were utilized for other statistics. P < 0.05 was considered to be statistically significant.

Analysis of differentially expressed circRNAs
In total, 13617 circRNAs were analyzed by microarray in 4 pairs of IH and adjacent skin tissues (S2 Table). Box plot demonstrated the intensity of all the samples after normalization (Fig 1A). Differentially expressed circRNAs between IH and adjacent skin tissues were acquired based on fold change filtering (Fig 1B). Finally, 234 up-and 374 down-regulated circRNAs (S3 Table) were identified by fold change ! 2.0 and P-value < 0.05 (Fig 1C). Hierarchical clustering was performed to visualize the differential circRNAs ( Fig 1D). The top 25 up-and downregulated circRNAs were displayed in Tables 2 and 3. The classification schematic revealed that up-regulated circRNAs included 211 exonic, 12 intronic, and 11 intragenic, while the down-regulated circRNAs included 256 exonic, 65 intronic, 36 intragenic, 8 intergenic, and 9 antisense (Fig 2A). In addition, we summarized the localization of circRNAs on human chromosomes ( Fig 2B). The data showed that the expression of circRNA in IH is much different from that of adjacent skin tissues.

qRT-PCR verification of selected circRNAs
On the basis of microarray raw intensity and downstream microRNAs, we selected 8 circRNAs from the most significant circRNAs for further verification. The expression of selected cir-cRNAs was substantiated by qRT-PCR in 10 pairs of IH and adjacent skin tissues including the previous microarray group. As a result, six of them (circRNA_100709, circRNA_102116, circRNA_051239, circRNA_100933, circRNA_102039, and circRNA_104310) had the same change direction with microarray results (Fig 3). The other two circRNAs (circRNA_023016 and circRNA_001654) showed a similar expression level between IH and skin tissues. Finally, two up-regulated circRNAs (circRNA_100933 and circRNA_100709) and one down-regulated circRNA (circRNA_104310) were validated with statistical significance.

GO and pathway analysis of target genes
To gain further insight into the function of target genes, GO and pathway analysis was conducted using DAVID and KOBAS. The top five significant GO terms of each subgroup (BP, CC, MF) were displayed. Intriguingly, the hsa_circRNA_100933-targeted genes were primarily involved in the physiological processes including epithelial cell migration, epithelium migration, vasculature development, and tissue migration (Fig 6A). The top ten pathway terms were also displayed (Fig 6B), which encompassed several angiogenesis-related pathways, such as signaling by EGFR; EGFR interacts with phospholipase C-gamma and Rap1 signaling pathway. In addition, the significant GO terms of hsa_circRNA_104310 included growth, developmental growth, cellular lipid metabolic processes, and intracellular signal transduction ( Fig 7A). Moreover, the HIF signaling pathway was among the top ten pathways of hsa_-circRNA_104310 (Fig 7B) that played a vital role in the pathogenesis of IH.

Discussion
The present study, for the first time, examined the expression of circRNAs in IH and identified 234 up-and 374 down-regulated circRNAs. A majority of the deregulated circRNAs were exonic; the intronic type was more regular in down-regulated circRNAs than those who are up-regulated. In addition, the expression of hsa_circRNA_100933, hsa_circRNA_100709 and hsa_circRNA_104310 was confirmed by PCR. CircRNAs are a large group of stable ncRNAs, which are widely found in various cells and exhibit a critical regulatory role in gene expression Red represents up-regulation in IH, yellow represents downregulation in IH, blue represents miRNAs, green represents mRNAs. [35]. In 2013, Memczak et al. [36] successfully identified more than 1980 kinds of circRNAs using high-throughput technology. Since then, the function of circRNAs in the tumor formation was gradually recognized. The characteristics of circRNAs include high stability, wide variety, transcription from exons, rich MREs, sequence specificity, and conserved across species [37]. These features indicate that circRNAs play a critical role in both transcriptional and post-transcriptional levels, and may serve as a new diagnostic and therapeutic target of IH. The functional study of circRNAs, which designates them to act as miRNA sponges, is an emerging research area of physiological processes and tumor formation. The common circRNA, CDR1as/ciRS-7, which harbors 71 miR-7 binding sites [38], was shown as the inhibitor of miR-7 in embryonic zebrafish midbrain [36] and islet cells [39,40]. The present study predicted the interaction of circRNAs and miRNAs. Among them, miR-137 targeted by cir-cRNA_100933, has been reported to act as tumor suppressor in various tumors [41,42]. Previous study demonstrated that miR-216a-3p, one of the target miRNAs of circRNA_104310, inhibited gastric and colorectal cancer cell proliferation and migration [43,44]. Moreover, hsa-520d-3p and hsa-miR-520a-5p targeted by hsa_circRNA_001654 and hsa_miR-518f-3p targeted by hsa_circRNA_102116 were implicated as biomarkers in IH [14], suggesting the potency of circRNAs as biomarker candidates in IH.
As circRNA and mRNA can compete for the same miRNA binding sites according to ceRNA theory [45], we constructed two circRNA-miRNA-mRNA networks with an integration of miRWalk prediction and mRNA profile(GSE78811) to gain further insight into the function of circRNAs in IH. As a result, 100 target genes of hsa_circRNA_100933 and 94 target genes of hsa_circRNA_104310 were acquired. Moreover, GO and pathway analysis suggested a number of angiogenesis and vascular development-related physiological processes and pathways. The hsa_circRNA_100933-targeted genes were primarily involved in neovascularization-related biological processes including epithelial cell migration, epithelium migration, vasculature development, and tissue migration. EGFR signaling and Rap1 signaling pathways that included the top 10 pathway terms, played a major role in vasculogenesis and angiogenesis of several vascular diseases [46,47]. Among the top 10 pathways of hsa_circR-NA_104310-targeted genes, HIF-1 signaling pathway has been demonstrated as an essential pathogenic mechanism in IH, which could regulate angiogenesis via VEGF signaling [48]. The bioinformatics data indicated that circRNAs may be involved in molecular mechanisms underlying neovascularization and may play a vital role in the mechanism of IH. Nevertheless, the function of circRNAs and the circRNA-miRNA-mRNA network in IH requires further investigations.

Conclusion
This is the first study to prolife circRNAs expression in IH. Furthermore, the circRNA/ miRNA interaction was predicted, and circRNA-miRNA-mRNA networks were constructed using bioinformatics tools. The GO and pathway analysis identified many vasculogenesis-and angiogenesis-related pathways, which may provide new insights into the pathogenesis of IH. However, the function of circRNAs in IH requires further in vitro and in vivo studies.