Weighted gene co-expression network analysis identified underlying hub genes and mechanisms in the occurrence and development of viral myocarditis.

BACKGROUND: Myocarditis is an inflammatory myocardial disease, which may lead to heart failure and sudden death. Despite extensive research into the pathogenesis of myocarditis, effective treatments for this condition remain elusive. This study aimed to explore the potential pathogenesis and hub genes for viral myocarditis.
METHODS: A weighted gene co-expression network analysis (WGCNA) was performed based on the gene expression profiles derived from mouse models at different stages of viral myocarditis (GSE35182). Functional annotation was executed within the key modules. Potential hub genes were predicted based on the intramodular connectivity (IC). Finally, potential microRNAs that regulate gene expression were predicted by miRNet analysis.
RESULTS: Three gene co-expression modules showed the strongest correlation with the acute or chronic disease stage. A significant positive correlation was detected between the acute disease stage and the turquoise module, the genes of which were mainly enriched in antiviral response and immune-inflammatory activation. Furthermore, a significant positive correlation and a negative correlation were identified between the chronic disease stage and the brown and yellow modules, respectively. These modules were mainly associated with the cytoskeleton, phosphorylation, cellular catabolic process, and autophagy. Subsequently, we predicted the underlying hub genes and microRNAs in the three modules.
CONCLUSIONS: This study revealed the main biological processes in different stages of viral myocarditis and predicted hub genes in both the acute and chronic disease stages. Our results may be helpful for developing new therapeutic targets for viral myocarditis in future research. 2020 Annals of Translational Medicine. All rights reserved.

Introduction

Myocarditis is an inflammatory myocardial disease, which is mainly caused by viral infections, such as those caused by coxsackieviruses, adenoviruses, influenza viruses, cytomegaloviruses, and human immunodeficiency virus (1). Pathologically, myocarditis is characterised by infiltration of the myocardium with mononuclear cells (2). Most myocarditis patients exhibit mild symptoms. However, a subset of patients experience fulminant myocarditis, leading to heart failure, arrhythmia, fulminant haemodynamic collapse and sudden mortality (3-5), whereas about one third of viral myocarditis patients develop a chronic form of the disease, which may lead to dilated cardiomyopathy (DCM) (1).

As we know, direct injury caused by the virus and indirect injury mediated by the host’s immune response underlie the pathogenesis of viral myocarditis (6). The progression of viral myocarditis is characterised by 3 stages: (I) an acute stage triggered by viral entry and replication, (II) a subacute stage characterised by inflammatory cell infiltration, and (III) a chronic stage characterised by cardiac remodelling (7). However, the mechanism underlying fulminant myocarditis and how it progresses into dilated cardiomyopathy still require further elucidation.

Conventional molecular biology techniques can uncover the specific function of a gene, but do not reveal interaction between genes in biological networks (8,9). To our knowledge, there is still a lack of systematic biological network analyses of hub genes in viral myocarditis. Weighted gene co-expression network analysis (WGCNA), a widely used systems biology approach, is a powerful tool that is used to define correlation patterns between genes (10-15). Here, we constructed a co-expression network using expression data from the Gene Expression Omnibus (GEO) database. Genes with similar expression patterns were clustered into the same modules. The relationships between different stages of the myocarditis and modules were calculated to identify highly related modules. The turquoise module was positively correlated with the acute disease stage, while the brown and yellow modules were positively and negatively correlated with the chronic disease stage, respectively. The hub genes and miRNAs identified in these modules could serve as biomarkers and therapeutic targets for viral myocarditis.

Methods

Data collection

The heart RNA expression profiles from mouse models of viral myocarditis were downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/; GSE35182) (16). This dataset was generated on the GPL6246 platform of the Affymetrix Mouse Gene 1.0 ST Array. To explore the occurrence and development of viral myocarditis, we analysed the RNA expression profiles in an acute infection group (coxsackievirus B3 infection for 10 days), chronic infection group (coxsackievirus B3 infection for 90 days) and normal group (PBS injection).

Data pre-processing

The R version 3.6.1 software was used for data processing and analysis. The raw data in CEL format were read using the Affy package version 1.62.0 (17) and processed using the robust multiarray average (RMA) algorithm () (18). Next, we used the impute.knn function in the impute package version 1.58.0 to impute the individual missing values in the raw data. This function used the Euclidean metric to find k nearest genes that were similar to the expression profiles of each gene with a missing value, and estimated the missing value through the expression values of nearest genes (19,20). After processing, the 6,680 genes exhibiting the most significant expression changes (the top 25% of rank genes with the largest variance) were selected for further analysis.

Figure 1

Data pre-processing and construction of a co-expression network. (A) Box plot of the expression data before RMA normalisation. (B) Box plot of the expression data after RMA normalisation. (C) Analysis of the scale-free fit index for various soft-threshold powers, the red line was set at 0.90. (D) Analysis of the mean connectivity for various soft-threshold powers. (E) The cluster dendrogram of genes in the co-expression network. RMA, robust multiarray average; N, Normal group; A, Acute infection group; C, Chronic infection group.

Co-expression network construction

The WGCNA package version 1.68 in R software was used to construct the co-expression network (21,22). An appropriate soft threshold was selected to ensure scale-free topology (R2>0.9). The topological overlap matrix (TOM) was constructed to measure the network connectivity of the genes. Genes with similar patterns were clustered into the same modules (minimum size =30) using average linkage hierarchical clustering. The relationships between phenotypes and modules were calculated to identify highly related modules. Finally, the highly correlated modules were analysed to explore their potential roles. In addition, the gene expression profiles of the highly correlated modules were visualised using the R software.

Functional annotation

The functional annotation analyses of the highly correlated modules were done using the clusterProfiler package version 3.12.0 on R (23,24). The P value was adjusted by the Holm–Bonferroni method. The adjusted P value cutoff was set at 0.05 and the q-value cutoff was set at 0.2. An adjusted P value less than 0.05 was considered to be significant, and the identified significant analyses were sorted by gene counts.

Identification of hub genes

Hub genes are defined as the genes that exhibit the greatest association with a disease. The IC for each gene was calculated by summing the connection strengths with other genes in the same module. For a given gene, the higher the IC, the stronger its relationship with other genes, and the more important it is within the module. Therefore, hub genes are the genes that possess high IC values. The interaction networks between hub genes and other genes in the module were visualised using Cytoscape version 3.7.2 (25).

Prediction of potential miRNA targeting networks

The top 100 genes ranked by IC were selected from each module and a miRNA targeting network was predicted by miRNet analysis (www.mirnet.ca) (18). The network was visualised using Cytoscape version 3.7.2.

Statistical analysis

WGCNA and functional annotation were performed in R (version 3.6.1) with default test statistics and cutoff values as specified in individual Methods sections. P<0.05 was considered statistically significant.

Results

Construction of weighted co-expression network and identification of key modules

We used the 6,680 genes that exhibited the greatest differences in expression to construct weighted gene co-expression networks. A total of 18 samples were used as input for the hierarchical clustering analysis. The hclust function was then used to cluster the samples and visualise outliers (18). Because this analysis did not reveal any outliers, all the samples were included in the downstream analyses (Figure S1). Before constructing the weighted co-expression matrix, a soft-threshold β was calculated to ensure a scale-free topology. When the soft-threshold β was equal to 14 the independence degree rose to 0.964. Therefore, co-expression gene modules were constructed using β=14 ().

Next, the co-expression network was constructed using a one-step method and the correlation matrix was constructed to calculate the correlation efficiency between genes. Genes with similar expression patterns were clustered into the same module. A total of 12 modules were identified (). Next, module eigengenes (MEs), components of a module representative of the gene expression profiles in that module, were calculated for each module. Subsequently, the correlation between MEs and viral myocarditis at different stages was calculated, and a correlation coefficient greater than, or equal to, 0.9 or less than, or equal to, −0.9 (| r | ≥0.9) was considered to be significant. In acute viral myocarditis, the turquoise module (r=0.9; P=5×10−7) was identified as the statistically significant module. In chronic viral myocarditis, the brown (r=0.92; P=5×10−8) and yellow (r=−0.95; P=4×10−9) modules were identified as the statistically significant modules (). These modules were then selected for further analyses. To evaluate the correlation between the different stages of viral myocarditis and the statistically significant modules, we calculated the values between module membership (MM) and gene significance (GS). The turquoise module (cor =0.94; P<1×10−200) showed a high positive correlation with acute viral myocarditis (), whereas the brown module (cor =0.9; P<1×10−200) showed a high positive correlation with chronic viral myocarditis. A high negative correlation (cor =−0.96; P<1×10−200) was observed between the yellow module and the chronic disease stage ().

Figure 2

Identification of key modules. (A) Module-trait relationships in the constructed network. The upper figure in each row represents the correlation with the different stages of myocarditis while the lower figure represents the P value. (B) The MM versus GS plot of the most positively related modules in the acute phase. (C) The MM versus GS plot of the most positively related modules in the chronic phase. (D) The MM versus GS plot of the most negatively related modules in the chronic phase. MM, module membership; GS, gene significance.

Genes expression of key modules at the different stages of viral myocarditis

To further explore the role of key modules in the occurrence and development of viral myocarditis, the gene expression in these modules at the different stages of viral myocarditis were visualised using a heatmap. Gene expression in the turquoise module was significantly increased in the acute infection group (), indicating that genes in this module play an important role in the acute stage of viral myocarditis. In chronic viral myocarditis, gene expression in the brown module was significantly increased, suggesting that the genes in this module play an important role in chronic viral myocarditis (). In the yellow module, gene expression continued to fall reaching its lowest level in the chronic phase, indicating that genes in this module are important for the development of viral myocarditis ().

Figure 3

Gene expression profiles of key modules. (A) Gene expression of the turquoise module in each sample. (B) Gene expression of the brown module in each sample. (C) Gene expression of the yellow module in each sample. In (A), (B), and (C), the upper part is a heatmap (the red and green colours correspond to high and low expression values, respectively), while the lower part is the module eigengene (a representative of the gene expression profiles). N, Normal group; A, Acute infection group; C, Chronic infection group.

Functional annotation of the key modules

As the turquoise, brown and yellow modules were identified as key modules, functional annotation was adopted to characterise them. Gene Ontology (GO) term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were executed using clusterProfiler. To avoid inaccurate enrichment results caused by too many genes, we selected the top 25% genes ranked by IC in each module for enrichment analyses.

The functional enrichment of GO term categories includes biological process (BP), cellular component (CC), and molecular function (MF). Analysis of BP in the turquoise module revealed that the genes were enriched in T cell activation, regulation of immune effector processes, negative regulation of immune system processes, positive regulation of cytokine production, and lymphocyte mediated immunity. Regarding the CC and MF categories, the genes were mainly associated with the lysosome (CC), lytic vacuole (CC), enzyme activator activity (MF) and GTPase activity (MF) (). The most enriched BPs in the brown module were regulation of actin filament-based processes, muscle organ development, blood circulation, circulatory system process, and regulation of actin cytoskeleton organization. In the CC and MF categories, genes were mainly associated with the extracellular matrix (CC), anchoring junction (CC), enzyme activator activity (MF), and actin binding (MF) (). The most enriched BPs in the yellow module were negative regulation of phosphorylation, positive regulation of catabolic process, epithelial cell proliferation, regulation of vasculature development, and autophagy. In the CC and MF categories, genes were mainly associated with transport vesicles (CC), cell-cell junction (CC), and ion channel binding (MF) (). Detailed information regarding the GO analyses is listed in .

Figure 4

GO and KEGG pathway enrichment analysis of the genes in the key modules. (A) GO enrichment analysis of the genes in the turquoise module. (B) GO enrichment analysis of the genes in the brown module. (C) GO enrichment analysis of the genes in the yellow module. (D) KEGG pathway enrichment analysis of the genes in the key modules. GO, gene ontology; BP, biological process; CC, cellular component; MF, molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Table 1

Gene ontology (GO) enrichment analyses of turquoise module

GO category	GO term	ID	Adjusted P	Count
BP	T cell activation	GO:0042110	1.51E-31	65
BP	Regulation of immune effector process	GO:0002697	4.68E-34	64
BP	Negative regulation of immune system process	GO:0002683	3.56E-28	60
BP	Positive regulation of cytokine production	GO:0001819	4.73E-26	56
BP	Lymphocyte mediated immunity	GO:0002449	2.22E-25	54
BP	Adaptive immune response based on somatic recombination of immune receptors built from immunoglobulin superfamily domains	GO:0002460	1.17E-24	54
BP	Response to virus	GO:0009615	6.48E-33	52
BP	Positive regulation of defense response	GO:0031349	3.27E-26	52
BP	Regulation of cell-cell adhesion	GO:0022407	2.33E-24	51
BP	Leukocyte cell-cell adhesion	GO:0007159	2.07E-27	50
MF	Enzyme activator activity	GO:0008047	5.90E-05	26
MF	GTPase activity	GO:0003924	2.30E-06	23
MF	Amide binding	GO:0033218	1.69E-04	22
MF	Peptide binding	GO:0042277	1.31E-04	20
MF	Ubiquitin-like protein ligase binding	GO:0044389	2.28E-04	20
CC	Lysosome	GO:0005764	3.35E-19	47
CC	Lytic vacuole	GO:0000323	3.35E-19	47
CC	Actin cytoskeleton	GO:0015629	3.36E-06	29
CC	Membrane raft	GO:0045121	6.49E-08	28
CC	Membrane microdomain	GO:0098857	6.49E-08	28

BP, biological process; CC, cellular component; MF, molecular function; Adjusted P, adjusted P value.

Table 2

Gene ontology (GO) enrichment analyses of brown module

GO category	GO term	ID	Adjusted P	Count
BP	Regulation of actin filament-based process	GO:0032970	9.90E-05	14
BP	Muscle organ development	GO:0007517	6.17E-04	13
BP	Blood circulation	GO:0008015	2.63E-03	12
BP	Circulatory system process	GO:0003013	2.63E-03	12
BP	Regulation of actin cytoskeleton organization	GO:0032956	1.43E-03	11
BP	Eye development	GO:0001654	2.63E-03	11
BP	Actin filament organization	GO:0007015	2.63E-03	11
BP	Visual system development	GO:0150063	2.63E-03	11
BP	Sensory system development	GO:0048880	2.63E-03	11
BP	Ameboidal-type cell migration	GO:0001667	2.63E-03	11
MF	Enzyme activator activity	GO:0008047	2.67E-03	11
MF	Actin binding	GO:0003779	4.36E-03	10
MF	DNA-binding transcription activator activity RNA polymerase II-specific	GO:0001228	4.44E-03	10
MF	GTPase regulator activity	GO:0030695	2.67E-03	9
MF	Nucleoside-triphosphatase regulator activity	GO:0060589	2.67E-03	9
CC	Extracellular matrix	GO:0031012	1.09E-02	10
CC	Anchoring junction	GO:0070161	1.09E-02	8
CC	Sarcomere	GO:0030017	9.14E-03	7
CC	Contractile fiber part	GO:0044449	9.14E-03	7
CC	Myofibril	GO:0030016	9.56E-03	7

BP, biological process; CC, cellular component; MF, molecular function; Adjusted P, adjusted P value.

Table 3

Gene ontology (GO) enrichment analyses of yellow module

GO category	GO term	ID	Adjusted P	Count
BP	Negative regulation of phosphorylation	GO:0042326	1.33E-02	13
BP	Positive regulation of catabolic process	GO:0009896	1.33E-02	12
BP	Epithelial cell proliferation	GO:0050673	1.33E-02	12
BP	Regulation of vasculature development	GO:1901342	1.33E-02	11
BP	Positive regulation of cellular catabolic process	GO:0031331	1.33E-02	11
BP	Autophagy	GO:0006914	1.34E-02	11
BP	Process utilizing autophagic mechanism	GO:0061919	1.34E-02	11
BP	Negative regulation of protein phosphorylation	GO:0001933	1.68E-02	11
BP	Blood circulation	GO:0008015	3.58E-02	11
BP	Small GTPase mediated signal transduction	GO:0007264	3.77E-02	11
MF	Ion channel binding	GO:0044325	2.29E-02	7
CC	Transport vesicle	GO:0030133	1.48E-02	10
CC	Cell-cell junction	GO:0005911	4.02E-02	10
CC	Adherens junction	GO:0005912	1.48E-02	9
CC	Anchoring junction	GO:0070161	1.48E-02	9
CC	Cell projection membrane	GO:0031253	3.07E-02	8

BP, biological process; CC, cellular component; MF, molecular function; Adjusted P, adjusted P value.

Next, we investigated the enriched pathways in the key modules. In the turquoise module, the enriched pathways were mainly associated with Herpes simplex virus 1 infection, Epstein-Barr virus infection and the phagosome. The brown module was enriched in the phosphatidylinositol 3-kinase (PI3K)/protein kinase B (AKT) and Ras-proximate-1 (Rap1) signaling pathways. The yellow module was enriched in the transcriptional misregulation in cancer (). Detailed information regarding the KEGG pathways analyses is shown in .

Table 4

KEGG pathway analyses of the key modules

Module	Pathway	ID	Adjusted P	count
turquoise	Herpes simplex virus 1 infection	mmu05168	1.37E-08	41
turquoise	Epstein-Barr virus infection	mmu05169	1.43E-16	40
turquoise	Phagosome	mmu04145	3.75E-15	34
turquoise	Influenza A	mmu05164	1.35E-14	32
turquoise	Cell adhesion molecules (CAMs)	mmu04514	1.99E-13	31
brown	PI3K-Akt signaling pathway	mmu04151	4.32E-02	8
brown	Rap1 signaling pathway	mmu04015	4.32E-02	6
brown	Regulation of actin cytoskeleton	mmu04810	4.32E-02	6
brown	ABC transporters	mmu02010	9.10E-04	5
yellow	Transcriptional misregulation in cancer	mmu05202	5.87E-03	9

KEGG, Kyoto Encyclopedia of Genes and Genomes; Adjusted P, adjusted P value.

In the turquoise module, integrin beta 2 (Itgb2, IC =610), tropomyosin 3, gamma (Tpm3, IC =607), Rho GTPase activating protein 30 (Arhgap30, IC =604), protein kinase C delta type (Prkcd, IC =596), and CD300C molecule 2 (CD300c2/AF251705, IC =596) were identified as hub genes (). In the brown module, hemecintin1 (Hmcn1, IC =186), pyruvate dehydrogenase phosphatase catalytic subunit 1 (Pdp1, IC =159), ATP-binding cassette, sub-family A (ABC1) member 8b (Abca8b, IC =158), Rab GTPase-activating protein 1-like (Rabgap1l, IC =157), and membrane-associated guanylate kinase WW and PDZ domain-containing protein 3 (Magi3, IC =157) were identified as hub genes (). In the yellow module, solute carrier family 9 (sodium hydrogen exchanger) member 3 regulator 2 (Slc9a3r2, IC =193), zinc finger and Broad-Complex, Tramtrack and Bric a brac (BTB) domain containing 16 (Zbtb16, IC =180), G protein-coupled receptor 56 (Gpr56, IC =179), proline and serine rich 2 (Proser2, IC =176), and midnolin (Midn, IC =175) were identified as hub genes ().

Figure 5

Interaction network of hub genes and potential miRNA-targeted regulatory networks. (A,B,C) Interaction network of the hub genes in the key modules, red nodes represent the hub genes, green nodes represent other related genes in the module. (D,E,F) Regulatory network of predicted miRNAs (Top 4 predicted miRNAs ranked by degree), red nodes represent miRNAs, green nodes represent genes.

Prediction of potential miRNA-target regulatory networks

Here, we used miRNet, a tool that integrates data from 11 different miRNA databases, to predict regulatory miRNAs of the key modules. In the turquoise module, 141 miRNAs were predicted, of which the top 4 (ranked by degrees) were mmu-mir-9-5p (degree =5), mmu-mir-466i-3p (degree =5) mmu-mir-362-3p (degree =4), and mmu-mir-329-3p (degree =4) (). In the brown module, 180 miRNAs were predicted, of which the top 4 miRNAs were mmu-mir-340- 5p (degree =12), mmu-mir-301b- 3p (degree =9), mmu-mir-124-3p (degree =7), and mmu-mir-19b-3p (degree =7) (). In the yellow module, 122 miRNAs were predicted, of which the top 4 were mmu-mir-34b-5p (degree =8), mmu-mir-149-5p (degree =6), mmu-mir-124-3p (degree =5), and mmu-mir-340-5p (degree =5) (). Surprisingly, mmu-mir-340-5p and mmu-mir-124-3p were predicted in both the brown and yellow modules, suggesting that these 2 miRNAs may be important in chronic viral myocarditis.

Discussion

It has been proposed that high-performance “omics” and bioinformatics approaches can provide unbiased information on potential genes that are important in virus pathogenesis or host defense (1). Here, we used WGCNA to identify stage-specific gene co-expression modules and uncover novel hub genes. Exploring the occurrence and development of viral myocarditis using a biological network approach may be helpful in identifying key targets. To the best of our knowledge, this is the first gene co-expression network study to analyse viral myocarditis.

In the acute infection stage, gene expression was significantly increased in the turquoise module. Moreover, the most enriched GO terms for BP were associated with antiviral and immune-inflammatory activation. It has been reported that the virus replicates significantly after infection, causing early virus-induced cardiomyocyte injury. Neutralizing antibodies appear around 1 week post-infection and play critical roles in limiting further viral replication in the heart. After the release of virus progeny into the interstitium, natural killer (NK) cells and macrophages migrate to the sites of injury and release a large amount of proinflammatory cytokines (1,7). Consistent with our analysis, this process primarily involves anti-viral and immune responses.

Based on the intramodular connectivity, we identified five hub genes (Itgb2, Tpm3, Arhgap30, Prkcd, CD300c2) in the turquoise module. Itgb2, also known as CD18, is a member of the integrin family that is critical for mounting effective immune responses (26,27). CD11/CD18 heterodimers (Mac-1) play critical roles in leukocyte trafficking, immune synapse formation, and co–stimulation (28). Itgb2 mutations causes leukocyte extravasation deficiency that impairs their migration to sites of inflammation (29). Previous studies show that Itgb2 may be involved in the regulation of B. burgdorferi-induced carditis and autoimmune carditis by recruiting macrophages (30,31). However, the role of Itgb2 in viral myocarditis needs to be further verified. Tpm3 modulates myosin-actin interaction for cardiac contraction (32). Previous studies have shown that CVB3 infection triggers a strong inflammatory response and production of autoantibodies to cardiac antigens. Significant increases in anti-myosin heavy chain, anti-actin, and anti-tropomyosin antibodies were seen in infected mice as early as day seven post-infection. Cardiac myosin-specific autoantibodies may have an immunopathogenic role in viral myocarditis (33,34). Arhgap30 is a RhoA- and Rac1-specific Rho GAP that modulates the cytoskeleton organisation and cell adhesion (35). It has been reported that it promotes p53 acetylation (36). Furthermore, an increase of acetylated p53 is related to apoptosis in CVB3-infected cardiomyocytes (37). Therefore, the increased expression of Arhgap30 may promote apoptosis in the early stage of viral myositis. Prkcd is a critical regulator of the inflammatory response in cancer, diabetes, ischaemic heart disease, sepsis, and neurodegenerative diseases (38), and its activation is associated with myocardial ischaemia-reperfusion injury (39). Furthermore, it has been previously shown that Prkcd participates in myocardial fibrosis in autoimmune myocarditis (40,41). CD300c2 (also MAIR-II, LMIR2, or CLM-4) is an Ig-like receptor expressed on macrophages, monocytes, and a subset of B cells (42,43). It has been reported that CD300c2 deficiency in a mouse sepsis model could reduce monocyte migration and suppress proinflammatory cytokine production (44,45). Furthermore, CD300c2 contributes to bleomycin-induced inflammatory responses (43). Importantly, the predicted hub genes indicate that that virus-mediated immune responses and inflammatory damage may play a major role during acute viral myocarditis.

In the chronic disease stage, gene expression was significantly increased in the brown module, while it was significantly decreased in the yellow module. In the brown module, the most enriched GO terms for BP were mainly associated with the cytoskeleton and the regulation of actin. In the yellow module, the most enriched GO terms of BP were mainly associated with phosphorylation, cellular catabolic process, vascular regulation and autophagy. Thus, the main BP changed significantly in the chronic disease stage, compared with the acute disease stage. Previous studies have shown that the chronic stage of viral myocarditis mainly includes myocardial fibrosis and cardiac remodelling, but the specific mechanism is unclear (1,2). Therefore, our analysis may offer a new direction in the exploration of viral myocarditis.

In the brown module, we identified five hub genes (Hmcn1, Pdp1, Abca8b, Rabgap1l, Magi3). The protein encoded by Hmcn1, also known as Fibulin-6, increases significantly in myocardial infarction and participates in the regulation of cardiac remodelling after myocardial infarction (46,47). Thus, Fibulin-6 might regulate fibrosis and cardiac remodelling in myocarditis. Pdp1 resides in the mitochondrial matrix and is responsible for dephosphorylation and reactivation of pyruvate dehydrogenase, thereby modulating the utilisation of carbohydrates in mammals (48). In heart failure, nicotinamide adenine dinucleotide (NADH)-dependent mitochondrial complex I activity and oxidative phosphorylation have been found to be defective, while pyruvate dehydrogenase complex activity is significantly enhanced (49). Therefore, Pdp1 may play an important role in the regulation of cardiac energy metabolism in chronic viral myocarditis. Abca8b is a membrane-associated protein belonging to the ATP-binding cassette transporters superfamily. Abca8b acts as a cholesterol transporter and may, therefore, be involved in atherosclerosis (50). Inhibition of miR-92a-3p has been is reported to regulate cardiomyocyte metabolic switching through Abca8b, thereby improving function and recovery of endothelial cells after acute myocardial infarction (51). Rabgap1l can modulate the migration of fibroblasts (52). Magi3 is associated with various diseases, including inflammatory bowel disease, thyroid disease, breast cancer, and colon cancer (53-56). Magi3 interacts with the Beta-1 and Beta-2 adrenergic receptors, selectively regulating the signal transduction of certain receptors (57,58).

In the yellow module, Slc9a3r2, Zbtb16, Gpr56, Proser2, and Midn were identified as hub genes. The protein encoded by Slc9a3r2, also known as NHERF2, plays an important role in cell signal transduction, ion transport, membrane targeting and transport of receptors (59,60). The protein encoded by Zbtb16, also known as PLZF, is a transcription factor containing 9 zinc fingers in the carboxyl terminal area. In spontaneously hypertensive rats, reduced expression of PLZF is associated with an improvement of cardiac hypertrophy and fibrosis (61). PLZF is an important angiotensin II receptor 2 binding protein and PLZF knockout mice are resistant to angiotensin-induced cardiac hypertrophy and fibrosis (62). Furthermore, PLZF may modulate congenital heart disease (63). Gpr56 is involved in a variety of BPs, including tumorigenesis, nervous system development, regulation of islet beta cells and muscle regulation of mechanical overload (64-69). In addition, GPR56 is reported to promote cardiomyocyte hypertrophy induced by angiotensin II, leading to cardiac hypertrophy (70). Proser2 is rich in serine and proline residues that maybe serve as phosphorylation sites. Although, Proser2 has been proposed as a breast cancer biomarker (71), its functions remain unclear. Midn is involved in the regulation of genes related to neurogenesis during mouse development (72), and promotes the expression of Parkin E3 ubiquitin ligase, whose down-modulation is associated with Parkinson's disease (73). Midn has been reported to interact with glucokinase through an N-terminal ubiquitin-like domain, resulting in a significant reduction in glucokinase activity and glucose-induced insulin secretion (74). In conclusion, the hub genes in the brown and yellow modules indicate that myocardial fibrosis, abnormal energy metabolism, cellular function, and cytoskeleton changes may be involved in chronic viral myocarditis.

Previous studies have shown that miRNAs also play key roles in viral myocarditis (75-79). Therefore, we predicted potential miRNA-targeted regulatory networks in the three key modules. Of these, we found that mmu-mir-340-5p and mmu-mir-124-3p may play an important role in chronic myocarditis.

In our study, the hub genes were identified according to intramodular connectivity. Generally speaking, genes with higher intramodular connectivity are often at the core position in the co-expression module, and play an important role in the development of the disease. However, some “hub genes” may show high intramodular connectivity due to the regulation of their upstream hub genes. These regulated “hub genes” may be significantly associated with the disease, but they are not the primary drivers of disease. Therefore, our results need further verification. Furthermore, rare interactions might have been obscured due to the small sample size available for analysis.

Conclusions

Our results reveal the main biological changes in different stages of viral myocarditis, and predict the important hub genes and miRNAs in both the acute and chronic disease stages. These findings provide novel insights into the pathogenesis of viral myocarditis, which might be helpful for developing new therapeutic targets in the future.

The article’s supplementary files as