Exploration and Validation of Potential Biomarkers and Therapeutic Targets in Ferroptosis of Asthma

Abstract

Purpose

Asthma is a chronic inflammatory airway disease involving multiple mechanisms, of which ferroptosis is a form of programmed cell death. Recent studies have shown that ferroptosis may play a crucial role in the pathogenesis of asthma, but no specific ferroptosis gene has been found in asthma, and the exact mechanism is still unclear. The present study aimed to screen ferroptosis genes associated with asthma and find therapeutic targets, in order to contribute a new clue for the diagnosis and therapy of asthma.

Methods

Ferroptosis-related differentially expressed genes (FR-DEGs) in asthma were selected by the GSE41861, GSE43696 and ferroptosis datasets. Next, the FR-DEGs were subjected by GO and KEGG enrichment, and the mRNA-miRNA network was constructed. Then, GSEA and GSVA enrichment analysis and Immune infiltration analysis were performed, followed by targeted drug prediction. Finally, the expression of FR-DEGs was confirmed using GSE63142 dataset and RT-PCR assay.

Results

We found 13 FR-DEGs by the GSE41861, GSE43696 and ferroptosis database. Functional enrichment analysis revealed that the 13 FR-DEGs were enriched in oxidative stress, immune response, ferroptosis, lysosome, necrosis, apoptosis etc. Moreover, our results revealed the mRNA-miRNA network of the FR-DEGs and identified candidate drugs. Also, immune infiltration revealed that ELAVL1, CREB5, CBR1 and NR1D2 are associated with the immune cells and may be potential targets in asthma. Finally, 10 FR-DEGs were validated by the GSE63142 database. It was verified that 7 FR-DEGs were differentially expressed by collecting asthma patients and healthy controls.

Conclusion

This study ultimately identified 7 FR-DEGs for the diagnosis and therapy of asthma. These 7 FR-DEGs contribute to oxidative stress and immune responses. This study provides potential therapeutic targets and biomarkers for asthma patients, shedding further light on the pathogenesis of asthma as well as providing new insights into the treatment of asthma.

Keywords:

• asthma
• bioinformatics
• ferroptosis
• GSEA
• GSVA
• immune microenvironment

Introduction

Asthma is an immune-related chronic airway inflammatory disease with major pathological changes including airway inflammation, smooth muscle dysfunction and airway remodeling with progressively increased morbidity and mortality in recent years.¹ According to the 2022 Global Initiative for Asthma (GINA) guidelines,² the common clinical phenotypes of asthma are allergic asthma, non-allergic asthma, adult-onset (late-onset) asthma, asthma with persistent airflow limitation, and asthma with obesity, and treatment guided by phenotype is clinically important, but no specific pathological features between clinical phenotype and response to treatment have been identified. The etiology of asthma is intricate, while current treatments work well to control symptoms in many patients, a large number of patients continue to suffer from worsening or treatment-related side effects.³ Therefore, it is critical to identify potential therapeutic targets for asthma, so as to offer new directions for the diagnosis, treatment and prognosis of asthma.

Ferroptosis is an iron-dependent and non-apoptotic form of programmed cell death that is symbolized by the formation of lipid peroxides from unsaturated fatty acids after catalysis by ester oxygenase or divalent iron,⁴ which results in cell death once excess lipid peroxides are accumulated.⁵ During the morphology of ferroptosis, the cell membrane and nucleus are intact, but the matrix is reduced or absent, and there is an increase in membrane density and a decrease in mitochondrial volume, which distinguishes it from other forms of cell death or apoptosis.⁶ The specific mechanism of ferroptosis is currently unknown and may be intricately linked to a number of biological pathways, the most well-known of which are the iron, lipid, and amino acid metabolic pathways.^7–9

Ferroptosis has recently been demonstrated to be essential in many diseases.^10–14 As a result, ferroptosis can be used to discover new diagnostic markers for diseases, treatments, and means of prevention, and ferroptosis has become a hot study. There are few studies on ferroptosis and asthma, and the only studies that have been conducted suggest that ferroptosis plays a very important role in the pathogenesis of asthma,^15,16 but the specific mechanisms involved have not been uncovered. Some studies have shown that in the pathogenesis of asthma, reactive oxygen species (ROS) dependent lipid peroxidation and destruction of iron homeostasis can promote ferroptosis in bronchial epithelial cells, and ferrostatin-1, an ferroptosis inhibitor, has been found to reverse the process of lipid peroxidation and ferroptosis.¹⁷ It has also been found in other studies that increased levels of lipid peroxidation and ROS are accompanied by increased levels of lung ferroptosis in asthma models. Ferrostatin-1 can inhibit the ferroptosis process in asthma. These findings suggest that ferroptosis play an important role in asthma, inhibiting ferroptosis process may be a potential target for the treatment of asthma.^16,18,19 In addition to ferrostatin-1, liproxstatin-1 and quercetin has also been found to alleviate asthma airway inflammation by inhibiting the ferroptosis process.^20,21 Wang found that ferroptosis could be inhibited in asthma by inhibiting M1 macrophage polarization, thereby improving airway inflammation associated with ferroptosis.²¹ These studies suggest that ferroptosis plays an important role in the pathogenesis of asthma by regulating oxidative stress and immune response, but in the study of ferroptosis in asthma it is limited to conventional markers of ferroptosis, and no asthma-specific ferroptosis genes have been identified. Therefore, there is an urgent need to further explore the ferroptosis genes related to asthma. Further study of the relationship between asthma and ferroptosis can identify asthma-specific ferroptosis genes, which can be used as the regulator of ferroptosis in asthma and provide new therapeutic targets for asthma patients.

In view of the lack of studies on ferroptosis in asthma and the specific mechanism is unclear, we explored the ferroptosis genes related to asthma through bioinformatics to lay the foundation for the subsequent mechanism research. This is the first study to link ferroptosis-related genes to asthma, which can reflect asthma pathology. The study screened out genes that are differentially expressed in asthma by using the GSE41861 and GSE43696 databases in the Gene Expression Omnibus (GEO) database, and further extracted co- FR-DEGs as target genes. Subsequently, functional analysis of FR-DEGs in Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome pathways was performed to determine the pathogenesis of these genes in asthma. Constructed mRNA-miRNA network and predicted the target drug of the genes. Then, Gene Set Enrichment Analysis (GSEA) and Gene Set Variation Analysis (GSVA) enrichment analysis and immune landscape analysis were performed, followed by targeted drug prediction. The FR-DEGs of asthma were further validated in the microarray dataset GSE63142. Finally, differential expression levels of these ferroptosis genes were verified in asthma and healthy individuals by RT-PCR.

Materials and Methods

Data Collection and Analysis

The FerrDb (http://www.zhounan.org/ferrdb/current; last accessed August 6, 2022) database contained 565 ferroptosis-related genes (Table S1). The mRNA expression profile datasets GSE41861, GSE43696 and GSE63142 in asthma were downloaded from GEO (https://www.ncbi.nlm.nih.gov/geo/; last accessed August 12, 2022) database. A total of 138 samples (91 asthma and 47 controls) are from the GSE41861 dataset; The GSE43696 dataset contains 58 samples (38 asthma and 20 controls); GSE63142 dataset with 83 samples (56 asthma and 27 controls). FR-DEGs were screened from GSE41861, GSE43696 and FerrDb databases in asthma, and finally the FR-DEGs in asthma was verified in the GSE63142 database.

Table 1 shows the specific information of the three datasets (GSE41861, GSE43696 and GSE63142), and Figure 1 shows the flow chart of this study.

Table 1 Details of GEO Asthma Date

Accession	Platform	Species	Asthma	Non-Asthma	Gene
GSE41861	GPL570	Human	91	47	mRNA
GSE43696	GPL6480	Human	38	20	mRNA
GSE63142	GPL6480	Human	56	27	mRNA

Abbreviation: GEO, Gene Expression Omnibus.

Figure 1 The overall protocol of this study.

Abbreviations: GO, differentially expressed genes; FR-DEGs, Ferroptosis-related differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; GSEA, Gene Set Enrichment Analysis; GSVA, Gene Set Variation Analysis.

Differentially Expressed Genes

Microarray data and annotation files from the GSE41861 and GSE43696 datasets were downloaded, additionally, the annotation files were used to annotate the probes. Then “GEO2R” packages including “GEO query” and “R/limma” packages were used to analyze the differentially expressed genes (DEGs), and the final determined DEGs were conditioned on P value less than 0.05. The OECloud tools at https://cloud.oebiotech.com was used to create volcano maps and box line plots, whereas the “R” package was used to create Principal component analysis (PCA) plots and Venn diagrams.

From the GSE41861 and GSE43696 databases, we first extracted 272 frequently DEGs in normal and asthma samples (Table S2). Then, using the FerrDb database, we downloaded the list of ferroptosis-related genes and intersected them with 272 co-DEGs to obtain 13 ferroptosis genes between asthma and non-asthma.

GO, KEGG and Reactome Pathway Analyses of FR-DEGs

The potential functions connected to FR-DEGs were examined using the OECloud tools at https://cloud.oebiotech.com; last accessed November 14, 2022. These analyses included GO and KEGG enrichment. The species were restricted to Homo sapiens, and statistical significance was determined by an “adjusted P value less than 0.05”. These were the three criteria for the GO terms: biological process (BP), cellular component (CC), and molecular function (MF). Reactome (https://reactome.org/; last accessed October 18, 2022) is a database of biological processes and bodily responses. Finally, the results of enrichment analysis are represented by circle diagram and bubble diagram.

Construction of mRNA-miRNA Network

The ENCORI (https://starbase.sysu.edu.cn/; last accessed October 19, 2022) database was used to predict interactions between miRNAs and differentially expressed mRNAs. Then, in order to explore the potential role between miRNA and FR-DEGs’ mRNA, a miRNA-mRNA regulatory network was established. The regulatory network was visualized with the help of Cytoscape (last accessed October 19, 2022).

Single-Gene Gene Set Enrichment Analysis (GSEA)

This analysis is implemented in the GSEA package in R (last accessed October 18, 2022). In the GSEA enrichment analysis, we used the KEGG predefined signaling pathway set to detect the enrichment pathway of each FR-DEG in the GSE43696 dataset. We calculated Pearson correlation coefficients between FR-DEGs and other genes in the dataset, which was done to explore the possible signaling pathways of FR-DEGs. Table S3 contains the particular enrichment results for each FR-DEG.

Gene Set Variation Analysis (GSVA)

GSVA were performed using the OECloud tools (https://cloud.oebiotech.com; last accessed August 29, 2022). In the GSVA²² enrichment analysis, we also used the KEGG predefined signaling pathway set. To determine whether distinct metabolic pathways are enriched between samples, GSVA performed unsupervised classification of samples by converting the expression matrix of genes between samples into the expression matrix of gene sets between samples. Simultaneously, we applied the oebiotech cloud platform (https://cloud.oebiotech.com/) to analyze the difference in GSVA score of the asthma and non-asthma. Eventually p <0.05, |t|>2 is significant, in the condition of metabolic pathways were activated in asthma was t>0, pathways were restrained when t<0. The final metabolic pathway enrichment results were collected in Table S4.

Immune Infiltration Analysis

CIBERSORT is a method that uses gene expression data to characterize the abundance of cells in mixed cell populations.²³ In this study, the Cibersort package (last accessed September 28, 2022) in R was used to predict the proportion of 22 types of immune infiltrating cells in the asthma immune microenvironment using the GSE41861 database. The final results are shown in Table S5.

Targeted Drug Prediction

This study used the DGIdb (https://www.dgidb.org/; last accessed November 2, 2022)²⁴ database to predict the target drugs for FR-DEGs and to provide directions for potential drug therapy for asthma.

Validation of Ferroptosis Genes

Finally, we validated the differential expression of 13 FR-DEGs in the GSE63142 dataset and performed receiver operating characteristic (ROC) curve analysis, expressed as the area under the ROC curve (AUC). The AUC of FR-DEGs with good diagnostic efficacy was greater than 0.6.

The Expression of FR-DEGs in Asthma Patients Was Verified

Patient population: Thirty patients with asthma (meeting the criteria for asthma diagnosis) were recruited from July 2022 to January 2023 at the Second Hospital of Shanxi Medical University, and these patients were compared with thirty healthy volunteers. The Ethics Committee of the Second Hospital of Shanxi Medical University approved the study, approval number (2022) YX (079), and received informed consent for each subject. Peripheral blood was collected and peripheral blood mononuclear cells were collected by density gradient centrifugation.

RNA isolation Quantitative reverse transcriptase-polymerase chain reaction (qRT-PCR): Total RNA from Peripheral Blood Mononuclear Cell was isolated using the RNA Isolation Kit (China, Beijing, Mei5bio MF044-01) according to the manufacturer’s protocol. RNA purity and concentration were measured using Nanodrop 2000 (ThermoFisher Scientific, USA). Next, cDNA was synthesized using the cDNA Reverse Transcription Kit (China, Nanjing, Vazyme R323-01). The specific primers were designed based on the sequence of the gene’s mRNA coding region and synthesized by Sangon Biotech. Quantitative PCR was performed using SYBR Green PCR (China, Nanjing, Vazyme Q711-02) and the details of the primers used are shown in Table 2. The relative expression of mRNA was calculated by 2−ΔΔCt method with the normalization to β-actin.

Table 2 The Primers Used in This Study

Gene	Primer-F	Primer-R
ADAMTS13	GAGGCAGAGGCAGAGGCAGAG	CCCGATGTTGAGGTTGGTGAGC
CBR1	CAGAAGTTCCGCAGTGAGACCATC	GGACAGAACGGTGACGCCAATC
ELAVL1	CATGACCCAGGATGAGTTACGAAGC	TAGTTCACAAAGCCATAGCCCAAGC
EPAS1	ACAACAACTGCCCTCCTCACAATAG	TGGCGGCTCAGGAAGGTCTTG
CREB5	TCACCAGACCTCGCCACATCC	ATCCTCGTCTACCACCCTTCGC
LRRFIP1	GGCAGCAGAAGGAGGTAGAAGAGAG	GAGAGGAAGTCCCACCCAGAGAG
MT1G	TGGACCCCAACTGCTCCTGTG	CACACTTGGCACAGCCCACAG
NOX1	CAGCAGGGGACTGGACAGAAAATC	CTCCAACCAGCACAGCCACTTC
NR1D2	TCATGCTTGCGAAGGCTGTAAGG	AGCGACATTGCTGACATCTGTTCC
VEGFA	GGTGCCCGCTGCTGTCTAATG	ACAAACAAATGCTTTCTCCGCTCTG
β-actin	CCTGGCACCCAGCACAAT	GGGCCGGACTCGTCATAC

Statistical Analysis

This study was analyzed in R. The Student’s t-test was used for comparison between the asthma and controls. DEGs were analyzed using the OECloud tools at https://cloud.oebiotech.com, and the final identified DEGs were statistically significant with p-values less than 0.05. GO and KEGG enrichment analysis and GSVA analysis were performed using the OECloud tools at https://cloud.oebiotech.com. Reactome pathway was analysed by Reactome website (https://reactome.org/). The ENCORI (https://starbase.sysu.edu.cn/) database was used to predict interactions between miRNAs and mRNAs. GSEA analysis is analyzed in the GSEA package in R. CIBERSORT analysis is performed in the Cibersort package in R. Targeted drug prediction was performed using DGIdb (https://www.dgidb.org/). Pearson correlation analysis was used for correlation analysis between FR-DEGs and immune infiltrating cells, and the AUC was used to evaluate the diagnostic efficacy of FR-DEGs in asthma. Cytoscape software was used for visualization of FR-DEGs-targeted drug networks and miRNA-miRNA networks. Statistical analysis was performed using SPSS 24.0 and GraphPad Prism 6. The measurement data were tested for normality, normally distributed data were expressed as mean ± standard deviation, and t-test for independence was used for comparison between two groups, and qualitative data were expressed as N (%); The Mann–Whitney U-test was used for comparison between two groups of non-normally distributed continuous variables. P-values less than 0.05 were statistically significant.

Results

Differentially Expressed Genes in Asthma

The “GEO2R” and “R/limma” packages were used to perform normalized and differential gene expression analysis on the GSE41861 and GSE43696 datasets that were taken from the GEO database. The results demonstrated that the boxplot distribution trends for the two datasets are essentially linear (Figure 2A and B). Additionally, PCA demonstrated good data reproducibility (Figure 2C and D).

Figure 2 The expression matrices of the GSE43696 and GSE41861datasets were normalized. (A) Boxplot of asthma and non-asthma in the GSE41861dataset; (B) Boxplot of asthma and non-asthma in the GSE43696 dataset. The green boxplot represents the asthma sample and the blue boxplot represents the control sample; (C) PCA diagrams of asthma and non-asthma in the GSE41861dataset; (D) PCA diagrams of asthma and non-asthma in the GSE43696 dataset.

Abbreviation: PCA, principal component analysis.

According to the adjusted P value<0.05, 4588 asthma-related DEGs (2325 up-regulated and 2263 down-regulated) were screened in the GSE41861 dataset, and 854 asthma-related DEGs (400 up-regulated and 454 down-regulated) were screened in the GSE43696 dataset. Volcano plots for GSE41861 and GSE43696 datasets are shown in Figure 3A and B. Then 272 common DEGs were identified in these 2 datasets, as shown in the venn plot (Figure 3C), of which 168 were down-regulated and 104 were up-regulated (Table S2). Then, we downloaded ferroptosis using the FerrDb database and intersected them with 272 co- DEGs to obtain 13 FR-DEGs between asthma and non-asthma, as shown in the venn plot (Figure 3D). Among these 13 differential genes, 11 down-regulated genes and 2 up-regulated genes were included (Table 3). Therefore, we finally identified ADAMTS13, CBR1, ELAVL1, EPAS1, CREB5, LRRFIP1, MT1G, NOX1, NR1D2, VEGFA, ACADSB, PAQR3 and STAT3 as FR-DEGs in asthma.

Table 3 13 FR-DEGs Has Been Identified as Downregulated or Upregulated in Asthma

Ensembl ID	Gene Name	Gene Type	Gene Start (bp)	Gene End (bp)	Chr no.	Status	P value	P adj
ENSG00000160323	ADAMTS13	Protein coding	133,414,337	133,459,386	9	Down	0.00107	0.04166977
ENSG00000159228	CBR1	Protein coding	36,070,024	36,073,164	21	Up	0.000286	0.02013983
ENSG00000066044	ELAVL1	Protein coding	7,958,573	8,005,641	19	Down	0.00088	0.0374346
ENSG00000116016	EPAS1	Protein coding	46,297,407	46,386,697	2	Down	0.00116	0.04364847
ENSG00000146592	CREB5	Protein coding	28,299,321	28,825,894	7	Up	0.000527	0.02799859
ENSG00000124831	LRRFIP1	Protein coding	237,527,587	237,781,643	2	Down	0.000689	0.03282089
ENSG00000125144	MT1G	Protein coding	56,666,730	56,668,065	16	Down	0.000151	0.0142877
ENSG00000007952	NOX1	Protein coding	100,843,324	100,874,359	X	Down	0.00127	0.04555435
ENSG00000174738	NR1D2	Protein coding	23,945,286	23,980,617	3	Down	0.000404	0.02426236
ENSG00000112715	VEGFA	Protein coding	43,770,209	43,786,487	6	Down	0.000897	0.03791496
ENSG00000196177	ACADSB	Protein coding	123,009,006	123,058,290	10	Down	0.00131	0.04632139
ENSG00000163291	PAQR3	Protein coding	78,887,076	78,939,438	4	Down	0.00133	0.04669752
ENSG00000168610	STAT3	Protein coding	42,313,324	42,388,442	17	Down	0.000572	0.02926981

Notes: The columns indicate the Ensembl ID, Gene name, Gene type, Gene start (bp), Gene end (bp), Chromosome number, Status=Up/Downregulated, p-value, p adj respectively. FR-DEGs were selected with p adj<0.05.

Figure 3 FR-DEGs of the GSE41861 and GSE43696 datasets. (A) The volcano plot presenting the DEGs in asthma and non-asthma in the GSE41861 dataset; (B) The volcano plot presenting the DEGs in asthma and non-asthma in the GSE43696 dataset. Red, upregulated genes; green, downregulated genes; black, no significant differences; (C) Venn diagram of DEGs, DEGs in the GSE41861 and GSE43696 datasets of asthma and non-asthma; (D) Venn diagram of FR-DEGs, the intersection of common DEGs of asthma and non-asthma in GSE41861 and GSE43696 datasets and ferroptosis database genes; (E) Heat maps of 13 FR-DEGs from the GSE41861 dataset; (F) Heat maps of 13 FR-DEGs from the GSE43696 dataset. The color of the small squares represents the size of the gene expression, with red representing high-expression genes and blue representing low-expression genes.

GO, KEGG and Reactome Pathway Enrichment Analysis of FR-DEGs

GO, KEGG and Reactome analyses were carried out to elucidate the biological functions and pathways connected to the 13 FR-DEGs. According to the findings, the GO categories was mainly enriched in surfactant homeostasis, intracellular receptor signaling pathway, monocyte differentiation, oxidoreductase activity acting on NAD(P)H, quinone or similar compound as acceptor, cellular response to vascular endothelial growth factor stimulus, angiogenesis and carbonyl reductase (NADPH) activity etc (Figure 4A and B). The KEGG pathway with the greatest enrichment were AMPK signaling pathway, pathways in cancer, HIF-1 signaling pathway, fatty acid degradation and growth hormone synthesis, etc (Figure 4C). Reactome pathway indicated that the interleukin-4 and interleukin-13 signaling, interleukin-2 family signaling and interleukin-10 signaling, regulation of gene expression by hypoxia-inducible factor, negative regulation of MAPK pathway, interleukin-7 signaling and cellular response to hypoxia were enriched (Figure 4D).

Figure 4 Enrichment functional analyses of the 13 FR-DEGs. (A) Bubble plot of GO enrichment analysis of the 13 FR-DEGs; (B) Circle diagram of GO enrichment analysis of the 13 FR-DEGs. (C) Bubble plot of KEGG enrichment analysis of the 13 FR-DEGs. (D) Bubble plot of Reactome pathway analysis of the 13 FR-DEGs.

Abbreviations: BP, biological process; CC, cell composition; MF, molecular function; KEGG, Kyoto encyclopedia of genes and genomes.

Construction of the FR-DEGs Network

We predicted the target miRNAs of FR-DEGs using the ENCORI tool. We identified 178 target miRNAs of 13 FR-DEGs. Based on the predicted 191 nodes and 370 edges, we constructed the mRNA-miRNA network by Cytoscape software (Figure 5). ADAMTS13 is regulated by 3 miRNAs, CBR1 is regulated by 11 miRNAs, ELAVL1 is regulated by 6 miRNAs, EPAS1 is regulated by 49 miRNAs, CREB5 is regulated by 31 miRNAs, LRRFIP1 is regulated by 40 miRNAs, MT1G is regulated by 7 miRNAs, NOX1 is regulated by 6 miRNAs, NR1D2 is regulated by 57 miRNAs, VEGFA is regulated by 68 miRNAs, ACADSB is regulated by 29 miRNAs, PAQR3 is regulated by 27 miRNAs and STAT3 is regulated by 36 miRNAs respectively (Table S6). Finally, using the Miranda database, mRNA-miRNA nucleic acid binding was predicted (Table S6).

Figure 5 MiRNA-mRNA regulatory network of FR-DEGs. The network included 191 nodes (178 miRNAs and 13 genes) and 370 edges.

Single‐Gene GSEA Enrichment Analysis

To more enrich and fully explore the role of these 13 FR-DEGs in asthma, we investigated the correlation between all genes and FR-DEGs in the GSE43696 dataset by single-gene GSEA enrichment analysis. Figure 6A–M shows the top ten enriched pathways for each FR-DEG.

Figure 6 Single-gene GSEA enrichment analysis in 13 FR-DEGs. (A) ACADSB, (B) ADAMTS13, (C) CBR1, (D) CREB5, (E) ELAVL1, (F) EPAS1, (G) LRRFIP1, (H) MT1G, (I) NOX1, (J) NR1D2, (K) PAQR3, (L) STAT3, (M) VEGFA.

By single-gene GSEA enrichment analysis, the results showed that these 13 FR-DEGs were mainly enriched in apoptosis, ferroptosis, lysosome, necroptosis, peroxisome, phagosome, proteasome, ribosome, amino acid metabolic pathway (“Glutathione metabolism”, “2-Oxocarboxylic acid metabolism”, “Alanine, aspartate and glutamate metabolism”, “Oxidative phosphorylation”, “Fatty acid metabolism”, “Porphyrin metabolism”, “Cysteine and methionine metabolism”, “Pyruvate metabolism”, “Citrate cycle (TCA cycle)”, “Arachidonic acid metabolism” and “Amino sugar and nucleotide sugar metabolism”), immune response (“Th1 and Th2 cell differentiation”, “Neutrophil extracellular trap formation”, “T cell receptor signaling pathway”, “B cell receptor signaling pathway” and“Th17 cell differentiation”) and various disease pathways (“Asthma”, “Acute myeloid leukemia”, “Influenza A”, “Gastric cancer”, “Tuberculosis” and “Hepatitis B and C”). In addition, most FR-DEGs were enriched in “JAK-STAT”, “Foxo”, “IL-17”, “NF-kappaB”, “Notch”, “PPAR”, “TGF-beta”, “Toll-like receptor”, “HIF-1”, “VEGF” and “Wnt” signaling pathways.

GSVA Enrichment Analysis

Based on the levels of each gene’s expression in combination with GSVA, we were able to identify the pathways that were differentially activated between the asthma and non-asthma groups (Figure 7A and B). The results showed that in the asthma group, the “FC-GAMMA-R-mediated-phagocytosis”, “Peroxisome”, “RNA- degradation”, “Cell-cycle”, “Long-term-potentiation”, “Regulation-of-autophagy”, “Amino-sugar-and-nucleotide-sugar-metabolism” and the signaling pathway (“VEGF”, “Calcium”, “T-cell-receptor” and “GNRH”) was activated, while the “ABC-transporters”, “Galactose-metabolism”, “Phenylalanine-metabolism”, “Glycerolipid-metabolism”, “Histidine-metabolism”, “Phenylalanine-metabolism”, and“Notch-signaling-pathway” was inhibited.

Figure 7 GSVA pathway analysis. High‐ and low‐expression groups in asthma-vs-controls (A) and controls -vs- asthma (B).

Immune Infiltration Analysis

Previous functional enrichment results showed that FR-DEGs were associated with immune response pathways. Studies had shown a strong connection between asthma and the immune microenvironment.^25,26 Consequently, we next used the GSE41861 database in the Cibersort package to predict the proportion of 22 immune infiltrating cells in the immune microenvironment of asthma. The results are shown in Figure 8A, the percentages of dendritic cells activated, plasma cells, T cells CD4 memory activated and mast cells activated in asthma were higher than controls, while macrophages M2, T cells regulatory (Tregs) and T cells follicular helper were lower. Moreover, Pearson analysis suggested that plasma cells were positively and negatively correlated with ADAMTS13 (p<0.01, r=0.356) and MT1G (p<0.01, r=−0.332,) respectively, Tregs had positive and negative correlations with ELAVL1 (p<0.01, r=0.374,) and CREB5 (p<0.01, r= −0.317,) respectively, macrophages M2 had positive correlations with ADAMTS13 (p<0.01, r=0.368,) and ELAVL1 (p<0.01, r= 0.349), respectively, mast cells activated had positive and negative correlations with CBR1 (p<0.01, r=−0.379,) and NR1D2 (p<0.01, r=−0.368,) respectively, NOX1 had positive correlations with T cells CD4 memory activated (p<0.01, r=0.307)(Figure 8B; Table S7).

Figure 8 Immune infiltration analysis. (A) Immune landscape analysis between asthma and non-asthma samples which implemented the CIBERSORT algorithm. (B) Pearson correlation analysis between 13 FR-DEGs and immune cells. *P<0.05, **P≤0.01, ***P≤0.001, ****P≤0.0001.

Targeted Drug Prediction

We predicted the targeted drugs of FR-DEGs using the DGIdb database (Table S8) and further demonstrated the FR-DEGs-drug network diagram using Cytoscape software (Figure 9). We investigated 71 drugs that target FR-DEGs, including 23 for STAT3, 5 for CBR1, 2 for ELAVL1, 3 for EPAS1, NOX1 targeted 2 drugs and VEGFA targeted 37 drugs. Unfortunately, we failed to predict the targeted drugs for ADAMTS13, CREB5, LRRFIP1, MT1G, NR1D2 and PAQR3. Moreover, the structural formulas of the 71 drugs listed above were retrieved using the DrugBank database. A total of 44 drug structures were retrieved. The structural information of the targeted drugs was revealed. The detailed results were illustrated in (Figure S1).

Figure 9 Prediction of gene‐targeted drugs (VEGFA, STAT3, NOX1, ELAVL1, EPAS1 and CBR1), The drugs may target genes through the DGIdb database.

GSE63142 Confirmed Differentially Expression and Diagnostic Efficacy of the FR-DEGs

Finally, we verified the expression of FR-DEGs with the GSE63142 dataset. The results showed that the expression of the 10 FR-DEGs (ADAMTS13, CBR1, CREB5, ELAVL1, EPAS1, LRRFIP1, MT1G, NOX1, NR1D2 and VEGFA) were consistent with the GSE41861 and GSE43696 dataset. The expression of CBR1 (p<0.001) and CREB5 (p≤0.001) in asthma patients was greater than controls, while the ADAMTS13 (p=0.01), ELAVL1 (p=0.003), EPAS1 (p<0.001), LRRFIP1 (p=0.015), MT1G (p<0.001), NOX1 (p<0.001), NR1D2 (p<0.001) and VEGFA (p≤0.001) were lower in asthma patients (Figure 10A–J).

Figure 10 Expression levels of the 10 FR-DEGs in asthma and non-asthma in the GSE631421 dataset. (A) ADAMTS13’ expression in the GSE631421 dataset; (B) CBR1’ expression in the GSE631421 dataset; (C) ELAVL1’ expression in the GSE631421 dataset; (D) CREB5’ expression in the GSE631421 dataset; (E) EPAS1’ expression in the GSE631421 dataset; (F) LRRFIP1’ expression in the GSE631421 dataset; (G) MT1G’ expression in the GSE631421 dataset; (H) NOX1’ expression in the GSE631421 dataset; (I) NR1D2’ expression in the GSE631421 dataset; (J) VEGFA’ expression in the GSE631421 dataset. *P<0.05, **P≤0.01, ***P≤0.001, ****P≤0.0001.

We created ROC curves for these 10 FR-DEGs in the GSE63142 database. The AUC of ADAMTS13 was 0.709 (95% CI: 0.595–0.823), CBR1 was 0.757 (95% CI: 0.652–0.863), CREB5 was 0.747 (95% CI: 0.634–0.861), ELAVL1 was 0.72 (95% CI: 0.601–0.839), EPAS1 was 0.722 (95% CI: 0.610–0.834), LRRFIP1 was 0.668 (95% CI: 0.546–0.790), MT1G was 0.851 (95% CI: 0.763–0.938), NOX1 was 0.812 (95% CI: 0.717–0.908), NR1D2 was 0.813 (95% CI: 0.717–0.910), and VEGFA was 0.729 (95% CI: 0.607–0.850) (Figure 11A–J). The results showed that these 10 FR-DEGs had good diagnostic roles in asthma.

Figure 11 Diagnostic ROC curves of 10 FR-DEGs in the GSE631421 dataset. (A) ROC curve of ADAMTS13; (B) ROC curve of CBR1; (C) ROC curve of ELAVL1; (D) ROC curve of CREB5; (E) ROC curve of EPAS1; (F) ROC curve of LRRFIP1; (G) ROC curve of MT1G; (H) ROC curve of NOX1; (I) ROC curve of NR1D2; (J) ROC curve of VEGFA.

Abbreviation: ROC, receiver operating characteristic.

Experimental Verification of Ferroptosis Genes in Asthma

Based on the asthma patients and healthy controls, the mRNA levels of 10 FR-DEGs were then measured. The general characteristics of the asthma and control are summarized in Table 4. The results showed that the mRNA levels of CBR1 and CREB5 were elevated in asthma patients compared to healthy controls, while the mRNA levels of ELAVL1, EPAS1, NOX1, NR1D2 and VEGFA were decreased. However, the expression levels of ADAMTS13, LRRFIP1, MT1G showed no significant difference between the two groups (Figure 12A–J).

Table 4 Characteristics of Asthma and Control in This Study

Characteristics	Total (n=60)	Asthma (n=30)	Control (n=30)	P-value
Age (years)	42.9±11.4	44.9±12.9	40.9±9.6	0.18
Male (n, %)	27 (45.0%)	14 (46.7%)	13 (43.3%)	-
BMI (kg/m^2)	24.4±3.0	24.2±2.9	24.7±3.1	0.58
FEV1 (L)	2.8±0.9	2.1±0.7	3.5±0.5	<0.01
FEV1/FVC (%)	79.0±16.8	64.8±8.6	93.2±9.0	<0.01
IgE (ng/mL)	55.0±49.4	87.5±51.3	22.5±12.2	<0.01

Note: Numeric and categorical variables were reported as mean ± SD and N (%), respectively.

Abbreviation: BMI, body mass index.

Figure 12 Validation of mRNA levels of 10 FR-DEGs in asthma patients and healthy controls. (A) mRNA levels of ADAMTS13; (B) mRNA levels of CBR1; (C) mRNA levels of ELAVL1; (D) mRNA levels of CREB5; (E) mRNA levels of EPAS1; (F) mRNA levels of LRRFIP1; (G) mRNA levels of MT1G; (H) mRNA levels of NOX1; (I) mRNA levels of NR1D2; (J) mRNA levels of VEGFA. **P≤0.01, ***P≤0.001, ****P≤0.0001.

Discussion

Asthma is a common and frequent disease that seriously endangers human health. In recent years, the incidence and mortality of asthma are on the rise. The main pathological changes of asthma include chronic airway inflammation, smooth muscle dysfunction, airway remodeling and mucus hypersecretion.^27,28 Therefore, asthma treatment can be achieved by modulating inflammation caused by abnormalities in the immune system,^15,29 restoring airway smooth muscle (ASM) function and reducing airway remodeling.^30,31 Although current treatments are effective in relieving clinical symptoms in patients with asthma, a large number of patients continue to suffer from worsening or treatment-related side effects. Therefore, there is a need to investigate more effective and safe treatment options,³ provided that the molecular biological process of asthma pathogenesis must be deeply explored to lay the theoretical foundation for subsequent experimental studies. The current study found that the regulation of airway inflammation in asthma is related to ferroptosis,^16,18 but studies on asthma and ferroptosis are rare. Therefore, this study explored the pathogenesis of ferroptosis in asthma using bioinformatics analysis, providing important evidence for the treatment and prevention of asthma.

Ferroptosis is a form of programmed non-apoptotic cell death caused by the accumulation of lipid peroxides and iron metabolism.⁴ Cell ferroptosis occurs when polyunsaturated fatty acids (PUFAs) are catalyzed by iron-containing lipoxygenases or divalent iron to form ROS and lipid peroxides, which accumulate in excess and eventually cause cell death.⁴ The lipid repair enzyme glutathione peroxidase 4 (GPx4) is an important protective factor of the ferroptosis process, which inhibits ferroptosis by reducing the production of lipid peroxides and ROS by attenuating lipid hydroperoxidation.^32,33 Current studies have found that ferroptosis plays a new role in the pathogenesis of airway inflammation in asthma, but the specific mechanism remains to be further studied. Tang’s study found elevated airway inflammatory indicators of asthma, ROS levels, markers of ferroptosis, and lipid peroxides in mice with house dust mite induced asthma,¹⁸ suggesting that ferroptosis plays a role in the pathology of asthma. In Zhao’s study, it was demonstrated that there is a relationship between ferroptosis and autophagy in the pathogenesis of asthma, that autophagy isolates ferritin and thus increases intracellular divalent iron levels,³⁴ and that the small scaffolding protein Phosphatidylethanolamine binding protein 1 (PEBP1) binds to iron-containing lipoxygenase and then facilitates the ferroptosis process by forming lipid peroxides. Further studies found that Ferr-1, a ferroptosis inhibitor, reduces airway inflammation in asthma by inhibiting the process of lipid peroxide production and iron release.¹⁶ Thus, it can be concluded that ferroptosis has a potential but integral role in asthma. Our study enriches the mechanism of ferroptosis in asthma through bioinformatics. The results suggest that FR-DEGs regulate asthma through oxidative stress, lysosomes, phagosomes, peroxisomes, proteasomes, autophagy, various immune response pathways, and molecular signaling pathways, but the specific mechanisms need further investigation.

The study screened out 13 FR-DEGs including ADAMTS13, CBR1, ELAVL1, EPAS1, CREB5, LRRFIP1, MT1G, NOX1, NR1D2, VEGFA, ACADSB, PAQR3 and STAT3. Among the 13 FR-DEGs s, ADAMTS13, CREB5, MT1G, and STAT3 are ferroptosis suppressor genes, while ELAVL1, EPAS1, NOX1, NR1D2, ACADSB and PAQR3 are ferroptosis driver genes. Through GO, KEGG and Reactome pathway analysis, it was finally concluded that 13 FR-DEGs were mainly enriched in oxidative stress and immune response pathways, which further confirmed that FR-DEGs may be closely associated with the pathogenesis of asthma.

The results of GSEA and GSVA analysis showed that these genes were enriched in ferroptosis, lysosomes, necrosis, autophagy, peroxisomes and phagosomes apoptosis, which were closely related to oxidative stress. Most of these genes were enriched in the immune signaling pathway (“Th1 and Th2 cell differentiation”, “Neutrophil extracellular trap formation”, “T cell receptor signaling pathway”, “B cell receptor signaling pathway” and “Th17 cell differentiation”), which is closely related to the immune response. Both oxidative stress and immune response are important pathogenesis of airway inflammation in asthma.^35,36 GSVA analysis revealed that these genes promote the ferroptosis process by activating “Regulation-of-autophagy”. Xue’s study³⁷ found that by enhancing the binding of the autophagy receptor for GPX4 protein degradation to GPX4, the autophagic degradation of GPX4 occurs, which eventually leads to ferroptosis, and the study confirmed that ferroptosis can be promoted by promoting autophagy. Zhao’s study identified a co-regulatory molecule PEBP1 of ferroptosis and autophagy, PEBP1 is also known as a rheostat between ferroptosis and autophagy, promoting ferroptosis in asthma Human airway epithelial cells (HAECs).³⁴ The current research on ferroptosis and autophagy in asthma is very limited, and more research is needed to explore. Furthermore, by GSVA analysis, we found that these specific ferroptosis genes could promote ferroptosis through activation of ‘Peroxisome’. Peroxisomes promote ferroptosis by providing substrate PUFA-ePLs for lipid peroxidation. Peroxisome gene knockout can reduce the abundance of peroxisomes, thereby reducing the inhibitory effect of GPX4 on ferroptosis, and also promote the ferroptosis process.³⁸ The “T cell receptor signaling pathway” is also involved in ferroptosis. Both in vivo and in vitro experiments, ferroptosis can promote the proliferation and differentiation of CD4 T cells and the production of cytokines by Th1 and Th17 cells through the T-cell receptor (TCR) signaling pathway, leading to inflammatory response.³⁹ GPX4 can inhibit Treg ferroptosis to maintain cellular immune balance. When GPX4 is deleted, Treg induces ferroptosis under TCR/CD28 co-stimulation.⁴⁰ No studies have yet investigated how ferroptosis plays a role in asthma through these pathways, and further studies are needed. The results of GSVA analysis could provide new ideas for the study of ferroptosis mechanism in asthma.

Immune homeostasis regulated by immune cells and various inflammatory factors have an important role in maintaining normal lung function, and these immune cells include T and B cells, neutrophils, mast cells, monocytes, eosinophils and dendritic cells (DC).¹⁵ Our results indicated that plasma cells, dendritic cells, mast cells and T cells CD4 memory were higher in asthma group, while macrophages M2, T cells follicular helper and Tregs were lower than the normal group by immune infiltration analysis. The result demonstrated that in asthma, the expression of CREB5 and CBR1 was high, and the expression of ELAVL1 and NR1D2 was low. Immune infiltration analysis showed that ELAVL1 was positively correlated with Tregs, CREB5 was negatively correlated with Tregs, NR1D2 was negatively correlated with mast cells, and CBR1 was positively correlated with mast cells. Consequently, ELAVL1, CREB5, NR1D2 and CBR1 may be latent targets for the treatment of asthma. Recently, it was discovered that ferroptosis is essential for immune homeostasis related to chronic inflammatory diseases.⁴¹ In neutrophil asthma, an excessive accumulation of ROS causes oxidative damage in addition to ferroptosis.²⁰ These 13 FR-DEGs except LRRFIP1, MT1G and VEGFA were all associated with neutrophil extracellular traps (NETs), NETs formation was upregulated in asthma. As a result, ferroptosis may also modulate neutrophil resistance to pathogens by affecting NETs formation, this can lead to the development of asthma.

In our results, it was shown that the function of Tregs might be induced by activation of ELAVL1 to alleviate asthma airway inflammation. The most well-known RNA-binding protein (RBP), ELAVL1, is crucial for the biological processes of cell cycle development, immune response, inflammation, and apoptosis.⁴² However, according to other studies, ELAVL1 increased the levels of pro-inflammatory cytokines in CD4+T cells,⁴³ and ELAVL1 inhibitors could improve the pro-inflammatory response,⁴⁴ which seems to be contrary to our results, possibly because of the different cells selected. The cAMP response element-binding (CREB) protein family includes CREB5. CREB5 plays crucial roles in apoptosis, gene regulation and cell proliferation.^45,46 Su found that Fn1 deposition in renal fibers could be regulated by targeting CREB5, which was closely related to ferroptosis.⁴⁷ Additionally, in Liu’s study, CREB5 was identified as a regulator of asthma exacerbation and could be targeted to treat asthma by inhibiting CREB5.⁴⁸ Therefore, we hypothesized that CREB5 may regulate asthma through the process of ferroptosis, which is consistent with our findings. NR1D2 is a nuclear hormone receptor associated with circadian clock disorder, which has attracted much attention in cancer research.⁴⁹ However, there are few studies on NR1D2 in asthma, and only one study on primary HAECs used to mimic asthma with circadian symptoms found that NR1D2 exhibited rhythmic expression,⁵⁰ the mechanism of which needs to be further investigated. Notably, the majority of other FR-DEGs have been documented in cancer, but no link has been found with asthma, and more researches are needed to ascertain their role in asthma.

We also used the DGIdb database to reveal drugs that might target FR-DEGs, and a total of six targeted drugs were queried (STAT3, CBR1, ELAVL1, EPAS1, NOX1, VEGFA). Among the five CBR1 target drugs retrieved, daunorubicin and doxorubicin were confirmed to be an inhibitor of CBR1. Daunorubicin is used in the treatment of leukemia, and doxorubicin is used in the treatment of acute leukemia and various cancers, but their application in asthma has not been reported. Nevertheless, a case report shows that a female patient with moderate to severe asthma who had breast cancer gradually improved her asthma symptoms after induction of chemotherapy and no longer required inhaled medication.⁵¹ Androstanolone, a steroid hormone generated by the adrenal glands, is a targeted drug for ELAVL1. Androstanolone is frequently known as the “mother hormone” since the body converts androstenedione into two steroid hormones (androgen and testosterone). It has been reported that androgen could reduce airway inflammation by inhibiting inflammatory factors produced by TH2 cells and promoting Treg cell function.^52,53 Asthma airway inflammation was inhibited by androgen receptors via the MAPK1 and MAPK14 pathways.⁵⁴ Testosterone could prevent the development of allergen-induced asthma by blocking mast cell activation.⁵⁵ Whether our predicted gene‐targeted drugs can play a role is uncertain, therefore, the selected drugs can be prospectively studied.

Finally, we verified the expression of FR-DEGs using the GSE63142 dataset. The results showed that CBR1 and CREB5 were highly expressed, while ACADSB, ELAVL1, EPAS1, LRRFIP1, MT1G, NOX1, NR1D2 and VEGFA were low expressed in asthma patients. The results of the validation dataset were consistent with the previous conclusion. We also used the GSE63142 dataset to create ROC curves. These 10 genes were of great value in the diagnosis of asthma. Except for LRRFIP1, the AUC values of the 9 FR-DEGs were all greater than 0.7, demonstrating that 10 FR-DEGs could be used as diagnostic markers and therapeutic targets for asthma. Again, we eventually experimentally validated that 7 asthma-associated ferroptosis genes were consistent with the results predicted by bioinformatics.

Limitations

Firstly, ferroptosis was poorly studied in asthma. At present, several studies on ferroptosis in asthma has involved only a few conventional markers of ferroptosis. Although our study eventually predicted and verified 7 differentially expressed genes associated with ferroptosis in asthma, there is still insufficient experimental evidence. Therefore, we will design animal and cellular experiments for further research, which will take a long period. Secondly, this study only found that the FR-DEGs may play a role in asthma by participating in oxidative stress and immune response pathways through functional analysis. Therefore, further molecular biological experiments are needed to confirm how these ferroptosis genes play a role in oxidative stress and immune response of asthma and to further verify the mechanism of action of these ferroptosis genes in the development of asthma. Finally, the diagnostic value of these ferroptosis genes in asthma also needs to be demonstrated in larger populations, furthermore, clinical trials related to targeted ferroptosis in asthma need to be conducted to provide practical insights into the potential therapeutic role of targeted ferroptosis in asthma management.

Conclusion

Ultimately, this study screened 7 FR-DEGs in asthma by bioinformatics to provide clues for the pathogenesis of ferroptosis in asthma. These FR-DEGs have good diagnostic properties in asthma and can be used as potential biomarkers specific to asthma. These results not only identify potential biomarkers and therapeutic targets specific to asthma, but also deepen the understanding of the pathogenesis of asthma. The present study provides novel insights into the mechanistic study of ferroptosis in asthma, and the next step is to experimentally demonstrate the mechanism of action of these FR-DEGs in asthma as well as the pathways to find new therapeutic avenues for asthma.

Abbreviations

FR-DEGs, Ferroptosis-related differentially expressed genes; DEGs, Differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; GSEA, Gene Set Enrichment Analysis; GSVA, Gene Set Variation Analysis; GEO, Gene Expression Omnibus; PCA, Principal component analysis; MF, Molecular function; CC, Cellular component; BP, Biological process; ROC, Receiver operating characteristic; HAECs, Human airway epithelial cells; ASM, Airway smooth muscle; PUFAs, Polyunsaturated fatty acids; ROS, Reactive oxygen species; GPx4, Glutathione peroxidase 4; RBP, RNA-binding protein; NETs, neutrophil extracellular traps; PEBP1, Phosphatidylethanolamine binding protein 1; TCR, T-cell receptor; DC, Dendritic cells.

Data Sharing Statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Ethical Approval

All procedures performed in studies involving human participants were in accordance with the ethical standards of the Ethics Committee of the second Hospital of Shanxi Medical University (approval number 2022 YX 079) and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards.

Disclosure

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We gratefully acknowledge contributions from all authors and GEO databases.

Additional information

Funding

This study was funded by key Research and Development Projects of Shanxi Province (201903D421066) and Basic Research Project of Shanxi Province (202203021211029).