Identification of immune feature genes and intercellular profiles in diabetic cardiomyopathy

Abstract

BACKGROUND

Diabetic cardiomyopathy (DCM) is a multifaceted cardiovascular disorder in which immune dysregulation plays a pivotal role. The immunological molecular mechanisms underlying DCM are poorly understood.

AIM

To examine the immunological molecular mechanisms of DCM and construct diagnostic and prognostic models of DCM based on immune feature genes (IFGs).

METHODS

Weighted gene co-expression network analysis along with machine learning methods were employed to pinpoint IFGs within bulk RNA sequencing (RNA-seq) datasets. Single-sample gene set enrichment analysis (ssGSEA) facilitated the analysis of immune cell infiltration. Diagnostic and prognostic models for these IFGs were developed and assessed in a validation cohort. Gene expression in the DCM cell model was confirmed through real time-quantitative polymerase chain reaction and western blotting techniques. Additionally, single-cell RNA-seq data provided deeper insights into cellular profiles and interactions.

RESULTS

The overlap between 69 differentially expressed genes in the DCM-associated module and 2483 immune genes yielded 7 differentially expressed immune-related genes. Four IFGs showed good diagnostic and prognostic values in the validation cohort: Proenkephalin (Penk) and retinol binding protein 7 (Rbp7), which were highly expressed, and glucagon receptor and inhibin subunit alpha, which were expressed at low levels in DCM patients (all area under the curves > 0.9). SsGSEA revealed that IFG-related immune cell infiltration primarily involved type 2 T helper cells. High expression of Penk (P < 0.0001) and Rbp7 (P = 0.001) was detected in cardiomyocytes and interstitial cells and further confirmed in a DCM cell model in vitro. Intercellular events and communication analysis revealed abnormal cellular phenotype transformation and signaling communication in DCM, especially between mesenchymal cells and macrophages.

CONCLUSION

The present study identified Penk and Rbp7 as potential DCM biomarkers, and aberrant mesenchymal-immune cell phenotype communication may be an important aspect of DCM pathogenesis.

Key Words: Diabetic cardiomyopathy; Immune feature genes; Proenkephalin; Retinol binding protein 7; Immune cell infiltration; Intercellular communication

Core Tip: In this study, we utilized bulk RNA sequencing (RNA-seq) data and machine learning techniques to identify and validate four immune feature genes associated with diabetic cardiomyopathy (DCM). Notably, retinol binding protein 7 and proenkephalin showed significantly elevated expression in cardiomyocytes and interstitial cells in DCM, as confirmed by single-cell RNA-seq and molecular experiments, highlighting their robust diagnostic potential. Furthermore, single-cell RNA-seq data revealed abnormal cellular phenotype transformations and communications in DCM, particularly between fibroblasts and macrophages.

INTRODUCTION

Cardiovascular disease continues to be the leading cause of mortality globally, with diabetes mellitus representing a significant risk factor[1,2]. Diabetic cardiomyopathy (DCM), a severe chronic complication of diabetes that profoundly affects the cardiac structure and function, culminates in heart failure without the presence of common causes like coronary artery disease, hypertension, and valvular heart disease[3]. DCM manifests as alterations in heart size, shape, and function, alongside an increase in fibrous tissue (cardiac fibrosis)[4,5]. Jia et al[6] reported the occurrence of heart failure in individuals with diabetes ranges between 19% and 26% and found that nearly one-third of asymptomatic patients had evidence of subclinical DCM.

The development and progression of DCM involve complex and multifactorial mechanisms, including hyperglycemia, insulin resistance, inflammation, immune activation, oxidative stress, and mitochondrial dysfunction[6-9]. The heart consists of cardiomyocytes as well as various non-muscle cells, including fibroblasts, blood vessel cells, neurons, and immune cells[10]. A comprehensive understanding of the intricate adaptive processes triggered by cardiac injury remains a considerable challenge. However, inflammation and immune cell signaling are key components of this process. Recent research has highlighted the significant role of inflammation and immune cell signaling in the pathology of DCM. Accumulating evidence suggests the involvement of cardiac myocytes, immune cells, and endothelial cells in the development of cardiac fibrosis[11-13]. Cardiomyocyte death triggers immune pathways that, initiate inflammatory responses[14]. Perturbations in the innate and adaptive immune systems play roles in the occurrence and progression of DCM[13,15], which is evidenced by imbalances in immune cell populations, enhanced generation of pro-inflammatory molecules, and elevated oxidative stress[10,16,17].

Many questions are unresolved in the molecular study of DCM. Despite advances in our understanding of immune dysregulation in DCM, many molecular and cellular characteristics of this process are not clear. Investigating these aspects is crucial for elucidating the immunopathology of DCM and identifying potential therapeutic targets[13]. According to a review by Lorenzo-Almorós et al[18], the prevalence of DCM varies greatly, influenced by the specific population studied and the diagnostic criteria applied. Therefore, the present study investigated immune dysregulation in DCM and identified key molecular characteristics associated with immune-mediated pathology.

Machine learning algorithms have been used to analyze large-scale omics datasets, such as gene expression profiles, proteomics data, and metabolic profiles, to identify molecular signatures associated with DCM. For example, Hathaway et al[19] used machine learning techniques to identify key genes and pathways involved in the pathogenesis of DCM, which provided insights into potential therapeutic targets. Single-cell sequencing has facilitated the identification and characterization of various heart cell types, such as cardiomyocytes, fibroblasts, immune cells, and endothelial cells. This technology allows for the examination of gene expression profiles at the single-cell level, helping researchers uncover cell-specific transcriptional changes and dysregulated pathways in DCM. We used the above multidimensional approach to analyze immune gene expression, profile immune cell populations, and assess intercellular communication patterns in DCM to advance our understanding of the intricate molecular and immunological mechanisms that drive the progression of DCM.

MATERIALS AND METHODS

Data sources and processing

Independent bulk RNA sequencing (RNA-seq) (GSE5606, GSE6880, GSE4745, and GSE197999) and single-cell RNA-seq (scRNA-seq) datasets (GSE213337) were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). GSE5606 and GSE68803 were combined into a comprehensive dataset containing 20 samples as discovery cohorts, and GSE4745 and GSE197999 were merged as validation cohorts. Gene expression was adjusted for each sample using the "ComBat" method in the R (version 4.2.1) package "sva" to eliminate differences due to batch effects. For the scRNA-seq dataset, GSE213337 contained DCM and a healthy control sample, and low-quality cells were filtered out (cells with gene counts between 200 and 2500 were retained). Two samples were then combined into one integrated dataset and subjected to batch effect removal and dimensionality reduction using the "Harmony" algorithm in the “Seurat” package (version 4.3.0) for subsequent analysis. ProteomicsDB (https://www.proteomicsdb.org/) contains thousands of samples from a variety of biological sources, such as tissues, body fluids, or cell lines, and it was used to explore tissue and cellular expression patterns of target genes.

Weighted gene co-expression network analysis and identification of disease-related modules

Weighted gene co-expression network analysis (WGCNA) was performed using the "WGCNA" package (version 1.71). Initially, the scale-free network properties were evaluated to determine a soft threshold. The correlation matrix was then converted into an adjacency matrix using this power value, and this matrix provided the foundation for building a co-expression network. A topological overlap matrix (TOM) was created with a minModuleSize value of 30 and a mergeCutHeight value of 0.2 to quantify similarity. Dendrograms with hierarchical clustering were created to display genes as determined by TOM (1-TOM). Module eigengenes (MEs) were calculated as the first principal component to summarize the gene expression profile within each module. Correlations between MEs and clinical features revealed modules associated with the disease of interest. Strong relationships between the modules and the genes of interest were shown by the correlation between gene significance (GS) and module membership (MM). Subsequently, normalization of the gene expression microarray data and differentially expressed genes (DEGs) were screened through the "LIMMA" package. DEGs were defined as P < 0.05 and |LogFC| > 1. The "ggVolcano" package was utilized to create volcano graphs visualizing the DEGs.

Identification of immune feature genes using machine learning algorithms

A list containing 2483 immune-related genes (IRGs) was acquired from the ImmPort database (https://www.immport.org/). Intersections were taken between the DEGs in the disease-related modules identified by WGCNA and the list of IRGs to identify differential IRGs (DIRGs). Least absolute shrinkage and selection operator (LASSO) and support vector machine recursive feature elimination (SVM-RFE), were used to select the most disease-related identification of immune feature genes (IFGs) and determine their intersections. LASSO and SVM-RFE effectively reduce the number of features and select the most relevant ones for predictive modeling or data analysis. Feature genes with non-zero coefficients were selected by using lambda.min to identify optimal λ values, and the LASSO regression model was constructed using the R package "glmnet". The R package "caret", which includes the SVM-RFE algorithm, is a powerful tool for feature selection. The "rfe" function in this package allows us to construct RFE models based on the "svmRadial" method to select feature genes.

Construction of the gene expression regulatory network

The NetworkAnalyst program (https://www.networkanalyst.ca/) was utilized to generate transcription factor-gene and miRNA-gene regulatory networks for IFGs. The transcription factor network was constructed with JASPAR as the retrieval target dataset, and the miRNA regulatory network was constructed using miRTarBase 8.0. The nodes and edges contributed by both networks were unionized into a single network, which was modified by Cytoscape software (version 3.10.1).

Nomogram model construction, receiver operating characteristic analysis, and decision curve analysis

Using the "rms" R package, a prognostic nomogram was created to assess the expected risk of disease based on the feature gene expression levels. Harrell’s concordance index (C-index) was performed to assess the performance of the nomogram. The calibration curve compares the actual probability of occurrence with the predicted probability, which was constructed by repeating bootstrap self-sampling 100 times (B = 100) in the R package "PredictABEL" to obtain reliable estimates of calibration performance for validating the nomogram model. The R package "pROC" was applied to assess the diagnostic value of IFGs in DCM. This approach facilitated the determination of the optimal gene expression cut-off value, the generation of receiver operating characteristic (ROC) curves, and the calculation of the area under the curve (AUC). A multivariable logistic regression model was employed to construct a combined diagnostic ROC curve. Decision curve analysis (DCA) attempts to construct models that predict maximum net benefits. Both simple models with single feature genes as predictors and complex models with multiple genes were constructed using the "decision_curve" function in the "rmda" package.

Single-sample gene set enrichment analysis

Single-sample gene set enrichment analysis (ssGSEA) assesses immune cell infiltration by comparing a gene expression to a reference immune cell gene set and was performed using the "GSVA" package (version 1.46.0), which outputs the ssGSEA score for each sample[20]. The correlation matrix was computed using the function rcorr in the package "Hmisc" for the feature genes and immune cell populations.

Human ventricular cardiomyocyte (AC16 cell) culture and high glucose exposure

AC16 cells were purchased from Procell (Procell, China) and cultured under standard conditions in Dulbecco's modified Eagle medium (DMEM) containing 10% fetal bovine serum (Procell, China) at 37 °C in an atmosphere of 5% CO₂ and 95% air. Cell passages ranging between 2 and 20 were utilized in the study. Before experimentation, the cells were cultured in DMEM with 5.5 mmol/L glucose for a minimum of four passages. The cells were exposed to high glucose (HG) conditions via treatment with 33 mmol/L D-glucose and 0.3 mmol/L palmitate (Sigma-Aldrich, New Zealand) for 24 to 48 hours to simulate an in vitro diabetic environment. Cells were plated at 3 × 10⁵ cells per 60-mm dish for RNA and protein extraction. Samples were harvested 24 and 48 hours following HG conditions.

Real time-quantitative polymerase chain reaction

As directed by the manufacturer, total RNA was isolated using the TRNzol Universal Kit (DP424, TIANGEN, Beijing, China). The HiScript III All-in-One RT SuperMix Reagent (R333-01, Vazyme, Nanjing, China) was used to reverse-transcribe one microgram of RNA template following the manufacturer's instructions. Using certain primers and the ChamQ Universal SYBR quantitative polymerase chain reaction Master Mix (Q711-02, Vazyme, Nanjing, China), real-time polymerase chain reaction (PCR) was carried out. Using an Applied Biosystems 7500 real-time PCR machine, the samples were cycled under the following conditions: 42 cycles of denaturation at 95 °C for 20 seconds and extension at 65 °C for 1.5 minutes followed by an initial denaturation step at 95 °C for 3 minutes. Using GAPDH as a reference gene, the expression of target genes was calculated through the 2^-ΔCq method. The primer sequences below were specifically used: Proenkephalin (Penk): CGGTTCCTGACACTTTGCACT (forward primer) and CACATTCCATTACGCAAGCCA (reverse primer); retinol binding protein 7 (Rbp7): CTCAGCGGTACTTGGACCC (forward primer) and CGAGTGGCAAAGTCAATACCT (reverse primer); and GAPDH: GGAGCGAGATCCCTCCAAAAT (forward primer) and GGCTGTTGTCATACTTCTCATGG (reverse primer).

Western blotting assay

Lysed AC16 cell samples were placed in RIPA buffer with a cocktail of 4% protease inhibitors (HY-K0021, Med-ChemExpress, Shanghai, China). The protein concentration was estimated using a BCA Protein Assay Kit (PA115, TIANGEN, Beijing, China), and an equal amount of protein (30 µg) was loaded onto a 10%-15% SDS polyacrylamide gel and separated for 1.2 hours at 90 V. The membranes were blocked with 5% nonfat milk and then left to overnight incubation with primary antibodies. Using Pro-Light HRP chemiluminescent reagent (PA112, TIANGEN, Beijing, China), PVDF membranes were observed with a ChemiDoc XR-S system (Bio-Rad, United States) after being treated with the relevant horseradish peroxidase-conjugated secondary antibodies. The protein band density was determined using ImageJ software. A 1:100 dilution of Rbp7 primary antibody was purchased from Thermo Fischer Scientific (MA5-24514, United States), and a 1:10000 dilution of GAPDH (ab181602, Abcam, Shanghai, China) was applied as an internal reference for expression value normalization. A 1:200 dilution of the Penk primary antibody was purchased from ABclonal (A6302, Wuhan, China), and a 1:2000 dilution of β-tubulin (10068-1-AP, Wuhan, China) was used as an internal reference.

scRNA-seq computational pipelines and analysis

In the integrated dataset, linear principal component analysis targeted variable genes, using the top 20 principal components for UMAP-based nonlinear dimensionality reduction. These principal components were employed to build a shared nearest neighbor (SNN) graph based on the "FindNeighbors" function. Clusters were then formed using the "FindClusters" function applied to the SNN graph. Identification of cell types was achieved by detecting DEGs in each cluster with the "FindAllMarkers" function of Seurat. Verification of cluster labels was manually conducted using established databases, including CellMarker2.0[21] and PanglaoDB[22]. Gene expression data for each identified cell type were retrieved using the "Fetchdata" function, facilitating further statistical analysis and visual representation.

The cell-cell communication analysis relied on two R packages, "celltalker" (version 0.0.7.900) and "cellChat''[23] (version 1.6.1). The "celltalker" package provided access to intercellular communication networks via the Ramilowski pairs dataset. Circle plots were generated to display the top 3 significant interactions for each cell type, which were determined based on specified criteria (false discovery rate < 0.05). The “cellChat” package was further used for the analysis of communication patterns by setting the ligand-receptor set to CellChatDB.mouse. The function "computeCommunProb" inferred the cell-cell interactions, and "computeCommunProbPathway" further inferred the communication at the signaling pathway level. The aggregated cell-cell communication network was constructed using the "aggregateNet" function. Network centrality scores were computed using the "netAnalysis_computeCentrality" function, and the number of patterns of interaction was determined using the "selectK" function, which is presented as river plots by "net-Analysis_river".

RESULTS

DCM-related gene modules

To identify DCM-related gene modules, we performed WGCNA on an integrated dataset of two independent bulk RNA-seq datasets (GSE5606 and GSE6880) (Table 1). The soft thresholding method emphasizes strong correlations between genes and penalizes weak correlations, which helps construct scale-free networks. We selected an optimal power value of 5 for the construction of DCM scale-free networks (Figure 1A). Gene clustering resulted in the identification of 17 distinct modules, represented in different colors, where each module is comprised of genes that strongly correlated with each other but not with genes in other modules (Figure 1B). After excluding the gray module (genes in the module typically do not align with well-defined biological functional categories), module correlation analysis showed that the green module had the most significant strong positive association with DCM (Figure 1C). Correlation analysis between GS and green MM (581 genes) showed a correlation coefficient of 0.96 (P < 0.0001) (Figure 1D). These findings indicate that the genes in the green module are crucially involved in the biological processes associated with DCM.

Figure 1 Weighted gene co-expression network analysis identified diabetic cardiomyopathy-related modules in the discovery cohort. A: Soft thresholding power analysis enabled the determination of the scale-free fit index for network topology; B: Cluster dendrogram and module assignment for modules from weighted gene co-expression network analysis in diabetic cardiomyopathy (DCM); C: Correlations of module eigengenes and DCM; D: Correlations between gene significance and green module membership. DCM: Diabetic cardiomyopathy.

Table 1 Detailed information on the datasets used in this study.

Accession number	Method	Origin	Samples		Aim
			DCM	CTL
GSE5606	Bulk RNA-seq	Rattus norvegicus	7	7	Merged as discovery cohorts
GSE6880	Bulk RNA-seq	Rattus norvegicus	3	3
GSE4745	Bulk RNA-seq	Rattus norvegicus	12	12	Merged as validation cohorts
GSE197999	Bulk RNA-seq	Rattus norvegicus	5	5
GSE213337	scRNA-seq	Mus musculus	1	1	Single cell analysis

RNA-seq: RNA sequencing; DCM: Diabetic cardiomyopathy; CTL: Control.

Penk, glucagon receptor, Rbp7, and inhibin subunit alpha as IFGs in DCM

To understand the immunological molecular pathology of DCM, we identified IFGs associated with this disease. First, we extracted the gene profile of the most correlated green module and calculated the DEGs (Figure 2A). We screened a total of 69 DEGs, including 37 up-regulated genes (Figure 2B) and 32 down-regulated genes in the green module. The 69 DEGs were intersected with the list of 2483 IRGs to obtain 7 DIRGs (Figure 3A). To further identify the most representative IFGs among the 7 DIRGs, machine learning models were constructed based on the LASSO and SVM-RFE algorithms. The LASSO regression model identified five genes with non-zero coefficients, which suggests that these genes [Penk, glucagon receptor (Gcgr), Rbp7, inhibin subunit alpha (Inha), and Angptl4] are strongly associated with DCM (Table 2). SVM-RFE is a feature selection technique that uses SVM to repeatedly construct a model and remove features with the smallest ranking criterion to identify the most critical features for disease classification or prediction. The SVM-RFE model optimized for the best predictive performance highlighted four genes (Penk, Inha, Gcgr, and Rbp7) as having the best prediction performance, where the accuracy and kappa coefficient of the model were 1.0 (Figure 3B and Table 2). Therefore, these four IFGs may play essential roles in the immunopathological characteristics of DCM. Among these characterized IFGs, the expression of Penk and Rbp7 was up-regulated, and Gcgr and Inha was down-regulated (Figure 3C).

Figure 2 Identification of differentially expressed genes in the green module. A: Volcano plot of differentially expressed genes (DEGs) in the green module; B: Heatmap of the 37 upregulated DEGs in the green module.

Figure 3 Identification of immune feature genes. A: The differentially expressed genes in the green module were intersected with the list of immune-related genes from the ImmPort database to obtain 7 differential immune-related genes; B: Machine learning algorithm based on least absolute shrinkage and selection operator (LASSO) and support vector machine recursive feature elimination (SVM-RFE) to select immune feature genes (IFGs) of diabetic cardiomyopathy (DCM). LASSO identified and selected 5 genes, whereas SVM-RFE identified and selected 4 genes; C: Box plot comparing the expression levels of the four common IFGs identified by LASSO and SVM-RFE between DCM patients and control patients; D: Gene expression regulatory networks. Red indicates target genes, green indicates transcription factors, and blue indicates miRNAs. LASSO: Least absolute shrinkage and selection operator; SVM-RFE: Support vector machine recursive feature elimination; DCM: Diabetic cardiomyopathy; CTL: Control; Penk: Proenkephalin; Gcgr: Glucagon receptor; Rbp7: Retinol binding protein 7; Inha: Inhibin subunit alpha.

Table 2 Results of least absolute shrinkage and selection operator and support vector machine recursive feature elimination algorithms for the selection of immune-related feature genes.

LASSO		SVM-RFE
Variables	Coefficients	Variables	Accuracy	Kappa
Penk	4.895	Il12a	0.800	0.600
Il12a	0	Angptl4	0.900	0.800
Gcgr	-4.308	Rbp7	0.967	0.933
Rbp7	1.761	Sctr	0.900	0.800
Inha	-1.155	Penk	1.000	1.000
Sctr	0	Gcgr	1.000	1.000
Angptl4	2.382	Inha	1.000	1.000

LASSO: Least absolute shrinkage and selection operator; SVM-RFE: Support vector machine recursive feature elimination; Penk: Proenkephalin; Gcgr: Glucagon receptor; Rbp7: Retinol binding protein 7; Inha: Inhibin subunit alpha.

Regulatory network of IFG expression

To characterize the in vivo expression regulation patterns of the identified IFGs, we constructed transcription factor and miRNA regulatory networks. The identified joint network consisted of a total of 27 nodes (IFGs with the transcription factors and miRNAs identified as their regulators) and 30 edges (regulatory interactions between the nodes) containing 20 transcription factors and 3 miRNAs. Among these transcription factors, TFAP2A is a co-transcription factor for Gcgr, Rbp7, and Inha. Two miRNAs, miR-215-5p and miR-192-5p, regulated Penk, and miR-26b-5p regulated Rbp7 expression (Figure 3D). Therefore, targeting TFAP2A or specific miRNAs could modulate the expression of multiple DCM-related genes simultaneously, which provides a strategic approach for therapeutic intervention.

IFGs as potential biomarkers for DCM

To evaluate the ability of the IFGs to distinguish between disease and control samples, we used predictive and diagnostic models to validate their ability to be extrapolated to RNA-seq validation cohorts (GSE4745 and GSE197999). A logistic regression nomogram model was constructed using Penk, Gcgr, Rbp7, and Inha as independent variables and DCM as the dependent variable. The model exhibited excellent discriminatory ability, as indicated by a C-index of 1.0 (Figure 4A). The nomogram model was further validated using a calibration curve, which demonstrated strong calibration with a mean absolute error of 0.008. This result indicated good agreement between the predicted probabilities and observed outcomes (Figure 4B).

Figure 4 Clinical prediction models and diagnostic value for immune feature genes in the validation cohort. A: Nomogram predicting disease risk by proenkephalin (Penk), glucagon receptor (Gcgr), retinol binding protein 7(Rbp7), and inhibin subunit alpha (Inha); B: Calibration curves visually assessing the risk nomogram model calibration, a key aspect of model validity; C: Receiver operating characteristic (ROC) analysis of immune feature genes in diabetic cardiomyopathy (DCM) patients; D: ROC analysis of the diagnostic value of multivariable (four feature genes, unregulated Penk and Rbp7) in DCM patients; E: Decision curve analysis model construction was performed using a simple logistic regression model with Penk, Gcgr, Rbp7, and Inha as independent predictors and DCM as the outcome; a complex model (complex) was constructed by combining the four genes as predictors. Penk: Proenkephalin; Gcgr: Glucagon receptor; Rbp7: Retinol binding protein 7; Inha: Inhibin subunit alpha.

ROC curves assessed the discriminatory capacity of the four feature genes (Penk, Rbp7, Gcgr, and Inha), demonstrating their strong diagnostic potential with AUC values between 0.973 and 0.996 in DCM patients (Figure 4C). Logistic regression-based multivariable diagnoses also resulted in AUCs of 1.000 for the four genes and two up-regulated genes (Penk and Rbp7) in DCM (Figure 4D). DCA was employed to assess the clinical utility of various models in predicting disease outcomes. The decision curves for the four individual gene models and the complex model (Penk, Gcgr, Rbp7, and Inha) exhibited noticeable separation from the extreme curves representing "None" (no intervention) and "All" (treating all individuals) (Figure 4E). Within the threshold range of 0 to 1, the decision curve of the complex model consistently outperformed the curves of the four simple models and demonstrated a greater net benefit (Figure 4E).

Correlation between immune cell infiltration and IGF expression in DCM

SsGSEA allows for the estimation of the relative abundance of immune cell types in complex tissue samples based on gene expression data. To examine the immune cell infiltration characteristics in DCM, we used ssGSEA to obtain immune infiltration scores for disease and control samples. The ssGSEA results discovered statistically significant differences in the scores of five immune cell types, type 2 T helper cells, activated B cells, memory B cells, neutrophils, and myeloid-derived suppressor cells, between the disease and control samples. Type 2 T helper cells exhibited the most significant difference (P < 0.001), which indicated their potential involvement in DCM pathology (Figure 5A). Type 2 T helper cells are known for their role in humoral immunity and allergic responses and may contribute to inflammatory processes in DCM. Correlation analysis between the four selected IFGs and immune cells revealed significant positive correlations (P < 0.01) between Penk and Rbp7 and type 2 T helper cells, and Gcgr and Inha showed significant negative correlations (P < 0.01) with type 2 T helper cells (Figure 5B). These genes may be involved in pathways that regulate Th2 cell proliferation or function, which suggests a potential mechanism for the influence of inflammation on the pathology of DCM.

Figure 5 Single-sample gene set enrichment analysis was used to analyze immune cell infiltration in diabetic cardiomyopathy patients and its correlation with feature genes. A: Box plot of immune scores for diabetic cardiomyopathy patients and control; B: Heatmap of the correlation between feature genes and immune cell infiltration. ^aP < 0.05, ^bP < 0.01 ^cP < 0.001. DCM: Diabetic cardiomyopathy; CTL: Control; Penk: Proenkephalin; Gcgr: Glucagon receptor; Rbp7: Retinol binding protein 7; Inha: Inhibin subunit alpha.

Expression patterns of IFGs and intercellular profiles determined by scRNA-seq of diabetic mouse hearts

To explore the expression patterns of the four identified IFGs, scRNA-seq analysis was conducted on the hearts of diabetic mice (Table 1). The integrated dataset, which consisted of DCM and healthy control samples, contained a total of 30905 genes and 12593 cells and allowed for robust comparative analysis. The cellular analyses showed that Penk and Rbp7 significantly elevated in DCM, which suggested their active involvement in the disease process. Conversely, Inha and Gcgr were expressed at low levels or undetectable, which indicated their minimal direct involvement under the conditions studied (Figure 6A). Additional tissue expression data from the ProteomicsDB database confirmed the substantial expression of Penk and Rbp7 in the heart, particularly Rbp7, which reinforced the RNA-seq findings and highlighted the potential biological relevance of these genes in cardiac tissues (Figure 6B). Molecular experiments on AC16 cells demonstrated minimal Penk expression in untreated cells, with significant upregulation of Penk and Rbp7 after a 48 hours HG treatment (Figure 6C and D). This response was consistent with the increased cardiac stress and metabolic changes observed in DCM, which supported the disease relevance of Penk and Rbp7 expression observed in the scRNA-seq analysis.

Figure 6 Expression levels of immune feature genes according to single-cell RNA sequencing data. A: Violin plots comparing the expression levels of immune feature genes in control (CTL) (5833 cells) and diabetic cardiomyopathy (6760 cells) samples, with each point representing one cell; B: Tissue-specific expression of proenkephalin and retinol binding protein 7 in the ProteomicsDB repository; C: Relative gene expression levels in AC16 cells were examined using real time-quantitative polymerase chain reaction. Each point represents one independent biological experiment. HG-24 and HG-48, AC16 cells exposed to high glucose for 24 hours and 48 hours, respectively; D: Western blot showing the protein expression levels of the two genes. β-Tubulin and GAPDH were used as endogenous CTLs. DCM: Diabetic cardiomyopathy; CTL: Control; Penk: Proenkephalin; Gcgr: Glucagon receptor; Rbp7: Retinol binding protein 7; Inha: Inhibin subunit alpha.

We used UMAP dimensionality reduction and graph-based clustering and identified 14 cell clusters (resolution = 0.5) (Figure 7A) and 9 cell types (Figure 7B), which were annotated based on their highly expressed genes (DEGs) (Figure 7C). Analysis of feature gene expression patterns across these cell types revealed that Penk was primarily expressed in fibroblasts, and Rbp7 was enriched in endothelial cells and CMs. Notably, Rbp7, which is crucial for vascular function, was also detected in pericytes (Figure 7D). We also examined intercellular dysregulation. Using dimensionality reduction, we observed increased fibroblast numbers and reduced macrophage numbers (Figure 7E and F). DCM enhanced intercellular communication, primarily via increased fibroblast signaling (Figure 8A), which highlighted the strong interaction between fibroblasts and macrophages (Figure 8B). Signaling analysis revealed the involvement of the macrophage migration inhibitory factor (MIF) signaling pathway (Figure 8C). Weighted network analysis indicated that macrophages were significant MIF receivers, and pericytes were significant sender proteins (Figure 8D and E). Macrophages received MIF signals and communicated with fibroblasts via the growth arrest-specific pathway (Figure 8F). Global communication analysis identified four patterns that showed how sender (outgoing) and target (incoming) cells coordinated with specific signaling pathways, which may be crucial for understanding the pathological landscape of DCM (Figure 8G).

Figure 7 Analysis of the expression patterns of immune feature genes and cell-cell communication. A: UMAP dimensionality reduction embedding was performed on the integrated dataset of single-cell RNA sequencing data from all profiled samples (n = 12593 cells). The embedding is color-coded by inferred cluster identity and annotated with manual cell type labels; B: Clusters with similar marker gene expression modules were merged and clustered into the same cell type; C: Dot plot showing cell-specific markers; D: Violin plots showing statistically significant feature gene (proenkephalin and retinol binding protein 7) expression values for each cell type, grouped by the donor of origin [control (CTL) and diabetic cardiomyopathy (DCM)]; E: UMAP for cell clusters grouped by the donor of origin (CTL or DCM); F: Bar chart representing the cell counts of different cell types in the two groups. DCM: Diabetic cardiomyopathy; CTL: Control; Penk: Proenkephalin; Gcgr: Glucagon receptor; Rbp7: Retinol binding protein 7; Inha: Inhibin subunit alpha.

Figure 8 Analysis of cell-cell signaling pathway networks and communication patterns in diabetic cardiomyopathy. A: The circle plot shows the cell-cell communication analysis of the control (CTL) and diabetic cardiomyopathy (DCM) groups, with the ligand color in blue and the receptor color in red; B: The circle plot depicting the aggregated cell-cell communication network in DCM illustrates either the number of interactions or the total interaction strength (weights) between any two cell groups; C: A signaling role analysis was conducted on the aggregated cell-cell communication network across all signaling pathways to identify the signals contributing the most to outgoing or incoming signaling in specific cell groups. This was visualized using a combined heatmap, where each row corresponds to a signaling pathway; D: Heatmap visualizing the cell-cell communication network specifically based on the migration inhibitory factor (MIF) signaling pathways. The intensity of the color in each cell reflects the strength of the interaction between the corresponding cell groups. Each row corresponds to a sender cell group, and each column represents a receiver cell group; E: Heatmap displaying the computed centrality scores (importance) of cell groups in the intercellular communication network (dominant senders, receivers, mediators, and influencers) inferred from the MIF signaling pathway; F: Chord diagram illustrating the significant signaling pathways between selected source cell groups and target cell groups (source: Endothelial cells, fibroblasts, and pericytes, target: Macrophages) (left panel) in the intercellular communication network. The chords connecting the source and target arcs represent the significant signaling pathways between the corresponding cell groups; G: River plot to visualize associations of latent patterns (outgoing and incoming) with cell groups and signaling pathways. The specific communication patterns are defined by different colors in the network. MIF: Migration inhibitory factor; DCM: Diabetic cardiomyopathy; CTL: Control.

DISCUSSION

The present study integrated transcriptome sequencing data and applied machine learning algorithms to analyze immune gene datasets and identified four IRGs as the most relevant to DCM. Analysis of single-cell sequencing data and in vitro modeling using AC16 human cardiomyocytes exposed to HG further validated the upregulation of Penk and Rbp7, which highlighted their potential importance in DCM immunopathology. Single-cell data also suggested a regulatory role for phenotypic transformation and signaling communication between the extracellular mesenchyme and immune cells outside of myocardial working cells, particularly fibroblasts and macrophages, in the pathomechanism of DCM. The workflow of this study is summarized in Figure 9.

Figure 9 The workflow of this study. DCM: Diabetic cardiomyopathy; RNA-seq: RNA sequencing; WGCNA: Weighted gene co-expression network analysis; DEGs: differentially expressed genes; IRGs: Immune-related genes; DIRGs: Differential immune-related genes; LASSO: Least absolute shrinkage and selection operator; SVM-RFE: Support vector machine recursive feature elimination; ROC: Receiver operating characteristic; DCA: Decision curve analysis; ssGSEA: Single-sample gene set enrichment analysis.

We investigated dysregulated IRGs in DCM and identified the upregulation of Penk and Rbp7 as disease signatures. Single-cell expression pattern analysis did not reveal which class of immune cells was enriched based on the limited availability of current DCM single-cell data, which led to low-resolution heart cell clustering. Despite recognizing this limitation, the included single-cell datasets and exploration of public databases provided evidence that both genes were expressed in cardiomyocytes. Our molecular experiments further confirmed this result.

Penk encodes the precursor peptide enkephalin, which is an opioid-like peptide that is transiently expressed in immune cells and may regulate immune cell activity and immune responses[24]. An endogenous opioid peptide was also found within the central and autonomic nervous systems of various organs, such as the heart. Recent meta-analysis findings highlight a notable association between Penk and diabetes, indicating its potential as a predictive biomarker for the disease[25]. Our DCM diagnostic and predictive models constructed based on gene expression data consistently emphasized the value of these factors as biomarkers. Previous evidence also suggested that human bone marrow mesenchymal stromal cells expressed and secreted Penk under inflammatory conditions, which led to elevated levels of the anti-inflammatory and immunomodulatory factor IL-10[26]. Another molecule, Rbp7, is a target gene of PPARγ and plays a role in antioxidant activity in the vascular system, which suggests that it is a protective factor in cardiovascular health[27]. Rbp7 is a putative RNA-binding protein (RBP) implicated in cell cycle arrest and transmission competence[28]. RBPs play a key role as mediators and regulators of RNA-protein interactions, influencing a wide range of genetic, epigenetic, and metabolic activities in both immune and non-immune cells. They are essential for numerous cellular functions related to immunity and homeostasis[29]. However, the knowledge of these two factors in diabetes is limited, and this analysis did not address how gene dysregulation in cardiomyocytes mediated the pathological process of DCM. Characterization of the expression of the gene in the heart in a well-rounded set of immune cells and mechanistic exploration of its dysregulated expression in cardiomyocytes could provide fresh perspectives on the pathogenesis of DCM.

Recent studies highlighted the importance of immune dysregulation in DCM pathogenesis, which emphasized the role of immune cells and inflammatory mediators in driving cardiac dysfunction[30,31]. New research points to the role of non-immune cells, including endothelial and fibroblast cells, in controlling the inflammatory response and encouraging cardiac fibrosis in DCM[32,33]. These findings underscore the complexity of immune-mediated mechanisms in DCM and the need for comprehensive approaches to unravel its pathophysiology. Our study also revealed differences in immune cell infiltration in DCM, particularly alterations in T-helper (Th) 2 cells. Th cell subsets, including Th2 cells, play a crucial role in directing immune responses via the production of signature cytokines. In the context of diabetes, T lymphocytes, including Th and regulatory T cells, modulate cardiac inflammation[34,35]. Th1/Th2 cytokine imbalances are notable features of diabetes complications and emphasize the significant role of Th cells in DCM pathology[36-40].

The scRNA-seq analysis of diabetic mouse hearts provided insights into enhanced intercellular communication, particularly between fibroblasts and other cell populations, especially macrophages. The MIF signaling pathway plays an essential role in the communication patterns among these cells, and macrophages are the main targets of signaling and vascular pericytes are the source of signaling. These findings elucidate the complex interactions between different cell types and signaling pathways. Cardiac fibrosis, a key pathological aspect of DCM, results from complex interactions among immune cells, fibroblasts, and other host-derived cells[13,41,42]. Diverse cell types including cardiomyocytes, endothelial cells, and pericytes play roles in promoting fibrosis through various mechanisms, such as the production of proteases involved in matrix metabolism, the secretion of fibrotic mediators and matrix cell proteins, and the transition to a fibroblast phenotype[43]. Under the influence of various cytokines, monocytes/macrophages differentiate into myofibroblasts, which produce inflammatory mediators and pro-fibrotic factors to establish a detrimental pathological cycle[5,44,45]. Although not yet reported in DCM, similar phenomena have been observed in models of myocardial infarction and diabetic nephropathy[46,47]. The role of cardiac pericytes is not clear because these cells are not well characterized, and their bidirectional differentiation into myofibroblasts and endothelial cells is controversial[48]. Cardiac pericytes may contribute to diabetes-associated fibrosis via the secretion of inflammatory mediators and fibroblast-activating growth factors[42,49]. Macrophages activate TGF-β signaling in pericytes via the secretion of amphiregulin, which promotes pericyte differentiation into myofibroblasts and participates in tissue repair[50]. Further examination of the MIF signaling-mediated macrophage-pericyte axis in DCM may provide insights into how these cells contribute to immune surveillance and effector functions in this disease[51].

Overall, the present study identified dysregulated Penk and Rbp7 as potential DCM biomarkers, but their exact mechanism of action in DCM must be further examined. Aberrant mesenchymal-immune cell phenotypic transformation and communication may be essential for DCM pathogenesis. Our analysis was limited by the heterogeneity of single-cell data from mice, which led to low-resolution heart cell clustering and limited characterization of immune cell types enriched in Penk and Rbp7. Our study did not investigate the precise mechanisms by which dysregulated genes in cardiomyocytes mediate the pathological process of DCM in vivo. Future research should focus on elucidating the molecular mechanisms of the role of Penk and Rbp7 in DCM and further investigating the functional consequences of aberrant mesenchymal-immune cell communication in disease progression.

CONCLUSION

Our study elucidated the immunological molecular mechanisms underlying DCM and proposed potential diagnostic and prognostic biomarkers. WGCNA and machine learning identified Penk and Rbp7 as promising candidates with high diagnostic and prognostic value for DCM. Our findings also highlight the involvement of type 2 T helper cells in immune cell infiltration associated with DCM. scRNA-seq analysis suggested abnormal cellular phenotype transformation and signaling communication, particularly between mesenchymal cells and macrophages, which suggests a crucial role for aberrant cell communication in DCM pathogenesis. Overall, the present study offers valuable insights into the immunological mechanisms of DCM and suggests avenues for future research and potential therapeutic targets.

ACKNOWLEDGEMENTS

We thank everyone who provided public datasets and software to make this work possible.