Screening and Identification of Hub Genes in the Corticosteroid Resistance Network in Human Airway Epithelial Cells via Microarray Analysis

Background and Objective: Corticosteroid resistance is a major barrier to chronic obstructive pulmonary disease (COPD), but the exact mechanism of corticosteroid resistance in COPD has been less well studied. Methods: The microarray dataset GSE11906, which includes genomic and clinical data on COPD, was downloaded from the Gene Expression Omnibus (GEO) database, and the differentially expressed genes (DEGs) were identified using R software. Gene set enrichment analysis (GSEA) and Kyoto Encyclopedia of Genes (KEGG) were utilized to enrich and analyze the gene cohort related to the response to steroid hormones, respectively. The Connectivity Map (CMap) database was used to screen corticosteroid resistance-related drugs that might exert a potential therapeutic effect. STRING was used to construct a protein-protein interaction (PPI) network of the gene cohort, and the CytoHubba plug-in of Cytoscape was used to screen the hub genes in the PPI network. The expression levels of hub genes in cigarette smoke extract (CSE)-stimulated bronchial epithelial cells were assayed by quantitative real-time PCR and western blotting. Results: Twenty-one genes were found to be correlated with the response to steroid hormones. In the CMap database, 32 small-molecule compounds that might exert a therapeutic effect on corticosteroid resistance in COPD were identified. Nine hub genes were extracted from the PPI network. The expression levels of the BMP4, FOS, FN1, EGFR, and SPP1 proteins were consistent with the microarray data obtained from molecular biology experiments. Scopoletin significantly restrained the increases in the levels of AKR1C3, ALDH3A1, FN1 and reversed the decreases of phosphorylated GR and HDAC2 caused by CSE exposure. Conclusion: The BMP4, FOS, FN1, EGFR, and SPP1 genes are closely correlated with CSE-induced glucocorticoid resistance in airway epithelial cells. Scopoletin may be a potential drug for the treatment of glucocorticoid resistance caused by CSE.


INTRODUCTION
Chronic obstructive pulmonary disease (COPD) is estimated to become the third leading cause of death worldwide by 2020. Environmental exposure can cause COPD, and the primary factor is cigarette smoking (López-Campos et al., 2016). Chronic inflammation in the airways and lung parenchyma is a characteristic of COPD, and various inflammatory cytokines are upregulated in the lungs (López-Campos et al., 2016;Barnes, 2019). Although glucocorticoids are the front-line antiinflammatory therapies for chronic inflammation, most patients with COPD have a poor response (Rossios et al., 2012). The molecular mechanisms of steroid resistance are due to the decreased expression of the glucocorticoid receptor (GR), which is the receptor for corticosteroids (Barnes, 2017). Several molecular pathways are involved in steroid resistance and lead to the downregulation of GR, and these pathways include the PIK3CD/Akt, AP-1, P38 MAPK, and JNK pathways (Barnes, 2010), whose activity and expression are upregulated in the lungs of smokers or patients with COPD (Kirkham and Barnes, 2013;Barnes et al., 2015;Fischer et al., 2015). Histone deacetylase-2 (HDAC2) can be recruited by corticosteroids, which results in the effective repression of activated inflammatory genes. Cigarette smoking markedly suppresses the activity and expression of HDAC2 induced by oxidative/ nitrative stress and results in inflammation showing resistance to glucocorticoids (Barnes, 2009;Li et al., 2012). However, the effects of the current treatments for corticosteroid resistance are poor. Therefore, finding new targets for treating corticosteroid resistance has potential clinical value.
With the recent advancement of gene chips, genome-wide expression profiles can be obtained using microarray platforms and analyzed using R language and unbiased bioinformatic tools to identify novel hub genes (Sun et al., 2019). Exporting microarray data to modern pathway profiling software could reveal novel targets for treatment (Wei et al., 2019;Xie et al., 2019). Gene set enrichment analysis (GSEA) can help analyze and interpret microarray platform data (Powers et al., 2018). The power of GSEA, which summarizes microarray data of gene expression changes into gene sets, provides a list of genes that are significantly enriched based on their correlation with biological functions and biological processes. This method can enrich genes that are significantly associated with a pathway and phenotype by enabling the detection of gene sets, and its practicality for many applications has been demonstrated (Subramanian et al., 2007;Powers et al., 2018). This method has been successfully used to identify key pathways involved in pancreatic carcinoma (Zhou et al., 2018), reveal seven pathways associated with methyladenosine-related genes in hepatocellular carcinoma , and elucidate the role of the TGF-β signaling pathway in preeclampsia (Wu et al., 2019). These results show the potential use of GSEA for studying COPD-related genes and extracting genes involved in corticosteroid resistance from the whole genome.
Here, we applied a bioinformatics method to reveal novel genes associated with corticosteroid resistance in patients with COPD. The "limma" package was used to identify the differentially expressed genes (DEGs), and a GSEA was performed to reveal corticosteroid resistance-related genes. A protein-protein interaction (PPI) network was established using String (https://string-db.org/), and the hub genes were extracted using the Cytoscape software. Furthermore, the corticosteroid resistance-related genes were verified by quantitative real-time polymerase chain reaction (qPCR) and western blot (WB) analyses.

Microarray Data Collection
The raw microarray dataset of COPD (GSE11906), which was obtained using the Affymetrix human genome U133 plus 2.0 array, was downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The tissue source of these samples was the bronchial airway epithelium. The mRNA expression data were normalized by log2 transformation to obtain the COPD/normal ratio in the GSE11906 dataset.

GSEA
The "limma" package of R software (version 3.5.1) was used to screen GSE11906 to identify the DEGs, and mRNAs with a P-value < 0.05 were considered to be genes with significant differential expression. For the identification of genes associated with the response to steroid hormones, the DEGs were analyzed by a GSEA using Gene Ontology (GO) gene sets accessed from MSigDB. Java GSEA software was used for GSEA, the gene sets from COPD vs normal samples were compared, and the gene sets with an FDR <0.1 and a nominal p < 0.05 were regarded as significantly enriched. The distinct expression profiles of genes associated with the response to steroid hormones are shown as a heatmap.

Functional Enrichment Analysis
To further research the functional implications of the DEGs associated with the response to steroid hormones, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses. The results from the KEGG enrichment analyses were based on the criterion P-value < 0.05. The significant enrichment results were visualized using R software.

PPI of Corticosteroid Resistance-Related Genes and Identification of Hub Genes
After screening for genes associated with the response to steroid hormones based on module connectivity, the genes that were found to be enriched using the GSEA software were uploaded to the Search Tool for the Retrieval of Interacting Genes (STRING) database to study the interactions among these genes, as well as NR3C1/GR and HDAC2. The interaction data for the abovementioned genes were downloaded from STRING and analyzed using Cytoscape (version: 3.5.1) to construct a PPI network. Furthermore, Cytoscape software was used to analyze the PPI network to identify the hub genes (Yi et al., 2018;Zhou et al., 2018). CytoHubba in Cytoscape was used to select the hub genes in the PPIs, and the analysis was performed using the default parameters.

Screening of Potential Drugs for Corticosteroid Resistance Treatment Based on a Connectivity Map
The CMap database was used to screen small molecules with therapeutic potential in corticosteroid resistance. The DEGs associated with corticosteroid resistance were uploaded to the CMap database and analyzed, and a negative connectivity score suggested that the drug might reverse the signature biology of human cells and show its potential therapeutic value . The small molecules were identified based on the following criteria: P-value < 0.05 and a negative connectivity score.

BEAS-2B Cell Culture and Treatment
The human bronchial epithelial cell line BEAS-2B (ATCC ® CRL-9609 ™ , Shanghai Cell Bank, Chinese Academy of Sciences) was maintained in RPMI 1640 medium (Gibco, New York, NY, United States) supplemented with 10% fetal bovine serum (Gibco) at 37°C under 5% CO 2 conditions. BEAS-2B cells were treated with cigarette smoke extract (CSE) for 24 h and collected for the extraction of RNA and proteins for qPCR and WB assays, respectively.

CSE Preparation
CSE was prepared according to the protocol from our previous study (Li et al., 2012). The smoke of 10 unfiltered cigarettes was slowly dissolved in 20 ml of RPMI-1640 medium. The pH of the CSE solution was adjusted to 7.4, and the CSE solution was filtered through a 0.22-μm leach to remove bacteria. The optical density was measured at a wavelength of 320 nm. The CSE solution was used in the experiments within 1 h of extraction. A CSE concentration of 0.1% was used for the study based on the results obtained using the Cell Counting Kit-8 (CCK-8) assay.

Statistical Analysis
The qPCR and WB data are presented as the means ±SDs. Three repeated independent experiments were performed from each group. Statistical comparisons were performed using the SPSS 18.0 software (Chicago, IL, United States). A two-tailed Student's t-test was used for comparison between two groups. One-way analysis of variance was conducted to evaluate significant differences among multiple groups. The results with p < 0.05 were considered statistically significant.

Identification of Differentially Expressed Genes
Thirty-three COPD samples from long-term smokers and 72 normal samples from nonsmokers were included in this study.
The main clinical features of the participants are listed in Table 2.
A total of 335 DEGs, including 187 upregulated and 148 downregulated genes, were obtained from GSE11906 using the limma package. Volcano plots of the DEGs are shown in Figure 1 (fold change ≥1.5 and P-value ≤ 0.05).

GSEA of DEGs
GSEA of the DEGs was conducted using GSEA software based on the following filter conditions: P-value < 0.05 and FDR q-value < 0.05. The GSEA plots of the genes enriched in the response to    steroid hormones are shown in Figure 2. The 21 genes related to the response to steroid hormones that were found to be enriched by the GSEA are shown in Table 3, and the expression of these 21 genes is shown in the heatmap diagram presented in Figure 3.

KEGG Enrichment Analysis
To study the functional implications of the 21 above-mentioned genes, we performed a KEGG functional enrichment analysis. The results showed that these genes were significantly enriched in  the immune system, signal transduction, cytokine signaling in the immune system, fluid shear stress, atherosclerosis and pathways in cancer (Figure 4).

Construction of the PPI Network and Identification of Hub Genes
The 21 above-mentioned genes, as well as NR3C1 and HDAC2, were uploaded to STRING to study their interactions, and PPIs were plotted using Cytoscape based on the interaction data obtained from STRING ( Figure 5A). In the PPI network of the 21 genes related to the response to steroid hormones in COPD and NR3C1, the top nine hub genes with the highest criticality were EGFR, FOS, SPP1, FN1, CCL2, IL1B, NQ O 1, BMP4, and C3 ( Figure 5B). Therefore, we speculated that these nine hub genes might play an important role in corticosteroid resistance in smoke-induced COPD.

Verification of Hub Genes
To verify the hub genes extracted from the PPI network, we determined their mRNA and protein expression levels in CSEtreated human bronchial epithelial cells using qPCR and western blot analyses. In addition, we verified the expression of ALDH3A1 and AKR1C3, which were significantly upregulated in patients with COPD and belonged to the core enrichment genes according to the microarray analysis. The qPCR results indicated that the expression of all hub genes, as well as ALDH3A1 and AKR1C3, was consistent with the data from the selected microarray datasets ( Figure 6). The WB data showed that the expression levels of the C-FOS, EGFR, GR, P-GR, and HDAC2 proteins were significantly downregulated in the CSE group. The expression levels of the BMP4, FN1, SPP1, ALDH3A1, and AKR1C3 proteins in BEAS-2B cells increased after stimulation with 0.1% CSE (p < 0.05; Figure 7). No significant differences in the expression of the CCL2, IL1B, C3, and NQO1 proteins were found between the CSE and control groups.

Prediction and Verification of Potential Agents for the Treatment of CSE-Induced Corticosteroid Resistance
Small-molecule compounds that might alleviate corticosteroid resistance were predicted through CMap analysis. The results suggest that 32 compounds might exert a therapeutic effect on corticosteroid resistance induced by oxidative stress in COPD ( Table 4). The most promising small-molecule compounds are scopoletin, harmol, AG-012559, AH-23848, perhexiline, and fluorometholone. We chose scopoletin as a potential agent for its ability to inhibit the activation of nuclear factor-κB (NF-κB) against COPD-associated oxidative stress (Leema and Tamizhselvi, 2018). A scopoletin concentration of 100 μM was chosen as the final intervention concentration according to the results of our CCK-8 assay, and the compound was added to the medium 2 h before CSE treatment. The results indicated that scopoletin significantly restrained the increases in the levels of AKR1C3, ALDH3A1 and FN1 caused by CSE exposure. In addition, scopoletin obviously reversed the decreases in P-GR and HDAC2 levels (Figure 8).

DISCUSSION
The imbalance of pulmonary oxidants and antioxidants caused by long-term exposure to cigarette smoke is thought to be the main driving mechanism of COPD (Kirkham and Barnes, 2013). Several studies have demonstrated that oxidative stress and inflammation in COPD are not only local but also interdependent and closely related processes in systemic inflammation in the lungs (Sundar et al., 2013). Cigarette smoke-mediated oxidative stress and inflammation induce a number of kinase signaling pathways that are considered to be associated with steroid resistance (Liao et al., 2020). Bronchial/ lung epithelial cells are the first barrier against inhaled cigarette smoke. Oxidative stress induced by cigarette smoke promotes the release of various cytokines and chemokines in epithelial cells. Epithelial damage and repair play a central role in the pathogenesis of COPD (Heijink et al., 2014). Previous studies have shown that the stimulation of lung epithelial cells with cigarette smoke can result in corticosteroid insensitivity. An important molecular mechanism leading to glucocorticoid resistance has been shown to be a result of decreases in the expression and activity of GR and HDAC2 caused by cigarette smoke exposure (Barnes, 2017).
Our study extracted DEGs from the GEO database of COPD at the whole gene level, and 21 enriched genes correlated with the response to steroid hormones were identified from a GSEA of these DEGs. This gene cohort included NQ O 1, ALDH3A1, AKR1C3, NR2F1, ABHD2, NR0B1, SPP1, BMP4, CCL2, DEFB1, FN1, CA2, NR4A3, IL1B, NR2F6, FOSB, FOS, EGFR, GHR, C3, and NTS. The nine core genes in the gene cohort, which were considered to make a major contribution to the enrichment score of the gene set, were NQO1, ALDH3A1, AKR1C3, NR2F1, ABHD2, NR0B1, SPP1, BMP4, and CCL2. We also analyzed the 21-gene cohort through the KEGG database. The results showed that the genes in the cohort were enriched in the immune system, signal transduction, cytokine signaling in the immune system, fluid shear stress, atherosclerosis, and pathways in cancer. These results indicate that the functions of the 21 genes are closely correlated with the immune response. Furthermore, we constructed a PPI network to assess the relationships among the 21 genes identified in the cohort, as well as HDAC2 and GR.
A previous study demonstrated that miR-183 regulates LPSinduced oxidative stress in rat hippocampal neurons by regulating the expression of FN1 to reduce oxidative damage (Xie et al., 2020). Nicotine and cigarette smoke obviously stimulate the expression of FN1 in lung fibroblasts and mouse lung tissues (Roman et al., 2004;Zhao et al., 2019). EGFR plays an important role in cancer progression and regulates inflammation and oxidative stress (Fang et al., 2016). EGFR was shown to be activated by the meprin alpha metalloproteinase to mediate oxidative stress in macrophages (Wang et al., 2017). Additionally, EGFR suppresses the activation of transcription factors, such as NF-κB, and induces proinflammatory gene transcription in macrophages. CCL2 and IL1B are common cytokines involved in chronic inflammation in COPD (Sethi   Barnes, 2016). Cigarette smoke can induce macrophages to release inflammatory mediators, including CXCL1, CXCL8, and CCL2. The increased release of IL1 leads to an increase in the number of macrophages in the lungs of smokers and patients with COPD (Barnes, 2016). Cigarette smoke-induced reactive oxygen species (ROS) can induce IL1B generation by activating Toll-like receptors and inducing immune responses (Zuo et al., 2014). C3 activation is closely related to oxidative stress in mice. Mice with transient focal cerebral ischemia with lower ROS levels show a better neurological outcome, and the inhibition of C3 activation in mice leads to an anti-inflammatory effect (Yang et al., 2013). SPP1 is a molecule with a complex function that decreases the oxidative stress effect and has been implicated in inflammation in the kidney (Trostel et al., 2018). Additionally, SPP1 might exert proinflammatory effects in vascular endothelial cells induced by cigarette smoke and might contribute to inflammation in smokers (Bishop et al., 2012). NQO1 protein is highly expressed in various cells, such as vascular endothelial cells, epithelial cells, and adipocytes. A number of studies have shown that NQO1 exerts antioxidant effects by reducing the level of ROS (Jo et al., 2016). However, no significant upregulation was observed at the protein level in CSEtreated BEAS-2B cells. FOS is a subunit of AP-1, which is a proinflammatory transcription factor sensitive to oxidative stress that interacts with GR proteins (Frenkel et al., 2015). The overexpression of AP-1 can prevents the binding of GR to glucocorticoid response elements and other proinflammatory  (Pujolsa et al., 2009;Jacques et al., 2010). Previous studies have shown that glucocorticoid resistance in asthma is correlated with the overexpression of AP-1 (Jacques et al., 2010). BMP4 is a new proinflammatory marker that stimulates the production of ROS and induces monocyte adhesion (Tian et al., 2012). BMP4 was found to be upregulated in the airway epithelium of asymptomatic smokers and patients with COPD. The accumulation of BMP4 aggravates airway remodeling by changing the phenotypes of basal stem/progenitor cells (Zuo et al., 2019). AKR1C3 is believed to be closely related to the metabolism of steroids. It can be activated by oxidative stress factors and leads to radiotherapy resistance in lung cancer tissues (Xie et al., 2013). ALDH3A1 was found to be upregulated in lung tissue after exposure to carcinogenic aldehydes (Patel et al., 2008).
To study whether the nine hub genes are altered by CSE stimulation in BEAS-2B cells, we determined the expression of these genes using qPCR and WB assays. Our results demonstrated that the changes in the mRNA expression of the nine genes were consistent with the results from the analysis of the microarray data. The expression levels of the BMP4, FOS, FN1, EGFR, SPP1, ALDH3A1, and AKR1C3 proteins were consistent with the mRNA results. The results demonstrate that CSE could induce glucocorticoid resistance by regulating these genes. The inconsistency between the mRNA and protein levels of CCL2, IL1B, C3, and NQO1 may be related to posttranscriptional modification or post-translational protein processing, specific mechanism needs further study.
CMap is a database covering diseases, gene expression patterns, and interactions between small-molecule compounds (Gao et al., 2019;Wei et al., 2020). This database is an effective genome-based tool for screening new chemopreventive drugs and is commonly used to find small-molecule compounds that exert therapeutic effects on certain diseases (Wei et al., 2020;Yu et al., 2020). We used CMap for drug discovery and identified smallmolecule compounds as potential drugs for corticosteroid resistance treatment. Our results showed that 32 smallmolecule compounds might exert potential therapeutic effects. Some compounds have been found to exert antioxidative stress effects in other fields. For example, in yeast cells, harmol and harmalol have a significant protective effect against oxidative stress induced by H 2 O 2 and paraquat. The antigenotoxic and antimutagenic effects of harmol result from its ability to scavenge hydroxyl radicals (Moura et al., 2007). Scopoletin efficiently quenches oxidative stress by stabilizing the nuclear factor erythroid 2-related factor 2 (Nrf2)/antioxidant responsive element (ARE) pathway through enhancement of the nuclear translocation and phosphorylation of Nrf2 (Narasimhan et al., 2019). We chose scopoletin as a candidate compound that might alleviate steroid resistance based on the results of the CMap Frontiers in Pharmacology | www.frontiersin.org October 2021 | Volume 12 | Article 672065 9 database analysis. The results of treatment indicated that scopoletin significantly restrained the increases in the AKR1C3, ALDH3A1, and FN1 levels caused by CSE exposure and reversed the decreases in the P-GR and HDAC2 levels. Thus, scopoletin may contribute to the reduction of CSE-induced glucocorticoid resistance in airway epithelial cells. However, the specific molecular regulatory mechanism and how these genes interact to regulate corticosteroid resistance still need to be clarified in future experiments.
In summary, we extracted and validated seven genes that were correlated with corticosteroid resistance and found that they interacted with GR and HDAC2. However, the study has some limitations. First, we have no direct evidence showing the mutual regulatory relationships among the hub genes, and further studies are needed to validate these relationships. Second, animal model validation is needed for future mechanistic research.

CONCLUSION
In this study, we extracted 21 DEGs by integrating the microarray data of COPD and identified a set of enriched genes related to the response to steroid hormones. It was found that the BMP4, FOS, FN1, EGFR, and SPP1 genes were closely correlated with CSEinduced glucocorticoid resistance in airway epithelial cells. In addition, scopoletin was found to be a potential drug for the treatment of glucocorticoid resistance caused by CSE. We hope this research can contribute to the development of corticosteroid resistance therapy for COPD in the future.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
ZH conceived the ideas. JD designed and supervised the research. GP and NM performed the research and wrote the manuscript.
LG and FC collected and analyzed the data. JB revised the manuscript. All authors read and approved the final version of the manuscript.