Dissecting the cellular landscape and transcriptome network in viral myocarditis by single-cell RNA sequencing.

Coxsackievirus B3 (CVB3)-induced myocarditis is commonly employed to study viral pathogenesis in mice. Chronically affected mice may develop dilated cardiomyopathy, which may involve the mediation of immune and nonimmune cells. To dissect this complexity, we performed single-cell RNA sequencing on heart cells from healthy and myocarditic mice, leading us to note significant proportions of myeloid cells, T cells, and fibroblasts. Although the transcriptomes of myeloid cells were mainly of M2 phenotype, the Th17 cells, CTLs, and Treg cells had signatures critical for cytotoxic functions. Fibroblasts were heterogeneous expressing genes important in fibrosis and regulation of inflammation and immune responses. The intercellular communication networks revealed unique interactions and signaling pathways in the cardiac cellulome, whereas myeloid cells and T cells had upregulated unique transcription factors modulating cardiac remodeling functions. Together, our data suggest that M2 cells, T cells, and fibroblasts may cooperatively or independently participate in the pathogenesis of viral myocarditis.

Introduction

Myocarditis is a significant clinical entity in young infants and adolescents (Bejiqi et al., 2019; Tschope et al., 2020). Although the disease is spontaneously resolved in most affected individuals, ∼20% of those affected develop chronic myocarditis that can lead to dilated cardiomyopathy (DCM) (Tschope et al., 2019). More recently, the term myocarditis has been designated as inflammatory cardiomyopathy to describe the occurrence of myocarditis in association with cardiac dysfunction (Tschope et al., 2020). Furthermore, it is not uncommon to detect low-grade inflammation in the hearts of healthy individuals, as has been suggested by a study involving accidental deaths in which heart infiltrates were detected in ∼1–9% of autopsies (Feldman and McNamara, 2000). Consistent with this finding, the presence of heart infiltrates in the sudden deaths of young athletes has raised a question as to the underlying mechanisms (Harris et al., 2021). Therapeutically, due to a lack of effective treatment options, ∼50% of DCM patients undergo heart transplantation, and children with acute myocarditis only have a ∼60% likelihood of transplantation-free survival (Towbin et al., 2006). This is complicated by the finding that myocarditis can result from multiple triggers, whose disease-inducing abilities are complex in nature.

Viruses are the major causative agents of myocarditis (Lasrado and Reddy, 2020; Lasrado et al., 2020; Tschope et al., 2020). Because it is difficult to study the pathogenic mechanisms of viral myocarditis in humans, animal models are commonly employed. However, infection models for all viruses are not available or feasible for routine experimentation—except for enteroviruses that include group B Coxsackieviruses. Thus, Coxsackievirus B3 (CVB3)-induced myocarditis is commonly employed to investigate the pathogenic mechanisms of viral myocarditis.

Of various rodent species, mice are highly susceptible to CVB3 infection (Fairweather and Rose, 2007), and their major histocompatibility complex (MHC) haplotypes influence the disease outcome. Although H-2 (IAb)-bearing C57Bl/6 mice develop acute infection and are resistant to the development of chronic myocarditis, A/J (IAk) and BALB/c (IAd) mice develop chronic disease, making them suitable to study the pathogenic mechanisms that may involve both viral and host factors (Blyszczuk, 2019; Fairweather and Rose, 2007; Lasrado and Reddy, 2020). We routinely use A/J mice that are highly susceptible to a human isolate of the Nancy strain of CVB3 (Fairweather and Rose, 2007; Gangaplara et al., 2012). Upon infection, animals develop myocarditis in two phases in continuum: acute or viral phase lasting approximately 14–18 days, followed by a chronic or nonviral phase in which cardiac dysfunctional features are manifested (Lasrado and Reddy, 2020). During the acute phase, CVB3 can directly lyse cardiomyocytes, as it is a cardiotropic and lytic virus, and activation of innate immune cells contribute to tissue damage through the production of various inflammatory cytokines such as interleukin (IL)-1, IL-6, tumor necrosis factor (TNF)-β, soluble IL-1 receptor-like 1 also called, sST2 and interferons (IFNs) (Fairweather et al., 2003; Nishimura et al., 2018). As the infection is established, T helper (Th) responses appear to play a critical role in developing chronic myocarditis/DCM, but their underlying mechanisms are complex. For example, Th1 responses mediate protective functions in acute myocarditis and inhibit replication of virus (Fairweather et al., 2004, 2005). On the other hand, the Th2 cells could attenuate acute myocarditis (Frisancho-Kiss et al., 2007; Liu and Huber, 2011) but can also promote cardiac remodeling events, leading to chronic myocarditis/DCM associated with cardiac fibrosis (Abston et al., 2012). Conversely, Th17 cells may contribute to the development of acute myocarditis and cardiac remodeling/DCM (Fairweather et al., 2012). These observations suggest that a combination of multiple Th subsets may be involved in the chronic progression of viral myocarditis. Identifying their molecular pathways is thus critical to understanding viral disease pathogenesis.

Translationally, the chronic nature of CVB3 is particularly relevant to humans because signatures of enteroviruses have been identified in DCM patients, as indicated by the detection of virus-reactive antibodies and the viral genomic material (Fujinami et al., 2006; Lasrado et al., 2020; Tschope et al., 2020). These findings raise the question of whether residual virus, if any, or reactivation of viral nucleic acid, if at all possible, may potentially contribute to DCM pathogenesis. But conclusive evidence is lacking to support these notions. Another possibility is autoimmunity, as autoantibodies have been detected in both DCM patients and CVB3 infection models in mice (Caforio et al., 2008; Fujinami et al., 2006; Lasrado et al., 2020; Neumann et al., 1991; Tschope et al., 2020). In our studies, by creating MHC class II dextramers for various cardiac antigens, we have demonstrated that CVB3 infection leads to the generation of pathogenic autoreactive T cells with multiple antigen-specificities localized in both lymphoid and nonlymphoid organs with the potential for them to be recirculated back into the heart under inflammatory conditions (Basavalingappa et al., 2020; Gangaplara et al., 2012).

Nonetheless, if autoimmunity is a key underlying mechanism for DCM pathogenesis, affected patients should be responsive to immune therapies, but mixed successes have been achieved in clinical trials (Tschope et al., 2019). Furthermore, if self-reactive, immune cells have free access to the heart muscle to be able to cause tissue destruction. From the standpoint of immune defense mechanisms, the arrival of immune cells to damaged cardiac tissue is expected to be beneficial to the host, but their detrimental effects cannot be discounted if their functions are dysregulated. In addition, numerous cardiac resident cells—specifically, cardiomyocytes, fibroblasts, and smooth muscle cells (SMCs)—and tissue-resident immune cells, such as macrophages and dendritic cells (DCs), may be severely affected by inflammatory responses, leading to alterations in cardiac functions. Collectively, many cell types may potentially participate in the cardiac dysfunction that culminates in cardiac remodeling events. To dissect this complexity, we used single-cell RNA sequencing (scRNAseq) to define the cardiac cellulome and its transcriptome profiles during the postinfectious phase of myocarditis in A/J mice. The data revealed detection of mainly myeloid cells, T cells, and fibroblasts in the heart infiltrates from myocarditic mice, which may have a role in developing chronic myocarditis and DCM.

Results

Using the mouse model of viral myocarditis induced with CVB3, we analyzed cellular infiltrates in hearts to identify genes, transcription factors (TFs), and signaling pathways that contribute to disease progression in the viral pathogenesis.

Myeloid cells, T cells, and fibroblasts are the major enriched cell types in the hearts of myocarditic mice

To elucidate cellular compositions and diversity in their transcriptome profiles, we performed scRNAseq using heart cells from both mice infected with CVB3 and healthy mice (Figure 1A). Single-cell suspensions obtained from myocarditic and healthy mice were stained with annexin-V and propidium iodide (PI) for sorting viable cells (annexin-V-PI-) by flow cytometry (Figures 1A and S1). After confirming the viability (100%), 16,000 cells from the healthy and myocarditic groups were loaded into the 10x genomics chromium 3′expression system, and their libraries were sequenced for downstream analysis (Figure 1A).

Figure 1

Phenotypic characterization of heart infiltrates in CVB3-infected mice

(A) Schematic representation of the experimental approach. Hearts harvested from A/J mice infected with or without CVB3 on day 21 post-infection were enzymatically digested to obtain single-cell suspensions. Single-cell libraries were prepared and sequenced after sorting the viable cells by flow cytometry as described in STAR Methods. Raw data were subjected for downstream analysis to characterize cellular distributions.

(B) Mapping of cell clusters. Uniform Manifold Approximation and Projection (UMAP) visualization of cells from healthy (9,734 cells) and myocarditic mice (13,251 cells) using Seurat identified 26 different clusters after unsupervised clustering.

(C) Distribution of cells between treatment groups. UMAP projections showing the relative distribution of cell clusters in healthy and myocarditic mice. (D) Prediction of cell types. Cell types were predicted using singleR R package (Borcherding et al., 2021); normalized correlation values for predicted immune cell phenotypes are shown. Cluster of columns based on Ward.D2 distance between normalized correlation values across all pure immune cell populations in the Immgen database (Heng et al., 2008).

(E) Identification of cell types using select lineage markers. Using the canonical markers, major cell types in both healthy and myocarditic mice were assigned.

(F) Annotating cell types. Twelve major cell types were identified and annotated based on the expression pattern of canonical cell markers.

(G) Relative distribution of cell types. The relative distribution of cell types scaled by the total number of cells per condition is shown. Red indicates myocarditis, and blue indicates control. Significance test on two-proportion Z-tests with p values correct for multiple comparisons using the Benjamini-Hochberg method; ∗<0.05, ∗∗<0.01, ∗∗∗<0.001, and ∗∗∗∗<0.0001.

A combined total of 22,985 cells from healthy (n = 9,734) and myocarditic mice (n = 13,251) were analyzed using the Seurat R package (Butler et al., 2018), and an unbiased clustering yielded 26 cell clusters (Figure 1B) along a single uniform manifold approximation and projection (UMAP). In general, cells from both control and myocarditic hearts were present in most clusters (Figure 1C). We next annotated each cluster by using two approaches: (1) correlation of mouse gene signatures in the Immune Genome Project (ImmGen) database (Heng et al., 2008) and (2) expression of canonical markers, with average gene expression in the clusters. We used the SingleR R package (Aran et al., 2019) to assign cell types to each cluster based on correlations between gene expression in the cluster and gene expression in purified mouse cell populations in the ImmGen database. Figure 1D shows the predicted proportions of cells identified as natural killer (NK) T cells, T cells, Tγδ cells, innate lymphoid cells (ILCs), NK cells, fibroblasts, stromal cells, neutrophils, macrophages, and monocytes that were present in multiple clusters. In contrast, B cells, DCs, endothelial cells (ECs), and basophils were restricted to one cluster each. Next, using canonical markers for indicated cell types (Table S1), as shown with two examples for each cell type (Figure 1E), we noted that T cells, fibroblasts, neutrophils, and erythroid cells were present at relatively higher proportions than B cells, NK cells, myeloid cells, ECs, and SMCs. Based on the above two approaches, cluster annotations were then made (Figure 1F). By comparing only the relative proportions of each cell type in both groups, we found that T cells (CD4 and CD8), fibroblasts, and myeloid cells were significantly enriched in myocarditis, whereas neutrophils, NK cells, ILCs, and B cells were reduced in myocarditic mice compared to healthy mice (Figure 1G, Table S2). Although ECs and SMCs were also elevated in myocarditic mice (Figure 1G), their number was relatively low (Table S2). Nonetheless, detection of the predominant populations of cells described above in myocarditic mice, notably myeloid cells, T cells, and fibroblasts raised questions as to their significance in CVB3 pathogenesis.

Transcriptome analysis of myeloid cells led to the detection of predominantly monocytes and macrophages with pathways involved in immune metabolism and inflammation

To understand the contribution of myeloid cells in post-infectious myocarditis, we re-clustered the myeloid cells separately. As shown in Figure 1B, myeloid lineage cells were scattered in five clusters—12, 14, 15, 19, and 25. Further unbiased subclustering led us to identify six distinct subpopulations (Figure 2A), with cellular proportions varying between healthy and myocarditic mice (Figure 2B, left panel). Cells in clusters 0 and 1 were significantly elevated, followed by clusters 4 and 2 in myocarditic mice (Figure 2B, right panel). By using canonical myeloid cell markers, cytokines, chemokines, and other molecules as shown (Figures 2C and S2A), we identified monocytes, macrophages, cDC1, moDC, and a discrete subset of cells that expressed predominantly mitochondrial genes (Mt-high) (Figure 2D, left panel). Although DCs contained two subpopulations, no differential gene expression (DGE) was noted in their transcriptomes. Furthermore, we used Slingshot (Street et al., 2018) to infer the cell lineage and pseudo-time trajectory, which indicated the branching of monocytes into other cell types (Figure 2D, left panel). The top immune-related genes used to construct this model are shown (Figure 2D, right panel), with Msrb1, an anti-inflammatory selenoprotein (Rose and Hoffmann, 2015), making the largest contribution to this developmental trajectory. Evaluation of cell-cycle phases revealed monocytes from healthy mice in the G2M and G1 growth phases, whereas monocytes, macrophages, and moDCs from myocarditic mice were in the S and G1 phases, suggesting that the growth stages of the latter cells might represent their activation status (Figure 2E).

Figure 2

Distribution and characterization of myeloid cells in heart infiltrates

UMAP of myeloid cells identified in heart infiltrates (A) and their distributions are shown in healthy controls and myocarditic mice (B); relative proportions of cells are indicated by cluster in the bar plot, with red indicating myocarditis and blue indicating control. By using select markers for various cell types (C), individual subsets were then identified (D); the cell trajectory is shown using the slingshot method with cluster 0 as the origin, and the ranked bar graph indicates the top immune-related genes across conditions. The bar graph represents cell-cycle phases by cluster between healthy and myocarditic mice (E). Percentage difference (Δ Percent of Cells) and log-fold change based on the Wilcoxon rank-sum test results for differential gene expression comparing myocarditis versus healthy controls in clusters 0, 3, 1, 4, respectively, are shown in (F). Genes highlighted in red or blue have adjusted p values < 0.05. qPCR analysis of selected transcripts in heart tissues obtained from myocarditis and control mice is shown (G). Mean ± SEM fold-change gene expression values representing three samples per group, each containing n = 3 to 5 mice, are shown. Mann-Whitney test was used to determine significance between groups. ∗p ≤ 0.05, ∗∗p ≤ 0.01, ∗∗∗p ≤ 0.001.

To understand the functional role of myeloid cells in viral myocarditis, we compared gene expression profiles between groups, leading us to note 354 upregulated and 389 downregulated genes overall (Table S3). Among various subtypes of myeloid cells (Figure 2D), differences in gene expression patterns were noted in monocytes representing clusters 0, 3, and 4 and in macrophages in cluster 1 (Figure 2F). We then sought to understand the significance of differentially expressed genes in each cluster by focusing on genes whose log-fold change (logFC) expressions were ≥0.25 with adjusted p value < 0.05, as well as a greater than 10% difference in the percentage of cells expressing the gene (Δ Percent of Cells) in myocarditic mice as compared with controls. Based on these criteria, in monocyte cluster 0 of myocarditis, we saw the upregulation of M2 macrophage marker genes Ccl24, Arg1, Gatm, and Chil3, which play roles in antiinflammatory functions and fibrosis/tissue repair (Murray et al., 2014; Pope et al., 2005). Similarly, expression of Tgfbi and S100a4, which mediate tissue repair and survival of cardiomyocytes (Schneider et al., 2007; Schwanekamp et al., 2017), may signify inhibitory functions of this cluster. Ccl2 and Ccl9, markers of the M1 phenotype implicated in myocardial infarction (Dusi et al., 2016; Morimoto and Takahashi, 2007), were also upregulated, thus suggesting a mixed phenotype for monocyte cluster 0 in myocarditis. Likewise, monocytes of cluster 3 in myocarditis had upregulation of Ccl24 and Ccl2, in addition to Lgals1 and Wfdc17 that counter inflammation (Karlstetter et al., 2010; Sundblad et al., 2017), suggesting that blood-derived monocytes arriving at inflamed hearts may take part in the reparative process. Similar analysis of the monocytes in cluster 4 revealed upregulation of mainly Fcgr2b, an inhibitory Fc receptor; Ms4a6d, a suppressor of IL-1b via NLRP3 activation; and Gatm, an activator of arginine metabolism, indicating their anti-inflammatory roles. Finally, transcripts in the macrophages of cluster 1 of myocarditis demonstrated upregulation of mainly Ccl24, Chil3 (M2 markers), and those involved in the modulation of macrophage functions (Fcgr2b) and inflammation and/or cardiac remodeling (Mt1, Wfdc17, and Socs3) (Carow and Rottenberg, 2014; Duerr et al., 2016; Karlstetter et al., 2010) (Figure 2F), suggesting a role for these molecules in viral myocarditis. We validated a few of the upregulated genes, namely Ccl24, Chil3, and Gatm that are M2 markers, and Mt1 that has a role in immune regulation, by qPCR, and their roles remain to be determined in CVB3 pathogenesis. As shown (Figure 2G), heart tissues obtained from myocarditic mice had significantly increased expression of all four genes compared with controls. It is possible that lack of expression of these molecules may lead to enhanced disease severity in CVB3 infection.

We next performed gene ontology (GO) analysis using the upregulated genes of myeloid cell populations in myocarditic mice. The analysis revealed prominence of metabolic pathways (oxidative phosphorylation, ATP, and TCA cycle), in addition to inflammatory responses, leukocyte migration and activation, hypoxia, and antigen-presentation functions (Figure S2B). Similar analysis corresponding to three clusters in myocarditis (0, 1, and 3) showed inflammatory and wound healing pathways to be highly prominent in the monocytes of cluster 0, whereas pathways related to cardiac muscle contraction, neutrophil degranulation, and hypoxia were prominent in cluster 3 (Figure S2C). In contrast, macrophages of cluster 1 in myocarditis were related to mainly signaling molecules and negative regulators of inflammation and metabolic pathways (Figure S2C). Overall, scRNAseq analysis revealed M2 cells to be the major myeloid cells, with Ccl24 being the major transcript to be upregulated in 3 out 4 clusters in myocarditis (0, 1, and 3) (Figure 2F), implying that Ccl24 may be critical for M2 functions in myocarditic mice. However, because Ccl24 can regulate inflammation and fibrosis in addition to activation and chemotaxis of immune cells and fibroblasts, Ccl24 may have a multitude of functions in viral pathogenesis. Such functionalities may be relevant to both resident macrophages and those derived from the recently migrated blood monocytes, as both cell clusters had upregulation of Ccl24.

Th17 cells form a dominant fraction of T cells in myocarditic mice

In dissecting the complexity of T cells, we identified nine subclusters (Figure 3A) by utilizing various phenotypic markers and unique transcriptional signatures (Figures 3B and S3A). Cells in each cluster varied, in that myocarditic mice had significantly higher proportions of T cells in clusters 1 (CD4+Th17), 2 (CD8+CTL), and 4 (CD4+Tregs), followed by 7 (CD4+Tcf4+), 6 (CD4+Ki-67+), and 5 (CD4+CD8+Ccr1+), whereas proportions of CD4+naive T cells were higher in healthy mice, as expected (Figure 3C). However, no apparent differences in cell-cycle status were noted between the two groups, as cells in most clusters in both healthy and myocarditic mice were in the growing (G1) phase with minor variations between groups (Figure S3B). Although the predominance of T helper (Th) 17 cells and cytotoxic T lymphocytes (CTLs) may indicate their effector functionalities, detection of Treg cells was not expected, and their presence indicates that Treg cells have a major role in the postinfectious phase of myocarditis. The significance of the presence of Tcf4+CD4 T cells is not clear, but Tcf4 can promote NF-kB activation (Liu et al., 2016). Detection of CD4+Ki-67+ cells indicate their proliferating status, and the presence of CD4+CD8+ cells expressing Ccr1 may suggest a possible role for them in disease mediation.

Figure 3

Analysis of T cell clusters reveals Th17 cells, CTLs, and Tregs to be the dominant fraction in myocarditis

Nine clusters were identified corresponding to healthy and myocarditic mice (A), using select markers for the indicated T cell types (B), and their proportions relative to the total number of cells per condition are shown. Red indicates myocarditis, and blue indicates control. (C). Based on the mRNA expression pattern and log-fold change for the top eight genes in each cluster, a heatmap was created, in which dot size equates to the percentage of cells expressing the gene, whereas color corresponds to the expression level (D). Percentage difference (Δ Percent of Cells) and log-fold change based on the Wilcoxon rank-sum test results for differential gene expression comparing myocarditis versus healthy controls in CD4+ and CD8+ T cells (E) and Th17, Treg, and CTLs (F), respectively, are shown. Genes highlighted in red or blue have adjusted p values < 0.05. Cell trajectory using the slingshot method (Street et al., 2018) with clusters 0 and 3 as the origin is shown; the top two genes with divergent expression at the root (Igfbp4) and at the terminal branches of the trajectory (Nkg7) are indicated in the UMAP (G). The ranked bar graph indicates the top immune-related genes in the order of their abundance in expression. GO analysis with enriched GO terms for various pathways (H) is shown for CD4+ (top panel) and CD8+ T cells (bottom panel). qPCR analysis of selected transcripts in myocarditis and control mice is shown (I). Mean ± SEM fold-change gene expression values representing three samples per group, each containing n = 3 to 5 mice, are shown. Mann-Whitney test was used to determine significance between groups. ∗p ≤ 0.05.

We next identified the top eight differentially expressed genes in each T cell cluster across both healthy and myocarditic mice (Figure 3D). To understand their significance in viral myocarditis, we performed DGE on CD4 and CD8 T cells (Figures 3E and S3C, Table S3), which revealed the detection of upregulated genes that mediate various functions. For example, in myocarditic mice, we noted the upregulation in CD4 T cells of Ccl5, which is a target for NF-kB activation (Huang et al., 2009), whereas Cxcr6 mediates recruitment of Th17 cells and CD8 T cells (Butcher et al., 2016; Sato et al., 2005). Strikingly, detection of Nkg7, which is implicated in cytotoxic functions (Ng et al., 2020), may mean that a proportion of CD4 T cells infiltrating hearts in viral myocarditis may have a cytotoxic function that has not been investigated thus far. Likewise, enhanced expression of Pdcd1 and Id2 may indicate the existence of potential checkpoints for T cell functions, as they are implicated in T cell exhaustion (Jubel et al., 2020) and plasticity of Treg cells (Hwang et al., 2018), respectively. Interestingly, DGE in CD8 T cells of myocarditic mice revealed increased Ccl5, S100a6, S100a4, Pdcd1, Cxcr6, and Bcl2a1b similar to CD4 T cells, implying their common functionalities in both CD4 and CD8 T cell subsets (Figure 3E). Although expression of Gzmb validates the identity of CD8 T cells, Ly6a, a regulator of memory T cell development, is also expressed in CD4–CD8–Treg cells (Zhang et al., 2002) and may have a role in viral myocarditis.

Because T cell infiltrates in myocarditic mice were predominantly Th17 cells, Treg cells, and CTLs (Figure 3C), we sought to analyze gene expression profiles in these subsets, expecting to identify genes of interest that may be involved in the viral pathogenesis. Th17 cells in myocarditic hearts when compared with those in controls (Figure 3F, top left panel) showed elevated expression of genes that were also upregulated in the complete CD4 T cell DGE analysis (Figure 3E). Although expression of Cxcr6 was expected in Th17 cells, upregulation of Nkg7, which mediates cytotoxic function (Ng et al., 2020), was not expected (Figure 3F). In addition, we also noted the increased expression of S100a4, a calcium-binding gene that has been previously shown to have increased expression during viral myocarditis with tissue repair functions in the myocardium (Schneider et al., 2007; Wang et al., 2019). Likewise, we also noted upregulation of other calcium-binding genes in the Th17 cells, such as S100a6 and S100a11 whose role in CVB3 myocarditis has not been investigated thus far. The Treg cells in cluster 4 of myocarditic mice (Figure 3F, top right panel) had upregulation of their known markers Lag3, Ly6a, Tnfrsf9, Izumo1r, Id2, Cst7, and fructose-bisphosphate aldolase, along with increased expression of Nkg7 and Pfn1, which mediate cytotoxicity (Ng et al., 2020; Young et al., 1986), indicating that the cytotoxic Treg cells may be critical to maintaining homeostasis in the local inflamed heart milieu. Similarly, CTLs of cluster 2 in myocarditic mice (Figure 3F, bottom panel) indicated expected upregulation of Gzmb, CD8 co-receptors, and Cd3 complex proteins; elevations of Icos and Pdcd1 may indicate activation status of CTLs (Dong et al., 2001; Jubel et al., 2020). We next built cell trajectories using Slingshot (Figure 3G). All T cell clusters had a common root-point from CD4+naive cells (cluster 0) together with CD8+ CM cells (cluster 3), branching into CD4+Th17 (cluster 1), which gave rise to terminal branches of other T cell populations (Figure 3G). We noted two genes with divergent expressions at the root (Igfbp4) and at the terminal branches (Nkg7) (Figure 3G). Igfbp4 can mediate both proliferative and inhibitory functions, whereas Nkg7 is implicated in the cytotoxic functions (Ng et al., 2020). The finding that upregulated expression of Nkg7 was evident in the whole CD4+ T cell subset (Figure 3E), Th17 cells, and Treg cells in myocarditis (Figure 3F, top panel) suggested the possibility that Nkg7 may have a role in the effector functions of these cell types in viral myocarditis. Furthermore, gene set enrichment analysis (GSEA) indicated that the genes in the whole CD4 T cell subset were primarily involved in immune metabolism, T cell activation, Th1, Th2 and Th17 differentiation, apoptosis, and cytokine production, whereas T cell activation, cytotoxicity, apoptosis, and NF-κB signaling, and inflammatory-response-related pathways were apparent in CD8 T cells of myocarditic mice (Figure 3H). These observations were also recapitulated individually in Th17 cells and CTLs of myocarditic mice (Figure S3D). Together, the data suggest prominent roles for Th17 cells, Tregs, and CTLs in the immune pathogenesis of viral myocarditis, with a possibility that the effector functions may be mediated mainly by cytotoxic functions regardless of T cell subsets. Whether such functionalities are beneficial for eliminating virus-infected cells or detrimental via tissue destruction should they be self-reactive is worth investigating in future studies. To that end, qPCR analysis of four upregulated genes, namely Nkg7 that is implicated in cytotoxicity and three genes involved in calcium binding (S100a4, S100a6, and S100a11) whose functions may also appear to be modulated in exhausted CD8 T cells (Wherry et al., 2007), showed expression of all the four transcripts to be significantly increased in hearts from myocarditic mice compared with controls (Figure 3I). Currently, data are lacking regarding the functionalities of Nkg7 and calcium-binding genes (S100a6 and S100a11) in CVB3 pathogenesis, and utilization of animals with a conditional knockout of these genes in T cells may be helpful to ascertain their role in viral myocarditis.

Fibroblasts in myocarditic mice contained various subtypes with unique transcriptome signatures that can influence the functionalities of immune cells

We analyzed the fibroblast population based on the known markers (Figure 1E) and re-clustered it into 12 distinct populations as indicated in the UMAP (Figure 4A). We analyzed the relative proportion of cells in each cluster and noted that clusters 1, 5, 6, and especially 8 were present in greater numbers in myocarditic mice than in healthy mice (Figure 4B). In addition, we noted that fibroblasts in myocarditic mice tended to be in the S phase compared with healthy mice, and such a trend was more evident in cluster 8 than others (Figure S4A). We next analyzed gene expression profiles of the top eight differentially expressed genes in all 12 fibroblast populations, revealing a few noteworthy findings. Expression of enriched genes, although unique in all clusters, showed overlapping patterns in clusters 5, 6, and 8 (Figure 4C). We noted three genes, Postn (periostin), Ltbp2 (latent TGF-β binding protein 2), and Thbs4, of interest in cluster 5. Although Postn, an ECM protein, is critical for tissue development and regeneration and plays a role in wound healing and ventricular remodeling in myocardial events (Shimazaki et al., 2008), Ltbp2 and Thbs4 have roles in myocardial fibrosis in DCM (Pang et al., 2020) or the fibrotic process (Frolova et al., 2012). Similarly, three transcripts, Apoe, Tgfbi, and Mfap4, were found interesting in cluster 6. They mediate regulation of T cell and macrophage functions, including inflammation and oxidative stress (Bonacina et al., 2018), cardiac fibrosis in concert with Postn (Schwanekamp et al., 2017), and ventricular remodeling and cardiac function, respectively (Wang et al., 2020). As to cluster 8, we noted that a few genes were uniquely expressed with varied functions. These include Apoe, Tgfbi; Wif1, which modulates cardiomyocyte differentiation (Buermans et al., 2010) necessary for cardiac remodeling in myocardial infarction and also promotes DCM (Lu et al., 2013); Clu, myocardial injury marker (Swertfeger et al., 1996); and Npy, which has a role in cardiac remodeling and heart failure (Widiapradja et al., 2017). These data indicated that clusters 5, 6, and 8 were mainly involved in the cardiac remodeling process in myocarditis. Although there also were more fibroblasts in clusters 1 and 3 in myocarditic mice, their transcriptional profiles were unique to each (Figure 4C). The genes expressed in cluster 1 include Mt1 and Mt2, which can modulate inflammation and support cardiac remodeling as demonstrated in ischemic cardiomyopathy (Duerr et al., 2016); Thbs1, a promoter for Tgfbi activation in fibrosis (Murphy-Ullrich and Suto, 2018); Cxcl1, which facilitates recruitment of neutrophils and nonhematopoietic cells to the site of injury and regulates immune and inflammatory responses (Dusi et al., 2016); and IL-6, an inducer of MMP1 that can mediate tissue remodeling (Fontes et al., 2015). Regarding cluster 3, expressions of C3 (activator of complement system) and Apod (putative marker of non-proliferating and senescent fibroblasts) (Provost et al., 1991) were elevated. These observations suggest that upregulated genes in the corresponding clusters have a role mainly in cardiac remodeling events and immune activation.

Figure 4

Analysis of fibroblasts reveals new cardiac fibroblasts involved in cardiac remodeling during viral myocarditis

(A) The UMAP of isolated fibroblasts across groups; their relative distributions in healthy and myocarditic mice are shown in UMAPs and cluster-wise in bar plots, with red indicating myocarditis and blue indicating control (B). Panel (C) represents the heatmap of the top eight markers shown with log-fold change, where dot size and color represent the percentage of cells expressing the gene and expression levels, respectively. Panel (D) indicates differential expression of the top 10 genes in myocarditic mice compared with healthy mice with respect to clusters 1, 5, 6, and 8. Percentage difference (Δ Percent of Cells) and log-fold change are based on the Wilcoxon rank-sum test results. The panel (E, top panel) shows localization of Wif1 expression in fibroblast cluster 8, and its correlation with other markers Vim, Cyb5a, and Fxyd6 (E, bottom panel). Slingshot analysis along the UMAP (F) shows the root and terminal branches of the cell fates of fibroblast populations. Branches terminating in Wif1+ fib and Cilp+ fib cell clusters promote cardiac remodeling. GO analysis with enriched GO terms for various pathways (G) is shown for the whole fibroblast population. qPCR analysis of selected transcripts in myocarditis and controls with mean ± SEM fold-change gene expression values representing three samples per group, each containing n = 3 to 5 mice, are shown (H). Mann-Whitney test was used to determine significance between groups. ∗p ≤ 0.05.

To further determine the specific roles of potential transcripts that might contribute to CVB3 pathogenesis, we performed DGE in fibroblast clusters 1, 5, 6, and 8 of myocarditis and control (Figure 4D). These analyses revealed increased expression of complement proteins C3 and C4b; Serpina3n, an accelerator of wound healing (Hsu et al., 2014); Tmem176b, which controls DC maturation and also has a role in fibrosis (Etikala et al., 2017); Tmem176a, a negative regulator of DCs (Picotto et al., 2020); and the anti-viral protein Ifitm (Diamond and Farzan, 2013) in all clusters of myocarditic mice (Figure 4D). Uniquely, however, Sbno2, a regulator of pro-inflammatory cascade (El Kasmi et al., 2007), and Ly6a, also called Sca-1, which has a known function in heart failure (Chen et al., 2018), was increased in cluster 1 of myocarditic mice, as opposed to Gsta3, an inhibitor of TGF-β-induced epithelial and mesenchymal transition and fibronectin expression (Xiao et al., 2016), and Socs3 (Carow and Rottenberg, 2014) in cluster 5; Apod (Provost et al., 1991) in cluster 6; and Apoe in cluster 8 of myocarditic mice (Figure 4D).

Because there were more fibroblasts in cluster 8, which also expressed injury marker Wif1, we sought to correlate its expression with other genes (Figure 4E, top panel). This analysis revealed a positive correlation with Vim. However, unexpectedly, two other genes (Cyb5a and Fxyd6) also showed positive correlation with Wif1 expression (Figure 3E, bottom panel); whether their expression can also be used as injury markers remains to be investigated. We next explored the origins of cardiac remodeling-associated fibroblasts by analyzing transcriptional activation using Slingshot. The analysis suggested that Cilp+ cluster 5 had its origins in cluster 7, whereas Wif1+ cluster 8 arose from cluster 6, which in turn had its origins in cluster 7 (Figure 4F, left panel), with Tmsb4x, having roles in the repair of human heart muscle, being the essential transcriptional activator for this cell trajectory (Figure 4F, right panel). By performing GO analysis on all populations, we noted the upregulation of functions associated with inflammatory responses and other immune signaling pathways (Figure 4G), including cytotoxicity, ossification, ECM protein synthesis, and immune/inflammatory regulatory networks, among others (Figures 4G and S4C), suggesting that fibroblasts can perform diverse functions in CVB3 viral myocarditis. Finally, we validated the differential expression of four genes of interest namely, Wif1, Cy5ba, Fxyd6, and Thbs4 by qPCR, and expectedly, their expression was elevated in myocarditic mice (Figure 4H). Of these, although Wif1 and Thbs4 have been implicated in the development of fibrosis in other models (Farbehi et al., 2019; McLellan et al., 2020), their functional role remains to be investigated in CVB3 infection. Our data may support a notion for their fibrogenic properties. Although elevated expressions of Cy5ba and Fxyd6 could suggest their importance as potential biomarkers of fibrosis/myocarditis, determination of their functional roles in CVB3 pathogenesis requires additional studies.

NK cells, ILC2, and ILC3 cells formed a major component of ILCs, but their numbers were low in myocarditic mice

We analyzed ILCs pooled from healthy and myocarditic mice, and the use of canonical markers allowed us to dissect ILCs into five distinct populations (Figure 5A). As shown in Figure 5B, three markers of NK cells (Nkg7, Klrd1, and Gzma) were consistently expressed in clusters 0, 1, and 3. Conversely, ILC markers (Gata3, Ltb4r1, Csf2, and Il17a) were expressed in clusters 2 and 4. By evaluating the relative proportion of cells present in each cluster, we observed that the Gzma+ NK cells, Tcf7+ NK cells, and Gata3+ ILCs were significantly reduced in myocarditic mice as compared with healthy controls (Figure 5C). Next, we evaluated the DGE in NK cell clusters and identified the top 15 differentially expressed genes between groups (Table S3). Figures 5E and S5A show upregulation of genes involved in NK cell development (Ltb and Thy1) (Kupz et al., 2013) and NK cell function (Gzmc). Genes implicated in cell migration, invasion, and fibrosis (S100a4 and S100a6) (Schneider et al., 2007) and anti-viral function (Ifi27L2a) (Diamond and Farzan, 2013) were also upregulated in myocarditis. Another gene, Tmem176b, may be an interesting candidate with a role in cardiac fibrosis, as upregulation of this gene has been noted in pulmonary fibrosis (Etikala et al., 2017). Finally, the GO analysis found that all NK cell clusters had a role in the activation and regulation of T cells and neutrophils and anti-viral responses in myocarditis (Figure 5F), with similar functions being noted in Gzma+ NK cells (cluster 0) in addition to the apoptotic signaling pathways (Figure S5C). Together, the data suggest a role for infiltrating NK cells in myocarditic mice in postinfectious myocarditis, but with no or minimal contribution from ILCs. Because the GO analysis revealed a role for NK cells to mediate anti-viral responses, NK cells might act in concert with T cells to eliminate the residual virus during the postviremic phase of CVB3 infection.

Figure 5

Analysis of ILCs in heart infiltrates of myocarditic mice

(A) The UMAP of ILCs across groups; the canonical markers used to identify the five clusters of ILCs are shown in (B). After their distributions corresponding to healthy and myocarditic mice were identified (C, top panel), the proportion of cells in each cluster relative to the total number of cells per condition was determined, with red indicating myocarditis and blue indicating control (C, bottom panel). The panel (D) indicates cell-cycle assignments across all ILC and NK cell clusters in healthy controls and myocarditic mice, Χ2 test p value < 2 × 10−16; in panel (E), the percentage of differentially expressed transcripts (Δ Percent of Cells and log-fold change) in NK cells in myocarditic mice is shown. Genes highlighted in red or blue have adjusted p values < 0.05. Panel (F) indicates enriched GO terms concerning pathways upregulated in NK cells of myocarditic mice.

Neutrophils infiltrating the myocarditic hearts predominantly modulate inflammatory responses

By utilizing the known markers, we subclustered the neutrophils obtained from myocarditic and healthy mice into six populations (Figures 6A andS6A). All clusters had similar proportions of neutrophils in both groups, barring Ly6g+ cluster 2 and Ccl5+ cluster 6. Evaluation of the top eight differentially expressed genes in each cluster revealed Camp, an antimicrobial peptide, and Ly6g, a marker of myeloid cells/granulocytes, as highly expressed genes in Ly6g+ cluster 2, whereas Ccl5, a chemokine expressed by nonneutrophils (Dusi et al., 2016), was identified in cluster 6. Similarly, increased expressions of Lyz2, a chemotactic for neutrophils, and Ccl6, a chemokine expressed by neutrophil and macrophage lineages (Dusi et al., 2016), were seen in Ccl6+ cluster 0. In addition, Fgl2, a regulator of immune and inflammatory responses (Zheng et al., 2018), and Itgax or CD11c, a marker of DCs that triggers a respiratory burst in neutrophils, were found in Gm2a+ cluster 1 (Figure 6C). Interestingly, cluster 3 neutrophils contained genes involved in the regulation of inflammation via inhibition of NF-kB, including Nfkbia, Nfkbie, Tnfaip3, and Tnf. Likewise, Ifit1+ cluster 4 had genes related exclusively to anti-viral activity (Ifi1, Ifit3b, and Oasl1) (Diamond and Farzan, 2013), whereas those in Mt-High cluster 5 were all mitochondrial genes (Mt-nd1/2/3 and Mt-atp6) that have been implicated in mitochondrial cardiomyopathy (El-Hattab and Scaglia, 2016). The data thus suggest heterogeneity in neutrophil populations that mediate varied functions in viral myocarditis.

Figure 6

Neutrophils mainly with pro-inflammatory functions were detected in myocarditic mice

(A) UMAP of neutrophils representing seven clusters across groups.

(B) UMAPs and bar plots demonstrating the relative distribution of neutrophils in healthy and myocarditic mice by cluster, with red indicating myocarditis and blue indicating control.

(C) Heatmap of the top eight markers by log-fold change. Dot size equates to the percent of cells in the cluster expressing the gene, and color corresponds to the expression level.

(D) Percentage differences (Δ Percent of Cells) versus log-fold change of the differentially expressed genes in myocarditic mice relative to the healthy group are indicated for clusters 0,1, 2, and 3.

(E) Enriched GO terms with respect to pathways upregulated in neutrophils of myocarditic mice.

Through DGE analysis in the myocarditic mice, we also noted a pattern that had similarities in several clusters. For example, prominently upregulated genes in Ccl6cluster 0 of myocarditis included Wfdc17, negative regulator of inflammation; anti-viral proteins Ifitm1, Ifitm2, and Lrg1, a neutrophil granule protein and modulator of myelopoiesis that also has a role in wound healing; and Tspo, a mitochondrial membrane protein overexpressed in inflammatory processes (Figure 6D). Similar trends were noted in Ly6g+cluster 2 in myocarditis, except that Il1b was upregulated, which promotes cardiac fibrosis and remodeling in myocardial infarction (Frangogiannis, 2015) (Figure 6D). Likewise, Clec4, an interacting partner of FcRIγ, was increased in cluster 1 of myocarditic mice (Figure 6D). We next performed GO analysis, leading us to note that the pathways involved in the regulation of various immune and inflammatory responses (Il17, NF-kB, TNF-signaling), cell death, cytokine signaling (Il1b, Il4 signaling), as well as modulation of the viral life cycle, were upregulated in myocarditis (Figure 6E). Several of these pathways also overlapped across individual neutrophil clusters, especially in cluster 1 of myocarditic mice (Figure S6B); Cluster 0 predominantly showed pathways of interferon functions, and cluster 2 showed neutrophil chemotaxis and inflammatory response. In sum, the data showed that neutrophils participate in the pathogenic inflammatory and cardiac fibrosis process, with Il1b being the major driver of this process. This observation is important to note because myocarditis is increasingly being considered as Il1b-mediated disease process, as neutralization of Il1b renders the mice more resistant to the development of CVB3 myocarditis (Wang et al., 2014).

Other cell types detected in myocarditic mice included B cells and, to a lesser extent, ECs, basophils, and SMCs

As presented in Table S2, scRNAseq analysis revealed the identification of B cells, basophils, ECs, erythroid cells, and SMCs in varied proportions. Because of low number, we did not investigate the transcriptomes of basophils, ECs, and SMCs (Table S2). However, by using known markers of B cells, we noted five distinct cell populations (Figure S7A), and no significant differences were noted in the relative proportions of cells in each cluster in myocarditic hearts (Figure S7B). In addition, by analyzing the differentially expressed transcripts, we found no significant differences between groups (Figure S7C), suggesting that the low proportion of B cells that accumulated in the myocarditic mice appeared not to play a major role in the immune pathogenesis of CVB3 infection. This was a surprising finding because immune complexes for various self-antigens have been demonstrated in the heart tissues (Kaya et al., 2012), but it does not mean that the B cells producing autoantibodies should also co-exist with the infiltrates. Alternatively, B cells producing antibodies to self- or viral antigens may be localized in the draining lymphoid tissues, and their products (antibodies) could infiltrate heart tissues.

Ligand-receptor analysis revealed diverse intercellular communications during myocarditis

To determine the intercellular communication networks between various cell types, we used CellChat (Jin et al., 2021), which utilizes the known structural composition of 2,021 validated ligand-receptor interactions and membrane-bound co-receptors deposited in CellChatDB (Jin et al., 2021). Although communications were evident between all cell types detected in both myocarditic (Figure 7A) and healthy (Figure S8A) mice, the number of ligand-receptor pairs involved in these interactions differed between groups as indicated by dense intercellular communication networks. Notable among these, especially in myocarditis, were the following (ligand-receptor): B cells-CD8 T cells; basophils-CD8 T cells and myeloid cells; ECs-CD8 T cells; CD4 T cells-CD8 T cells; CD8 T cells-neutrophils, myeloid cells, and NK cells; fibroblasts-SMCs and myeloid cells; and ILCs-CD8 T cells (Figures 7B and S9). CellChat also allowed us to evaluate the differential number of interactions and their strength in the myocarditic cardiac cellulome, as compared with healthy cells, based on the upregulated ligand-receptor pairs. We made a few observations: (1) the interaction strength and number of signals sent from fibroblasts to CD4 and CD8 T cells and SMCs were enriched in myocarditis (Figures 7C and S8B). (2) Myeloid cells appeared to have stronger and more interactions with CD4 T cells and SMCs. (3) The number of SMC interactions between most of the cell types and their relative strength were enriched. However, the analyses revealed no indication of autocrine networks for fibroblasts and myeloid cells (Figures 7A, 7C, and Figure S8B). GO enrichment analysis of upregulated interactions between cell types revealed that the ligand-receptor pairs involved in fibroblasts, CD4 and CD8 T cells, were associated with positive regulation of proliferation of epithelial cells and SMCs, adhesion and migration of cells, and inflammatory cytokine responses (IL-1b, IFN-γ, and TNF-α) (Figures 7D and 7E and Table S4). In combination with fibroblasts, ligand-receptor pairs involved in myeloid cells were also associated with cell adhesion and anti-inflammatory and wounding and regeneration processes. However, fibrosis-associated processes were mainly restricted to the receptor-ligand pairs in the fibroblasts, suggesting a prominent role for them in the development of cardiac fibrosis in viral myocarditis.

Figure 7

Intercellular communication between cardiac cell types in myocarditis

(A) Circle plot showing the intercellular communication between major cardiac cell types, using CellChat R workflow. The lines originating from a cell type indicate ligands being broadcast, with these lines connecting to the cell type where the receptors are expressed. The thickness of the line is proportional to the number of unique ligand-receptor interactions, with loops representing autocrine circuits.

(B) A detailed view of ligand and cognate receptor interaction for major cell types.

(C) Heatmap of differential interactions between myocarditic and healthy mice in the cell-cell communication network. The top-colored bar plot indicates the sum of column values (incoming signaling), and the right bar plot indicates the sum of row values (outgoing signaling). Red indicates increased signaling in myocarditis, and blue indicates decreased signaling.

(D) GO terms enriched in a set of genes that encode ligands upregulated in myocarditis. GO terms are ordered by their frequency of significant enrichment in different cardiac cell populations.

(E) GO terms enriched in a set of genes that encode receptors upregulated in myocarditis. GO terms are ordered by their frequency of significant enrichment in different cardiac cell populations. Dot size indicates the number of enriched genes, with colors indicating −log10 (adjP) value.

(F) All significant signaling pathways ranked based on their relative information flow within the inferred networks between healthy and myocarditic mice. The top signaling pathways colored red are more enriched in control mice; those colored black are equally enriched in control and myocarditic groups, and the blue colored pathways are more enriched in myocarditic mice

We next identified 66 signaling pathways associated with the ligand-receptor pairs we had identified and noted these to be enriched in either myocarditic (EGF, CD52, CXCL, CCL, MHC-I, LCK, ncWNT, AGRN, OSM, NOTCH, VTN, PDGF, GRN, and CD45) or healthy mice (MK, THY1, PTN, CD80, CD23, CDH1, NECTIN, CLEC, EPHB, VISFATIN, IL2, and MPZ, among others), whereas a few others were enriched equally in both (Figure 7F). Notably, most of the pathways enriched in myocarditic mice, such as EGF, LCK, ncWNT, OSM, NOTCH, and PDGF, have been implicated in CVB3 pathogenesis or atherosclerosis/cardiac fibrosis. For example, in our data, we observed that the EGF signaling involved in pathogenic plaques and remodeling of blood vessels (Makki et al., 2013) was upregulated in myocarditis, with signals being sent from the ILCs and SMCs to the fibroblasts (Figure S10). LCK signaling of the Src family p56lck, which is essential for CVB3 replication and pathogenesis (Liu et al., 2000), had signaling interactions from NK cells and CD4 T cells to CD8 T cells, with additional strong autocrine signaling on CD8 T cells (Figure S10). For the ncWNT pathway, which plays a role in cardiac remodeling (Bergmann, 2010), signals were transmitted from SMCs to ECs, with autocrine signaling on SMCs (Figure S10). OSM, which is upregulated in DCM patients (Kubin et al., 2011), included signals sent from basophils, myeloid cells, and neutrophils to SMCs (Figure S10). NOTCH signaling, which plays a role in repairing the myocardium (Ferrari and Rizzo, 2014), was secreted from SMCs in an autocrine fashion and to the ECs and neutrophils (Figure S10). The profibrotic signaling in PDGF was secreted by the ECs and SMCs to fibroblasts, in an autocrine network, as well. Finally, CD45 signaling was strongly secreted from almost all cells to the myeloid cells (Figure S10). Overall, the intercellular communication patterns in the cardiac cellulome and the signaling pathways being prominently upregulated in myocarditic mice may contribute to and play a critical role in the pathogenic progression of viral myocarditis. Because interactions of most of the cells were between a soluble ligand and a membrane-bound receptor, determining the significance of localization of cells is difficult. Alternatively, soluble factors may be involved in cell-cell interactions as might occur with cytokines.

TFs enriched in myocarditis modulate cardiac remodeling functions of target genes

In order to identify the TFs that regulate the differential expression of genes in different cell clusters, we employed single-cell regulatory network inference and clustering (SCENIC) (Aibar et al., 2017). This analysis revealed upregulated expression of four TFs in both CD4 and CD8 T cells (Elf1, Ets1, Irf7, and Stat1) (Figure 8A). Functionally, Elf1 is known to regulate anti-viral responses, especially of type I IFNs (Seifert et al., 2019), and Ets1 controls expression of cytokines and chemokines (Russell and Garrett-Sinha, 2010), whereas Stat1 is involved in the signaling cascades of both type I and type II IFN responses. We then analyzed co-expression patterns of these TFs with the upregulated genes in the CD4 and CD8 T cell subsets (Figure 8B, Table S5). Among all TFs, we noted increased interactions mainly between Ets1 and target genes involved in cardiac ischemia, remodeling, and heart failure (Ccr5, Ccl5, Cxcr3, Ccr2, Cxcr6, and S100a4) (Dusi et al., 2016; Schneider et al., 2007). As shown in Figure S11A, the target genes of Elf1 and Ets1 can mediate TCR signaling and T cell activation among other functions (cell adhesion, chemotaxis, and inflammation). Similarly, Stat1-targeted genes facilitate antigen-presentation/cross-priming and immunoproteasome functions in myocarditic mice. Extending these observations to myeloid cells, we noted upregulation of five TFs implicated in various functions (Figure 8A). These include Irf5, Mafb, Maff, Mef2c, and Rara. The putative co-expression patterns showed increased interactions of Irf5 and Mafb with each other, and with target genes involved in fibrosis and remodeling (Tgfbi, Fn1, Ccl24, Ccr2, Ccr5, Ccl2, S100a4) and M2-specific phenotype (Gatm, Arg1, Mrc1) with Irf5, and Mafb-targeted genes mediating various processes, such as antigen presentation, leukocyte activation, response to IFN-γ, and cytokine production or inflammation (Figure S11B). We also extended these analyses to other cell types of interest, neutrophils and NK cells. We noted upregulation of one TF in neutrophils (Bcl3, a regulator of cell proliferation) and one in NK cells (Egr1, a regulator of growth, cell survival, and cell death) in myocarditic mice, but no apparent differences were noted between groups in fibroblasts (Figure S12). Overall, as T cells and myeloid cells form a major component of cellular distribution in myocarditic mice and upregulate multiple TFs that can control various immunological processes involving many of the upregulated target genes in the respective cell types, the data point to a major role for them in the immune pathogenesis of CVB3 infection. To that end, using mice deficient for selected TFs such as Ets1 in T cells and Irf5 in myeloid cells or utilization of antagonists for these TFs may allow us to determine their roles in the mediation of viral myocarditis.

Figure 8

Analysis of myocarditis-specific transcription factors and their target genes

(A) Heatmaps showing the transcription factors being enriched in either myocarditis or control hearts for CD4, CD8, and myeloid cells using SCENIC.

(B) Network plot showing the upregulated target genes for the indicated transcription factors in T cells and myeloid cells. Lines indicate the interaction between TFs and their target genes, and colors indicate the interaction weight range.

Discussion

In this report, we have described the cellular complexity that occurs during the post-infectious phase of viral myocarditis induced with CVB3. Previous reports have delineated the cardiac landscape and intercellular communication in healthy mice (Skelly et al., 2018) and have investigated the fibrosis and cardiac remodeling process in angiotensin II (AngII) mouse models of fibrosis (McLellan et al., 2020). To our knowledge, this is the first report to comprehensively dissect the cardiac cellulome and the intercellular communication networks in viral myocarditis. By using whole heart cells, we captured a majority of immune cells, including fibroblasts, and a fraction of the ECs and SMCs (Figures 1F and 1G), but not cardiomyocytes. This is because the isolation protocol that we adopted does not yield live cardiomyocytes (Skelly et al., 2018). In addition, the pore size limit of ∼40 μm in the 10x genomics microfluidics chip would also exclude the cardiomyocytes. Enrichment of the ECs, SMCs, and cardiomyocytes requires specialized protocols utilizing single nuclei isolation, which were not used in our studies. After ascertaining the identities of various cell types, myeloid cells, T cells, and fibroblasts were significantly increased in myocarditic mice compared with healthy mice (Figure 1G).

In myeloid cells, we detected mainly monocytes and macrophages, and all subclusters of myeloid cells were branched from monocytes, with Msrb1 being the top immune-related gene driving this differentiation of myeloid cells (Figure 2D). Msrb1 is a selenoprotein promoting anti-inflammatory cytokine expression and has been found to be upregulated in mouse models of cardiac stress (Rose and Hoffmann, 2015). The fact that this gene drives the differentiation of myeloid cells could indicate the anti-inflammatory role of these cells in the postinfectious myocarditis phase. In addition, by analyzing the transcriptomes, we noted the upregulation of several genes that have roles in M2 macrophages and anti-inflammatory functions. Importantly, upregulation of Ccl24 was consistently noted in both monocytes and macrophages (Figure 2F). Ccl24, also called eotaxin-2, was shown to promote pathogenic fibrosis in skin, lung, and liver models of fibrosis in mice (Pope et al., 2005; Segal-Salto et al., 2020); it also was produced by F4/80+ macrophages in inflamed hearts, which may facilitate eosinophils in eosinophilic myocarditis (Diny et al., 2016). Although eosinophils are not commonly reported in viral myocarditis, Ccl24 expression may be necessary to recruit monocytes or facilitate their conversion to M2 cells, which may be critical to repair damaged cardiac tissue or participate in cardiac fibrosis. Nonetheless, it is to be noted that Ccl24 generally promotes cell trafficking and regulates inflammatory and fibrotic activities by interacting with Ccr3. Several cell types such as T cells, monocytes, epithelial cells, ECs, and fibroblasts express Ccl24 that can also induce the activation of Ccl3 cells that include immune cells and fibroblasts (Ablin et al., 2010). Thus, enhanced upregulation of CCL24 in CVB3 infection may have a wide range of chemotactic and immunoregulatory functions including fibrosis. Furthermore, we noted an increased expression of Mt1 in the macrophage cluster 1 (Figure 2F) that can modulate inflammation and cardiac remodeling in ischemic cardiomyopathy (Duerr et al., 2016). Our finding of macrophages expressing Mt1 may contribute to cardiac remodeling events during postinfectious phase of viral myocarditis. In addition, the upregulated genes of myeloid cells were found to have a prominent role in metabolic pathways, especially oxidative phosphorylation, ATP, and the TCA cycle (Figures S2B and S2C); the dependency of M2 cells on these pathways has been demonstrated in macrophages (Viola et al., 2019). Taken together, the myeloid cell populations noted in postinfectious myocarditis may primarily be involved in the reparative process in affected animals.

By investigating the T cell landscape, we noted few observations that offer new insights into myocarditis pathogenesis. T cells mainly consisted of Th17 cells, CTLs, and Tregs in myocarditic mice (Figures 3A and 3C). Detection of Th17 cells was not surprising, as IL-17 blockade can ameliorate the severity of CVB3 myocarditis (Xie et al., 2012), but their antigen specificity remains unknown. This is critical because IL-17-deficient mice develop acute myocardial inflammation in the setting of autoimmune myocarditis, but they chronically develop DCM. It may be that a proportion of these T cells are specific to cardiac antigens, as we have demonstrated with MHC class II tetramers and dextramers (Basavalingappa et al., 2020; Gangaplara et al., 2012). However, the appearance of CTLs was expected because of their disease-protective roles in virus infections, and indeed, CD8-deficient animals were previously shown to be highly susceptible to CVB3 infection (Opavsky et al., 1999). Nonetheless, detection of Treg cells was not expected, and their infiltration may be necessary to achieve immune homeostasis by suppressing the ongoing inflammation.

Furthermore, DGE analysis revealed identical transcriptomes between CD4 and CD8 T cells, with Ccl5, Nkg7, and S100a4 being prominent (Figure 3E), the latter two of which have been implicated in cytotoxic and fibrotic functions (Ng et al., 2020; Schneider et al., 2007). Although expression of Nkg7 in CD8 T cells can be related to cytotoxic functionality, its expression in CD4 T cells may suggest the possibility that CD4 T cells may function as cytotoxic CD4 T cells. We noted the upregulation of Nkg7 in Th17 cells, which was a new finding in myocarditis (Figure 3F). Recent reports have shown an upregulated Nkg7 expression in CD4+ T cells in patients treated for visceral leishmaniasis, and Nkg7¯¯ mice infected with Leishmania donovani and Plasmodium berghei had reduced inflammation (Ng et al., 2020). These observations suggest that Nkg7 expression may promote pro-inflammatory functions of Th17 cells. In the same Th17 cluster, we noted the differential expression of two Ca2+ binding S100a proteins (S100a6, S100a11) that can promote apoptosis and generation of ROS (Gonzalez et al., 2020). Whether such mechanisms of these S100a proteins, especially modulation of Ca2+ homeostasis, are relevant to T cells is worthy to investigate in the future studies. Treg cells in the Treg cluster revealed the expression of various genes (Nfg7, Prfn1) implicated in cytotoxicity (Figure 3F). However, this observation may raise the question of whether Treg cells can contribute to the persistence of virus infection by killing virus-reactive CTLs or by suppressing autoreactive T cell responses via cytotoxicity. The former possibility is unlikely because adoptive transfer of Treg cells can suppress the development of CVB myocarditis (Shi et al., 2010). Nonetheless, expression of Ly6a, the marker of double-negative Tregs, was not expected in CD8 T cells (Figure 3F), suggesting whether a proportion of CD8+ Ly6a+ T cells can function as Treg cells in viral myocarditis.

Although fibroblasts form a major part of the cardiac cellulome and are implicated in the acute phase of CVB infection (Figures 1F and 1G, Table S2), their role in the postinfectious myocarditis phase has not yet been studied in detail (Lindner et al., 2014; Yu et al., 2013). Our scRNAseq analysis revealed 12 distinct clusters that express overlapping genes (Figures 4A–4C). Several of these (Wif1, Npy) (Lu et al., 2013) have been previously shown to be involved in wound healing, ventricular/cardiac remodeling, and myocardial fibrosis in DCM pathogenesis (Ltbp2, Thbs4, and Tgfbi) (Frolova et al., 2012; Pang et al., 2020; Schwanekamp et al., 2017). However, a subset of genes (Mt1, Mt2, Cxcl1) (Duerr et al., 2016; Dusi et al., 2016) that have immune-modulatory roles and promote angiogenesis (Vegfd) were also detected in various clusters. By evaluating the top eight genes in all clusters, we were able to identify two categories of subclusters with one gene in each that mediated cardiac remodeling events and anti-viral and immune activation functions (Figure 4C). These functionalities of fibroblasts could be further supported by evaluating the DGE of various transcripts that had roles in wound healing, regulation of inflammation, and immune responses (Figure 4D). Our data also revealed detection of two transcripts (Cyb5a and Fxpyd6) (Figures 4E and 4H) that were not previously associated with the known cardiac injury marker Wif1, suggesting the possibility of their use as cardiac injury markers. Overall, although the composition of fibroblasts was complex, their transcriptome profiles revealed a diverse role in fibrosis, immune activation, and inflammatory functions.

Among other cells, NK cells and ILCs were low in number (Figure 1G, Table S2), but several genes with known functions were noted in NK cells, including Tmem176b, which has a role in fibrosis. Likewise, DGE analysis in the neutrophils revealed transcripts that have varied functions, such as negative regulation of inflammation, anti-viral response, antigen-presentation functions, and wound healing (Figure 6D). Notably, Ly6g+ cluster 2 had increased Il1b, indicating that neutrophils might release Il1b to promote pathogenic fibrosis and cardiac remodeling, eventually leading to DCM, thus providing more evidence to previous hypotheses (Frangogiannis, 2015; Rivadeneyra et al., 2018). However, the IL-17 signaling pathway appeared to be dominantly influenced by genes expressed in neutrophils, raising a question as to their reactivity to Th17 cytokines (Figures 6E and S6B). On the one hand, Th17 cells facilitate neutrophil chemotaxis, but on the other, neutrophils do not respond to IL-17, as they lack expression of IL-17Rc (Pelletier et al., 2010). It may be that the genes expressed in neutrophils may modulate the IL-17 signaling pathway in non-neutrophils. Unexpectedly, B cells formed a minor fraction in the heart infiltrates, and their number was significantly low in myocarditic mice (Figure 1G, Table S2). Analysis of their transcriptomes also did not reveal significant differences between myocarditic and healthy mice (Figure S7C), suggesting that B cells appear not to have a major role in chronic myocarditis.

Constructing the intercellular communication networks within the cardiac cellulome using ligand-receptor interactions in myocarditic mice indicated strong associations between CD4 and CD8 T cells, CD8 T cells and neutrophils, myeloid cells and NK cells, and fibroblasts and SMCs and myeloid cells, among others (Figures 7A and 7B). Likewise, the number of receptor-ligand pairs and interaction strengths were relatively high between fibroblasts and CD4 and CD8 T cells and between myeloid cells and CD4 T cells and SMCs (Figures 7C and S8B). In a virus-free setting of autoimmune myocarditis, it has been suggested that the Th17 cells-fibroblasts-monocytes/macrophage axis may be critical for the development of inflammatory cardiomyopathy (Wu et al., 2014), which may not be relevant to viral myocarditis. Based on our data, we noted a dependency between CD4 and CD8 T cells in their functionalities, and CD8 T cells may potentially influence the effects of other innate cells, including fibroblasts. Such a possibility can be expected with virus infections, as CD8 T cells form an important component of anti-viral responses. Furthermore, by analyzing the signaling molecules in relation to receptor-ligand interactions, we noted distinct pathways important in reparative processes (Figure 7F). Increased EGF, LCK, ncWNT, and NOTCH signaling (Figures 7F and S10), which have been reported to be necessary for CVB3 replication and cardiac fibrosis, could indicate that these pathways could be promoting the anti-inflammatory reparative process in the myocardium, as well as cardiac remodeling in the postinfectious myocarditis phase (Bergmann, 2010; Ferrari and Rizzo, 2014; Liu et al., 2000; Makki et al., 2013). OSM’s role in cardiomyocyte de-differentiation has been found to be increased in mouse models of myocardial infarction and in end-stage heart failure patients (Kubin et al., 2011). Even though OSM plays a protective role during the acute phase of myocardial damage, its prolonged expression, along with the increased infiltration of macrophages during the chronic phase, may promote functional deterioration and loss of cardiac contractility, leading to DCM and heart failure. PDGF, which has been known to be involved in fibrosis and upregulated in CVB3 infection, showed increased signaling in the postinfectious myocarditis phase. Profibrotic Tgfbi was upregulated in monocyte cluster 0 and Mfap4+ in fibroblast cluster 6 and is known to promote PDGF signaling, thus suggesting a role for it in the cardiac fibrosis and remodeling phase of post-infectious myocarditis. Overall, our data indicate that T cells, myeloid cells, and fibroblasts may mainly contribute to the progression of viral myocarditis in chronically infected mice and that the upregulated signaling pathways could be targeted for therapeutic purposes in DCM.

Finally, in our efforts to identify the potential TFs that might regulate various cellular functions, we identified Ets1 in T cells and Irf5 and Mafb in myeloid cells, which have varied roles (inflammation, differentiation of monocytes/macrophages and cardiac morphogenesis) as potential target candidates TFs (Figure 8A). In an AngII mouse model of cardiac fibrosis, a recent report showed that deletion of Ets1 from ECs reduced cardiac fibrosis and hypertrophy (Xu et al., 2019). Upregulated expression of Ets1 in our data, which was regulating the expression of target genes implicated in cardiac fibrosis and heart failure, may mean that targeting this TF could be beneficial in attenuating the transition to DCM. Although the roles for Irf5 and Mafb have not been defined clearly in myocarditis, our data indicate their important role in controlling the expression of M2-specific genes and transcripts involved in cardiac fibrosis/remodeling/heart failure during the post-infectious phase of myocarditis. As these TFs can act as activators or repressors of various genes, determining the myocarditis phenotype in mice deficient for TFs may permit us to evaluate their roles in viral myocarditis and develop strategies for therapeutic interventions.

In summary, we have described the cellular compositions and their transcriptome profiles in heart infiltrates using the mouse model of viral myocarditis. Although T cells, myeloid cells, and fibroblasts formed a major component, the proportion of B cells was low. Although CVB3 infection has been extensively used to understand the pathogenesis of myocarditis, a long-standing question remains as to the underlying mechanisms of chronic myocarditis, with autoimmune theory as one possibility. Reports indicate detection of autoantibodies in CVB3 myocarditis, but their pathogenic role remains inclusive. We have been investigating the role of autoreactive T cells and, using MHC class II dextramers, have reported the appearance of pathogenic autoreactive T cells with specificities for multiple antigens (Basavalingappa et al., 2020; Gangaplara et al., 2012). Our scRNAseq data also point to a role for T cells, but in-depth analysis of their role can be investigated in antigen-specific (dextramer+) T cells by scRNAseq analysis. Likewise, detection of M2 cells was not surprising because of their role in reparative functions, and detection of fibroblasts was also expected due to their role in the formation of fibrosis. In adjuvant-induced myocarditis, it was recently shown that neutrophils, Th17 cells, Treg cells, and γδ T cells were found predominantly during induction, acute, subacute, and chronic (myopathy) phases, whereas macrophages were found in all phases (Hua et al., 2020). But their transcriptome signatures were different as compared with our study. For example, Hif1a-producing Th17 cells were found critical in the inflammatory response in the above model that was missing in our analysis. Such variations are expected between models involving viruses as in our studies versus virus-free systems described earlier. In addition, our analysis involved a single time point measurement that may not truly reflect the infiltrations of complete progression of disease. More importantly, induction of autoimmune myocarditis involves the use of two doses of complete Freund’s adjuvant (CFA) (Smith, 1999). Because CFA is a powerful immune-stimulating adjuvant, a proportion of the immune-stimulating genes might have been influenced by nonspecific effects of CFA. Of note, we could not study the compositions of nonimmune cells such as cardiomyocytes, SMCs, pericytes, and ECs in our studies, and the use of their purified populations along with spatial transcriptomics may yield new insights into viral pathogenesis. For example, although low in number, we investigated the transcriptome profile of ECs, leading us to note upregulation of TFs, which have a role in apoptosis of ECs and inflammation (Atf3, Irf1, and Stat6), with a corresponding downregulation of those involved in EC survival and angiogenesis (Myc and Nfe212) in myocarditic animals as compared with controls (data not shown). Overall, our scRNAseq analysis offers a new dimension to understanding the post-infectious phase of viral myocarditis and associated pathogenic cardiac remodeling. Likewise, analysis of cell types using tissues derived from specific areas of heart could yield new insights. For example, scRNAseq analysis of noncardiomyocyte fraction from myocardial infarction model led to the identification of myofibroblast subtypes expressing pro- or anti-fibrotic signatures (Farbehi et al., 2019). Similarly, transcriptomics of infarcted area in the ischemic injury model led to the determination of the significance of beta-2 microglobulin in the activation of fibroblasts in response to ischemic damage (Molenaar et al., 2021). Such studies could also lead to the identification of various cardiac-derived fibroblasts that may help in predicting the disease outcome of postmyocardial infarction (Forte et al., 2020). The scRNAseq analysis may also facilitate to the determination of relevance of immune cells in settings that do not involve infectious agents. For example, transcriptomic profiles in the pressure-overload-driven heart failure model resulted in the identification of innate and adaptive immune cell infiltrations with a role for oncostatin M in pro-inflammatory macrophages and PD-1 in Treg cells (Martini et al., 2019).

Limitations of the study

Our studies had three major limitations. (1) We could not perform the scRNAseq analysis at multiple time points in CVB3 infection to relate the disease phenotypes with the molecular signatures of different cell types during the disease progression. Such an analysis, in combination with the use of tools such as echocardiography to assess the functional states of hearts, may be helpful to relate the disease stage with scRNAseq data. (2) Our cell isolation procedures did not permit us to evaluate single cell analysis in nonimmune cells such as cardiomyocytes, SMCs, pericytes, and ECs. Because CVB3 is a cardiotropic virus that causes lysis of cardiac cells, determination of transcriptome analysis in the above cell types may lead to the identification of new molecular pathways of pathogenic or therapeutic importance. (3) Although, our data revealed a role for several genes, we did not investigate their functional significance in the pathogenesis of viral myocarditis that may require the use of genetically altered mice or pharmacological interventions.

STAR★Methods

Key resources table

Resource availability

Lead contact

Requests for information or reagents should be directed to the lead contact, Jay Reddy (nreddy2@unl.edu).

Materials availability

This study did not generate new unique reagents.

Experimental model and subject details

Mice

Six-to-eight-week-old, and nine-11-week-old male A/J mice (H-2a) were procured from Jackson Laboratory (Bar Harbor, ME) and maintained according to the Institutional Animal Care and Use Committee’s guidelines of the University of Nebraska-Lincoln (protocol #: 1904), Lincoln, NE. Infection studies were performed based on biosafety level 2 guidelines. When animals were found to have persistent clinical signs, such as failure to move when physically touched or prodded, or failure to eat or drink, they were euthanized using a carbon dioxide chamber as recommended by the Panel on Euthanasia of the American Veterinary Medical Association.

Method details

Virus propagation and infection

The Nancy strain of CVB3 was procured from the American Type Culture Collection (ATCC, Manassas, VA, USA), and the virus was titrated in Vero cells (ATCC). The adherent Vero cells were grown to 80 to 90% confluence in 75cm2 flasks in EMEM/10% fetal bovine serum (FBS) and were later infected with CVB3 with multiplicity of infection 1 in EMEM containing no FBS. After incubation at 37⁰ C for 1 hour with gentle intermittent rocking, maintenance medium (EMEM/2% FBS) was added. Based on the cytopathic effect of virus during the next 1 to 2 days, supernatants containing virus were harvested. After determining 50% tissue culture infective dose (TCID50) values based on the Reed-Muench method, the virus stocks were aliquoted and preserved at −80°C. To infect mice, virus stock diluted in 1x PBS to contain 10,000 TCID50 in 200 μl was administered intraperitoneally (i.p.). Animals were monitored closely, cages were changed once in 2 days, and body weights were taken daily until termination. In addition, an alternative food and fluid source, trans gel diet (ClearH2O, Portland, ME, USA), was placed on the cage floor as needed.

Heart single-cell preparation

Single-cell suspensions from mouse hearts were prepared as previously described (Pinto et al., 2013, 2016). Briefly, male CVB3-infected mice and their age-matched healthy control mice were euthanized on day 21 post-infection using 2% CO2 in an asphyxiation chamber as per the IACUC guidelines. The hearts were perfused using perfusion buffer (1 × DPBS with 0.8 mM CaCl2, 5 ml/min for 5 minutes) until the liver was completely blanched and appeared pale yellow/brown in color. Next, hearts were isolated, their atria and valves were removed, and the whole heart was minced to ∼0.5-1 mm cubes using surgical scissors. Minced heart tissue was digested in 3 mL of digestion buffer (2 mg/mL Collagenase IV [Worthington Biochemical], 1.2 units/mL Dispase II [Sigma-Aldrich], in perfusion buffer) for ∼45 minutes at 37°C using a rotating holder, with tissue suspension triturated once in every 15 minutes with 1000 μl wide-bore micropipette tips. The cell suspension was filtered through a 70 μm nylon filter mesh to remove residual undigested tissue pieces. The filtrate was then diluted in ∼15 mL perfusion buffer, and the cells were pelleted at ∼200 × G for 20 minutes with no centrifuge brakes engaged. Cell supernatant was then aspirated, and the pellet was re-suspended in ∼15 mL of 1 × HBSS (Sigma-Aldrich) + 0.8 mM CaCl2. The cells were pelleted again as described above. In order to remove unwanted debris, a debris removal kit was used as per the manufacturer’s guidelines (Miltenyi Biotec, San Diego, CA). The final debris-free cell pellet was resuspended in 1000 μL of 2% FBS in RPMI medium for downstream cell sorting using flow cytometry.

Flow cytometry and sorting

Freshly prepared cardiac single-cell suspensions from healthy and myocarditic hearts were subjected to surface staining with Annexin V (Biolegend, San Diego, CA) and propidium iodide (PI; Biolegend). In brief, cells were first washed with Annexin V binding buffer, followed by surface staining with Annexin V (1:200, vol/vol) and PI (1:100, vol/vol) at room temperature (25°C) for 15 minutes in the dark. Cells were then resuspended and sorted by flow cytometry (FACSAria II, BD Biosciences, San Jose, CA). Only singlets that were viable and non-apoptotic (Annexin V- PI-) were sorted and collected in tubes containing RPMI with 2% FBS.

Sample processing and sequencing

Two replicates, with n=7 mice per treatment group, were used for heart sample processing. Approximately 16,000 cells were loaded onto a single channel of the 10X Genomics Chromium controller (10X Genomics, Pleasanton, CA), with a target recovery of ∼10,000 cells using the chromium v2 and v3 single-cell reagent kit. After the generation of single-cell gel bead-in-emulsions, cDNA was synthesized using a C1000 Touch Thermal Cycler (Bio-Rad Laboratories, Hercules, CA) and amplified for 11 cycles as per the manufacturer’s protocol. Quality control (QC) and quantification were performed using the Agilent 2100 bioanalyzer (Agilent Technologies, Santa Clara, CA) per the manufacturer’s guidelines. Amplified cDNA (50 ng) was used to construct 3’ expression libraries, and the libraries were pooled and run on an Illumina HiSeq 4000. Each lane consisted of 150 base-pair, paired-end reads. The Illumina basecall files were converted to FASTQ format, and these files were aligned to the murine genome (mm10) using the CellRanger v3.0.2 pipeline as described by the manufacturer. Across aligned cells, the mean number of reads per cell was 39,923, with an average of 95.3% of reads mapped to the mm10 genome.

Single-cell data processing and analysis

Initial processing of cells isolated from the heart in myocarditis run 1 (n=2,617), myocarditis run 2 (n=10,618), control run1 (n=1,528), and control run 2 (n=8,201) was performed using the Seurat R package (v3.0.2) (Butler et al., 2018; Stuart et al., 2019). Samples were normalized using the sctransform approach (Hafemeister and Satija, 2019) with default settings. The transformed data was then formed into a single data set using canonical correlational analysis and mutual nearest neighbors (MNN) as described by Stuart et al. (Stuart et al. (2019). Dimensional reduction to form the uniform manifold approximation and project (UMAP) utilized the top 30 calculated dimensions and a resolution of 0.6. The schex R package (v1.1.5) was used to visualize mRNA expression of lineage-specific or highly differential markers by converting the UMAP manifold into hexbin quantifications of the proportion of single cells with the indicated gene expressed. Default binning was set at 80 unless otherwise indicated in the figure legend.

Cell type identification utilized the SingleR (v1.0.1) R package (Aran et al., 2019) with correlations of the single-cell expression values with transcriptional profiles from pure cell populations in the Immgen (Heng et al., 2008). In addition to correlations, canonical markers for cell lineages were utilized and are available in Table S1. Differential gene expression utilized the Wilcoxon rank-sum test on count-level mRNA data. For differential gene expression across clusters or subclusters, the FindAllMarkers function in the Seurat package was used, employing the log-fold change threshold >0.25, minimum group percentage = 10%, and pseudocount = 0.1. Differential comparisons between conditions utilized the FindMarkers function in Seurat without filtering and a pseudocount = 0.1. Multiple hypothesis correction was reported using the Bonferroni method. Cell cycle regression was performed in Seurat using the CellCycleScoring function and genes derived from Nestorowa et al. (2016). Genes were isolated by calling cc.genes.updated.2019 in R and then converting into murine nomenclature. Gene set enrichment analysis was performed using the escape R package (v1.0.1). Differential enrichment analysis was performed using the getSignificance function in escape that is based on the limma R package linear fit model.

Cell trajectory analysis

Cell trajectory analysis used the Slingshot (v1.6.0) R package with default settings for the slingshot function, and using the UMAP embeddings from the subclustering for each cell type. Ranked importance of genes was calculated using the top 300 variable genes, and rsample (v0.0.9) and tidymodels (v0.1.0) R packages were used to generate random forest models based on a training data set of 75% of the cells. The rand_forest function in the parsnip (v0.1.1) R package was used, with mtry set to 200, trees to 1400, and minimum number of data points in a node equal to 15 across all cell types.

GO and pathway enrichment analysis of DEGs

GO pathway enrichment analysis of myocarditis-related DEGs was performed by Metascape (http://metascape.org/gp/index.html) (version 3.5) (Zhou et al., 2019). Results were visualized using the ggplot2 R package (version 3.2.1). Single-cell normalized enrichment scores were calculated using the escape (v1.0.1) R package (Borcherding et al., 2021). From this analysis, differentially expressed ligand and receptor between myocarditic and healthy controls for indicated cell types were extracted to use for the size/count of the dot plot. Differential gene set enrichment utilized the Welch’s T test with the Bonferroni adjustment for multiple hypothesis correction comparing individual cells in myocarditis versus healthy controls.

Intercellular communication analysis

Cell-cell interactions based on the expression of known ligand-receptor pairs in different cell types were inferred using CellChat (version 1.0.0) R package (Jin et al., 2021). We used the default settings to predict major signaling interactions of cells and how these cells and signals coordinate various functions. In brief, we followed the workflow recommended in CellChat and loaded the normalized counts into CellChat and applied the preprocessing functions identifyOverExpressedGenes, identifyOverExpressedInteractions, and projectData with default parameters set. For the analysis of ligand-receptor interactions, the functions computeCommunProb, computeCommunProbPathway, and aggregateNet were applied using default parameters. Finally, we classified signaling pathways and depicted conserved and context-specific pathways between myocarditis and healthy hearts.

Analysis of TF regulatory network

TF regulatory network analysis was performed using SCENIC (Aibar et al., 2017) (version 1.1.2.2) with default parameters and a co-expression method set to top 50 results per target. Murine mm9 TFs were downloaded using RcisTarget (version 1.6.0) as a reference. Enriched TF-binding motifs predicted candidate target genes (regulons), and regulon activity was inferred by RcisTarget. The resulting AUC enrichments for individual cells were attached to the Seurat object, and the median by cluster and condition were visualized using the pheatmap R package (v1.0.12). Transcription factor regulons were concatenated across the results from the motif enrichment step and graphed using the igraph (v1.2.6) R package.

RNA isolation and real-time quantitative PCR

Hearts harvested from myocarditis and control mice were stored at −80°C. Approximately 20-30 mg of tissue were transferred to the RLT buffer and homogenized with a FastPrep96 system as recommended (Lysing Matrix D 1.4-mm ceramic beads; MP Biomedicals, Irvine, CA). RNA was isolated using the RNeasy kit (Qiagen, Hilden, Germany), and samples were treated with deoxyribonuclease (DNase) I and quantified using the NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA). One-step qPCR was performed using the TaqManTM Fast Virus one-step Master mix (Applied Biosystems). The real-time quantitative PCR analysis was done using the TaqMan Fam gene expression assays that included amplifications for various target genes namely Ccl24 (Mm00444701_m1), Mt1 (Mm00496660_g1), Chil3 (Mm00657889_mH), S100a4 (Mm00803372_g1), S100a6 (Mm00771682_g1), Gatm (Mm01268677_m1), Nkg7 (Mm00452524_g1), S100a11 (Mm07297383_g1), Wif1 (Mm00442355_m1), Cyb5a (Mm00518027_m1), Fxyd6 (Mm00445583_m1), Thbs4 (Mm03003598_s1), and the housekeeping gene glyceraldehyde-3-phosphate dehydrogenase (Gapdh; Mm99999915_g1) using the CFX96 Touch Real-time PCR detection system (BioRad). Expression of target genes was normalized to Gapdh using the 2−(ΔΔCt) method.

Quantification and statistical analysis

Statistical Analyses were performed in R (v3.6.3). Two-sample significance testing utilized Welch’s T-test, with significance testing for more than three samples utilizing one-way analysis of variance (ANOVA) with Tukey honest significance determination for correcting multiple comparisons. Two-proportion Z-tests were performed using the total number of cells in each condition as the number of trials and without a prior for proportion. Mann-Whitney test was used to determine the significance between groups for qPCR data. Graphs were prepared by GraphPad prism software v9.0.2 (GraphPad Software, Inc., La Jolla, CA).

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
Annexin V-APC	BioLegend	Cat# 640919
Propidium Iodide-PE	BioLegend	Cat# 421301
Bacterial and virus strains
CVB3-Nancy	ATCC	Cat# VR-30
Chemicals, peptides, and recombinant proteins
Type IV Collagenase	Worthington Biochemical	Cat# LS004188
Dispase II	Sigma-Aldrich	Cat# D4693
Debris Removal Solution	Miltenyi Biotec	Cat# 130-109-398
Lysing Matrix D 1.4 mm ceramic beads	MP Biomedicals	Cat# 116913050-CF
DNase I, Amplification Grade	Invitrogen	Cat# 18068015
Annexin V Binding Buffer	BioLegend	Cat# 422201
Critical commercial assays
RNeasy Total RNA Isolation Kit	Qiagen	Cat# 74104
TaqManTM Fast Virus 1-step Master mix	Applied Biosystems	Cat# 4444432
TaqManTM Gene Expression Assay - Ccl24	Applied Biosystems	Mm00444701_m1
TaqManTM Gene Expression Assay – Mt1	Applied Biosystems	Mm00496660_g1
TaqManTM Gene Expression Assay – Chil3	Applied Biosystems	Mm00657889_mH
TaqManTM Gene Expression Assay – S100a4	Applied Biosystems	Mm00803372_g1
TaqManTM Gene Expression Assay – S100a6	Applied Biosystems	Mm00771682_g1
TaqManTM Gene Expression Assay – Gatm	Applied Biosystems	Mm01268677_m1
TaqManTM Gene Expression Assay – Nkg7	Applied Biosystems	Mm00452524_g1
TaqManTM Gene Expression Assay – S100a11	Applied Biosystems	Mm07297383_g1
TaqManTM Gene Expression Assay – Wif1	Applied Biosystems	Mm00442355_m1
TaqManTM Gene Expression Assay – Cyb5a	Applied Biosystems	Mm00518027_m1
TaqManTM Gene Expression Assay – Fxyd6	Applied Biosystems	Mm00445583_m1
TaqManTM Gene Expression Assay – Thbs4	Applied Biosystems	Mm03003598_s1
TaqManTM Gene Expression Assay – Gapdh	Applied Biosystems	Mm99999915_g1
Chromium™ Single Cell 3’ Library & Gel Bead Kit v2	10X Genomics	Cat# 120237
Chromium™ Single Cell A Chip Kit, 16 rxns	10X Genomics	Cat# 1000009
Chromium™ Single Cell 3’ Library & Gel Bead Kit v3	10X Genomics	Cat# 1000075
Chromium™ Single Cell B Chip Kit, 16 rxns	10X Genomics	Cat# 1000154
Chromium™ i7 Multiplex Kit	10X Genomics	Cat# 120262
Deposited data
Raw sequencing data files for scRNA sequencing	This paper	GEO: GSE174458
Experimental models: Cell lines
Vero cells	ATCC	Cat# CCL-81
Experimental models: Organisms/strains
Mouse: A/J	The Jackson Laboratory	Cat# 000646
Software and algorithms
GraphPad Prism (v9.0.2)	GraphPad Software	https://www.graphpad.com/
FlowJo (v10.2)	BD Biosciences	https://www.flowjo.com/
Microsoft Excel	Microsoft	https://www.microsoft.com/en-us/microsoft-365/excel
RStudio (v3.6.3)	RStudio, Inc.	https://www.rstudio.com/
Cell Ranger (v3.0.2)	10x Genomics	https://support.10xgenomics.com/
Seurat (v3.0.2)	(Stuart et al., 2019)	https://satijalab.org/seurat/
schex (v1.1.5)	Saskia Freytag	https://github.com/SaskiaFreytag/schex
SingleR (v1.0.1)	(Aran et al., 2019)	https://github.com/dviraran/SingleR
escape (v1.0.1)	(Borcherding et al., 2021)	https://github.com/ncborcherding/escape
Slingshot (v1.6.0)	(Street et al., 2018)	https://github.com/kstreet13/slingshot
Metascape (v3.5)	(Zhou et al., 2019)	http://metascape.org/gp/index.html
ggplot2 (v3.3.5)	Tidyverse	https://github.com/tidyverse/ggplot2/
CellChat (v1.0.0)	(Jin et al., 2021)	https://github.com/sqjin/CellChat
SCENIC (v1.1.2.2)	(Aibar et al., 2017)	https://github.com/aertslab/SCENIC
igraph (v1.2.6)	https://igraph.org/	https://github.com/igraph/rigraph
pheatmap (v1.0.12)	Raivo Kolde	https://github.com/raivokolde/pheatmap
Other
Immunological genome project (ImmGen)	(Heng et al., 2008)	https://www.immgen.org/
Scripts for scRNA seq processing	This paper	https://github.com/ncborcherding/Lasrado_Myocarditis