Introduction
Coronary artery disease (CAD), also termed coronary arteriosclerosis, is one of the most common types of heart disease (1). Morbidity and mortality of CAD has increased in recent years, reducing the quality of life of patients and continuing to present an important socioeconomic problem (2). Early diagnosis of CAD is difficult, and the mechanism of its onset and progression is complicated (3). Coronary artery bypass graft surgery and drug treatments are the primary treatment strategies for CAD (4). Although advances have been achieved with regards to CAD treatment, CAD is a health burden that remains to be solved. Therefore, it is important to investigate the mechanisms of CAD progression, and explore potential methods of CAD diagnosis and treatment.
A previous study demonstrated that CAD progression may be driven by immune factors, traditional risk factors, such as high blood pressure, diabetes, hyperlipidemia and smoking, and other novel risk factors; for example, high blood pressure is involved in the cardiovascular outcome of patients with diabetes and CAD (5). Apolopoprotein B is a novel CAD-associated protein that has been identified to stimulate the proliferation of coronary artery smooth muscle cells and promote their movement into the subendocardial layer to enhance the progression of CAD (6). Interleukin-18 (IL-18) is an independent predictor of the cardiovascular events in patients with CAD (7). Additionally, increasing evidence demonstrated that microRNAs (miRNAs/miRs) are crucial in CAD progression (8). The overexpression of miR-1 downregulates B-cell lymphoma 2 (Bcl-2) expression levels by targeting the 3′-untranslated region of Bcl-2 in cardiac muscles and is, thus, closely associated with ischemic injury (9). miR-210 expression targets caspase-8-associated protein 2 in ischemic preconditioning and may contribute to the survival of stem cells; therefore, protecting from ischemic injury (10). Although numerous risk miRNAs and crucial genes have been associated with CAD progression, the mechanism of CAD remains largely unknown.
In previous studies, the molecules in different cell types were investigated. Due to its interactive and dynamic properties, blood composition it is often closely associated with alterations that occur during the progression of disease and responses to injury (11). Therefore, whole blood cells may be useful samples for CAD disease research, and may be an alternative to tissue biopsy (12). Additionally, adhesion of circulating leukocytes was confirmed to be an important step for the development of CAD (13). Furthermore, the platelets of CAD patients may be easily activated when coronary blood flow is increased (14). These samples are important for CAD research. The use of computational methods may allow for a more thorough investigation of the interactions between the molecules in different cell types, and thus lead to the identification of novel factors that contribute to CAD progression. Zhang et al (15) used GSE12288 microarrays to identify growth factor receptor-bound protein 2 and heat shock protein family A (Hsp 70) member 8 as the key genes for CAD development. Chen et al (16) used GSE28858 microarrays to analyze the key miRNAs (miR-545 and miR-585) associated with CAD. In addition, Hua et al (17) identified the CAD-associated miRNA clusters using the same microarray data. The present study aimed to elucidate the key genes and miRNAs associated with CAD progression.
Materials and methods
Data resources and preprocessing
The gene expression profiles of GSE20680 (18) and GSE12288 (19) were downloaded from the gene expression omnibus (GEO) database in NCBI (National Center for Biotechnology Information; http://www.ncbi.nlm.nih.gov/geo/) based on the platforms of GPL4133 Agilent-014850 Whole Human Genome Microarray 4×44K G4112F (Feature Number version) and GPL96 [HG-U133A] Affymetrix Human Genome U133A Array, respectively. The dataset of GSE20680 contained 143 CAD and 52 control samples, and that of GSE12288 included 110 CAD and 112 control samples. In addition, the miRNA expression profile data of GSE28858 (20) was comprised of 12 samples from patients with premature CAD and 12 age- and gender-matched healthy control samples. It was downloaded from the GEO database in NCBI based on the platform of GPL8179 Illumina Human v2 MicroRNA expression beadchip.
The gene profile data of GSE20680 was preprocessed using Agilent Feature Extraction software (version 9.5.3.1; Aglient Technologies, Inc. Santa Clara, CA, USA) (21). The CEL file data of GSE12288 was preprocessed using the robust multi-array analysis method from the affy package in R (22). If a gene had several probes the mean expression value was selected. Additionally, miRNA IDs from the preprocessed expression matrix of GSE28858 were transformed into the miRNA symbols.
Differentially expressed gene (DEG) screening and enrichment analysis
The DEGs in CAD samples were compared with the control samples from the two gene expression profile datasets using a t-test in the limma package in R software (23). P<0.05 and a log2 fold-change of 0.1 were selected as thresholds to indicate a statistically significant difference.
In addition, the significant biological functions and pathways of the screened DEGs in GSE20680 and GSE12288 were analyzed using Database for Annotation, Visualization, and Integrated Discovery (DAVID) (24) from the Gene Ontology (GO) (25) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (26) databases with P<0.05.
Protein-protein interaction (PPI) network construction and module selection
Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) is a database of known and predicted protein interactions that may aid in the comprehensive description of cellular mechanisms and functions (27). The PPI network of the selected DEGs was constructed using the STRING database. Interaction pairs with a PPI score of 0.7 were selected for the construction of the final network.
The modules from the constructed PPI network were selected using the ClusterOne plugin in Cytoscape software (version 2.8) (28). In addition, the significant interaction pathways of DEGs with P<2.727×10−9 were analyzed using DAVID with P<0.05 indicating a statistically significant difference.
Regulatory network construction
The transcriptional associations between transcription factors (TFs) and target genes are of great biological significance, and may aid in the analysis of numerous physiological activities (29). The TFs and target genes from the selected DEGs in the two profile datasets were analyzed based on the information of TF-target genes stored in the UCSC database (30). Also, the regulatory network of TFs-target genes was constructed using the Cytoscape software (version 2.8) (31).
Enrichment analysis of common DEGs
The screened DEGs that appeared in the two datasets (GSE20680 and GSE12288) were considered common DEGs. The significant biological functions and pathways of the selected common DEGs were analyzed using DAVID (24) in GO (25) and KEGG (26) database, respectively. P<0.05 was selected as the cut-off criteria for including a statistically significant difference.
miRNAs screening and regulatory network construction of miRNA-targets
The differentially expressed miRNAs in CAD samples from the GSE28858 dataset were screened and compared with the control samples using the t-test in the limma package in R (23). P<0.05 and a log2 fold-change of 0.1 were selected as the thresholds for indicating statistically significant differences.
In addition, miRecords (32) and MirWalk (33) are two databases that integrate miRNA-target interactions with the experimental validated target genes of miRNAs. The target genes that are regulated by the selected differentially expressed miRNAs were predicted based on the miRecords and MirWalk databases. Genes that are present in one of the two or in both databases were selected for inclusion in the current study.
CAD-associated miRNA-target selection and enrichment analysis
Genes that appeared in the predicted miRNA-target interactions and in the CAD-associated dataset from the Comparative Toxicogenomics Database (CTD) (34) were confirmed to be the CAD-associated genes. The significant biological functions and pathways of the predicted miRNA-targets were analyzed using DAVID (24) in GO (25) and KEGG (26) databases, respectively with P<0.05 indicating a statistically significant difference.
Analysis of CAD-associated DEGs and miRNAs
To investigate the associations between DEGs (GSE20680 and GSE12288) and differentially expressed miRNAs (GSE28858) in the CAD samples, the total genes and miRNAs were integrated to screen for CAD-associated differentially expressed miRNA-target genes. Additionally, the significant functions of miRNAs were analyzed using DAVID (24) in the GO (25) database with P<0.05 indicating a statistically significant difference. A crosstalk network of miRNAs involved in the same biological processes was constructed, P<0.0001 indicating a statistically significant difference.
Results
DEG screening and enrichment analysis
The GSE20680 dataset was comprised of 270 DEGs (167 upregulated and 103 downregulated), and the GSE12288 dataset was comprised of 2,268 DEGs (534 upregulated and 1,734 downregulated).
The enriched significant GO terms and KEGG pathways of DEGs in GSE20680 and GSE12288 are indicated in Tables I and II, respectively. The upregulated DEGs in GSE20680 were involved in cell chemotaxis and positive regulation of the defense response, while the downregulated DEGs were involved in the positive regulation of B cell activation (Table IA). Additionally, the significant pathways of upregulated DEGs included oxidative phosphorylation, cardiac muscle contraction and metabolic pathways, while the downregulated genes were enriched in primary immunodeficiency and gap junction pathways (Table IIA). Conversely, the upregulated DEGs in the GSE12288 dataset were involved in system and multicellular organismal processes, while the downregulated DEGs were involved in cell activation, activation of immune response and the immune system process (Table IB). In addition, the significant pathways of upregulated DEGs included neuroactive ligand-receptor interaction and dilated cardiomyopathy, while the downregulated genes were enriched in the phagosome and Fc γ R-mediated phagocytosis pathways (Table IIB).
PPI network construction and module selection
The PPI networks of DEGs in the two datasets were constructed. The PPI network of DEGs in GSE20680 contains 68 nodes and 60 interaction pairs (Fig. 1A) while the PPI network of DEGs in GSE12288 includes 1,558 nodes and 7,695 pairs. The four modules from the PPI network of DEGs in GSE12288 were selected for further analysis. The significant pathways of DEGs in the selected four modules are indicated in Table III. Genes in module 1 were involved in the spliceosome and RNA polymerase pathway, genes in module 2 were enriched in proteasome and oocyte meiosis pathways, and genes in module 4 participated in oxidative phosphorylation and cardiac muscle contraction (Table III). Module 3 did not contain any enriched pathway genes.
Regulatory network construction
No TF target genes were obtained from the GSE20680 dataset. However, a total of 3,400 TF target genes were obtained from the GSE12288 dataset, including 10 TFs [CCAAT enhancer-binding protein β; CCCTC-binding factor (zinc finger protein); interferon regulatory factor 1 and 3 (IRF1 and 3), jun D proto-oncogene, MYC-associated factor X (MAX); RAD21 cohesin complex component; retinoid X receptor α; signal transducer and activator of transcription 3 (STAT3); YY1 transcription factor] and 1,747 target genes. The expression level of these 10 TFs in CAD samples were analyzed compared with the control samples (Fig. 1B). The expression levels of these TFs in CAD samples were lower than that in control samples, indicating that they were all downregulated genes.
In addition, regulatory networks between the 10 TFs and their target genes [IRF2 and cell death-inducing DFFA-like effector b (CIDEB)] were constructed. The results indicated that the number of downregulated genes regulated by the TFs was greater than the number of upregulated genes.
Enrichment analysis of common DEGs
A total of 41 common DEGs between GSE20680 and GSE12288 dataset were identified, including IRF2, fibrinogen-like 2 (FGL2), CIDEB and ribosomal protein S4, Y-linked 1 (RPS4Y1). The enriched GO terms and KEGG pathways of these common DEGs are indicated in Table IV. The significant GO terms of common DEGs included polysaccharide and carbohydrate derivative metabolic processes, as well as leukocyte mediated immunity (Table IV), while the enriched pathways of common genes were amino and nucleotide sugar metabolism, and protein digestion and absorption (Table IV).
miRNA screening and regulatory network construction of miRNA-targets
A total of 214 differentially expressed miRNAs (102 upregulated and 112 downregulated) in CAD samples were screened. Fig. 2A represents a heat map of miRNA expression levels. A total of 71 miRNAs were confirmed to regulate 455 target genes based on the miRRecords and MirWalk databases. Finally, 640 interactions between the 71 miRNAs and their target genes were determined.
CAD-associated miRNA target selection and enrichment analysis of miRNA targets
A total of 402 common genes involved in the expression of 69 miRNAs were identified based on a comparison between the 455 target genes that are regulated by the 71 miRNAs and the genes stored in the CTD database.
In addition, the significant GO terms and KEGG pathways of the 69 miRNAs are indicated in Table V. GO terms, including major histocompatibility complex (MHC) protein assembly, MHC class I protein complex assembly and peptide antigen assembly with MHC protein complex, were identified as significantly enriched in the miRNAs tested (Table V). Only one miRNA was identified to participate in the allograft rejection pathway (Table V).
The regulatory network between miRNAs associated with CAD and their target genes was also constructed (Fig. 2B). It was identified that miR-455-5p, miR-455-3p and miR-1257 may regulate numerous miRNAs and target genes, including pro-opiomelanocortin (POMC), toll-like receptor 4 (TLR4), IL10, activating transcription factor 6 (ATF6), and calreticulin (CALR).
Analysis of CAD associated DEGs and miRNAs
A total of 5 and 138 miRNA-target interaction pairs were obtained from the GSE20680 and GSE12288 dataset, respectively. However, crosstalk analysis of the 5 miRNA-target interaction pairs was not significant (Table VI). Additionally, the crosstalk network of miRNAs indicated that miRNAs were predominantly involved in MHC class I protein complex assembly, peptide antigen assembly with MHC protein complex, peptide antigen assembly with MHC class I protein complex, regulation of metanephric cap mesenchymal cell proliferation, positive regulation of metanephric cap mesenchymal cell proliferation, and MHC protein complex assembly (Fig. 2C).
Discussion
CAD is one of the most common types of heart disease, exhibiting increasing morbidity and mortality rates. CAD may reduce the quality of life of the patients and is an important socioeconomic problem (1,2). The mechanism of CAD progression is complicated; thus, it is important to investigate the disease mechanism, in addition to exploring different methods for CAD diagnosis and treatment. In the present study, microarrays were used to screen for CAD-associated genes and miRNAs. The results indicated that IRF2 and CIDEB, which were regulated by STAT3 and MAX, were common DEGs for CAD. In addition, miR-455-5p, miR-455-3p and miR-1257, which are involved in the MHC protein complex assembly pathway and peptide antigen assembly with MHC class I protein complex pathway, may regulate numerous miRNAs and target genes, including POMC, TLR4, IL10, ATF6 and CALR.
The results of the current study indicated that the CAD-associated gene POMC was a target of the upregulated miR-455-5p, which enriched in the MHC protein complex assembly pathway. POMC is a polypeptide hormone precursor that functions as a feeding suppressant and is similar to leptin (35). It has been determined that POMC neurons were targeted by leptin in the hypothalamus to promote the synthesis of α-MSH from POMC (36). Additionally, α-MSH acts on the melanocortin 4 receptor to induce seeding suppression, thus protecting the body from obesity (37). Logue et al (38) reported that obesity frequently led to fatal CAD. Therefore, POMC may be associated with CAD progression. Conversely, miR-455-5p has been identified to target scavenger receptor class BI and reduce high density lipoprotein cholesterol (HDL-C) uptake (39). Lower HDL-C is a predictor for CAD risk (40). Based on the current study, upregulation of miR-455-5p may reduce the progression of CAD by targeting the POMC gene via the MHC protein complex assembly pathway.
The present study demonstrated that the TLR4 gene was the common target of the downregulated miR-455-3p and the upregulated miR-455-5p. TLR4 is a member of the Toll-like receptor family, which is important for pathogen recognition and the activation of innate immunity (41). Otsui et al (42) determined that TLR4 was highly expressed in smooth muscle cells in patients with atherosclerotic arteries, and TLR4-mediated inflammatory activation of human coronary artery endothelial cells via lipopolysaccharide (43). Therefore, TLR4 may be involved in CAD progression. It has been demonstrated that upregulated miR-455-3p was involved in the acute myocardial infarction (44), while myocardial infarction was the pathological basis for ventricular remodeling in CAD (45). Therefore, miR-455-3p may be associated with CAD progression via myocardial infarction. Based on the results of the present study, downregulated miR-455-3p may inhibit CAD progression by targeting the TLR4 gene.
The current study indicates that CALR was the only target gene for the upregulated miR-1257, implying their respective importance in CAD progression. CALR is a multifunctional protein that acts as a major Ca2+-binding protein in the lumen of the endoplasmic reticulum (46). CALR is also associated with the myocardial hypertrophy. Overexpression of CALR may induce the dilated cardiomyopathy (47). Additionally, myocardial hypertrophy is one of the pathophysiological alterations that occur during CAD (48). Therefore, CALR may be involved in CAD development. Notably, the role of miR-1257 in CAD remains to be fully investigated. However, Kamiński et al (49) reported that miR-1257 is a cardiovascular disease-associated miRNA that has an A binding site. Based on the observations of the present study, it is possible that miR-1257 may be a key regulator of CAD progression by regulating the CALR gene.
The current study also indicated that the downregulated IRF2 and CIDEB genes were the common DEGs for CAD. IRF2 was regulated by STAT3 while CIDEB was regulated by MAX. IRF2 is a member of the interferon regulatory transcription factor family of proteins that have a transcriptional binding site for STAT3 (50). The roles of IRF2 and CIDEB in CAD have not been fully elucidated in previous studies. However, co-operative IRF1 (the homologue of IRF2) and IL-6 expression was associated with myocardial infarction (51). In addition, IRF1 inhibited the differentiation of T helper cells from CD4+ T cells in the peripheral blood in cases of acute coronary syndrome, indicating their involvement in the development of this syndrome (52). STAT3 was also reported to contribute to heart failure, which is associated with CAD (53). Therefore, based on the results of the current study IRF-2 may be important in CAD development regulated by STAT3 while CIDEB may be a novel factor that is regulated by MAX in CAD.
Additionally, the observations of the present study indicate that the selected significant miRNAs (miR-455-5p, miR-455-3p and miR-1257) were involved in the MHC protein complex assembly pathway and peptide antigen assembly with MHC class I protein complex pathway. A previous study revealed that the T cell receptor may only recognize and bind to the peptide fragments of MHC (54). Higher numbers of CD4+ T cells may promote the progression of atherosclerosis (55). In addition, the interaction between dendritic cells and T cells contributed towards the process of atherosclerosis (56). The gathered dendritic cells may secrete tumor necrosis factor-α to induce CD4+ T cells to produce tumor necrosis factor superfamily member 10 (TNFSF10, also known as TRAIL) (57). The TRAIL may combine with its receptors (TRAIL-R1 or TRAIL-R2), which are located on the vascular smooth muscle cells surface, and then induce the apoptosis of smooth muscle cells (58). Therefore, it is possible that miR-455-5p, miR-455-3p and miR-1257, may be important for the progression of CAD by participating in the MHC protein complex assembly pathway and peptide antigen assembly with MHC class I protein complex pathway.
The screened DEGs and TFs were enriched in various GO terms, including carbohydrate metabolic process, and KEGG pathways including cardiac muscle contraction and protein digestion and absorption. Therefore, free fatty acid metabolism was the key factor in CAD patients, which may regulate the coupling between carbohydrate oxidation and glycolysis (59). In addition, Fichtlscherer et al (60) also determined that certain critical miRNAs, such as miR-133 and miR-208a, were significantly enriched in cardiac muscle, and further participated in CAD disease. Therefore, the screened target genes and their associated TFs may participate in CAD development by being enriched in the aforementioned pathways.
In conclusion, the present study suggests that miR-455-5p reduces the progression of CAD by targeting POMC while miR-455-3p inhibits CAD by targeting TLR4. miR-1257 may be a key regulator for CAD by targeting CALR. Additionally, IRF-2, which is regulated by STAT3, may be important in CAD development while CIDEB, which is regulated by MAX, may be a novel factor in CAD progression. The current study may provide a basis for future research on the mechanism of CAD progression. There were however limitations to the present study. For example, the expression of the identified molecules and the machinery of the disease process should be verified in patients with CAD using western blot analyses and reverse transcription-quantitative polymerase chain reaction. Therefore, further experimental and clinical studies are required to confirm the results presented of the present study.