Identification of potentially relevant genes for myocardial infarction using RNA sequencing data analysis.

Myocardial infarction (MI) is a heart disease with high morbidity and mortality rates, thus it is critical to identify genes that serve roles during its pathogenesis. The objective of the present study was to identify potentially relevant genes during the progression of the disease. Blood samples from patients with MI and normal controls (n=3/group) were obtained, the RNA was extracted and cDNA libraries were established. RNA sequencing (RNA-seq) was performed on a HiSeq 2500 platform and fragments per kilobase of exon per million fragments mapped was utilized to calculate the gene expression value following preprocessing of the RNA-seq data. Electronic validation of several identified differentially expressed genes (DEGs) was performed on a Gene Expression Omnibus (GEO) dataset GSE59867 (390 cases and 46 healthy controls). Functional enrichment and protein-protein interaction network analysis was conducted for DEGs. A total of 977 DEGs, including 817 upregulated and 160 downregulated genes were identified in patients with MI. These DEGs were significantly enriched for 'positive regulation of the immune system process,' 'inflammatory response,' 'regulation of I-kappaB-kinase/NF-kappaB signaling' and 'TNF signaling pathway'. A protein-protein interaction network of the top 40 DEGs was used to identify high degree genes, including interferon induced protein with tetratricopeptide repeats 3 (IFIT3), MX dynamin like GTPase 1 (MX1), major histocompatibility complex, class II, DQ α1 (HLA-DQA1), RAR related orphan receptor A (RORA), prostaglandin D2 synthase (PTGDS), cysteine rich protein 2 (CRIP2), collagen type VI α 2 chain (COL6A2) and S100 calcium binding protein P (S100P). The results of validation in the GEO dataset were consistent with the sequencing analysis. A total of eight genes, including IFIT3, MX1, HLA-DQA1, RORA, PTGDS, CRIP2, COL6A2 and S100P may therefore be considered as potentially relevant genes in the pathology of MI.

Introduction

Plaques form in the coronary arteries during atherosclerosis, which impede blood flow. Rupturing of plaques leads to the exposure of arterial blood to pro-thrombotic substances (1). If the intravascular thrombus completely blocks the artery, this causes myocardial infarction (MI), a type of heart disease with a high morbidity rate and mortality rate of 5% (1). Under this circumstance, the left ventricle undergoes a number of changes in shape, size and function, which may lead to heart failure (2). There are multiple factors associated with MI, including sex, age, smoking and complications, including hypertension and diabetes mellitus (3). In primary care, the diagnosis of MI is recorded using Read codes (4) and previous studies have reported lists of Read codes that may be used to identify acute MI (5,6). Clinically, antithrombotic therapy or revascularization is the primary treatment method of type 1 MI, which is defined as myocardial ischemia caused by primary coronary events such as plaque rupture (7). Percutaneous coronary intervention in a timely manner is the primary treatment for acute ST segment elevation MI (7). In addition, various cardiac medications, including aspirin, β-blockers, lipid lowering agents and angiotensin converting enzyme inhibitors have been used in patients with acute MI (8,9).

Although some targeted treatments, including β-blockers and angiotensin converting enzyme inhibitors, have made considerable progress in treating heart failure in patients post-MI (10), further understanding of the potential pathological mechanism of MI is required. An improved understanding of the underlying genes and molecular mechanisms, which increase the risk of MI, is essential to improve its diagnosis, prevention and treatment. Previous reports have indicated some of the potential risk genes and micro-RNAs (miRNAs) in MI. Gao et al (11) suggested that C-C motif chemokine ligand 5, B-cell CLL/lymphoma 3, nuclear receptor subfamily 2 group C member 2, MYC associated factor X, general transcription factor IIIC subunit 2 and nuclear receptor coactivator 7 serve an important role in the development of acute MI. Several higher expression genes, including F-box and leucine rich repeat protein 4, enoyl-CoA hydratase domain containing 3, potassium voltage-gated channel subfamily E regulatory subunit 1, TAO kinase 2, aurora kinase B, ERG and FK506 binding protein 5 and lower expression genes, including myocardial infarction associated transcript, poliovirus receptor-related 3 and PZP were also demonstrated to be present in ST-segment elevation MI compared to non-ST-segment elevation MI (12). In addition, it is reported that elevated cardiac-specific miR-208a in plasma may be a potential biomarker for the early detection of MI (13).

RNA sequencing (RNA-seq) is a progressive method for transcription profiling that provides more precise information about gene expression (14). In the present study, RNA-seq was used to evaluate the gene expression between blood samples from patients with MI and a control group of healthy patients, and a large number of differentially expressed genes (DEGs) were identified. The biological processes, signaling pathways and protein-protein interactions (PPIs) of the identified DEGs were studied. The results were used to further understand the underlying molecular mechanisms of and the potentially relevant genes in the development of MI.

Materials and methods

Collection of patient samples

In the current study, patients with MI were enrolled according to the following conditions: i) Clinical history of ischemic chest pain; ii) dynamic evolution of electrocardiogram and iii) concentration change of myocardial necrosis related serum biomarkers. All patients with MI had no serious infections, endocrine disease, liver and kidney disease, severe neurological disease or blood system disease. Those patients with angina, acute pericarditis, acute pulmonary embolism, aortic dissection, myocarditis, cardiomyopathy, coronary artery malformations and multiple arteritis were excluded. For the normal control group, there was no statistical difference in age, gender, high density lipoprotein and low density lipoprotein compared with the MI group. A total of 3 patients with MI and 3 healthy control were enrolled between December 2015 and February 2016 in the Central Hospital of Tai'an of Shandong Province (Shandong, China). The MI group consisted of one female and two males ages between 46–59 years and the control group consisted of one female and two males aged between 45–56 years. Blood samples (10 ml) were obtained from all patients and then stored in −80°C freezer for subsequent RNA extraction. The clinical results of patients with MI and normal controls are presented in Table I. The current study was approved by the Ethics Committee of Central Hospital of Tai'an of Shandong Province and written informed consent was obtained from all patients.

Table I.

The clinical information of patients with MI and normal controls.

	Sample	Age	Sex	Height (m)	Weight (kg)	BMI	Hypertension history	Diabetes history
MI	Peripheral Blood	50	Female	1.6	68	26.5	No	No
MI	Peripheral Blood	46	Male	1.67	75	26.9	Yes	No
MI	Peripheral Blood	59	Male	1.71	68	23.4	No	No
Normal control	Peripheral Blood	45	Female	1.63	70	26.3	No	No
Normal control	Peripheral Blood	49	Male	1.65	70	25.7	Yes	No
Normal control	Peripheral Blood	56	Male	1.72	71	23.9	No	No

BMI, body mass index.

Patient sample collection, RNA isolation and sequencing

Total RNA were extracted, according to the manufacturer's protocol, using TRIzol reagent (Invitrogen; Thermo Fisher Scientific, Inc., Waltham, MA, USA) and an RNeasy mini kit (Qiagen AB, Sollentuna, Sweden), respectively. The quality and quantity of RNA was determined on a Nanodrop ND-2000 spectrophotometer (Thermo Fisher Scientific, Inc., Pittsburgh, PA, USA). 1.5% agarose gel electrophoresis was used to further verify RNA integrity.

The mRNA was purified from total RNA with a TruSeq RNA library preparation kit (Illumina, Inc., San Diego, CA, USA) and fragmented into sizes of ~200 bp using fragmentation buffer (Illumina, Inc.). The fragmented RNA was converted into double-stranded cDNA in a mixture consisting of reverse transcriptase, buffer, dNTPs and primers (Illumina, Inc.) for 5 min at 65°C followed by 60 min at 42°C and 5 min at 72°C. Following end repair, adapter ligation and gel purification, the ~200 bp cDNA was isolated into 300-nt fragments. A QIAquick PCR kit (Qiagen AB) was used to estimate the relative concentration of the library using a Bioanalyzer 2100 (Agilent Technologies, Inc., Santa Clara, CA, USA) according to the manufacturer's protocol. A HiSeq 2500 platform (Illumina, Inc.) was then applied to perform sequencing. SeqPrep 1.3.2–1 (github.com/jstjohn/SeqPrep) and Sickle V1.0 (github.com/najoshi/sickle) software were used to remove low quality sequences, including ambiguous nucleotides and adaptor sequences of the raw sequencing data. The cleaned sequencing reads were aligned to the University of California Santa Cruz (UCSC, http://genome.ucsc.edu) human reference genome (build hg19) through TopHat 1.3.1 (http://ccb.jhu.edu/software/tophat/index.shtml). The comparison of the reference genome reads with different samples is presented in Table II.

Table II.

Summary of comparison rates of the reads in samples from the control and MI groups with the reference genome.

Sample	Number of reads on the left	Multiple Comparison rate on the left, %	Multiple Comparison rate on the left, %	Number of reads on the right	Comparison rate on the right, %	Comparison rate on the right, %	Total comparison number, (pairs)	Total multiple comparison rate, %	Compare rate of overall consistency, %
Control group
Patient 1	13,635,629	89.0	18.3	13,635,629	86.7	18.3	11,247,785	18.5	81.6
Patient 2	13,441,835	87.4	3.0	13,441,835	85.2	3.0	10,818,632	2.9	79.6
Patient 3	13,679,592	90.1	22.3	13,679,592	88.0	22.3	11,538,272	22.5	83.8
MI group
Patient 1	13,608,241	87.9	10.3	13,608,241	85.6	10.3	11,029,492	10.3	80.1
Patient 2	13,470,345	88.8	21.0	13,470,345	86.6	21.0	11,119,789	21.2	81.8
Patient 3	13,892,369	87.9	8.0	13,892,369	85.7	8.0	11,258,627	8.1	80.1

MI, myocardial infarction.

Identification of DEGs

The fragment assembly and relative expression of the reads with the normalized RNA-Sequencing fragment counts were performed using Cufflinks software 2.2.1 (http://cufflinks.cbcb.umd.edu/). Fragments per kilobase of exon per million fragments mapped was used to determine the transcription abundance of genes. A paired t-test was performed for statistical analysis. A false discovery rate of <0.05 was determined to indicate statistically significant differences.

Electronic validation of DEGs in the Gene Expression Omnibus dataset

The Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo) database was used to validate the expression of the identified DEGs. In the current study, the mRNA profile of 390 MI cases and 46 normal controls was obtained from the GSE59867 dataset. The expression level of the selected targeted DEGs between MI cases and normal controls were compared and the difference of expression level was displayed using box-plots.

Functional enrichment analysis of DEGs

Gene Ontology (GO, geneontology.org/) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG, genome.jp/kegg/pathway.html) pathway enrichment analysis was performed to select DEGs using the Database for Annotation, Visualization and Integration Discovery (DAVID 6.8, david.abcc.Ncifcrf.gov/) online software. The cut-off value for significance was P<0.05.

Construction of PPI network

BioGrid database (thebiogrid.org) was used to retrieve the predicted interactions between the top 40 DEGs (20 upregulated and 20 downregulated) and proteins. Interplayed pairs, which have a confidence score >0.4, were selected to build the PPI network and were visualized using Cytoscape 3.3.0 (cytoscape.org/) software. A node in the PPI network denotes DEGs and the degree of a node represents the number of proteins that are involved in the node.

Results

DEGs associated with MI

A total of 977 DEGs, including 817 upregulated and 160 downregulated genes, were identified in samples from patients with MI. The top 10 upregulated DEGs were histone cluster 1 H2B family member O, cluster of differentiation 177, mast cell expressed membrane protein 1, LOC105376505, vanin 1, annexin A3, disheveled associated activator of morphogenesis 2, LOC105372578, interferon induced protein with tetratricopeptide repeats (IFIT) 1 and butyrophilin like 8. The top 10 downregulated DEGs were major histocompatibility complex, class II, DR β5, major histocompatibility complex, class II, DQ α1 (HLA-DQA1), major histocompatibility complex, class II, DR β1, collagen type VI α2 chain (COL6A2), cysteine rich protein 2 (CRIP2), family with sequence similarity 153 member A, tandem C2 domains, nuclear, protein kinase, Y-linked, pseudogene, SH3 and SYLF domain containing 1 and neural EGFL like 2. A heat map of the top 100 DEGs (50 upregulated and 50 downregulated) is presented in Fig. 1. Additionally, the electronic validation of eight DEGs in the Gene Expression Omnibus dataset (GSE59867) is presented in Fig. 2. MX1 was significantly upregulated in the MI group compared with control (P<0.05), while RORA was significantly downregulated in MI compared with control (P<0.001). IFIT3 (P>0.05) and S100P (P>0.05) were upregulated and PTGDS and HLA-DQA1 were downregulated in the MI compared with the control however, this difference was not significant. The expression trend of these genes was consistent with the RNA sequencing analysis. However, the expression of COL6A2 (P<0.01) and CRIP2 (P>0.05) were upregulated in the MI group compared with control, which was not in line with the RNA sequencing analysis.

Figure 1.

Heat map of clustering analysis of top 100 differentially expressed genes in myocardial infarction obtained using heatmap.2 (http://127.0.0.1:28428/library/gplots/html/heatmap.2.html).

Figure 2.

Electronic validation of eight differentially expressed genes in the Gene Expression Omnibus dataset (GSE59867). Electronic validation of MX1, IFIT3, S100P, RORA, PTGDS, COL6A2, CRIP2 and HLA-DQA1 in the control and case groups. MX1, MX dynamin like GTPase 1; IFIT3, interferon induced protein with tetratricopeptide repeats 3; S100P, S100 calcium binding protein P; RORA, RAR related orphan receptor A; PTGDS, prostaglandin D2 synthase; COL6A2, collagen type VI alpha 2 chain; CRIP2, cysteine rich protein 2; HLA-DQA1, major histocompatibility complex, class II, DQ α1.

Enrichment analysis of DEGs

GO enrichment analysis and KEGG pathway analysis demonstrated that the identified DEGs are significantly involved in the ‘positive regulation of the immune system process,’ ‘inflammatory response,’ ‘adaptive immune response,’ ‘regulation of I-kappaB/NF-kappaB signaling,’ ‘hematopoietic cell lineage’ and ‘TNF signaling pathway’. Figs. 3 and 4 indicate the results of GO enrichment analysis and KEGG pathway analysis, respectively. The KEGG map of the nuclear factor (NF)-κB and tumor necrosis factor (TNF) signaling pathways are presented in Figs. 5 and 6, respectively.

Figure 3.

GO enrichment analysis of differentially expressed genes in myocardial infarction. GO, Gene Ontology.

Figure 4.

Kyoto Encyclopedia of Genes and Genomes pathway analysis of differentially expressed genes in myocardial infarction. NF, nuclear factor; TNF, tumor necrosis factor; HIF-1, hypoxia-inducible factor 1; HTLV-1, human T-lymphotropic virus-1.

Figure 5.

DEGs enriched KEGG map of the NF-κB signaling pathway in myocardial infarction. Colored rectangle represents the DEG the red and green represents upregulation and downregulation, respectively. DEGs, differentially expressed genes; NF-κB, nuclear factor-κB; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Figure 6.

DEGs enriched KEGG map of the TNF signaling pathway in myocardial infarction. Colored rectangle represents the DEG, the red and green represents upregulation and downregulation, respectively. DEGs, differentially expressed genes; TNF, tumor necrosis factor; KEGG, Kyoto Encyclopedia of Genes and Genomes.

PPI network of DEGs

Following the integration of protein information in the BioGrid database with the top 40 DEGs (20 upregulated and 20 downregulated), a PPI network of the identified DEGs was established. The primary nodes with high degrees (≥10) were identified as follows: IFIT3 (degree, 66); MX dynamin like GTPase 1 (MX1; degree, 34); HLA-DQA1 (degree, 28); RAR related orphan receptor A (RORA; degree, 23); prostaglandin D2 synthase (PTGDS; degree, 15); CRIP2 (degree, 15); COL6A2 (degree, 12); and S100 calcium binding protein P (S100P; degree, 10) (Fig. 7).

Figure 7.

Protein-protein interaction network of the top 40 DEGs in myocardial infarction. Rectangle and ellipse represents the DEG and protein, respectively. Red represents an upregulated gene and blue represents a downregulated gene. DEGs, differentially expressed genes.

Discussion

MI is a major cause of mortality worldwide and affects patient health and quality of life (1). In the current study, RNA-seq technology was applied to identify DEGs to investigate the roles served by important genes in MI. A total of 977 DEGs, with 817 upregulated genes and 160 downregulated genes, were identified. Following GO annotation and KEGG pathway analysis, it was determined that all of the DEGs were enriched with the ‘positive regulation of the immune system process,’ ‘inflammatory response,’ ‘adaptive immune response,’ ‘regulation of I-kappaB/NF-kappaB signaling,’ ‘hematopoietic cell lineage’ and ‘TNF signaling pathway’. PPI analysis also identified several high-degree DEGs, suggesting that they serve an important role in MI.

There are a small number of immune cells in the adult mammalian heart and these immune cells expand in size following MI (15). Necrotic cardiomyocytes release danger signals, activating immune cells, which clear the matrix debris and dead cells, allowing the activation of cardiac repair (15). The extent of the inflammatory response may be used to identify the patient's outcomes, including heart function healing, acute myocardial rupture or dilatation (16). Post-MI inflammatory responses serve a role in cardiac repair and determine the infarct size and quality of repair (17). In the current study, DEGs were significantly enriched in the biological processes of GO annotation, which further demonstrated the role they serve in the immune system and inflammatory response during MI.

NF-κB activation serves a critical role in coronary microembolism-induced cardiac dysfunction and advanced heart failure (18). For example, NF-κB activation and the subsequent inflammatory response contribute to cardiac dysfunction and maladaptive left ventricular remodeling following MI (19). The TNF signaling pathway protects the heart by ischemic preconditioning; however, this protective effect depends on the concentration of TNF-α (20,21). Therefore, the appropriate concentration of TNF-α is a critical factor in the determination of the outcome of cardiac injury, particularly following acute MI (22). In the present study, the DEGs were also involved in the NF-κB and TNF signaling pathways, thus, these genes may serve roles via these pathways in the development of cardiac dysfunction in MI.

To further study the function of the identified DEGs in MI, a PPI network was constructed. In the PPI network, several hub nodes were identified, which suggested that the hub genes may serve crucial roles in the pathology of MI. IFIT3 is a gene associated with innate immunity (23). The expression of IFIT3 is increased in transgenic mouse myocytes (24) and is upregulated in animal models of pathological cardiac hypertrophy (25). In the current study, the expression of IFIT3 was upregulated in patients with MI, which suggests it may function as an immunity response gene in MI.

Chen et al (10) indicated that MX1 was associated with MI. In the present study, MX1 was upregulated in MI, which suggested that it may be involved in the pathology of the disease. HLA-DQA1 is a member of the MHC Class II family and is a downregulated gene in patients with cardiorenal syndrome (26). In the present study, the decreased expression of HLA-DQA1 in the blood of patients with MI was demonstrated. This suggested that HLA-DQA1 may serve a notable role in the pathology of MI.

RORA is a downstream effector molecule of circadian locomotor output cycles kaput genes (27,28). Results of genome-wide association studies (GWAS) or GWAS meta-analyses have demonstrated that RORA is one of the primary risk factors for cardiovascular metabolism (29). Additionally, RORA serves a role in the regulation of the inflammatory response (30). In the current study, the expression of RORA was downregulated, which suggested that RORA may be associated with MI.

PTGDS is involved in the biosynthesis of prostaglandin D2 and promotes cardiomyocyte survival (31,32). It serves a role in numerous diseases, including inflammatory, infectious and cardiovascular disease (33). Furthermore, PTGDS protein has been designated as a circulating marker in coronary artery disease (34,35). In the present study, the expression of PTGDS was decreased, which indicated that it may be a biomarker of MI.

CRIP2 is a cysteine-rich protein, which is expressed in the heart endothelium during development and also in the adult heart (36,37). CRIP2 is considered as a heart vascular marker (38). In the present study, CRIP2 expression was decreased in MI. However, the electronic validation result was not consistent with the sequencing analysis. This may be due to the size of the dataset used in the current study, thus larger datasets should perhaps be used for further investigations.

Cardiomyocyte death due to long-term ischemia induces infarction and robust adhesion molecule upregulation (39). In addition, COL6A2 is a cardiac fibroblast-enriched gene and encodes adhesion proteins (40,41). In the current study, the expression of COL6A2 was decreased in patients with MI, which indicated that COL6A2 may be involved in its development. However, it was upregulated in the results of electronic validation. This again may be due to the size of the dataset used in the current study and larger datasets should perhaps be used for further investigations.

Research has also demonstrated that S100P is associated with myocardial injury in patients and rat models of MI (42). The expression of S100P was increased in the serum of patients with acute coronary syndrome and was associated with infarct size (42). Similarly, the expression of S100P was increased in the current study, which further demonstrated its role in MI.

In conclusion, several potentially relevant genes, including IFIT3, MX1, HLA-DQA1, RORA, PTGDS, CRIP2, COL6A2 and S100P, in MI were identified in the current study and these may serve important roles in the development of MI. However, there are certain limitations to the present study. The number of blood samples used for RNA-seq was small and a larger sample size is required. This is a preliminary study investigating MI using bioinformatics analysis. In vitro gene expression validation (such as real-time quantitative polymerase chain reaction), animal models (including gene overexpression and gene knockdown), and epigenetics studies (such as DNA methylation) are required to further investigate the function of identified genes in MI. However, the current study provides a basis for increasing the understanding of the gene expression profile in MI.