Skip to main content

Epigenomic response to albuterol treatment in asthma-relevant airway epithelial cells



Albuterol is the first-line asthma medication used in diverse populations. Although DNA methylation (DNAm) is an epigenetic mechanism involved in asthma and bronchodilator drug response (BDR), no study has assessed whether albuterol could induce changes in the airway epithelial methylome. We aimed to characterize albuterol-induced DNAm changes in airway epithelial cells, and assess potential functional consequences and the influence of genetic variation and asthma-related clinical variables.


We followed a discovery and validation study design to characterize albuterol-induced DNAm changes in paired airway epithelial cultures stimulated in vitro with albuterol. In the discovery phase, an epigenome-wide association study using paired nasal epithelial cultures from Puerto Rican children (n = 97) identified 22 CpGs genome-wide associated with repeated-use albuterol treatment (p < 9 × 10–8). Albuterol predominantly induced a hypomethylation effect on CpGs captured by the EPIC array across the genome (probability of hypomethylation: 76%, p value = 3.3 × 10–5). DNAm changes on the CpGs cg23032799 (CREB3L1), cg00483640 (MYLK4-LINC01600), and cg05673431 (KSR1) were validated in nasal epithelia from 10 independent donors (false discovery rate [FDR] < 0.05). The effect on the CpG cg23032799 (CREB3L1) was cross-tissue validated in bronchial epithelial cells at nominal level (p = 0.030). DNAm changes in these three CpGs were shown to be influenced by three independent genetic variants (FDR < 0.05). In silico analyses showed these polymorphisms regulated gene expression of nearby genes in lungs and/or fibroblasts including KSR1 and LINC01600 (6.30 × 10–14 ≤ p ≤ 6.60 × 10–5). Additionally, hypomethylation at the CpGs cg10290200 (FLNC) and cg05673431 (KSR1) was associated with increased gene expression of the genes where they are located (FDR < 0.05). Furthermore, while the epigenetic effect of albuterol was independent of the asthma status, severity, and use of medication, BDR was nominally associated with the effect on the CpG cg23032799 (CREB3L1) (p = 0.004). Gene-set enrichment analyses revealed that epigenomic modifications of albuterol could participate in asthma-relevant processes (e.g., IL-2, TNF-α, and NF-κB signaling pathways). Finally, nine differentially methylated regions were associated with albuterol treatment, including CREB3L1, MYLK4, and KSR1 (adjusted p value < 0.05).


This study revealed evidence of epigenetic modifications induced by albuterol in the mucociliary airway epithelium. The epigenomic response induced by albuterol might have potential clinical implications by affecting biological pathways relevant to asthma.


Asthma is a leading worldwide biomedical concern, affecting up to 300 million people and related to one of each 250 deaths worldwide [1]. Treatment response to current therapies is heterogeneous, with 10% of patients not responding to available therapies in whom severe asthma remains unresolved [1]. It carries a high burden for patients, families, and healthcare systems, since non-responders suffer frequently from asthma exacerbations, have a lower quality of life, and increased mortality rates [1]. The impact of asthma across different regions and populations is not homogeneous [2], and low-income populations, such as Puerto Ricans, have reported the highest asthma prevalence, impaired lung function, and mortality [3]. Albuterol, a short-acting β2-agonist bronchodilator, is the most widely used drug to relieve asthma symptoms and the current first-line asthma drug in these populations. However, albuterol effectiveness, measured by the bronchodilator drug response (BDR), is the lowest in African Americans and Puerto Ricans in comparison with other populations such as Mexican Americans [4,5,6]. Despite this, these populations are prone to use albuterol exclusively in the management of their asthma, rather than more expensive maintenance medications, including inhaled corticosteroids, leading to the overuse of these short-acting bronchodilators. The molecular effects of this overuse on the airway epithelium which lines the airway lumen are unknown.

Genetics and environmental factors contribute to asthma pathophysiology. Since 2007, genome-wide association studies (GWAS) have provided new insights into the genetics of asthma and BDR [4,5,6,7,8]. However, identified genetic variants are unable to explain the heritability of asthma susceptibility and treatment response, estimated at over 60–65% [7, 9]. Epigenetic modifications regulate gene expression without modifying the underlying DNA sequence, being DNA methylation (DNAm) the most widely studied epigenetic marker [10]. Elucidating the role of epigenetic mechanisms in asthma has been of great interest since they might mediate the contribution of genetics and the environment to the disease [7, 10]. In epigenetic studies, nasal samples are one of the most commonly used as non-invasive proxies of the lung environment [11], as the nasal airway transcriptome has been demonstrated to reflect asthma status and is highly correlated with that of the bronchial airway cells [12].

Epigenome-wide association studies (EWAS) have sought epigenetic biomarkers to explain missing asthma heritability [7, 10, 11]. Epigenetics has also been implicated in lung function and albuterol response in children [13,14,15]. Our group recently identified DNAm markers associated with BDR in African American and Latino children with asthma that interplay with genetic variation and gene expression, with potential applicability to BDR classification [14]. Although epigenetic changes are modifiable, there is a lack of studies evaluating the impact of asthma therapies on DNAm. To our best, few studies have assessed the potential epigenomic response of corticosteroids [16, 17], but none have focused on albuterol.

We hypothesized that albuterol induces genome-wide DNAm changes in the airway epithelium and that repeated exposures to albuterol will durably reprogram the molecular response of these cells. This study aimed to characterize the effect of albuterol on DNAm on the mucociliary airway epithelium, and to assess the influence of genetic variation and asthma-related clinical characteristics (i.e., asthma severity, medication use, and BDR), and potential functional consequences on gene expression and biological pathways. Epigenomic modifications were evaluated in a discovery and validation study design using in vitro models of paired primary nasal mucociliary epithelial cultures exposed or not to repeated doses of albuterol.


Airway epithelial model of albuterol treatment

To model human airway epithelial molecular responses to albuterol, we established an in vitro nasal airway epithelial culture system. In particular, we generated air–liquid interface (ALI), mucociliary airway epithelia from nasal brushing-derived basal stem cells (Fig. 1). To mimic chronic, repeated albuterol use, we stimulated paired donor cultures with either 100 μM albuterol or mock-stimulus, twice daily, for five consecutive days (Fig. 1). This stimulation course was not associated with visible toxicity (by light microscopy) or loss of epithelial barrier function.

Fig. 1
figure 1

Outline of the laboratory experiments design. A Air–liquid interface (ALI) mucociliary airway epithelia were generated from basal stem cells obtained from nasal brushing. Paired samples were stimulated with albuterol (Trx) and vehicle control (MeOH). B Experiment timeline. Basal airway epithelial cells were differentiated for 24 days (yellow), followed by 5 days of either control vehicle or albuterol stimulation (green). Cultures were harvested after 48 h of rest from the last stimulation (blue)

We leveraged our repeated albuterol stimulation model to evaluate epigenomic modifications induced by this treatment in a large cohort of children from the Genes-environments & Admixture in Latino Americans (GALA II) study with nasal airway epithelial brushings. After 48 h post-last exposure to the stimulation course, cultures were harvested for DNA extraction. This timing was selected to allow acute methylation changes induced by stimulation to recover, allowing more persistent stimulation-associated epigenetic changes to be monitored. Epigenetic changes were evaluated through genome-wide profiling of DNAm using the Infinium MethylationEPIC microarray (Illumina, San Diego, CA). The DNAm data quality control (QC) of samples and CpGs is summarized in Additional file 1: Table S1, remaining 97 samples and 689,483 probes after the QC. The characteristics of samples passing the QC are summarized in Table 1. Briefly, participants from GALA II were Puerto Rican children and young adults and included similar proportions of healthy individuals, patients with mild asthma, and moderate-to-severe asthma.

Table 1 Characteristics of the airway cell donors

Genome-wide DNAm changes induced by albuterol treatment in nasal epithelia

Albuterol-induced DNAm changes in nasal samples were assessed through a paired linear regression model. Analyses were corrected for cell-tissue heterogeneity, while all other potential confounders intrinsic to the sample were corrected by pairing (e.g., age, sex, asthma status, and ancestry). After bias and genomic inflation correction, the quantile–quantile (Q–Q) plots did not show evidence of genomic inflation (λ = 1.09; Additional file 1: Figure S1). A total of 66 CpG sites were associated with albuterol treatment with a false discovery rate (FDR) < 5%, and 22 CpGs surpassed the genome-wide significance threshold established for the EPIC array (p < 9 × 10–8) [18] (Fig. 2, Table 2, Additional file 1: Table S2). These included genes involved in inflammation (TNFRSF21, IFNGR1, CSF3, and LTA4H), interaction with the cytoskeleton (FLNC, MYLK4, KSR1, IPP, SPTAN1, MACF1, and HSD17B12), cell adhesion (TGFBI, PCDH12, EPCAM, FAT4, and SPTAN1), host defense against airway microbial infections (BPIFA1), and lipids metabolism and transportation (PPAP2B, PLA2G6, ATP8B1, GLTPD2, and GRAMD1B). Albuterol exerted predominantly a hypomethylation effect on CpG sites captured by the EPIC array across the genome (probability of hypomethylation: 0.76, 95% Confidence Interval [CI] 0.65–0.85, p value = 3.3 × 10–5) (Additional file 1: Figure S2).

Fig. 2
figure 2

Manhattan plot of the EWAS of DNAm changes induced by albuterol treatment in nasal epithelial cells. The blue and red lines represent the false discovery rate (FDR) < 5% and the genome-wide (p < 9 × 10–8) significance thresholds, respectively. Gene annotation is represented for genome-wide significant CpGs, highlighting in boldface those CpGs that were validated

Table 2 Summary results of genome-wide CpGs associated with albuterol exposure in nasal cells

Validation in independent samples and bronchial epithelia

Genome-wide significant CpGs were attempted for validation in an independent subset of nasal epithelial cells exposed under the same treatment conditions. The hypomethylation effect on three CpG sites annotated to the CREB3L1 (cg23032799), MYLK4-LINC01600 (cg00483640), and KSR1 (cg05673431) genes was validated with an FDR < 5% (Table 2). In addition, these three CpGs from the nasal validation subset donors were cross-tissue evaluated in bronchial epithelial cells. The CpG cg23032799 (CREB3L1) showed a similar hypomethylation effect in bronchial epithelial cells exposed to repeated doses of albuterol (logFC = − 0.169, p = 0.030) (Additional file 1: Table S3).

Evaluation of asthma-related conditions on albuterol-induced DNAm changes

We examined whether the albuterol-induced DNAm changes on CpGs could be partially conditioned by asthma status, previous use of asthma medications, or BDR during sample collection. First, a meta-analysis of stratified analysis based on asthma status (non-asthma, mild asthma, and moderate-to-severe asthma) did not report significant evidence of heterogeneity among groups (0.058 ≤ Cochran’s Q p value ≤ 0.916 and 0 ≤ I2 ≤ 64.94, Additional file 1: Table S4 and Figure S3). Similarly, stratified analyses based on the previous use of (i) albuterol or (ii) any controller medication, showed only one and two CpGs, respectively, out of the 22 evaluated with a significant heterogeneous effect according to previous medication use (Cochran’s Q p value < 0.05, Additional file 1: Tables S5, S6 and Figure S3).

Second, we applied a linear regression model to examine the effects of basal BDR on the albuterol-induced changes on DNAm in cell assays. Relative changes in DNAm induced by albuterol were estimated as follows: ΔDNAm = (DNAmalbuterol−DNAmcontrol)/DNAmcontrol. Only asthma subjects had available BDR data. After adjusting for age, sex, ancestry, and tissue heterogeneity, DNAm changes on the CpG cg23032799 (CREB3L1) were associated with basal BDR. Specifically, lower basal BDR was associated with higher hypomethylation effect induced by albuterol (coefficient = 0.032, p = 0.004). However, this association did not remain significant after multiple comparisons adjustment (FDR = 0.09) (Additional file 1: Table S7).

Identification of SNPs regulating changes on DNAm and gene expression

Whole-genome sequencing (WGS) data were available for individuals from the GALA II study [19]. A cis-methylation quantitative trait locus (meQTL) analysis was conducted to identify single nucleotide polymorphisms (SNPs) involved in the genetic susceptibility of albuterol-induced DNAm changes. Normalized ΔDNAm was included as the dependent variable and the genotypes of common SNPs as predictors. meQTL analyses were corrected for age, sex, ancestry, tissue heterogeneity, and asthma status. Only the three genome-wide associated CpGs validated in independent samples were associated with multiple meQTLs after multiple comparisons correction (FDR < 0.05, Table 3, Fig. 3, Additional file 1: Table S8).

Table 3 Summary results of independent meQTLs
Fig. 3
figure 3

Violin plots of independent meQTLs in nasal samples identified for the genome-wide associated CpGs that were validated: A cg23032799 (CREB3L1), B cg00483640 (MYLK4-LINC01600), and C cg05673431 (KSR1). The normalized ΔDNAm values (y-axis) are plotted against the genotypes for each independent SNP

The functional consequences of independent meQTLs on gene expression in asthma-related tissues were evaluated through an in silico expression quantitative trait locus (eQTL) analysis. The SNP rs11038897, which regulates the DNAm changes on cg23032799 (CREB3L1), was found to regulate the gene expression of multiple genes in fibroblasts (ATG13, MADD, and C11orf49) and lung tissue (ATG13, MDK, and C11orf49). Additionally, the SNPs rs11038897 and rs6505279, which are meQTLs of the CpGs cg00483640 and cg05673431, respectively, regulate the expression of the genes where they are located in fibroblasts (LINC01600) and lung tissue (KSR1) (Additional file 1: Table S9).

Functional consequences of DNAm changes on gene expression

Bulk RNA-seq data were generated in all paired samples from the discovery phase. An expression quantitative trait methylation (eQTM) analysis was carried out to test for the association between changes in DNAm and gene expression of nearby genes to each CpG. Relative changes in gene expression data were estimated as ΔRNA-seq = (RNA-seqalbuterol−RNA-seqcontrol)/RNA-seqcontrol. The association between normalized ΔRNA-seq and ΔDNAm was tested through linear regression models adjusted for age, sex, asthma, ancestry, tissue heterogeneity, and batch effect from gene expression data. After multiple comparisons adjustment, we identified that hypomethylation at two CpGs was associated with increased gene expression patterns in the genes where they are located (Additional file 1: Table S10). Specifically, the CpGs cg10290200 and cg05673431 were associated with FLNC (coefficient: − 0.42, p = 0.001, FDR = 0.009) and KSR1 (coefficient: -0.30, p = 0.010, FDR = 0.029) gene expression, respectively. The CpG cg05673431 was also associated with increased expression of LGALS9 (coefficient: 0.27, p = 0.026, FDR = 0.039).

Gene-set enrichment analysis

We performed a gene-set enrichment analysis (GSEA) to identify biological pathways, molecular mechanisms, and drug signatures related to the epigenetic changes induced by albuterol and its functional consequences. Those CpGs that surpassed a p < 5 × 10–4 threshold in the main EWAS were selected to be annotated to genes, which were analyzed in the GSEA. As described in the Methods section, we only retained the results that remained significant after multiple comparison corrections (FDR < 0.05) and showed to be robust to varying the input p value thresholds for selecting the CpGs. We observed enrichment in multiple biological pathways and gene ontologies potentially related to asthma and the mechanism of action of albuterol. These included the tumor necrosis factor-alpha (TNF-α) signaling pathway mediated by the nuclear factor kappa B (NF-kB) transcription factor, interleukin-2 (IL-2) signaling pathway, regulation of cell proliferation, and components of the actin cytoskeleton. Additionally, we reported that genes affected by albuterol through DNAm changes are more likely to be regulated by numerous compounds and drugs than chance, including estradiol and trichostatin A (TSA) (Fig. 4, Additional file 1: Table S11).

Fig. 4
figure 4

Bar plots of the significant results identified in the enrichment analyses using the Enrichr tool. The −log10(q-value) for each term is represented in the x-axis. Only the top-three significant terms of the drug signature enrichment analyses are reported.

Regional DNAm changes induced by albuterol

We examined whether repeated stimulations with albuterol induce changes in DNAm on genomic regions. Differentially methylated regions (DMRs) were estimated using two independent software to minimize false positive results. Regional DNAm changes induced by albuterol were identified in nine DMRs (adjusted p value < 0.05, Table 4, Additional file 1: Tables S12, S13). The top hit was a genomic region of 103 bp annotated to the GLTPD2 gene (adjusted p value = 3.75 × 10–12). We also identified DMRs annotated to CREB3L1 (adjusted p value = 9.81 × 10–10), MYLK4 (adjusted p value = 2.04 × 10–7), and KSR1 (adjusted p value = 0.007).

Table 4 DMRs associated with albuterol treatment in nasal samples


In this study, we identified for the first time epigenomic changes induced by albuterol in mucociliary airway epithelia. Repeated stimulations with albuterol induced an overall genome-wide hypomethylation on CpGs captured by the EPIC array in nasal epithelia. Among 22 genome-wide significant CpGs, the effects on three CpGs annotated to CREB3L1, MYLK4-LINC01600, and KSR1 genes were validated in independent nasal samples. Albuterol induced regional DNAm changes in at least nine DMRs, including CREB3L1, MYLK4-LINC01600, and KSR1. The hypomethylation effect of the CpG near CREB3L1 was cross-tissue validated in bronchial epithelia and suggestively associated with BDR. Specifically, we observed that the hypomethylation effect on CREB3L1 was inversely associated with BDR, suggesting that albuterol might have further functional consequences in the gene regulatory landscape in poor responders to bronchodilators. We also reported meQTLs associated with the genetic susceptibility to these albuterol-induced DNAm changes, which regulate gene expression of nearby genes in lungs and/or fibroblasts. Additionally, hypomethylation on the CpGs cg10290200 and cg05673431 was associated with higher expression of FLNC and KSR1, respectively. Moreover, epigenetic loci modified by albuterol were located in genes enriched in asthma-associated processes (i.e., TNF-α, NF-kB, and IL-2 signaling pathways) and likely regulated by a potential asthma drug (i.e., TSA).

Albuterol triggers a rapid-onset bronchodilator effect in the airway smooth muscle and regulates airway inflammation by diverse pathways that have not been completely uncovered. The main known mechanism of action is mediated by its agonism on the β2-adrenergic receptor (ADRB2) [20]. Briefly, the activation of ADRB2 initiates a cascade of intracellular signals characterized by increasing cyclic adenosine monophosphate (cAMP) intracellular levels and resulting in the inhibition of intracellular Ca2+ via protein kinase A. The bronchodilator effect is a consequence of the inhibition of the Ca2+-dependent myosin light chain phosphorylation [20]. The main epigenetic markers reported in this study are annotated to genes potentially involved in this pathway. Our findings are relevant by providing new insights into potential mechanisms and epigenomic cell responses induced by albuterol that have not been described yet.

First, we reported albuterol-induced DNAm changes on CREB3L1 in nasal and bronchial epithelial cells, which are associated with genetic variation and the BDR of the donors. CREB3L1 encodes the cAMP response element binding protein 3-like 1, a member of the CREB/ATF transcription factor family. CREB3L1 is activated in response to increased cAMP levels and is involved in stress response, regulation of cell secretory capacity, extracellular matrix production, cell migration, and virus infection response [21,22,23,24]. Interestingly, we recently found CREB3L1 to be a transcription factor expressed in human mucus secretory cells [25]. Our finding of albuterol-induced downregulation in CREB3L1 methylation locus, suggests that airway secretory cell differentiation could be modulated by albuterol usage. DNAm on CREB3L1 participates in the toxic effects of inhaled silica nanoparticles on bronchial epithelial cells [26]. Although genetic variants in the intergenic region of PHF21A-CREB3L1 have been associated with post-bronchodilator lung function measurements [27], the potential role of DNAm has not been described. Nonetheless, the family-related CREB1 protein has been widely studied in the context of asthma. CREB1 participates in promoting epigenetic changes to switch on pro-inflammatory genes and is involved in persistent bronchial inflammation in asthma and the anti-inflammatory effect of inhaled corticosteroids [28]. Additionally, it participates in the regulation of ADRB2 expression and downregulates lung receptors after chronic exposure to β2-agonists [29]. CREB1 is one of the proteins involved in airway smooth muscle contraction mediated by cAMP and diacylglycerol kinase (DGK) [30]. Although the effects in vivo are uncertain, in vitro models have shown that β2-agonists increase CREB1 binding to the DNA in murine lung cells, while corticosteroids reduce it [31, 32]. Similarly, CREB3L1 expression is diminished in the rat hypothalamus after dexamethasone exposure, while increased levels of cAMP enhance its expression [33].

Second, we reported DNAm changes, influenced by genetic variation, on epigenetic markers annotated to MYLK4 and KSR1. Additionally, we also identified that albuterol induced hypomethylation on two CpGs associated with increased gene expression of FLNC and KSR1. MYLK4 encodes a member of the myosin light chain kinase (MYLK) family. MYLKs are involved in the phosphorylation of the myosin light chain and promote the contraction of smooth muscle cells leading to bronchoconstriction [34]. Although studies on MYLK4 are limited, MYLKs are known to participate in cell migration, invasion, and proliferation [34]. Genetic variation in MYLK1, the best-studied gene in the MYLK family, is associated with asthma risk and exacerbations and is a common genetic factor in diseases involving smooth muscle contraction and inflammation [35, 36]. The kinase activity of MYLK is activated by calmodulin when is bound to Ca2+ [37]. The kinase suppressor of ras 1 (KSR1) is a scaffold protein required for the association between Ca2+ and calmodulin, with a known regulatory effect on the mitogen-activated protein kinase (MAPK) cascade [38]. On the other hand, the filamin C (FLNC) is a protein that participates in the interaction between actin filaments and binding partners in muscle cells [39]. Although its role in the ADRB2 pathway has not been described, FLNC interacts with β-arrestin-2 [39], which is one of the main mediators involved in the side effects of bronchodilators via ADRB2 [20].

Third, enrichment analyses revealed that epigenomic modifications of albuterol could participate in IL-2, TNF-α, and NF-κB pathways, and that gene expression levels of affected genetic loci were more likely to be regulated by TSA and estradiol than expected by chance. IL-2 is a pro-inflammatory cytokine involved in the development of Treg cells and associated with eosinophil proliferation, airway inflammation and hyperreactivity, and impaired lung function in asthma patients [40,41,42]. IL-2 is involved in airway smooth muscle contraction by regulating the airway sensitization to leukotrienes [42], and IL-2-based therapies induce severe bronchoconstriction that can be reversible with albuterol under certain conditions [43]. Additionally, the effects of albuterol on airway inflammation might be mediated by IL-2 through the inhibition of TNF-α and NF-κB, and the activation of type 2 inflammation via STAT5 [44, 45]. NF-κB is a pro-inflammatory transcription factor widely involved in the pathogenesis of asthma [44] whose activation is inhibited by β2-agonists [46]. TNF-α is a cytokine with a regulatory role in inflammatory responses and the development of allergic diseases, particularly asthma [40]. TNF-α can activate NF-κB-mediated pathways and its expression is increased in patients with severe and steroid-refractory asthma, being considered a potential therapeutic target [47, 48]. Similarly, albuterol inhibits TNF-α expression, a mechanism involved in the agonism of the anti-inflammatory effect of corticosteroids and theophylline [49, 50]. Finally, TSA is an inhibitor of histone deacetylases that regulates inflammatory genes via NF-κB and has demonstrated in vivo anti-inflammatory properties in murine asthma models [51]. It has also been shown to reduce airway constriction by decreasing Ca2+ mobilization [52]. Additionally, previous studies have reinforced its potential therapeutic interest in asthma exacerbations by regulating genes implicated in ICS response and microbiome composition [53, 54]. On the other hand, estradiol is one of the main sex hormones in females, which has been suggested to be implicated in the sex differences in asthma throughout the lifespan [55]. Although this finding may suggest a differential in vitro effect of albuterol between cells derived from females and males, except for probe cg15833511, we did not observe evidence of significant heterogeneity in stratified analyses based on biological sex (Additional file 1: Table S14).

This study has several strengths. First, our experimental design allowed us to identify the causality between albuterol exposure and airway methylome alterations. The sample size in the discovery stage provided an 80% statistical power to identify modest DNAm changes (≥ 10%) with a p value < 1 × 10–6 (which approximately corresponds to a genome-wide significance at an FDR < 5%) [56]. Second, we applied a discovery, validation, and cross-tissue evaluation design that supported the robustness of our associations and their potential effects in different asthma-relevant tissues. Third, potential spurious results were minimized by correcting for bias and genomic inflation using a Bayesian method designed for EWAS and controlling for potential cell heterogeneity and unknown confounders. Fourth, we integrated epigenomic, transcriptomic, and WGS data mainly focusing on a diverse and underrepresented population (i.e., Puerto Ricans). Nevertheless, some limitations must be acknowledged. First, the sample size in the validation stage is limited, restricting the statistical power for the validation of markers identified in the discovery phase. Since independent donors were balanced for diverse relevant characteristics (i.e., age, sex, asthma status, body mass index, and ethnicity), we were unable to assess phenotype-specific effects with evidence of validation. Second, although it is recognized that the nasal airway reflects the lower airways, the effects on bronchial epithelial cells remained partially understudied in our study due to a limited sample size. Third, despite examining the potential influence of the prescription of asthma medications, the frequency of use and medication adherence were not available. Fourth, our study design did not allow us to identify potential clinical implications of epigenetic changes on treatment response, adverse reactions, and/or drug tolerance. Fifth, identified epigenetic markers in vitro may not reflect in vivo effects of albuterol when administered at a therapeutical dosage in asthma patients. Human organoids and in vivo animal models should allow further investigation of potential effects taking place in a therapeutic context. Sixth, the effect of albuterol on sex chromosomes DNAm remained unaddressed in this study.


We characterized for the first time the epigenomic response induced by albuterol treatment in airway epithelia. DNAm in 22 genome-wide significant CpGs and nine genomic regions were modified by albuterol in nasal epithelia from Puerto Rican children with different asthma statuses. Affected genes were enriched in asthma-relevant processes, including IL-2, TNF-α, and NF-κB signaling pathways, and were more likely to be regulated by Trichostatin A than chance. We observed evidence of genetic susceptibility to albuterol-induced changes and functional consequences on FLNC and KSR1 expression. Our findings provide novel insights into the potential role of CREB3L1, MYLK4, KSR1, and FLNC on epigenetic-mediated biological effects induced by albuterol.


Human subjects information

Basal airway epithelial cultures and then subsequently derived ALI cultures, used for the 100 donors' repeated-use albuterol stimulation model, were derived from nasal airway epithelial brushings collected from subjects recruited as part of the Genes-environments and Admixture in Latino Americans II (GALA II) childhood asthma study, which was approved by local institutional review boards (UCSF, IRB number 10–00889, Reference number 153543, NJH HS-2627). All subjects and their parents provided written informed assent and written informed consent, respectively [57]. Paired nasal and bronchial airway brushings used in generating the ALI albuterol validation study cultures were derived from subjects recruited as part of the Obese Asthma: Unveiling Metabolic and Behavioral Pathways study at CU-Anschutz and approved by local institutional review boards (University of Colorado, IRB Number 19-0510, 21-3959, and 16-2522, NJH HS-3110).

Albuterol treatment of mucociliary airway epithelial cultures

To model the in vivo airway epithelium, we generated ALI mucociliary epithelial cultures. For primary human nasal and bronchial epithelial samples, basal cells were cultured for expansion using a modified Schlegel method as previously described [58, 59]. Primary basal cells were seeded onto 6.5 mm, 0.4 μm pore transmembrane inserts (4 × 104 cells/insert) in PneumaCult Expansion Plus medium (StemCell Technologies) supplemented with Y-27632. Cultures were air-lifted upon confluence and basolateral media was then switched to PneumaCult ALI (PC-ALI; StemCell Technologies) media to stimulate mucociliary differentiation over the next 24 days.

Albuterol stimulations were completed on cultures from nasal samples from the GALA II study (n = 100 total; 30 and healthy controls, 33 with mild asthma, and 37 with moderate-to-severe asthma) and paired nasal and bronchial samples from the Obese Asthma study (n = 10 total; five healthy controls and five asthma patients). Exposures were performed starting at ALI D24 for epithelial cultures and stimulations were performed twice daily, once in the morning, and once in the evening for five consecutive days. For epithelial cultures, each morning fresh basolateral media was added to all ALI inserts, apical chambers were washed with warm PBS, and 20μl of ALI media supplemented with MeOH as vehicle control or ALI media supplemented with 100μM albuterol was added to the apical chamber of each ALI culture. For evening stimulations, apical chambers were washed once with warm PBS, and stimulations were performed as described above. PBS washes were performed prior to re-stimulations to prevent to accumulation of albuterol over time. Following 5 days of albuterol stimulations, ALI received fresh culture media without vehicle or albuterol and were allowed to rest for 48 h to allow for interrogation of epigenetic changes resulting from albuterol exposure. Samples were lysed using Zymo DNA/RNA Lysis Buffer (Zymo Research, Irvine, USA) supplemented with 40 mM DTT, and triplicates were pooled for DNA and RNA isolations as per manufacturer instructions.

Genome-wide DNAm profiling and quality control

Whole-genome DNAm data from 856,553 CpGs was generated using the Infinium MethylationEPIC microarray (Illumina) following the manufacturer's recommendations. QC of DNAm data was performed using the ENmix and ewastools R packages, as previously described [14, 60,61,62]. Briefly, samples with (1) a high proportion of bad-quality methylation data (> 5% of CpGs), (2) discordance between the reported sex and the epigenetic predicted sex, (3) a high missingness rate (> 5% of CpGs), and (4) potential cross-sample contaminated (based on the genotype distribution of control probes) were removed. Considering the paired approach of this study, samples were also removed if their paired sample was discarded. After QC, a total of 97 nasal cells were retained for the discovery stage, and all 10 nasal and bronchial epithelial cells were kept for the validation analyses. A QC of CpG sites was also conducted using the ENmix package [62]. Briefly, we corrected bisulfite intensities for background noise, dye and probe-type biases (oob and RELIC methods), and inter-array heterogeneity (quantile normalization) [63, 64]. Then, we excluded probes (1) with a high proportion of bad-quality methylation data (> 5% of samples), (2) located in sex chromosomes, (3) that non-specifically bind to a genomic position (cross-reactive) [65, 66], (4) with a multimodal distribution, and (5) potentially polymorphic probes. Probes with a multimodal distribution of beta values, which could capture other artifacts different than DNAm, were identified using the nmode function from ENmix [62]. Polymorphic probes were defined as probes containing a common SNP with a minor allele frequency (MAF) ≥ 5% at the CpG site (either at cytosine or guanine nucleotides) and, in the case of Type I probes, also those with a common SNP at the position where single base extension takes place during the methylation assay. We used available WGS data from GALA II to identify SNPs and estimate MAFs specifically in the analyzed individuals. In samples from independent donors (validation stage), we used the data from the Illumina manifest file v1.0 B4 [67]. DNAm was computed as beta values and transformed into M-values for better statistical performance. QC was conducted separately in the discovery and validation sets. Probes with a multimodal distribution were independently assessed for each cell type.

We used the ReFACTor algorithm as a free-reference method to control for unknown variation, mainly driven by cell heterogeneity [68]. The appropriate number of ReFACTor components to include in subsequent analyses (n = 3) was estimated for each dataset based on scree plots following the software developers’ recommendations. ReFACTor components were estimated with adjustment for cell condition (albuterol vs. control) and controlling for pairing. The sva R package was used to estimate surrogate variables for batch effect correction, but no significant variables in our dataset were detected [69].

Epigenome-wide association study

The discovery phase included the largest dataset consisting of 97 paired nasal epithelial cells from Puerto Rican children included in the GALA II study. We assessed the effect of albuterol treatment on whole-genome DNAm through paired robust linear regression models using the limma R package [70]. Models were corrected for ReFACTor principal components, while all other potential confounders intrinsic to the sample were corrected by pairing (e.g., age, sex, asthma status, ancestry). Considering that methylation levels were log-transformed into M-values, the log2(Fold-change) or logFC parameter from limma corresponds to the difference in the average M-value between control and treated samples. Negative and positive logFC values indicate the hypomethylation and hypermethylation effects of albuterol, respectively, and their values are correlated with the magnitude change in methylation values. We used the bacon R package to correct the effect sizes, standard errors, and p values controlling for bias and genomic inflation in epigenomic studies [71]. After correction, genomic inflation was further inspected by examining the Q–Q plots and the genomic inflation factor (λ). An FDR < 5% was used to declare significance and minimize false positive results, and the genome-wide significant threshold was established in p < 9 × 10–8, as previously estimated for the EPIC array [18]. Those significant CpGs with changes on DNAm (expressed as beta-values) < 2% were flagged as being potentially related to technical variations [72]. CpGs were annotated based on the Illumina manifest file v1.0 B4 and using the Genomic Regions Enrichment of Annotations Tool (GREAT) v4.0.4 [67, 73]. A one-sided binomial test was applied to assess whether albuterol showed predominantly a hypomethylation or hypermethylation effect. Genome-wide significant CpGs were followed-up in all subsequent analyses.

Validation in an independent subset of nasal and bronchial epithelial cells

DNAm changes induced by albuterol were attempted for validation in a subset of independent nasal epithelial cells from 10 donors treated under the same albuterol exposure conditions. An FDR < 5% was used to declare a significant association in this validation stage. Genome-wide CpGs that were validated were followed up in a cross-tissue evaluation analysis using bronchial epithelial cells. In these samples, we tested for the association of changes in DNAm and albuterol exposure following the same procedure as in nasal samples. A nominal p < 0.05 was used to declare a significant association.

Evaluation of asthma-related conditions on albuterol-induced DNAm changes

We examined whether albuterol-induced DNAm changes were conditioned by asthma status, biological sex, previous use of asthma medications, or BDR status during sample collection. Asthma status was classified according to the diagnosis and severity of the disease in non-asthma subjects, mild asthma patients, and moderate-to-severe asthma patients. Asthma severity was derived from the control of asthma symptoms (assessed using the Asthma Control Test and the Asthma Control Questionnaire) and the therapeutic step of each patient according to the Expert Panel Report-3 guidelines for diagnosis and management of asthma [74]. Mild intermittent and mild persistent asthma categories were classified as mild asthma, while moderate persistent and severe persistent asthma categories were grouped into moderate-to-severe asthma. Regarding asthma medications, we evaluated having a prescription of (1) short-acting beta-agonists or (2) any controller medication (i.e., inhaled corticosteroids, long-acting beta-agonists, combo medication, leukotriene receptor antagonists, oral corticosteroids, and/or theophylline). BDR was clinically characterized in asthma patients during sample collection as the maximum relative improvement in lung function (i.e., forced expiratory volume in the first second or FEV1) after albuterol inhalation [8].

We conducted stratified analyses according to asthma status, use of asthma medication, and biological sex and meta-analyzed them using METASOFT [75]. The heterogeneity among groups was examined using Cochran’s Q p value and I2. A Cochran’s Q p value < 0.05 or > 0.05 indicates the presence or absence of heterogeneity, respectively. The I2, which ranges between 0 and 100%, measures the proportion of variance as a consequence of the observed heterogeneity. The higher the I2 value, the greater the presence of heterogeneity.

On the other hand, linear regression models were applied to examine the effects of basal BDR on albuterol-induced DNAm changes in the cell assays. The relative change in DNAm induced by albuterol was estimated for each subject and CpG site as follows: ΔDNAm = (DNAmalbuterol-DNAmcontrol)/DNAmcontrol. The ΔDNAm was normalized by applying an inverse normal transformation to ensure a normal distribution. Considering the inter-individual comparison of this approach, models were corrected for age, sex, ancestry (the first three principal components), and cell heterogeneity (the first three ReFACTor components). ReFACTor components were re-estimated using the DNAm values from control samples adjusting for BDR, age, sex, and ancestry. Multiple comparisons were adjusted using an FDR < 0.05.

Methylation quantitative trait loci analyses

WGS data were generated for Puerto Ricans from GALA II within the Trans-Omics for Precision Medicine (TOPMed) Consortia at the New York Genome Center and Northwest Genomics Center. Briefly, paired-end reads (150 bp × 2) were generated on the HiSeq XTM Ten platform (Illumina) reaching a minimum mean genome coverage of 30x. All technical details are described on the TOPMed website [19]. The genetic susceptibility of albuterol-induced DNAm changes was investigated through a meQTL analysis. A total of 96 individuals with available epigenomic and WGS data were included in the analysis. For each CpG site, all SNPs located within a window of 500 kb (± 250 kb upstream and downstream) with a MAF ≥ 5% were evaluated in the cis-meQTL analysis. The association between normalized ΔDNAm and genotypes of these SNPs was tested through linear regression models using fastQTL that were adjusted for age, sex, ancestry, asthma status, and tissue heterogeneity [76]. An FDR < 0.05 at the probe level was used to correct for multiple comparisons considering that each CpG maps to multiple SNPs. Independent meQTL signals were identified by conditional regression models conditioned on the most significant meQTL using PLINK 1.9, as previously described [14]. Briefly, meQTL analyses were repeated adding the corresponding top-hit meQTL for each CpG as a covariate. All non-significant (p > 0.05) meQTLs were filtered out since they are dependent on the top-hit signal. If any meQTL remained significant, the new top-hit was considered as another independent meQTL and the process was repeated until no significant meQTLs were identified.

In silico expression quantitative trait loci analyses

The effects on gene expression of independent meQTL were inspected through in silico eQTL analyses using the public database Genotype-Tissue Expression (GTEx) portal. We examined whether these SNPs associate with the expression levels of nearby genes in different tissues related to the respiratory system. Summary statistics were extracted from the GTEx portal.

Expression quantitative trait methylation analyses

RNA normalization and library construction were performed using the KAPA mRNA Hyper Prep Kit (Roche, Basel, Switzerland), using the Beckman Coulter FXp automation system. Briefly, 200 ng of RNA was used as input, and Illumina Dual-Index adaptors from Integrated DNA Technologies were used to barcode libraries using 12 cycles of amplification. Paired-end sequences from pooled libraries were obtained using the NovaSeq 6000 platform (Illumina). Regarding the pre-processing of bulk RNA-seq data, raw reads were trimmed with from the BBMap package v38.79 with parameters “trimq = 10 qtrim = r minlength = 70 ref = adapters ktrim = r k = 23 mink = 11 hdist = 1 tpe tbo” [77]. Transcript counts were quantified with Salmon v1.5.2 with parameters “—seqBias —gcBias” [78] using human genome reference Gencode [] release 38 (GRCh38.p13). The Salmon output was imported into R and transcript-level counts aggregated to gene-level with tximport [79]. A total of 198 RNA-seq samples were used for analysis. One donor was removed due to the control and treatment labels being swapped. Genes were included if they were protein-coding and had an inferred count of at least 6 in at least 10% of the samples. Expression counts were transformed with the variance stabilizing transformation (VST) using DESeq2 v1.34.0 [80].

The association between albuterol-induced DNAm and gene expression changes was examined in a cis-eQTM analysis. Individuals with available DNAm and RNA-seq data (n = 95 paired samples) were included in the analysis. Changes in gene expression (ΔRNA-seq) were estimated and normalized following the same procedures as with DNAm data. We estimated surrogate variables of unknown sources of variation to correct for batch effect in RNA-seq data using the sva R package [69]. The significant number of surrogate variables was estimated using the svaseq function while adjusting for age, sex, ancestry, and tissue heterogeneity. To minimize overfitting in the eQTM analysis, only those surrogate variables not significantly correlated with any ReFACTor principal component were retained for the analyses (p > 0.05) (Additional file 1: Figure S4). Similar to ReFACTOr, surrogate variables were estimated using the gene expression data from the control samples. Then, we tested for the association between the normalized ΔDNAm of each CpG site and the normalized ΔRNA-seq of all genes whose transcription start site is located ± 250 kb upstream and downstream of each CpGs using MatrixEQTL [81]. Models were corrected for age, sex, asthma, ancestry, tissue heterogeneity, and RNA-seq batch effect. An FDR < 0.05 at the probe level was used to correct for multiple comparisons considering that each CpG maps to multiple genes.

Enrichment analyses

A GSEA was conducted to investigate whether the associated CpGs from the discovery EWAS were annotated to genes significantly involved in biological pathways, gene ontologies, traits, and drug signatures. Based on the EWAS summary statistics, a p < 5 × 10–4 threshold was used to select CpGs to include in the enrichment analyses. The GSEA was performed using the Enrichr database and an FDR < 0.05 was used to correct for multiple comparisons [82]. The GSEA was re-assessed by varying the p value thresholds for CpGs selection (i.e., p < 5 × 10–5 and p < 5 × 10–3). Only the terms that remained significant (p < 0.05) using these two alternative thresholds were retained to ensure the robustness of the findings. Analyses were carried out in January 2023, with the available version of the following gene-set libraries in Enrichr: BioPlanet 2019, Molecular Signatures Database (MSigDB), Gene Ontologies (GO) Biological Process 2019, GO Cellular Components 2019, Protein–Protein Interactions (PPI) Hub dataset, Transcription Factor PPIs, PheWeb 2019, and Drug Signature Database (DSigDB).

Differentially methylated regions

Regional DNAm analyses were conducted to identify DMRs associated with albuterol-induced changes in DNAm. We identified DMRs using two independent software: comb-p and DMRcate [83, 84]. DMRs were started at CpGs individually associated with p < 0.05 (comb-p) or FDR < 0.05 (DMRcate) and significant CpGs subsequently located in sliced windows of 1000 bp were included in the same region. All regions were restricted to contain at least two CpGs. In DMRcate the scaling factor for bandwidth was used as recommended by the developers (C = 2). Multiple comparisons were corrected intrinsically by each software using a Šidák-corrected p < 0.05 in comb-p and FDR < 0.05 in DMRcate. To minimize false positive results, only DMRs identified using the two software were retained. DMRs were annotated to the single nearest gene using GREAT v4.0.4 [73].

Availability of data and materials

All data necessary to evaluate the conclusions of this manuscript are reported in the main text and/or the supplementary information. Raw methylation data of samples included in the discovery phase (GALA II) are publicly available in the Gene Expression Omnibus (GEO) database under the accession number (GSE240155). TOPMed whole-genome sequencing data from GALA II are available in the database of Genotypes and Phenotypes (dbGaP) under accession number phs000920.v1.p1. The summary statistics of the epigenome-wide association study and the full report of the gene-set enrichment analysis are available in the Zenodo repository (


95% CI:

95% Confidence interval


β2-Adrenergic receptor


Air–liquid interface


Bronchodilator drug response


Cyclic adenosine monophosphate


Cytosine-phosphate-guanine dinucleotide


Diacylglycerol kinase


Differentially methylated region


DNA methylation


Expression quantitative trait locus


Expression quantitative trait methylation


Epigenome-wide association study


False discovery rate


Genes-environments and Admixture in Latino Americans


Gene-set enrichment analysis


Genotype-tissue expression portal


Genome-wide association study






Mitogen-activated protein kinase


Methylation quantitative trait locus


Myosin light chain kinase


Nuclear factor kappa B


Quality control

Q–Q plot:

Quantile–quantile plot


Single nucleotide polymorphism


Tumor necrosis factor-alpha


Trichostatin A


Whole-genome sequencing


  1. Trevor JL, Deshane JS. Refractory asthma: mechanisms, targets, and therapy. Allergy. 2014;69(7):817–27.

    Article  CAS  PubMed  Google Scholar 

  2. Global Asthma Network. The Global Asthma Report 2018 [Internet]. Auckland, New Zealand; 2018.

  3. Pino-Yanes M, Thakur N, Gignoux CR, Galanter JM, Roth LA, Eng C, et al. Genetic ancestry influences asthma susceptibility and lung function among Latinos. J Allergy Clin Immunol. 2015;135(1):228–35.

    Article  PubMed  Google Scholar 

  4. Drake KA, Torgerson DG, Gignoux CR, Galanter JM, Roth LA, Huntsman S, et al. A genome-wide association study of bronchodilator response in Latinos implicates rare variants. J Allergy Clin Immunol. 2014;133(2):370–8.

    Article  PubMed  Google Scholar 

  5. Naqvi M, Thyne S, Choudhry S, Tsai HJ, Navarro D, Castro RA, et al. Ethnic-specific differences in bronchodilator responsiveness among African Americans, Puerto Ricans, and Mexicans with asthma. J Asthma. 2007;44(8):639–48.

    Article  PubMed  Google Scholar 

  6. Spear ML, Hu D, Pino-Yanes M, Huntsman S, Eng C, Levin AM, et al. A genome-wide association and admixture mapping study of bronchodilator drug response in African Americans with asthma. Pharmacogenomics J. 2019;19(3):249–59.

    Article  CAS  PubMed  Google Scholar 

  7. Kabesch M, Tost J. Recent findings in the genetics and epigenetics of asthma and allergy. Semin Immunopathol. 2020;42(1):43–60.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Mak ACY, White MJ, Eckalbar WL, Szpiech ZA, Oh SS, Pino-Yanes M, et al. Whole-genome sequencing of pharmacogenetic drug response in racially diverse children with asthma. Am J Respir Crit Care Med. 2018;197(12):1552–64.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Park HW, Tantisira KG, Weiss ST. Pharmacogenomics in asthma therapy: Where are we and where do we go? Annu Rev Pharmacol Toxicol. 2015;55:129–47.

    Article  CAS  PubMed  Google Scholar 

  10. Sheikhpour M, Maleki M, Ebrahimi Vargoorani M, Amiri V. A review of epigenetic changes in asthma: methylation and acetylation. Clin Epigenetics. 2021;13(1):65.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Edris A, den Dekker HT, Melén E, Lahousse L. Epigenome-wide association studies in asthma: a systematic review. Clin Exp Allergy. 2019;49(7):953–68.

    Article  CAS  PubMed  Google Scholar 

  12. Poole A, Urbanek C, Eng C, Schageman J, Jacobson S, O’Connor BP, et al. Dissecting childhood asthma with nasal transcriptomics distinguishes subphenotypes of disease. J Allergy Clin Immunol. 2014;133(3):670-678.e12.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Herrera-Luis E, Li A, Mak ACY, Perez-Garcia J, Elhawary JR, Oh SS, et al. Epigenome-wide association study of lung function in Latino children and youth with asthma. Clin Epigenetics. 2022;14(1):9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Perez-Garcia J, Herrera-Luis E, Li A, Mak ACY. Multi-omic approach associates blood methylome with bronchodilator drug response in pediatric asthma. J Allergy Clin Immunol. 2023;151(6):1503–12.

  15. Cardenas A, Sordillo JE, Rifas-Shiman SL, Chung W, Liang L, Coull BA, et al. The nasal methylome as a biomarker of asthma and airway inflammation in children. Nat Commun. 2019;10(1):3095.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Kere M, Gruzieva O, Ullemar V, Söderhäll C, Greco D, Kull I, et al. Effects of inhaled corticosteroids on DNA methylation in peripheral blood cells in children with asthma. Allergy. 2020;75(3):688–91.

    Article  PubMed  Google Scholar 

  17. Wan ES, Qiu W, Baccarelli A, Carey VJ, Bacherman H, Rennard SI, et al. Systemic steroid exposure is associated with differential methylation in chronic obstructive pulmonary disease. Am J Respir Crit Care Med. 2012;186(12):1248–55.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Mansell G, Gorrie-Stone TJ, Bao Y, Kumari M, Schalkwyk LS, Mill J, et al. Guidance for DNA methylation studies: statistical insights from the Illumina EPIC array. BMC Genomics. 2019;20(1):366.

    Article  PubMed  PubMed Central  Google Scholar 

  19. National Heart L and BITO for PM (TOPMed) P. TOPMed whole genome sequencing methods: freeze 8 [Internet]. 2020.

  20. Wendell SG, Fan H, Zhang C. G protein-coupled receptors in asthma therapy: pharmacology and drug action. Pharmacol Rev. 2020;72(1):1–49.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Kamikawa Y, Saito A, Matsuhisa K, Kaneko M, Asada R, Horikoshi Y, et al. OASIS/CREB3L1 is a factor that responds to nuclear envelope stress. Cell Death Discov. 2021;7(1):152.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Fox RM, Hanlon CD, Andrew DJ. The CrebA/Creb3-like transcription factors are major and direct regulators of secretory capacity. J Cell Biol. 2010;191(3):479–92.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Vellanki RN, Zhang L, Volchuk A. OASIS/CREB3L1 is induced by endoplasmic reticulum stress in human glioma cell lines and contributes to the unfolded protein response, extracellular matrix production and cell migration. PLoS ONE. 2013;8(1): e54060.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Denard B, Seemann J, Chen Q, Gay A, Huang H, Chen Y, et al. The membrane-bound transcription factor CREB3L1 is activated in response to virus infection to inhibit proliferation of virus-infected cells. Cell Host Microbe. 2011;10(1):65–74.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Goldfarbmuren KC, Jackson ND, Sajuthi SP, Dyjack N, Li KS, Rios CL, et al. Dissecting the cellular specificity of smoking effects and reconstructing lineages in the human airway epithelium. Nat Commun. 2020;11(1):2485.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Zou Y, Li Q, Jiang L, Guo C, Li Y, Yu Y, et al. DNA Hypermethylation of CREB3L1 and Bcl-2 associated with the mitochondrial-mediated apoptosis via PI3K/Akt pathway in human BEAS-2B cells exposure to silica nanoparticles. PLoS ONE. 2016;11(6): e0158475.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Lutz SM, Cho MH, Young K, Hersh CP, Castaldi PJ, McDonald ML, et al. A genome-wide association study identifies risk loci for spirometric measures among smokers of European and African ancestry. BMC Genet. 2015;16:138.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Chiappara G, Chanez P, Bruno A, Pace E, Pompeo F, Bousquet J, et al. Variable p-CREB expression depicts different asthma phenotypes. Allergy. 2007;62(7):787–94.

    Article  CAS  PubMed  Google Scholar 

  29. Barnes PJ. Scientific rationale for inhaled combination therapy with long-acting β2-agonists and corticosteroids. Eur Respir J. 2002;19(1):182–91.

    Article  CAS  PubMed  Google Scholar 

  30. Yadav SK, Sharma P, Shah SD, Panettieri RA, Kambayashi T, Penn RB, et al. Autocrine regulation of airway smooth muscle contraction by diacylglycerol kinase. J Cell Physiol. 2022;237(1):603–16.

    Article  CAS  PubMed  Google Scholar 

  31. Peters MJ, Adcock IM, Brown CR, Barnes PJ. Beta-adrenoceptor agonists interfere with glucocorticoid receptor DNA binding in rat lung. Eur J Pharmacol. 1995;289(2):275–81.

    Article  CAS  PubMed  Google Scholar 

  32. Hancox RJ, Stevens DA, Adcock IM, Barnes PJ, Taylor DR. Effects of inhaled β agonist and corticosteroid treatment on nuclear transcription factors in bronchial mucosa in asthma. Thorax. 1999;54(6):488–92.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Greenwood M, Greenwood MP, Mecawi AS, Loh SY, Rodrigues JA, Paton JFR, et al. Transcription factor CREB3L1 mediates cAMP and glucocorticoid regulation of arginine vasopressin gene transcription in the rat hypothalamus. Mol Brain. 2015;8(1):68.

    Article  PubMed  PubMed Central  Google Scholar 

  34. Yang M, Zhang T, Zhang Y, Ma X, Han J, Zeng K, et al. MYLK4 promotes tumor progression through the activation of epidermal growth factor receptor signaling in osteosarcoma. J Exp Clin Cancer Res. 2021;40(1):166.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Gao L, Grant AV, Rafaels N, Stockton-Porter M, Watkins T, Gao P, et al. Polymorphisms in the myosin light chain kinase gene that confer risk of severe sepsis are associated with a lower risk of asthma. J Allergy Clin Immunol. 2007;119(5):1111–8.

    Article  CAS  PubMed  Google Scholar 

  36. Acosta-Herrera M, Pino-Yanes M, Ma SF, Barreto-Luis A, Corrales A, Cumplido J, et al. Fine mapping of the myosin light chain kinase (MYLK) gene replicates the association with asthma in populations of Spanish descent. J Allergy Clin Immunol. 2015;136(4):1116-1118.e9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Stull JT, Tansey MG, Tang DC, Word RA, Kamm KE. Phosphorylation of myosin light chain kinase: a cellular mechanism for Ca2+ desensitization. Mol Cell Biochem. 1993;127–128(1):229–37.

    Article  PubMed  Google Scholar 

  38. Parvathaneni S, Li Z, Sacks DB. Calmodulin influences MAPK signaling by binding KSR1. J Biol Chem. 2021;296: 100577.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Mao Z, Nakamura F. Structure and function of filamin C in the muscle Z-disc. Int J Mol Sci. 2020;21(8):2696.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Akdis M, Aab A, Altunbulakli C, Azkur K, Costa RA, Crameri R, et al. Interleukins (from IL-1 to IL-38), interferons, transforming growth factor β, and TNF-α: receptors, functions, and roles in diseases. J Allergy Clin Immunol. 2016;138(4):984–1010.

    Article  CAS  PubMed  Google Scholar 

  41. Park CS, Kim JW, Uh ST, Chung YT, Kim YH, Choi BW, Hue SH. Interleukin-2 and soluble interleukin-2 receptor in bronchoalveolar lavage fluid from patients with bronchial asthma. Chest. 1994;106(2):400–6.

    Article  CAS  PubMed  Google Scholar 

  42. Nag S, Lamkhioued B, Renzi PM. Interleukin-2–induced Increased airway responsiveness and lung Th2 cytokine expression occur after antigen challenge through the leukotriene pathway. Am J Respir Crit Care Med. 2012;165(11):1540–5.

    Article  Google Scholar 

  43. Suda T, Hashizume H, Aoshima Y, Yokomura K, Sato J, Inui N, et al. Management of interleukin-2-induced severe bronchoconstriction. Eur Respir J. 2007;29(3):612–3.

    Article  CAS  PubMed  Google Scholar 

  44. Ferrada MA, Gordon EL, Jen KY, He HZ, Lu X, Barone LM, et al. (R)-albuterol decreases immune responses: role of activated T cells. Respir Res. 2008;9(1):3.

    Article  PubMed  PubMed Central  Google Scholar 

  45. Kato A. Group 2 innate lymphoid cells in airway diseases. Chest. 2019;156(1):141–9.

    Article  PubMed  PubMed Central  Google Scholar 

  46. Farmer P, Pugin J. beta-adrenergic agonists exert their “anti-inflammatory” effects in monocytic cells through the IkappaB/NF-kappaB pathway. Am J Physiol Lung Cell Mol Physiol. 2000;279(4):L675-682.

    Article  CAS  PubMed  Google Scholar 

  47. Brown SD, Brown LA, Stephenson S, Dodds JC, Douglas SL, Qu H, et al. Characterization of a “high” TNF-α phenotype in moderate-to-severe asthmatic children. J Allergy Clin Immunol. 2015;135(6):1651–4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Berry MA, Hargadon B, Shelley M, Parker D, Shaw DE, Green RH, et al. Evidence of a role of tumor necrosis factor alpha in refractory asthma. N Engl J Med. 2006;354(7):697–708.

    Article  CAS  PubMed  Google Scholar 

  49. Pang L, Knox AJ. Synergistic inhibition by beta(2)-agonists and corticosteroids on tumor necrosis factor-alpha-induced interleukin-8 release from cultured human airway smooth-muscle cells. Am J Respir Cell Mol Biol. 2000;23(1):79–85.

    Article  CAS  PubMed  Google Scholar 

  50. Ezeamuzie CI, Shihab PK. Interactions between theophylline and salbutamol on cytokine release in human monocytes. J Pharmacol Exp Ther. 2010;334(1):302–9.

    Article  CAS  PubMed  Google Scholar 

  51. Toki S, Goleniewska K, Reiss S, Zhou W, Newcomb DC, Bloodworth MH, et al. The histone deacetylase inhibitor trichostatin A suppresses murine innate allergic inflammation by blocking group 2 innate lymphoid cell (ILC2) activation. Thorax. 2016;71(7):633–45.

    Article  PubMed  Google Scholar 

  52. Banerjee A, Trivedi CM, Damera G, Jiang M, Jester W, Hoshi T, et al. Trichostatin A abrogates airway constriction, but not inflammation, in murine and human asthma models. Am J Respir Cell Mol Biol. 2012;46(2):132–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Hernandez-Pacheco N, Vijverberg SJ, Herrera-Luis E, Li J, Sio YY, Granell R, et al. Genome-wide association study of asthma exacerbations despite inhaled corticosteroid use. Eur Respir J. 2021;57(5):2003388.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Perez-Garcia J, Espuela-Ortiz A, Hernández-Pérez JM, González-Pérez R, Poza-Guedes P, Martin-Gonzalez E, et al. Human genetics influences microbiome composition involved in asthma exacerbations despite inhaled corticosteroid treatment. J Allergy Clin Immunol. 2023;152(3):799–806.e6.

  55. Han YY, Forno E, Celedón JC. Sex steroid hormones and asthma in a nationwide study of U.S. adults. Am J Respir Crit Care Med. 2020;201(2):158–66.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Tsai PC, Bell JT. Power and sample size estimation for epigenome-wide association scans to detect differential DNA methylation. Int J Epidemiol. 2015;44(4):1429–41.

    Article  PubMed  PubMed Central  Google Scholar 

  57. Neophytou AM, White MJ, Oh SS, Thakur N, Galanter JM, Nishimura KK, et al. Air Pollution and lung function in minority youth with asthma in the GALA II (Genes-Environments and Admixture in Latino Americans) and SAGE II (Study of African Americans, Asthma, Genes, and Environments) studies. Am J Respir Crit Care Med. 2016;193(11):1271–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Reynolds SD, Rios C, Wesolowska-Andersen A, Zhuang Y, Pinter M, Happoldt C, et al. Airway progenitor clone formation is enhanced by Y-27632-dependent changes in the transcriptome. Am J Respir Cell Mol Biol. 2016;55(3):323–36.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Everman JL, Rios C, Seibold MA. Utilization of air–liquid interface cultures as an in vitro model to assess primary airway epithelial cell responses to the type 2 cytokine interleukin-13. Methods Mol Biol. 2018;1799:419–32.

    Article  CAS  PubMed  Google Scholar 

  60. Heiss JA, Just AC. Identifying mislabeled and contaminated DNA methylation microarray data: an extended quality control toolset with examples from GEO. Clin Epigenetics. 2018;10:73.

    Article  PubMed  PubMed Central  Google Scholar 

  61. R Core Team. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2020.

    Google Scholar 

  62. Xu Z, Niu L, Li L, Taylor JA. ENmix: a novel background correction method for Illumina HumanMethylation450 BeadChip. Nucleic Acids Res. 2016;44(3): e20.

    Article  PubMed  Google Scholar 

  63. Niu L, Xu Z, Taylor JA. RCP: a novel probe design bias correction method for Illumina Methylation BeadChip. Bioinformatics. 2016;32(17):2659–63.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Xu Z, Langie SAS, De Boever P, Taylor JA, Niu L. RELIC: a novel dye-bias correction method for Illumina Methylation BeadChip. BMC Genomics. 2017;18(1):4.

    Article  PubMed  PubMed Central  Google Scholar 

  65. Pidsley R, Zotenko E, Peters TJ, Lawrence MG, Risbridger GP, Molloy P, et al. Critical evaluation of the Illumina MethylationEPIC BeadChip microarray for whole-genome DNA methylation profiling. Genome Biol. 2016;17(1):208.

    Article  PubMed  PubMed Central  Google Scholar 

  66. McCartney DL, Walker RM, Morris SW, McIntosh AM, Porteous DJ, Evans KL. Identification of polymorphic and off-target probe binding sites on the Illumina Infinium MethylationEPIC BeadChip. Genom Data. 2016;9:22–4.

    Article  PubMed  PubMed Central  Google Scholar 

  67. Hansen KD. IlluminaHumanMethylationEPICanno.ilm10b4.hg19: Annotation for Illumina’s EPIC methylation arrays. R package version 0.6.0. 2016.

  68. Rahmani E, Zaitlen N, Baran Y, Eng C, Hu D, Galanter J, et al. Sparse PCA corrects for cell type heterogeneity in epigenome-wide association studies. Nat Methods. 2016;13(5):443–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The SVA package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. 2012;28(6):882–3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7): e47.

    Article  PubMed  PubMed Central  Google Scholar 

  71. van Iterson M, van Zwet EW, Heijmans BT, ’t Hoen PAC, van Meurs J, Jansen R, et al. Controlling bias and inflation in epigenome-and transcriptome-wide association studies using the empirical null distribution. Genome Biol. 2017;18(1):19.

    Article  PubMed  PubMed Central  Google Scholar 

  72. Campagna MP, Xavier A, Lechner-Scott J, Maltby V, Scott RJ, Butzkueven H, et al. Epigenome-wide association studies: current knowledge, strategies and recommendations. Clin Epigenetics. 2021;13(1):214.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. McLean CY, Bristor D, Hiller M, Clarke SL, Schaar BT, Lowe CB, et al. GREAT improves functional interpretation of cis-regulatory regions. Nat Biotechnol. 2010;28(5):495–501.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. National Asthma Education. Expert Panel Report 3 (EPR-3): guidelines for the diagnosis and management of asthma-summary report 2007. J Allergy Clin Immunol. 2007;120(5 Suppl):S94-138.

    Google Scholar 

  75. Han B, Eskin E. Random-effects model aimed at discovering associations in meta-analysis of genome-wide association studies. Am J Hum Genet. 2011;88(5):586–98.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Ongen H, Buil A, Brown AA, Dermitzakis ET, Delaneau O. Fast and efficient QTL mapper for thousands of molecular phenotypes. Bioinformatics. 2016;32(10):1479–85.

    Article  CAS  PubMed  Google Scholar 

  77. Bushnell B. BBMap: A Fast, Accurate, Splice-Aware Aligner. No. LBNL-7065E. Berkeley: Ernest Orlando Lawrence Berkeley National Laboratory; 2014.

    Google Scholar 

  78. Patro R, Duggal G, Love MI, Irizarry RA, Kingsford C. Salmon provides fast and bias-aware quantification of transcript expression. Nat Methods. 2017;14(4):417–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  79. Soneson C, Love MI, Robinson MD. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Research. 2016;4:1521.

    Article  PubMed Central  Google Scholar 

  80. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

    Article  PubMed  PubMed Central  Google Scholar 

  81. Shabalin AA. Matrix eQTL: ultra fast eQTL analysis via large matrix operations. Bioinformatics. 2012;28(10):1353–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  82. Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 2016;44(W1):W90–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  83. Pedersen BS, Schwartz DA, Yang IV, Kechris KJ. Comb-p: Software for combining, analyzing, grouping and correcting spatially correlated P-values. Bioinformatics. 2012;28(22):2986–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Peters TJ, Buckley MJ, Statham AL, Pidsley R, Samaras K, V Lord R, et al. De novo identification of differentially methylated regions in the human genome. Epigenetics Chromatin. 2015;8:6.

    Article  PubMed  PubMed Central  Google Scholar 

Download references


The authors would like to thank the patients, families, recruiters, healthcare providers, and community clinics for their participation in the study. We thank Sandra Salazar for her support as GALA II study coordinator. We thank the investigators and teams at the New York Genome Center and the Northwest Genomics Center for whole-genome sequencing sample preparation, quality control, data generation, and data processing. We gratefully acknowledge the studies and participants who provided biological samples and data for TOPMed. We also gratefully acknowledge the contributions of the investigators of the NHLBI TOPMed Consortium ( The authors thank Soren Germer, Michael C. Zody, Lara Winterkorn, and Catherine Reeves from NYGC, and Deborah A. Nickerson from NWGC for overseeing the production of whole genome sequencing data for GALA II. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. The data used for the in silico eQTL analyses described in this manuscript were obtained from the GTEx Portal on 16/06/22. The Genotype-Tissue Expression (GTEx) Project was supported by the Common Fund of the Office of the Director of the National Institutes of Health, and by NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS. JP-G would like to thank the International Mentoring Foundation for the Advancement of Higher Education (IMFAHE), the University of La Laguna, and the Spanish Ministry of Universities for financially supporting his international internship at the University of California, San Francisco.


Whole-genome sequencing (WGS) for the TOPMed (Trans-Omics in Precision Medicine) program was supported by the NHLBI. Whole-genome sequencing for "NHLBI TOPMed: Gene-Environment, Admixture and Latino Asthmatics Study" (phs000920, GALA II) was performed at the New York Genome Center (3R01HL117004-02S3) and Northwest Genomics Center (NWGC, HHSN268201600032I). Centralized read mapping and genotype calling, along with variant quality metrics and filtering, were provided by the TOPMed Informatics Research Center (3R01HL-117626-02S1). Phenotype harmonization, data management, sample-identity QC, and general study coordination were provided by the TOPMed Data Coordinating Center (3R01HL-120393-02S1, contract HHSN268201800002I). WGS of part of GALA II was performed by the New York Genome Center under The Centers for Common Disease Genomics of the Genome Sequencing Program (GSP) Grant (UM1 HG008901). The GSP Coordinating Center (U24 HG008956) contributed to cross-program scientific initiatives and provided logistical and general study coordination. GSP is funded by the National Human Genome Research Institute, the National Heart, Lung, and Blood Institute, and the National Eye Institute. This work was supported in part by the Sandler Family Foundation, the American Asthma Foundation, the Amos Medical Faculty Development Program from the Robert Wood Johnson Foundation, the Harry Wm. and Diana V. Hind Distinguished Professor in Pharmaceutical Sciences II, the National Heart, Lung, and Blood Institute of the National Institutes of Health (R01HL155024-01, R01HL117004, R01HL128439, R01HL135156, X01HL134589), the National Institute of Health and Environmental Health Sciences (R01ES015794, R21ES24844), the National Institute on Minority Health and Health Disparities (R01MD010443, R56MD013312), the Tobacco-Related Disease Research Program under Award Number (24RT-0025 and 27IR-0030). This work was also funded by the Spanish Ministry of Science and Innovation grant PID2020-116274RB-I00 awarded by MCIN/AEI/10.13039/501100011033. MP-Y was supported by the Ramón y Cajal Program by the MCIN/AEI/10.13039/501100011033 and the European Social Fund “Investing in your future” (RYC-2015-17205). JP-G was funded by fellowship FPU19/02175 from the Spanish Ministry of Universities. JV received funding from the Instituto de Salud Carlos III (ISCIII), Madrid, Spain (PI19/00141). MP-Y and JV received funds from the CIBER -Consorcio Centro de Investigación Biomédica en Red- (CIBERES), ISCIII, and European Regional Development Fund (CB06/06/1088).

Author information

Authors and Affiliations



JP-G, MP-Y, EGP, JLE, EZ, MAS, and EGB were involved in the conceptualization of the study; MAS, and EGB in funding acquisitions; JLE, CE, MAS, and EGB in project administration; KBB, VM, SS, DEW, FH, JRS, MAS, and EGB in resources; JP-G, MP-Y, EGP, JLE, CMM, EZ, MAS, and EGB in methodology; JP-G, MP-Y, EGP, JLE, CE, EZ, MAS, and EGB in investigation; JP-G, MP-Y, EGP, and JLE in data curation; JP-G, MP-Y, CMM, EZ, MAS, and EGB in formal analyses; MP-Y, JLE, JV, EZ, MAS, and EGB in supervision; JP-G and JLE in visualization; JP-G and MP-Y in writing—original draft; and EGP, JLE, KBB, JV, EZ, MAS, and EGB in writing—review and editing. All authors have read the manuscript and accepted the final version considered for publication.

Corresponding authors

Correspondence to Javier Perez-Garcia or Maria Pino-Yanes.

Ethics declarations

Ethics approval and consent to participate

The GALA II study was approved by the Human Research Protection Program (HRPP) Institutional Review Board (IRB) at the University of California, San Francisco (UCSF) (UCSF, IRB number 10–00889, Reference number 153543, NJH HS-2627). All parents and participants provided signed written consent and assent as appropriate. This work has been conducted in accordance with The Code of Ethics of the World Medical Association (Declaration of Helsinki).

Consent for publication

Not applicable.

Competing interests

Authors declare they have no competing interests or other interests that might be perceived to influence the interpretation of the article. No supporting institution may gain or lose financially through this publication. The funders of the study had no role in study design, data collection, data analysis, data interpretation, or writing of the report.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1

. Supplementary material (Supplementary Figures S1-S4 and Supplementary Tables S1-S14).

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Perez-Garcia, J., Pino-Yanes, M., Plender, E.G. et al. Epigenomic response to albuterol treatment in asthma-relevant airway epithelial cells. Clin Epigenet 15, 156 (2023).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: