Alterated gene expression in dilated cardiomyopathy after left ventricular assist device support by bioinformatics analysis

Introduction Heart transplantation is the best treatment for end-stage dilated cardiomyopathy (DCM). Left ventricular assist device (LVAD) support is becoming more prevalent and may delay heart transplantation. Gene expression of the left ventricular myocardium usually changes following LVAD implantation. In this study, we aimed to identify potential biomarkers to determine the prognosis of patients with DCM after receiving LVAD support. Methods We extracted microarray datasets from Gene Expression Omnibus (GEO), including GSE430 and GSE21610. There were 28 paired DCM samples in the GSE430 and GSE21610 profiles. Differentially expressed genes (DEGs) were identified at LVAD implantation and heart transplantation. DEGs were annotated according to Gene Ontology (GO) and analyzed according to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. A protein–protein interaction (PPI) network was constructed. The top 10 crucial genes were predicted using Cytoscape plugin CytoHubba in conformity with the network degree algorithm. The levels of gene expression and the diagnostic values of crucial genes were confirmed in the clinical datasets. Results The 28 DEGs were clustered into the GSE datasets. GO annotations and KEGG pathway enrichment analyses revealed that inflammation might be involved. They were associated with correlative inflammation. Combined with PPI networks, these results revealed CytoHubba's top 10 hub genes, including CCL2, CXCL12, CXCL1, CTGF/CCN2, CX3CR1, POSTN, FKBP5, SELE, AIF1, and BMP2. Among them, CCL2, CXCL12, FKBP5, and BMP2 might be considered prognostic and diagnostic biomarkers after LVAD support and have confirmed their validity in clinical datasets. The area under the curve of the four main hub genes was more than 0.85, indicating high diagnostic ability and good prognosis for patients with DCM with LVAD implantation. However, a significant effect of CCL2, CXCL12, FKBP5, and BMP2 expression was not observed on the left ventricular end-diastolic diameter (LVEDD), left ventricular ejection fraction (LVEF), cardiac index (CI), or support time of LVAD. Conclusion CCL2, CXCL12, FKBP5, and BMP2 could be potential gene biomarkers for patients with DCM after LVAD support. These findings provide critical clues for the therapeutic management of patients with DCM and LVADs. LVEDD, LVEF, CI, and support time of LVAD were not correlated with the expression of these hub genes.


Introduction
Dilated cardiomyopathy (DCM) is a common manifestation of end-stage heart disease characterized by left ventricle dilation and impaired contractility without secondary causes (1). High mortality of DCM suggests that heart transplantation is the best treatment for end-stage DCM. However, a shortage of heart donors increases waiting times for patients to undergo heart transplantation. Left ventricular assist device (LVAD) is becoming a bridge to heart transplantation or a destination therapy for patients with DCM and end-stage congestive heart failure (CHF) (2, 3). The use of LVAD significantly decreases mortality and improves the quality of life of patients with CHF (4,5). Some patients' symptoms of heart failure improve after implantation of a ventricular assist device, which can delay the time of heart transplantation or even prevent it.
However, due to the different support times of LVAD among individuals, we hypothesized that specific biomarkers could predict the prognosis after implantation of assist devices and intervene early. Gene expression in the left ventricular myocardium changed after LVAD implantation (6), and changes in genes may be another manifestation of CHF (7). Changes in gene expression could be considered as a method to detect the progression of the disease and cardiac remodeling after LVAD implantation. Some studies have investigated molecular levels after mechanical unloading. However, there has been little agreement on the molecular mechanisms of LVAD support. A better understanding of the gene expression patterns would provide a new perspective on the potential prognosis and treatment of DCM. Our study indicates that CCL2, CXCL12, FKBP5, and BMP2 might be related to clinical outcomes. These findings provide critical clues to the therapeutic strategies for DCM with end-stage CHF.

Data source
Gene expression profiles [GSE430 (6) and GSE21610 (7)] were obtained from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database (http://www. ncbi.nlm.nih.gov/geo/). The study participants were patients with DCM who did not experience secondary myocardial diseases, such as ischemia, valvulopathies, and congenital disorders. The GSE430 dataset consists of seven paired DCM and end-stage heart failure samples with LVAD support. The GSE21610 profile contains 21 paired samples of DCM and 9 paired samples from ischemic cardiomyopathy (ICM) with LVAD support. We analyzed only 28 paired DCM samples from the GSE430 and GSE21610 datasets. All data are freely available from GEO. The GSE430 dataset was based on the GPL96 platform [(HG-U133A) Affymetrix Human Genome U133A Array], while the GSE21610 dataset was based on the GPL570 platform [(HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array]. This study did not have a relation with any experiments conducted by any authors on humans or animals. No ethical approval was required, given the lack of human or animal subjects in this study.

Identification of differentially expressed genes
GEO2R, an online tool (http://www.ncbi.nlm.nih.gov/geo/geo2r/) with an R-based web (R 3.2.3, Biobase 2.30.0, GEOquery 2.40.0, limma 3.26.8), was used to analyze GEO data to screen out DEGs (A) Volcanic map of the GSE430 dataset; (B) volcanic map of the GSE21610 dataset. All genes were presented in volcano plots. Red ones were upregulated DEGs, green ones were downregulated DEGs, and black ones were the rest of the genes. The name of the top 10 upregulated and downregulated DEGs are shown in the volcanic plots. DEGs: differentially expressed genes.
Wei et al. 10.3389/fcvm.2023.1013057 Frontiers in Cardiovascular Medicine between pre-LVAD and post-LVAD groups in the dataset. A false discovery rate P value <.05 and |log2fold change (FC) |≥1 were set as the standard for identifying DEGs. Log2FC ≥1 for upregulated genes and log2FC ≤ −1 for downregulated genes were considered statistically significant. Overlapping DEGs were analyzed by a Venn diagram. Venn diagrams and heatmaps of DEGs were produced using TBtools, an integrative toolkit for interactive biological data (8). A free online data analysis tool, Sangerbox, was used to analyze volcano data (http://www.sangerbox.com/tool) (9), and the gene names of the top 10 upreregulated DEGs and downreregulated DEGs were indicated in the volcano plots.

GO and KEGG pathway analysis
GO annotation and KEGG pathway analyses were conducted to determine the functions of DEGs on David (http://david.ncifcrf. gov) (10,11), an online tool for gene annotation and analysis. Biochemical processes (BPs), cellular components (CCs), and molecular functions (MFs) were included in the GO annotation analysis. P value < .05 was used as an enrichment threshold for GO terms and KEGG pathway segments. Sangerbox was used to visualize GO enrichment and KEGG pathway analyses (9).

PPI network construction and hub gene screening
Protein-protein interaction (PPI) networks help construct a functional protein system by combining all the protein-coding genes into one huge biological network. In the study, by mapping DEGs into the database, we found important protein pairs with a minimum required interaction score of 0.4 and removed the hidden disconnected nodes. An evaluation of the PPI network of DEGs was conducted using the Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org/), which integrates associations for the PPI network (12-14). Using Cytoscape software version 3.9.1, the PPI network was visualized by calculating the degree of protein connectivity. CytoHubba is a plugin for Cytoscape that uses multiple topological algorithms to rank nodes in a PPI network. We used Cytohubba software to identify the central genes.The top 10 hub genes were selected using the degree method because they perform better in predicting essential proteins accurately. The hub genes were analyzed using the GeneMANIA online database (https:// genemania.org/).

Impact on the characterization of clinical data
The expression of the top 10 hub genes in the pre-and post-LVAD groups was analyzed using GraphPad Prism 8 software. We also studied the relationship between the expression of the top 10 hub genes and clinical information.

Statistical analysis
GraphPad Prism 8 software was used to perform the statistical analysis. Measures of continuous variation were expressed as means ± standard deviations (SDs), and categorical data were expressed as numbers or percentages. Continuous variables were compared by Student's t-test if the data were normally distributed; otherwise, they were compared by nonparametric statistics. Categorical data were analyzed by the chi-square test. Statistical significance was set at P < .05, and all statistical analyses were two-sided.

FIGURE 2
Venn diagram of overlapping of the differentially expressed genes in the datasets. Heatmap of the differentially expressed genes in the datasets.   Protein-protein interaction network of differentially expressed genes created by STRING.

Screening for DGES
The GSE430 and GSE21610 datasets were selected and screened for DEGs using GEO2R. By screening using the standard of significance (P value <.05 and |log2FC| ≥ 1), GSE430 identified 198 DEGs (76 upregulated and 122 downregulated) and GSE21610 identified 330 DEGs (151 upregulated and 179 downregulated). All genes are presented in volcano plots ( Figures 1A,B). Volcano analysis was performed using Sangerbox tools, an online data analysis platform that is free. The names of the top 10 upregulated and downregulated DEGs are shown in the volcano plots. There were 28 DEGs in the overlapping dataset ( Figure 2). A heatmap of DEGs in the datasets is shown (Figure 3). Network of differentially expressed genes created by Cytoscape.

GO and KEGG functional enrichment analyses
Functional enrichment analyses of the DEGs were performed using the Database for Annotation, Visualization, and Integrated Discovery (DAVID; http://www.david.org). In the GO analysis results, BP, MF, and CC appeared to be significantly enriched in DEGs. BP ( Figure 4A) analysis majored in the inflammatory response, G-protein-coupled receptor signaling pathway, negative regulation of cell proliferation, cell adhesion, cellular response to lipopolysaccharide, etc. In CC ( Figure 4B), these DEGs were associated with extracellular space, extracellular region, external side of the plasma membrane, cytoplasm, excitatory synapse, cell surface, and nucleus. In MF ( Figure 4C), the analysis indicated growth factor activity, chemokine activity, receptor binding, and CXCR chemokine receptor binding. KEGG pathway analysis showed that the DEGs were significantly enriched in viral protein interaction with cytokine and cytokine receptor, NF-kappa B signaling pathway, cytokine-cytokine receptor interaction, chemokine signaling pathway, lipid and atherosclerosis, rheumatoid arthritis, AGE-RAGE signaling pathway in diabetic complications, amoebiasis, TNF signaling pathway, and pathways in cancer ( Figure 4D). The GO enrichment and KEGG analyses were visualized using Sangerbox tools.

Building PPI networks and identifying hub genes
DEGs were evaluated based on the connectivity degree of the PPI network, which was built using the STRING database ( Figure 5). The top 10 genes were considered hub genes according to the degree algorithms in the Cytoscape soft. Among the significant genes, those with the highest connectivity were identified, and the relationships between the hub genes are shown in Figure 6. The top 10 hub genes were C-C motif chemokine ligand 2 (CCL2), C-X-C motif chemokine ligand 12 (CXCL12), C-X-C motif chemokine ligand 1 (CXCL1), connective tissue growth factor (CTGF/CCN2), C-X3-C motif chemokine receptor 1 (CX3CR1), periostin (POSTN), FK506binding protein 5 (FKBP5), selectin E (SELE), allograft inflammatory factor 1 (AIF1), and bone morphogenetic protein 2 (BMP2). The results of the top 10 hub genes' eight algorithms are presented in Table 1. Table 2 presents a description of the top 10 hub genes. The PPI network of the top 10 hub genes was constructed using GeneMania (http://genemania.org/) (Figure 7). This indicated a relationship between the top 10 hub genes and other related genes. CTGF is also known as cellular communication network factor 2 (CCN2). CNN2 represents the CTGF, as shown in Figure 7. Among them, CTGF, POSTN, FKBP5, SELE, AIF1, and BMP2 were upregulated, while CCL2, CXCL12, CXCL1, and CX3CR1 were downregulated.

Impact of hub genes on the characterization of clinical data
Patients with DCM obtained from the GSE430 and GSE21610 datasets received LVAD support and heart transplantion. The characteristics of clinical data are presented in Table 3. The expression of the top 10 hub genes was assessed between preand post-LVAD groups in both datasets (Figures 8, 9). Only CCL2, CXCL12, FKBP5, and BMP2 hub genes were significant in both clinical data (P < .05). A comparison of clinical data between the high and low expression levels of CCL2, CXCL12, FKBP5, and BMP2 hub genes was conducted. There was no statistical difference between the different expression levels of the four hub genes and left ventricular end-diastolic diameter (LVEDD), left ventricular ejection fraction (LVEF), cardiac index (CI), or LVAD support time (P > 0.05) ( Figure 10). However, LVEF and CI were significantly different between the GSE430 and GSE21610 datasets (P < 0.05) ( Table 3). CCL2, CXCL12, CX3CR1, FKBP5, SELE, and BMP2 demonstrated powerful diagnostic ability, with an area under the curve (AUC) > 0.85. AUCs for CCL2, CXCL12, and FKBP5 were >0.90 ( Figure 11).

Discussion
Heart transplantation is the best treatment for advanced heart failure caused by DCM. LVAD is increasingly being used as a bridge to transplantation, recovery, or destination therapy (2, 3, 15) and is a sustainable alternative to heart transplantation. Implantation of LVAD can slow down the progress of heart failure and improve the symptoms of heart failure. In addition, it would decrease mortality and improve the quality of life of patients (4,5). It has been proven that LVAD can lead to the reduction of ventricular diameter, reduction of LV mass, and improvement of cardiac ejection fraction (16,17), and patients may experience different pathophysiological processes after LVAD support. Mechanical unloading induced significant Frontiers in Cardiovascular Medicine changes in myocardial gene expression (18). Due to the different support times of LVAD in different individuals, we expect to find some biomarkers to predict the prognosis of heart failure after LVAD support and intervention as soon as possible. Some evidence indicates that the changes in gene expression in the left ventricular myocardium may be another manifestation of CHF following LVAD implantation (6,7). In this study, gene expression profiles of pre-and post-LVAD from DCM populations were compared using comprehensive bioinformatics.
All the patients underwent cardiac transplantation after LVAD support. We screened 16 upregulated and 12 downregulated DEGs using GO and KEGG analyses. GO analysis of DEGs revealed the specific mechanisms of myocardial function after LVAD support and indirectly suggested that some molecular mechanisms might be involved in the LVAD support. However, further studies are required to confirm and validate these findings. Gene ontology analysis suggested that several mechanisms of inflammatory response activation including inflammatory response, G-protein coupled receptor signaling pathway, negative regulation of cell proliferation, cell adhesion, FIGURE 7 Protein-protein interaction network of the top 10 hub genes. Inner circles represent the hub genes, and outer circles correspond to GeneMANIAproposed hub genes.   Frontiers in Cardiovascular Medicine pathway, chemokine signaling pathway, lipid and atherosclerosis, rheumatoid arthritis, AGE-RAGE signaling pathway in diabetic complications, amoebiasis, and TNF signaling pathway. KEGG enrichment was associated with inflammation-related pathways. An inflammatory cascade following LVAD implantation is essential. Functional enrichment analysis demonstrated that cytokine-related receptor interaction was a consequential pathway detected as a potential candidate for recent polygenic human adaptation (19). Grosman-Rimon et al. found that inflammatory markers increased over time with long-term LVAD support (20).
The results of GO and KEGG enrichment analyses in our study showed that some genes responsible for the production of inflammatory proteins were overexpressed after LVAD implantation. The top 10 hub genes were identified using Cytohubba in this study. These are CCL2, CXCL1, CXCL12, CTGF, CX3CR1, POSTN, FIGURE 10 Patients of LVEDD, LVEF, CI, and the support time of LVAD between high-expression groups and low-expression groups of CCL2, CXCL12, FKBP5, and BMP2 hub genes. High and rRed ones represented high-expression groups. The low and blue ones represented low-expression groups. LVEDD, left ventricular end-diastolic diameter; LVEF, left ventricular ejection fraction; CI, cardiac index.
Wei et al. 10.3389/fcvm.2023.1013057 Frontiers in Cardiovascular Medicine FKBP5, SELE, AIF1, and BMP2. We found that the expression levels of CCL2, CXCL1, CXCL12, and CX3CR1 were obviously downregulated and those of CTGF, POSTN, FKBP5, SELE, AIF1, and BMP2 were upregulated after LVAD implantation. The expression levels of CCL2, CXCL12, and FKBP5 were remarkably different between pre-and post-LVAD. The three main hub genes demonstrated powerful diagnostic abilities with an AUC >0.90. CCL2, CXCL12, and FKBP5 could be considered biomarkers for predicting and evaluating the effects of LVAD transplantation.
Our study shows that CCL2 was regarded as one of the marvelous hub genes. The expression level of CCL2 significantly decreased after LVAD implantation in our study. CCL2 is an inflammatory chemokine that plays an important role in the inflammatory response cascade (21). It could induce inflammation-associated cell aggregation and stimulate the activity of monocytes and basophils after LVAD implantation (22). Grosman-Rimon et al. found that markers of inflammation remained higher than before LVAD implantation at 3, 6, and 9 months after LVAD implantation (20). CCL2 may be involved in cardiac remolding and protects cardiomyocytes against cell death mediated by hypoxia. CCL2 showed a powerful diagnostic ability with an AUC > 0.90. The expression level of CCL2 seems to predict the severity of heart failure. CCL2 could be considered an ideal biomarker for LVAD support. Li et al. reported that CCL2 is a crucial hub gene in LVAD pathophysiology (23).
CXCL12 is an encoded protein and chemokine ligand. It binds to a G-protein-coupled receptor and CXC receptor 4 and plays a variety of important roles in many diverse cellular functions, including immune surveillance, inflammation response, tissue homeostasis, and tumor growth and metastasis (24)(25)(26)(27). In our study, we showed that CXCL12 might be involved in BP, including chemokine-mediated signaling pathway, positive regulation of monocyte chemotaxis, cell adhesion, response to hypoxia, response to peptide hormone, G-protein-coupled receptor signaling pathway, positive regulation of T cell migration, positive regulation of calcium ion import, blood circulation, and positive regulation of cell migration. The MF of CXCL12 includes growth factor activity, chemokine activity, receptor binding, and CXCR chemokine receptor binding. The AUC of CXCL12 was also >0.90 in our study, suggesting that it has a powerful diagnostic capability. At least three independent transcriptional studies of LVADs have detected CCL2 and CXCL12, indicating that mechanical unloading may influence gene expression regardless of the changes in cardiac function (28).
Multiple tissues express FKBP5 with different distributions, which plays an important role in cellular processes. FKBP5 exhibits peptide-prolyl isomerase activity and regulates protein folding. It influences the steroid receptor signal and NF-κB pathway. Stress may also result in the chemotaxis of proinflammatory cells and peripheral inflammation, possibly leading to a higher proinflammatory status and increased cardiovascular risk (29).
FKBP5 regulates glucocorticoid receptors. High levels of FKBP5 reduce effective glucocorticoid receptors and accelerate glucocorticoid resistance. Glucocorticoid resistance weakens the compensatory ability of the heart under stressful conditions. The expression of FKBP5 is fundamentally associated with stressrelated diseases. When exposed to high pressure, patients exhibit high FKBP5 expression and are prone to inflammation and acute cardiovascular events (29). The FKBP5 gene product is a potential biomarker and therapeutic agent for heart failure (30) and may be used to diagnose heart failure with an LVAD. FKBP5 displayed a high degree of diagnostic ability with an AUC of 0.9977.
Bone morphogenetic protein-2 (BMP2) is a low-molecularweight glycoprotein classified as a morphogen. It belongs to the transforming growth factor-β (TGF-β) superfamily. BMP2 plays an important role in vascular biology (31, 32), and increased BMP-2 levels can promote inflammation and atherosclerosis Frontiers in Cardiovascular Medicine through oxidative stress and endothelial dysfunction. BMP2 is involved in inflammation and fibrogenesis during heart remolding. BMP2 also demonstrated a high diagnostic ability in our study, with an AUC of 0.8701. LVAD support time was longer for high-expression CCL2 and CXCL12 groups than that for low-expression groups, while it was shorter for the high-expression FKBP5 and BMP2 groups in our study. However, there was no statistically significant difference in the duration between the groups with high and low expression levels (P > 0.05). LVEDD, LVEF, and CI showed no difference in hub gene expression levels.
The individuals included in our study were patients with DCM who underwent heart transplantation after LVAD implantation. LVAD unloading reduced the levels of inflammatory factor transcripts in patients with DCM. CCL2 and CXCL12 expression was low, whereas FKBP5 and BMP2 expression was high. This suggests reverse remodeling of the myocardium caused by mechanical unloading or deterioration of the patient's condition. These genes could be considered biomarkers for predicting and evaluating the effects of LVAD implantation in patients with DCM.
Despite discovering some potential genes associated with LVAD support at the end stage of CHF from DCM, the bioinformatics analysis of the study has some limitations. The number of patients investigated was small. We obtained only 28 paired myocardium microarrays from patients with DCM who underwent LVAD implantation and heart transplantation. It is not easy to obtain microarray datasets from the myocardium of patients with DCM who underwent LVAD implantation and heart transplantation from GEO datasets. This is an inherent limitation of bioinformatics analysis. The results from a bioinformatics analysis should be further validated using a larger sample cohort. In addition, the GEO dataset did not provide sufficient clinical information for this special population. Although the bioinformatic analysis has the aforementioned shortcomings, it can still provide a better understanding of LVAD pathophysiology from a gene expression perspective.
CCL2, CXCL12, FKBP5, and BMP2 could be potential gene biomarkers for patients with DCM after LVAD support. The expression level of hub genes was not statistically significant in LVEDD, LVEF, CI, and the support time of LVAD. Nevertheless, it may still be helpful to therapeutically manage DCM patients with LVADs.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.