Long non-coding RNA expression profiles in peripheral blood mononuclear cells of patients with coronary artery disease.

BACKGROUND: Coronary artery disease (CAD) is a leading cause of death and often presents as a complex systemic disease. The aim of the presents study was to determine the expression profiles of long non-coding RNAs (lncRNAs) in peripheral blood mononuclear cells (PBMC) of CAD patients and controls.
METHODS: The lncRNA expression profiling of PBMC obtained from 40 CAD patients and 10 non-obstructive coronary atherosclerosis (NOCA) subjects were analyzed using the Illumina Hiseq 4000 sequencer. Quantitative real-time polymerase chain reaction (qRT-PCR) was used to validate the differentially expressed lncRNAs. Gene ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of mRNA were conducted to predict biologic functions.
RESULTS: Our results indicated that lncRNAs were differentially expressed; 83 were upregulated lncRNAs and 114 were downregulated lncRNAs (P<0.05). A horizontal comparison of lncRNA expression indicated that the change in the expression profile of 48 lncRNAs was consistent with the degree of CAD. Six lncRNAs were validated using qRT-PCR, confirming the accuracy of the RNA sequencing analysis. GO analysis indicated that these dysregulated lncRNA transcripts were related to chromatin organization, cell and cell part, and protein heterodimerization activity. Pathway analysis indicated that these differentially expressed genes mainly included viral carcinogenesis, systemic lupus erythematosus and alcoholism.
CONCLUSIONS: Our preliminary findings indicate that lncRNAs could be used as potential biomarkers of subclinical cardiac abnormalities in the PBMC of CAD patients. However, further studies are needed to verify our findings and hypothesis. 2020 Journal of Thoracic Disease. All rights reserved.

Introduction

Cardiovascular disease is a leading cause of death, with more than 7 million people dying from coronary artery disease (CAD) annually (1,2). Although great progress in the field of percutaneous coronary intervention and stent implantation, the number of serious treatment-related adverse events and mortality continue to grow worldwide. It is important for physicians to identify individuals who are at high risk of CAD and to develop appropriate prevention strategies (3,4). There are a number of well-established CAD risk factors, including age, sex, smoking, alcohol use, and comorbidities, such as diabetes mellitus, dyslipidemia, hypertensive disease, and atherosclerotic vascular disease (5). Despite increased awareness of risk factors for CAD, the molecular pathogenesis of CAD is complex and not well understood (6-8). Presently, clinicians use several biomarkers for the differential diagnosis of CAD, mainly including Cardiac troponin I, Cardiac troponin T, creatine kinase-MB and heart fatty acid binding protein. The value of these biomarkers in clinical assessment of CAD risk and decision on CAD treatment is not satisfactory, which is partly due to the clinical symptoms of CAD varies greatly. Therefore, there is an urgent need to explore non-invasive biomarkers for the diagnosis and evaluation of CAD in the early stages (9-11).

Long non-coding RNAs (lncRNAs) are transcripts at least 200 nucleotides in length with little or no coding potential, but exhibit a variety of important biologic processes (12-15). In addition, it has been demonstrated that lncRNAs exist in human body fluid at measurable levels with high stability, and can accurately reflect the situation in vivo, indicating their potential use as biomarkers (16,17). These advances harbor the potential to better elucidate the complex pathological process of CAD (18-20). Several freely circulating lncRNAs have been related to distinct cardiovascular disorders, indicating that altered lncRNAs have a significant impact on the pathogenesis of CAD (19,21-23). CAD is usually caused by the abnormal activity of biological processes, such as cell motility, immune response, metabolism or inflammation response, and lncRNAs are directly or indirectly involved in it (7,9). Inflammation plays a crucial role in the development of CAD and peripheral blood mononuclear cells (PBMC) are one of the main type of immune cells. However, no studies to date have specifically focused on assessing the role of lncRNA on these aspects.

Therefore, the aim of the current study was to explore the expression of lncRNAs in PBMC of CAD patients by high throughput sequencing, and to identify lncRNAs as potential biomarkers for clinical and subclinical CAD. These findings could also provide a useful basis for exploring the role of lncRNAs in CAD progression. We present the following article in accordance with the MDAR reporting checklist (available at http://dx.doi.org/10.21037/jtd-20-3105).

Methods

Participants

A total of 40 CAD patients and 10 controls who underwent diagnostic coronary angiography at the Meizhou Peoples’ Hospital (Huangtang Hospital, Meizhou, Guangdong Province, China) between January 2016 and January 2018 were enrolled in the present study. CAD diagnosis was based on the coronary angiography, electrocardiogram and serological examination on the discharged medical records, according to the European Society of Cardiology/American College of Cardiology. The 50 participants were classified into 5 groups: stable angina (SA; n=10), unstable angina (UA; n=10), ST-segment elevation acute myocardial infarction (STEMI; n=10), non-STEMI (n=10), and non-obstructive coronary atherosclerosis (NOCA; n=10). NOCA on coronary angiography conformed was define as controls. Exclusion criteria were patients previously diagnosed with malignancies, stroke, valvular heart diseases, autoimmune diseases, or other malignant disease. To validate the lncRNA profiling data, another 120 patients and 30 controls were also included. Demographic information, medical history, and medication data of all participants were obtained by face-to-face or telephone interviews and recorded on the questionnaire. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The present study was approved by the Ethics Committee for Research of Meizhou Peoples’ Hospital (No. MPH-HEC 2015-A-23). Signed informed consent was obtained from each participant at the time of enrollment.

RNA isolation

A total of 10 mL ethylenediaminetetraacetic acid-anticoagulated venous blood samples were collected from all participants. The sample was centrifuged at 3,000 ×g for 15 min, and supernatants were again centrifuged at 10,000 ×g for 20 min. The plasma interlayer was carefully transferred, aliquoted, and stored at –80 °C until use. Total RNA was extracted using the miRNeasy Mini Kit (Qiagen, Hilden, Germany), following the manufacturer’s protocols. The total RNA concentration and protein contamination were assessed using the NanoPhotometer spectrophotometer (Implen, Los Angeles, CA, USA). RNA quality and integrity were estimated using the Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA).

LncRNA sequencing

Preparation for lncRNA sequencing library was based on the manufacturer’s standard protocols. In summary, ribosomal RNA was depleted from total RNA after enzymatic fragmentation by using the Ribo-Zero rRNA removal kit (EPICENTRE Biotechnologies, Madison, WI, USA). The resulting RNA pellet was then washed and recovered by alcohol precipitation. Libraries with high-strand specificity for sequencing were analyzed using the mRNA-ONLY Eukaryotic mRNA Isolation Kit (EPICENTRE Biotechnologies, USA). Samples for which the RNA integrity number was at least 8 were used to prepare a library of template molecules suitable for subsequent sequencing. The Illumina Hiseq 4000 sequencer (Illumina Inc., San Diego, CA, USA) was applied for sequencing, and a read length of 150 bp paired end was generated according to the standard protocol supplied by the manufacturer.

Mapping and identification of differently expressed genes

The sequencing raw reads were first converted into FASTQ format to filter out dirty raw reads, including sequence adaptors, reads containing ploy-N, and low-quality reads, prior to data analysis (24). The sequence reads were mapped onto the assembled genome, and transcript quantification and identification of differential gene expression was performed using the default parameters in Bowtie version 2.0.6 (http://bowtie-bio.sourceforge.net/index.shtml), TopHat version 2.0.9 (http://ccb.jhu.edu/software/tophat/index.shtml), and Cufflinks version 2.1.1 (http://cufflinks.cbcb.umd.edu/) (25). The read counts mapped to different isoforms of individual mRNA or lncRNA were pooled, and normalized and expressed as reads per kilobase per million mapped reads using the Cufflinks package. Differentially expressed genes were expressed as fold change with a false discovery rate (FDR) and corrected q-values. Only genes with an FDR <0.05 and |log2FC| >1 were identified as differentially expressed genes.

Functional enrichment analyses

Gene ontology (GO) enrichment analysis was performed to elucidate the biologic function of statistically significant genes, with GO annotations downloaded from GO (http://www.geneontology.org/) (26). Pathway analysis was used to analyze the potential functions in pathways where there was a significant enrichment of differentially expressed mRNAs based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg/pathway.html) (27). Fisher’s exact test was used in the GO and pathway analyses. Significance of data was set at P<0.05.

Quantitative real-time polymerase chain reaction (qRT-PCR)

LncRNA expression from the RNA sequencing data analysis was validated by qRT-PCR. A fixed volume of total RNA was subjected to reverse transcription using a PrimeScript RT Reagent Kit (Takara-Bio, Dalian, China) with a gDNA Eraser in accordance with the manufacturer’s protocol. The Premix Ex Taq II Reagent Kit (Takara-Bio, Dalian, China) and ABI 7500 Real-Time PCR thermocycler (Applied Biosystems, Carlsbad, CA, USA) were used for the qRT-PCR reactions. The primer sequences were listed in Table S1. Each reaction was conducted in triplicate. LncRNA expression was calculated using cycle threshold (Ct) values, and was internally normalized to glyceraldehyde 3-phosphate dehydrogenase using the 2−ΔΔCT equation method. All primer pairs are available upon request.

Statistical analysis

All statistical calculations were performed using SPSS version 19.0 (IBM, Armonk, NY, USA). Data for continuous variables are presented as means and standard deviations, and data for categorical variables are presented as numbers and percentages. The one-way analysis of variance comparison method was used to assess differences between different groups. P<0.05 indicated a statistically significant difference.

Results

Baseline clinical characteristics

A total of 160 CAD patients and 40 NOCA subjects were enrolled in the present. Characteristics of the study population are listed in . Statistically significant differences were not observed among the groups in terms of smoking, drinking, clopidogrel, and calcium blockers, as well as the level of systolic blood pressure (BP), diastolic BP, triglycerides, high-density lipoprotein cholesterol, and low-density lipoprotein cholesterol.

Table 1

Baseline clinical characteristics

Variable	NOCA	SA	UA	NSTEMI	STEMI	P value
Age (years)	59.9±9.3	63.7±8.9	66.1±11.4	68.8±11.1	67.1±9.9	0.002
Sex (male)	22 (55%)	31 (77.5%)	32 (80.0%)	35 (87.5%)	33 (82.5%)	0.007
Smoking	9 (22.5%)	16 (40%)	20 (50%)	16 (40%)	15 (37.5%)	0.155
Drinking	1 (2.5%)	2 (5%)	1 (2.5%)	0 (0%)	2 (5%)	0.662
Systolic BP (mmHg)	122.2±19.6	134.6±21.7	134.8±21.3	132.8±25.2	131.5±22.2	0.823
Diastolic BP (mmHg)	76.0±11.1	80.5±12.4	77.5±12.1	80.7±13.2	78.2±16.3	0.107
Hypertension	13 (32.5%)	28 (70%)	23 (57.5%)	23 (57.5%)	15 (37.5%)	0.004
Diabetes	2 (5%)	7 (17.5%)	14 (35%)	18 (45%)	16 (40%)	<0.001
Hyperlipidemia	6 (15%)	14 (35%)	10 (25%)	19 (47.5%)	14 (35%)	0.026
Aspirin, n (%)	8 (20.0%)	30 (75.0%)	31 (77.5%)	33 (82.5%)	34 (85.0%)	<0.001
Clopidogrel, n (%)	13 (32.5%)	22 (55.0%)	26 (65.0%)	21 (52.5%)	18 (45.0%)	0.052
ACEI/ARB, n (%)	19 (47.5%)	35 (87.5%)	34 (85.0%)	33 (82.5%)	30 (75.0%)	<0.001
β-blocker, n (%)	21 (52.5%)	31 (77.5%)	30 (75.0%)	35 (87.5%)	34 (85.0%)	0.002
PPI, n (%)	18 (45.0%)	17 (42.5%)	22 (55.0%)	28 (70.0%)	31 (77.5%)	0.004
Calcium blockers, n (%)	6 (15.0%)	11 (27.5%)	5 (12.5%)	5 (12.5%)	4 (10.0%)	0.277
TC, mmol/L	4.58±1.26	4.80±1.25	4.93±1.60	5.07±1.40	4.96±1.26	0.808
TG, mmol/L	1.27±0.39	2.02±1.87	1.60±1.17	1.99±1.20	1.74±1.29	<0.001
HDL-C, mmol/L	1.30±0.30	1.21±0.31	1.23±0.33	1.16±0.39	1.24±0.27	0.903
LDL-C, mmol/L	2.58±0.80	2.58±0.89	2.69±1.16	2.97±0.94	2.87±0.84	0.609

ACEI, angiotensin-converting enzyme inhibitor; ARB, angiotensin receptor blocker; BP, blood pressure; HDL-C, high-density lipoprotein cholesterol; LDL-C, low-density lipoprotein cholesterol; NOCA, non-obstructive coronary atherosclerosis; NSTEMI, non-ST-segment elevation acute myocardial infarction; PPI, proton-pump inhibitor; SA, stable angina; STEMI, ST-segment elevation acute myocardial infarction; TC, total cholesterol; TG, triglycerides; UA, unstable angina.

Quality assessment of lncRNA sequencing data

The cDNA sequencing libraries were constructed and sequenced through the Illumina HiSeq 4000 platform. Approximately 11.42 Gb sequence data were generated with the Q30, and averaged 94.14–95.32%, indicating high-quality data. Next, the clean reads were aligned to the human reference genome using TopHat version 2.0, where more than 93.96% of the total reads were useable reads that aligned uniquely to the reference genome with high uniformity. Detailed information on data quality and the mapping rate were summarized in Tables S2 and S3.

Four different algorithms including CPC (Coding Potential Calculator), CPAT (Coding Potential Assessment Tool), PFAM (Protein Family) and CNCI (Coding-Non-Coding Index) for calculating highly reliable potential non-coding transcripts were used. In the analysis, 2,925 candidate lncRNAs were detected, including 1,746 lincRNAs (59.7%) and 1,179 anti-sense lncRNAs (40.3%) (). Compared with lncRNAs, mRNA expression levels were generally higher (). In addition, significant differences in length, open reading frame length and exon number between mRNAs and lncRNAs were also observed ().

Figure 1

Identification of long non-coding (lncRNAs) and differentially expressed genes among coronary artery disease patients and control subjects. (A) Selection of lncRNAs using the Cufflinks package. (B) Screening of lncRNAs using CPC (Coding Potential Calculator), CPAT (Coding Potential Assessment Tool), PFAM (Protein Family) and CNCI (Coding-Non-Coding Index). Using the 4 tools to analyze the protein-coding potency of lncRNAs, 2925 lncRNAs were identified after the removal of putative protein-coding transcripts. (C) Distribution of the 3 types of lncRNAs.

Figure 2

Comparison of gene structure and expression analyses between long non-coding (lncRNAs) and mRNAs. (A) Expression level of mRNAs and lncRNAs. (B) Distribution of lncRNA and mRNA length. (C) Distribution of the length of open reading frames in lncRNAs and mRNAs. (D) Distribution of the number of exons in lncRNAs and mRNAs. FPKM, reads per kilobase per million mapped reads; ORF, open reading frame.

LncRNA expression profiles in PBMC

To specifically examine the entire repertoire of lncRNAs in PBMC, we determined the lncRNA expression profile from the PBMC of 40 CAD patients and 10 NOCA subjects. As shown in , the unsupervised hierarchical clustering of the lncRNA and mRNA expression profiles showed that the overall changes from normal to clinical CAD were separated (upregulated genes are in red and downregulated genes are in green). These findings indicate an underlying dynamic change in gene expression signatures according to stage at diagnosis during CAD development.

Figure 3

Hierarchical clustering of long non-coding (lncRNAs) and mRNAs in coronary artery disease patients and NOCA (non-obstructive coronary atherosclerosis) subjects. NOCA1–NOCA10, NOCA subjects; SA1–SA10, stable angina (SA) patients; UA1–UA10, unstable angina (UA) patients; NSTEMI1–NSTEMI10, non-ST-segment elevation acute myocardial infarction (NSTEMI) patients; STEMI1–STEMI10, ST-segment elevation acute myocardial infarction (STEMI) patients. Red and green indicate an increase and a decrease in expression level, respectively, across all samples. (A) lncRNA; (B) mRNA.

The Volcano plot of differential lncRNA expression in CAD subclinical patients and NOCA subjects are shown in . Of these, the expressions of 83 lncRNAs were up-regulated, whereas the expressions of 114 lncRNAs were downregulated (P<0.05). In addition, a horizontal comparison of lncRNA expression was further performed. As shown in , the change in the expression profile of 48 lncRNAs was found to be consistent with the degree of CAD, and could be used as potential biomarkers.

Figure 4

Volcano plot of differential long non-coding (lncRNA) expression. X-axis: log2 (fold change); Y-axis: −1×log10 (corrected q-value) for each probe. (A) stable angina patients (SA) vs. NOCA (non-obstructive coronary atherosclerosis) subjects; (B) unstable angina patients (UA) vs. NOCA subjects; (C) non-ST-segment elevation acute myocardial infarction (NSTEMI) patients vs. NOCA subjects; (D) ST-segment elevation acute myocardial infarction (STEMI) patients vs. NOCA subjects.

Figure 5

Hierarchical clustering of 48 long non-coding (lncRNAs) in coronary artery disease patients and non-obstructive coronary atherosclerosis (NOCA) subjects. Red and green indicate an increase and a decrease in expression level, respectively, across all samples. SA, stable angina; UA, unstable angina; NSTEMI, non-ST-segment elevation acute myocardial infarction; STEMI, ST-segment elevation acute myocardial infarction.

GO and KEGG pathway analyses

GO analysis was performed for gene and gene product enrichment based on biologic process, cellular component, and specific molecular function (). We found that dysregulated lncRNA transcripts are associated with chromatin organization (ontology: biologic process), cell and cell part (ontology: cellular component), and protein heterodimerization activity (ontology: molecular function). In addition, KEGG pathway analysis was further performed; the findings indicated that the pathways mainly included viral carcinogenesis, systemic lupus erythematosus and alcoholism ().

Figure 6

Gene ontology (GO) analysis of differentially expressed long non-coding (lncRNAs), which covers 3 domains: biologic process, cellular component, and molecular function. X-axis: GO terms of biologic process, cellular component, and molecular function. Green column indicates biologic process, red column indicates cellular component, and blue column indicates molecular function. Y-axis on the left: number of genes (lncRNAs).

Figure 7

Pathway analysis of differentially expressed long non-coding (lncRNAs). Pathway analysis is a functional analysis mapping genes to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and other pathway databases; the lower the P value, the more significant the pathway.

Evaluation of lncRNA expression via qRT-PCR analysis

To validate the expression levels of lncRNA identified through sequencing, 3 upregulated lncRNAs (RN7SKP203, SMIM25, and RP11-482G13.1) and 3 downregulated lncRNAs (RP11-763B22.4, XLOC_055653, and RN7SL471P) with only 1 transcript were randomly selected and evaluated. As shown in , there was consistency between the qRT-PCR and RNA sequencing analysis results in the expression of 6× lncRNAs detected in CAD patients compared with NOCA subjects.

Figure 8

Validation of RNA sequencing results using quantitative real-time polymerase chain reaction. Data are expressed as mean ± standard deviation (n=30). *P<0.05, **P<0.01, ***P<0.001 vs. non-obstructive coronary atherosclerosis (NOCA). Ns, no significance; SA, stable angina; UA, unstable angina; NSTEMI, non-ST-segment elevation acute myocardial infarction; STEMI, ST-segment elevation acute myocardial infarction.

Discussion

CAD is a leading cause of death, and often presents as a complex systemic disease. Despite significant advances in the prevention and treatment of CAD, including antiplatelet agents, anticoagulants, lipid-lowering therapy, and percutaneous coronary intervention, the prevalence of CAD is still increasing rapidly (2,28-30). The identification of novel biomarkers could improve the management of CAD patients (7,31). Several lncRNAs have been shown to be involved in the physiological and pathophysiological processes of CAD, and have the potential to be used as biomarkers (10,32). However, the majority of them remain largely unknown. In the present study, we performed a systematic sequencing analysis of differentially expressed lncRNAs in different clinical CAD subtypes. Our findings indicated that the expression levels of lncRNAs and mRNAs significantly changed in the PBMC between the 40 CAD patients and 10 NOCA subjects. Interestingly, the expression levels of some lncRNAs were altered depending on the stage of subclinical coronary atherosclerosis.

LncRNAs are non-coding transcripts with a length >200 nt and lack protein-coding capacity. They have been found to play an important role in various important biologic processes and diverse human diseases (6,33). Although the functions of lncRNAs involved in the physiological and pathological processes of CAD are not fully understood, the findings of several previously published studies indicate their beneficial effects on cardiovascular disease. LncRNA-p21 enhances p53 activity and promotes the progression of atherosclerosis, including vascular intimal hyperplasia, vascular smooth muscle cell proliferation, and apoptosis (34). Smooth muscle enriched lncRNA (SMILR) has been identified as a driver of vascular smooth muscle cell proliferation, which may be a novel therapeutic strategy for vascular pathologies (35). LncRNA antisense non-coding RNA in the INK4 locus (ANRIL) has been identified as a prototype of a cirRNA regulating ribosome biogenesis and might promote antiatherogenic cell function (36). The lncRNA metastasis-associated lung adenocarcinoma transcript 1 (MALAT1) has been found to be involved in regulating endothelial cell function and angiogenesis (37), whereas the smooth muscle and endothelial cell-enriched migration/differentiation-associated long non-coding RNA (SENCR) mediates a repressive effect on the migration of smooth muscle cells (38).

In the present study, we investigated the expression profiles of lncRNAs in PBMC in the CAD subtypes of a Chinese population. We identified 83 upregulated and 114 downregulated lncRNAs that were differentially expressed between CAD patients and NOCA subjects. The expression profiles of lncRNAs in PBMC were found to alter depending on CAD stage. Our findings support those of previously published studies, in that lncRNA MALAT1 and KCNQ1OT1 expression levels increased, whereas ANRIL and long intergenic non-coding RNA predicting cardiac remodeling expression levels decreased (19,39,40). The lncRNA myocardial infarction-associated transcript (MIAT) has been reported as a biomarker for acute myocardial infarction (AMI), and our findings indicated that the MIAT expression level in patients with myocardial infarction and NOCA subjects was equivalent. Nevertheless, there was altered expression in the setting of AMI, and more specifically, patients with ST had lower levels of MIAT compared to those with non-ST, our results are consistent with the literature (40). These results indicate that these differentially expressed lncRNAs might be implicated in the development and progression of CAD; however, the molecular regulation remains to be elucidated. The altered expression levels of other lncRNAs, such as SENCR and SMILR, which have been reported to be associated with CAD, did not seem to be statistically significant in our study (35,38,41). The probable reason for the heterogeneity of the results may be because the lncRNAs detected differed between studies, whereas the lncRNAs in the study were mainly from PBMC. Therefore, lncRNA changes in PBMC could be related to the development and progression of these diseases. Our results provide supplementary data in this field.

The qRT-PCR results were in agreement with the initial RNA sequencing data in the present study. Interestingly, the upregulated or downregulated expression profiles observed in 48 lncRNAs were in accordance with the subclinical CAD stages. We speculated that these lncRNAs and their co-expressed mRNAs may play an important in CAD development and progression. In addition, the findings of our study indicated the potential of these lncRNAs as candidate biomarkers in clinical use. However, future studies on the molecular mechanism of these lncRNAs are warranted.

Despite providing informative expression profiles of lncRNAs, the present study has several limitations that need to be acknowledged. The participants in our study were from a single center with a high geographic concentration of participants, and the data should be extrapolated to other regions and ethnic groups with caution. Therefore, these results need to be validated using a more extensive population. Also, lncRNA mechanisms in CAD have not been fully elucidated, and further research is warranted to elucidate their biologic functions.

Conclusions

Our preliminary findings indicate that lncRNAs could be used as potential biomarkers of subclinical cardiac abnormalities in the PBMC of CAD patients using the RNA sequencing approach. Our data might provide valuable insight into transcriptome and related regulatory mechanisms in CAD progression. However, further studies are needed to verify our findings and hypothesis.

The article’s supplementary files as