A 6-gene panel as a signature to predict recovery from advanced heart failure using transcriptomic analysis.

To the Editor:

Heart failure (HF) is a global health problem with a high mortality rate. The various stimuli associated with HF can cause maladaptive remodeling, gradually weakening cardiac functions. Left ventricular assist device (LVAD) can improve advanced HF patients' cardiac functions, termed “reverse remodeling”. However, only a small percentage of patients who have received LVAD have experienced reverse remodeling. Moreover, heart failure relapses in many patients after LVAD explantation, indicating that current biomarkers may not be sufficient to provide a complete view of pathological status. In this study, we identified a 6-gene panel (6GP; positive correlation: MYC, FOXO1, and ZFP36; negative correlation: LONRF1, FRZB, and NPPA), which demonstrated outstanding potential for predicting cardiac recovery, by combining weighted gene co-expression analysis (WGCNA) and differential gene expression analysis. 6GP was highly correlated with key molecular signatures related to cardiac recovery, including myogenesis, axonogenesis, and epidermis development.

The RNA-seq dataset GSE46224 contained quantitative expression data of 19,481 genes and was used as the discovery set. It was obtained from the ventricular tissue of eight patients with ischemic cardiomyopathy (ICM) and eight patients with non-ischemic cardiomyopathy (NICM) treated with a left ventricular assist device (LVAD). Samples were collected at two-time points: before LVAD implantation (pre-LVAD) and during heart transplantation (post-LVAD). Patients receiving LVAD treatment were divided into responders and non-responders. Criteria for dividing responders and non-responders and the detailed distribution of responders and non-responders in ICM and NICM patients were described in Supplementary Data section 1.2. The validation set GSE21610 (9 ICM patients, 21 NICM patients, and 8 healthy donors) had a probe set of 20,161 genes.

GSE46224 dataset was used for WGCNA. After construction of the co-expression network, 14 co-expression modules were obtained, among which the “brown” module had the highest correlation with the target clinical trait: post-LVAD (R2 = 0.76 and P-value = 3 × 10−8) (Fig. 1A). Therefore, it was identified as the LVAD-induced co-expression module.

Figure 1

Establishment of the 6-gene panel (6GP). (A) The association between traits and modules. Colors in the heatmap denote the strength and direction of correlations between co-expression modules and corresponding clinical traits (red and blue denote positive and negative correlations, respectively). (B) The heat map of differentially expressed genes (DEGs) between responders and non-responders in GSE46224. The red, white, and black blocks indicate high, medium, and low expression levels. The vertical axis represents the significant DEGs, while the horizontal axis represents the post left ventricular assist device treatment (post-LVAD) samples. Color legend on the top of the heatmap denotes the two comparing groups (light red and cyan denote responders and non-responders, respectively). (C) Receiver operating characteristic (ROC) curve of the signature score of 6GP in classifying responders and non-responders in GSE46224. (D) The signature score of 6GP in patients with ischemic cardiomyopathy (ICM), non-ischemic cardiomyopathy (NICM), and the whole GSE46224 dataset including both ICM and NICM (Whole). Statistical analysis was based on Student's t-test. (E) Gene set enrichment analysis plot of the post-LVAD group in GSE21610 demonstrating the physiological effects of 6GP. Red curves denote the running sum of the enrichment score. The red-to-blue panel in each subgraph represents a gene list sorted by correlations with the signature score of 6GP (from highly positive to highly negative).

To select pivotal modulators that determine patients' recovery status, we selected 20 hub genes from the LVAD-induced co-expression module based on their connectivity in the co-expression and protein–protein interaction network (Supplementary Data Section 2.3). The hub genes were further screened by the least absolute shrinkage and selection operator (LASSO) with a well-established leave-one-out cross-validation method. Three of the hub genes: MYC, FOXO1, and LONRF1, were identified as key modulators that determine the outcomes of LVAD treatment by LASSO. LONRF1 was the most influential variable and negatively correlated with the response variable. In contrast, FOXO1 and MYC were highly positive variables that prompt patients to become responders (Supplementary Data Section 2.4).

To find out more genes related to cardiac recovery and improve the effectiveness of the gene panel, we further conducted a differential gene expression analysis between responders and non-responders in GSE46224. Compared with the non-responders, 181 up-regulated and 86 down-regulated differentially expressed genes (DEGs) were identified in the responders (Fig. 1B).

Among the 181 up-regulated DEGs in responders of GSE46224, 13 were also identified in the LVAD-induced co-expression module by WGCNA. To select genes that can be added to the gene panel from the 13 DEGs, another dataset from the mice model: GSE107568 was applied. Among the 13 up-regulated DEGs, ZFP36 had the highest area under the receiver operating characteristic curve (AUC) in classifying mechanically unloaded samples in GSE107568 (AUC = 0.889) and was added into the predictive gene panel. CXCL2 and GPR183 also passed the screening by GSE107568 but were excluded due to their pro-inflammatory effects (Supplementary Data Section 2.5).

From the 86 down-regulated DEGs in responders, we further identified two important genes: NPPA and FRZB. They were up-regulated in patients with heart failure compared with healthy donors and were persistently elevated in non-responders to LVAD compared with responders. In another clinical dataset, GSE21610, their AUC was also higher than 0.8 when classifying heart failure patients and healthy donors (Supplementary Data Section 2.6). Therefore, three additional genes - ZFP36, NPPA, and FRZB – were added to the predictive gene panel.

Collectively, combining WGCNA and differential gene expression analysis, we obtained the predictive 6GP. To further evaluate the correlation between 6GP and the rehabilitation of patients with advanced HF, we performed scoring of 6GP signature based on the “singscore” method. According to the previous analysis, patients with higher levels of MYC, FOXO1, ZFP36, and lower LONRF1, NPPA, FRZB levels would recover better. Transcriptomes with higher concordance to this pattern would be given higher scores.

The AUC of scoring of 6GP was 0.95 when distinguishing responders and non-responders in GSE46224, indicating that the 6GP has good potential in the classification of LVAD responders and non-responders (Fig. 1C). It is worth noting that although there were differences in gene expression between ICM and NICM, scores of responders were higher than in non-responders across all patient subtypes, indicating that the classifying ability of this signature was robust across both types of cardiomyopathy (ICM and NICM) (Fig. 1D).

To further verify that 6GP reflects the presence of cardiac recovery, we performed the same 6GP signature scoring on the post-LVAD samples in an independent dataset, GSE21610. Gene set enrichment analysis (GSEA) revealed that the 6GP signature score in the GSE21610 dataset was positively correlated with molecular signatures related to myogenesis, axonogenesis, epidermis development, cell–cell adhesion, and adrenergic signaling. At the same time, unfolded protein response was suppressed (Fig. 1E). Although GSE21610 did not provide post-LVAD echocardiography data, we could infer that the high score of 6GP was associated with myocardial regeneration, better contractility, active development of blood vessels, and reduced oxidative stress.

Overall, this study derived a 6GP to optimize therapies and recovery assessing strategies from co-expression analysis and differential gene expression analysis. Among the 6GP, we revealed the positive correlations between MYC, ZFP36, FOXO1, and patients' responses. Meanwhile, we found that LONRF1, an LVAD-induced co-expression module hub gene, was negatively correlated with patients' responses. Therefore, further studies on LONRF1 may shed light on improving the recovery rate of LVAD treatment. The other two pivotal genes: NPPA and FRZB, were validated biomarkers for HF. The signature score of 6GP achieved high performance in classifying responders and non-responders to LVAD. GSEA based on another dataset: GSE21610 revealed positive correlations between the signature score of 6GP and the regeneration of myocardium, blood vessel development, resistance to oxidative stress, and contractility of the heart.

This correlation pattern provides a general idea of the mechanisms of cardiac recovery from advanced heart failure and a candidate list of drug targets and biomarkers for optimizing disease management. Further studies with larger cohorts and comprehensive clinical records would better confirm our findings.

Author contributions

Conceived and designed the study: SJZ, BLK. Performed the study: SJZ. Analyzed the data: SJZ. Wrote the paper: SJZ, BLK. All authors have read the manuscript.

Conflict of interests

Authors declare no conflict of interests.

Funding

City University of Hong Kong (No. 9610430), which is funded by the Research Grants Council (RGC); - Research Talent Hub (RTH) 1–5.