Research Paper Volume 12, Issue 3 pp 2798—2813

Functional lncRNA-miRNA-mRNA networks in rabbit carotid atherosclerosis

Yingnan Wu1, , Feng Zhang2, , Rui Lu1, , Yanan Feng1, , Xiaoying Li1, , Shuang Zhang1, , Wenying Hou3, , Jiawei Tian1, , Xianchao Kong4, *, , Litao Sun5, *, ,

  • 1 Department of Ultrasound, The 2nd Affiliated Hospital of Harbin Medical University, Harbin 150081, Heilongjiang, China
  • 2 Department of Ultrasound, The First Affiliated Hospital of Xiamen University, Xiamen 361003, Fujian, China
  • 3 Department of Ultrasound, Xuanwu Hospital Capital University, Beijing 100053, China
  • 4 Department of Gynecology and Obstetrics, The 2nd Affiliated Hospital of Harbin Medical University, Harbin 150081, Heilongjiang, China
  • 5 Department of Ultrasound, Shenzhen University General Hospital, Shenzhen 518055, Guangdong, China
* Equal contribution

Received: April 9, 2019       Accepted: January 19, 2020       Published: February 11, 2020
How to Cite

Copyright © 2020 Wu et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.


Atherosclerosis is one of the most common clinical cardiovascular disorders. Accumulating evidence indicates that lncRNAs exert critical functions in atherosclerosis; however, their functional roles and regulatory mechanisms remain unclear. In this study, we induced atherosclerotic plaques in three rabbit carotid arteries through an atherogenic diet and balloon injury; three age-matched rabbits were fed normal chow and served as controls. We thoroughly investigated the RNA (mRNA, lncRNA and miRNA) expression profiles in atherosclerotic rabbit carotid models with deep RNA sequencing. We identified several significantly differentially expressed RNAs. The corresponding lncRNA-miRNA-mRNA network was constructed, and the significantly dysregulated network was selected. Furthermore, Gene Ontology and  Kyoto Encyclopedia of Genes and Genomes analyses indicated that the mRNAs in the network were involved in leukocyte activation, cell proliferation, cell adhesion molecules and cytokine-cytokine receptor interaction. After rigorous screening, we obtained a differentially expressed lncRNA-miRNA-mRNA interaction network associated with atherosclerosis. In the network, XLOC_054118 and XLOC_030217 upregulate the CHI3L1, SOAT, CTSB and CAPG genes by competitively binding to the miRNA ocu-miR-96-5p. XLOC_062719 and XLOC_063297 upregulate CTSS, CTSB and EDNRA genes by competitively binding to the miRNA ocu-miR-185-5p.


Atherosclerosis (AS), a chronic disorder affecting the blood vessel walls, is characterized by an imbalance between the inflammatory response and lipid metabolism [1]. AS is the common pathological basis of multiple cardiovascular and cerebrovascular disorders and has become the main cause of mortality and long-term morbidity worldwide [2]. Currently, AS medications primarily reduce plasma cholesterol concentrations and blood pressure, effectively reducing tissue damage caused by the disease [3]. However, AS-related mortality and morbidity remain high. Genetic factors represent a major determinant of AS risk [4]. Thus, there is a need to identify the molecular pathways responsible for AS development and to develop valuable therapeutic medicines and prognostic biomarkers for AS.

Long noncoding RNAs (lncRNAs) were recently recognized to be broadly transcribed in various eukaryotic genomes ranging from those of nematodes to those of humans [5]. LncRNAs comprise transcripts > 200 nucleotides in length and regulate many biological mechanisms [6], such as genomic imprinting and chromatin modification [7]. Additionally, accumulating evidence has indicated that lncRNAs have extensive and complex impacts on the development and progression of AS. For example, the expression levels of the lncRNAs Zfas1, SNHG6 and GAS5 distinctly increased in individuals with atherosclerotic plaques compared with the expression levels in normal controls [8]. Li et al. [9] found that lncRNAs regulated AS-related processes in endothelial cells, macrophages, smooth muscle cells, and lipid metabolism. However, the roles of lncRNAs and their functions in AS have mostly been unreported. Salmena et al. [10] suggested a complicated posttranscriptional gene expression regulatory network termed competing endogenous RNAs (ceRNAs), wherein circular RNAs (circRNAs), lncRNAs, and other noncoding RNAs act as molecular sponges to inhibit mRNA function by sharing at least one microRNA (miRNA) recognition element (MRE). LncRNAs with similar sequences to targeted miRNAs serve as ceRNAs to modulate the level of protein-coding genes and to modulate cell biology by sponging miRNAs [11]. LncRNA-related ceRNAs serve as novel posttranscriptional regulators and have been reported to have crucial roles in various disorders. The lncRNA XIST was markedly overexpressed in gastric carcinoma and acted as a ceRNA to meditate EZH2 expression by naturally sponging miRNA 101 [12]. The lncRNA HOST2 functioned as a molecular sponge of the miRNA let-7b to downregulate let-7b expression, consequently affecting epithelial ovarian carcinoma [13]. However, the mechanisms of lncRNA-related ceRNA in AS remain largely unknown. Therefore, it is necessary to study lncRNA-miRNA-mRNA competitive regulatory networks to comprehensively understand the impact of ceRNA crosstalk on AS.

In this study, miRNA, mRNA and lncRNA expression profiles were constructed in carotid atherosclerotic rabbit models with deep RNA sequencing (RNA-seq). Furthermore, we investigated differentially expressed (DE) profiles of miRNA, mRNA and lncRNA in AS and detected lncRNA-miRNA-mRNA networks in carotid atherosclerotic rabbit models. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed to explore the potential regulatory functions of lncRNA. Then, functional lncRNAs were detected in lncRNA-miRNA-mRNA networks with the highest function in AS. Our results provide valuable resources to develop therapeutic pharmaceutical targets and molecular diagnostic tools.


Establishment of atherosclerotic rabbit carotid artery models

By combining endothelial injury and an atherogenic diet, we thickened and roughened the intima and induced multiple irregular hypoechoic plaques in rabbit right common carotid arteries as indicated by two-dimensional (2D) ultrasound at the 12th week (Figure 1A1C). In contrast, the carotid artery intima of rabbits fed normal chow remained clear, smooth and continuous with no obvious changes (Figure 1D1F). Hematoxylin and eosin (H&E) was used to stain vascular tissue from sites with obvious plaques in the AS group, and tissues were imaged using 2D ultrasound. Thickened intimae, foam cell deposition, smooth muscle cell proliferation and varying degrees of plaques on rabbit carotid artery intimae were observed (Figure 2A2C). However, the intimae of the control group were smooth and intact, and no obvious foam cell deposition or smooth muscle cell proliferation was observed (Figure 2D2F).

2D-ultrasound images of rabbit carotid arteries at the 12th week. (A–C) 2D-ultrasound images reveal that obvious atherosclerotic plaques formed on rabbit carotid arterial intima as the arrows show. (D–F) 2D-ultrasound images demonstrate that carotid arteries intima of rabbits treated with normal chow still maintain smoothness as the arrows show.

Figure 1. 2D-ultrasound images of rabbit carotid arteries at the 12th week. (AC) 2D-ultrasound images reveal that obvious atherosclerotic plaques formed on rabbit carotid arterial intima as the arrows show. (DF) 2D-ultrasound images demonstrate that carotid arteries intima of rabbits treated with normal chow still maintain smoothness as the arrows show.

H&E staining of rabbit carotid arteries. (A–C) H&E-stained vasculature shows plaques of varying degrees in rabbit carotid arteries in the AS group as the arrows show, ×40. (D–F) H&E-stained arteries show that no evident abnormal changes are detected in rabbit carotid arteries of the case group, ×40.

Figure 2. H&E staining of rabbit carotid arteries. (A–C) H&E-stained vasculature shows plaques of varying degrees in rabbit carotid arteries in the AS group as the arrows show, ×40. (DF) H&E-stained arteries show that no evident abnormal changes are detected in rabbit carotid arteries of the case group, ×40.

Overview of lncRNA-seq, mRNA-seq and miRNA-seq in the rabbit carotid arteries

Overall, 594,089,506 raw reads were generated, comprising 298,318,754 reads for the carotid atherosclerotic animal models and 295,770,752 reads for the normal controls. From the raw data, reads with adapter, poly-N, and low-quality sequences were removed. In total, we obtained 572,726,694 clean reads in the two different groups, including 288,533,106 reads for the AS samples and 284,193,588 reads for the normal control samples. Furthermore, 33 annotated lncRNAs and 1,622 novel lncRNAs were identified and subsequently analyzed; these lncRNAs included 1,496 lincRNAs (90.4%) and 159 antisense lncRNAs (9.6%) that have not been reported in the past. Sequence length analyses demonstrated that the lncRNA transcripts mostly ranged in length from 200 to 1,000 bp. Furthermore, the length of the ORF in the lncRNA transcripts ranged from 38 to 545, and the number of exons ranged from 2 to 10. Overall, 20,588 protein-coding transcripts were identified. Detected mRNAs were subsequently utilized for further analysis.

As for miRNA-seq in the rabbit carotid arteries, a total of 48,096,306 raw reads were produced and contained 25,287,408 reads from the atherosclerotic rabbit carotid artery samples and 22,808,898 reads from the normal controls. Reads with poly-N sequences, 5’ adapter contamination, poly C, G, T, or A sequences, and low-quality reads were removed as were those without a 3’ adapter or insert tag; finally, 47,233,534 clean reads were obtained from the raw data, including 24,808,717 reads from the AS samples and 22,424,817 reads from the normal control samples. Overall, we found 12,998 known miRNAs and 398 novel miRNAs. These miRNAs were further analyzed.

DE analysis of AS samples versus normal controls

First, compared with the normal control samples, the AS samples had 73 significantly and DE lncRNA transcripts (q < 0.05, Supplementary Table 1): there were 23 upregulated transcripts and 50 downregulated transcripts in the AS rabbits compared with those in the normal rabbits. Additionally, 50 significantly and DE miRNAs were identified between the two groups (q < 0.05, Supplementary Table 2): 29 DE miRNAs were upregulated, whereas 21 were downregulated. Fragments per kilobase of exon per million fragments mapped (FPKM) were applied to evaluate the mRNA transcript expression level. Finally, we detected 1,099 DE mRNA transcripts (q < 0.05, Supplementary Table 3), including 512 upregulated and 587 downregulated transcripts. The top 10 significantly and DE lncRNAs, miRNAs, and mRNAs based on q-values are summarized in Tables 13.

Table 1. Top 10 significantly DE lncRNA transcripts between AS and control rabbits.


Table 2. Top 10 significantly DE miRNA transcripts between AS and control rabbits.


Table 3. Top 10 significantly DE mRNA transcripts between AS and control rabbits.


Construction of the lncRNA-miRNA-mRNA network

RNA transcripts can effectively interact with one another based on the ceRNA hypothesis. RNAs in ceRNA can compete for the same MREs to modulate one another. LncRNAs can absorb miRNAs by binding to miRNAs and subsequently exhibiting a miRNA sponge function. Prediction results from bioinformatic analyses indicated that 746 lncRNAs were targeted by 382 miRNAs; 701 lncRNAs acted as decoys for 369 miRNAs. Based on the lncRNA-miRNA and miRNA-mRNA interaction pairs, a lncRNA-miRNA-mRNA network was constructed. The network consisted of 15,241 lncRNA-miRNA relationship pairs and 42,732 miRNA-mRNA relationship pairs (Supplementary Table 4).

Functional annotation: GO and KEGG

RNA interactions exert various functions expressed in associated mRNA genes. GO and KEGG analyses were conducted on the mRNA genes of the lncRNA-miRNA-mRNA network. Briefly, 424 significant GO-BP terms were observed and enriched (q < 0.05, Supplementary Table 5) in the GO enrichment analysis. The top three terms included immune system process (GO: 0002376), cell activation (GO: 0001775), and biological adhesion (GO: 0022610). Some cognition-related terms were also visualized, such as cell adhesion (GO: 0007155), immune response (GO: 0006955), T cell activation (GO: 0042110), leukocyte activation (GO: 0045321), and cell proliferation (GO: 0008283) (Figure 3). Through the KEGG survey, 33 significant KEGG terms were clustered (p < 0.05, Supplementary Table 6). We observed several cognition-associated terms, including cell adhesion molecules (CAMs) (ocu04514), T cell receptor signaling pathway (ocu04660), cytokine-cytokine receptor interaction (ocu04060), ECM-receptor interaction (ocu04512), and natural killer cell-mediated cytotoxicity (ocu04650) (Figure 4). Overall, the lncRNA-miRNA-mRNA network is involved in the etiopathogenesis of AS from various aspects.

Top 30 terms of biological process in GO analysis of mRNA genes in the lncRNA-miRNA-mRNA network.

Figure 3. Top 30 terms of biological process in GO analysis of mRNA genes in the lncRNA-miRNA-mRNA network.

Top 20 KEGG terms in KEGG analysis of mRNA genes in the lncRNA-miRNA-mRNA network.

Figure 4. Top 20 KEGG terms in KEGG analysis of mRNA genes in the lncRNA-miRNA-mRNA network.

Analysis of the DE lncRNA-miRNA-mRNA network in AS

We further established three restrictions to determine the most likely link between the lncRNA-miRNA-mRNA network and AS. A selective analysis of the lncRNA-mRNA network was conducted, wherein the lncRNAs, miRNAs, and their target genes were significantly dysregulated between the AS rabbits and control group (q < 0.05). Next, the concentrations of the chosen lncRNAs, miRNAs, and their target genes in the carotid arteries of rabbits were determined. The selected dysregulated lncRNA-miRNA-mRNA networks included only one situation (Figure 5): lncRNA (upregulated in AS rabbits)-miRNA (downregulated in AS rabbits)-mRNA (upregulated in AS rabbits) (Supplementary Table 7). These target genes in the selected triple pairs should be related to AS. The filtered lncRNA-miRNA-mRNA network demonstrated high interactions and could be used to explore the biological mechanism of AS. Intriguingly, only four lncRNAs, XLOC_054118, XLOC_030217, XLOC_062719 and XLOC_063297 were identified in our selected data. CHI3L1, SOAT, CTSB and CAPG shared the common miRNA ocu-miR-96-5p with XLOC_054118 and XLOC_030217 as ceRNAs. CTSB, CTSS and EDNRA shared the common miRNA ocu-miR-185-3p with XLOC_062719 and XLOC_063297 as ceRNAs. The results are listed in Table 4.

The dysregulated lncRNA-miRNA-mRNA network in AS rabbits. The network was based on lncRNA-miRNA and miRNA-mRNA interactions. LncRNA (up in AS rabbits)-miRNA (down in AS rabbits)-mRNA (up in AS rabbits). The red circles represent mRNAs. The red triangles represent lncRNAs. The green circles represent miRNAs.

Figure 5. The dysregulated lncRNA-miRNA-mRNA network in AS rabbits. The network was based on lncRNA-miRNA and miRNA-mRNA interactions. LncRNA (up in AS rabbits)-miRNA (down in AS rabbits)-mRNA (up in AS rabbits). The red circles represent mRNAs. The red triangles represent lncRNAs. The green circles represent miRNAs.

Table 4. The lncRNA-associated ceRNA network that is most likely involved in AS pathogenesis.

lncRNA idstatusq valuemiRNA idstatusq valueTranscript idmRNA idstatusgene_nameq value
XLOC_063297Up0.001644ocu-miR- 185-5pdown0.035479ENSOCUT00000011082ENSOCUG00000011085upCTSS0.00164462
XLOC_054118Up0.001644ocu-miR- 96-5pdown0.022014ENSOCUT00000008207ENSOCUG00000008204upSOAT10.00164462


Atherosclerosis is a chronic inflammatory disorder involving various immune cells at lesion sites and is the main cause of cardiovascular disorders [14]. Efforts have been made to find novel therapeutic targets and biomarkers of AS; however, the focus has been limited to mRNA [15]. Recently, lncRNAs and circRNAs have received attention as novel diagnostic markers for a number of disorders, including carcinomas [16]. In our previous study, we performed a comprehensive analysis of the circRNA expression pattern and circRNA-miRNA-mRNA network in the pathogenesis of AS in rabbits [17]. However, the expression profiles and functions of lncRNAs in AS remain unknown. LncRNAs, miRNAs and mRNAs participate in large-scale ceRNA crosstalk through MREs. This crosstalk has exciting influences for posttranscriptional gene regulation in various physiological and pathophysiological processes.

In this study, we established a rabbit model of the progression and regression of AS through balloon injury followed by an atherogenic diet. Then, we applied a deep RNA-seq analysis to explore the changes in the expression levels of mRNAs, lncRNAs and miRNAs in the progression and regression of AS. First, we detected 73 significantly and DE lncRNA transcripts between the AS samples and the normal control samples with 23 upregulated trasnscripts and 50 downregulated transcripts. Furthermore, 50 significantly DE miRNAs were identified between the two groups: 29 DE miRNAs were upregulated, whereas 21 were downregulated. Finally, we detected 1,099 DE mRNA transcripts, including 512 upregulated transcripts and 587 downregulated transcripts.

Many of these genes are well-known AS-related genes, such as MARCO, which is highly expressed in AS patients [18]. Heme oxygenase-1 (HO-1) has been reported to function as an intrinsic protective factor against atherosclerotic lesion formation by inhibiting lipid peroxidation in rabbits [19]. TLR2, 4, and 8 mRNAs are overexpressed in rabbit aortas after a high cholesterol diet, and their expression is correlated with inflammatory and biochemical markers and the further progression of AS [20]. Pryshchep et al. [21] found that TLR8 expression was very low in the aorta but high in the carotid arteries. Our findings suggested that TLR2, 4, and 8 were also highly expressed in the carotid arteries of atherosclerotic rabbits. VEGFs and other growth factors and cytokines can accelerate neointimal formation and AS by influencing monocyte activation, adhesion, and migration and by enhancing vascular permeability [22]. CXCL8/G31P, an IL-8 analog, can inhibit the formation of atherosclerotic plaques in the coronary artery [23]. Depletion of miR-34a facilitated endothelial cell growth and blocked apoptosis in AS by upregulating Bcl-2 [24]. MiR-497 expression was negatively correlated with apelin protein expression in atherosclerotic lesions [25]. In our experiment, ocu-miR-34a-5p and ocu-miR-497-5p were also highly expressed in atherosclerotic rabbits.

LncRNAs with sequences similar to those of their target miRNAs serve as ceRNAs to modulate the level of protein-coding genes [10, 11]. According to this ceRNA hypothesis theory, to fully understand the impact of lncRNA-related ceRNA crosstalk on AS, a global miRNA-lncRNA and protein-coding mRNA triple network of AS rabbits and controls was constructed by calculating the Pearson correlation coefficient of miRNA versus lncRNA and that of miRNA versus mRNA. GO and KEGG analyses were performed for the genes in this network and indicated that many enriched terms were associated with AS, including leukocyte activation (GO: 0045321), cell proliferation (GO: 0008283), CAMs (ocu04514) and cytokine-cytokine receptor interaction (ocu04060). The development of atherosclerotic lesions is driven by chronic inflammatory and proliferative processes. Swapnil V. et al. [26] demonstrated that the selective activation of leukocyte GPR120/FFAR4 by n-3 PUFAs, which are atheroprotective in humans, resulted in decreased leukocyte inflammation and AS. The proliferation of vascular smooth muscle cells (VSMCs) in the media layer, stimulated by growth factors from different sources, is an essential step in the formation of plaques [27]. Vascular cell adhesion molecule-1 (VCAM-1), an inflammatory and atherosclerotic marker, is an adhesion molecule abundantly expressed by smooth muscle cells in atherosclerotic lesions and injured arteries. VCAM-1 facilitates monocyte infiltration into atherosclerotic vascular walls and increases VSMC migration, further exacerbating AS [28]. An et al. [29] found that IL-8, a proinflammatory cytokine, interacted with its receptor CXC chemokine receptor 2 on neutrophils leading to the formation of neutrophil extracellular traps to aggravate AS progression in vivo. Thus, we predicted that these lncRNAs may be correlated with AS by regulating gene expression.

We applied rigorous restrictions to select the greatest probable lncRNA-related-ceRNA network that participated in the origin, development and changes in AS. First, we selectively analyzed the lncRNA-mRNA network by examining the likely significantly dysregulated RNAs. Next, the concentrations of the selected lncRNAs, miRNAs, and their target genes were determined. Third, the target genes in the chosen triple pairs should be related to AS. Finally, six qualified triple pairs were selected (Figure 6). We discovered that XLOC_054118 and XLOC_030217 were ceRNAs of ocu-miR-96-5p targeting CHI3L1, SOAT1, CAPG and CTSB. Macrophages, which are the primary cells of innate immunity, have vital effects in each stage of AS [30]. CHI3L1 has been detected to to be secreted from differentiated macrophages in early-stage AS lesions [31]. SOAT1-/- regulated the suppression of inflammatory molecules, including TNF-α, IL-6 and IL-1b, and alleviated the development of atherosclerotic lesions, exhibiting a positive effect [32]. Association analysis revealed that the rs6886 polymorphism in CAPG was linked with carotid intima-media thickness, which is a validated marker of AS; this mechanism could be associated with differential macrophage migration ability and the inflammation process [33]. Levels of CTSB mRNA and protein in atherosclerotic lesions of apoE-deficient mice have been observed to increase; CTSB immunoreactivity levels were highest in areas next to the lumen and in macrophages [34]. Interestingly, we found that the above four genes are associated with macrophages and involved in the inflammatory response, which further affects AS. Pordzik J [35] reported that as a platelet miRNA, miR-96 could be exploited as a biomarker in inflammatory disorders, such as cardiovascular disorders. Therefore, to the best of our knowledge, this is the first report that implies the lncRNAs XLOC_054118 and XLOC_030217 function as molecular sponges of the miRNA ocu-miR-96-5p to upregulate macrophages and inflammation-associated gene expression consequently participating in AS. Further experiments should be conducted to acquire detailed information on the pathway of genes and function of this network.

The dysregulated lncRNA-miRNA-mRNA network most likely involved in AS pathogenesis. The circles represent mRNAs with blue indicating inflammation and yellow indicating angiogenesis. The red triangles represent lncRNAs. The green squares represent miRNAs.

Figure 6. The dysregulated lncRNA-miRNA-mRNA network most likely involved in AS pathogenesis. The circles represent mRNAs with blue indicating inflammation and yellow indicating angiogenesis. The red triangles represent lncRNAs. The green squares represent miRNAs.

We also found that XLOC_062719 and XLOC_063297 were ceRNAs of ocu-miR-185-3p targeting CTSB, CTSS and EDNRA. Inflammation and angiogenesis broadly interact by sharing several cellular and molecular mediators. A key molecule in both processes is CTSB [36]. In vivo, CTSB inhibition retards tumorigenesis and vascular growth [37]. In contrast, Im et al. [37] found CTSB to be antiangiogenic in vitro, limiting endothelial tube formation through the degradation of endogenous VEGF-A. EDNRA mediates vasoconstriction and cell proliferation [38]. Research on EDNRA gene polymorphism showed that rs5333 was significantly associated with mean carotid intima-media thickness in men [39]. Microvessels may increase atherosclerotic plaque growth and lesion instability [40]. Recent studies have suggested that CTSS may contribute to the degradation of arterial ECM, which provides a path for neovessel growth [41]. CTSS deficiency reduces microtubule formation in vitro and microvessel growth in vivo. Interestingly, we found that these three genes are associated with angiogenesis and proliferation, which further affects AS. K Shan et al. showed that miR-185-5p is a posttranscriptional regulator of inflammation. RNCR3 is significantly upregulated, which alleviates miR-185-5p repression, thereby upregulating the level of the miR-185-5p target gene, KLF2. This regulatory loop maintains an associated balance in endothelial function to resist proatherogenic stress [42]. Therefore, we hypothesized that XLOC_062719 and XLOC_063297 functioned as molecular sponges of the miRNA ocu-miR-185-5p to upregulate angiogenesis and proliferation-associated gene expression, consequently participating in AS. Additional research is required to fully reveal the pathway of targeted genes and the involvement of these three lncRNAs in regulating angiogenesis in AS.


The results of the present study have not been experimentally validated and require confirmation. Additionally, the sample sizes were also small (three control samples and three AS samples).

Future directions

As more functional lncRNAs are identified in AS, lncRNAs will be increasingly and progressively used as diagnostic and therapeutic markers. An RNAi therapeutic strategy for AS could be to interfere with the lncRNA expression program that underlies AS [15]. Its potential for treating AS is supported by the discovery of lncRNAs and their effect on gene expression profiles. The use of lncRNA strategies could open novel therapeutic avenues for AS. However, large groups of essential, important lncRNAs remain unreported due to technical or cognitive limitations. More research is required to identify the novel mechanisms and functions of lncRNAs.


This is the first study to comprehensively analyze dysregulated lncRNA-miRNA-mRNA networks in AS. We showed that functional lncRNAs are involved in regulating carotid atherosclerotic rabbit models. In the network, XLOC_054118 and XLOC_030217 upregulate the CHI3L1, SOAT, CTSB and CAPG genes by competitively binding to the miRNA ocu-miR-96-5p. Furthermore, XLOC_062719 and XLOC_063297 upregulate the CTSS, CTSB and EDNRA genes by competitively binding to the miRNA ocu-miR-185-5p. These four lncRNAs could act as prospective clinical markers associated with the development of AS. However, additional work is needed to uncover the underlying mechanisms of lncRNAs in AS.

Materials and Methods


Six New Zealand white adult male rabbits (2.5-3.5 kg) were purchased from the Model Animal Center of the Second Affiliated Hospital of Harbin Medical University. Fan et al. reported the feasibility and validity of rabbit models to study human AS [43]. All animals received humane care in compliance with the Guidelines for the Care and Use of Laboratory Animals. All procedures were conducted on the basis of the guidelines established by the National Institutes of Health. Our study protocol was approved by the Medical Ethics Committee on Animal Research of the 2nd Affiliated Hospital of Harbin Medical University (Ethics No. KY2016-090).

Carotid atherosclerotic animal models

In order to induce carotid atherosclerosis, three rabbits were fed an atherogenic diet including 10% lard (Shandong Shiyuantianjiaji Factory), 3% yolk powder (Shandong Shiyuantianjiaji Factory), and 1% cholesterol (Shanghai Lanji Technology) combined with basal feed as the AS group. After one week, the three rabbits were intramuscularly anesthetized with ketamine (35 mg/kg), xylazine (5 mg/kg) and acepromazine (0.75 mg/kg). Following the methods of a previous study [44], the right common carotid arteries (CCA) were injured with a 2F Fogarty balloon catheter (Boston Scientific, Temecula, California), which was gently advanced into the CCA through the external carotid artery. The balloon was gradually inflated at 2 atm and retracted. This procedure was repeated three times in each rabbit, and anesthesia was maintained by isoflurane inhalation. Then, the balloon catheter was removed, the incision was closed with a suture, and the rabbits were allowed to recover while on the atherogenic diet. The carotid intimae of these rabbits were monitored weekly by two-dimensional ultrasound until significant plaques appeared at the 12th week. Three rabbits were sacrificed by air embolism. From the right carotid artery with the most obvious plaque (based on ultrasonic detection), the carotid atherosclerotic plaque and intima were excised and immediately preserved in liquid nitrogen. Three age-matched rabbits were fed a basal diet served as the control group, and the same tissue samples of the carotid artery as in the AS group were excised at the 12th week. The obtained arteries were numbered, marked and fixed with a 4% paraformaldehyde fixative and embedded in paraffin for subsequent H&E staining. Continuous cross-sections with a thickness of 3 μm were stained with H&E and observed by light microscopy (Olympus, BX41, Tokyo, Japan).

RNA extraction and library preparation for sequencing

NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, Ipswich, MA, USA) was used to prepare lncRNA and mRNA libraries following the manufacturer’s recommendations. Paired-end reads of 150 bp were generated with an Illumina HiSeq X platform. After quality control, paired-end clean reads were mapped to the reference genome (OryCun2.0) by TopHat v2.0.9. Cufflinks ( was used to assemble and annotate the transcripts. According to the annotation of the genome sequence (OryCun2.0), known lncRNAs and mRNAs were identified. We used the remaining transcripts to screen for putative lncRNAs according to the following criteria: (1) exon number ≥ 2; (2) sequencing coverage ≥ 3; (3) length ≥ 200 bp; and (4) identification in more than one sample. The transcripts meeting the above criteria were further filtered by removing known non-lncRNA transcripts. Next, CPC (0.9-r2) [45] and Pfam-scan (v1.3) [46] were applied to evaluate the transcripts that passed the filters for coding potential. Only transcripts without coding potential were classified as novel lncRNAs.

NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, Ipswich, MA, USA) was used to prepare small RNA library preparation following the manufacturer’s recommendations. Single-end reads of 50 bp were generated with an Illumina HiSeq 2500 platform. After quality control, the clean reads were aligned to a reference sequence (OryCun2.0) using Bowtie [47]. Mapped reads were used to identify known miRNAs applying miRBase 20.0 [48]. MiREvo [49] and miRDeep2 [50] were utilized to predict novel miRNAs by investigating characteristic hairpin structures, minimum free energy and Dicer cleavage sites.

Identification and clustering analyses of DE lncRNAs, miRNAs and mRNAs

We used the DE Seq R package (1.8.3) to identify DE lncRNAs, microRNAs and mRNAs (Benjamini & Hochber method corrected p-value <0.05) for AS rabbits and controls. We applied the R package “pheatmap” to analyze the expression of lncRNAs, microRNAs and mRNAs with unsupervised hierarchical clustering. The expression of each RNA type was normalized for unsupervised hierarchical clustering. Expressions of lncRNA and mRNA were normalized using the following formula: FPKm =lg10 FPKM + 1. Expression of miRNA was normalized with the following formula: TPm =lg10 TPM + 1 (TPM, transcripts per kilobase million). Next, the degree of similarity between the expression profiles of samples was measured using the Euclidean distance (using the R package, “complete”).

Prediction of lncRNA and miRNA target genes

Correlating expression levels between lncRNAs and mRNAs were used to assess potential trans roles of lncRNAs (acting on non-neighboring genes). Based on the expression correlation coefficient (Pearson correlation absolute value greater 0.8), we examined the trans role of lncRNAs in coding genes. To predict miRNA targets, we searched for targets in the 3′ UTR of gene models. For genes lacking a predicted 3′ UTR, regions 1,000 bp downstream of the stop codon were included. Miranda was conducted to perform the prediction (free energy <10 kcal/mol and score >140) [51].

GO and KEGG enrichment analyses

We utilized the GO database to analyze the lncRNA-miRNA-targeted genes. GO analysis was conducted with the GO-seq R package. Briefly, 424 GO-BP terms were significantly enriched (q < 0.05). As a database resource, KEGG can comprehend high-level utility functions of biological systems. KOBAS software was used to check the statistical enrichment of lncRNA target or dysregulated genes in the KEGG pathways. Thirty-three KEGG terms were significantly enriched (q < 0.05).

Construction of the lncRNA-miRNA-mRNA network

LncRNAs predicted to function as miRNA targets or decoys by Fan's methods were selected [52] to construct the lncRNA-miRNA-mRNA network. Next, we calculated the Pearson correlation coefficient value between a microRNA and its target mRNA and subsequently selected strongly correlated miRNA-mRNA pairs (with values greater than 0.8 or less than -0.8) to define the miRNA-mRNA relationships. To construct the lncRNA-miRNA-mRNA network, each DE RNA hub must be in either a miRNA-mRNA pair or a lncRNA-miRNA pair. The hubs in the lncRNA-miRNA-mRNA network consisted of microRNAs, lncRNAs acting as microRNA decoys, lncRNAs acting as microRNA targets, and mRNAs acting as microRNA targets. We used Cytoscape (version 3.4.0) [53] to visualize this network.

Conflicts of Interest

There are no conflicts of interest to declare.


This work was supported by the National Natural Science Foundation of China (No. 81671689) and the Natural Science Foundation of Heilongjiang Province (H2017021).


  • 1. Weber C, Noels H. Atherosclerosis: current pathogenesis and therapeutic options. Nat Med. 2011; 17:1410–22. [PubMed]
  • 2. Husain K, Hernandez W, Ansari RA, Ferder L. Inflammation, oxidative stress and renin angiotensin system in atherosclerosis. World J Biol Chem. 2015; 6:209–17. [PubMed]
  • 3. Steffens S, Veillard NR, Arnaud C, Pelli G, Burger F, Staub C, Karsak M, Zimmer A, Frossard JL, Mach F. Low dose oral cannabinoid therapy reduces progression of atherosclerosis in mice. Nature. 2005; 434:782–86. [PubMed]
  • 4. Kovacic S, Bakran M. Genetic susceptibility to atherosclerosis. Stroke Res Treat. 2012; 2012:362941. [PubMed]
  • 5. Guttman M, Amit I, Garber M, French C, Lin MF, Feldser D, Huarte M, Zuk O, Carey BW, Cassady JP, Cabili MN, Jaenisch R, Mikkelsen TS, et al. Chromatin signature reveals over a thousand highly conserved large non-coding RNAs in mammals. Nature. 2009; 458:223–27. [PubMed]
  • 6. Ponting CP, Oliver PL, Reik W. Evolution and functions of long noncoding RNAs. Cell. 2009; 136:629–41. [PubMed]
  • 7. Pandey RR, Mondal T, Mohammad F, Enroth S, Redrup L, Komorowski J, Nagano T, Mancini-Dinardo D, Kanduri C. Kcnq1ot1 antisense noncoding RNA mediates lineage-specific transcriptional silencing through chromatin-level regulation. Mol Cell. 2008; 32:232–46. [PubMed]
  • 8. Chen L, Yao H, Hui JY, Ding SH, Fan YL, Pan YH, Chen KH, Wan JQ, Jiang JY. Global transcriptomic study of atherosclerosis development in rats. Gene. 2016; 592:43–48. [PubMed]
  • 9. Li H, Zhu H, Ge J. Long Noncoding RNA: Recent Updates in Atherosclerosis. Int J Biol Sci. 2016; 12:898–910. [PubMed]
  • 10. Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language? Cell. 2011; 146:353–58. [PubMed]
  • 11. Chen DL, Ju HQ, Lu YX, Chen LZ, Zeng ZL, Zhang DS, Luo HY, Wang F, Qiu MZ, Wang DS, Xu DZ, Zhou ZW, Pelicano H, et al. Long non-coding RNA XIST regulates gastric cancer progression by acting as a molecular sponge of miR-101 to modulate EZH2 expression. J Exp Clin Cancer Res. 2016; 35:142. [PubMed]
  • 12. Fang J, Sun CC, Gong C. Long noncoding RNA XIST acts as an oncogene in non-small cell lung cancer by epigenetically repressing KLF2 expression. Biochem Biophys Res Commun. 2016; 478:811–17. [PubMed]
  • 13. Gao Y, Meng H, Liu S, Hu J, Zhang Y, Jiao T, Liu Y, Ou J, Wang D, Yao L, Liu S, Hui N. LncRNA-HOST2 regulates cell biological behaviors in epithelial ovarian cancer through a mechanism involving microRNA let-7b. Hum Mol Genet. 2015; 24:841–52. [PubMed]
  • 14. Lusis AJ. Atherosclerosis. Nature. 2000; 407:233–41. [PubMed]
  • 15. Liu Y, Zheng L, Wang Q, Hu YW. Emerging roles and mechanisms of long noncoding RNAs in atherosclerosis. Int J Cardiol. 2017; 228:570–82. [PubMed]
  • 16. Quan J, Pan X, Zhao L, Li Z, Dai K, Yan F, Liu S, Ma H, Lai Y. LncRNA as a diagnostic and prognostic biomarker in bladder cancer: a systematic review and meta-analysis. Onco Targets Ther. 2018; 11:6415–24. [PubMed]
  • 17. Zhang F, Zhang R, Zhang X, Wu Y, Li X, Zhang S, Hou W, Ding Y, Tian J, Sun L, Kong X. Comprehensive analysis of circRNA expression pattern and circRNA-miRNA-mRNA network in the pathogenesis of atherosclerosis in rabbits. Aging (Albany NY). 2018; 10:2266–83. [PubMed]
  • 18. Chu Y, Lao W, Jin G, Dai D, Chen L, Kang H. Evaluation of the relationship between CD36 and MARCO single-nucleotide polymorphisms and susceptibility to carotid atherosclerosis in a Chinese Han population. Gene. 2017; 633:66–70. [PubMed]
  • 19. Liu D, He Z, Wu L, Fang Y. Effects of induction/inhibition of endogenous heme oxygenase-1 on lipid metabolism, endothelial function, and atherosclerosis in rabbits on a high fat diet. J Pharmacol Sci. 2012; 118:14–24. [PubMed]
  • 20. Kapelouzou A, Giaglis S, Peroulis M, Katsimpoulas M, Moustardas P, Aravanis CV, Kostakis A, Karayannakos PE, Cokkinos DV. Overexpression of Toll-Like Receptors 2, 3, 4, and 8 Is Correlated to the Vascular Atherosclerotic Process in the Hyperlipidemic Rabbit Model: The Effect of Statin Treatment. J Vasc Res. 2017; 54:156–69. [PubMed]
  • 21. Edfeldt K, Swedenborg J, Hansson GK, Yan ZQ. Expression of toll-like receptors in human atherosclerotic lesions: a possible pathway for plaque activation. Circulation. 2002; 105:1158–61. [PubMed]
  • 22. Bhardwaj S, Roy H, Heikura T, Ylä-Herttuala S. VEGF-A, VEGF-D and VEGF-D(DeltaNDeltaC) induced intimal hyperplasia in carotid arteries. Eur J Clin Invest. 2005; 35:669–76. [PubMed]
  • 23. Qin Y, Mao W, Pan L, Sun Y, Fan F, Zhao Y, Cui Y, Wei X, Kohama K, Li F, Gao Y. Inhibitory effect of recombinant human CXCL8(3-72)K11R/G31P on atherosclerotic plaques in a mouse model of atherosclerosis. Immunopharmacol Immunotoxicol. 2019; 41:446–54. [PubMed]
  • 24. Su G, Sun G, Liu H, Shu L, Liang Z. Downregulation of miR-34a promotes endothelial cell growth and suppresses apoptosis in atherosclerosis by regulating Bcl-2. Heart Vessels. 2018; 33:1185–94. [PubMed]
  • 25. Cui J, Ren Z, Zou W, Jiang Y. miR-497 accelerates oxidized low-density lipoprotein-induced lipid accumulation in macrophages by repressing the expression of apelin. Cell Biol Int. 2017; 41:1012–19. [PubMed]
  • 26. Shewale SV, Brown AL, Bi X, Boudyguina E, Sawyer JK, Alexander-Miller MA, Parks JS. In vivo activation of leukocyte GPR120/FFAR4 by PUFAs has minimal impact on atherosclerosis in LDL receptor knockout mice. J Lipid Res. 2017; 58:236–46. [PubMed]
  • 27. Meneghini BC, Tavares ER, Guido MC, Tavoni TM, Stefani HA, Kalil-Filho R, Maranhão RC. Lipid core nanoparticles as vehicle for docetaxel reduces atherosclerotic lesion, inflammation, cell death and proliferation in an atherosclerosis rabbit model. Vascul Pharmacol. 2019; 115:46–54. [PubMed]
  • 28. Li HY, Leu YL, Wu YC, Wang SH. Melatonin Inhibits in Vitro Smooth Muscle Cell Inflammation and Proliferation and Atherosclerosis in Apolipoprotein E-Deficient Mice. J Agric Food Chem. 2019; 67:1889–901. [PubMed]
  • 29. An Z, Li J, Yu J, Wang X, Gao H, Zhang W, Wei Z, Zhang J, Zhang Y, Zhao J, Liang X. Neutrophil extracellular traps induced by IL-8 aggravate atherosclerosis via activation NF-κB signaling in macrophages. Cell Cycle. 2019; 18:2928–38. [PubMed]
  • 30. Binder CJ, Hartvigsen K, Chang MK, Miller M, Broide D, Palinski W, Curtiss LK, Corr M, Witztum JL. IL-5 links adaptive and natural immunity specific for epitopes of oxidized LDL and protects from atherosclerosis. J Clin Invest. 2004; 114:427–37. [PubMed]
  • 31. Rehli M, Niller HH, Ammon C, Langmann S, Schwarzfischer L, Andreesen R, Krause SW. Transcriptional regulation of CHI3L1, a marker gene for late stages of macrophage differentiation. J Biol Chem. 2003; 278:44058–67. [PubMed]
  • 32. Tavori H, Su YR, Yancey PG, Giunzioni I, Wilhelm AJ, Blakemore JL, Zabalawi M, Linton MF, Sorci-Thomas MG, Fazio S. Macrophage apoAI protects against dyslipidemia-induced dermatitis and atherosclerosis without affecting HDL. J Lipid Res. 2015; 56:635–43. [PubMed]
  • 33. Burillo E, Recalde D, Jarauta E, Fiddyment S, Garcia-Otin AL, Mateo-Gallego R, Cenarro A, Civeira F. Proteomic study of macrophages exposed to oxLDL identifies a CAPG polymorphism associated with carotid atherosclerosis. Atherosclerosis. 2009; 207:32–37. [PubMed]
  • 34. Chen J, Tung CH, Mahmood U, Ntziachristos V, Gyurko R, Fishman MC, Huang PL, Weissleder R. In vivo imaging of proteolytic activity in atherosclerosis. Circulation. 2002; 105:2766–71. [PubMed]
  • 35. Pordzik J, Pisarz K, De Rosa S, Jones AD, Eyileten C, Indolfi C, Malek L, Postula M. The Potential Role of Platelet-Related microRNAs in the Development of Cardiovascular Events in High-Risk Populations, Including Diabetic Patients: A Review. Front Endocrinol (Lausanne). 2018; 9:74. [PubMed]
  • 36. Nakao S, Zandi S, Sun D, Hafezi-Moghadam A. Cathepsin B-mediated CD18 shedding regulates leukocyte recruitment from angiogenic vessels. FASEB J. 2018; 32:143–54. [PubMed]
  • 37. Joyce JA, Baruch A, Chehade K, Meyer-Morse N, Giraudo E, Tsai FY, Greenbaum DC, Hager JH, Bogyo M, Hanahan D. Cathepsin cysteine proteases are effectors of invasive growth and angiogenesis during multistage tumorigenesis. Cancer Cell. 2004; 5:443–53. [PubMed]
  • 38. Zhang L, Sui R. Effect of SNP polymorphisms of EDN1, EDNRA, and EDNRB gene on ischemic stroke. Cell Biochem Biophys. 2014; 70:233–39. [PubMed]
  • 39. Yasuda H, Kamide K, Takiuchi S, Matayoshi T, Hanada H, Kada A, Yang J, Miwa Y, Yoshii M, Horio T, Yoshihara F, Nakamura S, Nakahama H, et al. Association of single nucleotide polymorphisms in endothelin family genes with the progression of atherosclerosis in patients with essential hypertension. J Hum Hypertens. 2007; 21:883–92. [PubMed]
  • 40. Moulton KS. Plaque angiogenesis and atherosclerosis. Curr Atheroscler Rep. 2001; 3:225–33. [PubMed]
  • 41. Libby P, Schönbeck U. Drilling for oxygen: angiogenesis involves proteolysis of the extracellular matrix. Circ Res. 2001; 89:195–97. [PubMed]
  • 42. Shan K, Jiang Q, Wang XQ, Wang YN, Yang H, Yao MD, Liu C, Li XM, Yao J, Liu B, Zhang YY, J Y, Yan B. Role of long non-coding RNA-RNCR3 in atherosclerosis-related vascular dysfunction. Cell Death Dis. 2016; 7:e2248. [PubMed]
  • 43. Fan J, Kitajima S, Watanabe T, Xu J, Zhang J, Liu E, Chen YE. Rabbit models for the study of human atherosclerosis: from pathophysiological mechanisms to translational medicine. Pharmacol Ther. 2015; 146:104–19. [PubMed]
  • 44. Meng L, Lv B, Zhang S, Yv B. In vivo optical coherence tomography of experimental thrombosis in a rabbit carotid model. Heart. 2008; 94:777–80. [PubMed]
  • 45. Kong L, Zhang Y, Ye ZQ, Liu XQ, Zhao SQ, Wei L, Gao G. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007 (suppl_2); 35:W345–9. [PubMed]
  • 46. Punta M, Coggill PC, Eberhardt RY, Mistry J, Tate J, Boursnell C, Pang N, Forslund K, Ceric G, Clements J, Heger A, Holm L, Sonnhammer EL, et al. The Pfam protein families database. Nucleic Acids Res. 2012; 40:D290–301. [PubMed]
  • 47. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009; 10:R25. [PubMed]
  • 48. Griffiths-Jones S. miRBase: microRNA sequences and annotation. Curr Protoc Bioinformatics. 2010; 29:12.9.1–12.9.10. [PubMed]
  • 49. Wen M, Shen Y, Shi S, Tang T. miREvo: an integrative microRNA evolutionary analysis platform for next-generation sequencing experiments. BMC Bioinformatics. 2012; 13:140. [PubMed]
  • 50. Friedländer MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012; 40:37–52. [PubMed]
  • 51. Enright AJ, John B, Gaul U, Tuschl T, Sander C, Marks DS. MicroRNA targets in Drosophila. Genome Biol. 2003; 5:R1. [PubMed]
  • 52. Fan C, Hao Z, Yan J, Li G. Genome-wide identification and functional analysis of lincRNAs acting as miRNA targets or decoys in maize. BMC Genomics. 2015; 16:793. [PubMed]
  • 53. Smoot ME, Ono K, Ruscheinski J, Wang PL, Ideker T. Cytoscape 2.8: new features for data integration and network visualization. Bioinformatics. 2011; 27:431–32. [PubMed]