An extensive β1-adrenergic receptor gene signaling network regulates molecular remodeling in dilated cardiomyopathies

We investigated the extent, biologic characterization, phenotypic specificity, and possible regulation of a β1-adrenergic receptor–linked (β1-AR–linked) gene signaling network (β1-GSN) involved in left ventricular (LV) eccentric pathologic remodeling. A 430-member β1-GSN was identified by mRNA expression in transgenic mice overexpressing human β1-ARs or from literature curation, which exhibited opposite directional behavior in interventricular septum endomyocardial biopsies taken from patients with beta-blocker–treated, reverse remodeled dilated cardiomyopathies. With reverse remodeling, the major biologic categories and percentage of the dominant directional change were as follows: metabolic (19.3%, 81% upregulated); gene regulation (14.9%, 78% upregulated); extracellular matrix/fibrosis (9.1%, 92% downregulated); and cell homeostasis (13.3%, 60% upregulated). Regarding the comparison of β1-GSN categories with expression from 19,243 nonnetwork genes, phenotypic selection for major β1-GSN categories was exhibited for LV end systolic volume (contractility measure), ejection fraction (remodeling index), and pulmonary wedge pressure (wall tension surrogate), beginning at 3 months and persisting to study completion at 12 months. In addition, 121 lncRNAs were identified as possibly involved in cis-acting regulation of β1-GSN members. We conclude that an extensive 430-member gene network downstream from the β1-AR is involved in pathologic ventricular remodeling, with metabolic genes as the most prevalent category.


Introduction
In patients with heart failure who have a reduced left ventricular (LV) ejection fraction (HFrEF), competitive inhibition of β 1 -adrenergic receptors (β 1 -ARs) by beta-blockers (BBs) produces the largest reduction in mortality among cornerstone drug therapies, thought largely to be due to favorable effects on pathologic ventricular remodeling (1,2).The 2-to 3-month time course for detection of reverse remodeling (RR) indicates a biologic mechanism, presumably mediated by alterations in gene expression (1)(2)(3).Given the primary mechanism of action of BBs is to competitively inhibit β 1 -ARs, it is likely that a substantial portion of the gene expression changes associated with RR are the result of inhibition of chronic hyperadrenergic signaling with β 1 -ARs acting as a key node.This hypothesis is supported by work in the intact failing human heart (1)(2)(3), transgenic (Tg) mice overexpressing human β 1 -ARs (4,5), cultured cardiac myocyte exposure to β-agonists (6), and isoproterenol infusions in murine models (7).However, BB-associated RR in dilated cardiomyopathies causes numerous gene expression changes (8,9) that are unlikely to be primarily related to β 1 -AR signaling, considering the relatively small number of candidate genes identified as β 1 -AR regulated (3).Furthermore, the time course of gene expression changes and myocardial phenotypic associations has not been previously investigated.In addition, how such a β 1 -AR-driven gene network is regulated or coordinated has not been established.
We investigated the extent, biologic characterization, phenotypic specificity, and possible regulation of a β 1 -adrenergic receptor-linked (β 1 -AR-linked) gene signaling network (β 1 -GSN) involved in left ventricular (LV) eccentric pathologic remodeling.A 430-member β 1 -GSN was identified by mRNA expression in transgenic mice overexpressing human β 1 -ARs or from literature curation, which exhibited opposite directional behavior in interventricular septum endomyocardial biopsies taken from patients with beta-blocker-treated, reverse remodeled dilated cardiomyopathies.With reverse remodeling, the major biologic categories and percentage of the dominant directional change were as follows: metabolic (19.3%, 81% upregulated); gene regulation (14.9%, 78% upregulated); extracellular matrix/fibrosis (9.1%, 92% downregulated); and cell homeostasis (13.3%, 60% upregulated).Regarding the comparison of β 1 -GSN categories with expression from 19,243 nonnetwork genes, phenotypic selection for major β 1 -GSN categories was exhibited for LV end systolic volume (contractility measure), ejection fraction (remodeling index), and pulmonary wedge pressure (wall tension surrogate), beginning at 3 months and persisting to study completion at 12 months.In addition, 121 lncRNAs were identified as possibly involved in cisacting regulation of β 1 -GSN members.We conclude that an extensive 430-member gene network downstream from the β 1 -AR is involved in pathologic ventricular remodeling, with metabolic genes as the most prevalent category.
In RR responder interventricular septum biopsied from patients nonischemic dilated cardiomyopathy (NDC), there were 6,909 genes with responder versus nonresponder changes in expression, to which Tg mouse mRNA expression was compared.In TG1, there were 225 Tg genes with expression changes opposite to RR responder versus nonresponder changes, or 14.3% of the Tg total changes (Supplemental Table 2A).In human RR LVs, 61 of these genes were downregulated while upregulated in TG1 mice, and 164 were upregulated coupled with downregulation in TG1 (Supplemental Table 2A).For TG2, a total of 243 genes were antithetically changed in RR, or 29.0% of the 837 total changes (Supplemental Table 2A).Of these, 63 (14 downregulated and 49 upregulated) were also present in TG1.Therefore, for pharmacologically concordant changes, 28.0% of TG1's and 25.9% of TG2's gene expression changes were shared with the counterpart Tg mouse group.The total number of individual genes with pharmacologically concordant significant expression changes across all Tg mouse models and RR was 405 (176 upregulated and 229 downregulated, Supplemental Table 2A).

Literature curation
There were an additional 25 curated genes with expression changes opposite to those in RR human hearts that were not identified as changed under ≥ 2 experimental conditions in Tg mice (Supplemental Table 2B and Supplemental Table 4).All the curated genes were present on the mouse array, and 17 were identified by same-direction changes in at least 2 independent studies, while 8 were through the combination of 1 published study and a single change in the Tg β1-AR overexpressors (Supplemental Table 4).These 14 downregulated and 11 upregulated genes were then added to the 405 from the Tg mouse data, resulting in a total β1-GSN population of 430 genes -190 downregulated and 240 upregulated in RR LVs (Supplemental Table 2C).

Gene expression changes from baseline to last observation in patients with dilated cardiomyopathy
The number of genes with P < 0.05 mRNA abundance changes from baseline in NDC RR responders versus nonresponders at last observation carried forward (LOCF, n = 39 at 12 months, 8 at 3 months) is given in Supplemental Table 5.With the 3 platforms used to measure mRNA abundance in 2 different cohorts, by at least 1 platform measurement, there were 2,975 unique genes upregulated and 3,934 downregulated.Concordance between at least 2 platforms was 6.5% for upregulated and 8.2% for downregulated genes.Compared with the total number of changed expression genes, concordance was markedly higher in the 430 β 1 -GSN genes, 25% in the 240 upregulated genes (P <0.0001), and 35% in the 190 downregulated genes (P < 0.0001) (Supplemental Table 5).For the Any Platform (RT-PCR, microarray, or RNA-Seq) analysis, downregulated genes outnumbered upregulated (P = 0.016), but upregulated were more numerous than downregulated (P = 0.025) in the β 1 -GSN (Supplemental Table 5).

Gene ontology using general classification tools
To investigate biologic function, we performed a gene ontology analysis on the 430 β 1 -GSN RR genes (Figures 1 and 2), including separate analyses with the 190 downregulated (Supplemental Methods, Supplemental Figures 1 and 2) and 240 upregulated genes (Supplemental Figures 3 and 4).The -log scores of the enrichments of the β 1 -GSN genes versus the entire human genome were displayed as waterfall plots for all pathways reaching a P < 0.05 after correcting for multiple comparisons.All 3 analyses (whole network,

R E S E A R C H A R T I C L E
JCI Insight 2023;8(16):e169720 https://doi.org/10.1172/jci.insight.169720downregulated, and upregulated genes) revealed significant enrichments in some known cardiac pathologies of contractile, metabolic, immune, and extracellular matrix (ECM) pathways.The mechanistic diversity is easily visualized by a Reactome pathway map (Supplemental Figure 5).Despite these findings, the canonical and GO databases only cover 249 and 264 β 1 -GSN-responsive genes, respectively, and when combined, only 316 of the 430 total β 1 -GSN genes.

Human ventricular myocardial ontology
With approximately 27% of β 1 -GSN-responsive genes lacking annotation and much of the network lacking a cardiac-specific assignment, we utilized a biologic classification developed specifically for pathologic eccentric ventricular remodeling in the human left ventricle (9).With this ventricular myocardial ontology (VMO) framework, 96% of the β 1 -GSN-responsive genes were assigned to specific biologic categories (Table 2).

Cardiac myocyte contractile or pump function
One of the 2 components of pathologic eccentric remodeling that directly characterizes HFrEF is systolic contractile dysfunction (2).There were multiple mRNA expression changes of β 1 -GSN genes whose encoded proteins affect cardiac myocyte or chamber contractile function.In this group, there was no difference in upregulated versus downregulated genes (Table 2).However, changes in the apoptosis and contractile and associated protein categories would be expected to improve chamber contractile function (Table 2, Table 3, Supplemental Table 6, and Supplemental Figure 6A).
Ca 2+ handling.The SR calcium ATPase Serca 2a (ATP2A2) and the gene for its major regulator, phospholamban (PLN), as well as the ryanodine release channel (RyR2) were upregulated with RR, as previously described (3, 9) (Table 3, Supplemental Table 6, and Supplemental Figure 6A).Also upregulated was S100A1, whose encoded protein can affect Ca 2+ handling and is a positive regulator of contractility (10), and ASPH, whose encoded product is important for SR Ca 2+ homeostasis (11).Downregulated genes included SLC8A1 -which encodes the Na + /Ca 2+ exchanger (NCX) gene that is upregulated in the failing heart and in response to isoproterenol treatment (12) -and calnexin (CANX), encoding an ER associated molecular chaperone that is upregulated in the failing heart (13).Also downregulated were the gene for the Ca 2+ binding protein S100A11, and SYT12, whose encoded protein is involved in Ca 2+ mediation of neurotransmitter release (Table 3, Supplemental Table 6, and Supplemental Figure 6A).
Contractile and associated proteins.As previously reported (3,6,9,16), MYH6, which encodes the higher ATPase activity isoform of myosin heavy chain, meets criteria as a β 1 -GSN member and was upregulated in RR (Table 3, Supplemental Table 6, and Supplemental Figure 6A).MYL3, essential ventricular light chain 1, encoding a protein associated with an increased LV contractile state in isolated rat hearts (17), was similarly upregulated with RR.Also upregulated was myosin light chain kinase 4 (MYLK4), encoding a Ca 2+ /calmodulin-independent light chain kinase that phosphorylates ventricular regulatory light chain resulting in an increase in contractility (18), the striated muscle ubiquitin ligase gene TRIM63, and the muscle stress transducer gene LRRC39.Only 1 contractile protein gene was downregulated, the skeletal muscle isoform of α-actin (ACTA1).This isoform is expressed in fetal development in rodents and is upregulated in the failing human heart to become the majority isoform, but in isolated in vitro motility experiments, it is not associated with contractility differences compared with the cardiac isoform (19) (Table 3 and Supplemental Table 6).
Apoptosis.Seven genes with encoded proteins that can facilitate apoptosis were downregulated with RR (Table 3, Supplemental Table 6, and Supplemental Figure 6A), including CASP3 and DAP.The other 5 (RERE, TM2D2, FILIP1L, MIF, PHLDA3) play various roles in biologic processes that can lead to apoptosis.The only upregulated gene in this category was TRIAP1, which encodes an apoptosis inhibitor (20).

Pathologic chamber hypertrophy
The second major component of eccentric pathologic chamber remodeling is hypertrophy, for which the VMO classification has 4 categories (Table 3 and Figure 6A).In these 4 categories combined, downregulated genes outnumbered upregulated 56 to 13 (P <0.0001) (Table 3 and Supplemental Table 6).
Growth/hypertrophy.With RR, the growth/hypertrophy category exhibited an excess in downregulated genes, 13 versus 4 upregulated (P = 0.029) (Table 3).Notable known mediators of hypertrophy were downregulated, including IGF1, the translation regulator IGF2BP2, EDN1, and the hypertrophy biomarker natriuretic peptides NPPA and NPPB.Also downregulated were MUSTN1, which encodes a myogenesis factor, and 2 genes encoding regulators of cell growth, GRN and GPC1.TMEM43 encodes a protein necessary for normal ventricular contractile protein content, and EXT1 encodes an ER enzyme responsible for chain elongation of heparan sulfate, a cardiac myocyte growth regulator.Two of the downregulated genes, RCAN1 (encoding a calcineurin inhibitor) and BTG1 encode antiproliferative proteins, indicating that downregulation of genes in this category doesn't transfer 1:1 to hypertrophy regression (Table 3 and Supplemental Table 6).
Four genes were upregulated with RR in the growth/hypertrophy category, and 2 of them (MTSS1 and SESN1) encode growth inhibitors (Table 3 and Supplemental Table 6).The other 2 (KY and GAB1) encode proteins involved in normal cell growth responses; in mice, GAB1 gene ablation causes a dilated cardiomyopathy and mitochondrial damage (21).Cytoskeleton.Five genes assigned to the cytoskeleton category were downregulated with RR, including 2 RAC1 regulators (PAK3 and PAK4) whose enzyme gene products link this GTPase to cytoskeletal reorganization and regulation (Table 3 and Supplemental Table 6).In this group, there are also 2 genes (TWF1, DBN1) encoding actin binding proteins and a FHL1 encoding a LIM-only protein expressed in fetal development.The only upregulated gene in the cytoskeleton category was LDB3, which encodes a LIM domain binding protein that interacts with α-actinin-2.
Microtubule function/integrity.Microtubules are assigned a separate subcategory from cytoskeleton because of their known regulation by β 1 -AR signaling and their potential to modify contractile function (22).There were 5 upregulated and 2 downregulated genes (P = 0.26), and this upregulation trend differed from the downregulation in the other 3 Pathologic Chamber Hypertrophy subcategories (2 × 2 χ 2 P = 0.0002; Table 3 and Supplemental Table 6).Of the upregulated genes, 2 (DPYSL2, MAPT) encode proteins that promote microtubule assembly, 2 (FSD2, CEP85) protein products are involved in microtubule organization, and 1 (FYCO1) encoded protein participates in autophagic vesicular transport.Of the 2 downregulated genes, MAP1LC3A encodes a light chain subunit of microtubule-associated protein 1A or 1B that mediates cytoskeletal interactions, and the 1B gene (MAP1B) was itself downregulated.

Metabolism
Mitochondrial based functions.The number of changes in metabolic genes associated with RR is striking, and the majority are in the mitochondrial upregulation subcategory (47 of 67 upregulation and 53 of 83 total changes, both P <0.0001; Table 2, Supplemental Table 7, and Supplemental Figure 6B).Of the 47 upregulated mitochondrial genes, 46 encode proteins involved in either fatty acid transport, utilization,  7).
Fatty acid metabolism, including β-oxidation.There were 12 upregulated genes whose encoded proteins are critically involved in mitochondrial fatty acid β-oxidation, and none that were downregulated (P <0.0001; Supplemental Table 7 and Supplemental Figure 6B).These include genes encoding proteins involved in long-chain fatty acid β-oxidation, CPT1B, ACADVL, HADHA, HADHB, and ACAA2.There is also evidence for increased utilization of medium-chain fatty acids, with upregulation of ACADM.Other upregulated genes encoding enzymes that participate in fatty acid β-oxidation were DECRI and ECI1.Mitochondrial acyl-CoA dehydrogenases that utilize both long-chain fatty acids and branchedchain amino acids -ACAD8, ACAD10, ACADSB, and HADH -were also upregulated in RR.Beyond entry into and effectuation of β-oxidation, there were multiple other upregulated genes encoding proteins that facilitate fatty acid metabolism (ACSL1, ECHDC2, and ECHDC3; Supplemental Table 7).In addition to the 15 β 1 -GSN upregulated genes functioning in fatty acid metabolism, 4 other upregulated genes encode fatty acid transport or enzymatic activity in peroxisomes, and 1 (CD36) fatty acid transport in the sarcolemma (23).In all, 21 upregulated genes encode enzymes (n = 18) or transporters (n = 3) involved in fatty acid metabolism.In contrast, only 1 downregulated gene (FADS1; P < 0.0001 versus upregulated) encodes a protein involved in fatty acid metabolism.
Branched-chain amino acid metabolism.Four genes encoding enzymes that degrade branched-chain amino acids to form acetyl CoA were upregulated.AUH and MCCC1 are key enzymes in leucine degradation, and HIBADH and BCKDHA are involved in valine and branched-chain keto acid metabolism.There were no downregulated genes involved in metabolism of branched chain amino acids (P = 0.046; Supplemental Table 7).
Tricarboxylic acid cycle.Of the remaining upregulated genes encoding proteins important for mitochondrial function, 4 are enzymes critical for the tricarboxylic acid (TCA) cycle: CS, SUCLA2, IDH2, and IDH3B.Also upregulated was ACAT1, which encodes a protein involved in cholesterol biosynthesis but can also regulate acetyl-CoA entry into the TCA cycle.Additionally, 2 upregulated genes involved in anaplerotic replenishment of succinyl-CoA are included in this group, MUT and MMAB (Supplemental Table 7).Thus, 7 TCA associated β 1 -GSN genes were upregulated with RR, versus none downregulated (P = 0.008).Regulation of mitochondrial respiratory capacity (electron transport chain).Of the 47 mitochondrial genes that were upregulated with RR, 15 encode proteins that reside in the electron transport chain (Supplemental Table 7).NDUFS1, NDUFS2, NDUFS4, NDUFA5, and NDUFC2 are enzymatic subunits of respiratory chain complex I, SDHA and SDHD are 2 of the 4 enzymatic subunits of complex II, and CYC1 encodes a subunit of the cytochrome bc1 complex (complex III).Genes encoding enzymes necessary for coenzyme-Q (ubiquinone) synthesis, including COQ3 and COQ9, were upregulated, as was 1 subunit of complex IV -the COX6C subunit of cytochrome C oxidase.Electron transfer flavoprotein-ubiquinone oxidoreductase complex protein subunits ETFB and ETFDH, which are a critical link between the β-oxidation pathway and the respiratory chain, are also found in the upregulated group, as is NUDT13, which encodes an enzyme believed to be involved in regulation of mitochondrial NAD(P) + /NAD(P)H ratios, energy homeostasis, and redox control (24).Of additional interest is GBAS (also known as NIPSNAP2), which encodes a mitochondria-associated protein suspected to play a role in mitochondrial respiration (25).
Mitochondria-and cytosol-based functions -glycolysis.With RR, genes encoding PDH subunit B (PDHB) and the E2 subunit (DLAT) were upregulated (Supplemental Table 6A, Supplemental Table 7, and Supplemental Figure 6B).However, PDK2 was also upregulated; its encoded protein reversibly inactivates the PDH complex via phosphorylation of the E1 subunit.Further complicating the interpretation of inferred PDH activity is the fact that the gene whose encoded phosphatase dephosphorylates the E1 subunit, PDP1, was upregulated in RR.In addition, expression of PFKM, whose encoded protein is the rate-limiting step in glycolysis, was upregulated in RR while the "platelet" isoform, PFKP, was downregulated.Also downregulated was HK1, primarily expressed in mitochondria but also in the cytosol, whose encoded enzyme along with HK2 catalyze the phosphorylation of intracellularly transported glucose to glucose-6-phosphate (Supplemental Table 6B, Supplemental Table 7, and Supplemental Figure 6B).PGAM1 was another downregulated cytoplasmic glycolysis enzyme gene.PANK4, whose encoded protein can interact with PDK2 (26), is also included in the glycolysis subcategory.Therefore, of the 6 upregulated genes encoding constituents of glycolysis, gene products of 5 were mitochondria based and involved in the terminal phase of pyruvate dehydrogenase activity or regulation.In contrast, of the 3 downregulated glycolysis genes, 2 encoded enzymes localized in the cytosol (PFKP and PGAM1), and 1 was expressed in both mitochondria and cytosol (HK1) (Supplemental Table 6B and Supplemental Table 7).
Peroxisome-based functions.Very long-chain fatty acids are metabolized in the peroxisome, and ABCD2, a gene encoding a critical peroxisomal enzyme responsible for transmembrane import of fatty acids, was upregulated in RR (Supplemental Table 7 and Supplemental Figure 6B).Also upregulated in the peroxisomal compartment was PHYH, whose encoded protein catalyzes the first step in α-oxidation of the branched-chain fatty acid phytanic acid; SCP2, encoding a protein involved in the oxidation of branched chain fatty acids; and ECHI whose encoded protein functions in the auxiliary step of the β-oxidation pathway (Supplemental Table 7 and Supplemental Figure 6B).There were no downregulated genes whose encoded proteins are peroxisome based.
Sarcolemma, sarcoplasmic reticulum, or endoplasmic reticulum genes with metabolism functions.Of potential particular interest is the marked (by 12.6-fold) upregulation of DHRS7C (Supplemental Table 6A, Supplemental Table 7, and Supplemental Figure 6B), which encodes a short-chain dehydrogenase/reductase that is highly expressed in the endoplasmic reticulum (ER) of cardiac myocytes.DHRS7C gene expression has been previously shown to be downregulated by both β-and α 1 -agonists, as well as in HF (27).

Gene regulation, transcription, and translation
This category, which encompasses any process affecting protein gene product formation, was markedly biased toward upregulation in RR, 50 versus 14 downregulated (P <0.0001;Table 2).The gene regulation subcategory membership is listed in Supplemental Tables 6 and 8.
TFs.Among upregulated TFs was CREBZF (Supplemental Table 8), encoding a protein that binds cooperatively with ATF4 to transactivate cAMP response elements when activated by PKA phosphorylation (28).Cardiac overexpression of a dominant negative CREB/ATF bZIP construct produces dilated cardiomyopathy in mice (29).The downregulated gene NR4A1 encodes an orphan nuclear receptor TF that is inducible by isoproterenol and regulates the expression of multiple cardiac metabolic genes (30).Also downregulated was the nuclear receptor coactivator PNRC1, whose encoded protein facilitates activation of multiple nuclear hormone receptors (31).
DNA repair, stability and synthesis.The DNA repair, stability and synthesis category is biased toward upregulated genes, 6 versus 1 downregulated (P = 0.059, Supplemental Table 8) and included 2 whose encoded proteins target mitochondrial DNA (REVL3 and RRM2B).
Chromatin/histones.Of potential importance for gene network coordination, the bromodomain and extraterminal (BET) family protein BRD4, a nodal effector of pathologic myocardial remodeling (32) via superenhancer binding and chromatin regulation, was downregulated (Supplemental Table 8).

Other changes within VMO categories
None of the remaining categories in Table 2 exhibited statistically significant differences in up-or downregulation by 1 × 2 χ 2 test or had component changes that are likely directly related to pathologic remodeling.These categories and their gene membership are described in Supplemental Methods, Supplemental Table 6, Supplemental Table 9B, and Supplemental Table 10.
Relationships of β 1 -GSN gene expression to phenotypic characteristics.Microarray global gene expression measurements were available at months 0, 3 (46 subjects each), and 12 (39 subjects), enabling a temporal assessment of β 1 -GSN VMO category gene expression to RR phenotypic measures as overviewed in Supplemental  1 and Supplemental Table 11).In responders as well as the entire cohort (EC), ESV, LVEF, and RVEF were favorably changed at 3 months and then progressively improved at 12 months, with highly statistically significant linear trend tests.EDV exhibited the same pattern.Reduction in PWP reached statistical significance in responders at 12 months with a significant linear trend, but RAP was not changed.NYHA class exhibited a marked improvement at 12 months, with 31% of responders and 26% of the EC becoming asymptomatic, plus a shift from NYHA class III to II at 3 months in both groups.NE was not significantly changed in responders or the EC.HR exhibited the expected reduction from β-blockade at both 3 and 12 months in responders and the EC but, interestingly, was unchanged in nonresponders who treated with the same target doses of β-blockade (3).SBP exhibited a paradoxical small increase at both 3 and 12 months in responders and the EC, which is likely due to improved systolic function.BMI increases slightly at both 3 and 12 months in responders and the EC, while CrCl was unchanged in all groups.In contrast, in nonresponders, no linear trend tests were positive, and only 2, NYHA at 3 months and BMI at 12 months, were statistically significant compared with month 0.
Relationships between mRNA abundance and the phenotypic measures in Supplemental Table 11 were assessed by nonparametric permutation testing, generating a permuted Z score (Z p ) that reflects the strength of β 1 -GSN VMO category mRNA abundance-phenotypic parameter relationship relative to the null.The key component of the Z p is the within-VMO category averaged Rho value across all the gene members.To illustrate the type of correlations produced, Supplemental Figure 8 is a plot of the average normalized mRNA abundance for metabolism genes versus LVEF or PWP in the 46 study subjects, generating a composite Rho in contrast to the Z p calculation, which is based on the average Rho of individual gene mRNA abundance-phenotype correlations.The characteristics of these plots are discussed further in Supplemental Methods.
Heatmaps for mRNA abundance versus phenotypic measure Z p values and cluster analyses for each of the 21 gene expression categories and 12 phenotypic characteristics for the 3 time points are in Figure 3. Data in Supplemental Table 11 indicate that the study population was composed of 2 populations of ventricular myocardial RR related responses, which are substantially favorable in the 30 responders and absent in the 16 nonresponders.The Z p in the heatmap are the average of Spearman's rank correlation coefficients (Rho) in the β 1 -GSN ontology categories versus each phenotypic measure, standardized using control values from non β 1 -GSN genes.At month 0 (Figure 3A), for the 430 β 1 -GSN member mRNA abundance measurements versus the 12 phenotypic parameters, of the 252 possible VMO category-phenotype measures there were 152 that had a positive Z p value, and 100 that were negative (P = 0.001).For Z p ≥ or ≤1.96 the numbers were 40 and 5, respectively (P <0.0001).At 3 and 12 months (Figure 3, B and C), the positive/negative Z p counts were, respectively, 170/82 (P <0.0001) and 200/52 (P <0.0001), and the Z p ≥ 1.96 and ≤ -1.96 were 58/6 (P <0.0001) and 72/1 (P <0.0001) (Supplemental Table 12, Supplemental Methods).The number of Z p values ≥ 1.96 at all 3 time points is markedly greater than expected by random chance versus 2.5%, at months 0, 3, and 12, values were 40, 58, and 72 versus 3.8, 4.2, and 5.0 (P < 0.0001).In addition, the patterns of gene expression categories with Z p ≥ 1.96 within phenotypic measures are different across the 3 time points (χ 2 P = 0.003; Supplemental Table 12 and Figure 3) with the number of statistically significant Z p generally increasing over time (proportional trend test P = 0.0006; Supplemental Table 13).However, for phenotypic measures, the pattern of statistically significant Z p scores varies, which can be appreciated by inspection of Figure 4 and Supplemental Figure 9, A and B. Only 1 measure, HR, has a greater number of Z p ≥ 1.96 at month 0 compared with months 3 or 12 -9 versus 2 or 3, respectively (P = 0.019 by χ 2 , 0.026 by proportional trend; Supplemental Tables 13 and 14 and Supplemental Figure 9A).LVEF, RVEF, NE, PWP, and BMI (Figure 4 and Supplemental Figure 9, A and B) all exhibit increasing numbers of Z p ≥ 1.96 with time, with proportional trend tests of P < 0.05.In contrast, ESV and EDV have more statistically significant Z p values at 3 months versus either months 0 or 12 (Supplemental Tables 13 and 14, Figure 4, and Supplemental Figure 9, A and B).
The direct or inverse correlation of mRNA abundance within up-or downregulated genes on the LOCF analysis (Table 2) versus phenotypic measure is given in Figure 4 and Supplemental Figure 9, A and B, above the individual gene category bars, depicted as direct relationship (+), inverse relationship (−), or 0 (no statistically significant relationship).In all cases where a significant relationship was observed for both up-and downregulated genes (n = 48 of 170 possible combinations), respective inverse associations with phenotype measures were noted.For example, for PWP and metabolic genes, those that are upregulated on LOCF analysis have an inverse (−) relationship to PWP values, while downregulated genes have a direct (+) association, indicating that the expression of all genes in this category was associated with declining PWP values.For net association including both up-and downregulated genes, the color coding of the bar is inverse, meaning that upregulated genes correlated with declining PWP values were the dominant relationship.JCI Insight 2023;8(16):e169720 https://doi.org/10.1172/jci.insight.169720 Upset plots comparing Z p ≥ 1.96 at months 0, 3, and 12 are given in Figure 5. Here, the temporal course of the gene expression-phenotypic measure relationship is readily observed as is the aggregation of various gene category relationships to individual phenotypic measures.Note that the ordering by number of Z p ≥ 1.96 from most to least is different at each time point and that 3 and 12 months are very different from month 0. Month 12 has greater numbers of aggregated significant correlations than month 3, but month 3 has statistically significantly larger Z p values for ESV and EDV than month 0 or month 12.
Long noncoding RNAs.To investigate the regulation of β 1 -GSN-responsive genes beyond the known input of TFs and downstream cross-regulation, we considered a role for lncRNAs since they are capable of contributing to gene network coordination (34).The RNA-Seq data were aligned to the NONCODE database, and the entire spectrum of identified non-coding RNAs are shown in Figure 6A.Unannotated species are by far the majority (75.6%), with pseudogenes (8.49%) and antisense (7.71%) also constituting larger fractions than annotated lncRNAs (4.01%).A total of 30,861 unique annotated or unannotated lncRNAs were identified in at least 1 subject by NONCODE (Supplemental Table 15).
Figure 6B is a heatmap of annotated and unannotated lncRNAs that changed from baseline to 12 months/LOCF.Using a fold change of 1.5 and a P < 0.05 for a meaningful change, 3,446 lncRNAs changed between baseline and end of study, 2,278 in the 6 analyzed superresponders (R SR ) versus 1,168 in the 4 nonresponders (NR SR ).In R SR 1,118 lncRNAs were P < 0.05 for increased compared with baseline, and 1,160 were decreased.In NR SR , the number of increased lncRNAs was 690, compared with 478 decreased.None of the significantly changed lncRNAs that overlapped with a changed β 1 -GSN mRNA were annotated in GENCODE.
The cis-acting lncRNAs are positioned in close genomic proximity to the genes they regulate.For the changed lncRNAs in Figure 6B, we determined the closest coding gene and then asked if they were enriched in the β 1 -GSN versus the entire number of protein coding genes (PCG) in hg19 GENCODE (Figure 6C).In responders, 121 of the 430 (28%) member β 1 -GSN-responsive genes were the closest PCG to a changed lncRNA, compared with 315 of the total 19,500 PCGs (1.6%) (P = 6.3 × 10 -19 ).Of the 121 β 1 -GSN-responsive genes that were the closest PGCs, 68 were associated with downregulated and 53 with upregulated transcripts, and each set of changed genes was enriched relative to respective decreases and increases in lncRNAs (Figure 6C).In contrast, lncRNAs significantly changed in nonresponders were not enriched in β 1 -GSN genes: all, n = 23, P = 0.48; upregulated, n = 12, P = 0.50; downregulated, n = 11, P = 0.88.A listing of the 121 nearest proximity β 1 -GSN members in responders is in Supplemental Table 16.For the 68 downregulated β 1 -GSN genes associated with changed lncRNAs in responders, 13, 6, 1, and 1 were the closest gene to 2, 3, 4, and 5 different changed lncRNAs, respectively, and the remaining 47 were the closest gene to a single lncRNA (Supplemental Table 16).For the 53 upregulated transcripts, 5 were from the closest genes to each of 2 and 3 lncRNAs and 43 genes were closest to a single lncRNA.

Discussion
We have provided evidence of an extensive, 430-member gene network downstream from the β 1 -AR that participates in pathologic ventricular eccentric remodeling and its reversal.

Global summary of gene expression changes
In addition to standard gene enrichment analyses, we utilized a VMO classification specifically designed for ventricular chamber remodeling of the human heart (9).Similarities between the general canonical pathway and VMO analyses included the metabolism categories such as upregulation in β-oxidation of fatty acids (in agreement with Wiki, Reactome; https://www.gsea-msigdb.org/gsea/msigdb/collection_details.jsp), fatty acid metabolism in general (Hallmark/KEGG), TCA cycle (Wiki/KEGG), branched chain amino acids metabolism (KEGG), and oxidative phosphorylation/electron transport (Hallmark/ Wiki).Additional agreement occurred for downregulation in apoptosis (Hallmark) and fibrosis (Wiki).In addition, for biologic processes, there was approximate agreement with VMO for upregulation in force of contraction, branched chain amino acid metabolism and fatty acid β-oxidation, and downregulation in heart growth and collagen fibril organization.
By the VMO classification, the β 1 -GSN changes in gene expression were widespread across multiple biologic categories and consistent with the reversal of the structural and functional phenotypic characteristics of eccentric pathologic remodeling (2).Changes in mRNA expression occurred in key genes whose encoded proteins would positively influence contractility or chamber pump function (Ca 2+ handling, β-adrenergic signaling, contractile proteins, apoptosis) and would lead to regression of pathologic hypertrophy (downregulation in growth promoting factors, fibrosis/ECM promoting and cytoskeletal genes, upregulation in microtubular function).Some of these changes have previously been identified (3,9), but many have not been definitively connected to a β 1 -AR mechanism.Presumably most of the changes in gene expression in Tg mice and RR human LVs originated from β 1 -AR-expressing cardiac myocytes, but a myocyte origin is not a requirement for β 1 -GSN membership.For example, the marked number of fibrosis/ECM genes exhibiting downregulation in RR is likely mediated through fibroblasts, which do not contain β 1 -ARs, responding to the inhibition of microRNAs (35) and/or growth factor (36) secretion from cardiac myocytes.
Two large categories of changes likely underpinned the favorable changes in gene expression contributing to RR. Metabolism category genes were dramatically biased in favor of upregulation versus downregulation, with genes encoding mitochondrial proteins accounting for the majority of the imbalance.Improvement in depressed contractile function in the dilated, myopathic LV would require an enhanced and more efficient source of energy, which was likely supplied by upregulation in the expression of genes encoding fatty acid metabolism.Critical steps in fatty acid β-oxidation, respiratory chain function, branched chain amino acid metabolism, and metabolism with peroxisomal and membrane compartments were downregulated in β 1 -AR-overexpressing mice that ultimately develop a dilated cardiomyopathy (4,5) and then upregulated on reverse LV remodeling in NDC patients treated with β 1 -AR antagonist treatment.The dramatic upregulation of genes involved in mitochondrial metabolism supports the concept that RR involves normalization of metabolic changes that have been previously described in end-stage HF, a shift from fatty acid oxidation and oxidative metabolism toward glycolysis and ketone utilization in the human heart (37).In pathologic ventricular remodeling neurohormonal signaling and mitochondrial function are generally linked (38), and the β 1 -GSN may be playing a major role in this connectivity.
The other large category of gene expression changes favoring upregulation versus downregulation was the gene regulation, transcription, and translation category.Extensive changes in gene expression would require multiple protein synthesis pathway changes, for which there was evidence at multiple levels, with 50 genes exhibiting upregulation versus 14 downregulated on RR.The upregulated genes included 15 that are involved in transcription regulation and 7 TFs, compared with only 3 and 2 in these respective categories.
Key to any proposed β 1 -GSN construct is the signaling architecture between the β 1 -AR receptor and nuclear gene regulation.β 1 -Adrenergic canonical and noncanonical signaling pathways that were likely involved through cross-regulation were altered in RR.With canonical β 1 -AR signaling, RR-associated mRNA changes were clearly biased in favor of restoration of the downregulated signaling present in Tg mice, with 13 changes in RR favoring increased cAMP/PKA signaling or improved neurotransmission versus 4 changes that would be predicted to decrease signaling.For potential GPCR cross-talk-mediated signaling, which -based  2 and Figure 3) versus non-β 1 -GSN VMO controls at months 0, 3, and 12. Blue, direct relationship of phenotypic measurement with net (including upregulated and downregulated genes) RNA expression; yellow, inverse relationship of net mRNA abundance changes with phenotypic measure; gray, no statistically significant relationship of net mRNA expression and phenotype measurement.The designations above the bars are as follows: first entry is for upregulated genes (top), second is for downregulated (bottom); (+), direct directional correlation with the phenotypic measure; (−), inverse directional correlation; 0, no directional correlation.on cognate gene expression changes -was extensive and involved multiple pathways and signaling effectors, phosphoinositide/phospholipase and non-β-adrenergic neurohormonal signaling exhibited numerically greater numbers of upregulated versus downregulated genes that would be expected to improve signaling.Moreover, the extensive involvement of small GTPases and their regulators could play a crucial role in transducing a single proximal signaling node into widespread changes in gene expression.

Phenotypic association and temporal course
Serial gene expression measurements conducted at baseline, 3, and 12 months allowed for temporal assessment of β 1 -GSN gene category-phenotypic associations during RR.We used a variant of nonparametric permutation testing for these analyses, a valid method for generating a randomly distributed control (39).Compared with null/control mRNA abundance-phenotypic correlations empirically generated by 10,000 permutations using non-β 1 -GSN genes, β 1 -GSN member values were correlated with phenotypic measurements at all time points and well above (by 10-to 14-fold) the statistically expected number of Z p ≥ 1.96 values based on the null.The numbers of statistically significant Z p increased progressively from months 0 to 12, and for ESV and EDV, month 3 had higher Z p scores than either month 0 or month 12. Z p patterns also changed over time; at month 0, HR, LVEF, and RVEF had the most numerous Z p ≥ 1.96 values, while at 3 months, it was ESV, EDV, LVEF, and RVEF, and at month 12, LVEF, PWP, RAP and RVEF had the most statistically significant Z p values.Within categories, the directionality of relationships between gene expression and phenotypic measures was acute fight or flight physiologic responses -namely, when deployed chronically, activation of an extensive gene network that is involved in maladaptive myocardial remodeling.

Methods
Model system identification of genes whose myocardial expression is regulated by β 1 -AR signaling Genes exhibiting mRNA expression changes in Tg mice overexpressing human ADRB1 389 polymorphic variants in the heart.This methodology is described in Supplemental Methods.
Literature manual curation of experiments demonstrating changes in myocardial mRNA or protein expression by subacute or chronic administration of β-agonists or antagonists.Evidence for β 1 -AR signaling was also sought by literature curation as described in Supplemental Methods.As for Tg mouse data, ≥ 2 statistically significant experiments or reports with same-directional change were required to qualify as β 1 -AR regulated.Curated evidence could also qualify with 1 statistically significant study, with a single P < 0.05 condition from the Tg mice that would not otherwise be eligible.

BB-associated RR in nonischemic patients with dilated cardiomyopathy
RR responders and nonresponders were identified in the BB Effect on Remodeling and Gene Expression (BORG, NCT01798992) study as previously reported (9) and further described in Supplemental Methods.Global gene expression was measured as mRNA abundance in RNA extracted from right-sided mid-distal interventricular septum EmBx using microarrays in the 47 patient BORG EC, and RNA-Seq in a 12-patient super-responder cohort (SRC) (9) (Table 1).Prespecified candidate gene mRNA abundance (n = 50) was also measured in EC RNA by reverse transcription PCR (RT-PCR) (Supplemental Methods).Endomyocardial biopsies, right heart catheterization, and SPECT imaging were performed at baseline (time 0), 3 months, and 12 months.

Detection of members of a β 1 -GSN operative in LV RR
We used a 3-tiered algorithm for selection of genes qualifying for membership in a β 1 -GSN consisting of (a) establishing that a gene alters its myocardial expression in response to changes in β 1 -AR input, (b) demonstrating these candidate genes change their expression during RR in response to β 1 -blockade in NDC patients, and (c) providing evidence of multigene network behavior (Supplemental Table 1, Supplemental Methods).
Gene ontology and enrichment analyses.A modified Fisher's exact test was used to determine pathway enrichment, with a Benjamini-Hochberg correction for FDR (48).Pathway databases were downloaded from the Broad Institute's MSigDB (https://gsea-msigdb.org/gsea/msigdb).Myocardial biology gene ontology assignment of members of the β 1 -GSN to a previously described VMO classification tailored to eccentric pathologic remodeling (9) was performed.

Specificity and temporal characteristics of the β 1 -GSN for effects on physiologic measures or clinical consequences of ventricular RR
In order to investigate the specificity of β 1 -GSN VMO categories to physiologic and clinical consequences of ventricular RR, the 21 VMO categories were compared by nonparametric permutation testing to 12 cardiac physiologic, pharmacodynamic, or clinical measures (phenotypic measures) for their relationships at baseline, 3 months, and 12 months.For these analyses, microarray identification of 19,673 transcripts in the 46 members of the entire cohort who had microarray measurements on biopsies taken at all time points were used to compare 430 β 1 -GSN member genes with 19,243 nonmembers.

Identification of lncRNAs potentially involved in β 1 -GSN regulation
In RNA extracted from the 6 SRC R SR and 4 age-and sex-matched NR SR , lncRNA was measured by RNA-Seq.Ribodepleted RNA-Seq reads were aligned via TopHat2 to NONCODE v4 and hg19.The resulting transcripts were then filtered by those that significantly changed in RR responders versus nonresponders (≥1.5 fold change, P < 0.05) using within-responder group analyses as described for mRNAs (Supplemental Methods).lncRNAs were annotated from NONCODE using a minimum of a 99% overlap with the most recent GENCODE annotation from hg19.The cis regulation was inferred computationally using proximity between a given lncRNA and a corresponding mRNA transcript.Overlapping lncRNA and mRNA

R E S E A R C H A R T I C L E
JCI Insight 2023;8(16):e169720 https://doi.org/10.1172/jci.insight.169720transcripts were identified using bedtools.Instances in which both transcripts significantly changed in expression were isolated as examples of potential cis regulation downstream of BB-mediated RR.

Statistics
Statistical path to membership in the β 1 -GSN.Tg mouse or curated source (Supplemental Methods) gene expression changes were matched to P < 0.05 opposite-direction changes for the same gene in RR human ventricular septal biopsy material.Using previously defined methodology (9), a gene's transcript abundance was considered changed in the EC if it was P < 0.05 by RT-PCR or microarray in responders versus nonresponders (R EC versus NR EC , or if the SRC if RNA-Seq data were P < 0.05 in R SR (n = 6) versus NR SR (n = 6).In addition, abundance change was considered statistically significant for a P < 0.05 in a paired analysis of LOCF versus baseline abundance in responders, provided that any same-direction change in nonresponders was P ≥ 0.10 (Supplemental Methods).The estimated net α and FDRs for Tg gene identification combined with human RR responder versus nonresponder statistical validation are given in the Supplemental Methods.
Statistical testing for β 1 -GSN membership.In the Tg ADRB1 389Arg or 389Gly mice the original (4, 5), and the confirmatory statistical tests performed in the current study 1-way ANOVA with a Benjamini-Hochberg false discovery rate of 5% in TG1 mice (4) and 1-way ANOVA with a Bonferroni correction in TG2 mice (5).The multiple comparison testing assumed a 2-tailed distribution, and P < 0.05 was required for statistical significance for both ANOVA and multiple comparisons.For analysis of changes in human myocardial mRNAs, Wilcoxon rank-sum or Wilcoxon signed rank tests were used (3,9).Categorical variables were compared using χ 2 1 × 2 configurations comparing the number of upregulated versus downregulated genes in a single category and 2 × 2 configurations comparing a single category or subcategory to all other categories in the data set.LVEF and LV volume changes were assessed by unpaired t tests.For temporal analysis of RR changes, 1-way ANOVA coupled with a Holm-Šídák test on paired values was used to assess changes from baseline at 3 and 12 months, followed by a linear trend test on unpaired values.
Temporal relationships of β 1 -GSN mRNA expression to phenotype measurements.A variant of nonparametric permutation testing (49,50) (summarized in Supplemental Figure 7 and Supplemental Methods) was employed to compare patterns of β 1 -GSN gene expression relative to cardiac and clinical phenotypic responses before and after RR.In separate analyses of data from months 0, 3, and 12, Spearman's correlation coefficients (Rho values) were generated from individual gene mRNA abundance-phenotype measurements within each β 1 -GSN or null (control) permuted VMO category and averaged, from which a Z p was calculated.The null effect was calculated using randomly selected gene sets obtained from the 19,243 microarray identified transcripts that were not members of the n = 430 β 1 -GSN.Further details of the permutation testing method are given in Supplemental Methods.P < 0.05 was considered significant for all measurements, including for the 3-month-old TG1 mouse mRNA analysis (4) where a Benjamini-Hochberg test for FDR correction at α = 0.05 was applied (Supplemental Table 2 and Supplemental Methods).Analyses were performed using the R statistical suite version (v.4.1) or GraphPad Prism.Heatmaps (9) and Upset plots (51) were generated as previously described.

Study approval
The Tg mice studies were approved by the University of Maryland School of Medicine IACUC (4) (TG1) and the University of Colorado Health Sciences Center Animal Care Committee (5) (TG2).The methods of anesthesia and euthanization for TG1 mice were CO 2 narcosis followed by cervical dislocation, and for TG2, a combination of xylazine and ketamine i.p. followed by cervical dislocation were used.
The clinical study (NCT01798992) was conducted according to Declaration of Helsinki guidelines, the protocol was approved by the University of Colorado and University of Utah IRBs, and all patients signed informed written consent.
Data are shown for nonischemic patients with dilated cardiomyopathy in the BORG trial(3), entire cohort of responders (R EC ) versus nonresponders (NR EC ), and superresponders (R SR ) versus SR subcohort nonresponders (NR SR ).C P < 0.05 versus NR EC or NR SR ; D P < 0.001 versus NR EC or NR SR .A W, White; B, Black; A, Asian; I, Islander or Native American; H, Hispanic.B Compared by Wilcoxon rank sum.R EC , entire cohort responder, LVEF increase by ≥5 absolute % at 3mos, or ≥ 8 absolute percentage at 12 months; NR EC , entire cohort nonresponder, LVEF change not meeting R EC criteria; R SR , responder in SR subcohort, LVEF change ≥ 10 absolute percentage; NR SR , SR subcohort nonresponder, LVEF change not meeting R SR criterion; LV, left ventricular; RV, right ventricular; SBP, systolic blood pressure; PWP, pulmonary wedge pressure; PAP, pulmonary arterial pressure; NYHA, New York Heart functional class; Δ, change from baseline; CrCl, creatinine clearance.LOCF, last observation carried forward (12 months or 3 months in subjects not completing study, ± SEM).Baseline data are ± SD.

Figure 1 .
Figure 1.Canonical pathway enrichment of the β 1 -GSN in reverseremodeled human left ventricles.-Log P values at the end of bars are from a modified Fisher's exact test to determine pathway enrichment, with a Benjamini-Hochberg correction for FDR.n = 47.

Figure 2 .
Figure 2. GO pathway enrichment of the β 1 gene signaling network in reverse-remodeled human left ventricles.-Log P values at the end of bars are from a modified Fisher's exact test to determine pathway enrichment, with a Benjamini-Hochberg correction for FDR.n = 47.

Figure 3 .
Figure 3. Heatmap and cluster analysis of the Z p comparing β 1 -GSN mRNA abundance-phenotype correlations with non-β 1 -GSN controls.A positive or negative Z p measures the degree of correlation that is respectively > or < controls/null, with an absolute Z p of > 1.96 (rounded to 2.0) considered statistically significant.Blue color intensity is the degree that the β 1 -GSN Z p exceeds the expected correlation in controls, red represents correlation less than expected, and white (0.0) is the null.The linear scale is suppressed at absolute values > 5.0.(A) Month 0 (baseline, n = 46), mRNA abundance of the 430 β 1 -GSN genes compared with values of 12 different cardiac and patient phenotypic characteristics in the 30 reverse remodeling responders.(B and C) Month 3 (n = 46) and month 12 (n = 39) are shown.

Figure 4 .
Figure 4. Temporal pattern of mRNA abundance-ventricular function and structure phenotypic relationships.Patients with nonischemic dilated cardiomyopathy were investigated at baseline (month 0; n = 46) on no beta-blocker treatment; n= 46 were investigated at 3 months and n = 39 at 12 months on beta-blockade.The y axes are Z p values ≥ 1.96 from nonparametric permutation testing of average Spearman's rank correlation Rho values of β 1 -GSN ventricular myocardial ontology (VMO) categories (Table2and Figure3) versus non-β 1 -GSN VMO controls at months 0, 3, and 12. Blue, direct relationship of phenotypic measurement with net (including upregulated and downregulated genes) RNA expression; yellow, inverse relationship of net mRNA abundance changes with phenotypic measure; gray, no statistically significant relationship of net mRNA expression and phenotype measurement.The designations above the bars are as follows: first entry is for upregulated genes (top), second is for downregulated (bottom); (+), direct directional correlation with the phenotypic measure; (−), inverse directional correlation; 0, no directional correlation.