Circulating MicroRNAs associated with Bronchodilator Response in Childhood Asthma

Rationale: Bronchodilator response (BDR) is a measure of improvement in airway smooth muscle tone, inhibition of liquid accumulation and mucus section into the lumen in response to short-acting beta-2 agonists that varies among asthmatic patients. MicroRNAs (miRNAs) are well-known post-translational regulators. Identifying miRNAs associated with BDR could lead to a better understanding of the underlying complex pathophysiology. Objective: The purpose of this study is to identify circulating miRNAs associated with bronchodilator response in asthma and decipher possible mechanism of bronchodilator response variation. Methods: We used available small RNA sequencing on blood serum from 1,134 asthmatic children aged 6 to 14 years who participated in the Genetics of Asthma in Costa Rica Study (GACRS). We filtered the participants into high and low bronchodilator response (BDR) quartiles and used DeSeq2 to identify miRNAs with differential expression (DE) in high (N= 277) vs low (N= 278) BDR group. Replication was carried out in the Leukotriene modifier Or Corticosteroids or Corticosteroid-Salmeterol trial (LOCCS), an adult asthma cohort. The putative target genes of DE miRNAs were identified, and pathway enrichment analysis was performed. Results: We identified 10 down-regulated miRNAs having odds ratios (OR) between 0.37 and 0.76 for a doubling of miRNA counts and one up-regulated miRNA (OR=2.26) between high and low BDR group. These were assessed for replication in the LOCCS cohort, where two miRNAs (miR-200b-3p and miR-1246) were associated. Further, functional annotation of 11 DE miRNAs were performed as well as of two replicated miRs. Target genes of these miRs were enriched in regulation of cholesterol biosynthesis by SREBPs, ESR-mediated signaling, G1/S transition, RHO GTPase cycle, and signaling by TGFB family pathways. Conclusion: MiRNAs miR-1246 and miR-200b-3p are associated with both childhood and adult asthma BDR. Our findings add to the growing body of evidence that miRNAs play a significant role in the difference of asthma treatment response among patients as it points to genomic regulatory machinery underlying difference in bronchodilator response among patients. Trial registration: LOCCS cohort [ClinicalTrials.gov number: NCT00156819], GACRS cohort [ClinicalTrials.gov number: NCT00021840]


Background
Asthma is a heterogeneous disease that affects over 300 million people worldwide (1). It is characterized by chronic airway in ammation and a clinical history of wheezing, coughing, chest tightness, and shortness of breath that varies with time and intensity, as well as expiratory air ow limitation that is variably reversible with inhaled bronchodilators (2).
Albuterol, a short-acting beta-2 agonist (SABA) and bronchodilator, is one of the most used asthma medications in both children and adults (3). β2-agonists promote bronchodilation by activating β2-adrenergic receptors (β2ARs) on airway smooth muscle cells, resulting in decreased bronchoconstriction via increases in cyclic adenosine monophosphate (cAMP) and protein kinase A (PKA) (4). This is a physiological response that involves interaction between several cells and tissues such as in ammatory (5) and various airway cell types: the epithelium (6); smooth muscle cells (7); and cells of the autonomic nervous system (8).
Bronchodilator response (BDR) assesses the difference in FEV1 before and after the administration of a SABA. SABAs have variable e cacy among patients and BDR testing can be used to assess such effectiveness (9). Thus, studying BDR may provide insight into both the pathophysiology and pharmacogenetics of asthma.
High variability in BDR among individuals and populations has been described and is in part due to environmental and genetic factors (10)(11)(12). Estimates of BDR heritability range between 31-92% (13)(14)(15), and genome-wide association studies (GWASs) and whole-genome sequencing studies have identi ed susceptibility genes for BDR in subjects with asthma. Such genes include those for the β2-adrenergic receptor (ADRB2) (16), adenylyl cyclase type 9 (ADCY9) (17), corticotrophin-releasing hormone receptor 2 (CRHR2) (18), arginase 1 (ARG1) (19), and Spermatogenesis Associated Serine Rich 2 Like (SPATS2L) (20). These also highlighted several biological pathways that are likely to be involved in BDR control (e.g., Erk1/2 signal transduction, PI3K/Akt signal transduction, and nitric-oxide (NO) signaling pathway) (21,22). Previously, pharmaco-metabolomics of bronchodilator response in asthma were also studied (23,24). MicroRNAs (miRNAs) are small non-coding RNA molecules that have sizes ranging from 18 to 22 nucleotides which act as posttranscriptional regulators of target gene expression (25) and have emerged as key regulators of epithelial cell and in ammatory processes (26). Circulating miRNAs have been suggested as biomarkers for a number of conditions, and have been shown to be important in a number of in ammatory-mediated processes (27,28), including asthma (29). We hypothesized that miRNAs regulate bronchodilator response, a key intermediate phenotype of asthma. To test this hypothesis, we examined the relation between serum miRNAs and bronchodilator response in children with asthma and then attempted to replicate the signi cant ndings in a cohort of adults with asthma.

Study Population
Subject recruitment and study procedures for the Genetics of Asthma in Costa Rica Study (GACRS) have been described in detail elsewhere (30,31). In brief, the GACRS included 1165 Costa Rican children with asthma aged 6 to 14 years who were recruited between February 2001 and July 2011. Asthma was de ned as physician-diagnosed asthma and having either at least two respiratory symptoms (wheezing, coughing, or dyspnea) or a history of asthma exacerbations in the previous year. Further, all participants had a high probability of having at least six great-grandparents born in Costa Rica's Central Valley, as determined by a genealogist based on each of the child's parents' paternal and maternal last names. Participants in the study completed a protocol that included a questionnaire on respiratory and general health that was slightly modi ed from one used in the Collaborative Study on the Genetics of Asthma (32

Replication Population
The Leukotriene modi er Or Corticosteroids or Corticosteroid-Salmeterol trial (LOCCS) (ClinicaTrials.gov-NCT00156819) has been previously described in detail (33). In summary, LOCCS enrolled 500 individuals with mild asthma to nd the best step-down therapy for those who were well-controlled on low-dose ICS. The trial was held between 2003 and 2005. Fluticasone 100 g twice daily or uticasone/salmeterol 100 g/50 g once daily provided better asthma control than montelukast alone, as measured by fewer treatment failures, fewer nocturnal awakenings, improved lung function, and higher asthma control questionnaire (ACQ) scores. The treatment was given in a double-blind fashion for 16 weeks. The time to treatment failure was the primary outcome. FEV1 and Bronchodilator assessment were measured at enrollment during the run-in period of ICS treatment. Participants in the LOCCS were mostly white, but there were a few participants from other racial or ethnic groups (e.g., Black, and Asian [see Table 1]).

Primary Outcome
Bronchodilator response (BDR) testing was performed according to American Thoracic Society criteria (34). BDR was calculated as the percent change in FEV1 in response to administration of 200 µg of inhaled albuterol, as (([post-BD FEV1 -pre-BD FEV1]/pre-BD FEV1) x 100). Percent-predicted FEV1 (ppFEV1) was computed using expected FEV1 formulae for age, sex, height, and race according to Hankinson et al. (35).

Sample Sequencing and Quality Control
We performed small RNA sequencing on serum from 1,134 GACRS samples of asthma patients and 390 samples of asthma patients in LOCCS. Both cohorts were sequenced following the same protocols (36). In brief, small RNA-seq libraries were prepared by using the Norgen Biotek Small RNA Library Prep Kit (Norgen Biotek, Therold, Canada) and sequenced on the Illumina NextSeq 500 platform. The ExceRpt pipeline was used for quality control (QC) of the RNA-seq data (37). The samples with less than 100k mapped reads were removed. miRNAs with less than ve mapped reads in at least 50% of subjects were removed. We used the guided Principal Component Analysis (gPCA) (38) package for the identi cation of batch effects in GACRS and LOCCS.

Identi cation of Differentially Expressed miRNAs and Statistical Approach
To decrease the intrinsic difference of BDR, high versus lowest quartiles of BDR were considered (named as high and low BDR group) and identi ed differentially expressed miRNAs between the high versus low BDR group using DESeq2 (39), which uses negative binomial regression, with a Benjamini-Hochberg false discovery rate (FDR) correction for multiple testing. A signi cance threshold of 10% FDR was used. The analysis was performed with adjustment for age, sex, use of inhaled corticosteroids (ICS) in the previous year and, baseline (pre-bronchodilator) percent predicted FEV1. Negative binomial regression was used to obtain estimates of effect size (betas and Odds Ratios) for a doubling of miRNA counts.
Top DE miRNAs were assessed for association with LOCCS high vs low BDR using DESeq2 and adjusted for age, sex, race, and baseline percent predicted FEV1.
Clinical and demographic features were compared using a Chi-square test for dichotomous variables and a t-test for continuous variables.

Functional Annotation of Differentially Expressed miRNAs
Putative target mRNA transcripts were identi ed for 11 DE miRNAs between high and low BDR group and two replicated miRNAs separately using the miRecords version 4 (40), TarBase   Similar trends were seen in the LOCCS replication cohort (Table 1B): both the groups were similar in terms of gender and race distribution but there was a signi cant (p-value = 0.01) trend in age distribution. The high BDR response group participants were younger than the low BDR response group. In this cohort also, the participants in high-response group had lower baseline ppFEV1 as compared to low-response group (Table 1B) though this was not signi cant for blood eosinophil count or airway hyperresponsiveness.  Fig. 2. These 11 miRNAs were tested for differential expression in two BDR-response groups in another independent asthma cohort of adults and identi ed two miRNAs (miR-1246 and miR-200b-3p) with same direction of effect at p-value < 0.05 (Table 2B).  Identi cation of Putative Targets and Functional Assessment of Differentially Expressed miRNAs We found 6,245 putative target mRNA transcripts for 11 DE miRNAs that were reported in databases using experimental approaches and performed functional pathway enrichment analysis of Reactome biological pathways (44) using the clusterPro ler R package (45). We also performed functional pathway enrichment analysis on each of the DE miRNA's targets separately. The top 30 enriched pathways are shown in Fig. 3. Regulation of cholesterol biosynthesis by SREBPs, ESR-mediated signaling, G1/S transition, RHO GTPase cycle, and Signaling by TGFB family pathways were among the top enriched pathways. It was also observed that miR-200b-3p target genes were also enriched in these pathways (Fig. 4).

Discussion
Asthma is a common disease of the airway, causing expiratory air ow limitation that is partially reversible with inhaled bronchodilators. Bronchodilators are the rst-line treatment for asthma, and they work by acting on beta-2-adrenergic receptors on airway smooth muscle (ASM) cells in the lower respiratory tract, allowing muscle relaxation and bronchodilation (46) as well as inhibits liquid accumulation and mucus section into the lumen (47)(48)(49). In this study, we tried to identify serum miRNAs as indicator of BDR in a cohort of children with asthma (GACRS) followed by replication study in an adult asthma (LOCCS) cohort. We found differences in ICS use, baseline FEV1, and PD20 among participants with high vs. low BDR in the GACRS cohort. Sometimes, BDR is greater in patients with lower starting, since they have more to gain from a bronchodilator FEV1 (regression beta = -0.41, pvalue < 10 − 15 ). We have attempted to correct for this effect by including baseline FEV1 as a covariate in our analysis, so that miRNAs associated with BDR should be more indicative of the airway's plasticity rather than the magnitude of lung function de cit. ICS use increases pre-BDR FEV1 (50,51) and would then decrease BDR (50), however, we noticed that ICS use was higher in patients with high BDR group. This may be due to confounding by indication, with patients with more serious disease and lower lung function requiring ICS therapy. Although PD20 differed by BDR response group, we did not include this as a covariate since it was anti-correlated with BDR, and inclusion in our model would result in decreased power.
We performed differential expression analysis using DeSeq2 to identify miRNAs associated with high vs low BDR and adjusting model for age, sex, use of inhaled corticosteroids (ICS) in the previous year and, baseline (pre-BDR) percent predicted FEV1. We found 11 miRNAs signi cantly associated with high vs low BDR in a study of Costa Rican children with asthma. In subjects with a high bronchodilator response, 10 of these 11 miRNAs were down-regulated, while one was up-regulated. Two of these miRNAs (miR-1246 and miR-200b-3p) were validated as being signi cant in the LOCCS cohort and regulated in the same direction, i.e., miR-1246 was down-regulated and miR-200b-3p was up-regulated in subjects with high BDR. These 11 miRNAs, and particularly miRs 200b-3p and 1246, may be indicative of general biological state that leads to difference in BDR.
The two replicated miRNAs were previously reported as potential biomarkers for respiratory diseases. miR-1246 has been reported to predict response to benralizumab in severe eosinophilic asthma (52), to distinguish healthy subjects from those with asthma (53), and to differentiate asthma from COPD or asthma-COPD overlap (ACO) along with two other miRNAs in a logistic regression model (54). Over-expression of miR-200-3p has been shown to reduce airway in ammation, mucus hypersecretion, and remodeling in asthma (55). Another study also suggested adenosine to inosine (A-to-I) edited sites in miR-200-3p in lower airway cells is associated with moderate-to-severe asthma (56). The putative target identi cation of these miRNAs revealed that miR-200b-3p regulates the expression of SPATS2L, a gene that was previously reported as a BDR gene (20) and miR-1246 regulates ADCY9, another BDR gene (17,57). DIANA-miTED: a microRNA tissue expression database (58) also shows that these replicated miRNAs, namely miR-200b-3p (42986.7 RPM) and miR-1246 (147.7 RPM), are expressed in the bronchus.
Both gene targets of all 11 DE miRNAs in aggregate and gene targets of the two replicated miRNAs separately were enriched in regulation of cholesterol biosynthesis by SREBPs, ESR-mediated signaling, G1/S transition, RHO GTPase cycle, and signaling by TGFB family pathways (Figs. 3 & 4).
Regulation of cholesterol biosynthesis by the SREBPs pathway promotes cholesterol accumulation through uptake (low-density lipoprotein receptor) and synthesis (e.g., hydroxymethylglutaryl coenzyme A reductase) in macrophages and other cells (59).
Recent ndings indicate that cholesterol tra cking and in ammation are associated in the lung (60-63). In the present study we found that the target genes of DE miRNAs were enriched in this pathway, which may indicate the role of BDR-responsive miRNAs in cholesterol tra cking and in ammatory response in asthma. This is of interest as there are studies that link BDR to the presence of in ammation.
RHO GTPase pathway is known to regulate many essential cellular processes, including actin dynamics, gene transcription, cellcycle progression and cell adhesion (64). We found that miR-200b-3p regulates the expression of ROCK2 gene that encodes Rhokinase, known to play a role in regulating mucus overproduction (65), airway smooth muscle (ASM) tone (66) and ASM cytoskeletal stiffness (67). Further, ROCK2 expression is increased in ASM and pulmonary blood vessels in human asthma (68). This indicates a possible role of miR-200b-3p in regulating bronchoconstriction. Further exploration of the mechanisms by which miR-200b-3p and its target gene ROCK2 affect BDR may be worth pursuing.
TGFB family pathway is known to play a role in epithelial shedding, mucus hyper-secretion, angiogenesis, airway hyperresponsiveness, ASMC hypertrophy and hyperplasia in an asthmatic mouse model (69-71). Previously, it has been reported that eosinophils constitute a major source of TGF-β in asthmatic airways (72,73). In this study, participants with a high bronchodilator response had a higher eosinophil count than those with a low bronchodilator response (Table 1). Previous investigations suggest that TGF-β1 may play a role in the development of resistance to bronchodilators in asthma by reducing the e cacy of β2-agonists and by inducing PDE4D gene expression in a Smad2/3-dependent pathway manner (74)(75)(76)(77). The DE miRNAs (miR-26b-5p, miR-378a-3p, miR-378i, miR-200b-3p, and miR-885-5p) were found to regulates the expression of TGFB1, TGFB2, TGFBR1, TGFBR2, and TGFBR3 genes encoding TGF-β and TGF-β receptor (Supplementary Table S1). Additionally, the target genes of DE miRNAs were found to be enriched in TGFB pathway and pathway associated with downregulation of SMAD2/3: SMAD4 transcriptional activity (Fig. 3), showing possible role of DE miRNAs in regulation of TGF-β associated pathway and thus involved in smooth muscle remodeling (Fig. 5). TGF-β works upstream of the RHO GTPase pathway, TGF-β activates RhoA/rho-associated protein kinase (ROCK), and the cross-talk between these two pathways promotes airway remodeling (78) and mucus formation (69).
Strengths of our study include leveraging a large cohort of childhood asthmatics with circulating miRNA sequencing data and careful spirometric evaluation. That two of the identi ed miRNAs were able to replicate in an adult asthma population, despite etiological differences between childhood and adult asthma, gives weight to their importance in determining e cacy of SABAs as rescue inhalers. Weaknesses of our study include the retrospective study design and inability to assess miRNA differences in airway smooth muscle cells. We anticipate that future work into ASM miRNAs would provide additional biological insight into differences in BDR.

Conclusion
In summary, we have identi ed differential expression of 11 miRNAs by bronchodilator response in children with asthma, that these miRNAs in uence biological pathways associated with in ammatory response, airway smooth muscle cell contraction and airway remodeling, and that two of these miRNAs were replicated in another cohort of adults with asthma. Our ndings add to the growing body of evidence that miRNAs play an important role in asthma treatment response differences among patients.

Declarations
Ethics approval and consent to participate The study was approved by the Institutional Review Boards of the Hospital Nacional de Niños (San José, Costa Rica) and Brigham and Women's Hospital (Boston, MA). The current work is covered by Brigham and Women's Hospital IRB# 2017P001799. All participants provided written informed consent to take part in the study.

Consent for publication
Not applicable

Availability of data and materials
Sequencing data is available in Gene Expression Omnibus (accession pending).

Competing interests
The authors declare no con ict of interest in relation to this manuscript.

Funding
This work was supported by R01 HL139634, R01 HL162570, R01 HL161362, and R01 HL127332.  Volcano plot for differential expression of miRNA between high and low BDR in the GACRS. Y-axis: represents the multiple testing corrected (FDR) p-value.

Figure 2
Clustered heatmap of all 11 differentially expressed miRs in the GACRS across conditions.DESeq2 normalized expression counts with shifted logarithm transformation was used. The heat map was created using unsupervised hierarchical clustering, and the distance metric was Pearson correlation. *Marked miRNAs were replicated in the LOCCS cohort. Reactome pathways enriched for 11 DE miRNAs in the GACRS at 5% FDR cut-off. The target genes were identi ed using Micro T-CDS, TarBase, and Target Scan databases. The pairwise similarities of the enriched terms calculated by the pairwise_termsim function using Jaccard's similarity index (JC) and the agglomeration method ward.Din is used for clustering in R. If a pathway was found to be enriched with a speci c DE miRNA's target gene, the DE miRNA name is written next to it. Figure 4 miRNA-target gene network between two replicated miRNAs. Nodes with different colors represent the genes in selected Reactome pathways.