The GSE19303 gene expression microarray and clinical features dataset was obtained from Gene Expression Omnibus (http://ncbi.nlm.nih.gov/geo/). The dataset contained a total of 81 endomyocardial biopsy samples, 40 baseline biopsies from patients with DCM, 33 of 40 patients had received immunotherapy (immunoadsorption with subsequent immunoglobulin substitution (IA/IgG) treatment), 33 follow-up biopsies of DCM patients collected 6 months after the treatment, and eight biopsies from individuals without DCM (Ameling et al., 2016). Among them, 40 baseline and 33 follow-up endomyocardial biopsy samples from DCM patients were utilized in our study and were divided into groups of different sexes. The male patient group contained 28 samples, and the female patient group included 12 samples. The platform for the gene expression profiles was GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array (Affymetrix; Thermo Fisher Scientific, Inc., Waltham, MA, United States).

We retrieved the expression matrix from the GEO database and preprocessed it by using the robust multiarray analysis (RMA) method (http://www.bioconductor.org/). After log2 transformation and quantile normalization of the expression data, we annotated the converted probe ID for each gene to a gene symbol utilizing hgu133plus2. db, org. Hs.eg.db, and the annotate package in Bioconductor (http://www.bioconductor.org/). If a gene’s symbol corresponded with the multiple probe IDs, the expression level of that gene was represented by the mean of the probes.

T-tests or paired t-tests were used to test for differences, and outlier samples were assessed using 1.5 times the interquartile range of the differences, and Shapiro–Wilk normality test was used for normality tests. A p-value less than 0.05 is considered statistically significant.

IOBR (Immuno-Oncology Biological Research) is a computational tool for immuno-tumor biology research (Zeng et al., 2021). Here, based on our expression profiles, we use the IOBR package in R to analyze the immune cell infiltration abundance of GSE19303 datasets. The xCell method was selected to calculate infiltration abundance of 64 kinds of immune cells, stem cells, and stromal cells in each sample (Aran et al., 2017). The DEGs of the ceRNA network were divided into the high-expression and low-expression groups by the median. The relationship between DEG expression and the fractions of immune cells was investigated by Wilcoxon test. The results were visualized using the ggplot packages in R software.

The Linear Models for Microarray Analysis (Limma) package in R software was applied to identify the differentially expressed genes (DEGs) in the male DCM samples compared with the female DCM ones, based on Student’s t-test. Adjusted p-values were calculated using the Benjamini–Hochberg method. The significant DEGs were selected with a threshold of p-value < 0.05 and fold change>1.5. We obtained the volcano plot utilizing the pheatmap package in R.

Moreover, we excluded the genes in the Y chromosomes and then compared the male and female DCM samples again using the Limma package. p-value < 0.05 and fold change>1.5 were selected to be the cutoff criteria of the significant DEGs. We obtained another volcano and heatmap plot. In order to plot the differentially expressed lncRNAs (DE-lncRNAs), we used the pheatmap package in R to construct a heatmap.

In order to explore the potential functions and pathways that may be altered by the DEGs, we applied the clusterProfiler package in R to perform the functional and pathway enrichment analyses of the identified DEGs. The Gene Ontology (GO; http://www.geneontology.org/) database was used to determine the biological processes (BPs) that the DEGs may be involved in. In addition, according to the modified Fisher’s exact test, the Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.jp/kegg/pathway.html) database was used for pathway enrichment analysis of the identified DEGs. The selection criteria for the significant GO terms and pathways were p < 0.05, and the number of enriched genes were (also called count) > 2. The plots were performed by the ggplot2 package in R.

GSEA software (version 3.0) and c2. all.v7.4. symbols.gmt subcollection were obtained from the GSEA website (http://software.broadinstitute.org/gsea/index.jsp) (Subramanian et al., 2005). We divided the samples into two groups by sex. The minimum gene set was 5 and the maximum gene set was 5,000, with 1,000 resampling. p-value < 0.05 (as needed) or FDR <0.25 (as needed) were considered statistically significant.

Potential interactions between DE-lncRNAs and DE-miRNAs and between DE-miRNAs and DEGs were predicted using DIANA (https://diana.e-ce.uth.gr/lncbasev3) (Karagkouni et al., 2020) and ENCORI databases (https://starbase.sysu.edu.cn/index.php) (Li et al., 2014), respectively. Only the lncRNA–miRNA and miRNA-DEG interactors, that had an opposite expression trend, were used to construct the ceRNA network. The immune gene list was obtained from the Immunology Database and Analysis Portal (IMMPORT) database (http://www.immport.org/) (Bhattacharya et al., 2014). The Venny online tool was used to analyze the overlapping genes (http://jvenn.toulouse.inra.fr/app/example.html) (Bardou et al., 2014).

Coexpression patterns in 29 immune-related DEGs were analyzed using Pearson’s correlation coefficient, and the results were visualized using the heatmap packages in R software. Cytoscape software (version 3.8.2, https://cytoscape.org) was used to develop the ceRNA network. For gene set functional enrichment analysis, we used the GO and KEGG annotations of genes in the R package org. Hs.eg.db (version 3.1.0) as the background to map the genes to the background set using the R package clusterProfiler (version 3.14.3) to perform enrichment analysis to obtain gene set enrichment results. p-value < 0.05 was considered statistically significant.

A protein–protein interaction (PPI) network, comprising 50 ceRNA network coexpression proteins, was constructed by GeneMANIA (http://genemania.org/) (Warde-Farley et al., 2010). These nodes represent genes that are closely related to the ceRNA network in terms of physical interactions, shared protein domains, predictions, colocalization, pathway, coexpression, and genetic interactions. We use NetworkAnalyst (version 3.0, https://www.networkanalyst.ca/) to carry out the heart (left ventricle)-specific PPI, TF-miRNA interactions, and protein–chemical interaction analysis on the ceRNA network coexpression module (Zhou et al., 2019). In these networks, the nodes represent individual genes/proteins/chemicals, while the edges which connect the nodes correspond to a known, curated interaction between a given pair of nodes.

All the statistical analyses were performed using R (version 3.6.4) or SPSS (version 19.0), and a p-value less than 0.05 is considered statistically significant. The plots were performed by R, Cytoscape software (version 3.8.2, https://cytoscape.org), or SangerBox tools (version 3.0, http://www.sangerbox.com/tool).

The RNA array data for a total of 40 DCM patients were acquired from the GES19303 dataset. The detailed baseline clinical features are listed in Supplementary Table S1. Among the 40 participants, 28 were male and 12 were female. DCM patient gender (male vs. female) was significantly correlated with age (52.2 ± 9.21 years vs. 45.42 ± 8.08 years, p = 0.032) and Epstein–Barr virus (EBV) infection (0/28 vs. 2/12, p = 0.027). However, gender was not significantly correlated with other clinical features such as LVEF, LVIDD, BMI, inflammation index, total virus infection, PVB19 infection, HHV6 infection, HSV1 infection, and IgG treatment.

Moreover, we investigated the association between IA/IgG treatment and clinical outcome features in DCM patients. In female patients, immunotherapy was significantly correlated with LVEF (32.00 ± 7.44% vs. 40.00 ± 7.93%, p = 0.002) and LVIDD (67.20 ± 9.40% vs. 64.40 ± 9.79%, p = 0.047). For male patients, immunotherapy was more significantly correlated with LVEF (34.70 ± 5.43% vs. 43.22 ± 8.88%, p = 3.74E-04) and LVIDD (70.74 ± 6.72% vs. 65.70 ± 10.30%, p = 1.417E-04) (Figures 1A,B). However, IgG treatment was not significantly correlated with the inflammation index (Table 1). These results suggest that IA/IgG immunotherapy could significantly improve the outcome of DCM especially in male patients.

To gain insight into potential target immune cells of DCM IA/IgG immunotherapy, we estimated the composition of the microenvironment in baseline DCM patients by using the xCell algorithm. Our result showed that the composition of the microenvironment of DCM was complex (Figure 2A). The top five abundant cell types were multipotent progenitors (MPPs), mesenchymal stem cells (MSCs), natural killer T cells (NKTs), immature dendritic cells (iDCs), and microvascular endothelial cells (mv endothelial cells) (Figure 2B). Moreover, when compared with the female DCM patients, male DCM patients had a significant positive correlation with abundance of B cells, memory B cells, effect memory CD8⁺ cells, and NK cells (Figure 2C). These results suggest that sex-differentiated microenvironments may contribute to differences in immunotherapy efficacy.

A total of 1,138 DEGs were finally screened from the comparison of male DCM samples with female DCM samples, including 579 upregulated and 556 downregulated DEGs (Figure 3A). Considering that genes in the Y chromosome is few in number and misleading, we removed the genes in the Y chromosome in the subsequent analyses. After excluding the genes in the Y chromosomes, we obtained 1,071 DEGs, of which 542 were upregulated and 529 were downregulated in the male group (Figure 3B). The top 50 significant positive and negative sex-related DEGs correlated with DCM are shown in the heatmap (Figure 3C). Interestingly, there were 68 differentially expressed lncRNAs (DE-lncRNAs) and five differentially expressed miRNAs (DE-miRNAs) in DEGs, and the top 50 significant DE-lncRNAs are shown in the heatmap (Figure 3D).

The upregulated and downregulated DEGs were processed separately for the GO and KEGG pathway analyses. The significantly enriched biological processes (BPs) were negative regulation of cellular process, negative regulation of programmed cell death, cardiovascular system development, ncRNA metabolic process, and RNA modification (Figure 4A). The significantly enriched cellular components (CCs) were the cytoskeleton, endoplasmic reticulum part, extracellular matrix, complex of collagen trimers, and RNA polymerase I transcription factor complex (Figure 4B). The significantly enriched molecular functions (MFs) were molecular function regulator, signaling receptor binding, extracellular matrix structural constituent, core promoter binding, and antiporter activity (Figure 4C). The significantly enriched KEGG were focal adhesion, tight junction, amebiasis, ECM–receptor interaction, and amino-acid metabolism (Figure 4D). These results suggest that there is a widespread impact of sex difference on the signaling transduction, transcriptional regulation, and metabolism.

To confirm the abovementioned results, we performed GSEA on the sex-related DEGs to identify Gene Ontology and signaling pathways that were differentially activated in DCM. The top 20 significant positive and negative enrichment pathways are shown in Figure 5A. GSEA showed that significant positively enriched pathways in male are aminoacyl tRNA biosynthesis, TP53 targets apoptotic, protein repair, RNA polymerase I transcription initiation, and TRAIL signaling. (Figure 5B). Gene sets related to collagens, TGF-beta receptor signaling, epithelial-to-mesenchymal transition, degradation of the extracellular matrix, and IL6 signaling showed enrichment in the female DCM patients (Figure 5C).

Among the 68 DE-lncRNAs and five DE-miRNAs, five DE-lncRNAs (ARHGEF7-IT1, LINC00632, LINC02135, TEX36-AS1, and X-inactive specific transcript (XIST)) were predicted to interact with four DE-miRNAs (miR-1-3p, miR-17-5p, miR-22-3p, and miR-146a-5p) by the DIANA database. The four DE-miRNAs were predicted to interact with a total of 7,747 genes by the ENCORI database (Figure 6A). To further construct the sex difference immune-related ceRNA network, the Venny method was used to analyze the intersection between DEGs, miR-related genes, and immune genes (Figure 6B). The coexpression plot among these 29 genes is presented in Figure 6C. These DE-lncRNA–miRNA–immune gene interaction pairs were integrated to construct the ceRNA network (Figure 6D). Furthermore, we extracted a central module from the ceRNA network, including two lncRNAs (XIST and LINC00632), three miRNAs (miR-1-3p, miR-17-5p, and miR-22-3p), and six mRNAs (Casitas B-Lineage Lymphoma Proto-Oncogene (CBL), C-X-C Motif Chemokine Ligand 12 (CXCL12), Estrogen Receptor 1 (ESR1), Insulin-like Growth Factor 1 Receptor (IGF1R), Interleukin 6 Cytokine Family Signal Transducer (IL6ST), and Stanniocalcin 1 (STC1)) (Figure 6E).

GO and KEGG enrichment analyses were conducted to investigate the functions of the ceRNA network. GO_BP analysis showed that ceRNA network was significantly enriched in programmed cell death, immune system process, and heart development (Figure 7A). GO_CC analysis showed that the ceRNA network was significantly enriched in the cell surface, endoplasmic reticulum, and collagen-containing extracellular matrix (Figure 7B). GO_MF analysis showed that the ceRNA network was significantly enriched in molecular function regulator, cytokine activity, and RNA polymerase II transcription factor binding (Figure 7C). The KEGG pathway enrichment analysis showed that ceRNA network was significantly enriched in cytokine–cytokine receptor interaction, Ras signaling pathway, and natural killer cell–mediated cytotoxicity (Figure 7D).

The protein–protein interaction (PPI) network was assembled based on the ceRNA network in the DCM cohort by GeneMANIA. The analysis showed that the 50 most significantly coexpressed genes play roles in peptidyl-tyrosine phosphorylation, leukocyte chemotaxis, and the vascular process in the circulatory system (Figure 8A). Next, the ceRNA network was assembled based on the heart (left ventricle)-specific data collected from the DifferentialNet database (Basha et al., 2018) by NetworkAnalyst (Figure 8B). The top five hub proteins were CBL, FYN Proto-Oncogene (FYN), Kinase Insert Domain Receptor (KDR), ESR1, and Histone Deacetylase 1 (HDAC1) (Supplementary Table S2). Furthermore, a graph of TF-miRNA coregulatory interactions of the ceRNA network was constructed based on the RegNetwork database (Liu et al., 2015) (Figure 8C). From this, the top five TFs identified were Myelocytomatosis Oncogene (MYC), Nuclear Factor Kappa B Subunit 1 (NFKB1), Specificity Protein 1 (SP1), MYC Associated Factor X (MAX), and Upstream Transcription Factor 1 (USF1) (Supplementary Table S3).

To gain insight into the potential target drugs based on our established ceRNA network, we examined protein–chemical interactions from the Comparative Toxicogenomics Database (CTD) (Davis et al., 2021) (Figure 8D). Excluding hazardous chemicals, the top five drugs were 4-fluorobenzoyl-TN-14003, polyphenon E, 3-(4-methylphenylsulfonyl)-2-propenenitrile, deoxyglucose, and arachidonic acid (Supplementary Table S4). 4-fluoro benzoyl-TN-14003(BKT140, motixafortide) is a high-affinity CXCR4 antagonist, which can inhibit the migration of prostate cancer cells (Peng and Kopecek, 2014). BKT140 has been tested in stimulating megakaryopoiesis and platelet production (Abraham et al., 2013). Polyphenon E is a green tea polyphenol preparation which possesses potent antioxidative and anti-inflammatory properties (Bornhoeft et al., 2012). 3-(4-methylphenylsulfonyl)-2-propenenitrile (BAY11-7,082) is an NF-kappaB inhibitor, which can protect the myocardial infarction heart from cardiac dysfunction in the mouse model (Martinez-Martinez et al., 2017). Moreover, BAY11-7,082 significantly reduced the TNF and IL-6 protein expression in atherogenesis (Vallejo et al., 2018). Deoxyglucose (2-DG) is a glucose molecule which cannot undergo further glycolysis. 2-DG can antagonize DOX-induced cardiomyocyte death, which is mediated through multiple mechanisms, including the preservation of ATP content, the activation of AMPK, and the inhibition of autophagy (Chen et al., 2011). Arachidonic acid (AA) is an essential fatty acid, which can be found in fish and certain plant oils. Higher in vivo circulating and tissue levels of AA were associated with lower risk of major cardiovascular events (Marklund et al., 2019). In addition, sex differences in the AA levels could be an important underlying mechanism for different effects of sex hormones and cardiovascular disease differences between males and females (Gerges and El-Kadi, 2021). Based on the previous research, these drugs show promising potential as novel therapies against DCM via the immune-related ceRNA network. However, further evaluation is still needed.

For further analysis, the DEGs of the ceRNA network were divided to the high-expression and low-expression groups by median. The correlation between immune cell abundance and DEGs was analyzed by the Wilcoxon test. CBL expression had a significant positive correlation with macrophages (p = 0.03), neutrophils (p = 0.008), and fibroblasts (p = 0.02). Moreover, CBL expression had a significant negative correlation with CD8+Tem (p = 0.0096) and conventional DC (cDC, p = 0.0056) (Figure 9A). CXCL12 expression had a significant positive correlation with M1 macrophages (p = 0.03), activated DCs (aDC, p = 0.0035), MSCs (p = 0.02), lymphatic endothelial cells (ly endothelial cells, p = 0.04), and mv endothelial cells (p = 0.0032). Furthermore, CXCL12 expression had a significant negative correlation with B-cells (p = 0.01), Treg (p = 0.03), M2 macrophages (p = 0.04), hematopoietic stem cells (HSC, p = 0.04), and common lymphoid progenitor (CLP, p = 0.0067) (Figure 9B). IL6ST expression had a significant positive correlation with Tregs (p = 0.0014) and myocytes (p = 0.03). Moreover, IL6ST expression had a significant negative correlation with iDC (p = 0.02), ly endothelial cells (p = 0.0044), and mv endothelial cells (p = 0.01) (Figure 9C). However, ESR1, IGF1R, XIST, and LINC00632 expression had no significant correlation with most immune cells in DCM.