Abstract
In 2019, the Global Initiative for Asthma treatment guidelines were updated to recommend that inhaled corticosteroid (ICS)/long-acting β2-adrenoceptor agonist (LABA) combination therapy should be a first-in-line treatment option for asthma. Although clinically superior to ICS, mechanisms underlying the efficacy of this combination therapy remain unclear. We hypothesized the existence of transcriptomic interactions, an effect that was tested in BEAS-2B and primary human bronchial epithelial cells (pHBECs) using formoterol and budesonide as representative LABA and ICS, respectively. In BEAS-2B cells, formoterol produced 267 (212 induced; 55 repressed) gene expression changes (≥2/≤0.5-fold) that were dominated by rapidly (1 to 2 hours) upregulated transcripts. Conversely, budesonide induced 370 and repressed 413 mRNAs, which occurred predominantly at 6–18 hours and was preceded by transcripts enriched in transcriptional regulators. Significantly, genes regulated by both formoterol and budesonide were over-represented in the genome; moreover, budesonide plus formoterol induced and repressed 609 and 577 mRNAs, respectively, of which ∼one-third failed the cutoff criterion for either treatment alone. Although induction of many mRNAs by budesonide plus formoterol was supra-additive, the dominant (and potentially beneficial) effect of budesonide on formoterol-induced transcripts, including those encoding many proinflammatory proteins, was repression. Gene ontology analysis of the budesonide-modulated transcriptome returned enriched terms for transcription, apoptosis, proliferation, differentiation, development, and migration. This “functional” ICS signature was augmented in the presence of formoterol. Thus, LABAs modulate glucocorticoid action, and comparable transcriptome-wide interactions in pHBECs imply that such effects may be extrapolated to individuals with asthma taking combination therapy. Although repression of formoterol-induced proinflammatory mRNAs should be beneficial, the pathophysiological consequences of other interactions require investigation.
SIGNIFICANCE STATEMENT In human bronchial epithelial cells, formoterol, a long-acting β2-adrenoceptor agonist (LABA), enhanced the expression of inflammatory genes, and many of these changes were reduced by the glucocorticoid budesonide. Conversely, the ability of formoterol to enhance both gene induction and repression by budesonide provides mechanistic insight as to how adding a LABA to an inhaled corticosteroid may improve clinical outcomes in asthma.
Introduction
As rescue medication, inhaled short-acting β2-adrenoceptor agonists provide acute symptom relief in individuals with asthma. Conversely, inhaled corticosteroids (ICS), acting via the glucocorticoid receptor (GR, NR3C1), reduce airway inflammation and generally control mild-to-moderate disease (Barnes, 2011; Oakley and Cidlowski, 2013). However, with more severe asthma, ICSs often provide inadequate control, and add-on therapy, typically a long-acting β2-adrenoceptor agonist (LABA), is recommended (Ross et al., 2015). Despite the ability of β2-adrenoceptor agonists to arrest and reverse bronchoconstriction, chronic use, especially of high-efficacy ligands, increases asthma mortality (Nelson et al., 2006; Cockcroft and Sears, 2013; Reddel et al., 2019). Explanations for this include the possibility that excessive β2-agonist use masks underlying disease progression and severity. Alternatively, β2-agonists may produce detrimental effects—for example, by enhancing expression of inflammatory genes. Certainly, β2-agonists increase the production of inflammatory mediators, particularly in the context of inflammatory stimuli (Ammit et al., 2002; Edwards et al., 2007; Strandberg et al., 2007; Holden et al., 2010). Moreover, transcriptomic analyses in human bronchial epithelial cells (HBECs) exposed to the LABAs indacaterol and salmeterol confirm increased expression of inflammatory genes (Yan et al., 2018).
Elevated mortality associated with chronic β2-agonist monotherapy is not seen in patients taking ICS/LABA combination therapy (Reddel et al., 2019). Accordingly, in 2019, the Global Initiative for Asthma (GINA) modified its treatment guidelines to recommend that ICS/LABA combination therapy, in which the LABA is formoterol, should be a first-line treatment option, even for patients with mild disease (https://ginasthma.org/reports/). Indeed, ICS/LABA combination therapy is superior to ICS alone, irrespective of dose, at improving lung function and reducing exacerbation rates (Newton and Giembycz, 2016). Such data raise the proposition that LABAs and ICSs may interact, possibly at a molecular level, to deliver superior clinical outcomes. Certainly, LABAs can cooperatively enhance glucocorticoid-driven transcription. For example, in BEAS-2B human airway epithelial cells, LABAs did not activate a simple glucocorticoid response element–dependent luciferase reporter, yet they enhanced (2- to 3-fold) the maximum glucocorticoid-induced response (Kaur et al., 2008; Joshi et al., 2015). Similar interactions occur for multiple genes (Kaur et al., 2008; Joshi et al., 2015; Rider et al., 2018). These include the regulator of G-protein signaling, RGS2, which attenuates signal transduction from G-protein–coupled receptors that act via Gq (Heximer, 2004; Kimple et al., 2011). In vivo, including in models of lung inflammation (Xie et al., 2012; Jiang et al., 2015; George et al., 2017; George et al., 2018), RGS2 is bronchoprotective (Holden et al., 2011). Moreover, in HBECs and smooth muscle cells, LABAs synergize with glucocorticoids to enhance and prolong RGS2 expression (Holden et al., 2011; Holden et al., 2014). This suggests therapeutic relevance.
Given that GR-dependent transactivation mediates anti-inflammatory effects of glucocorticoids (Newton and Holden, 2007; Clark and Belvisi, 2012; Newton et al., 2017; Oh et al., 2017), transcriptional cooperativity with LABA could be important to realize the enhanced therapeutic efficacy of ICS/LABA combination therapy. Equally, glucocorticoid-dependent repression of LABA-induced proinflammatory mediators would be beneficial (Ammit et al., 2002; Holden et al., 2010). However, as not all LABA- and glucocorticoid-induced genes show cooperative effects (Rider et al., 2018), detailed transcriptomic characterization is necessary in cells relevant to asthma pathogenesis.
Airway epithelial cells promote and regulate inflammation (Knight and Holgate, 2003), are a direct target of inhaled therapy, and are central to the anti-inflammatory effects of glucocorticoids in the airways (Klaßen et al., 2017). Furthermore, the ICS budesonide induces expression of multiple genes in vivo in the airways of healthy individuals and individuals with asthma (Kelly et al., 2012; Leigh et al., 2016). Many such genes are also induced in primary HBECs (pHBECs) and the HBEC line BEAS-2B (Mostafa et al., 2019). Because pHBECs and BEAS-2B cells express functional β2-adrenoceptors (Kaur et al., 2008; Rider et al., 2018; Yan et al., 2018), they were selected for transcriptomic analysis using formoterol and budesonide. These represent clinically relevant LABA and ICS, which are identified by GINA as the “preferred” or reference compounds, respectively (https://ginasthma.org/reports/). In each case, responses to formoterol and budesonide in epithelial cells were representative of other members of their class in modulating reporter activity and gene expression (Kaur et al., 2008; Rider et al., 2011, 2015; Yan et al., 2018). Following the hypothesis-free determination of differentially regulated genes by each treatment condition, gene ontology (GO) analyses were used to provide insight as to therapeutic relevance.
Materials and Methods
Cell Culture, Compounds, and Cell Treatments.
The human bronchial epithelial cell line BEAS-2B (American Type Culture Collection (ATCC), Manassas, VA), was cultured to confluence in Dulbecco’s modified Eagle’s medium/Ham’s F-12 medium supplemented with 14 mM NaHCO3, 2 mM l-glutamine, and 10% fetal calf serum (all Invitrogen; Burlington, ON). pHBECs from normal nonsmokers were cultured from airway brushings obtained at bronchoscopy (Leigh et al., 2016) or were isolated from nontransplantable normal human lungs obtained through a tissue retrieval service at the International Institute for the Advancement of Medicine (Edison, NJ) (Hudy et al., 2010). In each case, protocols and consent for the bronchoscopy study were approved by the Conjoint Health Research Ethics Board at the University of Calgary and Alberta Health Services (Research Ethics Board IDs: 23241 and 15-0336, respectively). The pHBECs were grown as submersion culture in bronchial epithelial cell growth medium (Lonza, Morristown, NJ) containing SingleQuots supplements (Lonza). BEAS-2B and pHBECs were cultured at 37°C in 5% CO2/95% air and, once confluent, were incubated in serum-free, or supplement-free, medium overnight prior to experiments. All cells tested negative on routine mycoplasma contamination testing. Budesonide (22RS 16α, 17α-butylidenedioxypregna-1,4-diene-11β,21-diol-3,20-dione) and the active RR version of formoterol fumarate [(R′, R′)-(+/−)-N-(2-hydroxy-5-(1-hydroxy-2-((2-(4-methoxyphenyl)-1-methylethyl)amino)ethyl)phenyl) formamide, (E)-2-butendioate (2), dihydrate] (gifts from AstraZeneca; Mölndal, Sweden) were dissolved in DMSO (Sigma-Aldrich) as stocks at 10 mM. Final DMSO concentrations on cells were ≤0.1%. Based on prior data that established maximally effective concentrations (Kaur et al., 2008; Rider et al., 2011; Holden et al., 2014; BinMahfouz et al., 2015; Mostafa et al., 2019), the BEAS-2B cells were either not treated or treated with formoterol (10 nM), budesonide (100 nM), or both combined for 1, 2, 6, and 18 hours. Likewise, pHBECs were treated with formoterol (10 nM), budesonide (100 nM), or both combined for 6 hours. In each case, treatments were conducted in serum-free or supplement-free medium to reduce the impact of signaling and transcriptional response due to mediators present in serum or other supplements. As was previously reported (Rider et al., 2011), this helps to ensure optimal responses to glucocorticoid and LABA.
Microarray Analysis.
Total RNA was extracted (RNeasy mini kit; Qiagen), and RNA quality was assessed on a 2100 Bioanalyzer using RNA 6000 Nano LabChips (Agilent Technologies, Santa Clara, CA). First and second strand synthesis was performed with GeneChip 3ʹ IVT Express kits (Affymetrix, Santa Clara, CA) and in vitro transcription-generated biotin-labeled amplified RNA. After purification and fragmentation, hybridization to PrimeView microarrays (Affymetrix) was for 16 hours prior to washing in a GeneChip Fluidics Station 450 and scanning with a Scanner 3000 G7. Robust multiarray averaging, quantile normalization, and median polishing on logged probe set intensity values were performed in Transcriptome Analysis Console software version 4.0 (Affymetrix). Transcriptome Analysis Console software was then used to identify differentially expressed genes between paired samples (treated vs. untreated) using the eBayes ANOVA method, as is appropriate for microarray data (Smyth, 2004), to produce descriptive statistics for data categorization. Data files were deposited with the The National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) under accession numbers GSE115830 and GSE161805 for BEAS-2B and pHBEC data, respectively.
When genes were represented by multiple probe sets, only those with the greatest overall change were retained. Four independent experiments in BEAS-2Bs and pHBECs from six healthy individuals were analyzed for gene expression changes due to formoterol, budesonide, or budesonide plus formoterol. In each case, data were expressed as fold of untreated for each time point, and genes showing ≥2- or ≤0.5-fold (each P ≤ 0.05) relative to untreated for any treatment were considered differentially expressed. Unless otherwise stated, these criteria are henceforth used to define induced or repressed gene expression. Supra- and infra-additive or enhanced repressive effects of the formoterol-plus-budesonide combination, compared with either treatment alone, were estimated using one-way ANOVA followed by Tukey’s post hoc test. Significant combinatorial effect is determined for genes showing adjusted P value of 0.05 or less. Positive or negative values of the difference between means define whether the significant combination effects were supra- or infra-additive, respectively. In each case, the calculated P values for differentially regulated genes compared with the untreated group or for the identification of combinatorial effects are purely descriptive. These are presented for the purpose of ranking and grouping purposes and should not be construed as testing a specific hypothesis.
Functional Annotation Analysis.
GO analyses were performed using the functional annotation chart and functional annotation clustering tools within the Database for Annotation, Visualization, and Integrated Discovery (DAVID) version 6.8 (Huang et al., 2009). GO categories for analysis were restricted to molecular function, biologic process, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway terms using the DAVID default cutoff, P value/EASE score (PEASE) ≤ 0.1, for enrichment. Multiple testing correction of enrichment P values, Benjamini (PB) ≤ 0.05 (or, where stated, ≤0.01), were used to highlight the more robustly enriched terms. For functional annotation analyses, the criterion of three or more genes per term was also required. In addition, and where specifically stated, less conservative search criteria using the DAVID default category settings were conducted to increase GO term coverage. For functional annotation clustering, enrichment scores representing −log10 of the geometric mean of all enrichment P values of each term in cluster were, as recommended (Huang et al., 2009), taken as significant when ≥1.3 [i.e., −log10(0.05)]. For clarity, GO terms are referred to using quotation marks.
Graphical Presentation and χ2 Testing.
Graphs were produced using GraphPad Prism version 6.01 software (GraphPad Software Inc., La Jolla, CA) or the R package ggplot2. The R package pheatmap was used to produce heat maps. χ2 testing was performed using Prism version 6.01 and assumes 18,843 unique mRNAs represented on the PrimeView arrays.
Results
Glucocorticoid- and LABA-Regulated Transcriptomes in BEAS-2B Cells.
Over an 18-hour time frame in BEAS-2B cells, formoterol (10 nM), budesonide (100 nM), or budesonide plus formoterol modified the expression of 267, 777, and 1165 genes, respectively (Fig. 1A; Supplemental Table 1). With formoterol, this was primarily represented by transient gene induction, with maximal effects occurring at 2 hours, and relatively few mRNAs were repressed (Fig. 1, A and B). By contrast, budesonide produced a modest 1- to 2-hour response that was primarily due to gene induction and that progressively increased to record the highest gene count at 18 hours (Fig. 1, A and B). Budesonide also promoted significant but delayed gene repression, such that at 18 hours the number of repressed transcripts greatly exceeded those that were induced (Fig. 1B). Similar profiles of gene expression changes were produced by budesonide plus formoterol, except that the numbers of genes and genes with peak expression at each time were generally increased relative to either treatment alone (Fig. 1, A and B). This was most apparent with the genes showing peak induction at 6 hours and those showing peak repression at 6 and 18 hours (Fig. 1B).
Overview effects of formoterol and/or budesonide on gene expression in BEAS-2B cells. BEAS-2B cells were either not stimulated (NS) or treated with formoterol (Form; 10 nM), budesonide (Bud; 100 nM), or both together (Bud + Form, B+F) for the indicated times. Cells (N = 4 separate experiments) were harvested for microarray analysis. Normalized data were expressed as fold relative to time-matched NS, and genes showing ≥2-fold or ≤0.5-fold (P ≤ 0.05) were considered to be differentially expressed (Supplemental Table 1). (A) The number of differentially expressed genes showing up- or downregulation at each time is shown for each treatment. The overall number of up- or downregulated genes at any time with each treatment is shown (right-hand panel). (B) For each treatment, the number of genes showing peak expression change from NS (+ or −log2 fold) were grouped according to the time of peak expression and the direction of change (up/down). Expression data for each group were summarized using local polynomial regression (solid line) ±95% confidence interval (CI; gray shading) of log2 fold over time (hours). For groups containing three or fewer genes, the regression line of each gene (dashed line) was plotted instead of CI. (C) The number of genes meeting the expression criteria for (i) induction (≥2-fold, P ≤ 0.05) or (ii) repression (≤0.5-fold, P ≤ 0.05) at any time with respect to each treatment is represented as Venn diagrams. Below each Venn, contingency tables assess the likelihood (χ2 test) that the number of genes (induced or repressed) in common between formoterol and budesonide may have occurred by chance. ****P = 0.0001. E, expected; O, observed.
Coregulation by Glucocorticoid and LABA.
With ∼20,000 genes in the human genome (Clamp et al., 2007; Ezkurdia et al., 2014) and 18,843 unique genes represented on PrimeView arrays, independent assortment for up- or downregulation by formoterol and budesonide should produce only one or two mRNAs regulated by both stimuli at each time. However, in common up- and downregulation by budesonide and formoterol was markedly over-represented (Fig. 1C; Supplemental Fig. 1A).
Scatterplots of the formoterol- and budesonide-regulated mRNAs also revealed over-representation of formoterol-induced genes that independently showed repression by budesonide (Supplemental Fig. 1A). Indeed, 19.3% of formoterol-induced mRNAs were independently repressed by budesonide (Supplemental Fig. 1B), an effect that was significant and consistent with the 39.2% of formoterol-induced mRNAs that failed the criteria for induction by formoterol plus budesonide (Fig. 1Ci).
In addition to coregulation, budesonide plus formoterol induced or repressed 196 and 198 mRNAs, respectively, that failed to meet the criteria for regulation by either treatment alone (Fig. 1C). These mRNAs represented ∼one-third of the total budesonide plus formoterol induced or repressed transcriptome and suggest extensive gene-dependent interactions between each treatment.
LABA-Induced Genes and the Effect of Glucocorticoid.
Although formoterol-induced genes showed variability in their mRNA kinetics (Fig. 2A), positive and negative effects of budesonide (defined as ≥2/≤0.5-fold further change at any time) were produced at each time of peak formoterol-induced expression (Fig. 2B). Overall, 58.5% of formoterol-induced genes were modulated by budesonide, and repression (36.8%) was the dominant effect (Fig. 2C). This was also consistent with 39.2% of formoterol-induced genes not achieving inducibility by budesonide plus formoterol (Fig. 1Ci). Conversely, many of the 46 formoterol-induced genes that were enhanced by budesonide revealed evidence of independent inducibility by budesonide (Fig. 2B). Indeed 37 of these fell within those genes showing independent formoterol and budesonide inducibility (Fig. 1Ci).
Formoterol- and budesonide-induced gene expression and effects of combination. (A) The 212 genes induced (≥2-fold, P ≤ 0.05) at any time by formoterol in BEAS-2B cells were grouped according to the time of peak expression (log2 fold) and are presented as heat maps. Red = induced, blue = repressed, white = log2(1), and color intensity reflects increasing change in log2 fold. (B) Expression of the genes showing peak formoterol-induced expression at each time was summarized using local polynomial regression (solid line) ± 95% confidence interval (CI; gray shading) of log2 fold over time (hours). For each time of peak formoterol-induced expression, the genes were stratified according to the greatest effect (+1, 0, or −1) of budesonide occurring at any time: +1 = ≥2-fold, 0 = >0.5/<2-fold, and −1 = ≤0.5-fold produced by budesonide addition when compared with formoterol-induced expression. Red line = formoterol, blue line = budesonide, and purple line = budesonide + formoterol. (C) Graphical representation showing the gene count for formoterol-induced genes meeting the criteria (+1, 0, or −1) in (B) for the effect of budesonide (Bud) addition. (D) The 370 genes induced (≥2-fold, P ≤ 0.05) at any time by budesonide in BEAS-2B cells were grouped according to the time of peak expression (log2 fold) and presented as heat maps (red = induced, blue = repressed). (E) Expression of genes showing peak budesonide-induced expression at each time was summarized as in (B). For each time of peak budesonide-induced expression group, the genes were stratified according to the greatest effect (+1, 0, or −1) of formoterol occurring at any time: +1 = ≥2-fold, 0 = >0.5/<2-fold, and −1 = ≤0.5-fold produced by formoterol addition when compared with budesonide-induced expression. Line colors are as in (B). (F) Graphical representation showing the gene count for budesonide-induced genes meeting the criteria (+1, 0, or −1) in (E) for the effect of formoterol (Form) addition. (G) The 196 genes induced (≥2-fold, P ≤ 0.05) at any time by budesonide + formoterol, but not by either formoterol or budesonide alone in BEAS-2B cells, were grouped according to the time of peak expression and summarized as in (B). Line colors are as in (B).
Glucocorticoid-Induced Genes and the Effect of LABA.
The prevalence of delayed induction and considerable variability in kinetics was apparent among the 370 budesonide-induced genes (Fig. 2D). Further modulation (≥2/≤0.5-fold further change at any time) by formoterol occurred for 20.5% of these transcripts, with 3.0% and 17.5% being repressed or enhanced, respectively (Fig. 2, E and F). As described above, independent induction by each stimulus was clearly apparent, and 29 of the formoterol-enhanced mRNAs showed independent inducibility by formoterol (Fig. 2, E and F). Furthermore, an additional 20 budesonide-induced mRNAs showed significant, but <2-fold, induction by formoterol. Thus, in BEAS-2B cells, 75% of the budesonide-induced mRNAs that were enhanced by formoterol were also, at least modestly, induced by formoterol alone.
Genes Only Upregulated by LABA plus Glucocorticoid.
Of the 196 mRNAs uniquely induced by budesonide plus formoterol, 85 showed peak expression at 6 hours, and overall, 70.9% of these genes showed peak induction at 6 or 18 hours (Fig. 2G). Furthermore, 133 or 168 of these genes were modestly (fold >1 but <2), yet significantly (P ≤ 0.05), induced by formoterol or budesonide alone (at any time), and 115 genes were significantly upregulated by both treatments (Supplemental Table 1). Thus, low-level inducibility may drive combinatorial effects.
LABA/Glucocorticoid Interactions: Supra- and Infra-additive Effects.
To explore combinatorial interactions, the effect (i.e., fold - 1) due to each of formoterol and budesonide were summated and then compared with the effect of budesonide plus formoterol. Applying this to the 46 genes upregulated in common by budesonide and formoterol revealed 20 genes for which the effect of budesonide plus formoterol was significantly (P ≤ 0.05) greater than the sum of the effects of budesonide and formoterol (i.e., supra-additivity) at one, or more, times (Fig. 3A; Supplemental Table 2). However, multiple other genes were induced in a simple additive manner, and 13 of these coregulated genes revealed significantly less than additivity (infra-additivity) at one or more times (Fig. 3A). For example, this was clearly apparent with SGK1 and KLF9, which were induced in common and showed peak expression at 1 or 2 hours, respectively. Thus, despite coregulation, supra-additivity was not inevitable with cotreatment; absence of interaction and negative interactions also occurred. A similar analysis for the 196 mRNAs only induced by budesonide plus formoterol revealed significant supra-additivity between budesonide and formoterol for 107 genes (Fig. 3B; Supplemental Table 2). Examples illustrating the dramatic nature of some such effects are shown (Fig. 3C).
Supra- and infra-additive effects of formoterol and budesonide on gene expression. (A) The 46 genes induced (≥2-fold, P ≤ 0.05) by both formoterol and budesonide at any time in BEAS-2B cells were grouped according to the time of peak expression, and expression was summarized using local polynomial regression (solid line) ± 95% confidence interval (CI; gray shading) of log2 fold over time (hours) (left panel). For each group of genes, heat maps were generated to show fold gene expression for each treatment and time (center panel). Red = induced, blue = repressed, white = 1, and color intensity reflects fold change. For each gene (right-hand main panel), the effect (fold-1) of budesonide and formoterol was summed (B − 1 + F − 1 = B + F − 2) and displayed next to the effect of budesonide + formoterol (BF − 1) as a heat map using the same scale as the prior panel. Differences between B + F − 2 and BF − 1 were compared by one-way ANOVA followed by Tukey’s post hoc test (far right panel). P values presented as heat maps are shown in which purple color intensity corresponds to increasing −log10(P value), whereas values below 1.3 (P < 0.05) are all white. (B) The 196 genes showing induced expression by budesonide + formoterol (2 ≥ fold, P ≤ 0.05), but not by budesonide or formoterol alone, were grouped, as in Fig. 2G, by time of peak expression. The effect of budesonide was summed with the effect of formoterol (B + F − 2) and displayed as a heat map next to the effect of budesonide + formoterol (BF − 1) (left panel) using the same scale as in (A). Similarly, significance was tested, and P values are presented as a heat map (right panel), as in (A). (C) For each time of budesonide + formoterol–induced peak expression, one gene showing budesonide+ formoterol, but not budesonide or formoterol alone, induced expression from (B) is shown. (D) Examples of formoterol-induced genes displaying significant supra- or infra-additivity (upper panels and lower panels, respectively) are shown for peak formoterol-induced expression at 1, 2, and 6 hours. Color coding is as in (C). (E) Examples of formoterol-induced genes displaying significant supra- or infra-additivity (upper panels and lower panels, respectively) are shown for peak formoterol-induced expression at 1, 2, and 6 hours. Color coding is as in (C). B, budesonide; BF, budesonide plus formoterol; F, formoterol.
Extending the above analysis to all 732 genes induced by any treatment revealed 259 mRNAs showing supra-additivity and 196 genes with infra-additive effects (all ≤0.05) (Supplemental Table 2). Of the 212 formoterol-induced mRNAs, 45 were increased in a supra-additive manner, whereas 111 showed infra-additivity when formoterol was combined with budesonide (Supplemental Fig. 2A; Supplemental Table 2). Figure 3D reveals striking positive enhancements (upper panels), as well as profound repression of IL6, IL11, or PDE7B (lower panels), by budesonide on formoterol-induced gene expression. For the 370 budesonide-induced mRNAs, 127 showed supra-additivity and 87 revealed infra-additivity (all P ≤ 0.05) with formoterol cotreatment (Supplemental Fig. 2B; Supplemental Table 2). As previously reported, this illustrates the extensive enhancement of budesonide-induced genes by formoterol (Altonsy et al., 2017; Rider et al., 2018), an effect which occurred at each time of peak budesonide-induced gene expression (Fig. 3E, upper panels; Supplemental Fig. 2B). Infra-additivity was also apparent at early, intermediate, and late times of peak budesonide-induced expression (Fig. 3E, lower panels; Supplemental Fig. 2B).
Finally, many mRNAs showed both positive and negative interactions, albeit at different times. This is illustrated by AREG (Fig. 3E), a budesonide-induced gene that was enhanced by budesonide plus formoterol at 1 and 2 hours, but at 6 hours, formoterol significantly reduced budesonide-induced expression. Similar effects were apparent for multiple other budesonide-, formoterol-, and formoterol plus budesonide–induced genes (Supplemental Fig. 2C). Although reflecting gene-specific control, this indicates the need for considerable caution when interpreting data derived from single time points.
LABA-Repressed Gene Expression and the Effect of Glucocorticoid.
Reduced expression of the 55 formoterol-repressed mRNAs was primarily restricted to the time of peak repression, and with 60% of transcripts not showing further modulation, the overall effect of budesonide was modest (Fig. 4, A–C). Thus, budesonide cotreatment further reduced expression of 15 formoterol-repressed genes and enhanced expression of seven genes (Fig. 4C).
Formoterol- and budesonide-repressed gene expression and effects of combination. (A) The 55 genes that were repressed (≤0.5-fold, P ≤ 0.05) at any time by formoterol in BEAS-2B cells were grouped according to the time of peak repression (log2 fold) and presented as heat maps [red = induced, blue = repressed, white = log2(1), and color intensity reflects increasing log2 fold change]. (B) Expression of the 2, 25, 24, and 4 genes showing peak formoterol-induced repression at 1, 2, 6, or 18 hours was summarized using local polynomial regression (solid line) ± 95% confidence interval (CI; gray shading) of log2 fold over time (hours). For each time of peak formoterol-induced repression group, the genes were further stratified according to the greatest effect (+1, 0, or −1) of budesonide occurring at any time: +1 = ≥2-fold, 0 = >0.5/<2-fold, and −1 = ≤0.5-fold produced by budesonide addition when compared with formoterol-repressed expression. Red line = formoterol, blue line = budesonide, and purple line = budesonide + formoterol. (C) Graphical representation showing the gene count for formoterol-repressed genes meeting the criteria (+1, 0, or −1) in (B) for the effect of budesonide (Bud) addition. (D) The 413 genes that were repressed (≤0.5-fold, P ≤ 0.05) at any time by budesonide in BEAS-2B cells were grouped according to the time of peak repression (log2 fold) and presented as heat maps (red = induced, blue = repressed). (E) Expression of the 3, 17, 132, and 261 genes showing peak budesonide-induced repression at 1, 2, 6, or 18 hours was summarized as in (B). For each time of peak budesonide-induced repression group, the genes were further stratified according to the greatest effect (+1, 0, or −1) of formoterol occurring at any time: +1 = ≥2-fold, 0 = >0.5/<2-fold, and −1 = ≤0.5-fold produced by formoterol addition when compared with budesonide-repressed expression. Red line = formoterol, blue line = budesonide, and purple line = budesonide + formoterol. (F) Graphical representation showing the gene count for budesonide-repressed genes meeting the criteria (+1, 0, or −1) in (E) for the effect of formoterol (Form) addition. (G) The 198 genes, which were repressed (≤0.5-fold, P ≤ 0.05) at any time by budesonide + formoterol, but not but by either formoterol or budesonide alone, in BEAS-2B cells were grouped according to the time of peak repression. The 25, 77, and 96 genes with peak repression at 2, 6, or 18 hours were summarized as in (B). Red line = formoterol, blue line = budesonide, purple line = budesonide + formoterol. (H) The 117 genes that showed enhanced repression by budesonide plus formoterol when compared with both budesonide and formoterol alone were grouped according to the time of peak repression. There were 22, 76, and 31 genes showing peak repression at 2, 6, and 18 hours, respectively. Fold data for each group were expressed as log2(fold) and are plotted as box and whiskers to show the effect of each treatment.
Glucocorticoid-Repressed Gene Expression and the Effect of LABA.
Heat maps of the 413 budesonide-repressed mRNAs highlight the few genes with rapid onset kinetics compared with the more prevalent repression at 6 or 18 hours (Fig. 4D). Indeed, 63.2% of these genes revealed repression that was maximal, or still increasing, at 18 hours. The effects of formoterol cotreatment were restricted to few genes (Fig. 4, E and F). Only 2.6% revealed further repression by formoterol, mostly due to independent repression by both treatments (Fig. 1Cii), whereas 8.7% showed budesonide-induced repression that was opposed by formoterol (Fig. 4, E and F). Of these latter 36 genes, 23 were formoterol-induced and 12 were significantly induced, albeit below the 2-fold cutoff (Supplemental Table 1). These represent formoterol-induced genes whose expression was repressed by budesonide.
Genes Only Downregulated by LABA plus Glucocorticoid.
The 198 mRNAs that were only repressed by budesonide plus formoterol revealed primarily late-onset effects, with 49% of genes showing peak repression at 18 hours (Fig. 4G). However, for each time of peak repression, effects of each treatment alone were also clearly apparent. Thus, significant formoterol- or budesonide-dependent repression that did not reach the ≤0.5-fold cutoff occurred for 142 or 186 genes, respectively, and 134 of these mRNAs were significantly repressed by both treatments.
Enhanced Repression by Formoterol plus Budesonide in Combination.
In the above sections, response additivity and either supra- or infra-additivity for formoterol and budesonide at maximally effective concentrations were used to define positive or negative cooperativity (Foucquier and Guedj, 2015). However, similar analysis is inappropriate for repressed genes where repression from fold = 1 can only approach zero. Although small repressive effects due to two stimuli may combine to give greater than additive repressive effects, more highly repressed mRNAs cannot achieve this. Instead, a highest single agent analysis (Foucquier and Guedj, 2015), which simply identifies genes that display more repression by formoterol plus budesonide compared with either treatment alone, was performed.
Of the 641 mRNAs repressed by any treatment (Fig. 1Cii), 117 showed repression by budesonide plus formoterol that was significantly greater than for budesonide and formoterol alone (Supplemental Table 3). Eight mRNAs were formoterol-repressed and 51 were budesonide-repressed, with seven mRNAs being repressed in common. The remaining 65 genes were all repressed by budesonide plus formoterol only. With 0, 22, 76, and 31 genes at 1, 2, 6, and 18 hours, enhanced repression was predominantly a delayed effect. Overall, 108 mRNAs were significantly repressed by formoterol alone, and 114 were repressed by budesonide alone. Furthermore, whereas modest repression produced by both treatments combined to give greater repression at 2 hours, at 6 and 18 hours, enhanced repression was predominantly repression due to budesonide, and this was enhanced by formoterol (Fig. 4H).
GO of LABA-Induced Gene Expression and Effect of Glucocorticoid.
As described for the LABA indacaterol (Yan et al., 2018), functional annotation of the 212 formoterol-induced genes produced enriched GO terms (PB ≤ 0.05) relating to positive and negative control of transcription, cytokines, signal transduction, inflammatory signaling, Wnt signaling, differentiation, and development (Table 1). Functional annotation clustering reinforced themes for transcriptional control, cell cycle/proliferation, growth factors, development, pluripotency, epithelial-mesenchymal transition, and cancer (Supplemental Fig. 3A). In addition, a strong inflammatory signal was apparent with clusters 4 and 9 (Fig. 5A; Supplemental Fig. 3A).
Functional annotation of the 212 genes induced by formoterol in BEAS-2B cells
The 212 genes induced ≥2-fold (P ≤ 0.05) by formoterol in BEAS-2B cells were subjected to functional annotation analysis in the DAVID. The 24 GO terms showing PB ≤ 0.05 (where 1.3 = −log100.05) are listed to show fold enrichment (Fold), PEASE and PB [both as −log10(P)], gene count (Count), and the fraction of those genes that showed further increased (≥2-fold = Fraction up) or decreased (≤0.5-fold = Fraction dn) expression at any time by budesonide. Expected fractions are as follows: Fraction up = 0.22, Fraction dn = 0.37. GO terms were manually curated into the indicated categories.
Functional aspects of formoterol-induced gene expression with effect of budesonide and overview GO analysis of genes modulated by formoterol, budesonide, and budesonide plus formoterol in BEAS-2B cells. (A) Functional annotation clustering of the 212 genes induced by formoterol resulted in 11 GO term clusters with enrichment scores ≥1.3 (Supplemental Fig. 3). (i) GO terms making up cluster 4 were depicted along with a bar graph showing the fraction of formoterol-induced genes that were further enhanced (≥2-fold) or repressed (≤0.5-fold) by budesonide. Expected fractions (enhanced = 0.22, repressed = 0.37) are indicated by dotted lines. Of the 21 genes in this cluster, three were further enhanced (≥2-fold) and 12 were repressed (≤0.5-fold) by budesonide cotreatment, making this cluster predominantly repressed by budesonide. (ii) Heat map showing effects of formoterol (Form), budesonide (Bud), and Bud + Form on the 21 genes within this cluster. Red = induced, blue = repressed, white = log2(1), and color intensity reflects increasing change in log2 fold. (B) Heap maps showing gene expression, as in (A), in response to formoterol, budesonide, and Bud + Form for all induced genes that were associated with the GO terms (i) “Pathways in cancer” and (ii) “cytokine activity” plus “Jak-STAT signaling pathway” with any treatment. (C) GO terms (biologic process, molecular function, KEGG) for the 267, 777, and 1165 genes modulated, either up (≥2-fold, P ≤ 0.05) or down (≤0.5-fold, P ≤ 0.05), by Form, Bud, or budesonide plus formoterol (B + F) in BEAS-2B cells were extracted using the DAVID. The number of individual genes within each functional grouping is indicated for each treatment. dev, development; diff, differentiation; prolif, proliferation.
GO analysis of the 46 formoterol-induced genes that were upregulated by budesonide showed enrichments (PEASE ≤ 0.1) for terms relating to transcription, development, and apoptosis. These genes are preferentially associated with the formoterol-induced gene clusters relating to development, growth factors, Wnt signaling, and cancer (Supplemental Fig. 3A). Conversely, GO annotation of the 78 formoterol-induced genes that were repressed by budesonide showed enrichment of terms related to inflammatory signaling, including “signal transduction” and “TNF signaling pathway” (with PB ≤ 0.05) and “inflammatory responses,” “positive regulation of GTPase activity,” “cytokine-cytokine receptor interaction,” and “Jak-STAT signaling pathway” (with PEASE ≤ 0.1). This suggests that budesonide had downgraded the inflammatory signature of formoterol. Repression of inflammatory responses was also indicated by the greater than expected fractions of budesonide-repressed genes within the formoterol-induced annotation (Supplemental Fig. 3A; Table 1). Thus, formoterol-induced inflammatory cytokines (IL6, IL11, IL20, LIF), chemokines (CCL2, CCL7, CXCL2), and receptors (ACKR3, CSF2RB), as per “cytokine-cytokine receptor interaction,” “Jak-STAT signaling pathway,” and “cytokine activity,” were repressed by budesonide (Fig. 5A).
GO of Glucocorticoid-Induced Gene Expression and the Effect of LABA.
Functional annotation identified 12 (PB ≤ 0.05) GO terms that were enriched for the 370 budesonide-induced genes (Table 2). “Positive regulation of apoptosis” was most significant, and various terms related to inflammation and stress, growth factor signaling, migration, differentiation, development, insulin responses, and “extracellular matrix organization” were also enriched. These themes were confirmed by functional annotation clustering, which produced clusters relating to TNF signaling/nuclear factor κB (NF-κB) and apoptosis/PI3K-Akt signaling, adhesion and extracellular matrix-receptor interaction, development, and morphogenesis (Supplemental Fig. 3B). Of note was that the 70 genes with peak expression at 1 or 2 hours revealed “negative” and “positive regulation of transcription from RNA polymerase II promoter” (14 genes each and 4.8- and 3.5-fold enrichments, respectively) as the most and third most significant terms (PB = 4.4 × 10−3 and 0.019, respectively).
Functional annotation of the 370 genes induced by budesonide in BEAS-2B cells
The 370 genes induced ≥2-fold (P ≤ 0.05) by budesonide were subjected to functional annotation analysis in the DAVID. The 12 GO terms showing PB ≤ 0.05 (where 1.3 = −log100.05) are listed to show fold enrichment (Fold), PEASE and PB [both as −log10(P)], gene count (Count), and the fraction of those genes that showed further increased (≥2-fold = Fraction up) or decreased (≤0.5-fold = Fraction dn) expression at any time by formoterol. Expected fractions are as follows: Fraction up = 0.18, Fraction dn = 0.03. GO terms were manually curated into the indicated categories.
With only 11 budesonide-induced genes repressed by formoterol, GO was unclear. Equally, no terms with PB ≤ 0.05 and only seven terms with PEASE ≤ 0.1 were identified for the 65 budesonide-induced, formoterol-enhanced genes. These included terms for apoptotic processes, development, “angiogenesis,” and “TNF signaling pathway,” a number of which appear in the budesonide-induced gene GO profile (Supplemental Fig. 3B; Table 2).
Comparative GO Analysis of LABA plus Glucocorticoid–Induced Genes.
Although the 196 genes induced only by formoterol plus budesonide produced no GO terms meeting PB ≤ 0.05, terms for transcriptional control, signaling, proliferation, development, angiogenesis, growth factors, migration, and signaling pathways (PI3K-Akt, notch, wnt, hippo) all showed PEASE ≤ 0.1. To explore the overall effects of combination treatment, GO terms showing PB ≤ 0.05 for at least one treatment were clustered in a comparative analysis (Supplemental Fig. 4). This identified the following groups of terms: A, enriched in all treatments; B, primarily enriched with budesonide, but maintained with the combination; and, C, modestly enhanced with formoterol or budesonide, but highly enriched with the combination (Supplemental Fig. 4). Collectively, these terms related to apoptosis, transcription, development, cell fate, angiogenesis, proliferation, migration, morphogenesis, differentiation, epithelial-to-mesenchymal transition, and similar effects. Furthermore, significance generally increased for all these terms with combination treatment and, as illustrated by “Pathways in cancer” (Fig. 5Bi), this was associated with increased gene counts (Supplemental Fig. 4). A fourth group of terms (group D) were most significantly enriched by formoterol alone (Supplemental Fig. 4). Examples include “cytokine activity” or “Jak-STAT signaling pathway,” found in formoterol-induced cluster 4, where proinflammatory genes induced by formoterol were repressed by budesonide (Fig. 5A). However, despite reduced enrichments with combination treatment, the gene count associated with each term was not necessarily reduced (Supplemental Fig. 4). For example, whereas budesonide reduced expression of many formoterol-induced genes in “cytokine activity” and “Jak-STAT signaling pathway” (Fig. 5Bii), expression of other formoterol-induced genes was maintained. Furthermore, budesonide also induced genes that were either unaffected by formoterol (GHR, IL12A, JAK2, PIK3R1) or were further enhanced by combination treatment (AREG, BMP2, GDF6, KITLG, MYC, VEGFA, WNT5A) (Fig. 5Bii). Thus, whereas many inflammatory genes induced by formoterol were repressed by budesonide, budesonide separately induced signaling, cytokines, and growth factor genes such that the overall quality of the response was modified with combination treatment.
GO of LABA-Repressed Gene Expression and Effect of Glucocorticoid.
Functional annotation of the 55 formoterol-repressed genes produced 20 GO terms with PEASE < 0.1, but all failed to reach PB ≤ 0.05. These terms are related to the positive and negative regulation of transcription, proliferation, differentiation, and apoptosis. Functional annotation clustering produced 4 clusters containing terms for transcriptional control, cancer, proliferation/differentiation, and stress responses (Supplemental Fig. 5A). However, this was insufficient for insight into the few genes showing further modulation by budesonide.
GO of Glucocorticoid-Repressed Gene Expression and Effect of LABA.
Inputting the 413 budesonide-repressed genes into the DAVID yielded 20 enriched GO terms (PB ≤ 0.05) (Table 3). These included terms for apoptosis, transcriptional control, cell migration, chemotaxis, angiogenesis, proliferation, growth factors, and related signaling pathways. These data were corroborated by functional annotation clustering, which produced 11 clusters with ≥1.3 enrichment scores (Supplemental Fig. 5B). Cluster 1 contained terms for transcriptional control and represented 58 genes, nearly half of which were sequence-specific transcription factors. This, along with nine additional transcription factors in cluster 4, confirms the extensive repressive effects of budesonide on transcription. Other clusters related to growth factor signaling, proliferation, and chemotaxis.
Functional annotation of the 413 genes repressed by budesonide in BEAS-2B cells
The 413 genes repressed ≤0.5-fold (P ≤ 0.05) by budesonide in BEAS-2B cells were subjected to functional annotation analysis in the DAVID. The 20 GO terms showing PB ≤ 0.05 (where 1.3 = −log100.05) are listed to show fold enrichment (Fold), PEASE and PB [both as -log10(P)], gene count (Count), and the fraction of those genes that showed increased (≥2-fold = Fraction up) or further decreased (≤0.5-fold = Fraction dn) expression at any time by formoterol. Expected fractions are as follows: Fraction up = 0.09, Fraction dn = 0.03. GO terms were manually curated into the indicated categories.
With only 11 budesonide-repressed mRNAs that were further repressed by formoterol, GO effects were modest (Supplemental Fig. 5B; Table 3). Conversely, various budesonide-repressed terms and clusters showed higher than expected fractions of the 36 budesonide-repressed genes that were enhanced by formoterol (Supplemental Fig. 5B; Table 3). These included clusters enriched with terms for proliferation, transcription, cytokine signaling, and chemotaxis. These represent functions associated with the budesonide-repressed genes that were enhanced by formoterol. This supports findings that inflammatory processes upregulated by formoterol were reduced by budesonide.
Comparative GO Analysis of Genes Repressed by LABA and Glucocorticoid in Combination.
The 198 genes that were repressed only by budesonide plus formoterol revealed a weak GO signature, which included GO terms (all PEASE ≤ 0.1) for transcriptional control and proliferation. All GO terms with PB ≤ 0.05 for at least one treatment were therefore combined in a comparative analysis (Supplemental Fig. 6). This highlighted both the modest GO signature due to formoterol-repressed genes as well as the 20 terms (PB ≤ 0.05) enriched for the budesonide-repressed genes Supplemental Fig. 6; Table 3). Although some budesonide-enriched terms were reduced in significance with combination treatment, other terms were increased. With the exception of three terms, containing formoterol-induced but budesonide-repressed genes (e.g., IL6, IL11, LIF), the number of repressed genes associated with each term increased from budesonide treatment to the combination. Thus, despite a modest repressed gene GO profile elicited by formoterol, effects of budesonide were largely maintained or enhanced by formoterol.
Overall GO of LABA- and Glucocorticoid-Modulated Genes.
As many GO terms were enriched in both induced and repressed gene lists, functional annotation was performed using combined lists of genes modulated (up or down) by each treatment in BEAS-2B cells to identify GO terms with PB ≤ 0.05 for at least one treatment (Supplemental Table 4). These were further filtered by PB ≤ 0.01 to give 86 GO terms that were manually curated into six functional groups (Supplemental Fig. 7A). Thus, transcriptional control terms, including positive and negative regulation, were highly enriched with all treatments. This group of terms contained 77, 157, and 238 genes with respect to formoterol, budesonide, and combination treatment, respectively, and represents a central feature of each response (Fig. 5C).
Numerous signaling and signal transduction terms, for example, the generic term “signal transduction,” were enriched with each treatment and showed the highest gene numbers with combination treatment (Supplemental Fig. 7A, Supplemental Table 4). Likewise, pathway-specific signaling terms, such as “PI3K-Akt signaling pathway” or “MAPK signaling pathway,” although containing fewer genes, revealed the same trend (Supplemental Fig. 7A, Supplemental Table 4). The greatest levels of significance were obtained for “TNF signaling pathway.” With 14, 20, and 27 genes regulated by formoterol, budesonide, and formoterol plus budesonide, respectively, this term illustrates the interplay between genes induced and repressed by each treatment. Thus, formoterol induced proinflammatory mRNAs (e.g., CCL2, CXCL2, EDN1, FOS, IL6, IL15, LIF), many of which (e.g., CCL2, FOS, LIF, IL6) were independently repressed by budesonide and/or showed repression upon combination treatment (Supplemental Fig. 7B). However, other “proinflammatory,” formoterol-induced mRNAs were relatively unaffected (CXCL2, IL15) or were enhanced (EDN1) by budesonide. Budesonide also induced mRNAs for TNFAIP3 and NFKBIA and other genes that negatively regulate NF-κB. In the case of NFKBIA, or inhibitor of κBα, budesonide-induced expression was unaffected by formoterol, whereas TNFAIP3, or A20, as previously reported (Altonsy et al., 2017), was induced by formoterol alone and more highly induced by combination treatment. Thus, downregulation of cytokines and chemokines by budesonide combined with the upregulation of repressors of NF-κB, especially in the context of LABA, may reduce proinflammatory effects. Furthermore, reduced expression of transcription factors (e.g., JUN, CREB3L1, CREB5) and signaling components (e.g., MAP2K3, MAP3K14, TRAF3) could also attenuate inflammatory processes. Conversely, induction of EDN1, a proinflammatory G-protein–coupled receptor agonist, or MAP3K8, an upstream MAPK kinase kinase, by formoterol and budesonide, and their further enhancement by the combination, implies complexity.
In terms of associated gene counts, these data show the effects of budesonide to be generally increased by combination treatment (Fig. 5C; Supplemental Fig. 7A). Signaling and transcriptional control were the second and third most highly represented groups by gene number, and the proliferation, differentiation, and development group of terms was represented by 354 genes with combination treatment (Fig. 5C). This represents the dominant effect, by gene and GO term number, of combination treatment. Nevertheless, cell adhesion and migration (171 genes) or inflammation and stress (142 genes) were also major components of the response to combination treatment and, in each case, showed increased gene numbers, and general significance, with combination treatment.
Formoterol- and Budesonide-Modulated Transcriptomes in pHBECs.
In pHBECs, formoterol (10 nM), budesonide (100 nM), or formoterol plus budesonide for 6 hours induced expression of 140, 94, and 243 genes and repressed 208, 141, and 370 genes, respectively (Fig. 6A; Supplemental Table 5). This revealed greater than predicted overlap between genes induced, or repressed, by both treatments (Fig. 6, B and C). With formoterol plus budesonide, 50.7% of formoterol-induced mRNAs were not induced (Fig. 6Ci). Furthermore, 17.9% of the formoterol-induced mRNAs showed ≥50% repression by budesonide, whereas 7.9% were further enhanced (≥2-fold) (Fig. 6D). Thus, repression was the main effect of budesonide on formoterol-induced gene expression. Conversely, few budesonide-induced mRNAs were repressed (≥50%) by formoterol, whereas 15.9% showed further enhancement (≥2-fold). Enhancement was therefore the main effect of formoterol on budesonide-induced expression (Fig. 6D). With combination treatment, 109 genes (44.9%) were induced but failed the induction criteria for either treatment alone (Fig. 6Ci). Of these, 74 were significantly induced by formoterol and 86 by budesonide, with 59 mRNAs induced by both treatments, albeit below 2-fold. Thus, modest inducibility by each treatment in fact characterizes those genes only significantly (≥2-fold) induced by combination treatment.
Effects of formoterol and budesonide alone, or in combination, on gene expression in pHBECs. pHBECs from six normal individuals were either not stimulated (NS) or treated with formoterol (Form; 10 nM) and budesonide (Bud; 100 nM) or both together (Bud + Form or B + F) for 6 hours prior to RNA extraction and microarray analysis. Following robust multiarray averaging and normalization, data were expressed as fold relative to NS, and the genes showing ≥2-fold or ≤0.5 fold (P ≤ 0.05) were considered to be differentially expressed (Supplemental Table 5). (A) The number of differentially expressed genes showing up- or downregulation for each treatment is shown. Overall gene numbers of modulated (up- and downregulated) genes are shown below the graph. (B) Expression of genes showing different expression by formoterol was plotted against the effect of budesonide alone, and expression of the genes showing differential expression to budesonide was plotted against the effect of formoterol alone. Gene numbers in each corner quadrant are shown along with expected numbers in brackets. (C) The number of genes meeting the expression criteria for (i) induction (≥2-fold, P ≤ 0.05) or (ii) repression (≤0.5-fold, P ≤ 0.05) at any time with respect to each treatment is represented as Venn diagrams. Beside each Venn diagram, contingency tables assess the likelihood (χ2 test) that the number of genes (induced or repressed) in common between formoterol and budesonide may have occurred by chance. ****P = 0.0001. (D) The 140 formoterol-induced genes were plotted to show the number of genes that were further induced (≥2-fold) or repressed (≤0.5-fold) with budesonide cotreatment, and the 94 budesonide-induced genes were plotted to show the number of genes that were further induced (≥2-fold) or repressed (≤0.5-fold) by formoterol cotreatment. (E). The 208 formoterol-repressed genes were plotted to show the number of genes that were induced (≥2-fold) or further repressed (≤0.5-fold) with budesonide cotreatment, and the 141 budesonide-repressed genes were plotted to show the number of genes that were induced (≥2-fold) or further repressed (≤0.5-fold) by formoterol cotreatment. E, expected, O, observed.
Turning to repression, 42 formoterol-repressed mRNAs were not repressed with budesonide plus formoterol, and eight (3.8%) of these were increased by ≥2-fold (Fig. 6, Cii and E). Expression of 19 (9.1%) formoterol-repressed genes was further reduced (≥50%) by budesonide cotreatment (Fig. 6E). A more polarized effect was apparent for the budesonide-repressed genes, in which only 15 failed to meet the criteria for repression with budesonide plus formoterol and just one gene (0.7%) revealed a ≥2-fold increase (Fig. 6Cii and E). Conversely, 38 (27%) of the budesonide-repressed genes showed ≥50% further repression with formoterol (Fig. 6E). Thus, although repression produced by formoterol was enhanced by budesonide, budesonide-dependent repression was more widely enhanced by formoterol. In addition, with formoterol plus budesonide, 162 genes reached the criteria for repression that was not reached with either treatment alone (Fig. 6Cii). Furthermore, with 125 genes significantly repressed by each treatment and 106 being significantly repressed by both treatments, modest repression by formoterol and/or budesonide predicts combinational effects (Supplemental Table 5).
Induced Genes Show Supra- and Infra-additive Effects of LABA/Glucocorticoid Combination in pHBECs.
Of the 328 genes induced by formoterol, budesonide, or the combination, 55 showed significant supra-additivity, whereas 35 showed significant infra-additivity (Supplemental Table 5).
Of those genes showing supra-additivity, five were induced genes common to both treatments, and 34 were genes induced only by formoterol plus budesonide. Similarly, of the 79 budesonide-induced genes that were not formoterol-induced, 11 genes showed significant supra-additivity, whereas only four revealed infra-additivity (Supplemental Table 5). Thus, enhancement of budesonide-induced gene expression by formoterol was more common than repression. This was compared with the 125 genes that were only induced by formoterol, of which, 6 revealed significant supra-additivity and 30 showed infra-additivity in the presence of budesonide. Thus, the major effect of budesonide on formoterol-induced genes was to reduce gene expression.
Enhanced Repression in pHBECs by Formoterol plus Budesonide.
Assessment of combinatorial effects between budesonide and formoterol on repressed genes was restricted to identifying mRNAs that showed more repression with the combination compared with each treatment alone. Of the 429 genes repressed by any treatment, 43 were significantly more repressed by budesonide plus formoterol compared with either treatment (Supplemental Table 5). In total, 41 of these mRNAs were significantly repressed by budesonide alone, and 38 were significantly repressed by formoterol alone. Thus, repression by each treatment may combine to produce enhanced repression with combination treatment (Fig. 6F).
GO Associated with Budesonide- and Formoterol-Induced Genes in pHBECs.
Functional annotation of the genes induced by formoterol, budesonide, or budesonide plus formoterol produced weak enrichments with just a single term, “cellular response to insulin stimulus,” in the budesonide-induced gene list, reaching PB ≤ 0.05. Numerous terms for each treatment met the lower (PEASE ≤ 0.1) level of significance, including terms for positive and negative transcriptional control and signaling. Formoterol-induced genes showed enrichment for development, differentiation, positive and negative regulation of proliferation, wounding responses, cancer, apoptosis, and metabolism. Similarly, the budesonide-induced list revealed terms for cell cycle, differentiation, and apoptosis, as well as signaling terms relating to metabolism, inflammation, and proliferation. Functional annotation for the 243 budesonide plus formoterol–induced genes also introduced terms for development, differentiation, proliferation, apoptosis, cancer, and many related pathways.
GO Associated with Budesonide- and Formoterol-Repressed Genes in pHBECs.
With the 208 formoterol-repressed genes, functional annotation yielded GO terms (PB ≤ 0.05) that related to protein phosphatase activity and “inactivation of MAPK activity,” plus a fifth term, “inflammatory response” (Table 4). Various terms for transcription, apoptosis, development, positive and negative proliferation, migration and chemotaxis, differentiation, and “positive regulation of GTPase activity” all met the lower PEASE ≤ 0.1 threshold.
Functional annotation of the 208 genes repressed by formoterol in pHBECs
The 208 genes repressed ≤0.5-fold (P ≤ 0.05) by formoterol in pHBECs were subjected to functional annotation analysis in the DAVID. The five GO terms showing PB ≤ 0.05 (where 1.3 = −log100.05) are listed to show fold enrichment (Fold), PEASE and PB [both as −log10(P)], and gene count (Count). GO terms were manually curated into the indicated categories.
The 141 budesonide-repressed genes produced 23 terms with enrichments of PB ≤ 0.05 (Table 5). Of these, “growth factor activity” was most significant, with other terms relating to proliferation, cell division, and angiogenesis. Likewise, “inflammatory response” was highly enriched, and this linked with “cytokine activity,” “immune response,” and “cytokine-cytokine receptor interaction.” Together with pathway terms for “NF-kappa B signaling pathway” and “TNF signaling pathway,” these terms suggest repression of inflammatory signals. However, “inactivation of MAPK activity” was also highly enriched, and since the term includes multiple phosphatases, this rather paradoxically suggests reduced control of MAPK activity.
Functional annotation of the 141 genes repressed by budesonide in pHBECs
The 141 genes repressed ≤0.5-fold (P ≤ 0.05) by budesonide in pHBECs were subjected to functional annotation analysis in the DAVID. The 23 GO terms showing PB ≤ 0.05 (where 1.3 = −log100.05) are listed to show fold enrichment (Fold), PEASE and PB [both as −log10(P)], and gene count (Count). GO terms were manually curated into the indicated categories.
With the 370 budesonide plus formoterol–repressed genes, 8 of the 22 GO terms showing PB ≤ 0.05 related to transcriptional control (Table 6). Other terms, including “growth factor activity,” correlated with cell proliferation, angiogenesis, and apoptosis, whereas “cytokine activity” and “inflammatory response” suggest repression of inflammation. However, “inactivation of MAPK activity” and “negative regulation of ERK1 and ERK2 cascade” again indicate complexity in determining net function.
Functional annotation of the 370 genes repressed by formoterol plus budesonide in pHBECs
The 370 genes repressed ≤0.5-fold (P ≤ 0.05) by formoterol plus budesonide in pHBECs were subjected to functional annotation analysis in the DAVID. The 22 GO terms showing PB ≤ 0.05 (where 1.3 = −log100.05) are listed to show fold enrichment (Fold), PEASE and PB [both as −log10(P)], and gene count (Count). GO terms were manually curated into the indicated categories.
Overall GO for Genes Modulated by Formoterol, Budesonide, or the Combination in pHBECs.
The weak GO signatures produced by the induced gene lists and overlap with the repressed gene GO prompted combined functional annotation analysis of all up- and downregulated genes. This produced 65 GO terms (with PB ≤ 0.05 for at least one treatment), which were manually curated into six functional groups (Supplemental Fig. 8; Supplemental Table 6). Positive and negative terms for transcriptional control were significantly enriched with all treatments, with greatest significance and gene number with combination treatment (Supplemental Fig. 8). However, despite almost 150 genes relating to transcriptional control with budesonide plus formoterol, this was exceeded by proliferation, differentiation, and development, which represented 208 genes within 20 GO terms and was the largest functional group for each treatment (Fig. 7A; Supplemental Fig. 8). Again, overall significance and gene number was generally greatest with combination treatment. For example, “growth factor activity” contained 10, 14, and 17 significantly induced or repressed genes for formoterol, budesonide, or formoterol plus budesonide, respectively (Fig. 7B; Supplemental Table 6). Visualization of gene expression changes revealed not only the extent to which many growth factors were repressed by formoterol and budesonide but also how repression can be enhanced by combination treatment (Fig. 7B). Conversely, although various factors (e.g., AREG, BDNF, EPGN) revealed formoterol-induced expression that was reduced by combination treatment, other genes (e.g., DKK1, INHBB) showed increased expression.
GO analysis of genes modulated by formoterol, budesonide, and budesonide plus formoterol in pHBECs. (A) GO terms (biologic process, molecular function, KEGG) for the 348, 235, and 613 genes modulated, either up (≥2-fold, P ≤ 0.05) or down (≤0.5-fold, P ≤ 0.05), by formoterol, budesonide, or formoterol plus budesonide in pHBECs were extracted using the DAVID. The number of individual genes within each functional grouping are indicated for each treatment. (B) Heat map showing effects of formoterol (F), budesonide (B), and budesonide + formoterol (BF) on the expression of all those genes, within the indicated GO terms, that were significantly modulated by any treatment. Red = enhanced fractions and blue = repressed. (C) Expression data for selected genes showing significant enhancement by formoterol. Microarray data for the indicated genes were plotted as log2 fold. Effects of formoterol (Form), budesonide (Bud), and budesonide + formoterol (B + F) are shown. Significance (ANOVA) is indicated where *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001. dev, development; diff, differentiation; prolif, proliferation.
Related to the effects on proliferation and growth were GO terms in the adhesion and migration group (Supplemental Fig. 8). Although representing fewer overall terms and genes, these showed the same trend toward greatest significance and number of genes with formoterol plus budesonide (Fig. 7A; Supplemental Fig. 8). This effect was less clear for the inflammation and stress group, in which enrichments were consistently most significant with budesonide treatment. Furthermore, although the number of genes associated with formoterol plus budesonide for each term was generally higher, or at least equal, to those obtained for budesonide alone, some terms, for example “Cytokine-cytokine receptor interaction” and “cytokine activity,” showed reduced gene numbers with combination treatment (Fig. 7B). For illustration, CXCL2 and AREG were induced by formoterol and repressed by budesonide, but in combination, these effects largely cancelled out (Fig. 7, B and C). Thus, multiple inflammatory genes, including cytokines (TSLP), chemokines (CXCL3, CXCL5, CXCL6), chemokine receptors (ACKR3 and CXCR4), and tyrosine kinase receptor ligands (EPHA4, EFNA1) were induced by formoterol, repressed by budesonide alone, and with the combination, were repressed relative to formoterol (Fig. 7, B and C). Furthermore, repression of many inflammatory genes, including cytokines (IL1A, IL1B, IL11, IL24, IL36G, TSLP) and chemokines (CXCL1, CXCL2, CXCL8), was widely observed with budesonide, and in combination with formoterol, repression was maintained or enhanced (Fig. 7B). Conversely, other inflammatory products, receptors, or growth factors, albeit fewer in number (see Fig. 7B), revealed budesonide- and or formoterol-induced expression that was maintained or even increased by combination treatment.
Finally, effects of formoterol and budesonide on GO for signaling generally conformed to the pattern of increasing significance and associated gene number with combination treatment (Fig. 7A; Supplemental Fig. 8). Furthermore, many of the signaling terms related to the other functional groupings, including cell proliferation and movement, differentiation, development, inflammation, and stress. More difficult to explain are GO terms for phosphatases and inactivation of MAPK pathways that were prominent among the genes repressed by formoterol and budesonide plus formoterol (Tables 4 and 6). Visualization of gene expression changes associated with “inactivation of MAPK activity” and “protein tyrosine phosphatase activity” revealed six dual-specificity phosphatases (DUSPs), which target MAPKs (Jeffrey et al., 2007), five protein phosphatases, and three protein tyrosine phosphatase receptors (Fig. 7B). Since a majority (13 of 19) of these genes showed formoterol-dependent repression that, with the exception of DUSP10, was maintained or enhanced with budesonide plus formoterol, enhanced protein phosphorylation is predicted. Conversely, as DUSP1 represses MAPK activity, and expression was increased by budesonide and enhanced by formoterol, cross-regulation may be important.
Discussion
The GINA guidelines now recommend ICS/LABA combination therapy as reliever and/or maintenance in all patients with asthma (https://ginasthma.org/reports/). This revised approach should reduce mortality in individuals with asthma who take insufficient ICS but overuse short-acting β2-adrenoceptor agonists for symptom relief. However, despite the clinical superiority of ICS/LABA combination therapies over ICS alone, a molecular basis remains unclear. We hypothesized involvement of transcriptional interactions, and this was tested in airway epithelial cells, a primary therapeutic target for these drugs. Widespread gene expression changes occurred in response to formoterol and budesonide, with genes showing in common up- and downregulation by both treatments being over-represented. Similarly, genes upregulated by formoterol but repressed by budesonide were also over-represented. Such data suggest biologic relevance. Indeed, the effect of formoterol plus budesonide was not simply a summation of the responses by each component; formoterol modified responses to budesonide, and budesonide modified responses to formoterol. Furthermore, ∼one-third of the genes modulated by formoterol plus budesonide only met the criteria for regulation with combination treatment. This outcome could not have been anticipated from the actions of each drug alone. Thus, effects of formoterol and budesonide alone, or in combination, are defined by gene populations, each with distinct “functional” signatures. Moreover, commonality between glucocorticoid-regulated genes, and their associated GO profiles in BEAS-2B, pHBECs, and the human airways post–budesonide inhalation (Leigh et al., 2016; Mostafa et al., 2019), suggest that ICS/LABA interactions are pertinent to asthma therapy.
Formoterol generated a gene induction profile that reflected rapid-onset (1 to 2 hours) cAMP signaling (Mayr and Montminy, 2001; Zhang et al., 2005). GO analysis indicated transcription and signaling as major effector responses in which multiple terms specified outputs associated with inflammation. This confirms studies showing β2-agonists promote inflammatory mediator expression in structural cells (Ammit et al., 2002; Edwards et al., 2007; Strandberg et al., 2007), including in HBECs (Yan et al., 2018). One example from pHBECs was TSLP, a cytokine associated with asthma pathogenesis (Mitchell and O’Byrne, 2017), which was upregulated by formoterol but then repressed upon budesonide cotreatment. Similar effects occurred with other cytokines, chemokines, growth factors, and/or their receptors. Furthermore, despite inflammatory stimuli augmenting proinflammatory effects of β2-agonists, repression by glucocorticoid persisted and suggests benefits for combination therapies in asthma (Ammit et al., 2002; Edwards et al., 2007; Holden et al., 2010).
Compared with formoterol, most budesonide-induced gene expression changes occurred after a lag of several hours. In BEAS-2B cells, expression of 70 induced genes peaked at 1 to 2 hours, whereas 300 mRNAs showed peak expression at 6–18 hours. Similarly, 393 of the 413 budesonide-repressed genes showed peak repression at 6–18 hours. Thus, the main responses to budesonide were delayed and are consistent with early GR-dependent transcription eliciting later-onset activation and repression (Chinenov et al., 2014; Newton and Holden, 2007; Reddy et al., 2009; Sasse and Gerber, 2014). Indeed, GO analysis of budesonide-induced genes at 1 to 2 hours showed positive and negative regulation of transcription among the most highly enriched terms. This is consistent with transcriptional activators and repressors being glucocorticoid-inducible in multiple cell types, including epithelial cell lines (BEAS-2B, A549), pHBECs, and the human airways post–budesonide inhalation (Chinenov et al., 2014; Himes et al., 2014; Leigh et al., 2016; Kan et al., 2019; Mostafa et al., 2019). Nevertheless, other control mechanisms also regulate downstream gene expression. For example, in BEAS-2B cells (this study), A549 cells, and the human airways, glucocorticoids induced expression of the mRNA destabilizing protein ZFP36 (Smoak and Cidlowski, 2006; King et al., 2009; Leigh et al., 2016). This will reduce inflammatory gene expression (Clark and Dean, 2016). However, in BEAS-2B cells, budesonide induced negative regulators of inflammatory signaling and gene expression, including DUSP1, NFKBIA, and TNFAIP3; other glucocorticoid-induced genes (JAK2, MAP3K8, IL12A, and EDN1) were proinflammatory. Thus, GO terms relating to inflammation may reflect mixed pro- and anti-inflammatory effects.
In BEAS-2B cells, numerous budesonide-induced genes related to growth factor responses, proliferation, differentiation, development, and migration. The strength of these associations (enrichment and gene number) was increased when induced and repressed gene lists were combined, suggesting that both induction and repression contribute to net biologic function. When formoterol and budesonide were used in combination, these GO enrichments and associated gene numbers generally increased relative to either treatment alone. This occurred in both BEAS-2B cells and pHBECs, strengthening the idea that glucocorticoid alone, or in combination with LABA, can modulate proliferation, differentiation, development, adhesion, and migration. Indeed, growth factor expression was generally suppressed by glucocorticoids in cell lines and pHBECs (Mostafa et al., 2019) and is consistent with inhibition of proliferation (Bird et al., 2015). However, GO terms linking to both positive and negative regulation of proliferation and apoptosis were enriched with formoterol and budesonide treatments. As these were enhanced by combination treatment, overall function may be complex to unravel.
As described above, the principal effects of glucocorticoid were delayed, with the greatest numbers of budesonide-induced or -repressed genes occurring at 6–18 hours. At these times, formoterol-modulated gene expression changes were past their peak and, therefore, modest. However, in combination with budesonide, formoterol markedly enhanced the number of genes showing late induction or repression. Explanations for this are multiple, but as β2-adrenoceptor (ADRB2) mRNA expression was not markedly induced by any of the treatments in BEAS-2B or pHBECs, and glucocorticoid had no significant effect on LABA-induced signaling or CRE reporter activity (Rider et al., 2018), these effects were likely not due to enhancement of β2-adrenoceptor signaling by glucocorticoids. More plausible are interactions between glucocorticoid- and cAMP-activated and/or -induced transcription factors. In this respect, formoterol does not generally enhance GR binding to target genes in BEAS-2B cells (Rider et al., 2018). However, cAMP activates multiple transcription factors (Sassone-Corsi, 2012). Thus, GR and cAMP-activated factors may interact to regulate specific genes. Indeed, the over-representation of coregulated genes supports this possibility. However, regulation by both budesonide and formoterol was not synonymous with positive cooperativity; many genes showed no more than simple additivity. For example, in BEAS-2B cells and airway smooth muscle cells, DUSP1 is independently induced by LABA and glucocorticoid, and these effects summate in combination (Kaur et al., 2008; Manetsch et al., 2012, 2013). This not only illustrates how glucocorticoid-induced gene expression can promote repression of transcripts whose expression requires MAPK activity (Clark et al., 2008; Shah et al., 2014) but also shows how repression could be enhanced by LABA. However, clinical data suggest more than simple additivity between LABAs and ICSs (Newton and Giembycz, 2016). Indeed, the current study identifies numerous genes that are regulated in a supra-additive manner by budesonide and formoterol. Furthermore, analysis of CRE and glucocorticoid response element reporters or gene expression in BEAS-2B cells shows that 10 nM formoterol and 100 nM budesonide, as used herein, are maximally effective and markedly exceed the affinities for their respective receptors (Kaur et al., 2008; Rider et al., 2011, 2018; Holden et al., 2014; Alexander et al., 2017; Yan et al., 2018). Thus, supra-additive responses produced at saturating concentrations of formoterol and budesonide meet a simple definition of synergy, one that is achieved by multiple genes. These included transcription factors, solute carriers, cell cycle regulators, and metabolic enzymes, for which functional roles in the context of asthma remain unclear. Although such genes require investigation, one synergistically upregulated gene that illustrates the therapeutic superiority of ICS/LABA combination therapy in asthma is, as mentioned above, the bronchoprotective gene RGS2. Aside from reducing proinflammatory effects of LABAs, ICS acting on the airway epithelium may mimic endogenous glucocorticoids in regulating stress responses, healing, and repair (Busillo and Cidlowski, 2013). Furthermore, glucocorticoid/LABA interactions may replicate physiologic stress responses whereby endogenous adrenaline and cortisol interact to control the consequences of damage or insult. Focused mechanistic studies are therefore essential to tease out the functional properties of these commonly used asthma therapies, especially when used in combination. As glucocorticoids promote lung maturation, regulate proliferation and apoptosis, and may stimulate epithelial differentiation (Bird et al., 2015), such questions are important. Excessive modulation of these effects may be undesirable in a ubiquitously used therapy, whereas promotion of healing and repair could be beneficial in asthma. Since cooperative interactions between LABAs and glucocorticoids are class-specific (Newton and Giembycz, 2016), the current findings should apply to all clinically used ICS/LABA combinations. Understanding these effects may enable improvements to the effectiveness of ICS/LABA combinations indicated at each stage of the GINA guidelines. Finally, since β2-agonist– and ICS-dependent gene expression changes are apparent in transcriptomic analyses of severe asthma (Djukanović, 2019; Weathington et al., 2019), the current analysis helps distinguish therapy-dependent changes in gene expression from disease phenotype. Since the current study focused on ICS/LABA interactions in vitro, further investigations are required in more physiologically and pathologically relevant systems. For instance, the influence of inflammatory stimuli, the presence of serum, and the use of primary cells from patients with mild, moderate, and severe disease should be assessed. Likewise, understanding the transcriptomic signature due to ICS/LABA in patients, for which pharmacokinetics, differences in severity, or endotypes of asthma may affect responsiveness, will be required for the identification of gene expression changes that associate with therapeutic benefit.
Authorship Contributions
Participated in research design: Mostafa, Rider, Giembycz, Newton.
Contributed new reagents or analytic tools: Leigh, Newtown.
Conducted experiments: Mostafa, Rider.
Performed data analysis: Mostafa, Rider, Wathugala, Newton.
Wrote or contributed to the writing of the manuscript: Mostafa, Rider, Wathugala, Leigh, Giembycz, Newton.
Footnotes
- Received August 16, 2020.
- Accepted December 7, 2020.
↵1 Current affiliation: Respiratory Medicine, University of British Columbia, Vancouver, British Columbia, Canada.
This study was supported by Canadian Institutes for Health Research (CIHR) project grants [MOP 125918, PJT 156310] and by a Natural Sciences and Engineering Research Council (NSERC) Discovery Grant [RGPIN-2016-04549]. M.M.M. received NSERC Postgraduate Scholarship-Doctoral, Elizabeth II Doctoral scholarship, and The Lung Association – Alberta and Northwest Territories studentship awards. Real-time polymerase chain reaction (PCR) machines was facilitated by an equipment and infrastructure grant from the Canadian Foundation for Innovation and the Alberta Science and Research Authority.
↵
This article has supplemental material available at molpharm.aspetjournals.org.
Abbreviations
- DAVID
- Database for Annotation Visualization and Integrated Discovery
- DUSP
- dual-specificity phosphatase
- FGF
- fibroblast growth factor
- GINA
- Global Initiative for Asthma
- GO
- gene ontology
- GR
- glucocorticoid receptor
- HBEC
- human bronchial epithelial cell
- ICS
- inhaled corticosteroid
- JAK
- Janus kinase
- KEGG
- Kyoto Encyclopedia of Genes and Genomes
- LABA
- long-acting β2-adrenoceptor agonist
- MAPK
- mitogen-activated protein kinase
- NF-κB
- nuclear factor κB
- pHBEC
- primary human bronchial epithelial cell
- PI3K
- phosphatidylinositol 3-kinase
- STAT
- signal transducer and activator of transcription
- Copyright © 2021 by The American Society for Pharmacology and Experimental Therapeutics
