Transcriptomic Analysis Identifies the Effect of Beta-Blocking Agents on a Molecular Pathway of Contraction in the Heart and Predicts Response to Therapy.

Over the last decades, beta-blockers have been a key component of heart failure therapy. However, currently there is no method to identify patients who will benefit from beta-blocking therapy versus those who will be unresponsive or worsen. Furthermore, there is an unmet need to better understand molecular mechanisms through which heart failure therapies, such as beta-blockers, improve cardiac function, in order to design novel targeted therapies. Solving these issues is an important step towards personalized medicine. Here, we present a comprehensive transcriptomic analysis of molecular pathways that are affected by beta-blocking agents and a transcriptomic biomarker to predict therapy response.

G-protein–coupled receptors are the most commonly targeted proteins of recently designed drugs in the cardiovascular field, making them a key component of pharmacogenomic investigations and genetic variability studies (1). One of the best studied G-protein–coupled receptors is the beta-adrenergic receptor (AR) (1), with β1- and β2-AR both being expressed on human cardiomyocytes. Stimulation of ARs, in particular the β1-AR, induces increased cardiac inotropy and chronotropy (1). The β3-ARs have been shown to have negative regulatory functions on inotropy and cardiac reserve 2, 3 through Gi coupling to cyclic guanine monophosphate-nitric oxide (4). Although beta-adrenergic stimulation is a major compensatory mechanism in the acute setting such as traumatic hypovolemia, it appears to worsen ventricular function and outcome in conditions with limited metabolic and physiological reserves, such as heart failure (HF) (1). Accordingly, beta-blocking agents were developed to partially antagonize beta-adrenergic “excess” of norepinephrine at the cardiomyocyte level (1).

Since the first discovery of beneficial effects of beta-blocker therapy in a small case series of 7 patients with HF by Waagstein et al. (5) in 1975, a sequence of large clinical trials 5, 6, 7, 8, 9 has confirmed clinical improvement with beta-blockade in HF and suggested various hypotheses on how beta-blockade in the failing heart may improve outcomes. These hypotheses included effects of beta-blockade through activation of myocardial contractile proteins and sarcoplasmic reticulum calcium-dependent ATPase (SERCA) activity, as well as alteration of gene expression 10, 11, 12.

In the field of gene expression analysis, Lowes et al. (10) made the important observation that functional improvement during beta-blocker therapy measured by improved ejection fraction (EF) (increase by 18.8 ± 1.8%) was associated with overexpression of SERCA and α-myosin heavy chain (MYH), whereas there was a decrease in expression levels of β-MYH. Yasumura et al. (11) made a similar observation in a clinical study, in which they treated patients with dilated cardiomyopathy for 4 months with beta-blockers. Improvement of EF during therapy with beta-blocking agents was associated with overexpression of SERCA and phospholamban.

While the previously mentioned investigational approaches used polymerase chain reaction (PCR) to gain valuable information about the molecular effects of beta-blocking agents on specific candidate genes, our group sought to expand the current knowledge by applying microarrays, a technology that evaluates expression levels of all genes in a given individual. Therefore, it allows the discovery of new genes and pathways in a more comprehensive approach 13, 14, 15, 16, 17, 18 that has proven useful to delineate diagnosis and prognosis in HF populations 15, 16, 18. While typically only genes of interest or candidate genes are investigated with polymerase chain reaction, microarray technology analyzes the entire transcriptome of about 30,000 genes in 1 experiment (19) and, therefore, provides a less biased approach of gene discovery.

Methods

Patient population

To identify genes that undergo expression changes during treatment with beta-blocking agents, we first analyzed endomyocardial biopsies (EMBs) obtained from patients with new-onset HF (n = 43), who were treated with beta-blockers (n = 30) versus alternative standard therapy (n = 13). This cohort has been previously described and analyzed for prognostic information (15). EMBs were obtained from a biorepository containing samples from patients with new-onset HF 15, 16, 18. In brief, transvenous EMBs were obtained from the right interventricular septum and immediately flash frozen in liquid nitrogen for storage and subsequent microarray analysis. All patients gave their written informed consent to participate in this study.

New-onset HF was defined as onset of clinical symptoms of HF within the past 6 months at the time of diagnosis. HF was diagnosed on the basis of clinical signs and symptoms as recommended by the American Heart Association and the American College of Cardiology (20). All patients had been treated with beta-blockers or alternative therapy for a maximum of 6 months at the time when biopsy samples were obtained. Metoprolol (tartrate and succinate) and carvedilol were used as beta-blocking agents at the maximum tolerated dose, whereas alternative standard therapy included angiotensin-converting enzyme inhibitors, aldosterone antagonists, and diuretic agents. Patients who received alternative therapy had contraindications to beta-blocking agents, such as severe chronic obstructive pulmonary disease, sinus bradycardia, or atrioventricular node conduction disorders.

The cohort was matched on the basis of age, sex, hemodynamic parameters, and medical therapy. To avoid different types of cardiomyopathy from being a possible confounding factor for gene expression analysis, only samples from patients with idiopathic dilated cardiomyopathy were investigated. Idiopathic dilated cardiomyopathy was a diagnosis of exclusion, after extensive workup including cardiac catheterization and special immunohistochemical stainings were performed (16).

After identifying genes that were affected by carvedilol or metoprolol, we sought to identify transcriptomic changes specifically related to improved outcomes. Therefore, in a case-control fashion on the basis of the same baseline parameters as for our first analysis, we investigated gene expression changes that were unique to patients who were responsive to beta-blocking agents. For this second analysis, within the previously described participants who received beta-blocking agents, we identified a group with good prognosis (n = 17) and a group with poor prognosis (n = 13) (Figure 1). Good prognosis was defined as event-free survival for at least 5 years after diagnosis of HF. Poor prognosis was defined as having an event within 2 years after symptom onset. Events included death, requirement for left ventricular assist device placement, or cardiac transplant (15).

Figure 1

Study Design

We analyzed a cohort of 43 patients with idiopathic dilated cardiomyopathy for molecular changes that are induced by beta-blocking agents. Among the patients on beta-blocking agents, we selected patients with poor (n = 13) versus good prognosis (n = 17) to identify gene expression changes in patients who improved during therapy with beta-blockers. Two-thirds of the data were used as a train set to develop the biomarker (poor prognosis n = 8; good prognosis: n = 11) and one-third of the data was used as an independent test set for validation of the molecular signature (poor prognosis: n = 5; good prognosis: n = 6).

Ribonucleic acid isolation from EMBs

We extracted and hybridized total RNA from EMBs as previously described 15, 16, 17. All EMBs were obtained at the time of diagnosis of new-onset HF and were immediately flash-frozen in liquid nitrogen. The average size of tissue samples was 2 mm. Samples were homogenized with the MM 301 Mixer Mill (Retsch, Inc., Newtown, Pennsylvania) (catalog no. 85120). Ribonucleic acid (RNA) extraction of total RNA was done with Trizol reagent in combination with the Micro-to-Midi Total RNA Purification System (Invitrogen, Carlsbad, California) (catalog no. 12183-018). Concentration and integrity of total RNA were measured with the Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, California). All RNA samples exhibited intact 28S and 18S ribosomal RNA on denaturing agarose gel electrophoresis, and the 260-/280-nm absorbance readings fell within the acceptable range of 1.8 to 2.1. Subsequently, RNA was pre-processed with the Ovation Biotin RNA Amplification and Labeling System (NuGEN Technologies, Inc., San Carlos, California) (catalog no. 2300-12).

Microarray analysis

Without prior amplification, samples were hybridized to the Human Genome U133 Plus 2.0 Array from Affymetrix (Santa Clara, California). Average background and noise of all chips registered within acceptable ranges, and hybridization efficiencies were similar for all samples.

Statistical analysis

Robust multiarray average was used for normalization of microarray data and significance analysis of microarrays (SAM) for identification of differentially expressed genes in patients taking beta-blocking agents (n = 30) versus alternative standard therapy (n = 13). SAM identifies statistically significant genes by assimilating a set of gene-specific t tests. This highly robust statistical algorithm provides a q value of statistical significance for each individual gene within a phenotype. The resulting set of genes was further processed with Meta Core pathway analysis from GeneGo Inc. (bioinformatics software, St. Joseph, Michigan). A major strength of GeneGo Metacore is to evaluate very comprehensively if genes are connected to a certain pathway through ligands or receptors within a signaling cascade, whereas gene ontology (GO) biological process terms will usually only list a direct involvement of a gene in a molecular process.

Organ- and species-specific pre-filtering was performed before network analysis to identify pathways that are truly inter-related within the human heart. To obtain a false discovery rate of a pathway, the p value was adjusted for multiple comparisons. For that purpose, the software applies a false discovery rate adjustment. This adjustment of the p value of the enrichment of a pathway takes into account the number of tests that are being performed during the enrichment analysis. In addition, a z score was calculated for each network, which reflects the saturation with genes from the experiment. A high z score indicates a network that contains a large amount of genes from the experiment, therefore suggesting that this particular network was significantly affected. The g score takes into account how many canonical pathways were involved to create the network by modifying the z score on the basis of the number of linear canonical pathway fragments contained within the network. A network with a high g score is saturated with objects from the investigated gene list and contains a large number of canonical pathway fragments.

To develop a marker that predicts responsiveness to beta-blockers, we applied Misclassification Penalized Posteriors (MiPP) to compare patients who have been on beta-blocking agents and had good prognosis (event-free survival for at least 5 years) versus patients who had poor prognosis.

Poor prognosis was defined as having an event, such as death, cardiac transplant, or left ventricular assist device placement within 2 years after diagnosis with HF. The MiPP package is an application in the R environment, which employs the libraries MASS for algorithms such as linear discriminant analysis and e1071 for support vector machine 21, 22. This software sequentially adds genes to a classification model on the basis of the Misclassification-Penalized Posteriors principle, which takes into account the likelihood that a sample belongs to a given class by using posterior probability of correct classification. Linear discriminant analysis uses a linear combination of features, which best separates 2 or more classes 21, 22.

In support vector machine algorithms for classification, the input data is plotted as 2 vectors in an n-dimensional space, and a virtual hyperplane is created that best separates the 2 phenotypes. This hyperplane is then used to classify samples with unknown phenotypes.

We developed a molecular signature in two-thirds of the data using 5-fold cross validation and performed additional validation in an independent sample that contained the other one-third of the data.

To evaluate if distinct models are generated from additional random splits, we performed 50 random divisions to develop individual classification models, which were then validated in 10 independent splits. Furthermore, we calculated mean sMiPP for every given gene model, an additional parameter for performance, which approximates 1 with increasing accuracy.

Results

Table 1 illustrates baseline variables of the initial set of samples in our first analysis including patients with new-onset HF treated with beta-blocking agents versus alternative standard therapy. There were no significant differences between the 2 groups (evaluated by Student t test and Fisher exact test). Importantly, the cohort was balanced between male and female subjects (70% vs. 62%). Hemodynamic parameters, including EF, left ventricular internal dimension-diastole, pulmonary artery pressure, and pulmonary capillary wedge pressure, were similar between groups. In both patient populations, a predominance of angiotensin-converting enzyme inhibitor and diuretic agent use was observed. A predictive molecular signature of therapy responsiveness was then developed in a train set, which contained two-thirds of the cohort treated with beta-blocking agents. One third of data was used as independent test set, in which accuracy of the biomarker was evaluated. Table 2 depicts baseline conditions of the train set (two-thirds of the cohort taking beta-blocking agents). The prognostic outcome in male versus female patients was comparable. Hemodynamic measurements of EF, left ventricular internal dimension-diastole, pulmonary artery pressure, and pulmonary capillary wedge pressure were again similar between both groups, as was medication use.

Table 1

Baseline Conditions in Patients With New-Onset Heart Failure on Beta-Blocking Agents Versus Alternative Therapy

	Beta-Blocking Agent(n = 30)	No Beta-Blocking Agent(n = 13)	p Value
Age, yrs	50 ± 5	45 ± 3	0.41
Male	21 (70)	8 (62)	0.73
LVEF, %	21 ± 2	29 ± 5	0.06
LVIDD, cm	6.3 ± 0.3	5.7 ± 0.6	0.31
PAP, mm Hg
Systolic	38 ± 3	37 ± 4	0.74
Diastolic	18 ± 2	17 ± 2	0.62
PCWP	16 ± 2	12 ± 2	0.12
Medications
ACE inhibitor	22 (73)	7 (54)	0.29
Aldosterone antagonist	8 (27)	0	0.08
Diuretic agent	21 (71)	11 (85)	0.45
IV inotropic therapy	0	0	1.00

Values are mean ± SD or n (%).

ACE = angiotensin-converting enzyme; IV = intravenous; LVEF = left ventricular ejection fraction; LVIDD = left ventricular internal dimension-diastole; PAP = pulmonary artery pressure; PCWP = pulmonary capillary wedge pressure.

Table 2

Baseline Conditions of the Train Set Containing Two-Thirds of the Population With New-Onset Heart Failure on Beta-Blocking Agents Who Had Good Versus Poor Prognosis

	Good Prognosis(n = 11)	Poor Prognosis(n = 8)
Age, yrs	42 ± 2	42 ± 4
Male	7 (64)	5 (63)
LVEF, %	20 ± 2	23 ± 3
LVIDD, cm	6.2 ± 0.3	5.8 ± 0.6
PAP, mm Hg
Systolic	35 ± 4	40 ± 4
Diastolic	15 ± 2	19 ± 3
PCWP	14 ± 2	16 ± 3
Medications
ACE inhibitor	9 (82)	5 (63)
Aldosterone antagonist	3 (27)	3 (38)
Diuretic agent	0	1 (13)
IV inotropic therapy	0	0

Values are mean ± SD or n (%).

Abbreviations as in Table 1.

Overexpressed genes and pathways in patients receiving beta-blocker therapy with either metoprolol or carvedilol

First, we analyzed EMBs from patients treated with beta-blockers (n = 30) versus alternative standard therapy (n = 13) to identify gene patterns that were overexpressed in patients who received beta-blocker therapy with metoprolol or carvedilol. Using SAM, we identified a total of 94 transcripts that were overexpressed (q value <5%, fold change [FC] >1.2) (Figure 2, Table 3) in patients treated with beta-blocking agents versus alternative therapy. The positive FC in patients on beta-blocking agents ranged from 1.14 to 2.92 and averaged 1.47. The q value for this analysis ranged from 0% to 4.1% and averaged 2.9%. Interestingly, no downregulated genes were detected. A heat map created by an unsupervised clustering approach illustrates a distinct set of genes between the 2 phenotypes (Figure 3). Overexpressed genes were further evaluated with GeneGo Metacore pathway analysis. We discovered 15 developmental and metabolic pathways (Table 4) that were activated in patients who had been treated with beta-blocking agents. Those pathways included metabolic processes involving phosphate, developmental processes, and regulation of apoptosis. The percentage values listed for each GO process in Table 4 represent the percentage of objects on that specific particular subnetwork that have been linked to that GO process. Similarly, the p value listed next to the percentage is calculated on the basis of the total number of objects on that particular subnetwork, the number of those objects on that subnetwork associated with the GO process, and of the objects associated with that GO process in the database. It uses the same calculation method as the other p values throughout MetaCore (hypergeometric distribution).

Figure 2

Significance Analysis of Microarrays Plot of Differentially Expressed Genes in Patients on Beta-Blocking Agents Versus Alternative Therapy in IDCM

There were 94 genes differentially expressed in patients treated with beta-blocking agents versus alternative therapy (q value <5%, fold change >1.2). No downregulated genes were detected in patients treated with beta-blocking agents. The 94 overexpressed genes are depicted in red. IDCM = idiopathic dilated cardiomyopathy.

Table 3

Differentially Expressed Genes in Patients Treated With Beta-Blockers Versus Alternative Standard Therapy for Heart Failure

Probe Set ID	Gene Symbol	Gene Title	Go Biological Process Term
204737_s_at	MYH7	MYH7, myosin, heavy chain 7	Cardiac muscle, beta
216265_x_at	MYH7	MYH7, myosin, heavy chain 7	Cardiac muscle, beta
214468_at	MYH6	MYH6, myosin, heavy chain 6	Cardiac muscle, alpha
209186_at	SERCA/ATP2A2ase	ATP2A2, ATPase	Ca++ transporting, cardiac muscle, slow twitch 2
212361_s_at	SERCA/ATP2A2ase	ATP2A2, ATPase	Ca++ transporting, cardiac muscle, slow twitch 2
212362_at	SERCA/ATP2A2ase	ATP2A2, ATPase	Ca++ transporting, cardiac muscle, slow twitch 2
239996_x_at	SERCA/ATP2A2ase	ATP2A2, ATPase	Ca++ transporting, cardiac muscle, slow twitch 2
1553992_s_at	NBR2	neighbor of BRCA1 gene 2 (nonprotein coding)	GTP binding
1554868_s_at	PCNP	PEST proteolytic signal containing nuclear protein	Cell cycle, protein ubiquitination
1556414_at	C21orf71	Chromosome 21 open reading frame 71	NA
1557383_a_at	NA	NA	NA
1560109_s_at	NA	NA	NA
1563498_s_at	SLC25A45	Solute carrier family 25, member 45	Transport, transmembrane transport
200821_at	LAMP2	Lysosomal-associated membrane protein 2	NA
201458_s_at	BUB3	Budding uninhibited by benzimidazoles 3 homolog (yeast)	Mitotic sister chromatid segregation, cell cycle, chromosome segregation
201534_s_at	UBL3	Ubiquitin-like 3	NA
201979_s_at	PPP5C	Protein phosphatase 5, catalytic subunit	Protein amino acid dephosphorylation, mitosis, response to morphine
202125_s_at	TRAK2	Trafficking protein, kinesin binding 2	Regulation of transcription from RNA polymerase II promoter
202932_at	YES1	v-yes-1 Yamaguchi sarcoma viral oncogene homolog 1	Protein amino acid phosphorylation, glucose transport, regulation of vascular permeability
202933_s_at	YES1	v-yes-1 Yamaguchi sarcoma viral oncogene homolog 1	Protein amino acid phosphorylation, glucose transport, regulation of vascular permeability
203212_s_at	MTMR2	Myotubularin-related protein 2	Protein amino acid dephosphorylation, protein tetramerization
203387_s_at	TBC1D4	TBC1 domain family, member 4	Regulation of Rab GTPase activity
203689_s_at	FMR1	Fragile X mental retardation 1	Negative regulation of translational initiation
204334_at	KLF7	Kruppel-like factor 7 (ubiquitous)	Regulation of transcription from RNA polymerase II promoter, axon guidance
205857_at	SLC18A2	Solute carrier family 18 (vesicular monoamine), member 2	Response to amphetamine, response to toxin, post-embryonic development
209337_at	PSIP1	PC4 and SFRS1 interacting protein 1	Interspecies interaction between organisms, regulation of transcription
210257_x_at	CUL4B	Cullin 4B	DNA repair, cell cycle
211552_s_at	ALDH4A1	Aldehyde dehydrogenase 4 family, member A1	Proline metabolic process, proline biosynthetic process, oxidation reduction
212008_at	UBXN4	UBX domain protein 4	Response to unfolded protein
212071_s_at	SPTBN1	Spectrin, beta, nonerythrocytic 1	Common-partner SMAD protein phosphorylation, actin filament capping
212214_at	OPA1	Optic atrophy 1 (autosomal dominant)	Inner mitochondrial membrane organization transport of mitochondrion, positive regulation of antiapoptosis
212598_at	WDFY3	WD repeat and FYVE domain containing 3	NA
212688_at	PIK3CB	Phosphoinositide-3-kinase, catalytic, beta polypeptide	Activation of MAPK activity, cell-matrix adhesion, calcium ion homeostasis
212764_at	ZEB1	Zinc finger E-box binding homeobox 1	Immune response, cell proliferation, regulation of mesenchymal cell proliferation, regulation of transforming growth factor beta receptor signaling pathway, negative regulation of epithelial cell differentiation
213169_at	SEMA5A	Sema domain, seven thrombospondin repeats (type 1 and type 1-like), transmembrane domain and short cytoplasmic domain, (semaphorin) 5A	Patterning of blood vessels, cell adhesion, cell-cell signaling, development
213251_at	SMARCA5	SWI/SNF-related, matrix-associated, actin-dependent regulator of chromatin, subfamily a, member 5	Chromatin silencing at rDNA, regulation of transcription from RNA polymerase II promoter, embryonic development
214666_x_at	IREB2	Iron-responsive element binding protein 2	Cellular iron ion homeostasis, post-embryonic development
214757_at	PMS2L2	Post-meiotic segregation increased 2-like 2 pseudogene	Mismatch repair
218082_s_at	UBP1	Upstream binding protein 1 (LBP-1a)	Angiogenesis, negative regulation of transcription
218197_s_at	OXR1	Oxidation resistance 1	Response to oxidative stress
218252_at	CKAP2	Cytoskeleton-associated protein 2	Apoptosis, cell cycle
218528_s_at	RNF38	Ring finger protein 38	NA
218658_s_at	ACTR8	ARP8 actin-related protein 8 homolog (yeast)	NA
219218_at	BAHCC1	BAH domain and coiled-coil containing 1	NA
220494_s_at	NA	NA	NA
220777_at	KIF13A	Kinesin family member 13A	Transport, microtubule-based movement
221273_s_at	RNF208	Ring finger protein 208	NA
221428_s_at	TBL1XR1	Transducin (beta)-like 1 X-linked receptor 1	Regulation of transcription, DNA-dependent
221472_at	SERINC3	Serine incorporator 3	Induction of apoptosis
221768_at	SFPQ	Splicing factor proline/glutamine-rich (polypyrimidine tract binding protein associated)	Alternative nuclear mRNA splicing, via spliceosome, DNA repair, response to DNA damage stimulus
222209_s_at	TMEM135	transmembrane protein 135	Protein folding, response to unfolded protein
222572_at	PDP1	Pyruvate dehydrogenase phosphatase catalytic subunit 1	Protein amino acid dephosphorylation
223282_at	TSHZ1	Teashirt zinc finger homeobox 1	Transcription, regulation of transcription, DNA-dependent, multicellular organismal development
224047_at	NA	NA	NA
224471_s_at	BTRC	Beta-transducin repeat containing	Ubiquitin-dependent protein catabolic process, signal transduction, regulation of I-kappaB kinase/NF-kappaB cascade, interspecies interaction between organisms, positive regulation of ubiquitin-protein ligase activity involved in mitotic cell cycle
224639_at	SPPL3	Signal peptide peptidase 3	Regulation of neuron apoptosis
224897_at	WDR26	WD repeat domain 26	NA
224945_at	BTBD7	BTB (POZ) domain containing 7	NA
224994_at	CAMK2D	Calcium/calmodulin-dependent protein kinase II delta	G1/S transition of mitotic cell cycle, regulation of cell growth, response to hypoxia, cardiac muscle contraction
225026_at	CHD6	Chromodomain helicase DNA binding protein 6	Chromatin assembly or disassembly, regulation of transcription, DNA-dependent, nervous system development
225328_at	NA	NA	NA
225538_at	ZCCHC9	Zinc finger, CCHC domain containing 9	NA
225603_s_at	C8orf83	Chromosome 8 open reading frame 83	NA
225785_at	REEP3	Receptor accessory protein 3	NA
225912_at	TP53INP1	Tumor protein p53 inducible nuclear protein 1	Induction of apoptosis, response to stress, cell cycle arrest
226004_at	CABLES2	Cdk5 and Abl enzyme substrate 2	Cell division, regulation of cell cycle
226134_s_at	NA	NA	NA
226558_at	LOC389834	Ankyrin repeat domain 57 pseudogene	NA
226771_at	ATP8B2	ATPase, class I, type 8B, member 2	ATP biosynthetic process
226797_at	MBTD1	Mbt domain containing 1	Transcription, chromatin modification
226799_at	FGD6	FYVE, RhoGEF and PH domain containing 6	Cytoskeleton organization, regulation of cell shape, actin cytoskeleton organization
226806_s_at	NFIA	Nuclear factor I/A	Transcription
226886_at	GFPT1	Glutamine–fructose-6-phosphate transaminase 1	Carbohydrate metabolic process
227434_at	WBSCR17	Williams-Beuren syndrome chromosome region 17	NA
227506_at	SLC16A9	Solute carrier family 16, member 9 (monocarboxylic acid transporter 9)	Transmembrane transport
227948_at	FGD4	FYVE, RhoGEF and PH domain containing 4	Apoptosis, cytoskeleton organization
227988_s_at	VPS13A	Vacuolar protein sorting 13 homolog A (S. cerevisiae)	Transport, protein localization
227990_at	SLU7	SLU7 splicing factor homolog (S. cerevisiae)	RNA splicing, mRNA processing
228045_at	NA	NA	NA
228151_at	NA	NA	NA
228242_at	N4BP2	NEDD4 binding protein 2	NA
228853_at	STYX	Serine/threonine/tyrosine interacting protein	Protein amino acid dephosphorylation
228948_at	EPHA4	EPH receptor A4	Protein amino acid phosphorylation, signal transduction, transmembrane receptor protein tyrosine kinase signaling pathway
229376_at	PROX1	Prospero homeobox 1	Cell fate determination, heart development, positive regulation of cell proliferation, negative regulation of viral genome replication, positive regulation of cell cycle, ventricular cardiac myofibril development, positive regulation of sarcomere organization, ventricular septum morphogenesis, positive regulation of heart growth
230831_at	FRMD5	FERM domain containing 5	NA
230894_s_at	NA	NA	NA
232429_at	NA	NA	NA
232871_at	NA	NA	NA
232892_at	C20orf166	Chromosome 20 open reading frame 166	NA
232909_s_at	BPTF	Bromodomain PHD finger transcription factor	Positive regulation of gene-specific transcription
235072_s_at	NA	NA	NA
235855_at	COX15	COX15 homolog, cytochrome c oxidase assembly protein (yeast)	Mitochondrial electron transport, cytochrome c to oxygen
236241_at	MED31	Mediator complex subunit 31	Transcription, regulation of transcription
236379_at	NA	NA	NA
236384_at	NA	NA	NA
238174_at	NA	NA	NA
238768_at	C2orf68	Chromosome 2 open reading frame 68	NA
240044_x_at	TNRC6B	Trinucleotide repeat containing 6B	Regulation of translation
241734_at	SRFBP1	Serum response factor binding protein 1	Regulation of transcription
243630_at	NDUFB1	NADH dehydrogenase (ubiquinone) 1 beta subcomplex, 1, 7kDa	Mitochondrial electron transport, NADH to ubiquinone
37577_at	ARHGAP19	Rho GTPase activating protein 19	Signal transduction

DNA = deoxyribonucleic acid; MAPK = mitogen-activated protein kinase; mRNA = messenger ribonucleic acid; NADH = nicotinamide adenine dinucleotide; rDNA = ribosomal deoxyribonucleic acid; RNA = ribonucleic acid.

Figure 3

Heat Map of Samples From All Patients With IDCM

When we applied this unsupervised clustering algorithm, samples from patients on beta-blocking agents were grouped in a classification tree (highlighted in a dotted square) that was very distinct from patients who received alternative therapy. Each column represents a patient sample, and each row corresponds to a gene. Down-regulated genes are depicted in red, whereas up-regulated genes are labeled in blue.

Table 4

Overexpressed Pathways in Patients Treated With Beta-Blocking Agents Versus Alternative Therapy

Network	GO Processes	p Value	z Score	g Score
PROX1, SMARCA5, PSF, LBP-1B, PCNP	Cellular metabolic process (85.7%; 1.139e-06)	1.08E-32	59.84	59.84
FALZ, TBCD4, NFIA, IRP2, LAMP2	Cellular iron ion homeostasis (8.9%; 4.641e-06)	1.02E-24	46.82	46.82
Ephrin-A receptor 4, PP5, Fodrin (spectrin), YES, FMR1	Developmental process (66.0%; 5.674e-09), phosphate metabolic process (34.0%; 3.201e-08), phosphorus metabolic process (34.0%; 3.201e-08)	3.67E-22	42.46	42.46
SMARCA5, PI3K cat class IA (p110-beta), Cul4/DDB/ROC1 E3 ligase, BUB3, KLF7	Chromatin assembly or disassembly (14.0%; 2.435e-06)	1.01E-17	35.8	35.8
VMAT2, KLF7, PSF, Cul4/DDB/ROC1 E3 ligase, CaMK II	Regulation of growth (31.0%; 2.324e-11), developmental process (76.2%; 2.531e-11) regulation of cell growth (26.2%; 7.286e-11)	1.84E-17	34.7	34.7
TBLR1 (DC42), FMR1, LBP-1B, beta-TrCP, P53DINP1a	Positive regulation of apoptosis (18.6%; 6.013e-05)	3.34E-15	30.81	30.81
Ephrin-A receptor 4, KLF7, FMR1, Frabin, Fodrin (spectrin)	Cellular developmental process (63.8%; 1.245e-15), neuron differentiation (36.2%; 5.415e-13)	5.89E-09	19.57	29.57
LBP-1B, LEDGF/p75, Y549 (GRIF1), TRAP/SMCC complex, LEDGF/p52	DNA metabolic process (23.7%; 3.963e-06), cellular biosynthetic process (57.9%; 6.829e-06)	2.36E-13	28.41	28.41
PROX1, PP5, PSF, LBP-1B, CaMK II delta	Negative regulation of cellular process (52.0%; 9.319e-12)	5.17E-13	26.93	26.93
YES, Cul1/Rbx1 E3 ligase, beta-TrCP, UBXD2, CKAP2	Interspecies interaction between organisms (22.9%; 6.687e-09), T cell activation (16.7%; 4.271e-08)	6.75E-11	23.04	23.04
Cul1/Rbx1 E3 ligase, TBLR1 (DC42), TCF8, beta-TrCP, Cullin 4B	Negative regulation of biological process (62.0%; 1.065e-15), organ development (62.0%; 1.615e-15)	7.28E-09	19.16	19.16
YES, CaMK II, Ephrin-A receptor 4, PI3K cat class IA (p110-beta), SRFBP1	Protein amino acid phosphorylation (45.8%; 1.541e-16)	7.28E-09	19.16	19.16
CaMK II delta, PP5, CaMK II, PSF, SP1	Response to stimulus (70.2%; 1.014e-10)	5.80E-07	15.44	15.44
VMAT2, H('+) = H('+), H('+), Serotonin + H('+) = Serotonin + H('+), H('+) cytosol	Synaptic vesicle amine transport (100.0%; 6.095e-05), monoamine transport (100.0%; 1.097e-03), response to amphetamine (100.0%; 1.463e-03)	6.63E-03	12.19	12.19
AL4A1, CaMK II, p53, FAK1, SHP-2	Positive regulation of cellular process (71.4%; 3.026e-15),	1.25E-03	8.57	8.57

Abbreviations as in Table 3.

Molecular signature to predict responsiveness to beta-blocking agents in patients with new-onset HF

MiPP classification software was used to identify a transcriptomic biomarker containing the minimal set of genes necessary to distinguish patients responsive to beta-blockers from nonresponsive patients (22). A set of 4 transcripts was identified, which predicted therapy responsiveness with very high accuracy (mean error: 0.04, mean sMiPP: 0.89). These included transcripts of the genes ARHGEF1, ALP1, ZNF404, and the transcript 239497_at (Affymetrix ID, not yet characterized). The 6 most robust molecular signatures predictive for response to beta-blockers are illustrated in Table 5. Among those 6 signatures, in particular, 2 were highly accurate, with a mean error of 0.04 and sMiPP ranging from 0.89 to 0.91. Detailed information on biological and molecular function of individual genes within those 2 molecular signatures is illustrated in Tables 6 and 7.

Table 5

Results From Misclassified Penalized Posteriors Classification Illustrating the Most Robust 6 Sets of Genes to Classify Patients Who Were on Beta-Blocker Therapy and Had Good Prognostic Outcome

Split	G1	G2	G3	G4	G5	G6	G7	mean ER	mean sMiPP
31	ALP1	ARFGEF1	ZNF404	239497_at	NA	NA	NA	0.04	0.89
46	C8orf47	1558458_at	PSEN2	TNFRSF14	ITGA2	NA	NA	0.04	0.91
44	SFRS16	MCF2	NA	NA	NA	NA	NA	0.07	0.80
7	SEMA3B	238953_at	PDE1A	NA	NA	NA	NA	0.08	0.82
11	SMOX	GP6PC3	ACSS1	GXYLT1	NA	NA	NA	0.08	0.79
41	231275_at	1565830_at	1564240_at	KCTD5	NA	NA	NA	0.08	0.82

ER = error; MiPP = Misclassified Penalized Posteriors.

Table 6

Detailed Information About Genes Identified in Split 31

Probe Set ID	Gene Symbol	Gene Title	Go Biological Process Term	Go Molecular Function Term
203055_s_at	ARHGEF1	Rho guanine nucleotide exchange factor (GEF) 1	Rho protein signal transduction, cell proliferation	Rho guanyl-nucleotide exchange factor activity, GTPase activator activity, protein binding
211618_s_at	ALPI	Alkaline phosphatase, intestinal	Metabolic process, phosphorylation	Magnesium ion binding, catalytic activity, alkaline phosphatase activity, protein binding, zinc ion binding
239043_at	ZNF404	Zinc finger protein 404	Transcription, regulation of transcription, DNA-dependent	Nucleic acid binding, DNA binding, zinc ion binding, metal ion binding
239497_at	NA	NA	NA	NA

Abbreviations as in Table 3.

Table 7

Detailed Information About Genes Identified in Split 46

Probe Set ID	Gene Symbol	Gene Title	Go Biological Process Term	Go Molecular Function Term	Go Cellular Component Term
1552389_at	C8orf47	Chromosome 8 open reading frame 47	NA	NA	NA
1558458_at	LOC401320	Homo sapiens, clone IMAGE: 4860560, mRNA	NA	NA	NA
204262_s_at	PSEN2	Presenilin 2 (Alzheimer disease 4)	Cell death, apoptotic program, protein processing, amyloid precursor protein catabolic process, positive regulation of catalytic activity	Protein binding, peptidase activity, hydrolase activity	Golgi membrane, kinetochore, integral to nuclear inner membrane, endoplasmic reticulum
209354_at	TNFRSF14	Tumor necrosis factor receptor superfamily, member 14 (herpesvirus entry mediator)	Apoptosis, immune response, signal transduction, cell surface receptor linked signal transduction	Receptor activity, transmembrane receptor activity, tumor necrosis factor receptor activity, protein binding	Cytoplasm, integral to plasma membrane
227314_at	ITGA2	Integrin, alpha 2 (CD49B, alpha 2 subunit of VLA-2 receptor)	Cell adhesion, cell-matrix adhesion, integrin-mediated signaling pathway, organ morphogenesis	Magnesium ion binding, receptor activity, calcium ion binding, collagen binding	Plasma membrane, integrin complex, integrin complex, integral to membrane

Abbreviations as in Table 3.

Among the genes of therapy responsiveness was Rho Guanine Nucleotide Exchange Factor, which is a key player in antiapoptotic function and in a metabolic pathway that regulates contraction and actin-myosin interaction. The pathway is illustrated in Figure 4.

Figure 4

Signaling Pathway of Rho Guanine Nucleotide Exchange Factor

Important functions of this pathway involve antiapoptotic actions and regulation of actin-myosin interaction.

Evaluation of recently suggested candidate genes affected by beta-blocking agents, which may be causative for improved heart function

Finally we tested a set of genes that had been recently related to improvement of EF during beta-blocker therapy (10). Lowes et al. (10) found MYH6 and SERCA2/ATP2A2ase to be overexpressed and MYH7 to be down-regulated in patients with an improvement of EF by at least 5% while on beta-blocking agents. When we evaluated those candidate genes in our population treated with beta-blocking agents, we similarly found overexpression of MYH6 (FC 1.3, q value = 0%) and SERCA2/ATPA2 (FC 1.4, q value = 0%) in patients on beta-blockers and good prognostic outcome versus poor outcome. However, there was also slight overexpression of MYH7 (FC 1.2, q value = 0%) in patients with good prognostic outcome.

Discussion

The major new finding of this study is the discovery of a molecular signature that identified patients with favorable outcomes receiving beta-adrenergic antagonist therapy for new-onset HF. We hypothesize that this signature can be used as a clinical biomarker. Given the mentioned limitations of this study, further prospective trials are necessary to test this hypothesis, in particular if this molecular signature exists in patients prior to initiation of beta-blocker therapy. We have previously employed transcriptomic-based biomarkers (TBBs) to establish diagnosis and prognosis in patients with new-onset HF 15, 16, 18. The present findings show the utility of this approach for pharmacogenomic personalized medicine, and as such, support the use of TBBs for advancing precision medicine.

The success of beta-blocking therapy in the treatment of HF has been attributed to molecular alterations of pathways in the heart, including reverse cardiac remodeling, as well as enhancement of contractile and calcium-channel proteins. Lowes et al. (10) discovered that treatment with beta-blockers leads to activation of a “fetal” gene program encoding SERCA and the α- and β-MYH. The goal of the present study was to use microarray analysis as a comprehensive approach to discover molecular pathways that are affected by beta-blocking agents and to identify transcriptomic markers of therapy responsiveness.

Among 94 overexpressed transcripts in patients taking beta-blocking agents, several genes were involved in apoptosis, such as cytoskeleton associated protein (CKAP)-2, serine incorporator (SERINC)-3, tumor protein p53 inducible nuclear protein (TP53INP)-1 and FYVE, RhoGEF and PH domain containing (FGD)-4. The effect on programmed cell death may lead to reduction of necrosis and inflammation in the myocardium and therefore contribute to the observed reduction in adverse remodeling in this drug category. Furthermore, those patients overexpressed oxidation resistance (OXR)-1, a response gene to oxidative stress, which may provide additional protection from tissue damage. Prospero homeobox (PROX)-1, also up-regulated during beta-blocker therapy, is a gene involved in cell fate determination, proliferation, and heart development, including ventricular cardiac myofibril development and positive sarcomere organization. Similar to PROX-1, other genes involved in heart development have been found to be overexpressed during therapy with beta-blocking agents, as mentioned in the previous text (10). In addition, upstream binding protein (UBP)-1, a gene involved in angiogenesis, was overexpressed in patients taking beta-blocking agents. Surprisingly, no down-regulated genes were identified. This is in disagreement with prior studies 10, 11, 12 and may be due to small sample size, which can lead to identification of only the very most significant gene expression changes with the robust statistical algorithm that was used in this study. Furthermore, prior investigated cohorts are slightly different from ours. Lowes et al. (10) as well as Yasumura et al. (11) investigated samples from patients with chronic HF, whereas our inclusion criteria was new-onset HF. Furthermore, Hamdani et al. (12) investigated biopsies of the left heart, whereas our study analyzed gene expression in the right heart. Asp et al. (23) had previously shown that the gene expression profile of the right atrium differs by about 2% of the genes covered on a microarray from the transcriptome of the left ventricle. Animal studies showed similar findings for right versus left ventricle 24, 25, 26.

When we performed MetaCore pathway analysis, we identified 15 molecular pathways that were activated during treatment with metoprolol or carvedilol. Pathways that were overexpressed with highest significance values were cellular metabolic processes involving phosphate metabolism, cellular iron homeostasis, and developmental processes.

We then went on to investigate if transcriptomic biomarkers derived from heart biopsies can be used to predict therapeutic responsiveness in patients with HF. After we applied MiPP, 2 very robust transcriptomic markers were identified. The first biomarker contained 4 transcripts, including Rho guanine nucleotide exchange factor (GEF)-1, a protein that plays a key role in actin-myosin interaction during contraction. It could be speculated that through the effect on actin-myosin coupling, beta-blocking agents improved contractile efficiency in those patients. Other biomarkers were alkaline phosphatase (ALPI), a gene involved in magnesium binding and phosphorylation; zinc finger protein (ZNF)-404; and the transcript 239497_at. The second biomarker consisted of 5 transcripts, including chromosome 8 open reading frame 47 (C8orf47), presenilin (PSEN)-2, tumor necrosis factor receptor superfamily, member (TNFRSF)-14, integrin alpha (ITGA)-2, and LOC 401320. Importantly, TNFRSF-14 has recently been shown by our group to be overexpressed in patients with good prognostic outcome (15). Therefore, its overexpression in patients on beta-blocking agents with good prognostic outcome may not solely be related to better response to therapy, but may confirm TNFRSF-14 as a generalizable prognostic biomarker independent of therapy.

Study limitations

A limitation of this study is that in order to develop a marker of therapy responsiveness, it was necessary to divide patients on beta-blocking agents into groups of poor versus good outcomes, which ultimately led to small sample size. Our samples were collected at 1 of the most active heart biopsy centers in the United States over a course of 10 years. Despite having collected a large biorepository of 350 EMB samples over a decade, only 18 patients of whom samples were taken had poor prognostic outcome. Of these 18 patients, only 13 were treated with beta-blockers, whereas 5 had contraindications. Therefore, significance values could be overestimated in this subgroup analysis. For the same reasons, the biomarker could not be validated in an additional cohort. Despite this limitation, this is the most comprehensive study to date investigating the effects of beta-blocking agents on the gene expression level and the first, to the best of our knowledge, that has identified a transcriptomic biomarker of therapy responsiveness in this drug class. Mechanistic animal studies will be necessary in the future to evaluate which of the changes during therapy were causative.

Also it warrants mention that it is difficult to establish a clear cut-off for each individual gene to determine which FC in gene expression has significant effect on its function downstream, as some genes are regulated more tightly than others. This issue has been previously discussed in the published data 15, 16, 27. Because part of the purpose of this study was gene discovery, we chose an inclusive approach with an FC cutoff of 1.2 and q value <5% to define significance.

Finally, it should be mentioned that some of the molecular differences that were observed in patients treated with beta-blocking agents versus alternative therapy may reflect changes secondary to comorbidities that prevented patients from being on beta-blockers. For example, chronic obstructive pulmonary disease, a comorbidity that limits the use of beta-blocking agents in patients with HF, may have led to right ventricular strain and consequently changes on the molecular level.

Conclusions

We are presenting results of the first comprehensive transcriptomic analysis that investigates the molecular effects of beta-blocking agents in patients with new-onset HF. Our data suggest a selection of genes that may be involved in better outcomes of patients with HF who are treated with beta-blockers. These findings have implications for the use of TBBs in pharmacogenomic drug development and precision medicine, which is currently a major unmet need in cardiovascular medicine.

COMPETENCY IN MEDICAL KNOWLEDGE: In the era of emerging precision medicine, it is essential to develop techniques that can help the physician to not only prescribe the correct medications for a given patient, but also to monitor and predict patient responsiveness. This is particularly important in disorders such as new onset heart failure due to cardiomyopathy, where there are a growing number of different classes of drugs from which to choose. One of the most successful classes of drugs employed for heart failure are the beta-adrenergic antagonists. In this paper, our group employed a comprehensive transcriptomic-based approach to predict a favorable clinical response to the use of the beta-adrenergic blockers. Transcriptomics employs a comprehensive micro-array that can measure the expression level of known transcripts. With this approach and a dataset of known outcomes in a group of patients with new-onset heart failure, we determined a gene signature that predicted a favorable outcome in patients taking beta-blocking agents. Using various bioinformatic classification algorithms, we developed the most parsimonious signature comprising the expression levels of 4 genes. These genes are overexpressed in patients with a good prognosis relative to patients who had a poor prognosis, defined as death, requirement for left ventricular assist device, or heart transplantation within 2 years. Thus, this transcriptomic-based biomarker can be potentially used by practitioners caring for these patients to help monitor patient's response to drugs. This approach is similar to other types of precision medicine approaches that have been developed to predict clinical trajectory independent of drug treatment or to enhance the diagnostic accuracy in myocarditis.

TRANSLATIONAL OUTLOOK: The transcriptomic-based biomarker described in this study can be potentially used by practitioners caring for these patients to help monitor patient's response to drugs. This approach is similar to other types of precision medicine approaches that have been developed to predict clinical trajectory independent of drug treatment or to enhance the diagnostic accuracy in myocarditis.