Functional lncRNA-miRNA-mRNA networks in rabbit carotid atherosclerosis

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 1A–1C). In contrast, the carotid artery intima of rabbits fed normal chow remained clear, smooth and continuous with no obvious changes (Figure 1D–1F). 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 2A–2C). 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 2D–2F).

Figure 1. 2D-ultrasound images of rabbit carotid arteries at the 12^th 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 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. (D–F) 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 1–3.

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.

Figure 3. Top 30 terms of biological process in GO 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.

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.

Limitations

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.

Animals

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 2^nd 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 12^th 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 12^th 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 (http://cufflinks.cbcb.umd.edu/) 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 =lg₁₀ FPKM + 1. Expression of miRNA was normalized with the following formula: TPm =lg₁₀ 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.