Whole blood transcriptomic profiling identifies molecular pathways related to cardiovascular mortality in heart failure.

AIMS: Chronic heart failure (CHF) is a systemic syndrome with a poor prognosis and a need for novel therapies. We investigated whether whole blood transcriptomic profiling can provide new mechanistic insights into cardiovascular (CV) mortality in CHF. METHODS AND
RESULTS: Transcriptome profiles were generated at baseline from 944 CHF patients from the BIOSTAT-CHF study, of whom 626 survived and 318 died from a CV cause during a follow-up of 21 months. Multivariable analysis, including adjustment for cell count, identified 1153 genes (6.5%) that were differentially expressed between those that survived or died and strongly related to a validated clinical risk score for adverse prognosis. The differentially expressed genes mainly belonged to five non-redundant pathways: adaptive immune response, proteasome-mediated ubiquitin-dependent protein catabolic process, T-cell co-stimulation, positive regulation of T-cell proliferation, and erythrocyte development. These five pathways were selectively related (RV coefficients >0.20) with seven circulating protein biomarkers of CV mortality (fibroblast growth factor 23, soluble ST2, adrenomedullin, hepcidin, pentraxin-3, WAP 4-disulfide core domain 2, and interleukin-6) revealing an intricate relationship between immune and iron homeostasis. The pattern of survival-associated gene expression matched with 29 perturbagen-induced transcriptome signatures in the iLINCS drug-repurposing database, identifying drugs, approved for other clinical indications, that were able to reverse in vitro the molecular changes associated with adverse prognosis in CHF.
CONCLUSION: Systematic modelling of the whole blood protein-coding transcriptome defined molecular pathways that provide a link between clinical risk factors and adverse CV prognosis in CHF, identifying both established and new potential therapeutic targets.

Introduction

The prognosis of patients with chronic heart failure (CHF) remains poor, despite substantial improvements in CHF diagnosis and treatment. Several prognostic risk scores have been developed for CHF, , yet the precise mechanistic links between such clinical risk scores and survival remain unclear. Furthermore, beyond natriuretic peptides, several circulating protein biomarkers that have been shown to be independently associated with prognosis in CHF , , remain to have the underpinning molecular nature of those relationships elucidated.

Whole blood transcriptomic profiling is a practical and powerful methodology to study genome‐wide biology in large numbers of patients. Despite CHF being a systemic syndrome, no study has investigated whether whole blood transcriptomic profiles can inform about the molecular processes underpinning the observed relationships between clinical or protein biomarkers and cardiovascular (CV) mortality in these patients. If the whole blood transcriptome can be shown to reflect important drivers of mortality, then modelling whole blood global gene expression may help identify clinically relevant pathways to target therapeutics development.

In this study, we report the first large‐scale whole blood transcriptomic profiles from patients with CHF. Our main aims were to (i) contrast such profiles between patients who survived or died, (ii) examine whether this provides new insights into molecular mechanisms through which clinical and protein biomarkers associated with adverse outcome might affect prognosis, and (iii) explore the utility of blood‐based transcriptomic signatures as a novel tool to aid in vitro drug‐repurposing analysis in CHF.

Methods

Participants

This was a nested case–control study undertaken within the BIOSTAT‐CHF study. Briefly, BIOSTAT‐CHF enrolled an index cohort of 2516 patients from 69 hospital centres in 11 European countries between 2010 and 2014. Patients were aged >18 years with symptoms of new‐onset or worsening CHF, confirmed either by a left ventricular ejection fraction (LVEF) of ≤40% or B‐type natriuretic peptide (BNP) level >400 pg/ml and/or N‐terminal pro‐B‐type natriuretic peptide (NT‐proBNP) plasma levels >2000 pg/ml, treated with either oral or intravenous furosemide ≥40 mg/day or equivalent at the time of inclusion. After enrolment and obtaining a blood sample, all participants had their heart failure medication optimized to guideline‐recommended or maximally tolerated doses. The median follow‐up of the cohort was 21 months, with an interquartile range of 15 months. The study complied with the 2008 Declaration of Helsinki and was approved by the relevant ethics committee in each centre. All participants gave written informed consent to participate.

For the current transcriptome analysis, participants were categorized into two groups: survivors (those participants who survived during the follow‐up period without hospitalization), and non‐survivors (those participants who died from a CV cause). Cause of death (CV or non‐CV) was adjudicated by a senior clinician (A.A.V.) based on reports of mortality events using previously reported criteria. There were 441 CV deaths (out of a total of 657 deaths), and 1437 participants survived without hospitalization. After the exclusion of participants with no usable RNA samples, 332 who died and 1061 who survived, were available to study. The two groups, those that survived or died, were matched for age and sex with an approximate ratio of 1:2, to match the number of arrays available (online supplementary Figure  ).

Generation of transcriptomic profiles

Whole blood was collected into PAXgene tubes (Qiagen) and RNA was isolated using TRizol® (Life Technologies), dissolved in 20 µl RNAse‐free water and profiled using the Affymetrix Human Transcriptomic Array 2.0 (HTA, Thermo Fisher Scientific) which contains ∼7 million probes to provide a comprehensive view of the transcriptome. Samples were processed following the manufacturer's protocol (Sutton Bonington Array facility, University of Nottingham). Briefly, 250 ng of RNA was processed to single‐stranded sense fragmented DNA using the WT PLUS Reagent Kit. Five µg of fragmented, end‐labelled sense stranded target cDNA was hybridized to each array, and scanned in batches of 24, using a 30007G scanner. The running order used a random list generator to ensure even distribution of cases and controls during profiling. Full details of the processing of the arrays and quality checks of the gene expression data are provided in online supplementary Appendix  . A total of 944 samples were successfully profiled, of whom 626 (66%) were survivors and 318 (34%) died of a CV cause.

Statistical and bioinformatics analysis

RNA expression values were derived for 17 748 protein‐coding genes (ENSG), then normalized and logarithm (base 2) transformed. We applied a multivariable linear model analysis with group (survival and non‐survival), age, sex and estimates of major white cell populations as covariates using the limma package in R, to identify differentially expressed ENSGs. Estimated white cell proportions were derived by applying ABIS, with absolute values being scaled to 100 prior to further analysis (online supplementary Appendix  ). P‐values were adjusted using the Benjamini–Hochberg false discovery rate (FDR) method; an FDR value of <5% was used. Gene ontology (GO) profiles (Biological Processes) were assessed using Metascape, with the detectable transcriptome used as the GO background to calculate enriched pathway statistics. The main drivers to GO enrichment were identified by removing associated processes, using the topGO hierarchical analysis package in R.

To assess whether there was a relationship between the observed transcriptomic differences between survivors and non‐survivors, and clinical and biochemical risk factors associated with CHF prognosis, we used a reproducible prognostic risk score for mortality developed in BIOSTAT‐CHF (BIOSTAT‐CHF risk score, BRS). The BRS model considered 42 clinical and biochemical variables, identifying five variables (more advanced age, higher blood urea nitrogen and NT‐proBNP, lower haemoglobin, and failure to prescribe a beta‐blocker) as the strongest independent predictor of mortality with a combined C‐statistic of 0.73. We undertook linear correlation analysis of gene expression with the BRS and compared the overlap of this set of genes with those differentially expressed between survivors and non‐survivors. The effect of adjusting for the BRS on the latter set of genes was also investigated.

To gain insight into potential protein drivers of the transcriptional changes, we correlated the transcriptomic changes with plasma levels of >150 circulating proteins available in the BIOSTAT‐CHF cohort, some of which have been reported to be associated with prognosis in the cohort. , The majority of proteins were measured using two commercial platforms (Olink Proseek Multiplex and Alere Luminex Panel), with a minority derived from routinely available hospital clinical assays or bespoke assays. Specifically, fibroblast growth factor 23 (FGF23) data were generated using the c‐terminal ELISA (Immutopics, Inc., San Clemente, CA, USA) while hepcidin (HEPC) was measured using the competitive ELISA developed by Kroot et al. To examine the relationship between each of the non‐redundant GO categories, and the protein biomarkers, we estimated the RV coefficient, a multivariate generalization of the squared Pearson correlation coefficient. Pearson correlation coefficients, adjusting for multiple testing (Bonferroni correction), were also calculated for the significant associated proteins and individual genes within the top‐ranked GO categories. For proteins most associated (RV coefficient >0.20) with at least one significant GO, we examined the relationship with CV mortality, using a Cox proportional hazards model.

Finally, to explore the potential of our approach to identify new drug treatments, we utilized the CMap‐L1000v1 in vitro drug signature database (https://clue.io/) as previously described. Briefly, to represent the prognosis‐related transcriptomics signature, we used the 60 most upregulated and 60 most down‐regulated genes that were also present in the CMap‐L1000v1 database. Choice of 60/60 genes is empirical and informed by our recent work on insulin resistance drug repurposing that yielded a drug list composed of ∼50% true‐positive drug classes and/or drug targets. The matching score was calculated by aggregating the transcriptional pattern across nine independent cell lines for >2500 compounds, many of which include Food and Drug Administration approved drugs. For each matching drug, the known protein targets were identified using PubChem and the small molecule suite. A network of significant pathway interactions was derived using metascape.org, using the CHF survival associated genes and the top ranked known protein target of each drug as input values, and the detectable transcriptome as background.

Results

Subjects

The two groups (survivors and non‐survivors) were, by design, closely matched for age and sex (Table  and online supplementary Figure  ). Patients were predominantly Caucasian; there was no difference in hypertension but there was a higher prevalence of both type 2 diabetes mellitus (T2DM) and CHF of ischaemic aetiology, in the non‐survivor group. Baseline use of CHF medication was similar in both groups, except for a greater use of beta‐blockers in the survivor group. The echo‐estimated LVEF was similar between groups (≈31%) as was heart rate. However, patients that died had worse symptoms, lower systolic and diastolic blood pressure and higher plasma concentrations of natriuretic peptides, and lower estimated glomerular filtration rate (eGFR), haemoglobin and body mass index (BMI), at baseline.

Table 1

Demographic, clinical and laboratory variables for study groups

Variable	Survivor (n = 626)	Non‐survivor (n = 318)	p‐value
Demographics
Male sex	75.7 (474)	74.5 (237)	0.748
Age (years)	71 (10.7)	71.4 (11)	0.584
BMI (kg/m2)	27.7 (6.2)	26.6 (7.3)	0.206
Clinical profile
NYHA class			<0.001
I	2.9 (18)	0.9 (3)
II	44.2 (277)	22.6 (72)
III	41.2 (258)	58.5 (186)
IV	8.6 (54)	14.2 (45)
LVEF (%)	31.6 (9.9)	31.3 (12.4)	0.747
Heart rate (bpm)	79.2 (20.2)	80.4 (19.4)	0.386
Systolic blood pressure (mmHg)	128.0 (21.7)	121.1 (22.3)	<0.001
Diastolic blood pressure (mmHg)	76.6 (13.7)	71.7 (12.3)	<0.001
HF history
Ischaemic aetiology	61.2 (340)	70.3 (206)	0.01
HF hospitalization in previous year	22.7 (142)	42.8 (136)	<0.001
Medical history
Hypertension	65.2 (408)	63.2 (201)	0.599
Diabetes mellitus	28.3 (177)	39.0 (124)	0.001
Medication at baseline
ACE inhibitors or ARB	74.8 (468)	69.2 (220)	0.081
Beta‐blockers	84 (526)	77 (245)	0.011
Mineralocorticoid receptor antagonist	51.3 (321)	50.9 (162)	0.977
Laboratory measurements
Haemoglobin (g/dl)	13.4 (1.8)	12.6 (1.8)	<0.001
Erythrocytes (million cells/µl)	4.5 (0.6)	4.4 (0.9)	0.31
Leucocytes (×109/L)	7.8 [6.6, 9.3]	7.8 [6.4, 9.6]	0.908
eGFR (ml/min/1.73 m2; CKD‐EPI)	62.0 [46.9, 76.7]	52.3 [34.8, 69.4]	<0.001
BNP (pg/ml)	174.7 [78.3, 371.1]	338.1 [176.2, 652]	<0.001
NT‐proBNP (ng/L)	2005 [949, 4642]	4121 [2332, 9275]	<0.001

Values are expressed as % (n), mean (standard deviation), or median [interquartile range].

ACE, angiotensin‐converting enzyme; ARB, angiotensin II receptor blocker; BMI, body mass index; BNP, B‐type natriuretic peptide; CKD‐EPI, Chronic Kidney Disease Epidemiology Collaboration equation; eGFR, estimated glomerular filtration rate; HF, heart failure; LVEF, left ventricular ejection fraction; NT‐proBNP, N‐terminal pro‐B‐type natriuretic peptide; NYHA, New York Heart Association.

Whole blood transcriptomic differences between survivors and non‐survivors

Laboratory‐measured total leucocyte count did not differ between those that survived or died (Table  ). However, counts modelled for individual blood cell types using specific gene expression markers (see Methods) identified significant reduction in T‐cells (−17%, p < 1 × 10−10) and B‐cells (−11%, p < 1 × 10−3) in the non‐survivor group, and an increase in neutrophils (+5%, p < 1 × 10−5) and basophils (+14%, p < 1 × 10−3; online supplementary Figure  ). After including the estimated cell subtypes as covariates in the analysis together with age and sex, we found that 1153 genes (6.5%) were differentially expressed (FDR <5% and >5% difference, 557 up‐regulated and 596 down‐regulated) between survivors and non‐survivors (Figure  and online supplementary Table  ). Up‐ and down‐regulated genes often belonged to the same biological pathways (online supplementary Figure  ) and mapped to several interrelated biological pathways (Figure  and online supplementary Table  ). The processes most responsible for the GO results were identified by removing the redundant topological features of the GO structure (online supplementary Figure ). This identified five main ‘driver’ GO categories, each with a robust fold enrichment (FE) and low FDR: adaptive immune response (FE = 2.3, FDR <0.1%), proteasome‐mediated ubiquitin‐dependent protein catabolic process (UPP; FE = 1.7, FDR <3%), T‐cell co‐stimulation (FE = 3.9, FDR <0.1%), positive regulation of T‐cell proliferation (FE = 3.0, FDR <0.1%), and erythrocyte development (FE = 4.3, FDR <0.1%).

Figure 1

Differential gene expression and biological process gene ontology analysis. (A) A volcano plot of the differential gene expression pattern between groups, calculated using limma and the following model: ENSG ∼ Age + Sex + log(sum T‐cells) + log(sum B‐cells) + log(neutrophils) + log(basophils) + group. (B) Global gene ontology (biological processes) analysis carried out using Metascape and the list of 1153 differentially regulated genes. The 17 748 detected protein‐coding genes were used as the background and the pre‐analysis settings were: threshold >2.0 and non‐adjusted p‐values threshold p < 0.0001. The plot x‐axis is p‐values (log base 10) while the displayed categories had a false discovery rate of 1.4% or better (see online supplementary Table ).

The differentially expressed genes in each of these pathways are shown in Figure  . Overall, the pattern of changes in gene expression indicates an impairment of the adaptive immune response, down‐regulation of T‐cell co‐stimulation and disruption of positive regulation of T‐cell proliferation in non‐survivors. In contrast, the UPP pathway and the erythrocyte development pathway were both activated in non‐survivors (see Discussion).

Figure 2

Individual gene expression responses within the topology adjusted top‐ranked gene ontology categories. Differential gene expression is presented across the five topology adjusted gene ontology categories identified using weighted Fisher's exact test, fold enrichment calculations and BH correction of Fisher's test statistics using the R package topGO. FC, fold change.

BIOSTAT‐CHF prognostic risk score and transcriptomic changes

Linear regression analysis of gene expression versus the BRS identified 1086 correlated genes, and these predominantly belonged to the same pathways as the genes differentially expressed between survivors and non‐survivors (online supplementary Figure  ). On an individual gene level basis, 390 genes were common to both lists. Inclusion of the BRS as a covariate in the analysis of the expression differences between survivors and non‐survivors removed all but one of the gene expression differences between groups. This suggests that the observed variation in the blood transcriptome is strongly related to clinical factors associated with prognosis.

Transcriptomic changes and circulating protein biomarkers

To examine whether circulating proteins represented potential drivers of the transcriptomic changes, and therefore potential intermediaries between the clinical prognostic variables and the observed gene expression differences, we modelled the relationships between >150 circulating proteins, and the five top‐ranked GO categories. While most demonstrated limited association (online supplementary Table  ), seven proteins demonstrated an RV coefficient of 0.2 or greater, with at least one of the five GOs (Table  ). These were, in order of association: FGF23, soluble ST2 receptor (sST2), adrenomedullin (ADM), HEPC, pentraxin‐3 (PTX3), WAP 4‐disulfide core domain 2 (WFDC2) and interleukin‐6 (IL‐6). Notably, six of these proteins positively covaried with each other, while HEPC demonstrated a negative correlation, especially with FGF23 (online supplementary Figure  ). These seven proteins also demonstrated a significant univariate association with mortality in the present nested case–control sample (Table  ), as previously reported for some of these proteins in the full BIOSTAT‐CHF cohort. , Notably, although T2DM was more common in non‐survivors at baseline, we did not observe any pathway‐level association with glycated haemoglobin (online supplementary Table  ).

Table 2

Top ontologies and RV coefficient estimates with circulating protein biomarkers

GO pathway (no. of genes)	FGF23	sST2	ADM	HEPC	PTX3	WFDC2	IL‐6
Adaptive immune response (56)	0.20	0.25	0.20	0.04	0.22	0.21	0.20
T‐cell co‐stimulation (14)	0.16	0.23	0.21	0.01	0.22	0.20	0.18
Positive regulation of T‐cell proliferation (18)	0.25	0.25	0.22	0.08	0.21	0.19	0.18
Erythrocyte development (11)	0.17	0.09	0.13	0.18	0.04	0.04	0.04
Proteasome‐mediated ubiquitin‐dependent protein catabolic process (47)	0.30	0.18	0.21	0.22	0.10	0.09	0.09
Univariate association with survival (HR)	1.64	1.70	1.92	0.84	1.57	1.88	1.43
	p < 1 × 10−34	p < 1 × 10−20	p < 1 × 10−22	p < 1 × 10−4	p < 1 × 10−9	p < 1 × 10−28	p < 1 × 10−12

They are listed in rank order from the RV analysis as follows: ADM, adrenomedullin; FGF23, fibroblast growth factor 23; HEPC, hepcidin; IL‐6, interleukin 6; sST2, soluble ST2 receptor (decoy receptor for IL‐33); PTX3, pentraxin‐3; WFDC2, WAP 4‐disulfide core domain 2. The table presents the RV coefficients – calculated between the expression of all genes identified within each of the top‐ranked gene ontology pathways and the levels of plasma protein biomarkers – for those proteins that demonstrated at least a 0.2 or > relationship with one of the five top‐ranked gene ontology pathways (all p < 0.001). The number of genes in each gene ontology pathway is listed in brackets. FE is the fold enrichment over the gene ontology database.

Co‐expression analysis

To further delineate the inter‐relationships between genes within each of the five GO categories and their relationship with the protein biomarkers, and how these might relate to prognosis, we created Pearson correlation heatmaps for each category, for survivors and non‐survivors separately (Figure  and online supplementary Figure  ). To facilitate identification of altered co‐expression relationships between survivors and non‐survivors we also created differential co‐expression plots (online supplementary Figures S7 and ). The direction of gene co‐expression was largely conserved between survivors and non‐survivors, partly speaking to the robustness of the methodologies utilized. However, these analyses highlighted several individual examples of altered gene co‐expression, e.g. the iron transporter, transferrin receptor was upregulated in non‐survivors (Figure  ) and correlated to IL6ST (gp130) expression only in non‐survivors (Figure  and online supplementary Figure  ). Gene expression within the UPP pathway was subject to numerous shifts in the correlative structure (Figure  and online supplementary Figure ) with genes like HECTD3 (up‐regulated in non‐survivors) becoming more positively co‐expressed with other members of the pathway (except for SMURF2 and RNF216). Likewise, the correlations between the seven protein biomarkers (FGF23, sST2, ADM, HEPC, PTX3, WFDC2 and IL‐6) tended to be stronger in non‐survivors, except for IL‐6 (Figure  vs. online supplementary Figure  ). The correlation between the seven individual proteins and individual genes varied substantially (online supplementary Figures and ) and this provided insight into some likely connections between the transcriptome profile and survival (see Discussion).

Figure 3

Top gene ontology pathway gene and protein biomarkers inter‐relationships. Gene expression was correlated with the top protein biomarkers using Pearson correlation coefficients for the two largest significant pathways ‘adaptative immune response’ and ‘proteasome‐mediated ubiquitin‐dependent protein catabolic process’ – the other three top ranked gene ontology pathways are presented in online supplementary Figure  . The protein values are enclosed by a black oblong box. Data are plotted separately (grey boxes) for survivors (n = 626) and non‐survivors (n = 318). Correlation values are represented by colour and are plotted for significant correlations (Bonferroni corrected threshold, p < 1.5 × 10−4).

Using mortality‐related transcript patterns to identify new drug targets

A search, using the combined 60 most up‐regulated and 60 most down‐regulated genes between survivors and non‐survivors that mapped to the iLINCS database, against >2500 in vitro drug signatures, found 45 matching compounds (online supplementary Table  ). This included 29 that reversed the mortality associated gene expression pattern, and these drugs had 47 known protein targets (<1 µM potency). Combined pathway analysis of the drug targets and the top regulated mortality associated genes established that the 47 proteins mapped to biological processes identified by the CHF mortality related transcriptome. This included the ‘immune response‐activating cell receptor signalling’ and ‘lymphocyte activation’ pathways (Figure  ).

Figure 4

Pathway level overlaps between survival‐related genes and the protein targets of drugs that reverse the survival‐related gene expression signature in vitro. A network of significant pathways (p = 10−18–10−4) was derived using metascape.org. The input genes were the 120 chronic heart failure (CHF) survival‐associated genes used to identify drugs that regulate the CHF signature in vitro (https://clue.io/) and the 47 protein targets of the 29 drugs (identified using PubChem and the small‐molecule suite) that reverse the CHF signature in vitro. There are two identical plots. The large plot, on the left‐hand side, presents the significant pathways for this combined gene list. Edges represent connected gene ontology biological processes (>0.3), and nodes within each cluster are coloured/named by their most statistically enriched gene ontology term, scaled in size by the total number of terms represented. The smaller plot, on the right‐hand side, is same network structure but now colour coded by input list membership. This identifies if the drug targets appear within CHF survival‐associated pathways (pathways common to Figure  ) or whether the drug target falls within a pathway more indirectly connected with the patient transcriptomic signature. Each node is presented as a pie‐chart, with the ‘slices’ coloured and scaled to indicate which gene list the terms originate from, and what proportion the lists contribute to the ontology groupings.

Discussion

Whole blood transcriptomic profiling is a potentially powerful methodology to study genome‐wide biology directly in patients. Here, in the first prospective study of this type in heart failure, we provide novel insights into molecular processes associated with adverse CV prognosis. Specifically, we show that such profiles differ at baseline between those patients with CHF who died of a CV cause (over a median follow‐up of 21 months) compared to those that survived, with more than 1 in 20 of the detectable protein‐coding transcriptome showing a difference. Although the overall fold differences in gene expression between survivors and non‐survivors were small, these differences are likely to be driven by a subset of immune cells and therefore the fold change values would be ‘diluted’ by other cell types contributing to the global RNA profile. Thus, the magnitude of the measured differences should not be equated with their biological significance. Indeed, the differences could be mapped to specific biological pathways, several of which have been postulated to be disrupted in CHF. Importantly, we show that the differentially regulated transcriptome is strongly related to clinical and protein predictors of survival. While these findings are correlative, they highlight potential mechanisms and intermediaries through which the risks related to clinical factors are mediated. Finally, we provide initial evidence that the mortality‐associated RNA pathways in CHF can be reversed in vitro by drugs, some already approved for clinical use, identifying them – and their targets – as being potentially relevant to drug development for CHF.

Altered biological pathways and prognosis in heart failure

In cross‐sectional analysis, lymphopenia has been associated with more severe NYHA class and poorer outcomes in CHF. , Our study provides further detailed evidence that CHF patients with a poorer prognosis not only have fewer T‐ and B‐cell, but they also have an altered immune cell phenotype. Specifically, T‐cell receptor (TCR) signal‐transducing molecules such as CD3 epsilon and gamma, as well the co‐stimulatory receptor CD28, are down‐regulated in non‐survivors. Further, a number of genes encoding protein adaptors of TCR signals (such as Lck, Zap70, PKC‐θ, ITK and SKAP1), as well as genes that act to amplifying TCR signals such as CD6 were all down‐regulated (Figure  ). Many of these ‘adaptative immune response’ pathway members, and those involved in TCR signalling in particular, were negatively correlated (Figure  ) with both arginase gene transcripts (ARG1 and ARG2); which were 14%–16% more abundant in those that subsequently died (online supplementary Table ). Arginase degrades arginine and decreased arginine, preceded by a measurable increase in arginase activity, is associated with poorer 6‐month outcomes in patients with ST‐elevation myocardial infarction. Arginase gene expression, induced by lactate production consequent to hypoxia, can limit T‐cell proliferation. A negative association between ARG2 and the pro‐survival IL‐7 receptor (IL7R/CD127) was observed in survivors (Figure  ) potentially linking arginine metabolism and IL7R signalling role in T‐cell survival and memory responses. These observations provide support for altered metabolism contributing to the altered T‐cell molecular profile, which our transcriptome‐wide analysis indicates is a central feature of those patients that died.

Increased net catabolic activity has long been associated with poorer prognosis in heart failure and our transcriptomics approach also identified that up‐regulation of the UPP pathway occurs in circulating blood cells (Figure  ). Siah E3 ubiquitin protein ligase 2 (SIAH2) was one of the most up‐regulated genes in the non‐survivor group (28%, p = 9 × 10−7) and it mediates degradation of heme oxygenase‐1, a cardioprotective enzyme involved with immunomodulation. HECTD3 – which has a role in both cardiac hypertrophy and T helper‐17 development  – was expressed to a greater extent in those that died, showing a stronger (positive) co‐expression (online supplementary Figure ) with other members of the activated UPP pathway (Figure  ). In contrast SMURF2 and RNF216 (two down‐regulated genes in those that died) were negatively associated with HECTD3 (Figure  ). This is interesting as SMURF2 is a E3 ubiquitin ligase targeting degradation of pro‐inflammatory TGFB/SMAD signalling, while RNF216 contributes to ubiquitin‐mediated protein degradation of TLR4. Presumably targeted acceleration of protein degradation, through UPP, is required to facilitate proinflammatory processes and these more specific events may represent more tractable therapeutic targets.

It is plausible that one driver of the catabolic profile in those that died was disrupted iron homeostasis. As discussed below, altered iron regulatory protein, HEPC (which regulates ferroportin capacity ) and induction of erythrocyte development revealed a number of key observations. For example, there was a strong positive association with the up‐regulated transferrin receptor (which imports iron into cells) and up‐regulated APFL (aprataxin PNK‐like factor) – a ubiquitous DNA damage response enzyme potentially responding to inappropriate iron accumulation. These observations indicate that molecular processes associated with cachexia in CHF, previously identified in cardiac and peripheral tissues, are also evident in the whole blood transcriptome where a pervasive role for disrupted iron homeostasis appears across several of the altered pathways.

Do the mortality‐associated protein biomarkers drive the transcriptional responses to chronic heart failure?

As the seven protein biomarkers associated with the transcriptomic changes co‐varied with each other (online supplementary Figure ), it is more challenging to identify which proteins are directly linked to the transcriptome differences, particularly as most have plausible mechanistic roles in CHF. For example, PTX3 can skew T‐cell differentiation towards the IL‐17 producing phenotype while ADM has a positive impact on the endothelial integrity and anti‐inflammatory properties. FGF23, a protein secreted that regulates calcium and phosphate homeostasis, is produced by osteocytes and cardiac cells, and promotes hypertrophy and fibrosis in a calcineurin‐dependent manner. WFDC2 (also known as HE4) suppresses osteopontin and promotes cell survival in an osteopontin‐ and interferon‐γ‐dependent manner. WFDC2 is also linked to increased ARG1 activity and immune suppression.

Exploration of single transcript and protein biomarker relationships does, however, provide several clues as to the order of events, particularly for the altered T‐cell biology and iron homeostasis (Figure  and online supplementary Figure  ). For example, HEPC links to both altered UPP activity and altered T‐cell biology. HEPC demonstrated a strong association with the UPP and ‘erythrocyte development’ pathways (Figure  and online supplementary Figure  ). Notably, the association between HEPC protein and expression of genes from the UPP pathway was more striking in the non‐survivors (Figure  B, e.g. GCLC, PSME4, SIAH2, TBL1XR1, YOD1 and ZER1). In contrast, mucolipin (MCOLN1) was negatively associated with HEPC and one of the few up‐regulated genes in the ‘adaptive immune response’ pathway. MCOLN1 has an important role in calcium‐mediated cardiac remodelling and is also implicated in iron homeostasis. Together, these examples, and those listed above, imply that alterations in iron status may be a key integrator for the immune‐related changes in transcriptome and catabolic pathway activation. sST2 was also up‐regulated in those that died, consistent with its relationship with outcome following myocardial infarction, where it subdues IL‐33 signalling. sST2 levels were positively correlated with the up‐regulated erythrocyte development pathway and of note, in the non‐survivor group, the increased ‘erythropoiesis development’ signature was accompanied by reduced haemoglobin (p < 1 × 10−6) and iron (Table  ). IL‐33 can be cardioprotective and is released upon tissue damage, stimulating T‐helper type 1 (Th1), Th2 and pro‐inflammatory Th17 cells, such that loss of IL‐33 through increased sST2 may help explain the observed immune‐related transcriptome profile. Overall, our analysis indicates that each of these seven survival‐related proteins could drive the transcriptome response, coupling immunological and iron homeostasis together with clinical risk factors.

Using mortality‐related transcript patterns to identify new drug treatments

Not all circulating protein biomarkers or regulated genes make relevant therapeutic targets. For instance, despite being modulated in CHF, targeting tumour necrosis factor‐α or IL‐6 pathways has not led to improved outcomes in heart failure, or been found to reduce major CV events. Notably, the present analysis did not rank either of these cytokines as being a dominant feature of the mortality‐related pathways. Instead, our analysis indicates that antagonizing FGF23 signalling and/or down‐regulating sST2 expression (or its binding to IL‐33) represent more enticing targets. Beyond these individual proteins, use of transcriptome signatures represents a promising alternative strategy to discover therapeutics and pathways that can modulate diseases. , , In fact, many of the protein targets of the 29 drugs that reversed the mortality‐related transcriptome in vitro, are not regulated in patients, but rather represent members of the same pathways regulated in those that died (Figure  ). Notably, several positively acting drugs are known to target heart failure‐related molecular processes (online supplementary Table ), and several have in vivo support for their potential utility in CHF, including a monoamine oxidase inhibitor which can be cardioprotective, forskolin which elevates cAMP and cardiac function, and the experimental epigenetic drug, chaetocin (a histone lysine methyltransferase inhibitor), which appears to be cardioprotective. Our findings encourage a further appraisal of the potential benefit in CHF of some of these drugs and pathways.

Limitations

Our study has several limitations. Despite being the largest whole blood genome‐wide transcriptomics study in CHF conducted in a very well‐characterized cohort, further large scale independent transcriptomic studies are required to validate our observations. An important component of our strategy was to adjust for differences in cell counts between survivors and non‐survivors. We modelled changes in individual white cell subtypes using the transcriptomics data and studies using fresh blood samples from patients with CHF are needed to fully explore the changes we predicted. Most patients in BIOSTAT‐CHF had heart failure with reduced ejection fraction and whether our findings extend to those with preserved ejection fraction remains to be determined. Finally, our results are correlative and experimental studies are required to confirm the direct involvement of the various pathways we have identified in adverse prognosis in CHF.

Conclusion

In conclusion, our analysis provides evidence that whole blood transcriptomics can identify molecular pathways that associate with CV mortality in CHF. These reflect clinical factors that are associated with poor prognosis and may be driven by specific intermediary protein factors. We demonstrate the potential utility of such transcriptomic profiling for identifying novel therapeutic targets for CHF.

Appendix S1. Supporting Information

Click here for additional data file.

Supplementary Tables

Click here for additional data file.