Transcriptomic Analysis of Inflammatory Cardiomyopathy Identifies Molecular Signatures of Disease and Informs in silico Prediction of a Network-Based Rationale for Therapy

Inflammatory cardiomyopathy covers a group of diseases characterized by inflammation and dysfunction of the heart muscle. The immunosuppressive agents such as prednisolone, azathioprine and cyclosporine are modestly effective treatments, but a molecular rationale underpinning such therapy or the development of new therapeutic strategies is lacking. We aimed to develop a network-based approach to identify therapeutic targets for inflammatory cardiomyopathy from the evolving myocardial transcriptome in a mouse model of the disease. We performed bulk RNA sequencing of hearts at early, mid and late time points from mice with experimental autoimmune myocarditis. We identified a cascade of pathway-level events involving early activation of cytokine and chemokine-signaling pathways that precede leucocyte infiltration and are followed by innate immune, antigen-presentation, complement and cell-adhesion pathway activation. We integrated these pathway events into a network-like representation from which we further identified a 50-gene subnetwork that is predominantly induced during the course of autoimmune myocardial inflammation. We developed a combinatorial attack strategy where we quantify network tolerance to combinatorial node removal to determine target-specific therapeutic potential. We find that combinatorial attack of Traf2, Nfkb1, Rac1, and Vav1 disconnects 80% of nodes from the largest network component. Two of these nodes, Nfkb1 and Rac1, are directly targeted by prednisolone and azathioprine respectively, supporting the idea that the methodology developed here can identify valid therapeutic targets. Whereas Nfkb1 and Rac1 removal disconnects 56% of nodes, we show that additional removal of Btk and Pik3cd causes 72% node disconnection. In conclusion, transcriptome profiling, pathway integration, and network identification of autoimmune myocardial inflammation provide a molecular signature applicable to the diagnosis of inflammatory cardiomyopathy. Combinatorial attack provides a rationale for immunosuppressive therapy of inflammatory cardiomyopathy and provides an in silico prediction that the approved therapeutics, ibrutinib and idelalisib targeting Btk and Pik3cd respectively, could potentially be re-purposed as adjuncts to immunosuppression.

Introduction

Inflammatory cardiomyopathy covers a group of diseases characterized by inflammation and dysfunction of the heart muscle, and often progresses to heart failure and death [reviewed in (1)]. Inflammatory cardiomyopathy can be triggered by infections (most commonly viral in the developed world), immune-mediated mechanisms and toxins (1, 2). Immunosuppressive therapy using agents such as prednisolone, cyclosporine and azathioprine is modestly effective in chronic virus-negative, giant cell and active autoimmune forms of inflammatory cardiomyopathy, and is recommended for such patients (2). A key challenge is to understand the molecular rationale underlying current immunosuppressive therapy for inflammatory cardiomyopathy, and to develop new therapeutic strategies that would improve patient outcomes.

Rodent models of autoimmune myocarditis have been extensively studied to understand the mechanisms of inflammatory cardiomyopathy [reviewed in (3)]. Cardiac inflammation, typically manifested by infiltration of myeloid and CD4⁺ T cells, occurs between 14 and 21 days after immunization, following which there is progressive fibrosis modeling clinical cardiomyopathy (3). A number of genes have been shown to be necessary for the development of inflammatory cardiomyopathy in this model and include the TLR-signaling adaptor Myd88 (4), cytokines Tnf (5), Il23a (6), Csf1 (7), Csf2 (8), Il6 (9), Il1r1 (9), Il17a (10), and Il4 (11), chemokines and their receptors Ccl2 (12), Ccl3 (12), Ccr1 (13), Ccr2 (14), and Ccr5 (12), and complement C3 and the complement receptor Cr2 (15). There is, however, no clear understanding of how these genes, and the pathways they function in, are interconnected, and of the most important drivers of myocardial inflammation. Neither is there a clear understanding of the molecular mechanism of action of immunosuppressive agents in clinical inflammatory cardiomyopathy.

We hypothesized that identifying the gene network activated during the course of murine experimental autoimmune myocarditis, a model of inflammatory cardiomyopathy, would provide a molecular rationale for informing new therapeutic strategies. In this study, we used sequential transcriptome profiling to identify the genetic network underlying inflammatory cardiomyopathy in a mouse model. We identified critical nodes that provide a molecular rationale for current immunosuppressive therapy, and also identified opportunities for future therapy including drugs that could be re-purposed for the treatment of inflammatory cardiomyopathy.

Materials and Methods

Mice

BALB/cAnNCrl mice (male, 6–8 weeks old), were purchased from Charles River UK and acclimatized for a week. Male mice were studied as they have a higher severity of myocarditis, reducing animal numbers needed for study (16). Animals were group housed in humidity and temperature-controlled rooms on a 12-h light-dark cycle, and provided food and water ad libitum. Mice (8–10 weeks) were immunized as described (17, 18). Peptide Ac-RSLKLMATLFSTYASADR-OH (19) (MyHC-α_614−629) was obtained at 95% purity from China Peptides. The peptide (1 mg) was dissolved in PBS (1 ml, containing 2%DMSO). Supplemented Complete Freund's Adjuvant (SCFA) was produced by adding 100 mg of powdered heat-killed M. tuberculosis strain H37Ra (MT) (BD, Cat # 231141-6X100 mg) to 25 mL Complete Freund's Adjuvant (Sigma Cat# F5881). SCFA (990 μL) and peptide (1,010 μL) were emulsified using 2.5 mL glass syringes. Control emulsions were prepared by emulsifying SCFA with PBS. Mice were immunized on day 0 and day 7 by injection of 200 μL SCFA + peptide/PBS emulsion (100 μg peptide). Control mice were injected on day 0 and day 7 with 200 μL SCFA + PBS emulsion. Injections of SCFA were performed under isoflurane inhalation anesthesia. All mice also received pertussis toxin (Sigma Cat# P7208, 500 ng/mouse, dissolved in 100 μL PBS) intraperitoneally on day 0. Mice were culled on day 10, 15 and 21 by anesthesia (4% isoflurane inhalation administered in a chamber) followed by cervical dislocation, and hearts were harvested. A portion of the heart was snap frozen in RNAlater (Ambion, Cat#AM7021) in liquid nitrogen and stored at −80°C for subsequent RNA extraction.

Cardiac Troponin Measurement

Serum levels of cardiac troponin-I were measured using an iSTAT1 (Abaxis) handheld analyzer and cardiac troponin cartridges (Abaxis) following the manufacturer's instructions.

Flow Cytometry

Cardiac tissue was digested with collagenase I, XI and hyaluronidase at 37°C for 1 h, and triturated through a nylon mesh to isolate a single cell population as described (14). Red blood cell lysis was performed by adding 10 mL of 1 × red blood cell lysis solution (ACK buffer) with 2 mM EDTA and incubated at room temperature for 3 min. Cells washed with HBSS supplemented with 0.2%(w/v) BSA and 1%(w/v) FCS, centrifuged at 4°C for 10 min at 300 g. Cells were then incubated in blocking solution containing 5% normal mouse serum, 5% normal rat serum, and 1% FcBlock (Cat#553142, BD Biosciences), in PBS on ice for 30 min. Cells were stained with 1:100 dilution of PE anti-mouse CD45 antibody clone 30-F11 PE conjugated (Biolegend, Cat#103106) for 30 min on ice. Cells were washed as described above and analyzed on an Attune Nxt Flow Cytometer with software version 3.1.2. Flow cytometry data were analyzed using FlowJo v10.7.1. The gating strategy used is described in Supplementary Figure 2.

RNA Extraction

Myocardial RNA was isolated from snap frozen myocardial samples using RNeasy Mini Kit (Qiagen, Cat#74104) following the manufacturer's instructions. Myocardial samples for RNA extraction were chosen to as far as possible satisfy the diagnostic criterion of myocarditis (presence of inflammation) as this may not be 100% penetrant in the BALB/c mouse experimental autoimmune myocarditis (EAM) model (20). Samples used are indicated in Supplementary Figure 1. Total RNA quantity and integrity were assessed using Quant-IT RiboGreen RNA Assay Kit (Invitrogen, Carlsbad, CA, USA) and Agilent Tapestation 2200 RNA Screentape. Purification of mRNA, generation of double stranded cDNA and library construction were performed using NEBNext Ultra II Directional RNA Library Prep Kit for Illumina (New England Biolabs) with minor modifications to manufacturer specifications. Amplified libraries were analyzed for size distribution using the Agilent Tapestation 2200 D1000. Libraries were quantified using Picogreen and relative volumes were pooled accordingly. Sequencing was performed as 75 bp paired end reads on a HiSeq4000 according to Illumina specifications.

RNA-seq Data Processing and Analysis

Samples were mapped to the GRCm38 genome using HISAT2 (https://ccb.jhu.edu/software/hisat2/index.shtml). Duplicate reads were identified and removed using Picard (https://broadinstitute.github.io/picard/). The mapped reads were then assigned to Ensembl genomic features defined in GRCm38.69 (http://primerseq.sourceforge.net/gtf.html) using feature Counts (http://subread.sourceforge.net/). Analysis was performed using edgeR (version 3.28.1) (https://www.bioconductor.org/packages/release/bioc/html/edgeR.html) and limma (version 3.42.2) packages for R (version 3.6.1). Genes with low read counts across all samples were removed. “Low” read counts were defined as not meeting the equivalent of 10 reads before normalization. If the number of samples where that threshold was met was at least as large as the smallest experimental group (4), then that gene was kept; otherwise the gene was excluded. For each timepoint (Day 10, Day 15, Day 21), the disease and control groups were compared to identify differentially expressed genes. The Benjamini–Hochberg correction was applied to the raw p-values to control the false discovery rate (FDR) in the final results.

Quality Check of RNA Sequences

The Ensembl annotations used for read assignment contained 41,388 genes. We observed at least one read assigned to 27,429 of them across the entire set of data. Under the rules for excluding marginally-expressed genes described above, the number of genes that remain varies by the comparison considered. At 10 days, 14,737 remain, at 15 days, 15,318 and at 21 days, 15,300. The assignment of reads to features for all the (non-outlier) samples are shown in Supplementary Figure 7. Different filters for different time point comparisons—genes with zero count, genes with counts too low in too many samples, and finally the genes taken forward for the analysis are presented in Supplementary Figure 7B.

Differential Expression, Gene Groupings and Pathway Analysis

Differentially expressed genes between myocarditis and control mice at each time point were defined based on FDR < 0.01 and at least two-fold change in expression level (either induced or repressed). Genes differentially expressed at each time point were compared to each other identifying seven major gene groups (see results), visualized using the UpSetR package (21). Gene groups were subjected to pathway enrichment analysis using the XGR package (22). Pathways and member genes (mouse) were obtained from KEGG (accessed on June 2020). The one-sided Fisher's exact test was used to test the enrichment of pathways for each of gene groups. We carried out such pathway enrichment analysis separately for KEGG Organismal Systems pathways and Environmental Information Processing pathways.

Network Analysis for Differential Genes Over Time

We extended our previous algorithm (heuristically solving a prize-collecting Steiner tree problem) (23) to identify a subnetwork through integrative analysis of time-course expression data and KEGG-merged network interaction data. The algorithm takes two inputs. The first is the score for genes quantifying the expression importance over time, this is, –log₁₀(FDR aggregated over 3 time points). The aggregation is based on Fisher's combined method implemented in the dnet package (23); in other words, a gene receives a higher score if it is more significantly changed over time. The second input is the gene interaction network (the parent network) merged from 11 KEGG pathways (illustrated in Figures 1C,D). The output is a gene subnetwork containing highly scored genes that are linked together through a few lowly scored genes as linkers. The subnetwork is highlighted in the parent network for which the layout is rendered using stress majorization (24).

FIGURE 1

Figure 1. Patterns of gene expression and enrichment of KEGG pathways in the heart over time following induction of experimental autoimmune myocarditis. (A) Bar plot showing the number of genes in each of seven gene groups that are classified by onset and persistence of expression. EPi, early-persistent induced; MPi, mid-persistent induced; MPr, mid-persistent repressed; MTi mid-transient induced; MTr, mid-transient repressed; LTi, late-transient induced; LTr, late-transient repressed. D, days from induction of EAM. (B) Molecular pathways significantly enriched for each of gene groups based on adjusted p < 0.05 (color-coded). The size of each individual dot is proportional to the number of member genes. Pathways are sourced from KEGG maps, and are displayed separately for Environmental Information Processing pathways (top) and Organismal Systems pathways (bottom). (C), Plot of fraction of induced genes (Y-axis) within the 4 pathways identified from the EPi group of genes over time (X-axis). (D) Plot of fraction of induced genes (Y-axis) within the 7 pathways identified from the MPi group of genes over time (X-axis).

Attackness and Combinatorial Attack

We used “attackness” to quantify the tolerance of the network to individual node removal, defined as the fraction of network nodes disconnected from the giant component after node removal. The attackness ranges from 0 to 1, with the higher value indicating the more vulnerable (critical) node for the network. Similarly, we introduced a new concept “combinatorial attack” to maximize the attackness for nodes in a specific combination removed, for example, identifying the optimal combination involving any two nodes.

Drug Repurposing Based on ChEMBL

We extracted current drug therapeutics, target genes and disease indications from the ChEMBL database (25) (v27), considering: (1) non-withdrawn drugs; (2) non-promiscuous therapeutic target genes with well-defined mechanism of action and explaining the efficacy of drugs in disease; and (3) drugs with the maximum phase for a gene given a disease indication.

Statistical Analysis

Summary statistics were calculated in GraphPadPrism, and statistical significance of differences between groups were evaluated using a one-way ANOVA with Dunnett's correction for multiple comparisons or using an unpaired t-test as indicated. The Benjamini–Hochberg correction was applied to the raw p-values to control the false discovery rate (FDR) in the RNASeq results.

Ethics Declarations

All animal procedures were approved locally by the University of Oxford Animal Welfare and Ethical Review Body and by the UK Home Office and carried out in accordance with the UK Animals (Scientific Procedures) Act 1986, under project license PPL P973A60F5. All procedures conformed to guidelines from Directive 2010/63/EU of the European Parliament on the protection of animals used for scientific purposes.

Data Availability

The datasets supporting the conclusions of this article are available in the Gene Expression Omnibus repository, with accession number GSE155423. Codes and associated data are provided as Supplementary Material.

Results

Experimental Autoimmune Myocarditis

We sought to determine how the cardiac transcriptome changes during autoimmune myocardial inflammation using an established model of experimental autoimmune myocarditis (EAM) in which heart-specific antigens are used to trigger myocarditis (3, 18). We induced EAM in mice by immunizing them with myosin heavy chain peptide and adjuvants (Complete Freund's Adjuvant and pertussis toxin, see schematic in Supplementary Figure 1A and methods). Control mice were generated by substituting PBS vehicle for myosin heavy chain peptide. We performed flow cytometry analysis (FACS) of single cell suspensions from hearts isolated at day 21 following immunization to determine the proportion of CD45⁺ (total leucocyte) cells as a quantitative measure of inflammation, in these mice, and also in naïve mice. We found that there was a significant increase in myocardial CD45⁺ cells in EAM mouse hearts (Supplementary Figures 1B, 2). All EAM hearts had elevated levels of CD45⁺ cells exceeding the 95% confidence interval of either the control or the naïve group, and 4 of 6 EAM mice were chosen for RNA sequencing. There was no significant difference in CD45⁺ cells between naïve and control mice (Supplementary Figure 1B). We next analyzed EAM mice at day 15 and day 10, again by FACS of single cell suspensions from hearts (Supplementary Figures 1C,D), and also by measuring cardiac troponin I (cTnI) (Supplementary Figures 1E,F). There was a statistically significant difference in myocardial CD45⁺ cells at day 15 between EAM and control mice (Figure 1C), and also a significant elevation of cTnI, indicative of myocardial damage (Supplementary Figure 1E). None of the EAM mice at day 10 showed elevated levels of either myocardial CD45% cells or cardiac troponin I (Supplementary Figures 1D,F). At each time point 4 mice in EAM and control groups were chosen for RNA sequencing.

RNA Sequencing Shows That Major Inflammatory Pathways Are Activated Early in the Course of Autoimmune Myocardial Inflammation

We identified genes that were differentially expressed between immunized and control groups at each time point using a false discovery rate (FDR) threshold of <0.01 and at least two-fold changes (either induced or repressed). Based on their temporal profiles and degree of persistence of expression over time, we then categorized these genes into seven major groups (Figure 1A, Supplementary Figure 3). For each gene group we next performed enrichment analysis using well-curated pathways based on the KEGG map (26) (Figure 1B). Twenty-two significantly enriched pathways (adjusted p < 0.05) were identified, with 11 pathways being identified from groups with persistent gene induction (Supplementary Figures 4, 5). Of these 11 pathways, four (cytokine-cytokine receptor, viral protein interaction with cytokine–cytokine receptor, TNF and chemokine signaling pathways) are activated early in the disease process (day 10) at a time point where there is no evidence of leucocyte infiltration (Figure 1C). There is substantial further activation of these four pathways at day 15 followed by a decline in activity by day 21, despite increased numbers of infiltrating leukocytes. Seven pathways are predominantly activated at day 15 (Figure 1D). These include antigen processing and presentation, cell adhesion molecules, complement and coagulation, NK cell mediated cytotoxicity, NF-kB signaling, NOD-like signaling and osteoclast differentiation. These results indicate that inflammatory pathways are activated early in the course of autoimmune myocardial inflammation, and molecular changes may precede the onset of leucocyte infiltration.

Integration of Gene Expression Changes With KEGG Molecular Interactions Identifies a 50-Gene Network in Autoimmune Myocardial Inflammation

Although the 11 pathways identified above are informative individually, it is unclear how genes from these pathways are networked as a whole, and which part of the network formed by these pathways is best reflective of our context-specific data (i.e., time-course expression data generated in this study). We reasoned that exploring the topology of this network in an unbiased manner could better understand disease mechanism. We therefore performed an integrated analysis of KEGG molecular interactions with gene expression changes at different time points (see online methods), identifying a subnetwork of 50 gene nodes based on their expression importance over time (Figure 2A). The cascade of subnetwork activation through the course of the disease process in this model of myocarditis is visualized in Figure 2B, where fold change over time and differential significance of individual genes are indicated respectively by color and size of individual gene nodes. Notably, interacting genes in this subnetwork are all induced with the changes peaking at the intermediate time point (day 15).

FIGURE 2

Figure 2. Gene subnetwork identification and temporal evolution in murine myocarditis model. (A), Overview of the gene subnetwork, with 50 gene nodes labeled by gene symbols and 75 connections between nodes (network edges). This subnetwork was identified from a parent network (in gray, obtained by merging 11 KEGG pathways identified in Figure 1) by scoring genes for their expression importance over time in EAM. Genes known to be required for the development of murine myocarditis are indicated in boxes. (B), Gene subnetwork illustrated at three different time-points with the same layout as shown in (A). Nodes at each time point are colored by log₂(fold change) and sized by –log₁₀(FDR) as indicated in the scale.

Modest Network Disconnection Is Achieved by Single Node Removal

One way of identifying the importance of a node is to quantify the tolerance of the network to node removal based on well-established percolation theory [reviewed in (27)]. A measure of overall network connectivity, and hence function, is the size of the largest connected, or “giant” component (28). The effect of node removal can be defined in terms of the fraction of network nodes disconnected from the giant component after node removal, a measure that we call “attackness” for simplicity. According to the percolation process, removal of a node with high attackness results in a large proportion of disconnected nodes. To assess the potential utility of this approach, we analyzed the 50-gene subnetwork to calculate attackness for each node (Figure 3A). We find that Traf2 and Nfkb1 have the highest degree of “attackness” (0.5 and 0.42, respectively).

FIGURE 3

Figure 3. Individual and combinatorial attack to identify vulnerable subnetwork nodes. (A) Individual attack strategy. Top panel: Plot of attackness (Y-axis) by node (X-axis, nodes ordered as in bottom panel). The “attackness” metric quantifies the susceptibility of the subnetwork to the individual node removed, and is defined as the fraction of nodes disconnected from the giant component after node removal. Bottom panel: Nodes in the 50-node subnetwork are illustrated as colored circles in descending order of attackness. Inserted panel (top): Graph of the 50-node subnetwork with the three nodes with highest attackness values indicated in pink. (B) Combinatorial attack strategy to identify node combinations that maximize network disconnection. Top panel: Plot of combinatorial attackness (Y-axis) by node combinations (X-axis, node combinations indicated in bottom panel). Bottom panel: Nodes removed in combination are indicated by connected colored circles. (C) Graph of the 50-node subnetwork with the nine nodes used most frequently in combinatorial attack indicated as pink circles.

Maximal Network Disconnection Is Achieved by Removal of Traf2, Nfkb1, Rac1, and Vav1

Most diseases, however, typically require a polypharmacological approach, where multiple targets are attacked (29). To minimize on-target adverse effects it is important to identify the smallest possible number of nodes that need to be targeted in order to have an adequate effect. To address this, we removed between two to five nodes in combination and assessed the effect of such removal on the fraction of network nodes disconnected from the giant component (i.e., the combinatorial attackness). The top-ranking combinations and their corresponding combinatorial attackness are shown in Figure 3B, with genes frequently involved in these top combinations highlighted in the subnetwork (Figure 3C). Removal of four nodes (Traf2, Nfkb1, Rac1, and Vav1) is predicted to disconnect 80% of nodes from the giant component. Combined attack on Nfkb1 and Rac1 is predicted to disconnect 56% nodes in the 50-node network.

Combinatorial Attack Identifies Additional Nodes Targeted by Approved Therapeutics

Drug repurposing is a strategy that repositions existing effective therapeutics with acceptable safety profiles to new indications (30). We used the ChEMBL database (25) to explore the evidence supporting potential therapeutics targeting products of the human gene equivalents of the mouse subnetwork. We identified 11 genes (BTK, CCR5, HCK, IL1B, ITGA4, ITGAL, ITGB2, ITGB7, LCK, PIK3CD, and PIK3CG) that could be targeted by drugs already in clinical use (Figure 4A). We also identified 7 preclinical (phased) drug targets. These include two genes (BCL2A1 and CCR2) targeted by drugs in clinical phase III, and five genes (BIRC3, CASP1, CCR1, CXCL10, and PRKCD) targeted by drugs in clinical phase II. We found the highest degree of support from clinical evidence (approved therapeutics), and also observed enrichment for preclinical (phased) drug targets, though at a less significant level (Figure 4B). We next calculated the attackness for these targets and found the lowest degree for all of them except for Lck (Figure 4C). This observation motivated us to explore the effect of a combinatorial attack strategy on the subnetwork. The top two-node, three-node and four-node combinations and their corresponding combinatorial attackness are shown in Figure 4D. We find that the tyrosine-protein kinase LCK inhibitor can be combined with any other inhibitor to achieve a small improvement in effect. A further improvement in effect may be observed by combining C–C chemokine receptor type 1 (CCR1) antagonist (Phase II), C–C chemokine receptor type 2 (CCR2) antagonist (Phase III), and C–C chemokine receptor type 5 (CCR5) antagonist (approved). The effect of this triplet combination can be further enhanced with addition of the LCK inhibitor.

FIGURE 4

Figure 4. Individual and combinatorial attack strategies to identify re-purposing candidates. (A) Dot plot showing 19 mouse genes in the subnetwork (X-axis, right) and their human homologs (Y-axis, left) mapped to currently approved and phase II/III drugs in human disease indications (indicated on the X-axis, and obtained from ChEMBL v25). Dots are color-coded by clinical phase status as indicated in the legend. Information on drugs and mechanism of action is indexed in the plot and detailed beneath the dot-plot. (B) Forest plot of approved or phased drug targets enriched in the subnetwork genes. The significance level (adjusted p-value), odds ratio and 95% confidence interval (represented by lines) calculated according to Fisher's exact test (one-sided). Also illustrated are gene lists overlapped. (C) Individual attack strategy. Top panel: Plot of attackness (Y-axis) by node (X-axis, nodes indicated by colored circles in the bottom panel). Bottom panel: 19 nodes within the 50-node subnetwork that have identifiable drugs. (D) Combinatorial attack strategy to identify node combinations that maximize network disconnection. Top panel: Plot of combinatorial attackness (Y-axis) by node combinations (X-axis, node combinations indicated in bottom panel). Bottom panel: Nodes removed in combination are indicated by connected colored circles. (E) Combinatorial attack strategy to identify additional node combinations that maximize network disconnection in the context of Nfkb1 and Rac1 node removal using state-of-the-art therapy (glucocorticoid and azathioprine). Top panel: Plot of combinatorial attackness (Y-axis) by node combinations (X-axis, node combinations indicated in bottom panel). Bottom panel: Nodes removed in combination are indicated by connected colored circles.

Combinatorial Attack in the Context of Established Therapy Identifies Btk and Pik3cd

As discussed above, established therapy with prednisolone and azathioprine, targeting Nfkb1 and Rac1, respectively, is predicted to disconnect 56% nodes in the 50-node network. We evaluated a repurposing strategy in the context of Nfkb1+Rac1 node removal. We find that the targeting Btk and Pik3cd in the context of established therapy will disconnect 72% of the nodes in the 50-node network, a substantial increase from which there is no further improvement by targeting another node (Figure 4E).

Discussion

The RNA sequencing studies reported here provide a resource with demonstrated utility in advancing our knowledge of the molecular pathology of autoimmune myocardial inflammation. Our work advances published studies of molecular pathology by integrating the temporal sequence of gene expression changes with well-established information on how these genes interact in pathways. Supporting our analysis, several genes within the 11 pathways identified in this study have already been shown to be critical for the development of EAM, including Ccl2 (12), Ccl3 (12), Ccr2 (14), Ccr5 (12), Il1r1 (9), Il17ra (31), and Tnf (5). The role of NK cells and complement pathways in myocarditis has also previously been elucidated (15, 32). Importantly, a number of genes and pathways are identified that are induced at early time points and persist through the course of inflammation. Our work advances the molecular pathology of autoimmune myocardial inflammation by integrating these genes and pathways into a single network and identifying a 50-gene molecular network that is induced over the time course of autoimmune myocardial inflammation. Elements within this molecular network have potential for exploitation as biomarkers or molecular signatures of disease. The diagnosis of myocarditis and inflammatory cardiomyopathy currently relies on the presence of leucocyte infiltration in endomyocardial biopsy specimens, which is the “gold-standard” (2). As such infiltrates may be patchy, the disease may go unrecognized.

The data from RNA sequencing analysis indicate that chemokine signaling pathway and cytokine–cytokine receptor interactions are elevated at D10, and increase further by D15 and at D21. While no evidence of myocardial inflammation is detected at D10, it is clearly detected at D15 and peaks by D21. Our results thus show that molecular changes in the myocardium appear before leucocyte infiltration occurs, and persist through to later stages. The delay between chemokine/cytokine gene expression and maximum accumulation of CD45⁺ cells is expected, and likely arises from the necessary subsequent stages of chemokine protein translation, post-translational modification, secretion, presentation on the endothelial luminal surface, binding to leucocyte chemokine receptor, activation of the leucocyte adhesion cascade, and leucocyte transendothelial migration that leads to tissue leucocyte accumulation (33, 34).

Alternative approaches that could be used for such molecular phenotyping include RT-qPCR, proteomic approaches such as western blotting or multiplexed immunoassay (e.g., Luminex xMAP) technologies, and single-cell RNA sequencing approaches to characterize the nature of the immune cell infiltrate. RNA sequencing is considered the “gold-standard” for whole transcriptome gene expression studies, and is highly correlated with alternative transcriptomic approaches that detect changes in expression of individual genes such as microarrays or RT-qPCR (35). Proteomic approaches such as western blotting or multiplexed immunoassays we consider impractical for the molecular characterization of multiple gene products due to the very small amounts of tissues that are typically obtained in myocardial biopsy specimens, and the requirement for highly specific and sensitive antibodies to the gene products of interest.

A limitation of our work is that the restricted amount of tissue available for FACS analyses (most tissue being used for RNA extraction), made more detailed studies of cellular subpopulations technically challenging. In future studies characterization of inflammatory subpopulations—for instance using single cell RNA sequencing—would be a desirable goal, as it would allow an understanding of the changes in cellular compartments, and also allow correlation with transcriptomic analyses. The application of single-cell RNA sequencing as a molecular diagnostic tool is indeed of great interest, and will likely provide complementary information to bulk RNA sequencing. Single-cell RNA sequencing is technically more challenging and substantially more expensive than bulk RNA sequencing, and currently may be impractical for diagnostic applications where very small amounts of tissue are typically available for study. Taken together, we suggest that molecular characterization of biopsied myocardium using bulk RNASeq may be a sensitive and cost-effective indicator of autoimmune myocardial inflammation in comparison to histopathology, and could complement histopathology in the evaluation of myocardial inflammation.

These studies also provide a molecular rationale for immunosuppressive therapy of inflammatory cardiomyopathy. Recommended immunosuppressive approaches for autoimmune and virus-negative inflammatory cardiomyopathy employ glucocorticoids, azathioprine and cyclosporine singly or in combination (1, 2). These immunosupressive approaches either directly or indirectly target three of these four critical nodes identified by the combinatorial attack strategy. The azathioprine metabolite 6-thio-GTP directly binds RAC1 to suppress VAV1 guanosine exchange activity on RAC1 and inhibits T-cell-APC conjugation, causing immunosuppression (36, 37). Glucocorticoids such as prednisolone repress NFKB activity through diverse mechanisms including increasing its cytoplasmic retention, and reducing its synthesis (38). Cyclosporine also inhibits NFKB activity through its action on the 20S proteasome (38). Combined attack on NFKB and RAC1 is predicted to disconnect 56% nodes in the 50-node network. These findings support the validity of the subnetwork definition coupled with node disconnection approach developed here to predict therapeutic targets, and suggests that the extent of network disconnection achieved predicts therapeutic response. These findings thus provide a molecular rationale for the immunosuppressive therapy recommended currently for inflammatory cardiomyopathy.

While the combinatorial attack strategy predicts critical network node/target combinations for drug therapy, it does not necessarily imply therapeutic success. Such a strategy may even result in the use of drugs with unacceptable on-target adverse effects. For instance, while TRAF2 has a key role in the subnetwork, its loss results in inflammatory and cardiac failure phenotypes which would invalidate it as a target (39, 40). Drug repurposing is a strategy that repositions existing effective therapeutics with acceptable safety profiles to new indications (30). As discussed above, established therapy with prednisolone and azathioprine, targeting Nfkb1 and Rac1, respectively, is predicted to disconnect 56% nodes in the 50-node network. Given that any new therapeutic can only be ethically tested in the context of combination with established therapy, we evaluated a repurposing strategy in the context of Nfkb1+Rac1 node removal. We find that the targeting Btk and Pik3cd in the context of established therapy will disconnect 72% of the nodes in the 50-node network, which predicts a substantial improvement in therapeutic efficacy. These in silico predictions need to be validated in appropriate studies, e.g., using genetic (e.g., knockout) or pharmacological (e.g., using agents such as Ibrutinib and Idelalisib, targeting BTK and PIK3CD, respectively) approaches in animal models of inflammatory cardiomyopathy.

In summary we have conducted a series of analyses (Supplementary Figure 6), starting with RNA sequencing to define the sequential evolution of the myocardial transcriptome during experimental autoimmune myocarditis. Combining the data with KEGG pathways allows identification of a myocarditis gene subnetwork that is induced at early, mid and late time points through the course of cardiac inflammation. We show that key vulnerable nodes in the temporally evolving subnetwork can be identified by a combinatorial attack strategy, the concept being based on well-established percolation theory. This approach is supported by the identification of the key nodes Nfkb1 and Rac1, in the myocardial inflammation subnetwork, which are well-established targets of glucocorticoids, cyclosporine and azathioprine, and which have proven efficacy in the therapy of autoimmune and virus-negative inflammatory cardiomyopathy. Our analysis provides an in silico prediction that agents such as Ibrutinib and Idelalisib, targeting BTK and PIK3CD, respectively, could be re-purposed for the therapy of virus-negative inflammatory cardiomyopathy. Combination therapy is increasingly adopted across a range of chronic immuno-inflammatory diseases with the aim of inducing sustained remission (41, 42). Our approach, combining sequential transcriptome profiling, in-depth pathway network analysis and combinatorial attack, adds a new strategy to the repository of network medicine analyses designed to aid development of drug repurposing and drug combinations (43, 44). We suggest that our approach would be applicable more broadly to other immuno-inflammatory disorders where a similar series of analyses could be applied not only in mouse models, but also in human disease.

Data Availability Statement

The datasets supporting the conclusions of this article are available in the Gene Expression Omnibus repository, with accession number GSE155423. Codes and associated data are provided as a series of online chapters at the project web interface (http://galahad.well.ox.ac.uk/Myocarditis).

Ethics Statement

All animal procedures were approved locally by the University of Oxford Animal Welfare and Ethical Review Body and by the UK Home Office and carried out in accordance with the UK Animals (Scientific Procedures) Act 1986, under project license PPL P973A60F5. All procedures conformed to guidelines from Directive 2010/63/EU of the European Parliament on the protection of animals used for scientific purposes.

Author Contributions

KS and GD performed and analyzed mouse and flow cytometry experiments. HF, JK, and SB conceived the network analysis strategy. HF designed all scripts for network analysis. BW and HL designed the RNA sequencing study and analyzed RNA sequencing data. SB, DC, RW, and JK conceived the study and obtained funding support. KS, HF, JK, and SB drafted the manuscript. We thank Dr. Ying-Jie Wang for his advice on FACS data analysis. All authors contributed to editing the manuscript and approved the final manuscript version.

Funding

This research was funded by British Heart Foundation Chair (CH/09/003/26631), Project Grant (PG/16/100/32632), Program Grant (RG/18/1/33351), and Oxford BHF Centre of Research Excellence (RE/13/1/30181) awards to SB. JK is supported by a Wellcome Trust Investigator Award (204969/Z/16/Z) and the NIHR Oxford Biomedical Research Centre. HF was supported by an EU/EFPIA Innovative Medicines Initiative (ULTRA-DD 115766). We thank the Oxford Genomics Centre at the Wellcome Centre for Human Genetics (funded by Wellcome Trust grant reference 203141/Z/16/Z) for the generation and initial processing of the sequencing data, and Dr. John Broxholme for running the RNA-Seq pipeline. Further support was received from Wellcome Trust Grant (090532/Z/09/Z) to core facilities at the Wellcome Centre for Human Genetics.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2021.640837/full#supplementary-material

References

1. Trachtenberg BH, Hare JM. Inflammatory cardiomyopathic syndromes. Circ Res. (2017) 121:803–18. doi: 10.1161/CIRCRESAHA.117.310221

CrossRef Full Text | Google Scholar

2. Caforio AL, Pankuweit S, Arbustini E, Basso C, Gimeno-Blanes J, Felix SB, et al. Current state of knowledge on aetiology, diagnosis, management, and therapy of myocarditis: a position statement of the European Society of Cardiology Working Group on Myocardial and Pericardial Diseases. Eur Heart J. (2013) 34:2636–48, 2648a–d. doi: 10.1093/eurheartj/eht210

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Blyszczuk P. Myocarditis in humans and in experimental animal models. Front Cardiovasc Med. (2019) 6:64. doi: 10.3389/fcvm.2019.00064

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Marty RR, Dirnhofer S, Mauermann N, Schweikert S, Akira S, Hunziker L, et al. MyD88 signaling controls autoimmune myocarditis induction. Circulation. (2006) 113:258–65. doi: 10.1161/CIRCULATIONAHA.105.564294

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Smith SC, Allen PM. Neutralization of endogenous tumor necrosis factor ameliorates the severity of myosin-induced myocarditis. Circ Res. (1992) 70:856–63. doi: 10.1161/01.RES.70.4.856

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Wu L, Diny NL, Ong S, Barin JG, Hou X, Rose NR, et al. Pathogenic IL-23 signaling is required to initiate GM-CSF-driven autoimmune myocarditis in mice. Eur J Immunol. (2016) 46:582–92. doi: 10.1002/eji.201545924

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Meyer IS, Goetzke CC, Kespohl M, Sauter M, Heuser A, Eckstein V. Silencing the CSF-1 axis using nanoparticle encapsulated siRNA mitigates viral and autoimmune myocarditis. Front Immunol. (2018) 9:2303. doi: 10.3389/fimmu.2018.02303

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Sonderegger I, Iezzi G, Maier R, Schmitz N, Kurrer M, Kopf M. GM-CSF mediates autoimmunity by enhancing IL-6-dependent Th17 cell development and survival. J Exp Med. (2008) 205:2281–94. doi: 10.1084/jem.20071119

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Eriksson U, Kurrer MO, Sonderegger I, Iezzi G, Tafuri A, Hunziker L, et al. Activation of dendritic cells through the interleukin 1 receptor 1 is critical for the induction of autoimmune myocarditis. J Exp Med. (2003) 197:323–31. doi: 10.1084/jem.20021788

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Baldeviano GC, Barin JG, Talor MV, Srinivasan S, Bedja D, Zheng D, et al. Interleukin-17A is dispensable for myocarditis but essential for the progression to dilated cardiomyopathy. Circ Res. (2010) 106:1646–55. doi: 10.1161/CIRCRESAHA.109.213157

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Diny NL, Baldeviano GC, Talor MV, Barin JG, Ong S, Bedja D, et al. Eosinophil-derived IL-4 drives progression of myocarditis to inflammatory dilated cardiomyopathy. J Exp Med. (2017) 214:943–57. doi: 10.1084/jem.20161702

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Goser S, Ottl R, Brodner A, Dengler TJ, Torzewski J, Egashira K, et al. Critical role for monocyte chemoattractant protein-1 and macrophage inflammatory protein-1alpha in induction of experimental autoimmune myocarditis and effective anti-monocyte chemoattractant protein-1 gene therapy. Circulation. (2005) 112:3400–7. doi: 10.1161/CIRCULATIONAHA.105.572396

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Futamatsu H, Suzuki J, Koga N, Adachi S, Kosuge H, Maejima Y, et al. A CCR1 antagonist prevents the development of experimental autoimmune myocarditis in association with T cell inactivation. J Mol Cell Cardiol. (2006) 40:853–61. doi: 10.1016/j.yjmcc.2006.03.432

CrossRef Full Text | Google Scholar

14. Leuschner F, Courties G, Dutta P, Mortensen LJ, Gorbatov R, Sena B, et al. Silencing of CCR2 in myocarditis. Eur Heart J. (2015) 36:1478–88. doi: 10.1093/eurheartj/ehu225

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Kaya Z, Afanasyeva M, Wang Y, Dohmen KM, Schlichting J, Tretter T, et al. Contribution of the innate immune system to autoimmune myocarditis: a role for complement. Nat Immunol. (2001) 2:739–45. doi: 10.1038/90686

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Myers JM, Fairweather D, Huber SA, Cunningham MW. Autoimmune myocarditis, valvulitis, and cardiomyopathy. Curr Protoc Immunol. (2013) Chapter 15:Unit 15.14.1–51. doi: 10.1002/0471142735.im1514s101

CrossRef Full Text | Google Scholar

17. Ciháková D, Rose NR. Myocarditis and other immunological models of cardiac disease. In: Ardehali H, Bolli R, Losordo DW, editors. Manual of Research Techniques in Cardiovascular Medicine. New York: John Wiley and Sons, Ltd. (2013). p. 197–202. doi: 10.1002/9781118495148.ch24

CrossRef Full Text | Google Scholar

18. Fontes JA, Barin JG, Talor MV, Stickel N, Schaub J, Rose NR, et al. Complete Freund's adjuvant induces experimental autoimmune myocarditis by enhancing IL-6 production during initiation of the immune response. Immun Inflamm Dis. (2017) 5:163–76. doi: 10.1002/iid3.155

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Sonderegger I, Rohn TA, Kurrer MO, Iezzi G, Zou Y, Kastelein RA, et al. Neutralization of IL-17 by active vaccination inhibits IL-23-dependent autoimmune myocarditis. Eur J Immunol. (2006) 36:2849–56. doi: 10.1002/eji.200636484

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Pummerer CL, Luze K, Grassl G, Bachmaier K, Offner F, Burrell SK, et al. Identification of cardiac myosin peptides capable of inducing autoimmune myocarditis in BALB/c mice. J Clin Invest. (1996) 97:2057–62. doi: 10.1172/JCI118642

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Conway JR, Lex A, Gehlenborg N. UpSetR: an R package for the visualization of intersecting sets and their properties. Bioinformatics. (2017) 33:2938–40. doi: 10.1093/bioinformatics/btx364

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Fang H, Knezevic B, Burnham KL, Knight JC. XGR software for enhanced interpretation of genomic summary data, illustrated by application to immunological traits. Genome Med. (2016) 8:129. doi: 10.1186/s13073-016-0384-y

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Fang H, Gough J. The “dnet” approach promotes emerging research on cancer patient survival. Genome Med. (2014) 6:64. doi: 10.1186/PREACCEPT-1248435140128963

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Gansner ER, Koren Y, North S. Graph Drawing by Stress Majorization. Berlin, Heidelberg: Springer (2005). p. 239–50.

Google Scholar

25. Mendez D, Gaulton A, Bento AP, Chambers J, De Veij M, Felix E, et al. ChEMBL: towards direct deposition of bioassay data. Nucleic Acids Res. (2019) 47:D930–40. doi: 10.1093/nar/gky1075

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. (2000) 28:27–30. doi: 10.1093/nar/28.1.27

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Stauffer D, Aharony A. Introduction to Percolation Theory. London, Bristol, PA: Taylor and Francis (1994).

Google Scholar

28. Braunstein A, Dall'asta L, Semerjian G, Zdeborova L. Network dismantling. Proc Natl Acad Sci USA. (2016) 113:12368–73. doi: 10.1073/pnas.1605083113

CrossRef Full Text | Google Scholar

29. Hopkins AL. Network pharmacology: the next paradigm in drug discovery. Nat Chem Biol. (2008) 4:682–90. doi: 10.1038/nchembio.118

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Pushpakom S, Iorio F, Eyers PA, Escott KJ, Hopper S, Wells A, et al. Drug repurposing: progress, challenges and recommendations. Nat Rev Drug Discov. (2019) 18:41–58. doi: 10.1038/nrd.2018.168

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Wu L, Ong S, Talor MV, Barin JG, Baldeviano GC, Kass DA, et al. Cardiac fibroblasts mediate IL-17A-driven inflammatory dilated cardiomyopathy. J Exp Med. (2014) 211:1449–64. doi: 10.1084/jem.20132126

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Ong S, Ligons DL, Barin JG, Wu L, Talor MV, Diny N, et al. Natural killer cells limit cardiac inflammation and fibrosis by halting eosinophil infiltration. Am J Pathol. (2015) 185:847–61. doi: 10.1016/j.ajpath.2014.11.023

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Vanheule V, Metzemaekers M, Janssens R, Struyf S, Proost P. How post-translational modifications influence the biological activity of chemokines. Cytokine. (2018) 109:29–51. doi: 10.1016/j.cyto.2018.02.026

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Graham GJ, Handel TM, Proudfoot AEI. Leukocyte adhesion: reconceptualizing chemokine presentation by glycosaminoglycans. Trends Immunol. (2019) 40:472–81. doi: 10.1016/j.it.2019.03.009

PubMed Abstract | CrossRef Full Text | Google Scholar

35. 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. doi: 10.1038/s41598-017-01617-3

PubMed Abstract | CrossRef Full Text

36. Tiede I, Fritz G, Strand S, Poppe D, Dvorsky R, Strand D, et al. CD28-dependent Rac1 activation is the molecular target of azathioprine in primary human CD4⁺ T lymphocytes. J Clin Invest. (2003) 111:1133–45. doi: 10.1172/JCI16432

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Poppe D, Tiede I, Fritz G, Becker C, Bartsch B, Wirtz S, et al. Azathioprine suppresses ezrin-radixin-moesin-dependent T cell-APC conjugation through inhibition of Vav guanosine exchange activity on Rac proteins. J Immunol. (2006) 176:640–51. doi: 10.4049/jimmunol.176.1.640

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Yamamoto Y, Gaynor RB. Therapeutic potential of inhibition of the NF-kappaB pathway in the treatment of inflammation and cancer. J Clin Invest. (2001) 107:135–42. doi: 10.1172/JCI11914

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Guo X, Yin H, Li L, Chen Y, Li J, Doan J, et al. Cardioprotective role of tumor necrosis factor receptor-associated factor 2 by suppressing apoptosis and necroptosis. Circulation. (2017) 136:729–42. doi: 10.1161/CIRCULATIONAHA.116.026240

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Shi JH, Sun SC. Tumor necrosis factor receptor-associated factor regulation of nuclear factor kappaB and mitogen-activated protein kinase pathways. Front Immunol. (2018) 9:1849. doi: 10.3389/fimmu.2018.01849

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Colombel JF, Sandborn WJ, Reinisch W, Mantzaris GJ, Kornbluth A, Rachmilewitz D, et al. Infliximab, azathioprine, or combination therapy for Crohn's disease. N Engl J Med. (2010) 362:1383–95. doi: 10.1056/NEJMoa0904492

CrossRef Full Text | Google Scholar

42. O'dell JR, Mikuls TR, Taylor TH, Ahluwalia V, Brophy M, Warren SR, et al. Therapies for active rheumatoid arthritis after methotrexate failure. N Engl J Med. (2013) 369:307–18. doi: 10.1056/NEJMoa1303006

CrossRef Full Text | Google Scholar

43. Cheng F, Kovacs IA, Barabasi AL. Network-based prediction of drug combinations. Nat Commun. (2019) 10:1197. doi: 10.1038/s41467-019-09186-x

CrossRef Full Text

44. Nabirotchkin S, Peluffo AE, Rinaudo P, Yu J, Hajj R, Cohen D. Next-generation drug repurposing using human genetics and network biology. Curr Opin Pharmacol. (2020) 51:78–92. doi: 10.1016/j.coph.2019.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar