Sequential analysis of myocardial gene expression with phenotypic change: Use of cross-platform concordance to strengthen biologic relevance

Abstract

Objectives

To investigate the biologic relevance of cross-platform concordant changes in gene expression in intact human failing/hypertrophied ventricular myocardium undergoing reverse remodeling.

Background

Information is lacking on genes and networks involved in remodeled human LVs, and in the associated investigative best practices.

Methods

We measured mRNA expression in ventricular septal endomyocardial biopsies from 47 idiopathic dilated cardiomyopathy patients, at baseline and after 3–12 months of β-blocker treatment to effect left ventricular (LV) reverse remodeling as measured by ejection fraction (LVEF). Cross-platform gene expression change concordance was investigated in reverse remodeling Responders (R) and Nonresponders (NR) using 3 platforms (RT-qPCR, microarray, and RNA-Seq) and two cohorts (All 47 subjects (A-S) and a 12 patient “Super-Responder” (S-R) subset of A-S).

Results

For 50 prespecified candidate genes, in A-S mRNA expression 2 platform concordance (CcpT), but not single platform change, was directly related to reverse remodeling, indicating CcpT has biologic significance. Candidate genes yielded a CcpT (PCR/microarray) of 62% for Responder vs. Nonresponder (R/NR) change from baseline analysis in A-S, and ranged from 38% to 100% in S-R for PCR/microarray/RNA-Seq 2 platform comparisons. Global gene CcpT measured by microarray/RNA-Seq was less than for candidate genes, in S-R R/NR 17.5% vs. 38% (P = 0.036). For S-R global gene expression changes, both cross-cohort concordance (CccT) and CcpT yielded markedly greater values for an R/NR vs. an R-only analysis (by 22 fold for CccT and 7 fold for CcpT). Pathway analysis of concordant global changes for R/NR in S-R revealed signals for downregulation of multiple phosphoinositide canonical pathways, plus expected evidence of a β1-adrenergic receptor gene network including enhanced Ca2+ signaling.

Conclusions

Two-platform concordant change in candidate gene expression is associated with LV biologic effects, and global expression concordant changes are best identified in an R/NR design that can yield novel information.

About the Authors:

Lee S. Toni

Roles Conceptualization, Data curation, Investigation, Methodology, Project administration, Validation, Writing – original draft, Writing – review & editing

Affiliation: Division of Cardiology, University of Colorado, Denver/Anschutz Medical Campus, Aurora, Colorado, United States of America

Ian A. Carroll

Roles Formal analysis, Investigation, Methodology, Software, Visualization, Writing – review & editing

Affiliations Division of Cardiology, University of Colorado, Denver/Anschutz Medical Campus, Aurora, Colorado, United States of America, ARCA biopharma, Westminster, Colorado, United States of America

Kenneth L. Jones

Roles Formal analysis, Investigation, Validation, Writing – review & editing

Affiliation: Department of Pediatrics, University of Colorado, Denver/Anschutz Medical Campus, Aurora, Colorado, United States of America

Jessica A. Schwisow

Roles Data curation, Investigation

Affiliation: Division of Cardiology, University of Colorado, Denver/Anschutz Medical Campus, Aurora, Colorado, United States of America

Wayne A. Minobe

Roles Data curation, Investigation, Project administration, Supervision, Validation, Writing – review & editing

Affiliations Division of Cardiology, University of Colorado, Denver/Anschutz Medical Campus, Aurora, Colorado, United States of America, University of Colorado Cardiovascular Institute Pharmacogenomics, Boulder and Aurora, Colorado, United States of America

Erin M. Rodriguez

Roles Data curation, Investigation

Affiliation: Division of Cardiology, University of Colorado, Denver/Anschutz Medical Campus, Aurora, Colorado, United States of America

Natasha L. Altman

Roles Investigation, Writing – review & editing

Affiliations Division of Cardiology, University of Colorado, Denver/Anschutz Medical Campus, Aurora, Colorado, United States of America, University of Colorado Cardiovascular Institute Pharmacogenomics, Boulder and Aurora, Colorado, United States of America

ORCID logo http://orcid.org/0000-0002-4190-057X

Brian D. Lowes

Roles Funding acquisition, Investigation

Affiliation: Division of Cardiology, University of Nebraska Medical Center, Omaha, Nebraska, United States of America

Edward M. Gilbert

Roles Funding acquisition, Investigation, Writing – review & editing

Affiliation: Division of Cardiology, University of Utah Medical Center, Salt Lake City, Utah, United States of America

Peter M. Buttrick

Roles Funding acquisition, Resources, Writing – review & editing

Affiliations Division of Cardiology, University of Colorado, Denver/Anschutz Medical Campus, Aurora, Colorado, United States of America, University of Colorado Cardiovascular Institute Pharmacogenomics, Boulder and Aurora, Colorado, United States of America

David P. Kao

Roles Data curation, Formal analysis, Investigation, Project administration, Software, Supervision, Validation, Writing – review & editing

Affiliations Division of Cardiology, University of Colorado, Denver/Anschutz Medical Campus, Aurora, Colorado, United States of America, University of Colorado Cardiovascular Institute Pharmacogenomics, Boulder and Aurora, Colorado, United States of America

Michael R. Bristow

Roles Conceptualization, Data curation, Formal analysis, Funding acquisition, Methodology, Project administration, Resources, Supervision, Visualization, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliations Division of Cardiology, University of Colorado, Denver/Anschutz Medical Campus, Aurora, Colorado, United States of America, ARCA biopharma, Westminster, Colorado, United States of America, University of Colorado Cardiovascular Institute Pharmacogenomics, Boulder and Aurora, Colorado, United States of America

ORCID logo http://orcid.org/0000-0002-8860-4698

Introduction

Heart failure is a consequence of structural and functional impairments of ventricular myocardium that reduce the ability of the heart to adequately supply blood to the body [1]. The current method of classifying heart failure is by whether left ventricular ejection fraction (LVEF) is reduced (“HFrEF”) or relatively preserved (“HFpEF), which generally subdivides heart failure into respective categories that respond to multiple types of drug or device therapies vs. those that do not [2]. Despite the established effectiveness of 8 drug and 2 device classes [3], HFrEF continues to increase in prevalence [4] and major adverse outcomes including mortality have not appreciably declined [5].

Current drug therapy of HFrEF is based on targeting final common pathway maladaptive mechanisms that provide short term compensation to impaired ventricular chamber contractile function, but have adverse effects long term [3]. Despite the relative success of this therapeutic approach, the molecular mechanisms underlying the progression of HFrEF as well as the molecular phenotypic signatures of effective therapies have not been established. Reduced to first principles, disease in any tissue is the result of altered gene expression that produces organ dysfunction, and therefore the causative mechanisms are discoverable by identifying abnormal patterns of gene expression directly related to the dysfunction. In HFrEF we have attempted to identify abnormal gene expression profiles at the mRNA [6–8] or miRNA [9] level by producing therapeutic interventional reversal of the structural and functional abnormalities in HFrEF ventricles, and identifying ventricular septal endomyocardial biopsy gene expression changes associated with dynamic left ventricular reverse remodeling. This approach, which we term “sequential analysis of myocardial gene expression with phenotypic change” (SAMGE-PC), has several potential advantages, including: investigating serial changes in subjects whose genetic and non-interventional therapeutic backgrounds remain constant before and after ventricular reverse remodeling; being able to link changes in gene expression to specific functional and structural changes in ventricular myocardium; having an internal control group consisting of reverse remodeling nonresponders; being able to investigate HFrEF at earlier stages than can be done by using starting material from explanted hearts removed from cardiac transplant recipients; and the inherent statistical advantage of within group repeated measures vs. cross-sectional comparisons [7].

In these studies, for mRNA expression we have used two types of quantitative PCR [6,8,9], two generations of microarrays [7–9] and most recently, RNA-Sequencing (RNA-Seq) [10]. Each of these three methodologies has been widely used for quantifying mRNA expression, and all have undergone extensive validation, typically in “cross-platform” analyses (11–13). However, the inherent methodological differences in these platforms as well as limitations in starting material means that in sampling human tissue in situ platform agreement may be lower than in model systems [14,15]. Using a repeated measures longitudinal design helps counter these limitations [16], as does the use of a restricted candidate module with fewer genes and higher probability of change, but these may not be adequate to address inherent methodological challenges. Another design feature with the potential to improve reproducibility and reduce false discovery is the deployment of multiple platforms in contemporaneous measurements, allowing for cross-platform validation [17–19]. The primary aim of the current report was to investigate the utility of cross-platform validation of gene expression changes in the setting of a SAMGE-PC study, with a secondary aim of evaluating the effects of downsizing sample size into a greater phenotypic change cohort.

Materials and methods

Study material and datasets

Clinical study.

The data are derived from the Beta-blocker Effect On Remodeling and Gene Expression study (“BORG”, NCT01798992), the full design of which is described elsewhere [8]. Briefly, eligible idiopathic dilated cardiomyopathy (IDC) patients with HFrEF (LVEF ≤40%, symptomatic heart failure) were treated with one of three combinations of β-adrenergic receptor (AR) antagonists that had β1-AR blockade in common [8]. RV distal septal endomyocardial biopsies were performed prior to and after 3 and 12 months treatment with these regimens as previously described [8]. In the current analysis the 12 month results are used unless data were unavailable (n = 8), in which case 3 month data were used by last observation carried forward (LOCF) [8]. The primary method of measuring LV remodeling was by ejection fraction (EF, a measure of chamber structure and function [3]) by radionuclide SPECT imaging [3,8] for determination of LVEF or RVEF. All patients signed written consent for this multi-center study, which was conducted at the University of Colorado Anschutz Medical Campus and the University of Utah Medical Center and approved by the Institutional Review Boards at both sites. The study was conducted according to the Declaration of Helsinki and included a Data and Safety Monitoring Board.

RNA extraction and mRNA abundance measurements.

RNA extraction was performed on RV distal septum endomyocardial biopsy samples as previously described [8]. All mRNA measurements were from the same RNA extraction, typically involving 2–5 separate biopsy samples from an individual patient. Extracted RNA was stored at -80°C until use. Gene expression measurements were conducted by 3 platforms: TaqMan RT-qPCR (for 50 candidate genes) [8]; microarray (for global as well as candidate gene expression); and RNA-Seq (global and candidate gene expression). Candidate genes were selected for the likelihood they may mediate myocardial remodeling in response to neurohormonal activation or inhibition, within paradigms of either fetal/adult expression patterns, β1-adrenergic receptor signaling or thyroid hormone mediated effects [8]. RT-qPCR and microarray analyses were performed as previously described [8]. The microarray method was the Affymetrix HG-U133 Plus 2.0 Microarray that contains sequences for 22,963 separate gene sequences, with which 19,672 separate mRNAs were identified, including for all 50 of the candidate genes. RNA-Seq was conducted as previously described [10], using an Illumina HiSeq2500 platform and ribosome depleted RNA. RNA-Seq identified 15,366 separate gene transcripts, including the 50 candidate genes. RT-qPCR Ct values were normalized to GAPDH gene expression. Normalized changes in RT-qPCR gene expression from baseline to LOCF were calculated with the ddCt. method as previously described [8]. Microarray gene expression data were normalized by log-scale robust multi-microarray analysis. RNA-Seq data transcript levels were quantified as fragments per kilobase of exon per million mapped reads (FPKM) as previously described [10].

Patient cohorts and reverse remodeling phenotypic groups.

Two patient cohorts were analyzed: the BORG IDC cohort of 47 subjects, (“All-Subjects” (A-S)); and an A-S derived subcohort of the largest degree of reverse remodeling responders and control nonresponders (“Super-Responder” restricted (S-R) cohort). The 3 LV phenotypes as defined by their LVEF change measured at LOCF are therefore Responders (R) and Nonresponders (NR) in the A-S cohort, and Super-Responders (SR) in the S-R cohort. SRs have their own NR analysis group, which is the same phenotype as A-S NRs.

A-S Cohort. This IDC cohort consisted of 31 reverse remodeling Responders (R, n = 31) defined by an LOCF increase in LVEF of ≥5 absolute % at 3 months or 8% at 12 months, and 16 Nonresponders (NR) [8]. In R the baseline LVEF/RVEFs in %±SD were 25.6± 8.2/27.0± 8.7, and in NR 27.8± 9.9/28.2/9.5 (R vs. NR P = 0.42 in LVEF, 0.71 in RVEF). The LOCF—baseline absolute % increase in LVEF (mean±SEM) was 21.2±1.8 in R, vs. 0.9±1.1% in NR (P <0.0001). The absolute % change in RVEF was 9.7±2.2 in R, vs. 5.2±3.8 in NR (P = 0.30).

S-R Cohort. In this cohort “Super Responders” (SR, n = 6) were defined by an increase in LVEF of ≥10 absolute % (mean±SEM increase by 31±4%, P = 0.035 vs. A-S cohort R), and NRs (n = 6) by age and gender-matched patients with no LVEF increase ≥5% (mean change 0±2%, P <0.0001 vs. R). The mean±SD baseline and LOCF LVEFs in the S-R cohort were respectively 32± 9% and 32±8% in NR, and 24±7% (P = 0.63 vs. A-S R baseline) and 53±6% (P <0.001 vs. A-S R LOCF) in SR. RVEF mean change±SEM in SR was 16.4±8.1%, vs. 7.5±6.7% in NR (P = 0.42).

Summary of gene analysis groups, platforms and patient cohorts.

Fig 1 summarizes the platforms, patient cohorts, reverse remodeling phenotypes and types of gene expression analyses (candidate or global) investigated. A data table roadmap is given in the Supporting Information (S1 Table).

[Figure omitted. See PDF.]

Fig 1. Summary of LVEF changes from baseline.

LVEF change from baseline in the 47 patients comprising the All-Subjects (A-S) Cohort, and the 12 Subjects in the Super-Responder (S-R) subcohort. Messenger RNA measurement platforms, Responder/Nonresponder reverse remodeling phenotype and types of gene expression changes measured in each platform are also given.

https://doi.org/10.1371/journal.pone.0221519.g001

Data analyses

Data analysis within reverse remodeling groups.

Changes in gene expression measured as mRNA abundance were assessed within R, SR or NR reverse remodeling groups as change from baseline, using paired analyses of baseline vs. LOCF values of dCt for RT-qPCR, fluorescence intensity units for microarray and FPKM for RNA-Seq. Thus, the primary analysis in each phenotypic group was conducted with the untransformed readouts from each platform. The R/NR analysis was conducted by comparing the subtracted (LOCF–baseline) values in Rs to NRs. Fold changes for RT-qPCR and microarray data were determined by averaging the 2^subtracted values (2-ddCt for RT-qPCR) from each patient. For RNA-Seq FPKM, which is not an exponentiated value as for RT-qPCR and microarray data, the within R (or SR) and NR groups average fold change was determined by subtracting log2 transformed baseline values from the LOCF values for each patient, and then converting the difference to the base 2 antilog, i.e. 2^((Log2 LOCF FPKM)—(Log2 baseline FPKM)) and averaging for all members of the group. Fold change for R/NR was calculated by dividing the R by the NR fold change.

Statistical tests.

Baseline to LOCF changes in mRNA abundance within Responder (R) and Nonresponder (NR) groups were assessed by the Wilcoxon signed-rank test. Between group changes (R/NR analysis) were determined by Wilcoxon rank sum test. LVEF and RVEF were analyzed by t-tests, as the data were not significantly non-normally distributed by the Anderson-Darling test. Potential differences in categorical variables were assessed by Chi-Square or Fisher’s exact test for small cell sizes, with a Benjamini-Hochberg [20] adjustment for multiplicity. Alpha was set at 0.05 (2-sided P <0.05) for selection of differentially expressed genes (DEGs) on changes from baseline within R or NR groups, and in R/NR analyses. Rank correlation analysis was performed using Spearman’s Rho. Statistical analyses were performed on the R statistical package [21].

Cross-platform or cross-cohort concordance.

Concordant change was defined as the same directional change at a P <0.05 in each compared platform. Total Cross-Platform Concordance (CcpT) in % for mRNA abundance changes (up plus down-regulated genes, “DEGs”) between platforms within the same cohort (A-S or S-R) and response (NR, R or R/NR) was calculated from the formula [22]where intersect is the number of genes with concordant changes.

When applied to 2 platforms measuring an identical set of genes’ mRNA expression this formula is identical to that of Shi et al [23], used to derive “POG”, or percentage of overlapping genes, from 2 lists of platform specific gene expression changes. The same formula can be used to derive cross-cohort performance (CccT) in the same platform.

Gene categorization and network analyses.

Genes exhibiting concordant changes were categorized by a subgroup of the authors (JAS, EMR, MRB) based on manual curation of literature and NCBI website information relevant to myocardium. Canonical pathways, biological functions and gene networks were generated through Ingenuity Pathway Analysis (IPA) [24] (https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis). “Superheat” maps were generated by the Barter-Yu graphical tool using R [25]. Column and row ordering of heatmaps was determined by hierarchical cluster solution for LVEF change tertile 3 subjects.

Source data and protocol warehousing.

Source data files from the 3 measurement platforms plus clinical data are warehoused in the GitHub repository (https://bristowlab.github.io/borg.html). Endomyocardial biopsy, RNA extraction and mRNA measurement protocols are available at https://www.protocols.io/view/borg-protocols-2uggetw/abstract.

Results

Cross-platform concordance

Candidate genes in the A-S cohort.

Responder (R) group analysis. For the 50 candidate genes, in the A-S cohort we compared mRNA measurements determined by microarray and RT-qPCR (Table 1). In the change from baseline Responder analysis, 23 (46%) of the 50 candidate gene module exhibited changes by RT-qPCR, and 24 (48%) by microarray. For downregulated genes, in RT-qPCR there were 18 changes (P = 0.0009 vs. upregulated) and 14 in microarray. The CcpT for all concordant changes in both platforms was 60% (Table 1). In the R analysis 33 (66%) of the 50 candidate genes exhibited changes by either RT-qPCR or microarray.

[Figure omitted. See PDF.]

Table 1. Comparison of significantly different (P <0.05) candidate genes (of 50 total mRNAs measured) detected by RT-qPCR or microarray in the A-S cohort.

https://doi.org/10.1371/journal.pone.0221519.t001

R/NR analysis. For the R/NR analysis (Table 1) 18 (36%) of the 50 candidate genes exhibited changes by either RT-qPCR or microarray (vs. R, P = 0.059 for RT-qPCR, 0.023 for microarray, and 0.003 for both platforms). For upregulated genes RT-qPCR identified 11 and microarray 6 changes, all of which were concordant with RT-qPCR. For downregulated genes, RT-qPCR identified 2 (P = 0.007 vs. upregulated) changes that were both concordant with microarray measurements, which had 7 changes. The R/NR CcpT was 62%.

Candidate gene expression measured by three platforms in the S-R cohort.

In the S-R cohort all three platforms generated measurements, shown in Table 2. Only one change was concordant in all 3 platforms, marked downregulation (by >90%) in the natriuretic peptide NPPB, a biomarker of ventricular myocardial hypertrophy, in the R/NR analysis. The Na+/H+ antiporter SLC9A1 (NHE1), which is involved in mediation of β1-adrenergic signaling of hypertrophy [26], exhibited two platform concordance between RNA-Seq and RT-qPCR for R, and RNA-Seq and microarray for R/NR.

[Figure omitted. See PDF.]

Table 2. Candidate genes in S-R cohort, mRNA expression fold changes from baseline within NR, R groups or comparison of R, NR fold changes (R/NR), changes with P <0.05.

https://doi.org/10.1371/journal.pone.0221519.t002

The relatively small, N = 12 S-R cohort provided an opportunity to assess sensitivity of change detection among the 3 platforms. Of the 50 candidate genes, RT-qPCR, microarray and RNA-Seq yielded respectively 16, 9 and 17 statistically significant changes from baseline (R + R/NR analyses), P = 0.15 for change detection among platforms. The CcpT value within specific platforms calculated as agreement with at least 1 other platform trended lower in microarray vs. respective values in RT-qPCR or RNA-Seq (in R 0% vs. 31% and 44%, P = 0.02; in R/NR 38% vs. 100% and 62%, P = 0.18).

Global gene expression in the S-R cohort, measured by microarray and RNA-Seq.

Microarray and RNA-Seq platform measurements of global gene expression were used to investigate change concordance in the S-R cohort. Microarray detected 19,672 individual gene transcripts and 15,364 (78%) of these were also detected by RNA-Seq, forming the index number of measurements for analysis. S2 Table gives the alpha 0.05 list of up- or down-regulated genes within R and R/NR analyses, which are summarized in the denominators for change and concordance calculations at the bottom of Table 3.

[Figure omitted. See PDF.]

Table 3. Cross-platform concordance of global gene expression changes, mRNA measurements by array and RNA-Seq in the S-R cohort, responder and R/NR datasets, changes at alpha 0.05.

https://doi.org/10.1371/journal.pone.0221519.t003

R analysis. In the R analysis microarray detected 573 genes with mRNA changes, or 2.9% of the transcripts measured (279 upregulated and 294 downregulated, P = 0.53). RNA-Seq identified 1724 changes (11.2% of the transcripts measured, P <0.0001 vs. microarray), consisting of 631 upregulated and 1093 downregulated genes (P <0.0001). Concordance analysis between microarray and RNA-Seq (Table 3) yielded a total of 27 genes in the R analysis (15 upregulated and 12 downregulated), for a CcpT for any gene expression change of 2.4%.

R/NR analysis. In the R/NR analysis microarray identified 1805 genes with changed expression (9.2% of total mRNAs measured, 869 upregulated and 936 downregulated, P = 0.11), while RNA-Seq detected 1456 changes (9.5% of total measured transcripts), of which 378 were upregulated and 1078 downregulated (P <0.0001). The R/NR dataset exhibited 297 changes that were concordant between microarray and RNA-Seq (92 upregulated and 205 downregulated (P <0.0001)) (Table 3), with a combined platform CcpT of 17.5% (P <0.0001 vs. A-S R/NR (PCR/microarray) CcpT of 62%). The within S-R CcpT comparison between candidate (6/16, 38%) and global changes measured by microarray/RNA-Seq yielded a P value of 0.036. Table 4 contains a biologic categorization of genes with concordantly changed expression for the microarray/RNA-Seq R/NR analysis in the S-R cohort, with the category assignments for each gene given in S3 Table. For all 297 genes with cross-platform concordantly changed expression, the disproportionate downregulation by a factor of 2.2 fold compared to upregulation is due to changes in the Growth/Hypertrophy (P = 0.026, Hochberg adjusted P = 0.078), Cell Homeostasis (respective Ps 0.005, 0.027), Fibroblast Growth/extracellular matrix (Ps 0.001, 0.008), Extracellular to Intracellular Signaling (Ps 0.022, 0.077), Cytoskeleton (Ps 0.041, 0.10), and Unclassified/unknown Function (Ps = 0.002, 0.016) categories. S3 Table is a more detailed biologic categorization of genes with concordant expression changes in S-R, to which 2 additional candidate genes that were concordantly changed by RT-qPCR and RNA-Seq (upregulation in ATP2A2 and downregulation in CASQ2, Table 2) have been added.

[Figure omitted. See PDF.]

Table 4. Biologic categories of concordant global gene expression changes measured by microarray and RNA-Seq in the S-R cohort, R/NR analysis (from Table 3).

https://doi.org/10.1371/journal.pone.0221519.t004

RNA-Seq detected more genes with changed expression in R than in R/NR (1724 vs. 1456, p <0.0001), while microarray identified more changes in R/NR than R (1805 vs. 573, p <0.0001). CcpT was 6.2 fold higher in R/NR compared to R (Table 3). The CcpT of 2.4% in R for global gene expression was matched by a 0% (no concordance in 10 changes) for microarray/RNA-Seq in S-R (Table 2). For R/NR CcpT was higher, 17.5% for global (Table 3, P <0.0001 vs. R) and 38% (6/16) for candidate gene microarray/RNA-Seq (Table 2, P = 0.053 vs. R, 0.036 vs. global R/NR).

Cross-cohort comparisons

In order to derive direct information on whether smaller cohorts of hyper-responders can substitute for larger cohorts comprised of smaller degrees of phenotypic response in longitudinal reverse remodeling gene expression studies, we compared measurements between the A-S and S-R cohorts using data from RT-qPCR for candidate genes and microarray for global gene expression.

Cross-cohort concordance for candidate gene expression changes in A-S vs. S-R assessed by RT-qPCR.

R analysis. Genes with p <0.05 changes in expression for RT-qPCR measurements, the protocol primary method for candidate gene mRNA change detection [8], are given for the 2 cohorts in Table 5. In the A-S cohort R analysis, 23 genes exhibited expression changes (5 upregulated and 18 downregulated, P = 0.002), while 13 changes were identified in the S-R cohort (3 upregulated, 10 downregulated, P = 0.037) (P = 0.059 vs. total changes in A-S). Thus, in the R analysis RT-qPCR identified more down-regulated genes in both the A-S and S-R cohorts. The two cohorts were concordant for 2 upregulated (ADRB1 and PLN) and 8 downregulated (GNAI2, SLC9A1, AGT, HK2, PDHX, CASQ2, SLC8A1, DMD, NPPB and IL6) genes (Table 5, P = 0.046). Within the R group the 10 concordantly changed genes in the A-S and S-R cohorts were from a total of 35 upregulated and downregulated genes, for a CccT of 56% (Table 5).

[Figure omitted. See PDF.]

Table 5. Cross-cohort candidate gene concordance of RT-qPCR in the A-S vs. S-R cohorts (changes of P < 0.05).

https://doi.org/10.1371/journal.pone.0221519.t005

R/NR analysis. Within the R/NR dataset of the A-S cohort 11 genes were upregulated while 2 genes were downregulated (Table 5, P = 0.007) (Table 5). In the S-R cohort 1 gene was upregulated while 2 were downregulated (total changes of 3 vs. 13 in A-S, P = 0.006). The 2 cohorts were concordant for 1 upregulated (ATP2A2) and 1 downregulated gene (NPPB), for a CccT of 25% (Table 5; P = 0.033 vs. R Group analysis).

Cross-platform comparisons between S-R and A-S cohorts. For candidate genes cross-platform concordance was compared between the S-R and A-S cohorts, using the 2 platform PCR/microarray comparison. Compared to the A-S cohort (Table 1), for the R analysis PCR/microarray CcpT in S-R (Table 2) was 0/14 = 0% (P = 0.01 vs. CcpT of 60% in A-S), and 2/11 = 18% for R/NR (P = 0.016 vs. CcpT of 62% in A-S). Compared to PCR/microarray the S-R PCR/RNA-Seq CcpT values were higher, 8/22 (36%) (P = 0.03 vs. S-R PCR/microarray) in R and 6/11 (55%, P = 0.08) in R/NR (Table 2). The S-R R and R/NR CcpT values for microarray/RNA-Seq were respectively 0/10 (P = 0.07 vs. PCR/RNA-Seq) and 6/16 = 38% (P = 0.38 vs. PCR/RNA-Seq) (Table 2). Thus, candidate gene cross-platform concordance for PCR/microarray was not as good in S-R as compared to the A-S cohort, but the higher 2 platform concordances of PCR/RNA-Seq vs. microarray containing CcpT combinations suggests part of the issue was microarray platform related.

Cross-cohort concordance for global gene expression changes as assessed by microarray.

R Analysis. The number of significantly different global genes based on microarray measurements in the S-R cohort is given in S2 Table, and for the A-S cohort in S4 Table. The R group analysis yielded a total of 2931 upregulated genes in the A-S cohort (represented in the denominator of total changes in S5 Table, while 2294 upregulated genes were detected in the S-R cohort (p <0.0001). A total of 3274 downregulated genes were identified in the A-S cohort while 2095 genes were detected in the S-R cohort (p <0.0001). There were only 60 upregulated and 42 downregulated concordant changes, for a CccT of 1.9% (S5 Table).

R/NR analysis. In the R/NR patient group microarray detected 850 upregulated genes in the A-S cohort (S4 Table), while 869 were detected in the S-R cohort (S2 Table). Downregulated genes were n = 1085 in the A-S cohort (S4 Table), and 936 in the S-R cohort (S2 Table, P = 0.0006 vs. A-S). In marked contrast to the R analysis, the R/NR analysis had much higher concordance in both up- (N = 391) and down- (N = 393) regulated genes, with a CccT of 42% (P <0.0001 vs. R analysis) (S5 Table).

Platform concordance of gene expression changes has biologic significance

The A-S cohort, which had LVEF changes ranging from -8 to 47.5 absolute %, created the opportunity to examine the relationship between reverse remodeling as a biologic effect and cross-platform concordance of gene expression changes. For changes from baseline, candidate gene microarray and RT-qPCR CcpT values were related to LVEF change tertiles (Fig 2, S6 Table). The first tile of LVEF change (Fig 2) corresponds to the NR group (LVEF change in absolute %±SD = 0.9±4.4%), while the 2nd and 3rd tiles consist of R group subjects that have respective LVEF increases of 13.3±4.0 and 29.6±6.1 absolute %. There were no differences in baseline LVEF values, which ranged from 27.8% in tile 1 to 23.6% in tile 3 (ANOVA P = 0.32). The tiles 2 and 3 LOCF LVEFs were respectively 40.9% and 53.1% (Fig 2), indicating phenotypic improvements into the HFmrEF (HF with “mid-range” LVEF) and normal ranges. The tile 1 LVEF remained in the lower range of HFrEF (P <0.0001 across the 3 LOCF tiles).

[Figure omitted. See PDF.]

Fig 2. LVEF change and CcpT by tertile.

(A) LVEF change by tertile, mean ± SD in absolute %; (B) Gene expression change from baseline at P <0.05, Overall Cross-platform Total % Concordance (CcpT) by tile. LVEF baseline → LOCF* values (mean in % ± SD): Tile 1, 27.8±10.3 → 28.7±10.3; Tile 2, 27.6±8.4 → 40.9±10.5; Tile 3, 23.6±7.7 → 53.1±6.7. ANOVA P values: baseline = 0.32, LOCF <0.0001.

https://doi.org/10.1371/journal.pone.0221519.g002

By RT-qPCR there were respectively 25, 12 and 20 up or downregulation changes progressing up the tertile (P = 0.026) (Fig 2). Microarray measurements yielded 2, 13 and 17 respective changes (P = 0.008) by tertile. For both platforms the number of changed genes by tertile was 27, 25 and 37 (P = 0.033) progressing from Tiles 1 to 3, but tile 1 was not significantly different from tile 3 (P = 0.083) and there is no obvious reverse remodeling vs. gene expression change relationship, with tile 2 having numerically fewer changes than tiles 1 or 3. However, RT-qPCR and microarray cross-platform concordance exhibited a direct relationship to reverse remodeling across the LVEF tiles, with respective CcpT values of 15%, 32% and 59% (overall P < 0.001) (Fig 2, S6 Table). Thus, in a set of candidate genes enriched for having a possible relationship to LV remodeling, concordance of gene expression changes increases directly and substantially with the degree of reverse remodeling.

Gene network/signaling pathway analyses

Pathway analyses of concordant R/NR changes.

Canonical pathways. To investigate the possible canonical signaling pathways involved in the concordance related biological effect on LV remodeling, for microarray/RNA-seq concordant changes in global gene expression in S-R we conducted an Ingenuity Pathway Analysis (IPA) [24] on the 299 R/NR concordant genes listed in S3 Table with the results shown in S7 Table. The methodological background for the IPA analyses is given in the Supporting Information. S3 Table indicates there were 92 upregulated and 207 downregulated genes based on 297 concordantly changed genes in Table 3 and an additional 2 in Table 2. Of the most negative z scores in S7 Table (values < -1.5), evidence for a net pathway inhibition, 3 of them were phosphoinositide related pathways containing many of the same molecules (Superpathway of Inositol Phosphate compounds, 3-phosphonositide Degradation and D-myo-inositol-5-phosphate Metabolism), belonging to the biologic categories (Table 4, S3 Table) of Growth/hypertrophy regulation (PTPRH, CDC25A), extracellular to intracellular signaling (PLCE1, PLCD4 (subcategory phospholipase C signaling), PPP1R1A (subcategory neurohormonal signaling); PIK3R2, INPP4B, MTMR3 (subcategory PI3 kinase signaling); and SIRPA (subcategory other signaling)). The other 2 pathways identified with z scores < -1.5, GP6 Signaling and Dendritic Cell Maturation contained collagen and phospholipase C or PI3 kinase signaling molecules. GP6 is a collagen receptor expressed in platelets, but the IPA is likely detecting extracellular matrix collagen transcript downregulation.

Only 6 pathways had positive z scores (indicating net activation), with the highest, calcium signaling, considered statistically significant with a z-score of 2.00 and an overlapping P value of < 0.01. Cardiac β-adrenergic signaling and protein kinase A signaling also had positive z-scores, in a clinical setting of high dose β1 adrenergic receptor blockade. Downstream of β1-adrenergic signal transduction is calcium signaling and the initiation and regulation of contractile function [27]. The 8 mRNA gene products grouped together within the calcium signaling Ingenuity canonical pathway in S7 Table are involved in calcium binding/uptake or release (GRIN2A, ATP2A2, CASQ1, CASQ2), or are directly or indirectly regulated or activated by calcium (RCAN1, MYL3, TPM3, TNNT1). In S3 Table they are members of the channels/solute exchanger/transporters (GRIN2A), Ca2+ handling or signaling (ATP2A2, CASQ1, CASQ2), contractile and associated proteins (MYL3, TPM3, TNNT1) or growth/hypertrophy regulation (RCAN1).

The second most positive z-score in the IPA was for Gα12/13 signaling (S7 Table), which has not been previously identified as having a molecular signature in the reverse remodeled failing human heart. GNA12 and GNA13 are closely related members of a class of pertussis-toxin insensitive G-proteins [28] that transduce signaling in lysophosphatidic acid (LPA), angiotensin II, endothelin-1, MAP kinase, PI3 kinase and other pathways involved in cardiac myocyte hypertrophy [28–31]. The Gα12/13 signaling module contained an upregulated LPA receptor (LPAR3), an upregulated MAP kinase kinase (RAF1), 2 downregulated cadherin genes (CDH2 and CDH17), a downregulated PI3 kinase regulatory subunit (PIK3R2) and an upregulated contractile protein gene, MYL3, the “essential” ventricular myosin light chain that is considered to be an “adult” expressed gene within the fetal-adult gene expression paradigm.

Noncanonical disease associations.

Fig 3 is the “tox function” IPA analysis for 4 domains relevant to ventricular chamber remodeling. Failure (A), Dysfunction (B), Cell Death (C), and Fibrosis (D) of the heart are all predicted to be inhibited by the noncanonical gene expression changes in the S-R cohort (Tables 2 and 3).

[Figure omitted. See PDF.]

Fig 3. Ingenuity pathway analysis of concordantly changed mRNAs from Table 3.

Colors indicate predicted effect of transcript on biological function or disease in center of circular diagram. (A) Dysfunction of heart. (B) Failure of heart. (C) Cell death of heart. (D) Dilated cardiomyopathy.

https://doi.org/10.1371/journal.pone.0221519.g003

Correlation of candidate with global gene expression changes in A-S cohort

In the A-S cohort we used the “Superheat” [25] graphical package to assess the extent of microarray-measured candidate gene expression changes compared to those in the 19,672 global genes (Fig 4), within the LVEF change tertiles shown in Fig 2. In Fig 4 the 19 genes that had concordant R or R/NR changes (Table 1) are designated by arrows, and the 19,672 gene transcripts are displayed in individual horizontal rows that cannot be distinguished in the resolution. The heat map is arranged by degrees of correlated fold change ranging from a Rho of +1.0 (two genes, a candidate and a global) with log2 fold changes perfectly positively correlated to -1.0 (2 genes perfectly negatively correlated) in the 47 individual patients. The candidate genes along the x-axis are arranged by degree of correlation with the 19,672 global genes identified on microarray in tile 3 (Fig 4C); closer spatial proximity of candidate gene columns suggests similar patterns of correlation across non-candidate genes. In the 3rd tile of LVEF change (Fig 4C, mean LVEF increase 29±6 absolute %) clear patterns of moderately positive or negative correlations are evident, which are not present in tiles 1 and 2 (Fig 4A and 4B; respective mean LVEF changes of 0.9±4 and 13±4 absolute %). In tile 3 (Fig 4C) there are 2 sets of candidate genes associated with reciprocal regulation of large numbers of global genes, which also contain a high percentage of concordant genes. Along the horizontal axis are candidate genes, ranging from GNAS to MYH7, which are positively (bottom part of the heat map) or negatively (middle-upper part) correlated of with large numbers of global genes’ expression. Of the 19 candidate genes in this region, 7 are concordantly changed, with 6 being upregulated including both ADRB1 and ADRB2. The positive and negative portions of the map encompass approximately 28% and 30% of the 19,672 genes displayed. The 19 candidate genes in this “hot” region (“Group 1”) have in common either being regulated by β1-AR signaling (n = 12; AGTR1, HNRNPD, MYL3, RYR2, ADRB2, PLN, ADRB2, PFKM, CTF1, ATP2A2, ADRA1A, MYH7), being involved in adrenergic receptor signal transduction (GNAS, GNAQ), being involved in Ca2+ handling or regulation (ATP2A2, SLC8A1, PLN, CANX), or being important in sarcomere structure (DMD, CTF1).

[Figure omitted. See PDF.]

Fig 4. Superheat maps.

(A) Spearman’s Rho of baseline vs. LOCF mRNA abundance, superheat map of 19,672 separate mRNA gene products (rows along Y-axis) vs. the 50 candidate gene mRNAs, patients in tile 1 of LVEF change (Fig 2). Arrows indicate candidate genes with concordant changes measured by RT-qPCR and microarray in either the R (Responder) or R/NR (Responder/Nonresponder) analysis; (B) Same as Fig 4A., except data are from tile 2 in Fig 2; (C) Same as Fig 4A., except data are from tile 3 in Fig 2.

https://doi.org/10.1371/journal.pone.0221519.g004

Beginning with TNF and continuing through SLC9A1 on the Fig 4C heat map are 17 candidate genes (“Group 2”) which, compared to the previous 19 candidate genes, reciprocally correlate with global genes (negatively correlated in the bottom part of the map and positively in the mid-upper region). Eight of these 17 were concordant, with 6 downregulated. This gene cluster contains 3 cytokines (TNF, IL1B and IL6), the thyroid hormone receptor THRA1, and several β1-adrenergically regulated genes (ADRBK1, MYH6, GNAI2, MYL4, HK2, CPT-1 and SLC9A1). Finally, the third region on the map beginning with PDK4 and proceeding through a total of 11 genes does not contain any obvious clustering of correlated mRNAs, similar to a 3 gene (ADCY5, CASQ2, MYL2) transition area between the two other regions. Of this total of 14 genes without correlative changes in expression, 4 (and CSRP3) were concordant, all downregulated. NPPA, NPPB, and ACTA1 are regulated by β1 adrenergic signaling. The other 7 members of this group consist of 3 metabolic genes (the β1-adrenergically regulated PDK4, PDHX and the reference gene GAPDH, 2 contractile proteins (TNNT2, TNNC1), and 2 neurohormonal genes (EDN1, AGT)).

Discussion

We used cross-platform concordance of gene expression changes as a method of increasing the reliability of detecting biologically meaningful molecular changes associated with reverse remodeling in an idiopathic dilated cardiomyopathy study population. Cross-platform concordance has been widely used as a method to validate measurement methodologies [17–19,22,23], but in model systems it is also correlated with biologic and/or therapeutic signals [22,32–34]. The use of cross-platform concordance was validated through a highly statistically significant relationship of concordant candidate gene changes with progressive reverse remodeling. We then extended the use of concordant gene expression changes to pathway analyses, which revealed both expected and unexpected results.

Cross-platform concordance of changes in gene expression

CcpT values, a standard measurement of expression change concordance, were in the 60% range for RT-qPCR and microarray measurements of candidate genes in the 47 patient A-S cohort, for both R and R/NR analyses. For the smaller (n = 12) S-R cohort that exhibited a greater degree of reverse remodeling, CcpT values were somewhat lower, particularly in the R analysis on 2 platform comparisons that included microarray data. The lower CcpT values in the S-R cohort likely reflect the difference in statistical power between the larger and smaller cohorts, plus possibly somewhat lesser accuracy of transcript detection by microarray compared to RT-qPCR or RNA-Seq.

In contrast to candidate genes in the A-S cohort, concordance of changes for global genes as assessed by microarray and RNA-Seq in the S-R cohort were extremely low in the change from baseline R analysis, with a CcpT of 2.4%. However, the microarray/RNA-Seq CcpT for candidate genes in the S-R cohort was 0%, suggesting that the smaller sample size S-R cohort was ill-suited for the R analysis. The R/NR analysis CcpT for global gene expression was much higher, 17.5%, but lower than the counterpart CcpT microarray/RNA-Seq value for candidate genes (38%, P <0.036). These data emphasize the value of comparing R vs. NR changes (R/NR analysis), with the implication that an in situ myocardial gene expression study needs to be large enough to generate both reverse remodeling responders and nonresponders. In addition, the data support an expectedly higher concordance for candidate genes, which are enriched for genes likely to change expression in reverse remodeling. The larger global gene CcpT for R/NR vs. R also supports the use of the R/NR analysis in pathway generation.

Cross-cohort concordance of changes in gene expression

Cross-cohort concordance, CccT, was measured in data generated by a single platform in a manner analogous to inter-platform concordance. For candidate gene mRNAs measured by RT-qPCR the A-S cohort tended to have (p = 0.059 in the R analysis) or had (p <0.01 in the R/NR analysis) more concordant changes, again indicating that the larger, A-S cohort had a statistical advantage that was not compensated by the greater degree of reverse remodeling in S-R. However, the A-S/S-R CccT of 56% in R and 25% in R/NR (Table 5) supports comparability of the S-R and A-S cohorts, not surprising since S-R is a subset of A-S.

In A-S and S-R changes in expression of global genes was assessed by microarray. In the R analysis the A-S cohort yielded proportionally more statistically significant changes than in S-R (by 41%), whereas in the R/NR analysis the number of changes was nearly equal in the 2 cohorts, with fewer than half the number of changes in both A-S and S-R. The CccT for the R analysis was extremely low, 1.9%. In contrast, the CccT for R/NR was markedly higher, 42%, and within the range for R or R/NR CccT for candidate genes. The basis for the marked difference in CccT values was an R/NR vs. R increase in concordant genes by 7.7 fold and a decrease in the number of either-cohort changed genes by 2.8 fold, with corresponding implied reductions in both Type 2 and Type 1 error in the R/NR analysis. These data provide support for the validity of the S-R cohort as representative of global gene expression changes in A-S, but only for an R/NR analysis.

Biologic significance of platform concordance for gene expression changes

In the A-S cohort an analysis of platform concordance of changes in gene expression by degree of LVEF change indicated that most of the changes from baseline in the 1st LVEF tertile, comprised entirely of Nonresponders, were non-concordant, and therefore likely to be false positives. In contrast, in patients undergoing reverse remodeling (tertiles 2 and 3) concordance was directly related to the degree of LVEF improvement and reached 59% in tile 3, where the increase in LVEF was by 29.6 absolute %. Candidate gene change analysis in the S-R cohort also provided support for the biologic significance of concordance. Two platform CcpT exhibited a strong trend for an increase from NR to R to R/NR, reaching a value in R/NR (58%) that was quite similar to the A-S R/NR value of 62% despite the lower power of the S-R cohort. A straightforward interpretation of these data is that when comparing organ phenotypes, platform concordance for gene expression measurements is directly related to the degree of phenotypic change; that is, concordant changes have biologic meaning.

Gene network/signaling pathway analyses for concordant S-R gene expression changes in microarray and RNA-Seq

The demonstration of an acceptable degree of platform and cohort concordance for the R/NR analysis in S-R created the opportunity to perform a network analysis, for purposes of hypothesis generation. In this dataset downregulation exceeded upregulation, and the biologic categories of gene changes were unsurprisingly led by the growth/hypertrophy, cell homeostasis, fibroblast/extracellular matrix and extracellular to intracellular signaling categories. These changes are all consistent with reverse remodeling with inactivation of genes involved in myocyte and interstitial cell growth.

Despite all study subjects receiving high doses of β-blocking agents, cardiac β-adrenergic signaling and protein kinase A signaling had positive z-scores. All study subjects were receiving target/guideline recommended doses of metoprolol succinate or carvedilol that produce high levels of blockade of β1-adrenergic receptors, and both agents have been shown to reverse adrenergically-mediated desensitization changes in failing hearts [35,36]. Although the low mRNA abundance ADRB1 gene was not concordantly upregulated in the R/NR global gene expression analysis in S-R (Table 3), it was by RT-qPCR in the R/NR analysis in A-S (Table 1), by PCR or RNA-Seq R analysis in S-R (Table 2), and by RNA-Seq in the R/NR analysis in S-R (Table 2). ADRB2 was concordantly upregulated in A-S R/NR (Table 1). Restoration of β-receptor signaling could have produced a net activation of downstream cAMP/PKA and calcium signaling via receptor constitutive activity or by adrenergic signaling intermittently overriding β-blockade. The cAMP/PKA canonical pathway contained 3 upregulated and 5 downregulated genes whose protein products have been shown to regulate cAMP levels or be modifiable by PKA phosphorylation. The canonical calcium signaling pathway featured upregulation of ATP2A2 (SR Ca2+ ATPase) downregulation of CASQ1 and CASQ2, the respective minor and major abundance calsequestrins expressed in the heart, and the calcium regulated contractile proteins MYL3, TNNT1 and TPM3. However, expression changes of other modulators of Ca2+ signaling in the heart were not included in the Ingenuity cluster, such as upregulation of the phosphatase 1 inhibitor PPP1R1A that would be expected to promote SR Ca2+ uptake and release by inhibiting phospholamban dephosphorylation. The activation or improvement in calcium signaling as the LV reverse remodels is supported by work in model systems [37].

There were also some interesting findings within the IPA canonical pathway analysis involving noncandidate genes that have not been previously reported in ventricular reverse remodeling in the human heart. Phosphoinositide signaling occupied 3 of the top 4 negative z score or inhibited pathways, 2 pathways containing numerous collagen transcripts were in the first 5 inhibited pathways, and Gα12/13 signaling was the second highest ranked of only 6 activated pathways. None of these pathways had z-scores with absolute values of >2.0, so the findings are only hypothesis generating. Neither phosphoinositide or Gα12/13 signaling have been prominently identified in our [38] or other large scale [39] human myocardial remodeling gene expression studies, but the role of phosphoinositide signaling in pathologic hypertrophy and heart failure have been reported in numerous animal model studies [40–45]. Importantly, the PI3 kinase signaling pathway interacts with β1-adrenergic receptor expression and signaling [46], providing biologic plausibility for phosphoinositide pathway gene expression changes in that in adult rat cardiac myocytes β1-adrenergic receptor signaling leads to increased PI3K activity which then attenuates adrenergically mediated positive inotropic effects [46]. Similarly, β1-adrenergic receptor signaling has been shown to be a powerful mediator of fibrosis and other forms of extracellular matrix deposition in model systems [40,47–50], but generally has not been considered a major therapeutic target of β-blockers in reverse remodeling. In the human heart [51] Gα12/13 transduces hypertrophic signals [52] for a number of non-β-adrenergic receptors such as endothelin-ETA and angiotensin II AT1, and Gα12 mRNA expression is upregulated in the transition of pathologic hypertrophy to failure in the Dahl salt-sensitive rat model [52]. Despite this evidence that upregulation of Gα12/13 signal transduction can produce pathologic hypertrophy, several lines of evidence indicate that this signaling module may be cardioprotective [53,54], and the contractile protein member (MYL3) is an "adult" gene within the fetal/adult expression paradigm [55] that regulates contractility [56].

Of note is that the IPA canonical pathway Cardiac Hypertrophy Signaling had a neutral activation z-score, with an overlapping P value of–log 1.82 (0.015) (S7 Table), despite the biologic category of Growth/hypertrophy regulation having 21 downregulated vs. 8 upregulated genes (S3 Table). This discrepancy is likely due to differences in how the two classification instruments are constructed; the biologic category (S3 Table) placed genes with demonstrated changes into groupings based on curation of the current myocardium relevant literature and NCBI websites, while the IPA inserted the same changes into a proprietary algorithm where genes are assigned to signaling modules.

Correlation of candidate with global gene changes in the A-S cohort

In a correlation analysis between changes in candidate and global gene expression in the A-S cohort two potential patterns were identified that may be candidates for extensive gene regulation networks. One group (Group 1) had in common β1-AR gene signaling, for which there is considerable previous evidence [8,57]. Group 2 included GNAI2, which regulates multiple signaling pathways involved in hypertrophy including LPA signaling [29,30]. Both these potential networks, which exhibited reciprocal regulation, appeared to contain an extensive number of genes likely numbering in at least the hundreds based on the size of the area on the heat map and assuming that only a small percentage (≤10%) of genes is imparting the color to the map. A third group of genes with no strong clustered correlations was also identified, with the natriuretic peptides NPPA/NPPB as the centerpiece.

Basis for the advantage of cross-platform concordance in human in situ studies of gene expression

We evaluated gene expression change concordance as a means to improve the likelihood of biologic relevance of differential gene expression associated with reverse remodeling of failing, eccentrically hypertrophied human LVs. In terms of failure to achieve change concordance, both technical/methodological and biological causes are possible. For cross-platform concordance in the A-S cohort the majority of the non-concordance was presumably technical/methodological related, most likely related platform sensitivity and specificity of transcript identification and quantitation. However, endomyocardial biopsy sampling and RNA extraction variability would also fall within this category. For biologic sources of non-concordance, within the A-S R and NR groups (i.e. the components of the R/NR analysis) variance in degrees of remodeling and the relationship between remodeling and change in individual patient gene expression are also potential sources. For the much smaller sample size S-R cohort, both types of causes would be accentuated, and for the cross-cohort comparisons biologic causes would be a relatively greater contributor to nonconcordance since the responder phenotypes were quantitatively different.

All 3 methodologies used for cross-platform measurements have been extensively validated against each other [11–13,22,23], and with the possible exception of microarray in the R analysis in the S-R cohort, they performed similarly in this study. Compared to RT-qPCR and RNA-Seq microarray may have been at some disadvantage for low abundance mRNAs, but nevertheless was able to detect changes in some of them such as β-adrenergic receptors. TaqMan or Real Time qPCR has long been considered the “gold standard” for RNA quantitation but is impractical for measurement of global gene expression. Microarrays revolutionized the measurement of genome-wide transcriptomics, and held sway as the dominant methodology for over a decade [58]. Due to its flexibility for measuring multiple RNA gene products and now being cost-effective, RNA-Seq has replaced microarrays as the most widely used methodology for identifying and quantifying multiple RNA species [59].

As the methods for quantifying mRNA abundance are different for the 3 technologies, the statistical probability of achieving a P <0.05 with the same directional change across 2 platforms is 0.052/2 or 0.00125. Therefore, some of the value of cross-platform concordance can be attributed to statistical stringency. However, the utility of a 2nd, confirming measurement in real time by different methodology likely confers value for concordance beyond setting a lower P value to minimize false discovery rate (FDR). In gene expression experiments that include large numbers of measurements problems with controlling FDR without having an adverse effect on sensitivity are well known [60,61], and contemporaneous measures of cross-platform concordance can be viewed as an alternative to using stringent P values or FDR adjusted q values to lower Type I error.

Limitations

The data presented and the conclusions reached apply only to a SAMGE-PC design and may not be applicable to cross-sectional studies, as the general advantages of longitudinal, repeated measures designs [62] certainly apply to gene expression studies [16]. The study used endomyocardial biopsies taken from the distal RV septum as starting material for RNA extraction. Although the septum is an interventricular shared wall whose gene expression tracks with LV [6,8,63] and RV [63] remodeling, tissue from the right side of the septum may not be as representative of LV molecular remodeling as LV biopsies. However, in this study only the LV demonstrated statistically significant reverse remodeling, and therefore the observed changes in gene expression are likely to be representative of LV biologic effects. For cost reasons we confined RNA-Seq to a 12 patient Super-Responder subcohort, and ideally all 47 patients would have had these measurements. The S-R subcohort was a subset of the larger A-S cohort, and therefore was not independent of the parent population. Finally, since RNA-Seq has become the favored method of measuring transcript abundance investigators may be reluctant to add another platform to their study design, and it is not clear if performing two separate RNA preparation and sequencing runs could substitute for mRNA abundance determination in a completely separate platform. Nevertheless, a cross-platform concordance design performed contemporaneously as opposed to for post hoc confirmation should be considered for in situ human gene expression studies, where starting material is limited and repeating experiments is untenable.

Conclusions

For measurement of myocardial gene expression in sequential studies that include a baseline control value, among RT-qPCR, microarray and RNA-Seq methods there was no convincing advantage of one platform over another. However, for measurements of both candidate and global gene expression the use of 2 platforms had an advantage over a single platform, allowing for the use of cross-platform concordance as a more robust measure of biologically-related changes in gene expression. The basis for the 2-platform concordance advantage is likely related to the derived alpha calculation of 0.00125, but also to the requirement that independent measurement technologies yield conventional significance levels. Measurement concordance is more valuable for global gene expression, where the false discovery rate is higher than for candidate genes preselected for potential involvement in the remodeling process, but concordance also adds value in a candidate/hypothesis driven setting. Multiple lines of evidence favored the use of the larger but less reverse remodeled A-S cohort over the smaller, more improved phenotype S-R subcohort, and global gene expression measurements were only available in S-R. Nevertheless, a 42% CccT between A-S and S-R for global R/NR changes indicated that for this analysis the S-R subcohort was representative of the A-S cohort. Finally, both IPA and biologic category analyses support an improvement in calcium signaling as the major gene expression correlate of reverse remodeling, consistent with the improvement in systolic function that is associated with improved LVEF as well as previous work in model systems.

IPA and Superheat map investigation of concordant gene expression changes in the R/NR global gene expression analysis in S-R identified the likely involvement of unexpected canonical pathways in the reverse remodeling process, as well as potentially large networks involving β-AR and possibly GNAI2 regulation.

The data support the advantage of sequential measurements of gene expression before and after organ phenotype modification, particularly if phenotype responders can be compared to nonresponders in an R/NR design.

Measuring human ventricular myocardial gene expression in situ by serial endomyocardial biopsy remains a research tool only, but has the potential, particularly by employing multi-platform measurements in larger cohorts undergoing reverse remodeling, to identify novel mechanisms important in myocardial failure and hypertrophy.

Supporting information

[Figure omitted. See PDF.]

S1 Text. List of abbreviations and acronyms.

https://doi.org/10.1371/journal.pone.0221519.s001

(DOCX)

S2 Text. Ingenuity Pathway Analysis (IPA) Background.

https://doi.org/10.1371/journal.pone.0221519.s002

(DOCX)

S1 Table. List of tables comparing platforms in A-S and S-R cohorts, or cohort comparisons.

https://doi.org/10.1371/journal.pone.0221519.s003

(DOCX)

S2 Table. Up or downregulated genes within the R and R/NR analyses, microarray or RNA-Seq measurements in the S-R cohort.

https://doi.org/10.1371/journal.pone.0221519.s004

(DOCX)

S3 Table. Biologic categories of 299 concordant gene expression changes identified in the R/NR analysis by microarray and RNA-Seq in the S-R cohort (Table 3), and RT-qPCR and RNA-Seq in the S-R cohort (S1 Table).

https://doi.org/10.1371/journal.pone.0221519.s005

(DOCX)

S4 Table. Up or downregulated genes within the R and R/NR analyses, microarray measurements in the A-S cohort.

https://doi.org/10.1371/journal.pone.0221519.s006

(DOCX)

S5 Table. Cross-cohort concordant global gene expression changes, mRNA measurements by microarray in the A-S and S-R cohorts, Responder and R/NR datasets.

(P <0.05 in each cohort).

https://doi.org/10.1371/journal.pone.0221519.s007

(DOCX)

S6 Table. Gene expression changes from baseline within LVEF change tertiles.

P <0.05 by Wilcoxon signed rank for change from baseline.

https://doi.org/10.1371/journal.pone.0221519.s008

(DOCX)

S7 Table. IPA canonical pathways with overlap P values < 0.05, microarray, RNA-Seq and RT-qPCR mRNA measurements in the R/NR analysis.

N = 299 concordant gene mRNA changes in the S-R cohort biologically classified in S3 Table.

https://doi.org/10.1371/journal.pone.0221519.s009

(DOCX)

Acknowledgments

The authors thank Rachel Rosenberg and Laura Hofstatter for their assistance in manuscript preparation, proofreading and document uploading.

Citation: Toni LS, Carroll IA, Jones KL, Schwisow JA, Minobe WA, Rodriguez EM, et al. (2019) Sequential analysis of myocardial gene expression with phenotypic change: Use of cross-platform concordance to strengthen biologic relevance. PLoS ONE 14(8): e0221519. https://doi.org/10.1371/journal.pone.0221519

References

1. Hunt SA, Abraham WT, Chin MH, Feldman AM, Francis GS, Ganiats TG et al. 2009 Focused Update Incorporated Into the ACC/AHA 2005 Guidelines for the Diagnosis and Management of Heart Failure in Adults: A Report of the American College of Cardiology Foundation/American Heart Association Task Force on Practice Guidelines: Developed in collaboration with the International Society for Heart and Lung Transplantation. Circulation. 2009;119: e391–479. pmid:19324966

2. Yancy CW, Jessup M, Bozkurt B,Butler J, Casey DR Jr, Colvin MM, et al. 2017 ACC/AHA/HFSA Focused Update of the 2013 ACCF/AHA Guideline for the Management of Heart Failure: A Report of the American College of Cardiology/American Heart Association Task Force on Clinical Practice Guidelines and the Heart Failure Society of America. J Am Coll Cardiol. 2017;70:776–803. pmid:28461007

3. Bristow MR, Kao DP, Breathett KK, Altman NL, Gorscan J, Gill EA, et al. Structural and Functional Phenotyping of the Failing Heart: Is the Left Ventricular Ejection Fraction Obsolete? JACC Heart Fail. 2017;5:772–781. pmid:29096787

4. Benjamin EJ, Virani SS, Callaway CW, Chamberlain AM, Chang AR, Cheng S, et al. Heart Disease and Stroke Statistics—2018 Update: A Report From the American Heart Association. Circulation. 2018;137:e67–e492. pmid:29386200

5. Gerber Y, Weston SA, Redfield MM, Chamberlain AM, Manemann SM, Jiang R, et al. A Contemporary Appraisal of the Heart Failure Epidemic in Olmsted County, Minnesota, 2000 to 2010. JAMA Intern Med. 2015;175:996–1004. pmid:25895156

6. Lowes BD, Gilbert EM, Abraham WT, Minobe WA, Larrabee P, Ferguson D, et al. Myocardial Gene Expression in Dilated Cardiomyopathy Treated with Beta-Blocking Agents. N Engl J Med. 2002;346:1357–1365. pmid:11986409

7. Lowes BD, Zolty R, Minobe WA, Robertson AD, Leach S, Hunter L, Bristow MR. Serial Gene Expression Profiling in the Intact Human Heart. J Heart Lung Transplant. 2006;25:579–588. pmid:16678038

8. Kao DP, Lowes BD, Gilbert EM, Minobe WA, Meyer LK, Ferguson DA et al. Therapeutic molecular phenotype of β-blocker-associated reverse-remodeling in nonischemic dilated cardiomyopathy. Circ Cardiovasc Genet. 2015;8:270–283. pmid:25637602

9. Sucharov CC, Kao DP, Port JD, Karimpour-Fard A, Quaife RA, Minobe WA, et al. Myocardial microRNAs associated with reverse remodeling in human heart failure. JCI Insight. 2017;2:1–16.

10. Leary PJ, Kronmal RA, Bluemke DA, Buttrick PM, Jones KL, Kao DP, et al. Histamine H2 receptor polymorphisms, myocardial transcripts and heart failure, From The Multi-Ethnic Study of Atherosclerosis and Beta-blocker Effect on Remodeling and Gene Expression Trial. Am J Cardiol. 2018;121:256–261. pmid:29191567

11. Arikawa E, Sun Y, Wang J, Zhou Q, Ning B, Dial SL et al. Cross-platform comparison of SYBR Green real-time PCR with TaqMan PCR, microarrays and other gene expression measurement technologies evaluated in the Microarray Quality Control (MAQC) study. BMC Genomics. 2008;9:328. pmid:18620571

12. Everaert C, Luypaert M, Maag JLV, Cheng QX, Dinger ME, Hellemans J, et al. Benchmarking of RNA-sequencing analysis workflows using whole-transcriptome RT-qPCR expression data. Sci Rep. 2017;7:1559. pmid:28484260

13. Peters TJ, French HJ, Bradford ST, Pidsley R, Stirzaker C, Varinli H, et al. Evaluation of cross-platform and interlaboratory concordance via consensus modelling of genomic measurements. Bioinformatics. 2019;35:560–570. pmid:30084929

14. Lucas SM, Heath EI. Current challenges in development of differentially expressed and prognostic prostate cancer biomarkers. Prostate Cancer. 2012;2012:640968. pmid:22970379

15. Zhang M, Yao C, Guo Z, Zou J, Zhang L, Xiao H, et al. Apparently low reproducibility of true differential expression discoveries in microarray studies. Bioinformatics. 2008;24:2057–2063. pmid:18632747

16. Tumor Analysis Best Practices Working Group. Nat Rev Genet. 2004;5:229–237. pmid:14970825

17. Detours V, Dumont JE, Bersini H, Maenhaut C. Integration and cross-validation of high-throughput gene expression data: comparing heterogeneous data sets. FEBS Lett. 2003;546:98–102. pmid:12829243

18. Siriwardhana C, Datta S, Datta S. Inter-platform concordance of gene expression data for the prediction of chemical mode of action. Biol Direct. 2016;11:67. pmid:27993158

19. Zhang XF, Ou-Yang L, Zhao XM, Yan H. Differential network analysis from cross-platform gene expression data. Sci Rep. 2016,6:34112. pmid:27677586

20. Hochberg Y, Benjamini Y. More powerful procedures for multiple significance testing. Statistics in Medicine. 1990;9:811–818. pmid:2218183

21. R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. 2018; URL https://www.R-project.org/.

22. Wang C, Gong B, Bushel PR, Thierry-Mieg J, Xu J, Fang H, et al. The concordance between RNA-seq and microarray data depends on chemical treatment and transcript abundance. Nat Biotechnol. 2014;32:926–932. pmid:25150839

23. Shi L, Tong W, Fang H, Scherf U, Han J, Puri RK, et al. Cross-platform comparability of microarray technology: intra-platform consistency and appropriate data analysis procedures are essential. BMC Bioinformatics. 2005;6 Suppl 2:S12.

24. Krämer A, Green J, Pollard J, Tugendreich S. Causal analysis approaches in ingenuity pathway analysis. Bioinformatics. 2014;30:523–530. pmid:24336805

25. Barter R, Yu B. Superheat: A graphical tool for exploring complex datasets using heatmaps. R package version 1.0.0.

26. Engelhardt S, Hein L, Keller U, Klämbt K, Lohse MJ. Inhibition of Na(+)-H(+) exchange prevents hypertrophy, fibrosis, and heart failure in beta(1)-adrenergic receptor transgenic mice. Circ Res. 2002;90:814–819. pmid:11964375

27. Bristow MR. Treatment of chronic heart failure with β-adrenergic receptor antagonists: a convergence of receptor pharmacology and clinical cardiology. Circ Res. 2011;109:1176–1194. pmid:22034480

28. Tilley DG. G protein-dependent and G protein-independent signaling pathways and their impact on cardiac function. Circ Res. 2011;109:217–230. pmid:21737817

29. Radhika V, Hee Ha J, Jayaraman M, Tsim ST, Dhanasekaran N. Mitogenic signaling by lysophosphatidic acid (LPA) involves Galpha12. Oncogene. 2005;24:4597–4603. pmid:15856019

30. Nishida M, Tanabe S, Maruyama Y, Mangmool S, Urama K, Nagamatsu Y, et al. G alpha 12/13- and reactive oxygen species-dependent activation of c-Jun NH2-terminal kinase and p38 mitogen-activated protein kinase by angiotensin receptor stimulation in rat neonatal cardiomyocytes. J Biol Chem. 2005;280:18434–18441. pmid:15743761

31. Kilts JD, Lin SS, Lowe JE, Kwatra MM. Selective activation of human atrial Galpha12 and Galpha13 by Galphaq-coupled angiotensin and endothelin receptors. J Cardiovasc Pharmacol. 2007;50:299–303. pmid:17878759

32. Ali-Seyed M, Laycock N, Karanam S, Xiao W, Blair ET, Moreno CS. Cross-platform expression profiling demonstrates that SV40 small tumor antigen activates Notch, Hedgehog, and Wnt signaling in human cells. BMC Cancer. 2006;6:54. pmid:16522205

33. Lin SH, Beane L, Chasse D, Zhu KW, Mathey-Prevot B, Chang JT. Cross-platform prediction of gene expression signatures. PLoS One. 2013;8:e79228. pmid:24244455

34. Siriwardhana C, Datta S, Datta S. Inter-platform concordance of gene expression data for the prediction of chemical mode of action. Biol Direct. 2016;11:67. pmid:27993158

35. Heilbrunn SM, Shah P, Bristow MR, Valantine HA, Ginsburg R, Fowler MB. Increased beta-receptor density and improved hemodynamic response to catecholamine stimulation during long-term metoprolol therapy in heart failure from dilated cardiomyopathy. Circulation. 1989;79:483–490. pmid:2537158

36. Iaccarino G, Tomhave ED, Lefkowitz RJ, Koch WJ. Reciprocal in vivo regulation of myocardial G protein-coupled receptor kinase expression by beta-adrenergic receptor stimulation and blockade. Circulation. 1998;98:1783–1789. pmid:9788834

37. Kranias EG, Bers DM. Calcium and cardiomyopathies. Subcell Biochem. 2007;45:523–537. pmid:18193651

38. Sweet ME, Cocciolo A, Slavov D, Jones KL, Sweet JR, Graw SL, et al. Transcriptome analysis of human heart failure reveals dysregulated cell adhesion in dilated cardiomyopathy and activated immune pathways in ischemic heart failure. BMC Genomics. 2018;19:812. pmid:30419824

39. Margulies KB, Matiwala S, Cornejo C, Olsen H, Craven WA, Bednarik D. Mixed messages: transcription patterns in failing and recovering human myocardium. Circ Res. 2005;96:592–599. pmid:15718504

40. Oudit GY, Crackower MA, Eriksson U, Sarao R, Kozieradzki I, Sasaki T, et al. Phosphoinositide 3-kinase gamma-deficient mice are protected from isoproterenol-induced heart failure. Circulation. 2003;108:2147–2152. pmid:12963636

41. Perrino C, Naga Prasad SV, Patel M, Wolf MJ, Rockman HA. Targeted inhibition of beta-adrenergic receptor kinase-1-associated phosphoinositide-3 kinase activity preserves beta-adrenergic receptor signaling and prolongs survival in heart failure induced by calsequestrin overexpression. J Am Coll Cardiol. 2005;45:1862–1870. pmid:15936620

42. Dent MR, Dhalla NS, Tappia PS. Phospholipase C gene expression, protein content, and activities in cardiac hypertrophy and heart failure due to volume overload. Am J Physiol Heart Circ Physiol. 2004;287:H719–27. pmid:15072958

43. McMullen JR, Amirahmadi F, Woodcock EA, Schinke-Braun M, Bouwman RD, Hewitt KA, Mollica JP, et al. Protective effects of exercise and phosphoinositide 3-kinase(p110alpha) signaling in dilated and hypertrophic cardiomyopathy. Proc Natl Acad Sci. 2007;104:612–7. pmid:17202264

44. Lin RC, Weeks KL, Gao XM, Williams RB, Bernardo BC, Kiriazis H, et al. PI3K(p110 alpha) protects against myocardial infarction-induced heart failure: identification of PI3K-regulated miRNA and mRNA. Arterioscler Thromb Vasc Biol. 2010;30:724–32. pmid:20237330

45. Ghigo A, Laffargue M, Li M, Hirsch E. PI3K and calcium signaling in cardiovascular disease. Circ Res. 2017;121:282–92. pmid:28729453

46. Leblais V, Jo SH, Chakir K, Maltsev V, Zheng M, Crow MT, et al. Phosphatidylinositol 3-kinase offsets cAMP-mediated positive inotropic effect via inhibiting Ca2+ influx in cardiomyocytes. Circ Res. 2004;95:1183–90. pmid:15539636

47. Engelhardt S, Hein L, Wiesmann F, Lohse MJ. Progressive hypertrophy and heart failure in beta1-adrenergic receptor transgenic mice. Proc Natl Acad Sci. 1999;96:7059–64. pmid:10359838

48. Bisognano JD, Weinberger HD, Bohlmeyer TJ, Pende A, Raynolds MV, Sastravaha A, et al. Myocardial-directed overexpression of the human β1-adrenergic receptor in transgenic mice. J Mol Cell Cardiol. 2000;32:817–30. pmid:10775486

49. Seeland U, Selejan S, Müller P, Lohse MJ, Böhm M. Interstitial remodeling in beta1-adrenergic receptor transgenic mice. Basic Res Cardiol. 2007;102:183–93. pmid:17122889

50. Dockstader K, Nunley K, Medway AM, Nelson P, Liggett S, Bristow MR, et al. Analysis of temporal expression of mRNAs and miRNAs in Arg- and Gly389 polymorphic variants of the β1-adrenergic receptor. Physiol Genomics. 2011;43:1294–1306 pmid:21954455

51. Kilts JD, Lin SS, Lowe JE, Kwatra MM Selective activation of human atrial Galpha12 and Galpha13 by Galphaq-coupled angiotensin and endothelin receptors. J Cardiovasc Pharmacol. 2007;50:299–303. pmid:17878759

52. Kong SW, Bodyak N, Yue P, Liu Z, Brown J, Izumo S, Kang PM. Genetic expression profiles during physiological and pathological cardiac hypertrophy and heart failure in rats. Physiol Genomics. 2005;21:34–42. pmid:15623566

53. Karliner JS, Honbo N, Summers K, Gray MO, Goetzl EJ. The lysophospholipids sphingosine-1-phosphate and lysophosphatidic acid enhance survival during hypoxia in neonatal rat cardiac myocytes. J Mol Cell Cardiol. 2001;33:1713–7. pmid:11549349

54. Lips DJ, Bueno OF, Wilkins BJ, Purcell NH, Kaiser RA, Lorenz JN, al. MEK1-ERK2 signaling pathway protects myocardium from ischemic injury in vivo. Circulation. 2004;109(16):1938–41. pmid:15096454

55. McNally EM, Buttrick PM, Leinwand LA. Ventricular myosin light chain 1 is developmentally regulated and does not change in hypertension. Nucleic Acids Res. 1989;17:2753–67. pmid:2717409

56. Morano I, Ritter O, Bonz A, Timek T, Vahl CF, Michel G. Myosin light chain-actin interaction regulates cardiac contractility. Circ Res. 1995;76:720–5. pmid:7728988

57. Vehlow C, Kao DP, Bristow MR, Hunter L, Weiskopf D, Görg C. Visual analysis of biological data-knowledge networks. BMC Bioinformatics. 2015;16:135. pmid:25925016

58. Lee NH, Saeed AI. Microarrays: an overview. Methods Mol Biol. 2007;353:265–300. pmid:17332646

59. Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009;10:57–63. pmid:19015660

60. Pawitan Y, Michiels S, Koscielny S, Gusnanto A, Ploner A. False discovery rate, sensitivity and sample size for microarray studies. Bioinformatics. 2005;21:3017–3024. pmid:15840707

61. Sweeney TE, Haynes WA, Vallania F, Ioannidis JP, Khatri P. Methods to increase reproducibility in differential gene expression via meta-analysis. Nucleic Acids Res. 2017;45:e1. pmid:27634930

62. Singer JD, Willett JB (2003). A Framework for Investigating Change over Time. In Applied Longitudinal Data Analysis: Modeling Change and Event Occurrence. Oxford University Press.

63. Lowes BD, Minobe WA, Abraham WT, Rizeq MN, Bohlmeyer TJ, Quaife RA, et al. Changes in gene expression in the intact human heart: down-regulation of alpha-myosin heavy chain in hypertrophied, failing ventricular myocardium. J Clin Invest. 1997;100:2315–2324, 1997. pmid:9410910