Systems Genomics Reveals microRNA Regulation of ICS Response in Childhood Asthma

Background: Asthmatic patients’ responses to inhaled corticosteroids (ICS) are variable and difficult to quantify. We have previously defined a Cross-sectional Asthma STEroid Response (CASTER) measure of ICS response. MicroRNAs (miRNAs) have shown strong effects on asthma and inflammatory processes. Objective: The purpose of this study was to identify key associations between circulating miRNAs and ICS response in childhood asthma. Methods: Small RNA sequencing in peripheral blood serum from 580 children with asthma on ICS treatment from The Genetics of Asthma in Costa Rica Study (GACRS) was used to identify miRNAs associated with ICS response using generalized linear models. Replication was conducted in children on ICS from the Childhood Asthma Management Program (CAMP) cohort. The association between replicated miRNAs and the transcriptome of lymphoblastoid cell lines in response to a glucocorticoid was assessed. Results: The association study on the GACRS cohort identified 36 miRNAs associated with ICS response at 10% false discovery rate (FDR), three of which (miR-28-5p, miR-339-3p, and miR-432-5p) were in the same direction of effect and significant in the CAMP replication cohort. In addition, in vitro steroid response lymphoblastoid gene expression analysis revealed 22 dexamethasone responsive genes were significantly associated with three replicated miRNAs. Furthermore, Weighted Gene Co-expression Network Analysis (WGCNA) revealed a significant association between miR-339-3p and two modules (black and magenta) of genes associated with immune response and inflammation pathways. Conclusion: This study highlighted significant association between circulating miRNAs miR-28-5p, miR-339-3p, and miR-432-5p and ICS response. miR-339-3p may be involved in immune dysregulation, which leads to a poor response to ICS treatment.


Introduction
Asthma is the most common chronic lung disease of childhood, with 44.3 percent of children in the United States reporting one or more asthma attacks in the previous year [1]. Inhaled corticosteroids (ICS) are the most effective and widely recommended controller medication for asthma. However, there are two key drawbacks: over 40% of patients do not respond well to ICS therapy, particularly those with severe disease, and ICS have dosedependent side effects [2][3][4]. Clinical asthma management would benefit from a method to identify patients who are poor responders to ICS before treatment begins, allowing other medications to be added or substituted for the ICS and avoiding the usual trial-and-error phase of treatment. Extracellular microRNAs (miRNAs) are miRNAs found outside of cells, in bodily fluids such as saliva, urine, breast milk, etc. They are released from cells through various mechanisms, such as passive leakage from damaged cells, active secretion in extracellular vesicles (EVs), or by binding to proteins such as Argonaute-2 (AGO2). The release of extracellular miRNAs is more localized and may affect cells that are in close proximity to the site of release. In contrast, circulating miRNAs are a type of extracellular miRNA that circulate in the bloodstream and have the potential to regulate gene expression in distant cells. They are mainly released by active secretion in extracellular vesicles or by binding to proteins such as high-density lipo-proteins (HDLs), apoptotic bodies, and microparticles [5]. Circulating microRNAs (miRNAs) have been suggested as biomarkers for a number of conditions and have been shown to be important in a number of inflammatorymediated processes [6,7], including asthma [8]. We hypothesized that circulating miRNA in serum would be associated with ICS response in children with asthma.
To distinguish between good and poor ICS responders, a precise quantitative definition of steroid response is necessary. Six clinical features were used by Clemmer et al. to establish a measure of the Steroid Responsiveness Endophenotype (SRE); while this method showed excellent performance, this measure is inherently longitudinal and necessitates surveillance during a period of ICS administration [9]. We have previously defined a Cross-sectional Asthma STEroid Response (CASTER) measure of ICS response, based on a combination of asthma control indicators and spirometry measures, which can be computed from data collected at a single time point [10]. Previously, analyses of gene or miRNA differential expression [11,12], metabolome, single-nucleotide polymorphism, and expression quantitative trait loci (eQTL) found multiple genes, miRNAs, metabolomes, or SNPs associated with asthma drug responsiveness [13][14][15][16].
The primary goal of this study was to identify circulatory miRNAs associated with ICS response in children with asthma and to interpret the mode of function of such miRNAs. We show that serum miRNAs were associated with ICS response in a cohort of Costa Rican children with asthma, and these findings were replicated in an independent cohort of North American children with asthma. Using transcriptomics of lymphoblastoid cell lines in response to steroids, we then constructed a steroid-responsive gene co-expression network, identified co-modulated genes clusters using a modified WGCNA method, and showed that these were associated with replicated miRNAs. The Genetics of Asthma in Costa Rica Study (GACRS) is a cross-sectional study that recruited from February 2001 to August 2008. It included 1165 asthmatic children between the ages of 6 and 14 years, each with a high probability of having six or more great-grandparents from the Central Valley of Costa Rica. A previous publication [17] described the GACRS protocols and assessments, which included FEV1, bronchodilator response (BDR), methacholine challenge, and questionnaires about prior medication use and exacerbations. GACRS did not specifically assess asthma severity. We used 2 major types of data from selected GACRS participants: clinical data (from which we calculated CASTER), and serum miRNA data.

Replication Cohort: The Childhood Asthma Management Program (CAMP)
CAMP is a multi-center, randomized, double-blinded, clinical trial of inhaled corticosteroids in 1041 children aged 5 to 12 years with mild-to-moderate persistent asthma who were followed for four years. In CAMP, asthma severity was defined at enrollment based on physician assessment (1 = no asthma, 2 = mild, and 3 = moderate asthma), and patients with mild and moderate asthma were enrolled in the cohort. A detailed definition is given in Saprio et al. 1999 [18]. CAMP participants predominantly self-identified as non-Hispanic white, but included small numbers of African American, Hispanic white, and other racial and ethnic groups (Table 1), which were grouped together in this analysis. The design [18] and outcomes [19] of the CAMP study have been previously reported ("other"). During the original trial, CAMP participants were randomized into 3 treatment arms: budesonide (an ICS), nedocromil, or placebo. We used 3 major types of data from selected CAMP participants: clinical data (from which we calculated CASTER), serum miRNA data, and transcriptomic profiling from lymphoblastoid cell lines.

CASTER: Cross-Sectional Asthma STEroid Response Measurement
Cross-sectional Asthma STEroid Response (CASTER) is a composite corticosteroid responsiveness phenotype or Steroid Responsiveness Endophenotype (SRE) measure. CASTER is calculated using five cross-sectional clinical measures: (1) oral steroid courses, (2) asthmarelated hospitalizations and ED (emergency department) visits, (3) pre-bronchodilator FEV1, (4) provocative concentration of methacholine required to effect a 20% reduction of FEV1 (PC20), and (5) bronchodilator response to albuterol (BDR), computed as a percentage change in FEV1 from baseline. In previous work, we showed that CASTER had good performance in childhood cohorts [10].

miRNA Sequencing and Quality Control
We performed small RNA sequencing on serum from 580 GACRS samples of asthma patients on ICS. Sequencing on 187 CAMP samples has been described previously [20]. Both cohorts were sequenced following the same protocols [20]. In brief, small RNA-seq libraries were prepared by using the Norgen Biotek Small RNA Library Prep Kit v2 (Norgen Biotek, Therold, ON, Canada) and sequenced on the Illumina NextSeq 500 platform. The ExceRpt pipeline was used for quality control (QC) of the RNA-seq data [21]. miRNAs with less than five mapped reads in at least 50% of subjects were removed. The data were normalized using DESeq2 [22].
Small RNA sequencing was completed on all available serum samples from GACRS [n = 1134]. Each sample produced an average of 15.8 million total reads. In total, 11.9 million reads per sample passed the initial quality control (QC). On average, 8.9 million reads per sample were mapped to the genome, which included 5.4 million miRNA sense sequences.
DESeq2 was used to perform normalization before association analysis. For GACRS cohort samples, sequencing was carried out in 2 batches, which may have introduced technical causes of discrepancy during preparation and handling, affecting the outcomes. Therefore, guided principal component analysis (gPCA) [23] was used to check for batch effects on mapped read counts per sample; the analysis showed that there was not a significant batch effect in normalized data (p-value = 0.41) (Supplemental Figure S1).

ICS-Response-Associated miRNA Identification
The CASTER measure of ICS response was computed for all subjects in GACRS and separately for CAMP. For more details, see [10], where partial least squares (PLS) distance and principal component (PCA) methods were used to classify good vs. poor ICS responders. We used mean (PCA-based) CASTER values from each cohort to dichotomize participants into poor and good responder groups. The patients having a CASTER value greater than the mean value of their cohort were considered to have a good response to ICS treatment, while those with a CASTER value less than the mean value were considered to have a poor response. Logistic regression with a Benjamini-Hochberg false discovery rate (FDR) correction for multiple testing was used to identify miRNAs associated with ICS response in GACRS using the GLM function in stats (v 4.1.2) R package. A significance threshold of 10% FDR was used. The analysis was adjusted for age and sex. Similarly, logistic regression was used to assess the association of serum miRNAs with ICS response using CASTER phenotype in CAMP while adjusted for age, sex, and race or ethnicity (non-Hispanic white, Hispanic, and others). Sensitivity analysis in CAMP was also performed, including asthma severity and genetic ancestry proportion as co-variates. A p-value < 0.05 was considered significant for replication in the CAMP cohort.
Logistic regression was also used to define a predictive model of ICS response. Three models were assessed in this study: Model1 included demographic variables (age, sex, race/ethnicity, height, weight, BMI), clinical variables (log10 IgE, vitamin, log10 Eosinophil count), smoking, and asthma severity; Model2 consisted of three miRNAs that were replicated; and Model3 combined variables from Model1 and Model2. These models were trained on the discovery cohort, GACRS, and were then evaluated on the replication cohort (CAMP) without parameter refitting. The area under the receiver operator characteristic curve (AUC) was used to measure predictive accuracy, and the AUC confidence interval was measured using 1000 bootstrap. The three models were compared using ANOVA.

CAMP Lymphoblastoid Cell Lines (LCL) Transcriptome Data
Some CAMP subjects also had gene array expression data available from a later timepoint. As previously described [14,15], subjects provided blood samples from which CD41 lymphocytes were isolated. These samples were immortalized and transcriptomically profiled under two treatment conditions: dexamethasone (DEX) (1026 mol/L) and a sham [ethanol] control. After 6 h, mRNA expression levels were measured with the Illumina Hu-manRef8 v2 BeadChip (Illumina, San Diego, CA, USA). Detailed protocols were previously described [14]. We performed quality control (QC) and filtering processes for the present study as described previously [15]. For this study, we looked at 88 samples of individuals on ICS for which serum miRNA data was also available.

Functional Annotation of Replicated miRNAs
The gene targets for three replicated CASTER-associated miRNAs were identified by multiMiR R package v 1.12 [24], where miRecords v 4, TarBase v 8, and miRTarBase v 7.0 databases were used to provide validated mRNA targets. These databases contain lists of validated miRNA target genes identified using experimental methods such as luciferase reporter assay, HITS-CLIP, CLASH, qRT-PCR, Western blot, degradome sequencing, immunocytochemistry, and others. The Database for Annotation, Visualization and Integrated Discovery (DAVID) v 2021 [25] was used to perform functional enrichment analysis of total unique targets of 3 replicated miRNAs, where gene ontology (molecular function and biological process) and pathway (KEGG and Reactome) datasets were considered. To show substantial enrichment of targeted genes for a pathway, we used a Bonferroni-adjusted p-value cut-off of 0.10 and a gene count of 3 or more.

Gene Expression Analysis
To identify differentially expressed genes in response to dexamethasone (DEX), gene expression from LCL cell lines was examined using the limma R package [26]. The association between three replicated miRNAs and differentially expressed genes in response to dexamethasone was then investigated using linear mixed modeling, using the lme4 R package [27]. We tested whether the direct targets of three miRNAs were over-represented among genes differentially expressed in response to DEX treatment using the hypergeometric test available in the stats R package. We identified comodulated genes (modules) using Weighted Gene Correlation Network Analysis (WGCNA) v 1.71 [28]. For the WCGNA input matrix, we used a modified version of Pearson correlation, following [29], to account for the correlation between the DEX and sham paired expression design. We then assessed the association of each module's eigengene (ME) to miRNA levels, where the module eigengene is defined as the first principal component of the module and represents the overall expression level of the module. To account for the within-pair correlation in data from the DEX-sham samples' paired design, we used a linear mixed-effects model (LMM) to test the association of a gene module to miRNA, with adjustment for age and sex using the lme4 R package [27]. We considered 10% FDR cutoff as a significance threshold.

GACRS
Small RNA-Seq data from the baseline serum were available for 1134 (97%) of the 1165 children in the GACRS. A total of 580 GACRS participants (51%) who were selfreported to use ICS in the prior 6 months and had enough data to calculate their CASTER phenotype were divided into two groups: those with CASTER values below the mean CASTER value (0.00403) were classified as poor responders (n = 379), and those with CASTER values above the mean were classified as good responders (n = 201) ( Table 1). Small RNA-Seq data on the baseline serum level were available for 492 (47%) of the 1041 children in the CAMP. A total of 187 CAMP participants (38%) who were on ICS and had enough data to calculate their CASTER phenotype were divided into two groups: those with CASTER values below the mean CASTER value (0.00547) were classified as poor responders (n = 116), and those with CASTER values above the mean CASTER value were classified as good responders (n = 71) ( Table 1).
In GACRS, older age (9.43 vs. 8.8) was associated with poor ICS response (p < 0.001), whereas in CAMP, there was a non-significant trend of younger age associated with poor ICS response. Total serum IgE level and eosinophil count were higher in the poor-responder group than in the good-responder group in both cohorts (Table 1). Further, the five constituent features of the CASTER measure differed from poor to good ICS responders, as expected, described presently. The poor-responder group patients had more courses of oral steroids (GACRS: 2.  (Table 1). Poor responders had a lower FEV1 than good responders in the GACRS (91.9 vs. 109), but not in CAMP. We noticed no significant differences by asthma severity at enrollment (in CAMP; unavailable in GACRS) or by genetic ancestry proportions.

Identification of miRNAs Associated with ICS Response
The sequenced serum miRNA datasets were subjected to-quality control (QC), filtration, and normalization, yielding 580 samples and 317 miRNAs for our discovery cohort (GACRS) and 187 samples and 257 miRNAs in our replication cohort (CAMP). Logistic regression showed that 36 of the 317 interrogated miRNAs were associated with ICS response at 10% FDR in the GACRS (Figure 1 and Supplemental Table S1). Of these miRNAs, 33 were associated with poor ICS response (odds ratio (OR) >1) and three miRNAs were associated with good ICS response (OR < 1). Of these 36 miRNAs, three (miR-28-5p, miR-339-3p and miR-432-5p) were validated as being significant in the CAMP cohort and in the same direction of effect (Table 2A,B, Supplemental Figure S5) at p < 0.05. All three replicated miRNAs were positively associated with poor ICS response (OR > 1, Table 2), which means that the increase in expression of these miRNAs leads to worse response to ICS treatment. In CAMP, the miRNA association finding remains consistent in sensitivity analysis adjusting by asthma severity status and genetic ancestry proportion.  In order to predict ICS response, three predictive models were assessed in this study. Each logistic regression model was trained on the GACRS cohort and tested on the CAMP cohort without parameter refitting. Model1, which incorporated demographic variables and clinical variables which were not constituents of the CASTER measure, had a 65% (CI: 57-73%) AUC in testing on the CAMP cohort (Supplemental Figure S2a). Model2, which included only the three replicated miRNAs, had a 64% (CI: 56-72%) AUC in the same cohort (Supplemental Figure S2b). Finally, Model3, which combined all variables from both Model1 and Model2, had a 72% (CI: 65-79%) AUC in the CAMP cohort (Supplemental Figure S2c). ANOVA was used to compare the models, and the results indicated a significant difference only between Model1 and Model3 (p = 0.0022).

Target Identification and Functional Enrichment Analysis of Replicated miRNAs
Next, we assessed the putative targets of these three miRNAs. We used the multiMiR R package to identify 1320 unique functionally validated gene targets for pathway and  In order to predict ICS response, three predictive models were assessed in this study. Each logistic regression model was trained on the GACRS cohort and tested on the CAMP cohort without parameter refitting. Model1, which incorporated demographic variables and clinical variables which were not constituents of the CASTER measure, had a 65% (CI: 57-73%) AUC in testing on the CAMP cohort (Supplemental Figure S2a). Model2, which included only the three replicated miRNAs, had a 64% (CI: 56-72%) AUC in the same cohort (Supplemental Figure S2b). Finally, Model3, which combined all variables from both Model1 and Model2, had a 72% (CI: 65-79%) AUC in the CAMP cohort (Supplemental Figure S2c). ANOVA was used to compare the models, and the results indicated a significant difference only between Model1 and Model3 (p = 0.0022).

Target Identification and Functional Enrichment Analysis of Replicated miRNAs
Next, we assessed the putative targets of these three miRNAs. We used the multi-MiR R package to identify 1320 unique functionally validated gene targets for pathway and ontology enrichment (Table 3, Figure 2A). We found enrichment of some common asthma pathogenesis-associated pathways such as PI3K-Akt, MAPK cascade, Wnt, Hippo, FoxO, and p53 signaling. We also discovered enriched ICS-response-related pathways such as estrogen receptor (ESR)-mediated, PI3K-Akt signaling (Supplemental Figure S3), and immune-response-and inflammation-associated pathways such as neutrophil degranulation and Interleukin-4 and -13 signaling ( Figure 2B).  Figure 2A). We found enrichment of some common asthm pathogenesis-associated pathways such as PI3K-Akt, MAPK cascade, Wnt, Hippo, FoxO and p53 signaling. We also discovered enriched ICS-response-related pathways such a estrogen receptor (ESR)-mediated, PI3K-Akt signaling (Supplemental Figure S3), and im mune-response-and inflammation-associated pathways such as neutrophil degranula tion and Interleukin-4 and -13 signaling ( Figure 2B).

In Vitro Steroid Response Lymphoblastoid Gene Expression
Paired mRNA gene expression data were available from immortalized lymphoblastoid cell lines treated with dexamethasone (DEX) or a sham vehicle treatment from some of the CAMP subjects. From these, we selected 88 paired DEX-sham arrays from individuals for which serum miRNA data was also available. A total of 3827 out of 4818 probes passing QC were differentially expressed in response to DEX (complete list by module in Supplemental Table S2). These genes showed significant enrichment of direct targets of the three replicated miRNAs (miR-28-5p, miR-339-3p, and miR-432-5p; p = 0.042). Further, association analysis between three miRNAs and DEGs showed 22 significant associations at 10% FDR (Supplemental Table S3 and Supplemental Figure S4).
WGCNA analysis was then conducted on the 88 paired DEX-sham expression arrays to identify associations between gene expression on ICS response and three miRNAs

In Vitro Steroid Response Lymphoblastoid Gene Expression
Paired mRNA gene expression data were available from immortalized lymphoblastoid cell lines treated with dexamethasone (DEX) or a sham vehicle treatment from some of the CAMP subjects. From these, we selected 88 paired DEX-sham arrays from individuals for which serum miRNA data was also available. A total of 3827 out of 4818 probes passing QC were differentially expressed in response to DEX (complete list by module in Supplemental Table S2). These genes showed significant enrichment of direct targets of the three replicated miRNAs (miR-28-5p, miR-339-3p, and miR-432-5p; p = 0.042). Further, association analysis between three miRNAs and DEGs showed 22 significant associations at 10% FDR (Supplemental Table S3 and Supplemental Figure S4).
WGCNA analysis was then conducted on the 88 paired DEX-sham expression arrays to identify associations between gene expression on ICS response and three miRNAs associated with ICS response. We identified 20 gene modules showing coordinated response of LCL cells to dexamethasone treatment. Next, we calculated the module eigengene (ME) for each of the 20 gene modules and used linear mixed-effects models to evaluate the association between each of the gene modules with the three miRNAs. We found that the black (OR = 0.99, FDR = 0.02) and magenta (OR = 0.99, FDR = 0.01) modules (Figure 3, Supplemental Table S4) were significantly associated with miR-339-3p. The functional enrichment analysis of these two modules showed that the magenta module is enriched with genes associated with adaptive immune response and inflammation pathways such as adaptive immune response, PIP3 activates AKT signaling, PTEN regulation, mTOR signaling, and glutathione metabolism. We also looked for an association between these 20 gene module eigengene values and clinical and demographic outcomes, but we did not find any significant association at the 5% FDR cut-off.

Discussion
MicroRNAs are emerging as pharmacogenomic biomarkers capable of identifying drug treatment response [30]. By enhancing the reliability of asthma treatment regimens, a biomarker for ICS response might subsequently help reduce asthma morbidity. In this study, we used CASTER, a composite corticosteroid responsiveness phenotype, to identify serum miRNAs as indicators of ICS response in two independent cohorts of children

Discussion
MicroRNAs are emerging as pharmacogenomic biomarkers capable of identifying drug treatment response [30]. By enhancing the reliability of asthma treatment regimens, a biomarker for ICS response might subsequently help reduce asthma morbidity. In this study, we used CASTER, a composite corticosteroid responsiveness phenotype, to identify serum miRNAs as indicators of ICS response in two independent cohorts of children with asthma. GACRS is a cross-sectional observational study with self-reported ICS use. CAMP, on the other hand, is a clinical trial in which all children in the ICS-arm received the same dose of the same ICS. Although GACRS is a cross-sectional study, and CASTER is designed to provide a cross-sectional measure of ICS response, this is still an inferred response phenotype that may differ from ICS response as revealed in a longitudinal cohort such as CAMP. We minimized such issues by also computing CASTER in the CAMP cohort, but some differences between the studies may have resulted in a lower number of replicated miRNAs.
We identified 36 miRNAs significantly associated with ICS response at 10% FDR cut-off; out of these, 17 miRNAs were found to have passed a more stringent 5% FDR cut-off, three of which (miR-28-5p, miR-339-3p and miR-432-5p) were found to be positively associated with poor ICS response in both cohorts (Supplemental Figure S5). Even at 5% FDR, both miR-28-5p and miR-339-3p remained significant, while miR-432-5p was borderline, at p = 0.058 (Table 2B). These three miRNAs showed some promise as potential pharmacogenomic biomarkers for ICS response with a 64% AUC in a test cohort when evaluated alone (Supplemental Figure S2b), providing a significant improvement to prediction based solely on clinical and demographic data (72% AUC). While we stress that this performance is with no parameter refitting in a completely independent cohort, and is therefore very encouraging, further work would be required to obtain more comprehensive models suitable for clinical use. miR-28-5p has been previously associated with inflammation in severe asthma [31]. miR-28-5p is a validated regulator of IL2RG (interleukin 2 receptor, gamma) (Supplemental Table S5), which is a pro-inflammatory cytokine receptor that plays an important role in glucocorticoid function in severe asthma [32,33]. miR-339-3p regulates the expression of NR3C1, which encodes an intracellular glucocorticoid receptor [34]. Fez et al. 2019 showed that miR-339-3p was up-regulated after 30 months of ICS treatment in COPD patients [35]. We are the first to report a link between miR-432-5p and asthma and ICS response, but this miRNA has previously been linked to cancer [36].
Functional enrichment analysis of their unique validated targets showed that these miRNAs regulate the expression of genes enriched in pathways of immune response and inflammation ( Figure 2B). This may indicate a possible role of these miRNAs in immune dysregulation. The previous literature has established a link between immune dysregulation and steroid insensitivity in patients [32,37,38]. We also discovered that these miRNAs regulate the expression of genes that were enriched in estrogen receptor (ESR)-mediated (Supplemental Figure S3) and PI3K-Akt signaling pathways (Supplemental Figure S3), which have been linked to steroid resistance in severe asthma [39]. According to Bi et al. 2020, inhibiting the PI3K-Akt signaling pathway with PI3K inhibitors mitigates glucocorticoid (GC) insensitivity in severe asthma by restoring HDAC2 activity and inhibiting the phosphorylation of nuclear signaling transcription factors [40]. It has previously been reported that estrogen inhibits glucocorticoid anti-inflammatory function. [41,42]. This suggests that these miRNAs indirectly contribute to glucocorticoid insensitivity by regulating the genes involved in the pathways. miR-339-3p is significantly associated with two DEX-responsive gene clusters ( Figure 3). The genes in these modules were found to be enriched in inflammation-related pathways: adaptive immune response, PIP3 activates AKT signaling, PTEN regulation, mTOR signaling, and glutathione metabolism (Figure 4). miR-339-3p showed a significant association (FDR = 0.02) with a DEX-responsive gene (PSMA7) (FDR = 0.0023) which participates in activation of NF-kappaB in B cells (R-HSA-1169091) (Supplemental Table S4). NF-kappaB is a pro-inflammatory transcription factor which is well-studied for its role in glucocorticoid response [43,44]. miR-339-3p also showed significant association with genes of adaptive immune system (CLTA, clathrin at FDR = 0.02 and RNF126, ring finger protein 126 at FDR = 0.04), mTOR-signaling pathway (EIF4EBP1, eukaryotic translation initiation factor 4E binding protein 1 at FDR = 0.04; CAB39, calcium binding protein 39 at FDR = 0.05), and HATs acetylate histones pathway (RUVBL2 at FDR = 0.06). The mTOR signaling pathway produces inflammatory or immune imbalance in asthma [45], and it has been reported that inhibition of mTOR activity can restore corticosteroid sensitivity in COPD [46]. The HATs acetylate histones pathway is critical in corticosteroid responsiveness. Corticosteroids suppress inflammatory genes in asthma by inhibiting HAT activity and recruiting HDAC2 to the activated inflammatory gene complex. HDAC2 activity and expression are reduced in severe asthma, which may account for the increased inflammation and resistance to corticosteroid action [47,48].These findings suggest that these three miRNAs, particularly miR-339-3p, play an important role in immune dysregulation, which results in poor response to ICS treatment. pathway produces inflammatory or immune imbalance in asthma [45], and it has been reported that inhibition of mTOR activity can restore corticosteroid sensitivity in COPD [46]. The HATs acetylate histones pathway is critical in corticosteroid responsiveness. Corticosteroids suppress inflammatory genes in asthma by inhibiting HAT activity and recruiting HDAC2 to the activated inflammatory gene complex. HDAC2 activity and expression are reduced in severe asthma, which may account for the increased inflammation and resistance to corticosteroid action [47,48].These findings suggest that these three miR-NAs, particularly miR-339-3p, play an important role in immune dysregulation, which results in poor response to ICS treatment.  The advantages of our study include a relatively large population for integrative pharmacogenomics in asthma, replication of the results in another childhood asthma cohort, use of a previously validated composite ICS response phenotype, and validation of the function of selected miRNAs using in vitro steroid response lymphoblastoid cell line gene expression data. Disadvantages of our study are the lack of mRNA expression data for the patient at the same time point as the LCL expression data, as well as racial and trial-design differences between our discovery and replication cohorts. While CASTER may be thought of as a measure of asthma severity, it would be one defined only in terms of ICS response: people doing poorly on ICS therapy are generally defined as having more severe asthma according to modern GINA (Global Initiative for Asthma) guidelines [49]. However, asthma severity was assessed at enrollment in CAMP in a treatment-independent manner, and this was not associated with the CASTER measure.

Conclusions
In conclusion, our findings show that serum miRNA expression differs among individuals according to response to ICS treatment. The miRNAs miR-28-5p, miR-339-3p, and miR-432-5p were associated with poor ICS response and may form part of a future biomarker of ICS response. Particularly, miR-339-3p may reduce ICS responsiveness through immune dysregulation.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/cells12111505/s1. Supplemental Methods; Figure S1: GACRS batch effect check. Figure S2: ICS response prediction using logistic regression in replication cohort (CAMP); (a) Model1, including age, sex, race/ethnicity, height, weight, BMI, log10 IgE, log10 Eosinophil, vitamin, smoking and severity; (b) Model2: three miRNAs (miR-28-5p, miR-339-3p, miR-432-5p); (c) Model3: variables of Model1 and Model2. Figure S3: Target genes of replicated miRs enriched in PI3K-AKT and estrogen signaling pathways; orange-colored boxes represent enriched target genes. Figure S4: Boxplot showing expression under DEX and sham condition for 22 differentially expressed genes associated with replicated miRNAs. Figure S5: BoxPlot showing expression in poor and good responder group for three replicated miRNAs in (A) GACRS and (B) CAMP. Table S1: miRNAs associated with ICS response in GACRS cohort. Table S2: Details of genes from LCL cell line used for differential expression analysis and WGCNA. Table S3: Association between differentially expressed genes and three ICS response miRs using a linear mixed model. Table S4: Details of genes from black and magenta modules. Table S5: List of unique validated gene targets of three ICS-response-associated miRs.