Gene expression profiling in whole blood of patients with coronary artery disease

CAD (coronary artery disease) is caused by multiple genetic and environmental factors and the interaction between them. Microarray analysis is a powerful technique for high-throughput global transcriptomic profiling of gene expression. In addition to identifying new target genes, this method also allows the clustering of genes according to similar patterns of expression and/or function from which gene expression fingerprints can be created. Gene expression fingerprints have become a useful tool in research and diagnosis, with multiple applications, for example in the classification of different tumour types [1], in defining the molecular pathogenesis of segmental glomerulosclerosis [2] in and detecting exposure to toxic substances [3].

Microarray studies of human disease are often limited by challenges in obtaining human tissues and by the lack of models that effectively capture clinically relevant disease features. Peripheral blood has become an attractive prime tissue for biomarker detection because of its critical role in immune response, metabolism, communication with cells and the extracellular matrix in almost all tissues and organs in the human body, as well as for the simplicity of sample collection [4–6]. The dynamic and interactive properties of blood give rise to the possibility that subtle changes occurring within the body, such as changes in association with a disease process or in response to an injury, may leave ‘footprints’ in blood. Therefore transcriptional profiling from whole-blood cells might provide an alternative to tissue biopsy in the search for biomarker genes of cardiovascular disease [7].

miRNAs (microRNAs) are small non-coding RNAs that bind mRNAs at their 3′-UTRs (untranslated regions), stimulating mRNA degradation or inhibiting protein translation [8]. Many miRNAs are up-regulated in response to cellular stress [9] and can modify essential cellular functions of proliferation, differentiation and programmed death [10–12]. Evidence indicates that miRNAs are also implicated in cardiovascular disease [13–15] and have the capacity to create cardiac pathology [16–18].

Given that miRNAs respond to acute changes in cell stress, we hypothesized that combining information from whole-blood miRNA profiles and mRNA signatures could improve the knowledge about the pathogenesis of cardiovascular disease. We examined this notion through comprehensive analyses of miRNA and mRNA expression levels in whole blood from patients affected by severe CAD, also attending a cardiac rehabilitation programme, in comparison with healthy subjects.

Patients undergoing elective CABG (coronary artery bypass graft) surgery (n=163) at the Western Infirmary, Glasgow, U.K. were enrolled between October 2006 and August 2008. These patients attended our clinical investigation unit 1 day before surgery. Control subjects (n=63) were recruited by placing an advertisement in the University of Glasgow Newsletter and in local sport centres. Although ethical and funding limitations prevented angiographic investigation of these individuals to absolutely discount the presence of CAD, these subjects presented as healthy at recruitment and were free of any evidence of CAD according to their patient history. Related studies [19–21] from our group have shown that our healthy control subjects present with normal ECG and no history of arterial hypertension, angina, CAD or peripheral artery disease. Furthermore, we have shown that control subjects have better aortic compliance and lower levels of vascular ROS (reactive oxygen species).

At 4–6 weeks after the CABG procedure, 37 out of the 163 patients were also included in a cardiac rehabilitation programme, consisting of 60 min exercise (15 min warming up, 30 min cardiovascular training and 15 min warming down), twice a week for 10 weeks. The initial 4–6-week period following surgery was to allow the patients to recover from the procedure to a degree where they were well enough to commence exercise, and also would allow some of the more acute responses to surgery (e.g. increases in inflammatory markers) to subside. Although some markers remained elevated [e.g. CRP (C-reactive protein)], this difference was not found to be significant with regard to control subjects. The exercise regimen was adapted to each individual's exercise capacity, assessed by validated step and walking tests [19,22,23] at inclusion and after completion of the rehabilitation programme. For microarray analysis, we selected 12 subjects with severe CAD and 12 control subjects; and ten patients attending the exercise programme examined before and after its completion. Details of these patients are given in Supplementary Tables S1 and S2 (at http://www.clinsci.org/cs/119/cs1190335add.htm).

The study complies with the Declaration of Helsinki and was approved by the local ethics committee. All participants gave written informed consent and used a standardized questionnaire to self-report the health status and medication taken.

Venous blood samples were drawn via antecubital venipuncture from each subject in a sitting position after fasting for 3 h, and routine biochemistry was assessed in plasma. An aliquot of whole blood (2.5 ml) for each subject was collected directly into PAXgene blood RNA tubes (PreAnalytiX) and stored at −20 °C within 24 h before RNA extraction. After collection of all samples, total RNA was isolated using the PAXgene blood RNA kit (PreAnalytiX), according to the manufacturer's protocol, and subjected to on-column DNase I treatment with RNase-free Dnase (Qiagen). RNA quantity and quality were determined using a NanoDrop Technologies ND-1000 spectrophotometer and an Agilent 2100 bioanalyser. A total of 400 ng of RNA (with an RNA integrity number >7.5) was used to generate cDNA (Illumina TotalPrep RNA Amplification kit; Ambion). Reverse transcription with T7 oligo(dT) primers was used to produce first-strand cDNA. The cDNA then underwent second-strand synthesis and RNA degradation by DNA polymerase and RNase H, followed by clean up. In vitro transcription technology, along with biotin UTP, was employed to generate multiple copies of biotinylated cRNA. The labelled cRNA was purified using a filter cartridge and quantified using a NanoDrop Technologies ND-1000 spectrophotometer. The integrity of cRNA was evaluated using an Agilent 2100 bioanalyser.

Small RNA species (including miRNA) were also extracted by modifying the PAXgene procedure by harvesting the flow-thorough of the PAXgene RNA column, and isolated using the miRNeasy Mini Kit (Qiagen), according to the manufacturers' protocols. Successful isolation of small RNAs was confirmed using Agilent 2100 bioanalyser analysis.

All mRNA and miRNA expression studies were performed on samples from individual patients (not pooled).

Microarray analysis of gene expression was performed on the Illumina Beadstation platform, using Illumina Sentrix humanref-6 beadchips. Labelled cRNA (1.5 μg) was used for hybridization to the array, according to the manufacturer's protocol. A maximum of 10 μl of cRNA was mixed with 20 μl of GEX-HYB hybridization solution. The pre-heated 30 μl assay sample was dispensed on to each array and incubated for 18 h at 58 °C. Following hybridization, samples were washed and scanned with a BeadArray Reader (Illumina). All data passed quality control analysis as assessed by the Illumina on-board software (BeadStudio 3.4) and Principle Component Analysis (Partek Genomics Suite). Quantile normalization of gene expression was performed (BeadStudio 3.4), and differential expression was determined by Rank Products analysis [24], where a 5% FDR (false discovery rate) cut-off value was used. Rank Products is a simulation method based on the average rank each gene has for a pairwise comparison difference compared with what would be expected under the null hypothesis of no differential expression. It has been specifically developed for the analysis of microarray experiments and was shown to be more powerful than Student's t tests [25]. In addition, in a comparison of 11 different methods for microarray analysis [26], Rank Products was shown to have performed well, in particular with small sample sizes. Significance was assessed using the FDR multiple testing correction method [27] with an FDR cut-off of 10% (gene expression) or 5% (miRNA expression). For probe sets with an FDR >10%, a further 1000 permutations were carried out to obtain more precise P values and FDR values. Molecular interactions between genes were mapped to common canonical pathways using the Pathway Explorer function within Ingenuity Pathway Analysis software (Ingenuity Systems; www.ingenuity.com).

miRNA microarray analysis was performed to determine differential expression in blood-borne miRNAs between (i) healthy individuals and patients with CAD (n=7), and (ii) patients with CAD before a cardiac rehabilitation exercise programme and at 6 weeks after exercise (n=8). All miRNA expression profiling was carried out by LC Sciences using a proprietary microfluidic array (Chip ID miRHuman_11.0_080411) based on the Sanger miRBase v11.0 database (http://www.sanger.ac.uk/Software/Rfam/mirna/), designed to detect 887 human miRNAs. Proprietary analysis employed a Student's t test to compare data. An α level of 10% was accepted for the Student's t test analysis in order to identify sufficient targets for further evalulation. miRNA data were subsequently re-analysed in-house using the Ranked Products method with an FDR value of 5% [24]. miRNAs achieving significance under both analyses, and which had a signal intensity >500 fluorescence units, were considered for further analysis.

Bioinformatic determination of downstream predicted targets for candidate miRNAs was performed as described previously by Selbach et al. [28]. Briefly, an unbiased search for sequence motifs (representing positions 2–7 of the miRNA seed region) in a library of genome-wide 3′-UTRs in human was undertaken for each miRNA. On this basis genes were assigned to either target gene or non-target groups. Gene expression was then compared between patients with CAD and controls, and patients with CAD pre- and post-exercise respectively, for gene targets and non-targets, and compared with unpaired (CAD compared with controls) or paired (post- compared with pre-exercise) Student's t tests (P<0.05).

qRT-PCR was performed to validate changes in selected genes by IPA software. Purified RNA was converted into cDNA using TaqMan® reverse transcription reagents in a total reaction volume of 20 μl (Applied Biosystems) containing: 10×TaqMan® reverse transcription buffer (2.0 μl), 25 nM MgCl₂ (4.4 μl), deoxyNTPs (4.0 μl), random hexamers (1.0 μl), RNase inhibitor (0.4 μl) and MultiScribe reverse transcriptase (0.5 μl), using the following cycling parameters: 25 °C for 10 min, 48 °C for 30 min and then 95 °C for 5 min. Appropriate ‘no template’ qRT-PCR controls were processed similarly. cDNA was diluted to a final volume of 80 μl by addition of nuclease-free water and stored at −20 °C until use. qRT-PCR was performed in 5 μl volumes, in a reaction buffer containing 2.5 μl of 2×TaqMan® Universal PCR Master Mix, 2 μl of cDNA template and 0.25 μl of each TaqMan® probe/primer mixture (Supplementary Table S3 at http://www.clinsci.org/cs/119/cs1190335add.htm). Amplification reactions were performed in triplicate, duplexed with GAPDH (glyceraldehyde-3-phosphate dehydrogenase) as an endogenous control. Assay plates were sealed, vortex-mixed, centrifuged (1 min at 850 g), and subjected to PCR (ABI Sequence Detection System 7900HT), using ABI-suggested cycling parameters: 50 °C for 2 min, 95 °C for 10 min and then 40 cycles of 95 °C for 15 s, followed by 60 °C for 1 min. Data were exported using a 0.2 threshold, with 3–15 cycles set as baseline. Cycle threshold (C_t) values were normalized by the geometric means of GAPDH. Triplicates were averaged to calculate an expression value for each sample.

Using the power calculation method described by Pounds and Cheng [29], a desired FDR of 0.05, a desired average power of 80% and assuming that 30 out of 833 miRNAs to be investigated have a θ value of 1.8 (with the remainder being assumed not to be differentially expressed), we calculated that a a sample size of 12 per group was needed.

Statistical analyses were performed using GraphPad Prism software. Values are expressed as means±S.D. Paired or unpaired Student's t tests were used to compare data as appropriate. A P value <0.05 (two-sided) was considered significant. For analysis of microarray expression data, a Ranked Products approach was used, with FDRs of either 10% or 5% as appropriate.

Between patients with CAD (n=12) and control subjects (n=12), we identified 365 differentially expressed genes: 175 genes were up-regulated and 190 genes were down-regulated. In patients with CAD studied across an exercise rehabilitation programme (n=10), we found 645 genes differentially expressed between the beginning and completion of the programme, with 196 genes up-regulated and 449 genes down-regulated.

Biological pathway analysis identified a number of canonical pathways as being significantly and consistently modulated across patients with CAD compared with control subjects (Figure 1a), and across the cardiac rehabilitation exercise regime (Figure 1b). From (i) the overall degree of modulation, (ii) the proportion of differentially expressed genes in each pathway, and (iii) biological relevance, further analyses were focused on mitochondrial dysfunction and oxidative phosphorylation pathways (also see Supplementary Table S4 at http://www.clinsci.org/cs/119/cs1190335add.htm).

Within these pathways, expression of NDUFB3 [NADH dehydrogenase (ubiquinone) 1β subcomplex 3], UQCRQ (ubiquinol-cytochrome c reductase, complex III subunit VII), COX7C (cytochrome c oxidase subunit VIIc), and ATP5I (ATP synthase, H⁺ transporting, mitochondrial F0 complex, subunit E) were greater in patients with CAD relative to control subjects. Completion of the rehabilitation period was associated with a down-regulation of the genes in these pathways: NDUFA1 [NADH dehydrogenase (ubiquinone) 1α subcomplex 1], CASP3 (caspase 3), COX7C, ATP5I and ATP5L (ATP synthase, H⁺ transporting, mitochondrial F0 complex, subunit G) (see Supplementary Figures S1 and S2 at http://www.clinsci.org/cs/119/cs1190335add.htm). qRT-PCR confirmed the up-regulation of NDUFB3 and COX7C in patients with CAD and the down-regulation of CASP3 and ATP5L following the rehabilitation programme (Figure 2).

The miRNAs hsa-miR-140-3p (median signal intensity, 738 compared with 383; P=0.017) and hsa-miR-182 (median signal intensity, 836 compared with 549; P=0.093) were differentially regulated between patients with CAD and control subjects (Figure 3a). Accordingly, global analysis of predicted targets within the gene expression dataset demonstrated significantly reduced expression of genes with target regions of hsa-miR-140-3p and hsa-miR-182 (Figures 3b and 3c). Within the post-surgery exercise groups, hsa-miR-92a (median signal intensity, 4712 compared with 6250 respectively; P<0.01) and hsa-miR-92b (median signal intensity, 1119 compared with 2163 respectively; P<0.01), which have identical downstream mature sequences, were differentially regulated between the beginning and completion of the rehabilitation programme (Figure 4a). Expression levels of predicted gene targets of hsa-miR-92a and hsa-miR-92b were significantly reduced (Figure 4b). Moreover, predicted targets of hsa-miR-92a and hsa-miR-92b include the genes NDUFA1 and CASP3. The observed reduced expression of these genes following the rehabilitation programme is consistent with the increased expression of hsa-miR-92a and hsa-miR-92b (Figure 2).

In the prsent study, we have successfully used gene expression profiling in whole blood to establish genes that are differentially regulated between patients with CAD and control subjects, and whose expression is modified following a cardiac rehabilitation programme in patients with CAD. We have, for the first time, identified differentially regulated miRNAs in whole blood in patients with CAD compared with healthy control subjects. Among the biological pathways containing the differentially expressed genes, we have identified those involved in mitochondrial dysfunction and oxidative phosphorylation as being the most consistently affected and modulated.

Studies on gene expression in humans are particularly difficult due to the limited availability of tissue from well-phenotyped patients, and even more so from healthy subjects. The use of blood as a ‘surrogate tissue’ that can be obtained with a minimally invasive procedure is therefore an attractive alternative. Although the majority of studies have been performed in mononuclear cells isolated from peripheral blood [30], there have been reports on gene expression profiling in whole blood in a number of human diseases, including Huntington's disease [4] and cancer [5,6]. Many different techniques have been proposed to handle peripheral blood samples prior to RNA isolation based on experimental design [31]. We utilized the PAXgene Blood RNA system as it provides a way to stabilize RNA immediately after sample collection and facilitates storage of the samples for a relatively long period of time without compromising RNA integrity [31–34]. The straightforward procedure of sample collection is an important condition for future use of whole-blood gene expression profiling in clinical practice. Most recently, whole-blood gene expression profiling has also been applied to cardiovascular disease. Sinnaeve et al. [35] showed a relationship between gene expression profiles in whole blood with a CAD index as a marker of the extent of disease. Our present study adds important new features to this area of research.

In face of the multiplicity of differentially expressed genes both between patients with CAD and controls, and between pre- and post-cardiac rehabilitation, we have used pathway analysis to examine patterns of expression rather than individual genes. Using this approach, we were able to interpret data on relatively small fold changes in gene expression, which may on their own be of limited importance but, in combination, represent multiple hits to specific pathways and may thereby affect, for example, mitochondrial function. The finding that we were unable to confirm all differences in gene expression from the microarray experiment with qRT-PCR may again be in keeping with the existence of multiple changes, but of relatively small magnitude, below the 2-fold limitation of qRT-PCR. This pattern of numerous small changes of expression could be reasonably expected in complex cardiovascular traits such as CAD. Of particular note is that the direction of these and other multiple ‘hits’ to mitochondrial pathways, although small in magnitude, reverse following exercise (Figure 2), indicating that this intervention may offer specific benefits in relation to reducing oxidative stress in individuals recovering from coronary surgery.

We have further corroborated our results by simultaneously examining miRNA expression. Again, interpretation of results, especially with regard to targets of particular miRNAs, is challenging [25], but the results support our main findings. The results on the integration of miRNAs with expression of their predicted targets does indeed show a significant generalized down-regulation of these targets. This is in keeping with a functional role of these molecules, although such a role of miRNA in blood remains to be elucidated. We have, however, shown that expression of hsa-miR-92a and hsa-miR-92b are increased following cardiac rehabilitation, which occurs in parallel with a confirmed reduction in expression of two of their predicted target genes. Hsa-miR-92 has been found to be modulated in the process of neointima formation in a rat model of vascular injury [36], but to the best of our knowledge has not yet been clearly implicated in the pathophysiology of human cardiovascular disease. Although acknowledging that hsa-miR-92a and hsa-miR-92b are likely to modify a plethora of potential targets, and that the specific interactions with NDUFA1 and CASP3 require further detailed examination, it is nonetheless intriguing to identify physiological candidates which are down-regulated in concert with increased expression of these miRNA. A further level of complexity of the regulation of miRNA expression and activity can be seen by the generation of the identical mature hsa-miR-92 miRNA from two different genetic loci. In relation to target prediction and validation, the present study has focused solely on that common mature sequence. However, differential expression of the two pri-miRNAs (potentially under the control of different up-stream regulators) may occur in a heterogeneous manner in different tissues, cells or compartments, offering even more exquisite control of their relative activities. The dynamic change of miRNA and gene expression towards a ‘healthier’ profile following cardiac rehabilitation is a unique feature of our present study, supporting further the relevance of the identified genes and miRNAs.

We acknowledge that our present findings only provide initial evidence for a role of the discovered genes and miRNAs. However, these findings are particularly attractive as targets for intervention. Oxidative stress has been implicated by a plethora of studies as playing a key role in a range of cardiovascular diseases [37–40]. Mitochondrial dysfunction, and particularly generation of ROS by mitochondria, is an important yet relatively novel mechanism contributing to the development of oxidative stress and cardiovascular disease [41,42]. We have recently shown that a specific mitochondria-targeted antioxidant, MitoQ₁₀, reduces the increased vascular superoxide generation in stroke-prone spontaneously hypertensive rats and leads to a reduction in BP (blood pressure) [43]. Our present results indicate that there may be further options for mitochondria-targeted therapies in human cardiovascular diseases including CAD.

With regard to limitations of the present study, we cannot exclude the possibility that the observations obtained in the cardiac rehabilitation group are not, at least to some extent, due to, for example, healing of the sternotomy. Nonetheless, significant changes were observed in both gene expression and miRNA profiling in the blood of these patients and, in particular, specific pathways were identified as being modified which were also shown to be changed in the CAD and control comparison. Moreover, in this latter comparison, these differences were determined in the patients with CAD prior to undergoing surgery, illustrating that dynamic changes in gene and miRNA expression can occur without the provocation caused by, for example, significant trauma such as surgery.

Our present study is also limited by the relatively small sample sizes due to the significant costs involved in microarray studies. However, the careful selection of patients and well-matched control subjects from a larger cohort of well-characterized individuals reduced the ‘noise’ in our expression profiling experiments. We have further addressed this issue by studying a second set of patients before and after completing the rehabilitation programme using paired statistical analysis. We are confident that, due to integration of these two datasets and further analysis, together with miRNA expression results, the identified genes and pathways are not only valid, but may also be clinically relevant. It should also be noted that the sample size in our study is equal to or greater than that in other human global gene expression profiling studies [44–46]. Clearly, the aim of the present initial study was to identify such targets. Subsequent validation of these targets in larger studies will be a necessary step to confirm their candidacy as biomarkers and/or functionally relevant molecules in blood.

In summary, mRNA and miRNA expression profiling in whole blood is a promising tool to discover genes that determine cardiovascular phenotypes. The identified genes and pathways are likely to play a key role in the pathogenesis and progression of chronic CAD. Moreover, we have identified biological pathways linked to vascular function which appear to be modifiable by an exercise intervention programme. Further elucidation of these mechanisms will contribute to our understanding of disease processes and ultimately lead to new therapeutic and preventative strategies.

This work was supported by a BHF (British Heart Foundation) Chair (to A.F.D.), a BHF Programme Grant [grant number RG/07/005/23633 (to A.F.D. and C.D.)], the European Union's Sixth Framework Programme InGenious HyperCare [grant number LSHM-CT-2006-037093 (to A.F.D. and C.D.)], and the Wellcome Trust Cardiovascular Functional Genomics Initiative [grant number 066780/2/012 (to A.F.D.)].

ATP5I

ATP synthase, H⁺ transporting, mitochondrial F0 complex, subunit E

ATP5L

ATP synthase, H⁺ transporting, mitochondrial F0 complex, subunit G

CABG

coronary artery bypass graft

CAD

coronary artery disease

CASP3

caspase 3

COX7C

cytochrome c oxidase subunit VIIc, CRP, C-reactive protein

FDR

False Discovery Rate

GAPDH

glyceraldehyde-3-phosphate dehydrogenase

miRNA

microRNA

NDUFA1

NADH dehydrogenase (ubiquinone) 1α subcomplex 1

NDUFB3

NADH dehydrogenase (ubiquinone) 1β subcomplex 3

qRT-PCR

quantitative real-time PCR

ROS

reactive oxygen species

UQCRQ

ubiquinol-cytochrome c reductase, complex III subunit VII

UTR

untranslated region

Chiara Taurino recruited the patients and performed the gene expression studies; William Miller performed the gene expression studies, analysed and interpreted the data, and drafted the manuscript; Martin McBride, John McClure and Raya Khanin analysed and interpreted the gene expression data; María Moreno performed the gene expression studies; Jane Dymott designed the study and the recruited patients; Christian Delles designed the study, interpreted the data and wrote the manuscript; and Anna Dominiczak designed the study and critically reviewed the manuscript.