Identification of hub genes associated with oxidative stress in heart failure and their correlation with immune infiltration using bioinformatics analysis

Data collection

Gene expression profiling (GSE36074) was obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The GSE36074 series on the platform GPL1261 (Affymetrix Mouse Genome 430 2.0 Array). The data included 7 TAC samples and five Sham samples.

Identification of oxidative stress-related differential genes

GEO2R online tool was used to screen the DEGs and oxidative stress-related genes (OSRGs) were extracted from the GeneCards website. DEGs were defined based on |Fold change |≥ 1 and adj.p-value of < 0.05. Then, we performed overlapping between DEGs and OSRGs to obtain OSRDEGs by the Draw Venn Diagram tool.

Enrichment analysis for OSRDEGs

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed to explore potential functions and crucial pathways of OSRDEGs using DAVID 6.8. A p-value of < 0.05 was considered as statistically different.

Screening of hub genes

The PPI network was established via the STRING database (http://string-db.org). Then, cytoHubba plugin in the Cytoscape was utilized to screen hub genes by the Degree algorithms. Hub genes were considered with degrees of ≥76.

Analysis of immune cell infiltration

CIBERSORTx (https://cibersortx.stanford.edu/) was applied to calculate the proportions of immune cells between the TAC and Sham samples. The expression profiles of 25 immune cell types in mouse tissues were obtained from Chen et al. (2017). Then, the percentages of each immune cell type were calculated and displayed in bar graph; Furthermore, we compared the levels of the 25 immune cell types in the TAC and Sham samples, and visualized the results using an online platform (https://www.bioinformatics.com.cn). To further understand the relationship between immune cells and hub genes, we conducted a Pearson correlation analysis and visualized the results through an online platform (https://www.bioinformatics.com.cn).

Animals

Ten specific pathogen-free male C57BL/6 mice (8–10 weeks, 20–23 g) were obtained from Yangzhou University. All mice were raised in Yangzhou University Comparative Medical Center with a controlled environment (26 °C, 55% humidity, 12 h light/dark cycle), and allowed access to feed and water ad libitum. The rearing environment was cleaned regularly by professionals. After acclimatization for seven days, mice were randomly allocated to either the Sham or TAC group. All mice were anesthetized with ketamine (80 mg/kg) and xylazine (5 mg/kg) by intraperitoneal injection and euthanized with cervical dislocation on the 28th day after TAC. All the animal experiments were approved by the Yangzhou University Ethics Committee (No. 20210398).

Transverse aortic constriction (TAC) mouse model

Mice were anesthetized with 80 mg/kg ketamine and 5 mg/kg xylazine, and when mice lost the response to foot squeeze, the surgical procedures were performed as previously described (Ma et al., 2016). Briefly, the transverse aortic arch was exposed by blunt dissection in the 2nd intercostal space, and ligated against a 26 G blunted needle with a 7–0 silk suture. Sham group mice underwent the same procedure without ligation of the aorta. No animals died after surgery. Cardiac function was analyzed by echocardiography on the 28th day after TAC. After echocardiography, the mice were euthanized, and their heart samples were collected for subsequent study.

Reverse transcription quantitative PCR (RT-qPCR)

The mRNA expression levels of hub genes were measured by RT-qPCR. Briefly, we used the RNA extraction reagent (cat. no. R401-01-AA; Vazyme, Nanjing, China) to isolate the total RNA; next, the RNA was reverse-transcribed into cDNA based on reverse transcription kit instructions (cat. no. RR036A; Takara, Dalian China). The primer sequences were listed in Table 1. The levels of hub gene mRNAs were normalized using GAPDH mRNA.

Immunohistochemistry (IHC)

IHC was performed as previously described (Kolkova et al., 2010). In brief, ventricular tissues were fixed in formalin and then embedded in paraffin, afterward, cut into 5-m thick sections. Lastly, the sections were stained with anti-Mapk14 (1:200, cat. no. AF6456; Affinity Biosciences, Jiangsu, China), anti-Hif1a (1:200, cat. no. #48085; Cell Signaling Technology, MA, USA), anti-Myc (1:100, cat. no. A11394; ABclonal, Wuhan, China), anti-Hsp90ab1 (1:200, cat. no. ET1605-56; Huabio, Hangzhou, China), and anti-Hsp90aa1 (1:200, Huabio, cat. no.ET1605-57; Huabio, Hangzhou, China) antibodies at 4 °C overnight. Lastly, the sections were incubated with the horseradish peroxidase-conjugated secondary antibody (1:500; cat. no. ab6721; Abcam, Cambridge, MA, USA) for 1 h at room temperature. Images of three representative fields in each section were captured using Leica DMi8 microscope (Leica Microsystems, German) at 200× magnification by a pathologist. For the reading of each antibody staining, a uniform setting for all the sections was applied. Then, Image-Pro Plus v6.0 software (Media Cybernetics Inc, Bethesda, MD) was used to measure the integrated optical density (IOD) of each image. Standard Optical Density was chosen in the Intense Calibration panel, and background was subtracted in the panel of Optical Density Calibration. The mean IOD (IOD/total area) was the result of cumulative optical density divided by total area.

Western blotting

The protein was extracted from heart tissues using RIPA lysis buffer; Then, the protein concentration was measured by the BCA protein assay kit (Beyotime Institute of Biotechnology, Jiangsu, China). Afterwards, the protein was separated on 10% SDS-PAGE, and transferred to PVDF membranes (Bio-Rad, Hercules, CA, USA), membranes were blocked with 5% skim milk powder (cat. no.BS102; Biosharp, Wuhan, China) for 2 h and then incubated with the anti-Mapk14 (1:1000, Affinity Biosciences, cat. no. AF6456, Jiangsu, China), anti-Hif1a (1:1000, cat. no. #48085; Cell Signaling Technology, MA, USA), anti-Myc (1:500, cat. no. A1309; ABclonal, Wuhan, China), anti-Hsp90ab1 (1:2000, cat. no. ET1605-56; Huabio, Hangzhou, China), anti-Hsp90aa1 (1:1000, cat. no. ET1605-57; Huabio, Hangzhou, China) and anti-Tubulin (1:1000, cat. no. AC030; ABclonal, Wuhan, China) primary antibodies at 4 °C overnight. Subsequently, membranes were washed with TBST and incubated with the HRP-labeled secondary antibody (1:3000; cat.no. S0001; Affinity Biosciences, Jiangsu, China) for 2 h. The bands were visualized using enhanced chemiluminescence reagents (Thermo Fisher Scientific, Waltham, MA, USA). Protein expression was quantified by densitometry using Image J (NIH, Bethesda, MD, USA).

Statistical analysis

Data were expressed as mean ± standard deviation. Statistical analysis was performed using GraphPad Prism (version 7.0). For two-group comparisons, a two-sample, unpaired Student t-test was used. P < 0.05 was considered as statistically different.

Identification of OSDEGs in HF

The flow diagram was provided in Fig. S1. The expression profile (GSE36074) was enrolled in this study, and the dataset included five Sham and seven TAC heart samples. OSRGs were screened from the GeneCards websites. A total of 1,918 DEGs were screened from the GSE36074 dataset (Fig. 1A). Moreover, a total of 899 OSRGs were found in GeneCards, and 136 OSRDEGs were screened by taking the intersection of OSRGs and DEGs (Fig. 1B). The clustering results of the OSRDEGs showed that samples from the same group clustered together (Fig. 1C).

Enrichment analysis

To gain a deeper understanding of the function of OSRDEGs, we conducted GO and KEGG enrichment analyses. Our GO analysis revealed that OSRDEGs were mostly associated with mitochondrion and response to oxidative stress (Fig. 2A). Additionally, our KEGG analysis demonstrated that OSRDEGs were related to HIF-1, FoxO, MAPK, and TNF signaling pathway (Fig. 2B).

Screening hub genes

To identify hub genes in HF, we constructed a PPI network of OSRDEGs using the STRING tool (Fig. 3A.). We then employed the Degree algorithms in Cytoscape to calculate each gene’s degree, ultimately identifying five hub genes (Mapk14, Hif1a, Myc, Hsp90ab1, and Hsp90aa1) (Fig. 3B).

Analyzes of immune infiltration

Immune infiltration in Sham and TAC groups was explored with the 25 immune cell types. Figure 4A showed the composition of the 25 immune cell types in each sample. As shown in Fig. 4B, compared with the Sham group, monocyte in the TAC group were relatively infiltrated less. To investigate the correlation between the hub genes and immune cells, we conducted Pearson correlation analyses. The result revealed that Mapk14 was positively correlated with Th2 Cells (R = 0.68, p = 0.015), while Hif1a and Hsp90ab1exhibited a negative correlation with Th2 Cells (R = − 0.68, p = 0.016; R = − 0.76, p = 0.004). Additionally, Myc exhibited a negative correlation with Monocyte (R = − 0.66, p = 0.019), whereas, Hsp90aa1 was negatively correlated with NK Resting (R = − 0.67, p = 0.017) (Fig. 5).

Validation of key genes

To validate the findings from our bioinformatics analysis, we established a transverse aortic constriction (TAC) mouse model. After 4 weeks of TAC, the cardiac function of mice was evaluated with echocardiography. The result showed that left ventricular ejection fraction (LVEF) in the TAC group was significantly decreased compared with the sham group (Fig. S2); Based on the results, we have successfully established a TAC-induced HF model in mice.

The RT-qPCR results showed that Mapk14 mRNA levels were remarkably reduced, while Hif1a, Myc, Hsp90ab1, and Hsp90aa1 mRNA levels were significant increased in the TAC group than in the Sham group (Fig. 6). Furthermore, we further verified those hub genes at the protein level by IHC (Figs. 7A–7B) and western blot (Figs. 8A–8B). It was observed that the TAC group showed lower expression of Mapk14 protein (p < 0.05), while Hif1a, Myc, Hsp90aa1 and Hsp90ab1 proteins were significantly upregulated compared to the Sham group (p < 0.05). Those results supported bioinformatic analysis.