Comprehensive analysis of circRNA expression pattern and circRNA-miRNA-mRNA network in the pathogenesis of atherosclerosis in rabbits

Successful established of atherosclerotic rabbit model

After one week of endothelial injury and high-fat diet, atherosclerotic plaques appeared. The atherosclerotic plaques of the right common carotid arteries were easily detected by the 8th week on the 2D-ultrasound images (Figure 1A-C). However, the carotid arteries of rabbits in the control group had no significant changes and the intima remained smooth (Figure 1D-F). Besides, hematoxylin and eosin stained vessels adjacent to the RNA samples revealed plaques of different severities in the carotid arteries of the case group (Figure 2A-C), while there were no obvious abnormalities in the control group (Figure 2D-F).

Identification of circRNAs in rabbit carotid arteries

In the present study, we used CIRCexplorer2 to systematically identify and annotate circRNAs. In total, we identified 9,418 circRNAs from 2 groups, including 9,336 exonic EcircRNAs (99.55%) and 82 intronic IcircRNAs (0.45%) with isoforms that were not previously reported. In addition, the circRNA transcripts were broadly distributed in all chromosomes (Figure 3A). Approximately 8.65%, 8.08% and 8.03% came from chr13, chr1 and chr2, respectively, whereas the percentages of circRNA from any other chromosome was less than 7%. As shown in Figure 3B, 6,321 and 7,976 circRNAs were detected in the AS group and control group, respectively, and 4,879 circRNAs were detected in both groups. We used transcripts per million reads (TPM) to estimate the expression level of the circRNA transcripts. Almost all circRNA transcripts were expressed at low levels (Figure 3C). Sequence length analysis revealed that circRNA transcripts were mostly 200-400 bp in length (Figure 4A). In addition, we found the number of exons in circRNA transcripts was mainly range from 2 to 6 (Figure 4B).

Identification of differentially expressed RNAs in atherosclerosis

To investigate the possible biological functions of circRNAs in AS, we analyzed miRNA, circRNA and mRNA expression profile data from AS samples and normal controls. RNAs with a fold-change >2 and p-value ≤ 0.05 were identified as significantly differentially expressed. RNA-seq analysis identified 813 mRNAs aberrantly expressed in AS, including 333 upregulated and 480 downregulated mRNAs (Supplementary Table S1). Meanwhile, 48 and 111 circRNAs were upregulated and downregulated, respectively, in AS samples compared to in controls (Supplementary Table S2). In addition, we found 55 upregulated miRNAs and 58 downregulated miRNAs in AS samples compared to controls (Supplementary Table S3). The top 10 dysregulated mRNAs, circRNAs and miRNAs based on p-value were summarized in Table 1, 2 and 3, respectively. Hierarchical clustering showed that DEmiRNA, DEcircRNA and DEmRNA expression patterns among samples were distinguishable (Figure 5). The data suggested that the expressions of miRNAs, circRNAs and mRNAs in AS are different from those in control samples.

Construction of the DEcircRNA-DEmiRNA-DEmRNA triple network

Increasing evidence indicates that circRNAs can sequester relevant miRNAs through MREs to post-transcriptionally regulate gene expression. To investigate ceRNA regulation in atherosclerosis and identify AS-related circRNAs, we identified putative DEmiRNA-DEcircRNA and DEmiRNA-DEmRNA interactions using the miRanda software. In total, we obtained 345 DEmiRNA-DEcircRNA and 2,874 DEmiRNA-DEmRNA interaction pairs. We next extracted miRNAs that paired with both circRNAs and mRNAs to construct the DEcircRNA-DEmiRNA-DEmRNA triple network, which included 81 miRNAs, 115 circRNAs, 399 mRNAs and 3,007 interaction pairs (Figure 6). The degree distribution of the network nodes was investigated. We observed a power-law distribution with a slope of -1.238 and an R-squared value of 0.808. This suggested that the network displayed scale-free characteristics typical of a biological network rather than a random network (Supplementary Figure S1).

To investigate the potential functional implication of the circRNAs, we performed functional enrichment analysis for all genes in the DEcircRNA-DEmiRNA-DEmRNA triple network based on gene ontology biological process terms. These genes were enriched for 41 GO biological terms (GOTERM-BP-FAT) (FDR cutoff =0.01), and Figure 7A shows the top 20 terms. The top ten enriched GO-BP terms, according to the FDR values, included cell adhesion, immune response, cell-cell adhesion, single organism cell adhesion, single organismal cell-cell adhesion, cell activation, regulation of cell activation, positive regulation of immune system process, leukocyte cell-cell adhesion and regulation of T cell activation. We also visualized and clustered the enriched GO-BP terms using the EnrichmentMap plug-in in Cytoscape (Figure 7B). We found that the enriched GO-BP terms were mainly clustered into cell adhesion and cell activation categories. Furthermore, these enriched GO terms have been implicated in many physiological and pathophysiological activities of atherosclerosis. Therefore, the results indicated that circRNAs act as ceRNAs to regulate mRNAs by competing for shared miRNAs in the pathogenesis of atherosclerosis.

Based on the above results, we extracted the genes from the cell adhesion and cell activation categories and their corresponding miRNAs and circRNAs from the DEcircRNA-DEmiRNA-DEmRNA triple network to form sub-network, respectively. In the cell adhesion sub-network, 112 circRNAs and 60 mRNAs competed for 70 miRNAs (Supplementary Figure S2A). Meanwhile, 111 circRNAs and 42 mRNAs competed for 69 miRNAs in the cell activation sub-network (Supplementary Figure S2B). The sub-network indicated that one mRNA might compete with multiple circRNAs. In particular, many circRNAs were predicted to participate in 2 sub-networks at the same time. Therefore, it is reasonable to predict that these circRNAs might contribute to the molecular regulation of atherosclerosis.

Construction of dysregulated circRNA-related ceRNA network for atherosclerosis

To identify miRNA-mediated DEcircRNA-DEmRNA interactions, a hypergeometric test was used to evaluate the significance of the number of shared miRNAs between circRNA and mRNA pairs. We defined circRNA-mRNA competing pairs if the circRNA and mRNA shared the same miRNA and had a p-value <0.05. The results of the dysregulated circRNA-related ceRNA network (DCCN) was comprised of 1,452 edges between 365 mRNAs and 112 circRNAs (Supplementary Table S4, Figure 8A). The degree of distribution of the nodes in the DCCN approximately followed the power-law distribution with a slope of -1.429 and an R-squared value of 0.807 (Figure 8B). The average degree of the circRNAs was 12.96 and was significantly higher than the average degree of mRNAs (3.98) (p = 2.2e-12, Figure 8C). This indicates that circRNAs have more interactions with other nodes than mRNAs in the DCCN. In addition, we found that circRNA nodes had more betweenness centrality compared to mRNA nodes (p=2.2e-12, Figure 8D), which revealed that circRNAs have more communication functions in the DCCN. These results suggested that, although the circRNAs are non-coding RNA, they exhibit more specific topological properties than mRNAs in the DCCN. These results also revealed that a part of circRNA nodes linked many mRNA nodes in the DCCN and acted as hubs. In the network, the maximum degree node was ocu-cirR-novel-18038 (degree=60), which is located in the FN1 gene and competes with many atherosclerosis-related mRNAs.

Topological analysis of the dysregulated circRNA-related ceRNA network for atherosclerosis

We considered three topological properties of the ceRNA network: the node degree, betweenness centrality (BC) and closeness centrality (CC). In general, a higher degree indicated that the node was a hub that participated in more ceRNA interactions. A higher BC value implied that the node acted as a bridge connecting different network modules. A higher CC value suggested that the node tended to be at the center of the network. Then, we ranked the topological features of all nodes and found that seven circRNAs (ocu-cirR-novel-18038, ocu-cirR-novel-18298, ocu-cirR-novel-15993, ocu-cirR-novel-17934, ocu-cirR-novel-17879, ocu-cirR-novel-18036 and ocu-cirR-novel-14389) were common nodes that were present in the top 10 of each index described above (Table 4).

By mapping these circRNAs into a ceRNA network, we found that these circRNAs and their neighbors formed a key module (Figure 9A). Surprisingly, the module contained 176 dysregulated genes and some of them were atherosclerosis-related genes. To explore the potential functional implications of the seven circRNAs, we performed functional enrichment analysis of gene ontology (GO) for the mRNAs in the module (Figure 9B). Interestingly, these genes were also enriched to the immune response, cell adhesion and cell activation terms. Furthermore, we found that these genes were annotated to the platelet activation category. Previous studies have suggested that platelet activation plays a critical role in atherosclerosis. The circulating activated platelets and platelet-leukocyte aggregates promote the formation of atherosclerosis [22]. Platelets may contribute to atherogenesis through the delivery of proinflammatory factors to leukocytes and endothelial cells [23]. Therefore, these seven circRNAs have potential use as diagnostic factors and new potential therapeutic targets of atherosclerosis. However, further specific studies are required to confirm the potential values of these circRNAs.

Animals and tissue preparation

All rabbits (New Zealand white male rabbits) used in our study were obtained from the model animal center of the 2nd Affiliated Hospital of Harbin Medical University. Fan et al. reported the feasibility and validity of rabbit models to study human atherosclerosis [42]. We obeyed the principles of laboratory animal care and made efforts to minimize suffering. All procedures were conducted on the basis of the guidelines established by the National Institutes of Health. The Medical Ethics Committee on Animal Research of the 2nd Affiliated Hospital of Harbin Medical University approved of our study protocol (Ethics No. KY2016-090).

Six New Zealand white male rabbits (2.5-3.0 kg) were randomly divided into two groups. The rabbits in the AS group (case group, n=3) were fed a high-fat diet [10% lard (Shandong Shiyuantianjiaji Factory), 1% cholesterol (Shanghai Lanji technology), and 3% yolk powder (Shandong Shiyuantianjiaji Factory)]. The rabbits in the control group (n=3) were fed with a standard diet. The high-fat diet in the AS group was initiated immediately after balloon injury and continued for 8 weeks. Endothelial injury caused by a 2 F Fogarty balloon catheter (Boston Scientific, Temecula, California) was performed to accelerate the formation of atherosclerotic plaque. Rabbits were anaesthetized by intramuscularly injected xylazine (5 mg/kg), ketamine (35 mg/kg) and acepromazine (0.75 mg/kg). Anesthesia was maintained with isoflurane inhalation through a mask during the procedure. As previously reported [43], the balloon-induced endothelial injury of the right common carotid artery was carried out by a balloon catheter introduced through right external carotid artery cut-down. The balloon was gently inflated with 0.3 ml of saline and then retracted. This procedure was performed three times in each rabbit, and then the balloon catheter was removed. The incision was closed with sutures, and the rabbits were allowed to recover.

The progression of atherosclerotic plaques was monitored by 2D-ultrasound weekly. After several weeks, the atherosclerotic plaques of the right common carotid arteries were obvious. Following euthanasia, the rabbits were killed with overdose of pentobarbital sodium. Then, we intercepted 1 cm of the right common carotid artery of each rabbit in the case group with ultrasound guidance where the atherosclerotic plaques were most obvious. Similarly, we intercepted 1cm of the right common carotid artery of each rabbit in the control group. These blood vessels were immediately transferred to liquid nitrogen for preservation. The right common carotid artery adjacent to the incision of each rabbit was swiftly removed. Each specimen was fixed with 4% paraformaldehyde fixative and embedded in paraffin for subsequent hematoxylin and eosin (H&E) staining. Serial cross-sections with a thickness of 3 μm were stained with H&E and observed by light microscopy (Olympus, BX41, Tokyo, Japan).

Total RNA isolation, library construction and sequencing

Total RNA was isolated from the right common carotid arteries from six samples, including three atherosclerotic sample and three control samples, using TRIzol reagent (Life Technologies, Beijing, China) according to the manufacturer’s recommended protocol. The RNA purity was verified with a NanoPhotometer® spectrophotometer (IMPLEN, CA). The RNA integrity was assessed using the RNA Nano 6000 Assay Kit for the Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

Approximately 8 μg of RNA per sample was used in the protocol to deplete ribosomal RNA according to the instructions in the Epicentre Ribo -zero™ rRNA Removal Kit (Epicentre, USA). The rRNA-depleted RNAs were used as input material for preparing the libraries following the manufacturer’s recommendations in the NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina (NEB, USA). The products were purified (AMPure XP system), and the library quality was assessed on the Agilent Bioanalyzer 2100 system. Sample clustering was performed on a cBot Cluster Generation System using the TruSeq PE Cluster Kit v3-cBot-HS (Illumina, San Diego, USA) according to the manufacturer’s instructions. Paired-end sequencing with 125-bp reads was performed on the Illumina HiSeq 2500 platform.

A total of 3 μg of RNA per sample was used as input material for the small RNA library. Sequencing libraries were generated using the NEBNextR Multiplex Small RNA Library Prep Set for IlluminaR (NEB, USA) following the manufacturer’s recommendations. Barcodes were added to map the sequences to each sample. The library quality was assessed on the Agilent Bioanalyzer 2100 system using DNA High Sensitivity Chips. Clustering of the barcoded samples was performed on a cBot Cluster Generation System using the TruSeq SR Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer’s instructions. After cluster generation, the libraries were sequenced on an Illumina HiSeq 2500 platform and 50-bp single-end reads were generated.

RNA-seq data analysis

The quality of the raw sequencing reads was estimated by the FastQC program (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), which included measures of base quality distribution, average quality of reads and adaptor check. Before alignment, we used Trimmomatic [44] to remove adaptor sequences, trim the 5’ and 3’ end low-quality sequences, and discard the reads containing the N base. This yielded high-quality clean reads for use in all of the downstream analyses.

The reference genome (OryCun2.0) and gene annotation files were downloaded from the Ensembl genome browser (ftp://ftp.ensembl.org/pub/release-92/fasta/oryctolagus_cuniculus/). The clean mRNA reads were aligned to the rabbit genome using the HISAT2 program [45], and the read counts for expressed transcripts were calculated with StringTie software [46]. To identify candidate circular RNAs, we mapped the circRNA-seq reads to the rabbit genome using the read alignment tool Tophat2 [47]. The CIRCexplorer2 [48] program with default parameters was used for transcript assembly, quantification and novel circRNA prediction. Meanwhile, miRDeep2 [49] with default parameters was used to compare the expression of known miRNAs against precursor and mature miRNA sequences in miRBase v22 [50]. Transcripts with read counts of less than 2 in all six samples were removed. The differentially expressed RNAs (DEmRNA, DEmiRNA and DEcircRNA) (P-value ≤ 0.05 and fold change (FC) > 2) were then identified by the ‘edgeR’ R package with default parameters. Hierarchical clustering of the samples was performed to generate an overview of the expression profile characteristics based on the CPM values of differentially expressed RNAs with the ‘pheatmap’ R package.

Identification of miRNA-target interactions

The miRNA-mRNA and miRNA-circRNA interactions in our study were predicted with miRanda (version 3.3a) [51], which predicted microRNA targets based on sequence complementarity, conserved target sites and free energy of formation. The input 3’UTR sequences of mRNA from the OryCun2.0 assembly were retrieved from the Ensembl database. To improve the reliability of miRNA target prediction, we used a maximum binding free energy of -20 kcal/mol for the miRNA-target interaction predicted in miRanda.

Identification of AS-related mRNA-circRNA interactions

We next extracted the miRNAs that were predicted to interact with both circRNAs and mRNAs. These miRNAs and their targets were used to construct a DEcircRNA-DEmiRNA-DEmRNA triple network. To identify competing circRNA-mRNA interactions, we used a hypergeometric test which enabled the evaluation of the significance of the shared miRNAs between each circRNA and mRNA. We considered a p-value < 0.05 as statistically significant. The p-value was calculated as follows:

$p-value=1-{\sum }_{i=0}^{r-1}\frac{C_t^i{C}_{m-i}^{n-i}}{C_m^n}$

where m is the number of miRNAs in the genome, and r represents the number of miRNAs shared between mRNAs and circRNAs. The number of miRNAs interacting with mRNAs and circRNAs is represented by t and n, respectively.

Network visualization and topological analysis

We used Cytoscape software (version 3.4.0) [52] to construct and visualize the network in this study. Several topological properties such as the node degree, betweenness and closeness were analyzed using the built-in NetworkAnalyzer tool. The degree of a node is the number of edges that link to this node. Betweenness is a measure of the centrality of the node in a network, which is the number of shortest paths from each node to all others that pass through the node. The closeness of a node is the average length of the shortest path between the node and all other nodes in the network. The Wilcoxon rank-sum test was used to test for significant differences in topological properties. The statistical calculations were performed with R software version 3.3.1. A p-value of less than 0.05 was considered statistically significant.

Functional enrichment analysis

Functional enrichment analysis at the GO level was performed using the DAVID bioinformatics resource (https://david.ncifcrf.gov/summary.jsp, version 6.8) [53,54]. In our work, GO terms for “Biological Process” (GOTERM-BP-FAT) with a threshold of FDR < 0.05 and p-value < 0.01 were considered significant functional categories using the whole rabbit genome as the background. Functional categories were visualized and clustered using EnrichmentMap [55], a plugin in Cytoscape.