A comparison of gene expression profiles in patients with coronary artery disease, type 2 diabetes, and their coexisting conditions

Background To support a hypothesis that there is an intrinsic interplay between coronary artery disease (CAD) and type 2 diabetes (T2D), we used RNA-seq to identify unique gene expression signatures of CAD, T2D, and coexisting conditions. Methods After transcriptome sequencing, differential expression analysis was performed between each disordered state and normal control group. By comparing gene expression profiles of CAD, T2D, and coexisting conditions, common and specific patterns of each disordered state were displayed. To verify the specific gene expression patterns of CAD or T2D, the gene expression data of GSE23561 was extracted. Results A strong overlap of 191 genes across CAD, T2D and coexisting conditions, were mainly involved in a viral infectious cycle, anti-apoptosis, endocrine pancreas development, innate immune response, and blood coagulation. In T2D-specific PPI networks involving 64 genes, TCF7L2 (Degree = 169) was identified as a key gene in T2D development, while in CAD-specific PPI networks involving 64 genes, HIF1A (Degree = 124), SMAD1 (Degree = 112) and SKIL (Degree = 94) were identified as key genes in the CAD development. Interestingly, with the provided expression data from GSE23561, the three genes were all up-regulated in CAD, and SMAD1 and SKIL were specifically differentially expressed in CAD, while HIF1A was differentially expressed in both CAD and T2D, but with opposite trends. Conclusions This study provides some evidences in transcript level to uncover the association of T2D, CAD and coexisting conditions, and may provide novel drug targets and biomarkers for these diseases. Electronic supplementary material The online version of this article (doi:10.1186/s13000-017-0630-7) contains supplementary material, which is available to authorized users.


Background
Type 2 diabetes (T2D) and coronary artery disease (CAD) often coexist and cause substantial public health and economic burden world-wide. CAD has long been established as a complication of T2D. The plaque formation in T2D patients may narrow the coronary arteries and thus predispose the occurrence of heart attack. It is assumed that there is an intrinsic interplay between T2D and CAD, in the form of shared etiology and pathophysiological mechanisms. The two diseases have shared common risk factors such as, age, gender, anthropometric, metabolic, socioeconomic and lifestyle variables, as well as psychosocial stress and environmental pollutant exposure. In addition, both diseases are characterized by a chronic inflammatory process [1,2] and disorders of the coagulation system [3].
T2D has been associated with increased risk of cardiovascular disease and death [4,5]. The underlying mechanisms may involve a complex interplay between genes predisposing to insulin resistance and those independently regulating lipid metabolism, coagulation processes and biological responses of the arterial wall [6]. The shared susceptibility regions (bin 9.3 and 6.5) were observed across T2D, obesity and CAD by Wu et al, suggesting the possibility of shared pathophysiology and risk through genetic pleiotropy [7], which may account for their frequent coexistence. However, there has been limited success in correlating T2D with CAD in terms of pathophysiologic changes up to now [8].
Peripheral blood gene expression profile has been used to reflect pathological conditions in a variety of diseases [9][10][11][12][13][14]. The increasing genomic information has provided an opportunity to better understand the complex biological processes of diseases, especially after the emergence of high throughput technologies. Previous studies have reported the gradual change in circulating gene expression profiles in patients with different extent of CAD [14].
Therefore, to support a hypothesis of an intrinsic interplay between T2D and CAD, we compared the peripheral blood gene expression profiles of T2D, CAD and coexisting conditions, to better understand the association of the three metabolic disorders.

Patients
All patients (2 T2D, 2 CAD, and 6 T2D + CAD) were recruited from the Third Municipal Hospital of Shijiazhuang City between March 2007 and December 2009. The diagnosis of T2D was according to World Health Organization criteria [15]. CAD was diagnosed with imaging techniques to detect flow-limiting coronary artery stenosis [16]. Patients met both of the above inclusion criteria were defined as T2D + CAD. The age-and race-matched patients (n = 7) attending the outpatient department were recruited as control during the study period. None of these patients had previous diagnosis of dyslipidemia, abnormal glucose tolerance, high blood pressure, or any illness. Demographic data and medication of the study population are summarized in Additional file 1: Table S1. The study was approved by the Institutional Review Board of the Third Municipal Hospital of Shijiazhuang City and all subjects provided written informed consent.

RNA isolation and sequencing
Peripheral blood mononuclear cells (PBMCs) were isolated from ethylene diamine tetraacetic acid (EDTA) anticoagulated whole blood using Ficoll-Hypaque gradients. Total RNA was extracted using a Trizol reagent (Invitrogen, Carlsbad, CA, USA). The quality and quantity of RNA were evaluated on a Nanodrop ND-2000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA). Isolation of messenger RNA (mRNA) was carried out using a TruSeq RNA library preparation kit (Illumina, San Diego, CA) according to the manufacturer's instruction. The products were subsequently fragmented into sizes of around 200 bp and subjected to double-stranded cDNA synthesis. A HiSeqTM 2500 platform (Illumina) was applied to perform sequencing.

Differential expression analysis
TopHat v1.3.1 software [17] was used to align raw sequencing reads to the UCSC human reference genome (Build hg19). The original alignment file was processed to measure transcript abundance using Cufflinks v1.0.3 software [18]. Transcript abundance of each gene was determined by calculation of Reads per kilobase of exon per million mapped reads (RPKM). The paired t-tests were performed to identify differentially expressed genes. P <0.05 was selected as the criteria for significant differences. Hierarchical clustering of differentially expressed genes was performed using the "pheatmap" function of the R/Bioconductor package [8].

Functional enrichment analysis of differentially expressed genes
Gene ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed to annotate the biological function of the differentially expressed genes using the online software GENECODIS [19]. A cut-off of FDR was defined at 0.05.

Protein-protein interactions (PPIs) network construction
To reveal the interactions of selected genes at molecular level, PPIs network was established based on the online database [20,21]. Biological General Repository for Interaction Datasets (BioGRID) (http://thebiogrid.org/) was used to construct PPI networks, and the distribution characteristics of selected genes in the PPI network were visualized using Cytoscape software [22]. Nodes in the PPI network represent proteins, while edges represent interactions between two proteins.

Differential expression analysis between CAD and control group
A total of 488 genes were significantly differentially expressed in CAD when compared with the control group, which include 400 up-regulated and 88 down-regulated genes. Signal transduction (GO: 0007165, P = 6.14E-11) and blood coagulation (GO:0007596, P = 8.96E-09) were significantly enriched. Chemokine signaling pathway was one of the most significantly enriched pathways (P = 6.14E-11)( Table 2).

Correlation among three disease states
A strong overlap of 191 genes was identified among genes that differentially expressed across CAD, T2D, and T2D + CAD groups when compared with the control group.
Additionally, pairwise Spearman's correlation coefficient was calculated to assess the correlation among three disease states using the differentially expressed genes in each disease. The results showed a significant correlation among CAD, T2D, and T2D + CAD (p <0.0001), suggesting the possibility of shared pathophysiology of the diseases. The correlation between T2D and T2D + CAD was the most significant in the peripheral blood gene expression profiles (Spearman's rho = 0.6757, p <0.0001).

Disease-specific PPI network
The pathology of T2D with CAD shares much in common with that of CAD and T2D. We selected 64 genes that were significantly differentially expressed in T2D or T2D with CAD, but were not significantly differentially expressed in CAD to construct T2D-specific PPI networks. The information of 7 genes wasn't available in BioGRID database. When removing duplicated edges, self-loops nodes, colocalization edges, finally 47 genes  were involved in T2D-specific PPI networks, including 1216 nodes and 1579 edges. The significant hub proteins contained TCF4 (Degree = 169), SKP1 (Degree = 164) and UBE2W (Degree = 75) (Fig. 1), suggesting their important role in the development of T2D. To construct CAD-specific PPI networks, we selected 69 genes that were significantly differentially expressed in CAD or T2D with CAD, but were not significantly differentially expressed in T2D. The information of 9 genes wasn't available in BioGRID database. When removing duplicated edges, self-loops nodes, colocalization edges, finally 54 genes were involved in T2D-specific PPI networks, including 943 nodes and 1085 edges. The significant hub proteins contained HIF1A (Degree = 124), SMAD1 (Degree = 112) and SKIL (Degree = 94) (Fig. 2), suggesting their important role in the development of CAD.

The verification of gene expression in GSE23561
, HIF1A, SMAD1, and SKIL, specifically differentially expressed in CAD, and also hub genes in the CAD-specific PPI networks, was selected to confirm the above results. With the provided expression data from GSE23561, HIF1A, SMAD1, and SKIL were all significantly upregulated in CAD, and the expression of SMAD1 and SKIL remained unchangeable in T2D as compared with normal control, which was consistent with our results. Differently, HIF1A was both differentially expressed in CAD and T2D, but with opposite trends (Fig. 3).

Discussion
In this study, we used RNA-seq to identify unique peripheral blood gene expression signatures of T2D, CAD, and coexisting condition. Previous literature have revealed that there is an intrinsic interplay between T2D and CAD, while the detailed mechanism remains unclear. Toward this end, we compared the gene expression profiles of T2D, CAD and coexisting condition to show the association in-between them, and tried to explain the shared pathophysiology.
Here, PBMCs are usually selected to monitor posttranslational modifications relevant to many diseases [24], updating the underlying molecular events of diseases. Furthermore, the feature of minimal invasion makes differentially expressed genes in PBMCs suitable as predictive biomarkers in clinical studies. Some studies that analyzed gene expression profile in peripheral blood cell of CAD or T2D have already been performed [25][26][27]. There were common features and characteristic differences between the current study and previous studies. T2D and CAD both manifest disordered coagulation system, local inflammatory process, or lipid-related disorder. The shared pathophysiology between T2D and CAD may be explained by the common genetic variant, such as CDKN2A/2B, ADIPOR1 [28] and TCF7L2 variants [29]. In our study, many differentially expressed genes from the individual comparisons of T2D, CAD and coexisting condition to control was also found to be overlapped among the three disorders, suggesting shared pathophysiology. Also Spearman's test for gene expression correlation revealed the significant correlation of the three disorders, among which the correlation between T2D and T2D + CAD was the most significant.
The overlapping genes were mainly enriched in GO terms of viral infectious cycle, anti-apoptosis, endocrine pancreas development, innate immune response, and blood coagulation, etc. Those overlapping genes were mainly involved in pathways of Ribosome (RPS3A, UBA52, RSL24D1, RPL7, RPL39, RPL21, RPS24, and RPS12). Ribosomal proteins are involved in cell growth and proliferation, differentiation and apoptosis. Ribosome biogenesis disruption could activate the p53 signaling Table 4 The significantly enriched pathways for the common differentially expressed genes in T2D, CAD, and T2D + CAD   KEGG ID  KEGG term  Count  FDR  Genes   hsa03010  Ribosome  8  0.00  RPS3A,UBA52,RSL24D1,RPL7,RPL39, pathway, resulting in cell cycle arrest and apoptosis. Moreover, it has been correlated to clinical manifestations in pathological states, such as cardiovascular diseases and metabolic disorders. RpL17 inhibited vascular smooth muscle cell growth, and limits carotid intima thickening in mice [30]. Animal study has shown that ribosomal protein S6kinase4 plays a crucial role in pancreatic β-cell function and glucose homeostasis [31]. Ribosomes dysfunction may be applied in the early diagnosis of chronic diseases, such as cancer and cardiovascular diseases.
Besides the common signatures among T2D, CAD and coexisting condition, we also analyzed the gene sets which functioned specially in the development of T2D or CAD. In the T2D-specific PPI networks, the significant hub proteins were TCF4 (Degree = 169), SKP1 (Degree = 164) and UBE2W (Degree = 75). TCF4, also named transcription factor 7-like 2 (TCF7L2), encodes a transcription factor involved in the Wnt signaling pathway. TCF4 stimulates the proliferation of pancreatic β-cells, regulates embryonic development of the Fig. 1 The T2D-specific PPI networks of 64 dysregulated genes. Red nodes indicate up-regulated genes in T2D, blue nodes indicate down-regulated genes in T2D. Pink nodes indicate genes interacting with the differentially expressed genes, and larger icons indicate hub proteins pancreatic mass, and induces the production of the insulinotropic hormone glucagon-like peptide-1 (GLP-1) in intestinal endocrine cells [32]. It plays a critical role in blood glucose homeostasis. Recently, numerous studies have demonstrated an association between TCF7L2 genotype and T2D [33][34][35][36][37]. In a meta-analysis by Wang J et al., the TCF7L2 rs7903146 polymorphism was found to be associated with increased T2D risk in the Chinese population [38]. There were also evidences for a strong interplay between TCF7L2 polymorphisms and CAD [29,39]. In another study on nine hundred subjects referred for cardiac catheterization for CAD diagnosis by Sousa AG et al., a significant association was identified between the TCF7L2 rs7903146 polymorphism and the prevalence and severity of CAD [40].
Interestingly, differentially expressed genes between T2D and normal controls were significantly enriched in Parkinson's disease, suggesting that a shared Fig. 2 The CAD-specific PPI networks of 69 dysregulated genes. Red nodes indicate up-regulated genes in CAD, and blue nodes indicate down-regulated genes in CAD. Pink nodes indicate genes interacting with the differentially expressed genes, and larger icons indicate hub proteins Fig. 3 The verification of mRNA expression of HIF1A, SMAD1, and SKIL in patients with CAD or T2D via GSE23561 pathophysiology of Parkinson's disease and T2D. The link has been confirmed in several epidemiological studies [41,42]. The mechanism behind this association may be that there was a reciprocal regulation between insulin and dopamine [43].
In case of CAD-specific PPI networks, the significant hub proteins were HIF1A, SMAD1 and SKIL. In the provided expression data from GSE23561, we found that SMAD1 and SKIL were specifically up-regulated in CAD with no change in T2D, while HIF1A was both differentially expressed in CAD and T2D, but with opposite trends. In CAD patients, oxygen supply was limited as a result of reduced blood flow brought about by atherosclerotic plaque formation and inflammatory processes taking place within the vascular endothelium [44]. As a consequence of oxidative stress in hypoxic conditions, hypoxia-inducible factor-1 alpha (HIF-1α) was produced to involve in adaptive and repair mechanisms. HIF-1α is a transcriptional factor encoded by the HIF1A gene, functioning in the preservation of oxygen homeostasis. HIF-1α may participate in the occurrence and progression of CAD through activating various genes such as VEGF, HO-1, and ET-1. HIF-1α mRNA was found to be markedly up-regulated in both monocytes and lymphocytes of CAD patients than that of controls, and the expression level of HIF-1α was highly correlated with severity of atherosclerosis and higher level of collateral score [45,46]. In addition, other studies have investigated the correlation between HIF1A polymorphism and CAD. A recent study showed that HIF-1α polymorphisms (rs11549465 and rs11549467) were associated with clinical type and formation of coronary collaterals [47]. The rs2057482 SNP of HIF1A was showed to be associated with increased susceptibility to premature CAD, which may be applied in clinical diagnostics as a susceptibility marker of premature CAD [48].
However, there were some evidences showing that the HIF1A was crucial for first phase insulin secretion and glucose homeostasis. Nagy et al. revealed that polymorphism (g.C45035TSNP, rs11549465) of HIF1A were associated with T1D and T2D in a Caucasian population [49]. Cheng et al. identified that HIF1A could play an important role in β cell function via binding to ARNT promoter in a mouse with β cell-specific Hif1a disruption [50]. Besides, hyperglycaemia appeared to down-regulated HIF-1α mRNA expression in ischaemic myocardium, and inhibited the defective response of HIF-1α to ischaemia [51,52], explaining a positive association between hyperglycaemia at the time of the event and subsequent mortality from myocardial infarction. The dual role of HIF1A connects the pathogenesis of T2D with that of CAD. For this study, we determined the up-regulation of HIF1A in CAD, but didn't detect the expression of down-regulation of HIF1A in T2D in RNA-seq results. The limited sample may explain the inconsistency between RNA-seq results and that of published studies. Further largesample studies are needed to confirm these results.
In our study, SMAD1 and SKIL were both significantly up-regulated in CAD, remained unchangeable in T2D, which were involved in SMAD pathway. SMAD1 polymorphism was reported to be associated with sudden cardiac arrest in CAD patients [53]. In a transgenic mice model with cardiac-specific overexpression of smad1, Masaki M found that transgenic mice had significantly smaller myocardial infarctions and fewer apoptotic deaths of cardiomyocytes after ischemiareperfusion (I/R) injury, suggesting a role of SMAD1 in cardioprotection against I/R injury [54]. SKIL, a component of the SMAD pathway, was showed to be involved in cardiac fibrosis [55,56]. There remained no reports on the association with CAD.

Conclusions
In summary, our data showed that the gene expression profile of T2D, CAD, and coexisting condition were all distinguishable from controls, and displayed common and specific gene expression pattern in each disordered state. To note, viral infectious cycle, antiapoptosis, endocrine pancreas development, innate immune response and blood coagulation were common biological processes among the three conditions. This study provides some evidences in the transcript level to show the association of T2D, CAD and coexisting condition. For this study, the number of sample for RNA-seq was small, which is a limitation of this study, so studies of large sample size need to be conducted to confirm this conclusion.

Additional file
Additional file 1: Table S1. Demographic data and medication of the study population. (XLSX 10 kb)

Abbreviations
BioGRID: Biological general repository for interaction datasets; CAD: Coronary artery disease; GEO: Gene expression omnibus; GO: Gene ontology; KEGG: Kyoto encyclopedia of genes and genomes; PPI: Protein-protein interaction; RPKM: Reads per kilobase of exon per million mapped reads; T2D: type 2 diabetes