Gene network analysis for identification of microRNA biomarkers for asthma

Abstract

Background

To date, reliable biomarkers for asthma have not been identified. MicroRNAs (miRNAs) are small, non-coding RNAs that negatively regulate post-transcriptional gene expression, and they are involved in various diseases, including asthma. MiRNAs may serve as ideal biomarkers due to their ability to regulate multiple pathways. This study aims to identify miRNA biomarker signatures for asthma.

Methods

We used the house dust mite (HDM) mouse model of allergic inflammation. Mice were phenotyped by assessing lung function, allergic response, airway inflammation, and remodeling. The miRNA signature profiles in serum and lung tissue were determined by small RNA sequencing, and data were analyzed using Qiagen CLC Genomics Workbench. To identify relevant gene targets, we performed mRNA sequencing, followed by miRNA-targets analysis. These miRNAs and targets were subject to subsequent pathway and functional analyses.

Results

Mice exposed to HDM developed phenotypic features of allergic asthma. miRNA sequencing analysis showed that 213 miRNAs were substantially dysregulated (FDR p-value < 0.05 and fold change expression >  + 1.5 and < − 1.5) in the lung of HDM mice relative to the control mice. In contrast, only one miRNA (miR-146b-5p) was significantly increased in serum. Target analysis of lung dysregulated miRNAs revealed a total of 131 miRNAs targeting 211 mRNAs. Pathway analysis showed T helper 2/1 (Th2/Th1) as the top significantly activated signaling pathway associated with the dysregulated miRNAs. The top enriched diseases were inflammatory response and disease, which included asthma. Asthma network analysis indicated that 113 of 131 miRNAs were directly associated with asthma pathogenesis.

Conclusions

These findings suggest that most dysregulated miRNAs in the HDM model were associated with asthma pathogenesis via Th2 signaling. We identified a panel of 30 miRNAs as potential biomarker candidates for asthma.

Introduction

Asthma is a prevalent chronic inflammatory disease of the airways that affects over 300 million people worldwide and 25 million people in the United States [1]. It is characterized by airway hyper-responsiveness (AHR), remodeling, and smooth muscle cell hyperplasia and hypertrophy [2]. Asthma is classified into multiple phenotypes and endotypes, including allergic, eosinophilic, and neutrophilic with allergic asthma the most prevalent subtype [3]. Heterogeneity in the origin of inflammation contributes to challenges in finding effective treatment and there is no cure or specific biomarkers for asthma. Thus, identifying reliable biomarkers to predict disease pathogenesis would fill a gap in our knowledge that could be useful in therapeutic development.

There is accumulating data supporting a role for microRNAs (miRNAs) as biomarkers for asthma. miRNAs are small, non-coding RNAs that negatively regulate post-transcriptional gene expression via mRNA destabilization or/and degradation [4, 5]. They are involved in various cellular processes such as cell proliferation, differentiation, migration, and apoptosis [6] and are dysregulated in various pathologies, including asthma [7,8,9,10,11,12]. They are additionally implicated in all key pathogenetic processes of asthma, including airway inflammation, airway hyper-responsiveness (AHR) and hypercontractility. For instance, in a mouse model of allergic asthma, down-regulation of miRNA-133a was associated with increased expression of interleukin 13, and AHR, indicating a potential role of this miRNA in inflammation and airway contractility [13]. Moreover, down-regulation of miRNA-145 alleviated eosinophilic inflammation, mucus hyper-secretion, Th2 production, and AHR [14] in a house dust mite (HDM) model of allergic inflammation. Recent studies have focused on the potential use of miRNAs as biomarkers in various inflammatory diseases, including asthma. For example, in asthma patients, circulating miRNAs were found to be associated with lung function (miR-126-3p, miR-1290, miR-146b-5p, and miR-206), asthma severity (miR-185-5p), and exacerbations (miR-146b-5p, miR-206, and miR-720) [10, 11, 15].

While the use of circulating miRNAs as potential biomarkers for asthma has been investigated, conflicting results that associate circulating miRNAs and asthma have impeded progress toward finding reliable diagnostic and prognostic biosignatures. For instance, Rodrigo-Muñoz et al., identified a panel of three circulating miRNAs (miR-185-5p, miR-320a, and miR-320b) as biomarkers for asthma [11]. In contrast, another group showed a distinct panel of four miRNAs (miRNA-16-5p, miR-223-3p, miR-570-3p, and miR-299-5p) differentially expressed in asthmatic patients relative to the healthy group and individuals with allergic rhinitis [12]. Subsequent study found a subset of seven miRNAs comprising five miRNA ratios as biomarker signatures for asthma [16].

Another confounding problem is that the expression profiles of miRNAs in response to physiological or pathological stimuli are tissue-specific. Therefore, lung tissue might be a more promising source of biomarker discovery for asthma than circulating miRNAs. In this study, we compared lung miRNA signatures to circulating miRNA profiles using a HDM preclinical model of asthma to address the hypothesis that the lung contains distinct miRNA signature patterns that can serve as unique biomarkers for asthma. We have identified a panel of 30 miRNAs as biomarker candidates for asthma. While 18 of the 30 miRNAs had been previously associated with asthma pathogenesis, we have found 12 novel miRNAs (miR-3473b, 7061-5p, 217-5p, 369-3p, 411-3p, 381-3p, miR-7224-3p, 491-5p, 3097-5p, miR-6540-5p, 33-3p, and 1943-5p) correlated with asthma features.

Materials and methods

Materials

Sterile endotoxin-free 0.9% saline solution (cat# IAX-900-003) was obtained from Innaxon, UK. House dust mite (HDM, Dermatophogoides pteronyssinus, cat # XPB70D3A2.5) extracts were purchased from Stallergenes (Greer Laboratories, Lenoir, NC). Methacholine (acetyl-β-methylcholine chloride) was from Sigma Life Sciences (St. Louis, MO). RNA extraction and cDNA library preparation kits were from QIAGEN (Germantown, MD).

Methods

Mice

Male and female BALB/c mice, aged 6–8 weeks old, were purchased from the Jackson Laboratory and housed in a pathogen-free environment in the laboratory animal medicine facility located in the Center for Molecular Medicine at the University of Nevada, Reno (UNR), School of Medicine. All experimental procedures were approved by the institutional Animal Care and Use Committee (IACUC).

HDM sensitization and challenge protocol

The HDM protocol was carried out as previously described [17]. Briefly, mice were sensitized with an intranasal administration of 50 μl of a 25 μg HDM solution or saline once daily for five consecutive days. This was followed by a daily challenge of HDM or saline for five consecutive days per week for 3 weeks.

Tracheostomy and lung function assessment

Mice were anesthetized by an intraperitoneal (i.p.) injection of a ketamine (90 mg/kg)/xylazine (10 mg/kg) cocktail. The anesthesia depth was monitored with a toe pinch and heart rate assessment. Mice were laid supine and hair removed from the chest and neck area using commercial hair removal product followed by wet cotton swabs to reduce irritation. The skin, muscles, and fat in the neck were removed, with stainless-steel surgical scissors, to expose the trachea. A small lateral incision was made to the upper trachea, and a surgical tube was inserted 1/3 deep into the trachea above the primary bronchi. The surgical tube was secured into the trachea using a sterile suture. Additional ketamine was often administered after the tracheostomy to ensure proper sedation before lung function assessment. Mice were mechanically ventilated and lung airway resistance and compliance were measured in response to increasing doses (0–12.5 mg/ml in PBS) of aerosolized methacholine using Buxco equipment and FinePoint software from Data Sciences International (St. Paul, MN).

Blood collection and serum preparation

Blood was collected after lung function analysis via cardiac puncture as previously described [18] using a 23G needle attached to a 1 ml syringe. After collection, the needle was removed from the syringe, and blood was transferred into serum gel with clotting activator tubes (Sarstedt, Newton, NC). The blood was clotted at room temperature (RT) for 30 min, followed by centrifugation at 2000×g for 10 min at 4 °C, and the serum was transferred to a new 1.5 ml tube. The serum was then centrifuged at 16,000×g at 4 °C for 10 min to remove additional debris and minimize gDNA contamination. The clarified serum was aliquoted into 200 μl aliquots.

Bronchoalveolar lavage and Luminex assays

After the cardiac puncture procedure, bronchoalveolar lavage (BAL) was performed on tracheostomized mice. A total of five lavages were performed per mouse using a 1 ml syringe with a 0.8 ml of sterile 2.6 mM EDTA-saline solution (BAL solution) per lavage. Cell pellets from all five lavages were collected by centrifugation (300×g at 4 °C for 5 min), combined into a single tube, resuspended in 500 μ of BAL solution, and subject to cell counting using a Coulter counter (Beckman Coulter Z series, Brea, CA). The supernatants from the first two lavages were pooled and used for Luminex assays. The Luminex multiplex assays were performed at Eve Technologies Corp (Calgary, AB, Canada) using 50 µl BAL fluid per sample with the Mouse Cytokine Array/Chemokine Array 32-Plex Discovery Assay® (MD31) according to the manufacturer's instructions.

Lung tissue collection and histology

The right lung lobes and post-caval lobe were collected and placed, at one lobe per tube, into 1 ml tubes containing RNAlater solution (Invitrogen, ThermoFisher Scientific, cat# AM7021, Carlsbad, CA) and stored at − 80 °C for downstream analysis. The left lung lobe was inflated and fixed with 10% neutral buffered formalin (Fisher Scientific, Kalamazoo, MI) at 25 cm above the mouse level. Fixed lungs were embedded in paraffin (FFPE), cut into five μM sections onto positively charged slides. Slides were stained with hematoxylin–eosin (H&E) and Masson’s trichrome according to the Reveal Biosciences’ standard protocol. The whole slide images were obtained using a Panoramic SCAN (3D Histech, Hungary). The quantitative analysis of whole slide images was performed using imageDX software (Reveal Biosciences, San Diego, CA). The inflammatory cells in the H&E-stained lung were determined and quantitated as the number of immune cells within the total image analysis area (mm²). The collagen and extracellular matrix (ECM) depositions in the Masson’s Trichrome-stained sections were evaluated as a percent of the entire image analysis area (mm²). The representative 20X photos were acquired from whole slide images using the CaseViewer software (3D Histech, Hungary).

Detection of serum HDM-specific IgE

HDM-specific IgE in the mouse serum was quantified using the mouse serum anti-HDM IgE antibody assay kit (cat# 3037, Chondrex, Redmond, WA) according to the manufacturer’s recommendations. Briefly, anti-HDM Mouse IgE antibody standards (diluted at a range of 0–50 ng/ml) and mouse serum samples (diluted at 1:10) were added in duplicate wells onto a 96-well ELISA plate precoated with anti-mouse IgE antibody. The plate was incubated for 2 h at RT and washed 3× with wash buffer. The plate was incubated with the biotinylated HDM detection antibody for 2 h at RT and washed 3×. The streptavidin peroxidase was added, incubated at RT for 30 min, and washed 3×. The plate was incubated with TMB (3,3′,5,5′-tetramethylbenzidine) substrate for 25 min at RT, the reaction was stopped with 2N sulfuric acid, and the plate was read at 450 mm.

microRNAs sequencing and bioinformatic analysis

microRNA sequencing (miRNA- seq) from serum and lung tissue was performed at the Qiagen Genomic Center (Hilden, Germany) according to the facility's established protocols. For serum samples, RNA was extracted from 200 μl of serum using the miRNeasy serum/plasma advanced kit, and the RNA quality control (QC) was determined by QPCR of serum endogenous miRNAs (such as miRNA-103a-3p and miR-191-5p), hemolysis markers (miR-23a and miR451a), and Spike-Ins control (UniSp6). For lung tissue samples, RNA was extracted with the RNeasy mini kit, quantified using a Nanodrop, and the RNA integrity was assessed using the Agilent TapeStation. miRNA libraries for serum and lung were prepared using Qiagen's miRNA Library Kit following the manufacturer's protocol. cDNA libraries with Unique Molecular Index (UMI) assignment were amplified, cleaned up, QC done by Bioanalyzer, and sequenced using single-end reads, 1 × 75 bp, to a depth of 12–14 million reads (NextSeq 550, Illumina). FASTQ files with quality scores (QS) > 30 were analyzed by CLC Genomics Workbench (version 12.0.2) using the QIAseq miRNA quantification workflow. The de-duplicated Unique Molecular Index (UMI) reads were grouped into unique search sequences and mapped to the miRNAs annotated in the miRbase v22. miRNA expression was calculated via counting the number of distinctive UMI reads corresponding to each miRbase-annotated miRNA and normalized to the sequencing depth (counts per million). Differential miRNA expression between control and HDM-treated mice was assessed using the “Exact Test” for two-group comparisons available as part of the EdgeR Bioconductor package in the CLC Genomic Workbench. The microRNA-Seq data have been deposited at the NCBI Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo) under series accession number GSE214937.

mRNA sequencing and bioinformatic analysis

Libraries were prepared from the same lung RNA samples used for miRNA sequencing using QIAseq Stranded mRNA Select Kit with poly-A enrichment. cDNA libraries were amplified, cleaned up, QC performed with a Bioanalyzer, and sequenced using paired-end 2 × 75 bp read, to a depth of 60 million reads (NextSeq 550, Illumina) FASTQ files with quality scores (QS) > 30 were analyzed, mapped, and quantified using the RNA-seq analysis tool from CLC Genomics Workbench (version 21.0.4). Reads were trimmed, mapped to the mouse genome version GRCm38, and annotated using the ENSEMBL Mus_musculus. GRCm38 version 98. Gene expression was calculated by counting number of reads mapping to the annotated gene loci, and the resulting expression values were normalized to transcripts per million (TPM). Differential gene expression between control and HDM-exposed mice was assessed using the “Exact Test” for two-group comparisons available as part of the EdgeR Bioconductor package in the CLC Genomic Workbench. The RNA-Seq data have been deposited at the NCBI Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo) under series accession number GSE214937.

miRNA targets identification and pathway analysis

miRNA target identification and pathway analysis were performed using Ingenuity Pathway Analysis (IPA, Qiagen). To determine the most significantly dysregulated molecules, miRNAs and mRNA datasets were filtered using the filter dataset tab on IPA. The cut-off values were log2 fold-change expression (log2 FC) > 0.58 or < − 0.58 and FDR p < 0.05. Next, the mRNA targets were identified using the miRNA target filter. It is known that one miRNA can target up to 200 mRNAs. Therefore, to reduce the list of mRNA targets, we applied filters to include diseases relevant to asthma pathogenesis, such as inflammatory disease, inflammatory response, and respiratory disease. Further, the pathways list was filtered to include those related to a disease-specific pathway, such as airway inflammation in asthma and cellular growth, proliferation, and development pathways (including Th1, Th2, and Th17 activation pathways). The resulting miRNA and mRNA targets were subject to the IPA core analysis [19].

cDNA synthesis and QPCR

The lung RNA samples from the same mice used for all the experiments listed above, including pulmonary function, lung histology, and miRNA/RNA sequencing, were utilized to validate the miRNA sequencing data. miRNA cDNA was synthesized from five ng of total RNA using TaqMan Advanced miRNA cDNA Kit (Cat# A28007, Thermofisher Scientific, Waltham, MA) according to the manufacturer’s recommendations. QPCR was performed on the QuantStudio3 Fast Real-time PCR system, 96-well, 0.1 ml plate with 2X Advanced TaqMan Fast Master Mix and specific miRNA probe-primer mix. miRNA expression was normalized to mmu-miR5121. The miRNA mmu-miR-5121 was chosen as reference miRNA because its expression in the dataset remains unchanged in response to HDM exposure.

Statistics

Mice characterization data were graphed and analyzed using Graph Pad Prism software. Data represent mean ± standard error mean (SEM). An unpaired t-test with Welch's correction was used for two-group comparison and p ≤ 0.05 is considered statistically significant. For the miRNA and RNA sequencing data, the false discovery rate (FDR) p-values were calculated using the Benjamini–Hochberg method to account for multiple testing issues, and FDR < 0.05 was considered significant. The networks, pathways, and functional analyses were done through IPA. Right-tailed Fisher’s Exact Test was used to calculate p-values, and p < 0.05 (− log₁₀ p-value > 1.3) was considered significant.

Results

HDM exposure increases airway resistance and induces allergic inflammatory response

To confirm the efficacy of HDM exposure on lung function, mice were anesthetized with ketamine/xylazine, tracheostomized, and airway resistance was measured in response to increasing doses of methacholine (0–12.5 mg/ml). HDM sensitized and challenged mice displayed a significantly increased lung resistance compared to the controls (Fig. 1A). These mice also showed a substantial elevation in the concentration of serum HDM-specific IgE relative to the saline controls (Fig. 1B). To analyze the effect of HDM sensitization and challenge on inflammation and remodeling, we assessed the total inflammatory cell count in the BAL fluid and performed histology in FFPE lung sections. The total BAL inflammatory cell counts were significantly higher in mice exposed to HDM compared to controls (Fig. 2A). Hematoxylin–eosin (H&E) staining revealed that HDM mice had a substantial increase in lung inflammation and thicker airways than the saline controls (Fig. 2B–D). However, the collagen deposition levels were unchanged between the two groups (Fig. 2E, F). These results suggest that HDM sensitization and challenge trigger airway inflammation and remodeling but do not affect ECM deposition.

Fig. 1

Assessment of airway hyper-responsiveness and IgE production. A Airway resistance was measured in response to increasing doses of methacholine, and B serum HDM-specific IgE was measured using ELISA. The graphs represent mean ± SEM, n = 16–18, *p < 0.05 and ****p < 0.0001 HDM-challenged mice versus alum control, two-way ANOVA followed by Dunnett's multiple comparisons (A) or unpaired t-test with Welch's correction (B)

Fig. 2

Total inflammatory cell counts in Bronchioloalveolar lavage (BAL) and lung histology. A Total BAL cell counts were performed by coulter counting. C and E are representative images of lung sections stained with hematoxylin and eosin (C), and mason trichrome (E). B and D are quantification of inflammatory cell count (B) and epithelial airway thickness (D) in H&E-stained lungs. F The quantification of collagen deposition in mason trichrome-stained lung sections. Data represent mean ± SEM, n = 12–18 per group, ***p < 0.001 and ****p < 0.0001 HDM-challenged mice relative alum control using unpaired t-test with Welch's correction

The effect of HDM sensitization and challenge on various inflammatory cytokines and chemokines was assessed via Luminex assays. The secretion of major Th2 cytokines, IL-4, 5, 6, 10, and 13, were significantly enhanced in HDM mice compared to the controls (Fig. 3A–E). In contrast, IL-1α levels were reduced in HDM mice relative to controls (Fig. 3F). Additionally, the production levels of chemokines, including eotaxin/CCL11, IP-10/CXCL10, KC/CXCL3, MCP-1/CCL2, MIG/CXCL9 and MIP-1β/CCL4, and were significantly higher in BAL fluid of HDM mice relative to saline controls (Fig. 4). As expected, the HDM model exhibited increased production of key inflammatory mediators. Taken together, these results validate the use of the HDM model for these studies. We also analyzed the phenotypic features of the HDM mouse model based on sex and found that female HDM-exposed mice exhibited a similar allergic inflammatory profile as the male HDM counterparts (data not shown). Therefore, both sexes were combined for all subsequent analyses.

Fig. 3

Quantification of cytokine production in the BAL fluid. The cytokine secretion in BAL fluid was analyzed via Luminex assays. A IL-4, B IL-5, C IL-6, D IL-10, E IL-13 and F IL-1α secretion levels. Data are expressed as mean ± SEM, n = 16–18, **p < 0.01 ***p < 0.001 and ****p < 0.0001 versus alum control, unpaired t-test with Welch's correction

Fig. 4

Quantification of chemokine secretion in the BAL fluid. The inflammatory chemokine secretion in BAL fluid was analyzed via Luminex assays. A Eotaxin/CCL11, B IP-10/CXCL10, C KC/CXCL3, D MCP-1/CCL2, E MIG/CXCL9 and F MIP-1β/CCL4, A–F chemokine quantities. Data are expressed as mean ± SEM, n = 16–18, **p < 0.01 ***p < 0.001 and ****p < 0.0001 versus alum control, unpaired t-test with Welch's correction

Dysregulation of miRNA expression in the HDM mouse model of allergic inflammation

Principal component analysis (PCA) of miRNA expression in lung samples showed two nicely separated clusters between mice treated with saline and HDM (Fig. 5A). In contrast, there was no clear clustering observed between saline and HDM mice in serum samples (Fig. 5B). In the lung tissue, a total of 345 miRNAs were significantly differentially expressed (FDR < 0.05 or − log₁₀ (FDR) > 1.3) in HDM mice relative to control (Fig. 5C), but only one miRNA (miR-146b-5p) was substantially upregulated (FDR < 0.05 or − log₁₀ (FDR) > 1.3) in serum (Fig. 5D, E). Thus, the rest of this study focuses on lung miRNA signatures. Further analysis of the lung miRNA profile indicated that 125 miRNAs were up-regulated and 88 miRNAs were downregulated by at least 1.5-fold change in HDM mice compared to saline control (Fig. 5E, F). This lung miRNA signature was the focus of the rest of this study (Fig. 5F and Table 1).

Fig. 5

miRNA signature profiles in the house dust mite (HDM) mouse model of allergic inflammation. A Principal component analysis (PCI) of lung samples, B PCI of serum samples, red dots represent mice exposed to HDM, and black dots are control mice. C, D Volcano plots of differentially perturbed miRNAs in the lung tissue (C) and serum samples (D), red dots indicate the significantly dysregulated miRNAs (FDR < 0.05), dotted line show the significance threshold (− log₁₀ (FDR) = 1.3). A − log₁₀ (FDR) > 1.3 is considered as statistically significant. E Venn diagram of significantly up-regulated miRNAs in lung (green circle), serum (orange circle) samples, and down-regulated miRNAs in lung (blue circle). The miRNAs included in this analysis have a log2 fold change expression cut off >  + 0.58 and cut off < − 0.58 in HDM vs saline mice. n = 14–16 per group, and FDR < 0.05 is considered statistically significant. The Venn diagram was designed using the online Venny 2.1.0, https://bioinfogp.cnb.csic.es/tools/venny/index.html, with slight modifications. F Heat map of the 213 miRNAs that were substantially dysregulated in the HDM mice vs. saline controls. Red squares show upregulated miRNAs and green represent downregulated miRNAs. The heat map was generated by scaling the data using Z-scores of each miRNA (normalized CPM were converted into z-scores)

Table 1 Fold change and FDR values of dysregulated microRNAs

Validation of dysregulated lung miRNA signatures by QPCR

The top ten downregulated miRNAs were miR-1264-5p and 3p, 1298-3p and 5p, 1912-3p and 5p, 448-3p, 12181-3p, 764-3p and 208-5p. The top ten upregulated miRNAs were miR-8090, 135b-3p and 5p, 7681-5p, 147-3p and 5p, 7062-5p, 215-3p and 5p, and 6978-3p (Table 1). Selected miRNAs (1264-3p, 1298-5p, 7681-5p, and 511-3p) among the top 25 dysregulated miRNAs were validation by QPCR. We noted that miR-1264-3p and 1298-5p are significantly decreased in HDM mice relative to control mice (Fig. 6A, B), while miR-7681-5p and miR-511-3p expression levels are substantially enhanced (Fig. 6C, D). These results are consistent with the miRNA-seq findings.

Fig. 6

Analysis and QPCR validation of selected differentially dysregulated miRNAs. QPCR validation of downregulated (A, B) and upregulated (C, D) miRNAs. Data are expressed as mean ± SEM, n = 14–16, **p < 0.01, ***p < 0.001 and ****p < 0.0001 versus alum control, unpaired t-test with Welch's correction

Target identification of dysregulated lung miRNAs and pathway analysis of target genes

miRNA target analysis releveled a total of 131 microRNAs, with conserved seed sequence homology between mouse and human, targeting 211 mRNAs. There were 78 miRNAs upregulated and 53 miRNAs downregulated in HDM mice relative to saline controls. The expression pairing between miRNAs and corresponding mRNA targets is displayed in Additional file 2: Table S1. These molecules (Additional file 2: Table S1) were subject to pathway analysis to identify relevant signaling pathways and biological functions. We found that 38 canonical pathways were substantially impacted (− log₁₀ (p-value) > 1.3) with significant predicted activity (Z-score ≥  + 2) or inhibition (Z-score ≤ − 2) (Additional file 1: Fig. SA). The top 15 significantly activated pathways associated with the target genes were identified using IPA and included T helper 2 and 1 (Th2 and Th1), high mobility group box protein1 (HMGB1), interleukin-17 (IL-17), T and B cell signaling, and RhoA signaling (Fig. 7).

Fig. 7

Pathway analysis. The 15-top significantly activated canonical pathways correlated with target mRNAs, with the highest z-score (> 2.8) as indicated by numbers between parentheses. All pathways shown are activated (z-score ≥  + 2). Dotted line shows the significance threshold (− log₁₀ (p-value) = 1.3). Pathways with -log₁₀ (p-value) > 1.3 are marked as statistically significant. SP (signaling pathway)

Functional analysis of dysregulated miRNAs and their targets

Using pathway analysis, we identified several pathologies significantly associated with miRNAs and corresponding gene targets (Additional file 1: Fig. SB). The five top impacted pathologies were inflammation, organismal injury and abnormalities, connective tissue disorders, and skeletal and muscular disorders (Fig. 8A). The anomalies with significantly increased (z-score > 2) biological functions were inflammatory response, organismal injury and abnormalities, and inflammatory disease (Fig. 8B–D). These biological functions were mainly related to the activation, accumulation, and chemotaxis of immune cells (Fig. 8B), lung damage, injury, fibrosis (Fig. 8C), and airway hyper-responsiveness (Fig. 8D). All these functions are associated with phenotypic features of asthma. Additionally, we noted that asthma was substantially enhanced in the top-enriched disease categories (Fig. 8B–D), indicating the correlation of the lung miRNA signatures with asthma pathogenesis.

Fig. 8

Identification of diseases and biological functions associated with miRNA signatures in the HDM mouse model. A Represent the five top diseases associated with the HDM model. B–D Show biological functions associated with inflammatory response (B), organismal injury and abnormalities (C), and inflammatory diseases (D). Z-score corresponding to each biological function is shown in the parentheses. Z-score > 2 indicate increased function. The dotted line shows the significance threshold, − log₁₀ (p-value) > 1.3 are considered as statistically significant

Identification of miRNA and mRNA targets associated with an asthma network

We performed further downstream analysis, based on published data found in the IPA knowledge base, to determine genes and targeting miRNAs associated with asthma (z-score = 2.41) in mice exposed to HDM. This analysis identified 45 genes and three miRNAs (miR-135-b-5p, miR-148a-3p, and miR-215-5p). These genes included cytokines (IL-4, IL-5, and IL-13), chemokines (CCL2), growth factors (TGFβ1), peptidases (MMP9), transcription regulators (FOXP3 and GATA3), transmembrane receptors (IL13RA2 and IL17RB), enzymes and Kinases (Fig. 9A, Table 2). To connect the genes in Fig. 9A to the targeting miRNAs, we utilized the build and grow tabs in the path designer tool. These genes were targeted by 113 distinct miRNAs (Fig. 9B, Table 3). These results may indicate that most miRNAs that were differentially dysregulated in HDM mice are involved in asthma pathogenesis.

Fig. 9

Molecules associated with an asthma network in the HDM mouse model. A Genes correlated with asthma pathogenesis. Nodes represent the molecules and lines indicate the relationship between two nodes. The description of the node shapes is displayed in the legends to the left side. The node color intensities show the expression level degrees between HDM mice relative to control. The downregulated molecules are in green and upregulated in pink or red. Dotted lines indicate relationships between these molecules and asthma, as predicted by the IPA findings. Red lines indicate predicted activation, yellow lines represent inconsistent findings in the literature, and grey lines mean the effect was unknown. B miRNA-gene network interactions. miRNAs are colored in red (upregulated) or green (downregulated), and mRNA targets are represented in the lower big circle of image 7B. Grey arrows indicate an association between miRNAs and their targets

Table 2 Molecules directly associated with asthma

Table 3 miRNAs associated with asthma network

Identification of miRNA biomarkers for asthma

To investigate which miRNAs, serve as potential biomarkers, we first chose genes in Fig. 9A that were used in previous studies as biomarkers for asthma [20,21,22,23,24]. We found nine gene candidates (Fig. 10), including interferon gamma (IFNγ), Th2 cytokines (IL-4, IL-5, IL-6, IL-10, and IL-13), monocyte chemoattractant protein-1 (MCP-1/CCL2), transforming growth factor beta 1 (TGFβ1), and matrix metallopeptidase 9 (MMP9). These genes were grouped into four biomarker categories for asthma diagnosis, prognosis, and response of therapy (efficacy and safety) (Fig. 10). All nine genes were significantly dysregulated in our data set. The expressions of IFNγ, IL-4, IL-5, IL-6, IL-10, IL-13, MCP-1/CCL2, and TGFβ1 were increased and MMP9 was downregulated (Fig. 10, Additional file 3: Table S2). Next, we identified the targeting miRNAs for each of the genes mentioned above, and have discovered a panel of 39 miRNA candidates (Additional file 2: Table S2). The diagnostic values of candidate miRNAs were assessed using the receiver operating characteristic (ROC) analysis (Table 4). We found 17 miRNAs (miRNA-26a-5p, 24-3p, 148a-5p, 146b-5p, 147-5p, 503-5p, 135b-3p, 130b-5p, 369-3p, 195b, 3473b, 671-5p, 34c-3p, 340-5p, 204-5p, 381-3p and miR-125a-5p) with excellent predictive power (AUC: 0.92–1.00), 13 miRNAs (miRNA-9-3p, 491-5p, 33-3p, 7061-5p, 1a-3p, 6540-5p, 193a-5p, 3097-5p, 217-5p, 7224-3p, 296-5p, 411-3p and 1943-5p) with good predictive power (AUC: 0.80–0.88), and 9 miRNAs with low predictive values (AUC < 0.80) (Table 4). Overall, we have identified 30 miRNAs, with good to excellent predictive power, that may serve as potential biomarker candidates for asthma diagnosis.

Fig. 10

The miRNA biomarker candidates for asthma. Biomarker analysis was performed using prior knowledge from IPA. There are four miRNA biomarkers categories (efficacy, safety, diagnosis, and disease progression), and each category is connected to the molecules in the network by solid black lines. Nodes represent the molecules and lines indicate the relationship between two nodes. Green nodes represent downregulated miRNAs, and the upregulated ones are in pink or red. All mRNA targets, except MMP9, are represented in V-shaped dark pink. Each group of miRNAs is pointing to its specific mRNA target. BM: Biomarker

Table 4 ROC analysis of candidate miRNA for asthma diagnosis

Correlation of 30 miRNA biomarker candidates with asthma parameters

Next, we performed correlation analysis between the 30 miRNA biomarker candidates and various asthma parameters, including airway hyper-responsiveness (AHR), total serum IgE levels, and inflammatory cell counts. We note that 17 miRNAs showed positive correlation with asthma and 13 miRNAs displayed negative association with asthma (Table 5). Further, we note that 22 miRNAs (miRNA-135b-3p, 147-5p, 3473b, 148a-5p, 7061-5p, 146b-5p, 369-3p, 204-5p, 130b-5p, 671-5p, 340-5p, 411-3p, 503-5p, 381-3p, 125a-5p, 24-3p, 491-5p, 26a-5p, 34c-3p, 195b, 33-3p and miR-1a-3p) were significantly associated with all asthma features including IgE levels, AHR and inflammation. Four miRNAs (miRNA-9-3p, 296-5p, 193a-5p and 3097-5p) were significantly correlated with inflammation and IgE levels, but not with AHR. Two miRNAs (miR-6540-5p and 1943-5p) were inversely correlated with AHR but not associated with IgE levels or inflammation. The miR-7224-3p was negatively associated with inflammation and AHR but not linked with IgE levels, and the miRNA-217-5p was positively associated with inflammation only (Table 5). These data indicate that 22 of the 30 miRNAs tested were associated with three features of asthma, and eight of the 30 miRNAs were correlated with at least one asthma parameter.

Table 5 Correlation of 30 miRNA biomarker candidates with the features of asthma in HDM exposed mice

Discussion

The primary goal of this study was to identify miRNA biomarkers for asthma using a house dust mite (HDM) mouse model of allergic inflammation. Although several common allergens are used to induce allergic inflammation, HDM remains the most clinically-relevant antigen, as it affects nearly 85% of asthma patients worldwide [25]. Even though HDM is a well-established model, the outcomes of the allergic response and associated molecular signatures, such as miRNA profiles, can vary significantly depending on the dose of the allergen, sensitization route, mouse strain, and protocol used. Our results show that HDM- exposed mice exhibited an increased airway resistance, prominent airway inflammation and thickening, elevated serum HDM-specific IgE, and increased production of key inflammatory mediators. These results are in line with previous studies and indicate that the HDM protocol used in this study elicited a potent allergic inflammatory response with airway remodeling that recapitulates key features of human asthma [17, 26].

We utilized miRNA sequencing followed by bioinformatic analysis to identify miRNA signatures in lung tissue and serum of HDM-exposed mice relative to the saline control group. The comparison between lung tissue and serum miRNA profiles revealed 213 miRNAs significantly dysregulated in lung tissue, with the log2 fold change expression cut off (> + 0.58 and cut off < − 0.58) and FDR p-values < 0.05. Unexpectedly, only miR-146b-5p was substantially upregulated in the serum. These results suggest that lung tissue is a more reliable source of miRNA-based biomarker discovery for asthma than circulating miRNAs. To the best of our knowledge, this is the second study to analyze the free circulating miRNA profiles in a murine model of allergic inflammation. The first group (Milger K. et al.) to profile plasma in a mouse model of inflammation identified 11 significantly dysregulated miRNAs with the fold change (FC) from 0.59 to 1.75 and p < 0.01 [16]. While we only found one miRNA significantly upregulated (FC = 2.5 and FDR p-value = 5.13E−13). In Milger K. et al.’s study, miRNA signatures analysis was done using a focus panel containing 179 known miRNAs in plasma samples from female mice exposed to HDM or PBS. However, our study used serum from male and female mice, and miRNA profiles were determined via deep sequencing. The discrepancy between Miller et al.’s results and ours may be due to the gender and experimental methodology differences between the two studies.

The lung is a complex organ comprised of various cell types, including epithelial, smooth muscle, fibroblast, and endothelial cells [27] that undergoes structural changes and inflammatory cell recruitment when exposed to allergens such as HDM. Thus, analysis of miRNA signatures from whole lungs likely reflects a pool of various miRNAs from specific lung cell types and infiltrating immune cells, such as eosinophils. This would probably be more beneficial in unveiling the accurate landscape of the miRNA signature profile for asthma rather than using a specific lung cell type. Based on this rationale, we decided to focus solely on lung miRNA signatures since miR-146b-5p was the only miRNA dysregulated in serum, and miR146b-5p was also included in the whole lung analysis.

Target analysis of lung miRNA signatures showed 131 microRNAs targeting 211 mRNAs. These miRNAs had conserved seed sequence homology with human miRNAs, indicating their relevance to human asthma. Pathway analysis demonstrated that the mRNA targets were implicated in various immune and inflammatory signaling pathways, and Th2 signaling was the most significantly enriched pathway. Functional analysis indicated that the dysregulated miRNAs and corresponding mRNA targets were involved in several inflammatory responses and diseases, specifically in immune cell activation, accumulation, and chemotaxis, lung damage, lung injury, lung fibrosis, airway hyper-responsiveness, and asthma. Interestingly, all these biological functions are involved in the phenotypic features of asthma. Taken together, these findings suggest that lung miRNA signatures identified herein are associated with Th2-mediated asthma pathogenesis, the most predominant feature observed in mice models of allergic inflammation, including the BALB/c mouse strain [28]. Therefore, these miRNA signatures would be specifically relevant for the diagnosis of patients with the Th2-high asthma phenotype but not those with other asthmatic phenotypes, such as Th2-low asthma.

The asthma network analysis unveiled 113 miRNAs targeting 45 genes. These genes included cytokines (IL-4, IL-5, and IL-13), chemokines (CCL2), growth factors (TGFβ1), peptidases (MMP9), transcription regulators (FOXP3 and GATA3), transmembrane receptors (IL13RA2 and IL17RB), enzymes and kinases that were significantly dysregulated with 41 genes were upregulated and four were downregulated in our dataset. These results confirm the implication of these genes in asthma and agree with previous human asthma studies [3, 29].

The miRNA-gene network analysis identified 39 miRNAs out of 113 miRNAs as potential biomarkers for asthma. ROC analysis showed that 30 of the 39 candidates had a substantial predictive power (AUC ≥ 80 for each miRNA) for asthma diagnosis. The association of the 30 miRNA candidates with asthma was performed via correlation analysis between these miRNAs and three asthma parameters, including airway hyperresponsiveness (AHR), total serum IgE levels, and BAL inflammatory cell counts. The correlation analysis indicated that 17 of the 30 miRNAs were positively correlated with asthma, and 13 of the 30 miRNAs were inversely linked to asthma. Further, 22 of the 30 miRNAs were associated with IgE levels, AHR, and inflammation and eight were correlated with at least one asthma parameters. The association of most of these miRNAs with several asthma features indicate that they may serve as viable candidates for asthma diagnosis.

Out of the 30 identified miRNAs, 18 (miR-135b-3p, 147-5p, 148a-5p, 146b-5p, 204-5p, 130b-5p, 671-5p, 340-5p, 9-3p, 296-5p, 503-5p, 125a-5p, 193a-5p, 24-3p, 26a-5p, 34c-3p, 195b, and miR-1a-3p) had been previously linked to asthma pathogenesis [12, 30,31,32,33]. For instance, miRNA-296-5p, 125a-5p, 195b, and 26a-5p were correlated with lung function [15, 34, 35]. The miRNA-296-5p expression was associated with AHR in response to increasing methacholine PC₂₀ in the childhood asthma management program cohort [34]. The miR-26a-5p expression was associated with FEV₁/FVC ratio in asthma patients [15]. The miR-34c-3p was downregulated in asthmatic patients and its overexpression inhibited airway remodeling via the insulin growth factor binding protein-3 (IGFBP-3) signaling [36]. Other miRNAs such as miRNA-671-5p, miR-135b-3p, 147-5p, 340-5p and miR-1a-3p were dysregulated in murine models of asthma [8, 37]. These findings agree with our results and confirm the implication of several of these miRNAs in human asthma. Surprisingly, among the 30 identified potential biomarkers, the miR-146b-5p was the only upregulated miRNA in common between serum and lung samples. This finding corroborates several studies that also found an increased miR-146b-5p expression in various murine models of allergic inflammation, as well as in asthma patients [8, 10, 11]. Our findings confirm previous studies and indicate that miRNA-146b-5p is a viable biomarker for asthma.

This study also identified 12 novel miRNAs (miR-3473b, 7061-5p, 217-5p, 369-3p, 411-3p, 381-3p, miR-7224-3p, 491-5p, 3097-5p, miR-6540-5p, 33-3p, and 1943-5p) implicated in asthma pathogenesis. Interestingly, some of these miRNAs, specifically miR-3473b, 217-5p, and 369-3p, display a role in human inflammatory diseases other than the lung. For instance, the miR-3473b has been shown to participate in post-stroke neuroinflammation injury by targeting the suppressor of cytokine signaling 3 (SOCS3), a crucial negative regulator of Th2-mediated allergic responses [38]. Additionally, the miR-3473b targets the triggering receptor expressed on myeloid cells 2 (TREM2) and Unc-51 like autophagy activating kinase 1 (ULK1) to regulate the inflammatory response, such as the secretion of TNF-α and IL-1β, in Parkinson's disease [39]. Next, the downregulation of miR-217 was associated with decreased endothelial cell proliferation and enhanced secretion levels of IL-6, IL-1β, ICAM-1, and TNF-α. Whereas its overexpression was parallel with a reverse phenomenon [40]. Further, the upregulation of miR-369-3p in dendritic cells was associated with a decrease in pro-inflammatory (TNFα, IL-6, IL-12, IL-1α, IL-1β) and an increase in anti-inflammatory (IL-10 and IL-1RA) cytokine secretions in response to lipopolysaccharides (LPS)-induced inflammatory response, suggesting a possible anti-inflammatory role of this miRNA [41]. Overall, these findings indicate that miRNAs 3473b, 217-5p, and 369-3p participate in allergic inflammation by potentially regulating the production of various inflammatory cytokines.

In addition to the inflammatory role of some miRNAs, other candidates have been involved in fibrosis or proliferation. For instance, miR-411-3p has been shown to alleviate pulmonary and skin fibrosis through modulation of TGF-β/Smad ubiquitin regulatory factor 2 (Smurf 2) signaling, whereas miR-491-5p inhibits the non-small cell lung cancer (NSCLC) cell proliferation and migration by targeting forkhead box P4 (FOXP4), a transcription factor involved in the regulation of various cancers [42,43,44]. However, the implication of the remaining seven novel miRNAs (miR-7061-5p, 381-3p, miR-7224-3p, 3097-5p, miR-6540-5p, 33-3p, and 1943-5p) in lung biology still unknown.

Similar studies have been done in ovalbumin models using either male or female mice [8, 37, 45]. However, our study has several strengths, including using the HDM mouse model, a most clinically-relevant model of allergic inflammation. Additionally, we surveyed an extensive number of miRNAs using male and female mice. Further, each of the 30 miRNA candidates had good predictive power (AUC ≥ 80%), making these miRNAs stand-alone biomarkers for asthma diagnosis. The weakness of the study is that we did not confirm our findings in human samples. Despite this, several miRNAs discovered herein were also implicated in human asthma, indicating the validity of HDM models in profiling lung miRNAs. Subsequent studies will be performed to analyze how these miRNAs regulate asthma. Further, we will test if these miRNAs can serve as biomarkers in clinical settings using lung biopsy samples from asthmatic and healthy individuals.

Conclusions/future directions

These results suggest that lung tissue is a more reliable source of miRNA-based biomarker discovery for asthma than circulating miRNAs. Overall, these findings indicate that mouse models of inflammation might not be ideal models for circulating miRNA studies. Regardless of this plausible specific pitfall, the HDM model used herein was very useful in recapacitating all key features of human asthma, leading to the discovery of unique lung miRNA signatures for allergic asthma. Although most of these miRNAs were previously implicated with asthma, we have identified 12 novel biomarkers. Subsequent studies will validate the efficiency of these miRNA in asthmatic patients.