Dian Chen1,2
Wenliang Wu1,2Lingling Yi1,2
Yuchen Feng1,2Chenli Chang1,2Shengchong Chen1,2Jiali Gao1,2Gongqi Chen1,2
Guohua Zhen1,2*- 1Division of Respiratory and Critical Care Medicine, Department of Internal Medicine, Tongji Hospital, Tongji Medical College, Huazhong University of Science and Technology, Wuhan, China
- 2Key Laboratory of Respiratory Diseases, National Health Commission of People’s Republic of China, and National Clinical Research Center for Respiratory Diseases, Wuhan, China
Background: Asthma is one of the most prevalent chronic respiratory diseases worldwide. Bronchial epithelial cells play a critical role in the pathogenesis of asthma. Circular RNAs (circRNAs) act as microRNA (miRNA) sponges to regulate downstream gene expression. However, the role of epithelial circRNAs in asthma remains to be investigated. This study aims to explore the potential circRNA-miRNA-messenger RNA (mRNA) regulatory network in asthma by integrated analysis of publicly available microarray datasets.
Methods: Five mRNA microarray datasets derived from bronchial brushing samples from asthma patients and control subjects were downloaded from the Gene Expression Omnibus (GEO) database. The robust rank aggregation (RRA) method was used to identify robust differentially expressed genes (DEGs) in bronchial epithelial cells between asthma patients and controls. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were used to annotate the functions of the DEGs. Protein-protein interaction (PPI) analysis was performed to identify hub genes. Three miRNA databases (Targetscan, miRDB, and miRWalk) were used to predict the miRNAs which potentially target the hub genes. A miRNA microarray dataset derived from bronchial brushings was used to validate the miRNA-mRNA relationships. Finally, a circRNA-miRNA-mRNA network was constructed via the ENCORI database.
Results: A total of 127 robust DEGs in bronchial epithelial cells between steroid-naïve asthma patients (n = 272) and healthy controls (n = 165) were identified from five mRNA microarray datasets. Enrichment analyses showed that DEGs were mainly enriched in several biological processes related to asthma, including humoral immune response, salivary secretion, and IL-17 signaling pathway. Nineteen hub genes were identified and were used to construct a potential epithelial circRNA-miRNA-mRNA network. The top 10 competing endogenous RNAs were hsa_circ_0001585, hsa_circ_0078031, hsa_circ_0000552, hsa-miR-30a-3p, hsa-miR-30d-3p, KIT, CD69, ADRA2A, BPIFA1, and GGH.
Conclusion: Our study reveals a potential role for epithelial circRNA-miRNA-mRNA network in the pathogenesis of asthma.
Introduction
Asthma is a heterogeneous disease and is characterized by chronic airway inflammation. Over 300 million people are suffering from asthma worldwide nowadays and the number is expected to reach more than 400 million by the year 2025 (To et al., 2012; Barcik et al., 2020; Global Strategy for Asthma Management and Prevention, 2020). Approximately 250, 000 confirmed deaths are reported annually worldwide due to respiratory failure during asthma exacerbations (Christiansen and Zuraw, 2019). Airway epithelial cells play a pivotal role in asthma pathogenesis, including airway inflammation, mucus overproduction, airway wall remodeling, and bronchial hyperresponsiveness (Gohy et al., 2020; Hellings and Steelant, 2020; Hammad and Lambrecht, 2021). However, the intrinsic molecular mechanisms of epithelial cells in asthma are still not fully clarified.
Circular RNAs (circRNAs), which comprise a large proportion of stable RNAs in eukaryotes, have been identified in large quantities owing to the widespread use of high-throughput RNA sequencing and the development of bioinformatics-based algorithms (Jeck et al., 2013; Memczak et al., 2013; Salzman et al., 2013; Wang et al., 2014; Gao and Zhao, 2018; Chen, 2020). CircRNAs are produced by the so-called backsplicing mechanism, a process in which a downstream 5′ donor site is covalently linked to an upstream 3′ acceptor site to form a stable closing RNA structure containing exon and/or intron sequences (Jeck et al., 2013; Memczak et al., 2013; Salzman et al., 2013; Kristensen et al., 2019; Chen, 2020). Notwithstanding a lack of 5′ and 3′ ends through the non-canonical splicing, circRNAs are generally thought to localize to the cytoplasm. Therefore, circRNAs might function as microRNA (miRNA) sponges and sequester miRNA away from mRNAs, thus indirectly regulate gene expression (Hsu and Coca-Prados, 1979; Jeck et al., 2013; Memczak et al., 2013; Chen, 2020). Such competing endogenous RNAs (ceRNA) mechanism has been considered a prominent function of circRNAs. Previous studies have revealed that circRNAs are involved in multiple biological processes including neuronal function (Piwecka et al., 2017; Kleaveland et al., 2018), pluripotency regulation (Yu et al., 2017), cell differentiation (Kristensen et al., 2018), and cancer progression (Wang et al., 2019a; Wong et al., 2020) via decoying miRNAs. MiRNAs constitute another major class of non-coding RNAs with approximately 22 nucleotides (Gebert and MacRae, 2019; Goodall and Wickramasinghe, 2021). Previous reports have demonstrated that miRNAs can bind mRNAs through partial complementarity and reduce gene expression by restraining mRNA translation and/or facilitating mRNA degradation (Yekta et al., 2004; Lewis et al., 2005; Gebert and MacRae, 2019; Goodall and Wickramasinghe, 2021). However, the role of circRNA-miRNA-mRNA networks in airway epithelial cells in asthma pathogenesis has not been reported yet.
Currently, data mining from Gene Expression Omnibus (GEO) database has been applied in asthma studies. However, the studies based on a single dataset with a limited sample size will lead to deviation and bias in data interpretation. Integrating multiple datasets to include more samples was a promising strategy to increase the credibility of the analysis. Using the robust rank aggregation (RRA) method can reduce the heterogeneity of different experiments with different analysis platforms.
In this study, we integrated five mRNA microarray datasets from the GEO database utilizing the RRA method to identify the robust DEGs in bronchial epithelial cells between asthma patients and healthy controls, based on which a potential circRNA-miRNA-mRNA regulatory network was constructed. Our study reveals a potential role for epithelial circRNA-miRNA-mRNA network in the pathogenesis of asthma.
Material and Methods
Microarray Datasets Acquisition and Processing
The mRNA and miRNA expression profiles of asthma patients were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) utilizing the getGEO function of the GEOquery package (version 2.58.0) in R software (version 4.0.3) (Barrett et al., 2012). The following terms were used to search the microarray studies: “Asthma”, AND “epithelial cells”, OR “bronchial brushings”, AND “mRNA expression”, OR “miRNA expression”, AND “Homo sapiens”, AND “Microarray”. The following eligibility criteria were used to include or exclude datasets and samples: 1) A minimum of 10 subjects in the dataset, containing at least four asthma patients and four controls; 2) Bronchial brushings were used for microarray analysis; 3) All subjects were older than 18 years old and were steroid-naïve; 4) Raw data were available in GEO.
The expression matrix file and related platform annotation document of each microarray dataset were downloaded from the GEO database to convert the names of microarray probes to the gene symbols. Probes matching with multiple gene symbols were eliminated and the mean values were calculated for gene symbols corresponding to multiple probes. After excluding samples with strong heterogeneity according to the hclust method, the differentially expressed genes (DEGs) and the differentially expressed miRNAs (DEMis) between asthma and control samples in each dataset were identified by using the “limma” (linear models for microarray data) package (version 3.46.0) (Ritchie et al., 2015) in R software with the cut-off criteria of |log2 fold change| > 0.5 and p-value < 0.05.
Robust Rank Aggregation Analysis
To integrate the outcomes of multiple microarray datasets and minimize the bias and inconsistencies, the RRA method was adopted to identify the robust DEGs. The up- and down-regulated DEGs were firstly ranked by expression fold changes in each dataset and were subsequently analyzed using the “Robust Rank Aggregation” R package (version 1.1) (Kolde et al., 2012). The score in the RRA result indicated the ranking degree of each gene in the final gene list. Genes with score < 0.05 and |log2 fold change| > 0.5 were considered as the significant robust DEGs.
Visualization of Chromosome Locations of Robust Differentially Expressed Genes
The “RCircos” R package (1.2.1) was utilized to visualize the expression patterns and the chromosomal positions of all robust DEGs.
Functional and Pathway Enrichment Analyses
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were used to investigate the biological process (BP), the cellular component (CC), the molecular function (MF), and the involved pathways of selected molecules, which were performed with the “clusterprofiler” R package (version 3.18.1) (Yu et al., 2012). The GO terms and KEGG pathways with p-value < 0.05 were considered statistically significant and further visualized via the “ggplot2” R package (version 3.3.3). The KEGG enrichment analysis of differentially expressed miRNAs was conducted using the miRPathDB v2.0 (https://mpd.bioinf.uni-sb.de/) database (Kehl et al., 2020).
DisGeNET Database Analysis
DisGeNET (http://www.disgenet.org) database is one of the largest available platforms of human disease-associated genes and variants through manually integrating the scientific literature (Piñero et al., 2017). For a given gene list, DisGeNET database can identify significantly correlated diseases.
Protein-Protein Interaction Network Construction and Clusters Analysis
All previously identified robust DEGs were uploaded to the STRING (version 11.0) database (https://www.string-db.org/) to construct the protein-protein interaction (PPI) network (Szklarczyk et al., 2021). Confidence > 0.4 was set as the screening criteria. The PPI network was subsequently reconstructed and visualized via the Cytoscape (version 3.8.2) (http://cytoscape.org/) software (Su et al., 2014). In the Cytoscape plot, each node represented a gene/protein/miRNA/circRNA, while the edge between nodes represented the interactions of molecules. The molecular complex detection (MCODE) plugin of the Cytoscape software was used to screen out significant clusters in the PPI network.
Hub Gene Identification
CytoHubba is another plugin of the Cytoscape software for ranking nodes in a network, which provides eleven topological analysis methods and six centralities to identify hub genes based on shortest paths, including Maximal Clique Centrality (MCC), Density of Maximum Neighborhood Component (DMNC), Maximum Neighborhood Component (MNC), Degree, Edge Percolated Component (EPC), BottleNeck, EcCentricity, Closeness, Radiality, and Betweenness (Chin et al., 2014). Considering the potential heterogeneity of the biological network, the RRA method was adopted to integrate the results of different analysis methods for catching essential proteins.
GeneMANIA Database Analysis
GeneMANIA (http://www.genemania.org) database was used to construct the PPI network and explore the putative functions of up-loaded genes (Warde-Farley et al., 2010). For a given query list, GeneMANIA analyzes target genes with functionally similar genes together to obtain regulatory networks.
CircRNA-miRNA-mRNA Network Construction
The Targetscan (http://www.targetscan.org/vert_72/) (Garcia et al., 2011), miRDB (http://mirdb.org/) (Chen and Wang, 2020), and miRWalk (http://mirwalk.umm.uni-heidelberg.de/) (Dweep et al., 2011) databases were used to predict the corresponding miRNAs of all 19 hub genes. The overlapping results of three databases were then intersected with the differentially expressed miRNAs of GSE142237. The ENCORI (http://starbase.sysu.edu.cn/index.php) database (Li et al., 2014) was used to predict the upstream circRNAs of the selected miRNA-mRNA pairs. The final ceRNA network was further processed with the Cytoscape software (version 3.8.2).
Statistical Analysis
The differential analysis was conducted by the “limma” package (version 3.46.0) in R version 4.0.3. Heatmap was used to reveal the logarithmic fold changes of robust DEGs in the RRA analysis. p < 0.05 was considered statistically significant.
Results
Subjects Characteristics of the Microarray Datasets Included in This Study
Five mRNA microarray datasets (GSE4302, GSE43696, GSE63142, GSE67472, and GSE41861) and one miRNA microarray dataset (GSE142237) derived from bronchial epithelial brushings were obtained from the GEO database. There were a total of 272 steroid-naïve asthma patients and 165 healthy controls in the five mRNA microarray datasets. The miRNA microarray dataset (GSE142237) included a total of eight asthma patients and four healthy controls. Only asthma patients without any steroid treatments were included for further analysis.
The workflow of the study was shown in Figure 1. Detailed information on the datasets mentioned above was shown in Table 1.
FIGURE 1. The whole study workflow. GEO, Gene Expression Omnibus; DEGs, differentially expressed genes; RRA, robust rank aggregation; PPI, protein-protein interaction.
TABLE 1. Characteristics of six microarray datasets included in the study.
Identification of Differentially Expressed Genes in Steroid-Naïve Asthma Patients
After base two logarithm conversions, all five mRNA microarray datasets were further analyzed using the “limma” package of the R software to identify DEGs in primary airway epithelial cells between steroid-naïve asthma patients and healthy controls. The cut-off criteria were as follows: p-value < 0.05 and |log2 fold change| > 0.5. The distribution of DEGs in five datasets was illustrated by the volcano plots in Figure 2. There were 851 DEGs (737 upregulated and 114 downregulated), 1218 DEGs (606 upregulated and 612 downregulated), 783 DEGs (345 upregulated and 438 downregulated), 1084 DEGs (650 upregulated and 434 downregulated), and 875 DEGs (348 upregulated and 527 downregulated) in GSE4302, GSE43696, GSE63142, GSE67472, and GSE41861, respectively.
FIGURE 2. Volcano plots of five mRNA microarray datasets. The upregulated genes were marked in red, while the downregulated genes were marked in blue. The gray dots represented genes with no significant difference. (A) GSE4302; (B) GSE43696; (C) GSE63142; (D) GSE67472; (E) GSE41861.
The RRA analysis was conducted to integrate five mRNA microarray datasets with minimal bias and a total of 127 integrated DEGs were identified, including 83 upregulated and 44 downregulated robust DEGs. Lower RRA scores in the results indicated higher gene ranks and the top 20 genes were shown in the lollipop chart (Figure 3A), only five of which were downregulated. The top 10 ranked genes in asthma were composed of eight upregulated genes [including CLCA1 (score = 4.99E-15), CPA3 (score = 1.21E-12), CST1 (score = 1.56E-11), PRR4 (score = 2.51E-11), SERPINB2 (score = 1.08E-08), POSTN (score = 3.57E-08), ITLN1 (score = 7.85E-08), UPK1B (score = 1.11E-07)] and two downregulated genes [including BPIFA1 (score = 3.52E-12), FHOD3 (score = 3.90E-08)]. Of note, CLCA1, SERPINB2, and POSTN, as the IL-13-responsive genes and airway epithelial gene signatures for type 2 status, were all at a higher rank in five microarray datasets. This suggested the significant roles of the three genes in asthma pathogenesis and the reliability of the RRA results (Woodruff et al., 2007; Woodruff et al., 2009). The top 10 up- and down-regulated DEGs were drawn on the heatmap (Figure 3B). The expression patterns and the chromosomal positions of all robust DEGs across the five microarray datasets were shown in Figure 3C. These DEGs were distributed in all chromosomes except chromosome Y. Chromosome three contained the most DEGs.
FIGURE 3. Identification of robust DEGs in the RRA analysis. (A) The lollipop chart showed the top 20 ranked genes identified by the RRA method. The red and blue dots represented the up- and down-regulated genes, respectively. Bigger dots represented higher ranks; (B) The heatmap showed the top 10 up- and down-regulated DEGs in five GEO series. Red represented high expression of DEGs in asthma patients, while blue represented low expression of DEGs in asthma patients. The numbers in the box indicated logarithmic fold changes in each dataset; (C) The circular heatmaps showed the chromosomal positions of all robust DEGs. The outer circle represented chromosomes, while the inner circle heatmaps represented logarithmic fold changes of all robust DEGs in five asthma microarray datasets.
Functional Annotation of Robust Differentially Expressed Genes in Asthma
GO annotation and KEGG pathway enrichment analyses of robust DEGs were performed to explore the biological classifications in asthma utilizing the clusterProfiler R package (Figures 4A–C). Three categories including biological process, cellular component, and molecular function were analyzed within up- and down-regulated DEGs, respectively. For 83 upregulated genes, the outcomes showed that changes in the biological process were mainly enriched in regulation of endopeptidase and peptidase activities, followed by response to tissue homeostasis, negative regulation of proteolysis. The significantly enriched entries for the cellular component part were collagen-containing extracellular matrix, cornified envelope, and external side of plasma membrane. Furthermore, endopeptidase and peptidase regulator activities, enzyme inhibitor activity, and cysteine-type endopeptidase inhibitor activity accounted for the majority of the molecular function terms (Figure 4A). In terms of 44 downregulated genes, the significantly enriched biological process terms were humoral immune response, response to drug, and pattern specification process. In the cellular component part, the downregulated genes were particularly enriched in tight junction, brush border membrane, and Z disc. Meanwhile, endopeptidase and peptidase regulator activities, enzyme inhibitor activity, and heme binding were mainly enriched in the molecular function group (Figure 4B). In addition, integrated DEGs were mainly involved in salivary secretion, metabolism of xenobiotics by cytochrome P450, IL-17 signaling pathway, and leukocyte transendothelial migration in KEGG pathway analysis (Figure 4C). The DisGeNET database was further used to identify DEGs related diseases. As shown in Figure 4D, the result indicated that robust DEGs participated in the progression of various diseases, such as Nasal Polyps, Allergic rhinitis disorder, Allergic asthma, and Atopic Dermatitis, which were all related to allergic reactions and chronic inflammation (Figure 4D). Taken together, the above results indicated that the robust DEGs were mostly associated with asthma-related functions.
FIGURE 4. Bar plots and bubble charts of functional annotations involved in asthma. GO enrichment annotations of upregulated DEGs (A) and downregulated DEGs (B) in three categories: BP, CC, and MF; (C) KEGG pathway enrichment analysis of all DEGs; (D) Enrichment analysis of all DEGs in DisGeNET database. GO, Gene Ontology; BP, biological process; CC, cellular component; MF, molecular function; KEGG, the Kyoto Encyclopedia of Genes and Genomes.
Protein-Protein Interaction Network Construction, Clusters Analysis, and Hub Gene Identification
In order to explore the potential protein-protein interactions in asthma, all 127 robust DEGs were uploaded to the STRING database for further analysis (http://string.embl.de/). After hiding the disconnected nodes, the Cytoscape software was adopted to visualize the network (Figure 5A). As shown in the final network, 77 nodes and 114 edges were obtained, including 57 upregulated and 20 downregulated genes. Three key clusters were identified from the whole network using the MCODE plugin (Figures 5B–D). GO enrichment analyses showed that the significantly enriched biological process terms of three clusters were regulation of myeloid leukocyte mediated immunity, T cell activation, and antibacterial humoral response, respectively (Figure 5E).
FIGURE 5. Protein-protein interaction (PPI) network construction, key clusters analyses, and hub genes identification. (A) The whole PPI network with all robust DEGs; (B) PPI network of Cluster 1; (C) PPI network of Cluster 2; (D) PPI network of Cluster 3. Red circles represented upregulated DEGs, while blue circles represented downregulated DEGs. (E) The significantly enriched entries for the biological process of three clusters; (F) The lollipop chart showed all hub genes identified by the RRA method. The red and blue dots represented the up- and down-regulated hub genes, respectively. Bigger dots represented higher ranks.
Hub genes were subsequently screened out utilizing the cytoHubba plugin, which investigates the most important nodes in the PPI network with several topological analysis algorithms. In order to improve the positive rate of hub gene identification, the RRA method was used to integrate the top 50 ranked genes calculated by ten different topological analysis algorithms (MCC, DMNC, MNC, Degree, EPC, BottleNeck, EcCentricity, Closeness, Radiality, and Betweenness) and a total of 19 genes were obtained accordingly (Figure 5F). The detailed information of all 19 genes was listed in Table 2, including full names, scores of the RRA method, direction, and primary functions. The upset diagram of the top 50 ranked genes from the ten algorithms was shown in Supplementary Figure 1. The GeneMANIA database was used to construct the regulatory network among these 19 hub genes with functionally similar genes and the result showed that these genes may be involved in the following functions: tissue homeostasis, secretory granule, serine-type peptidase and endopeptidase activities, serine hydrolase activity, regulation of protein processing, and humoral immune response (Figure 6).
TABLE 2. The detailed information of 19 hub genes.
FIGURE 6. Protein-protein interaction network of the 19 hub-genes target network. Protein-protein interaction network of the 19 hub-genes target network was constructed using GeneMANIA database. The colors of the edges in the network indicated different bioinformatics methods used, including co-expression, website prediction, shared protein domains, and co-localization. The colors of the nodes in the network indicated the functional enrichment analysis of the query gene list.
circRNA-miRNA-mRNA Network Construction
All 19 hub genes, including 14 upregulated genes (CST4, CTSG, CLCA1, CSTA, CPA3, KIT, SERPINB2, GGH, MUC5AC, POSTN, ADRA2A, TPSAB1, CD69, and AZGP1) and five downregulated genes (LTF, C3, MUC5B, BPIFA1, and SCGB1A1), were used for the circRNA-miRNA-mRNA network construction. MiRNA-mRNA analyses were performed with Targetscan, miRDB, and miRWalk databases. The intersection of miRNA results predicted by three databases was selected as the prediction result. Furthermore, GSE142237 was adopted to validate the prediction result. The volcano plot showed the distribution of DEMis of GSE142237 and there were 522 DEMis (184 upregulated and 338 downregulated) in asthma (Figure 7A). Generally, miRNA has a negative regulatory relationship with its targeted mRNA (Yekta et al., 2004; Lewis et al., 2005; Gebert and MacRae, 2019; Goodall and Wickramasinghe, 2021). Therefore, the miRNAs targeting upregulated hub genes further intersected with downregulated DEMis of GSE142237, while the miRNAs targeting downregulated hub genes intersected with upregulated DEMis. As shown in the Venn diagrams in Figures 7B,C and Table 3, there were 45 miRNA-mRNA relationships in total. The KEGG analyses of seven upregulated and 34 downregulated miRNAs were further conducted via the miRPathDB database (Supplementary Figure 2).
FIGURE 7. CircRNA-miRNA-mRNA network construction. (A) The volcano plot of differentially expressed miRNAs of GSE142237. The upregulated miRNAs were marked in red, while the downregulated miRNAs were marked in blue. The gray dots represented miRNAs with no significant difference. (B) The Venn diagram showed the intersection between the miRNAs targeting upregulated hub genes in the prediction result (red) and the downregulated miRNAs of GSE142237 (blue). (C) The Venn diagram showed the intersection between the miRNAs targeting downregulated hub genes in the prediction result (blue) and the upregulated miRNAs of GSE142237 (red). (D) The circRNA-miRNA-mRNA network. CircRNAs were marked as octagons, miRNAs were marked as ellipses, and mRNAs were marked as diamonds. The size of the nodes indicated the degree of the connection. Up-regulated molecules were marked in red, while down-regulated molecules were marked in blue. (E) The top 10 ranked ceRNAs identified by MCC algorithms. CircRNAs, miRNAs, and mRNAs were marked as octagons, ellipses, and diamonds, respectively. The colors of nodes represented the degree of the connection.
TABLE 3. 45 miRNA-mRNA relationships.
Previous reports have demonstrated that circRNA could function as miRNA sponge to prevent co-expressed mRNA from miRNA-mediated degradation (Piwecka et al., 2017; Yu et al., 2017; Kleaveland et al., 2018; Kristensen et al., 2018; Wang et al., 2019a; Wong et al., 2020). Therefore, the target mRNA has the same expression pattern as circRNA. Consequently, 45 miRNA-mRNA relationships were used to construct the potential circRNA-miRNA-mRNA network. The corresponding circRNAs were predicted with the ENCORI database, while hsa-miR-1293-CST4 and hsa-miR-3664-5p-CTSG failed to obtain targeted circRNAs. After cross-linking, only circRNAs that were able to regulate all remaining mRNAs by targeting corresponding miRNAs were retained to build the network. As was shown in Figure 7D, a ceRNA network, with three circRNAs, 27 miRNAs, and 12 mRNAs, was finally constructed. The detailed information of three circRNAs was listed in Supplementary Table 1. The top 10 ranked ceRNAs identified by MCC algorithm were shown in Figure 7E, including hsa_circ_0001585, hsa_circ_0078031, hsa_circ_0000552, hsa-miR-30a-3p, hsa-miR-30d-3p, KIT, CD69, ADRA2A, BPIFA1, and GGH.
Discussion
In the present study, five microarray datasets from the GEO database were downloaded and analyzed by the RRA method to identify robust DEGs in primary epithelial cells between asthma patients and healthy controls. The RRA method is a both computationally efficient and statistically stable algorithm to combine multiple lists of genes from several datasets, which has four key features: robustness to noise, ability to handle incomplete rankings, assigning a significant score to each element in the resulting ranking, and efficient to compute (Kolde et al., 2012). However, the RRA method has not been reported to be used in asthma research yet. Using the RRA method, our study systematically integrated multiple microarray datasets on asthma in GEO. Furthermore, GO and KEGG enrichment analyses were conducted to annotate the functions of robust DEGs. In addition, the PPI network was constructed and key clusters were selected. In order to identify hub genes, the RRA method was utilized again to integrate the results of ten cytohubba plugin algorithms and nineteen genes were obtained. Potential miRNA-mRNA pairs were predicted by three miRNA databases (Targetscan, miRDB, and miRWalk) and further validated by a miRNA microarray dataset (GSE142237) to increase the reliability. By using the ENCORI database, a circRNA-miRNA-mRNA regulatory network was finally constructed.
The final ceRNA network included three circRNAs, 27 miRNAs, and 12 mRNAs. KIT, CD69, ADRA2A, BPIFA1, and GGH were subsequently identified as hub genes using the MCC algorithm. Of note, BPIFA1 was among the top 10 ranked genes, while KIT, CD69, ADRA2A, and GGH ranked the 18th, the 20th, the 28th, and the 64th, respectively. Stem cell factor and its receptor, the KIT proto-oncogene receptor tyrosine kinase (henceforth known as KIT), is involved in mast cell development, migration, and function (Silva et al., 2006). Finotto and others found that the ligand of KIT, stem cell factor (SCF), played a critical role in a murine asthma model. Suppressing SCF expression in epithelial cells decreased various signs of lung inflammation (Finotto et al., 2001). In this study, KIT was also found to be significantly upregulated in bronchial epithelial cells. CD69 is a type II transmembrane receptor, an activation marker of eosinophils. Kwon et al. reported that oleoylethanolamide increased CD69 expression on purified eosinophils, thus playing a role in the pathogenesis of asthma by inducing eosinophilic airway inflammation (Kwon et al., 2021). Adrenoceptor Alpha 2A (ADRA2A) mediates the catecholamine-induced inhibition of adenylate cyclase through the action of G proteins. Yoshie et al. found that alpha-2 adrenoceptors existed in human airways and the overfunction of these receptors could cause intractable asthma (Yoshie et al., 1988). Bacterial permeability family member A1 (BPIFA1) is abundantly expressed in normal airway surface liquid and involved in the anti-inflammatory response. Thaikoottathil et al. found that BPIFA1 inhibited airway eosinophilic inflammation by reducing eotaxin-2 production in alveolar macrophages (Thaikoottathil et al., 2012), which was consistent with Schaefer’s research (Schaefer et al., 2019). γ-glutamyl-hydrolase (GGH) is a ubiquitously expressed enzyme that regulates cell proliferation, DNA synthesis, and repair. However, the relationship between GGH and asthma has not yet been characterized, which requires further investigation. Numerous studies have concentrated on the diagnostic functions and therapeutic targets of these regulatory molecules for patients with asthma. Cahill et al. reported that both airway hyperresponsiveness and mast cell counts were decreased in patients with severe asthma after treated with imatinib, a KIT inhibitor (Cahill et al., 2017). It was also reported that anti-CD96 mAb treatment could inhibit established airway inflammation as effectively as dexamethasone pretreatment in a mouse model of asthma (Wang et al., 2015). Sakai et al. found that the antagonist of ADRA2A might participate in the inhibition of the allergen provoked late asthmatic response (Sakai et al., 1995). However, there were no reports, so far, on the roles of BPIFA1, GGH, hsa-miR-30a-3p, hsa-miR-30d-3p, hsa_circ_0001585, hsa_circ_0078031, and hsa_circ_0000552 in asthma.
Extensive studies have revealed that miRNAs expressed in bronchial epithelial cells contribute to asthmatic pathogenesis. In this study, hsa-miR-30d-3p and hsa-miR-30a-3p were identified in the final top 10 ceRNAs. Previous bioinformatic analyses showed that hsa-miR-30d-3p was associated with non-small cell lung cancer and inhibited epidermal growth factor receptor-targeted medicine therapy (Wang et al., 2017; Pan et al., 2019). Hsa-miR-30d-3p has also been implicated as a novel biomarker for treatment monitoring of postmenopausal osteoporosis (Weigl et al., 2021) and cerebral ischemia-reperfusion injury (Jin et al., 2021). However, hsa-miR-30d-3p has not been reported in asthma yet. Hsa-miR-30a-3p was reported to regulate the tumorigenesis in various cancer, such as gastric cancer (Wang et al., 2019b), lung adenocarcinoma (Wang et al., 2020), and pancreatic ductal adenocarcinoma (Shimomura et al., 2020). Li and others reported that hsa-miR-30a-3p regulates eosinophil activity through targeting CCR3 in asthma (Li et al., 2020).
Hsa_circ_0001585, hsa_circ_0078031, and hsa_circ_0000552 were the three circRNAs finally identified in the ceRNA network. So far, there were no reports on these circRNAs. As shown in Supplementary Table 1, these circRNAs were mostly intergenic circRNAs with relatively long spliced lengths, bringing obstacles for research. However, the exceeding spliced length of these circRNAs could provide numerous miRNA response elements for miRNA to bind.
In summary, we included six microarray datasets (five mRNA datasets and one miRNA dataset) and utilized the RRA method to obtain robust DEGs and robust hub genes. Based on the prediction of three miRNA-related databases (Targetscan, miRDB, and miRWalk) and one circRNA-related database (ENCORI), an epithelial circRNA-miRNA-mRNA network was finally constructed and the top 10 epithelial ceRNAs were identified. This epithelial ceRNA network provides new clues for future study on airway epithelial cells in asthma.
Data Availability Statement
Publicly available datasets were analyzed in this study. This data can be found here: GEO database (https://www.ncbi.nlm.nih.gov/geo/).
Author Contributions
DC and GZ conceptualized the study design. DC, WW, LY, YF, and CC acquired data. DC, WW, and SC analyzed the data. DC, JG, and GC visualized the data. DC and GZ interpreted the results. DC and GZ wrote the manuscript with all authors providing feedback for revision. All authors read and approved the final report.
Funding
Supported by National Natural Science Foundation of China (grants 81670019, 91742108, 81800026), National Key Research and Development Program of China (2016YFC1304400).
Conflict of Interest
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 sincerely thank the scientists who shared their data on the public database. We thank Rong-bing Yi, Wei-feng Hong, Yu-zhou Xue, Lu Wang, Bo Liang for their suggestions and support.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2021.703307/full#supplementary-material
References
Barcik, W., Boutin, R. C. T., Sokolowska, M., and Finlay, B. B. (2020). The Role of Lung and Gut Microbiota in the Pathology of Asthma. Immunity 52, 241–255. doi:10.1016/j.immuni.2020.01.007
Barrett, T., Wilhite, S. E., Ledoux, P., Evangelista, C., Kim, I. F., Tomashevsky, M., et al. (2012). NCBI GEO: Archive for Functional Genomics Data Sets-Update. Nucleic Acids Res. 41, D991–D995. doi:10.1093/nar/gks1193
Cahill, K. N., Katz, H. R., Cui, J., Lai, J., Kazani, S., Crosby-Thompson, A., et al. (2017). KIT Inhibition by Imatinib in Patients With Severe Refractory Asthma. N. Engl. J. Med. 376, 1911–1920. doi:10.1056/nejmoa1613125
Chen, L.-L. (2020). The Expanding Regulatory Mechanisms and Cellular Functions of Circular RNAs. Nat. Rev. Mol. Cel Biol. 21, 475–490. doi:10.1038/s41580-020-0243-y
Chen, Y., and Wang, X. (2020). miRDB: an Online Database for Prediction of Functional microRNA Targets. Nucleic Acids Res. 48, D127–D131. doi:10.1093/nar/gkz757
Chin, C.-H., Chen, S.-H., Wu, H.-H., Ho, C.-W., Ko, M.-T., and Lin, C.-Y. (2014). CytoHubba: Identifying Hub Objects and Sub-Networks From Complex Interactome. BMC Syst. Biol. 8, S11. doi:10.1186/1752-0509-8-s4-s11
Christiansen, S. C., and Zuraw, B. L. (2019). Treatment of Hypertension in Patients With Asthma. N. Engl. J. Med. 381, 1046–1057. doi:10.1056/nejmra1800345
Dweep, H., Sticht, C., Pandey, P., and Gretz, N. (2011). miRWalk - Database: Prediction of Possible miRNA Binding Sites by "Walking" the Genes of Three Genomes. J. Biomed. Inform. 44, 839–847. doi:10.1016/j.jbi.2011.05.002
Finotto, S., Buerke, M., Lingnau, K., Schmitt, E., Galle, P. R., and Neurath, M. F. (2001). Local Administration of Antisense Phosphorothioate Oligonucleotides to the C-Kit Ligand, Stem Cell Factor, Suppresses Airway Inflammation and IL-4 Production in a Murine Model of Asthma. J. Allergy Clin. Immunol. 107, 279–286. doi:10.1067/mai.2001.113049
Gao, Y., and Zhao, F. (2018). Computational Strategies for Exploring Circular RNAs. Trends Genet. 34, 389–400. doi:10.1016/j.tig.2017.12.016
Garcia, D. M., Baek, D., Shin, C., Bell, G. W., Grimson, A., and Bartel, D. P. (2011). Weak Seed-Pairing Stability and High Target-Site Abundance Decrease the Proficiency of Lsy-6 and Other microRNAs. Nat. Struct. Mol. Biol. 18, 1139–1146. doi:10.1038/nsmb.2115
Gebert, L. F. R., and MacRae, I. J. (2019). Regulation of microRNA Function in Animals. Nat. Rev. Mol. Cel Biol. 20, 21–37. doi:10.1038/s41580-018-0045-7
Global Strategy for Asthma Management and Prevention (2020). Available at: https://ginasthma.org/wp-content/uploads/2020/06/GINA-2020-report_20_06_04-1-wms.pdf Accessed 2 April 2021.
Gohy, S., Hupin, C., Ladjemi, M. Z., Hox, V., and Pilette, C. (2020). Key Role of the Epithelium in Chronic Upper Airways Diseases. Clin. Exp. Allergy 50, 135–146. doi:10.1111/cea.13539
Goodall, G. J., and Wickramasinghe, V. O. (2021). RNA in Cancer. Nat. Rev. Cancer 21, 22–36. doi:10.1038/s41568-020-00306-0
Hammad, H., and Lambrecht, B. N. (2021). The Basic Immunology of Asthma. Cell 184, 1469–1485. doi:10.1016/j.cell.2021.02.016
Hellings, P. W., and Steelant, B. (2020). Epithelial Barriers in Allergy and Asthma. J. Allergy Clin. Immunol. 145, 1499–1509. doi:10.1016/j.jaci.2020.04.010
Hsu, M.-T., and Coca-Prados, M. (1979). Electron Microscopic Evidence for the Circular Form of RNA in the Cytoplasm of Eukaryotic Cells. Nature 280, 339–340. doi:10.1038/280339a0
Jeck, W. R., Sorrentino, J. A., Wang, K., Slevin, M. K., Burd, C. E., Liu, J., et al. (2013). Circular RNAs are Abundant, Conserved, and Associated With ALU Repeats. RNA 19, 141–157. doi:10.1261/rna.035667.112
Jin, D., Wei, W., Song, C., Han, P., and Leng, X. (2021). Knockdown EZH2 Attenuates Cerebral Ischemia-Reperfusion Injury via Regulating microRNA-30d-3p Methylation and USP22. Brain Res. Bull. 169, 25–34. doi:10.1016/j.brainresbull.2020.12.019
Kehl, T., Kern, F., Backes, C., Fehlmann, T., Stöckel, D., Meese, E., et al. (2020). miRPathDB 2.0: a Novel Release of the miRNA Pathway Dictionary Database. Nucleic Acids Res. 48, D142–D147. doi:10.1093/nar/gkz1022
Kleaveland, B., Shi, C. Y., Stefano, J., and Bartel, D. P. (2018). A Network of Noncoding Regulatory RNAs Acts in the Mammalian Brain. Cell 174, 350–362. doi:10.1016/j.cell.2018.05.022
Kolde, R., Laur, S., Adler, P., and Vilo, J. (2012). Robust Rank Aggregation for Gene List Integration and Meta-Analysis. Bioinformatics 28, 573–580. doi:10.1093/bioinformatics/btr709
Kristensen, L. S., Andersen, M. S., Stagsted, L. V. W., Ebbesen, K. K., Hansen, T. B., and Kjems, J. (2019). The Biogenesis, Biology and Characterization of Circular RNAs. Nat. Rev. Genet. 20, 675–691. doi:10.1038/s41576-019-0158-7
Kristensen, L. S., Okholm, T. L. H., Venø, M. T., and Kjems, J. (2018). Circular RNAs are Abundantly Expressed and Upregulated During Human Epidermal Stem Cell Differentiation. RNA Biol. 15, 280–291. doi:10.1080/15476286.2017.1409931
Kwon, E. K., Choi, Y., Yoon, I. H., Won, H. K., Sim, S., Lee, H. R., et al. (2021). Oleoylethanolamide Induces Eosinophilic Airway Inflammation in Bronchial Asthma. Exp. Mol. Med. 53, 1036–1045. doi:10.1038/s12276-021-00622-x
Lewis, B. P., Burge, C. B., and Bartel, D. P. (2005). Conserved Seed Pairing, Often Flanked by Adenosines, Indicates that Thousands of Human Genes are microRNA Targets. Cell 120, 15–20. doi:10.1016/j.cell.2004.12.035
Li, J.-H., Liu, S., Zhou, H., Qu, L.-H., and Yang, J.-H. (2014). starBase v2.0: Decoding miRNA-ceRNA, miRNA-ncRNA and Protein-RNA Interaction Networks From Large-Scale CLIP-Seq Data. Nucl. Acids Res. 42, D92–D97. doi:10.1093/nar/gkt1248
Li, X., Wang, B., Huang, M., and Wang, X. (2020). miR-30a-3p Participates in the Development of Asthma by Targeting CCR3. Open Med. (Wars) 15, 483–491. doi:10.1515/med-2020-0102
Memczak, S., Jens, M., Elefsinioti, A., Torti, F., Krueger, J., Rybak, A., et al. (2013). Circular RNAs are a Large Class of Animal RNAs With Regulatory Potency. Nature 495, 333–338. doi:10.1038/nature11928
Pan, L., Chen, Y., Ying, G., Wang, H., Jiang, C., and Li, W. (2019). Multiple MicroRNAs Synergistically Promote Tolerance to Epidermal Growth Factor Receptor-Targeted Drugs in Smoked Lung Cancer Therapies. J. Can. Res. Ther. 15, 876–881. doi:10.4103/jcrt.jcrt_208_18
Piñero, J., Bravo, À., Queralt-Rosinach, N., Gutiérrez-Sacristán, A., Deu-Pons, J., Centeno, E., et al. (2017). DisGeNET: a Comprehensive Platform Integrating Information on Human Disease-Associated Genes and Variants. Nucleic Acids Res. 45, D833–D839. doi:10.1093/nar/gkw943
Piwecka, M., Glažar, P., Hernandez-Miranda, L. R., Memczak, S., Wolf, S. A., Rybak-Wolf, A., et al. (2017). Loss of a Mammalian Circular RNA Locus Causes miRNA Deregulation and Affects Brain Function. Science 357, eaam8526. doi:10.1126/science.aam8526
Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res. 43, e47. doi:10.1093/nar/gkv007
Sakai, H., Dobashi, K., Iizuka, K., and Nakazawa, T. (1995). Protective Effect of an α2-Adrenoceptor Antagonist, Midaglizole, Against Allergen-Provoked Late Asthmatic Responses. J. Asthma 32, 221–226. doi:10.3109/02770909509089511
Salzman, J., Chen, R. E., Olsen, M. N., Wang, P. L., and Brown, P. O. (2013). Cell-type Specific Features of Circular RNA Expression. Plos Genet. 9, e1003777. doi:10.1371/journal.pgen.1003777
Schaefer, N., Li, X., Seibold, M. A., Jarjour, N. N., Denlinger, L. C., Castro, M.., et al. (2019). The Effect of BPIFA1/SPLUNC1 Genetic Variation on its Expression and Function in Asthmatic Airway Epithelium. JCI Insight 4. doi:10.1172/jci.insight.127237
Shimomura, H., Okada, R., Tanaka, T., Hozaka, Y., Wada, M., Moriya, S., et al. (2020). Role of miR-30a-3p Regulation of Oncogenic Targets in Pancreatic Ductal Adenocarcinoma Pathogenesis. Int. J. Mol. Sci. 21, 21. doi:10.3390/ijms21186459
Silva, C. A., Reber, L., and Frossard, N. (2006). Stem Cell Factor Expression, Mast Cells and Inflammation in Asthma. Fundam. Clin. Pharmacol. 20, 21–39. doi:10.1111/j.1472-8206.2005.00390.x
Su, G., Morris, J. H., Demchak, B., and Bader, G. D. (2014). Biological Network Exploration with Cytoscape 3. Curr. Protoc. Bioinformatics. 47, 8–24. doi:10.1002/0471250953.bi0813s47
Szklarczyk, D., Gable, A. L., Nastou, K. C., Lyon, D., Kirsch, R., Pyysalo, S., et al. (2021). The STRING Database in 2021: Customizable Protein-Protein Networks, and Functional Characterization of User-Uploaded Gene/measurement Sets. Nucleic Acids Res. 49, D605–D612. doi:10.1093/nar/gkaa1074
Thaikoottathil, J. V., Martin, R. J., Di, P. Y., Minor, M., Case, S., Zhang, B., et al. (2012). SPLUNC1 Deficiency Enhances Airway Eosinophilic Inflammation in Mice. Am. J. Respir. Cel Mol Biol. 47, 253–260. doi:10.1165/rcmb.2012-0064oc
To, T., Stanojevic, S., Moores, G., Gershon, A. S., Bateman, E. D., Cruz, A. A., et al. (2012). Global Asthma Prevalence in Adults: Findings From the Cross-Sectional World Health Survey. BMC Public Health 12, 204. doi:10.1186/1471-2458-12-204
Wang, H., Kanmangne, D., Li, R., Qian, Z., Xia, X., Wang, X., et al. (2020). miR-30a-3p Suppresses the Proliferation and Migration of Lung Adenocarcinoma Cells by Downregulating CNPY2. Oncol. Rep. 43, 646–654. doi:10.3892/or.2019.7423
Wang, H.-y., Dai, Y., Wang, J.-l., Yang, X.-y., and Jiang, X.-g. (2015). Anti-CD69 Monoclonal Antibody Treatment Inhibits Airway Inflammation in a Mouse Model of Asthma. J. Zhejiang Univ. Sci. B. 16, 622–631. doi:10.1631/jzus.b1400285
Wang, K., Chen, M., and Wu, W. (2017). Analysis of microRNA (miRNA) Expression Profiles Reveals 11 Key Biomarkers Associated With Non-Small Cell Lung Cancer. World J. Surg. Onc. 15, 175. doi:10.1186/s12957-017-1244-y
Wang, P. L., Bao, Y., Yee, M.-C., Barrett, S. P., Hogan, G. J., Olsen, M. N., et al. (2014). Circular RNA is Expressed Across the Eukaryotic Tree of Life. PLoS One 9, e90859. doi:10.1371/journal.pone.0090859
Wang, S., Zhang, Y., Cai, Q., Ma, M., Jin, L. Y., Weng, M., et al. (2019). Circular RNA FOXP1 Promotes Tumor Progression and Warburg Effect in Gallbladder Cancer by Regulating PKLR Expression. Mol. Cancer 18, 145. doi:10.1186/s12943-019-1078-z
Wang, Y., Wang, F., He, J., Du, J., Zhang, H., Shi, H., et al. (2019). miR-30a-3p Targets MAD2L1 and Regulates Proliferation of Gastric Cancer Cells. Onco. Targets Ther. Vol. 12, 11313–11324. doi:10.2147/ott.s222854
Warde-Farley, D., Donaldson, S. L., Comes, O., Zuberi, K., Badrawi, R., Chao, P., et al. (2010). The GeneMANIA Prediction Server: Biological Network Integration for Gene Prioritization and Predicting Gene Function. Nucleic Acids Res. 38, W214–W220. doi:10.1093/nar/gkq537
Weigl, M., Kocijan, R., Ferguson, J., Leinfellner, G., Heimel, P., Feichtinger, X., et al. (2021). Longitudinal Changes of Circulating miRNAs During Bisphosphonate and Teriparatide Treatment in an Animal Model of Postmenopausal Osteoporosis. J. Bone Miner Res. 36, 1131-1144. doi:10.1002/jbmr.4276
Wong, C. H., Lou, U. K., Li, Y., Chan, S. L., Tong, J. H., To, K.-F., et al. (2020). CircFOXK2 Promotes Growth and Metastasis of Pancreatic Ductal Adenocarcinoma by Complexing With RNA-Binding Proteins and Sponging MiR-942. Cancer Res. 80, 2138–2149. doi:10.1158/0008-5472.can-19-3268
Woodruff, P. G., Boushey, H. A., Dolganov, G. M., Barker, C. S., Yang, Y. H., Donnelly, S., et al. (2007). Genome-wide Profiling Identifies Epithelial Cell Genes Associated With Asthma and With Treatment Response to Corticosteroids. Proc. Natl. Acad. Sci. 104, 15858–15863. doi:10.1073/pnas.0707413104
Woodruff, P. G., Modrek, B., Choy, D. F., Jia, G., Abbas, A. R., Ellwanger, A., et al. (2009). T-helper Type 2-driven Inflammation Defines Major Subphenotypes of Asthma. Am. J. Respir. Crit. Care Med. 180, 388–395. doi:10.1164/rccm.200903-0392oc
Yekta, S., Shih, I. H., and Bartel, D. P. (2004). MicroRNA-Directed Cleavage of HOXB8 mRNA. Science 304, 594–596. doi:10.1126/science.1097434
Yoshie, Y., Iizuka, K., and Nakazawa, T. (1988). Inhibitory Effect of DG-5128 on Steroid-Dependent Asthma. J. Asthma 25, 321–324. doi:10.3109/02770908809071373
Yu, C.-Y., Li, T.-C., Wu, Y.-Y., Yeh, C.-H., Chiang, W., Chuang, C.-Y., et al. (2017). The Circular RNA circBIRC6 Participates in the Molecular Circuitry Controlling Human Pluripotency. Nat. Commun. 8, 1149. doi:10.1038/s41467-017-01216-w
Keywords: asthma, airway epithelial cells, microarray datasets, mRNA, microRNA, circular RNA, robust rank aggregation, integrated analysis
Citation: Chen D, Wu W, Yi L, Feng Y, Chang C, Chen S, Gao J, Chen G and Zhen G (2021) A Potential circRNA-miRNA-mRNA Regulatory Network in Asthmatic Airway Epithelial Cells Identified by Integrated Analysis of Microarray Datasets. Front. Mol. Biosci. 8:703307. doi: 10.3389/fmolb.2021.703307
Received: 30 April 2021; Accepted: 30 June 2021;
Published: 16 July 2021.
Edited by:
Michele Trabucchi, Institut National de la Santé et de la Recherche Médicale (INSERM), FranceReviewed by:
Luca Parca, Casa Sollievo della Sofferenza (IRCCS), ItalyQian Xu, Central South University, China
Copyright © 2021 Chen, Wu, Yi, Feng, Chang, Chen, Gao, Chen and Zhen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Guohua Zhen, Z2h6aGVuQHRqaC50am11LmVkdS5jbg==