Abstract
Background Bronchial thermoplasty is a nonpharmacological, device-based treatment option for a specific population of severe asthmatic subjects, but the underlying mechanisms are largely unknown. The purpose of this study is to identify potential altered pathways by bronchial thermoplasty using a transcriptomic approach.
Methods Patients undergoing bronchial thermoplasty were recruited to the study, and a bronchial brushing sample was obtained before each bronchial thermoplasty session and sent for RNA sequencing. A variance component score test was performed to identify those genes whose expression varied after bronchial thermoplasty sessions. Differential gene expression meta-analysis of severe asthmatic subjects versus controls was performed using public repositories. Overlapping genes were included for downstream pathway and network analyses.
Results 12 patients were enrolled in our study. A total of 133 severe asthma cases and 107 healthy controls from the public repositories were included in the meta-analysis. Comparison of differentially expressed genes from our study patients with the public repositories identified eight overlapping genes: AMIGO2, CBX7, NR3C2, SETBP1, SHANK2, SNTB1, STXBP1 and ZNF853. Network analysis of these overlapping genes identified pathways associated with neurophysiological processes.
Conclusion We have shown that bronchial thermoplasty treatment alters several gene networks that are important in asthma pathogenesis. These results potentially elucidate the disease-modifying mechanisms of bronchial thermoplasty and provide several targets for further investigation.
Abstract
Gene networks that are altered in the airway epithelium of asthmatic subjects after bronchial thermoplasty http://ow.ly/Ki6x30nbPSj
Introduction
In recent years, considerable debate has arisen over how to best phenotype and treat patients with severe asthma. While biological therapies are increasingly used to treat such patients, bronchial thermoplasty is an effective [1], nonpharmacological treatment option for a subset of severe asthmatic subjects whose symptoms are not controlled [2].
The overall modest effectiveness of bronchial thermoplasty for patients with moderate to severe asthma has been shown in multiple randomised clinical trials [3–9]. The end-points in these studies include improvement in asthma control, fewer asthma exacerbations, fewer hospital and emergency department visits, and a reduction in inhaled steroid dose for patients treated with bronchial thermoplasty. The high degree of heterogeneity in the response to bronchial thermoplasty suggests this therapy may be highly beneficial in specific asthma subsets. A better understanding of the mechanisms responsible for response to bronchial thermoplasty will provide better targeting of this therapy for the patients most likely to have a robust clinical response.
There are several proposed mechanisms that partly explain the effectiveness of bronchial thermoplasty [2]. It is known to reduce smooth muscle mass in bronchial biopsy samples as measured by morphometric analysis of α-actin staining. Importantly, this effect was observed in the contralateral or untreated airways, suggesting bronchial thermoplasty treatment has a global effect on the lung [10–12]. For example, bronchial thermoplasty may reduce global airway innervation that disrupts pathological neurophysiological processes in the lung. A recent study evaluated bronchial biopsy specimens from 15 severe asthmatic subjects before and 3 months after completion of three bronchial thermoplasty treatments [13]. The authors found that bronchial thermoplasty treatment decreased airway smooth muscle area, subepithelial basement membrane thickness, the number of neuroendocrine epithelial cells and the number of bronchial nerve fibres [13]. A similar study demonstrated that eosinophils and the concentration of inflammatory mediators, such as tumour necrosis factor-β and the chemokine ligand CCL5, were reduced in parallel with reduced smooth muscle mass after bronchial thermoplasty [11]. Another proposed mechanism is the reduction of smooth muscle contractile properties after exposure to high temperatures in airways treated with bronchial thermoplasty [14]. Collectively, these studies begin to identify a series of global effects on the lung following site-specific bronchial thermoplasty treatments.
Our hypothesis was that bronchial thermoplasty induces global effects on airway smooth muscle that disrupt airway smooth muscle/epithelial cell cross-talk. This disruption results in decreased inflammatory chemokine and cytokine expression measured at the airway epithelium. We further hypothesised that decreased signalling by the airway epithelium would decrease inflammatory cell influx into the airway, mitigating asthma pathogenesis. In support of our hypothesis, several studies have shown that cross-talk between epithelial cells and other cells can modulate inflammatory processes related to type 2 immune signalling [15] and IgE-mediated signalling through periostin [16].
The purpose of this study was to advance our understanding of the pathways and networks affected during bronchial thermoplasty through a transcriptomic approach by combining a meta-analysis of public repositories with a time-course longitudinal analysis of our bronchial thermoplasty dataset.
Methods
Study population
A total of 12 participants with severe asthma who underwent bronchial thermoplasty at the University of California Davis (Sacramento, CA, USA) were recruited for the study. The study was conducted between April 2013 and November 2015. Eligible subjects were adults (18–65 years of age) diagnosed with severe asthma with indications for bronchial thermoplasty defined by the US Food and Drug Administration. All of our subjects were using inhaled corticosteroids, an oral corticosteroid or both at the time of enrolment.
Airway cell sampling
Subjects recruited to the study consented to sampling of their airway using a disposable cytology brush with direct visualisation prior to their bronchial thermoplasty procedure. A standard 5 mm cytology brush (Olympus, Center Valley, PA, USA) was used to acquire brush samples targeting the tracheal wall prior to each bronchial thermoplasty session. The sample was placed in a 15 mL conical tube with 3 mL of RNAlater (ThermoFisher Scientific, Waltham, MA, USA) storage solution for gene expression analyses. Samples obtained with this technique are highly enriched for airway epithelial cells (>95% of all cells collected; data not shown).
RNA processing and sequencing from airway brush samples
RNA was isolated from the cytology brushes using the Qiagen RNeasy Micro Kit (Qiagen, Germantown, MD, USA). cDNA Libraries were prepared with the KAPA Stranded RNA-Seq Kit with the RiboErase Kit (Roche, Pleasanton, CA, USA) according to the manufacturer's protocol. Generation of cDNA libraries and sequencing was carried out by the DNA Technologies and Expression Analysis Core Laboratories at the University of California Davis Genome Center. Genes with <0.5 counts per million reads in all samples were filtered prior to the analysis. Data were transformed using a variance-stabilising transformation and then LOESS (locally estimated scatterplot smoothing) normalised. Gene annotation was based on Gencode genome assembly version GRCh38 (www.gencodegenes.org).
Systemic search for severe asthmatic subjects versus healthy controls datasets
Two public gene expression microarray repositories (Gene Expression Omnibus database (www.ncbi.nlm.nih.gov/geo) and ArrayExpress (www.ebi.ac.uk/arrayexpress)) were utilised for our study, which included all available datasets published prior to November 1, 2017. We included only those studies derived from adult human bronchial epithelium gene expression datasets that contained both severe asthmatic subjects and healthy controls. For datasets that contained a mixed cohort of asthmatic subjects, we selected subjects that met the criteria for severe asthma. We excluded datasets that met our initial inclusion criteria if asthma severity for each individual could not be determined. We reserved one dataset for validation purposes.
Statistical analyses
The schematic of our analysis work is shown in figure 1.
Modelling gene expression changes over sequential bronchial thermoplasty sessions
The longitudinal gene expression change over sequential bronchial thermoplasty sessions was tested with the variance component score test using the R package Time-Course Gene Set Analysis for RNA-Seq Data (tcgsaseq) [17]. This method is a linear mixed effect model to account for repeat measurements from the same participants and a time variable is included. The model was adjusted for age, sex and smoking status. We tested which genes demonstrated a change in gene abundance after bronchial thermoplasty, with the null hypothesis that there were no longitudinal changes in normalised gene expression over time. Compared with traditional two-time-point gene expression comparisons (e.g. time 2 versus time 1 and time 3 versus time 1), this method takes longitudinal trajectories into consideration with full usage of all gene expression information. With time-course modelling, we are able to analyse changes in gene expression patterns over time rather than differences between two discrete time-points. Due to the small sample size, 1000 permutations were performed to calculate the p-value and a Benjamini–Hochberg false discovery rate (FDR) correction was used for multiple comparisons. The threshold for a significant change in gene expression after serial bronchial thermoplasty sessions was set as an adjusted p-value <0.001.
Differential expression meta-analysis
Gene expression datasets were downloaded from public repositories for differential expression meta-analysis to be compared with our patient cohort. Datasets were then split into discovery and validation sets. Meta-analysis was run on the discovery dataset using the R package MetaIntegrator, as described previously [18]. In brief, the meta-analysis computes a Hedges’ g effect size for each gene in each dataset. The summary effect size is computed using a random effects model. The interdataset variation was estimated by the DerSimonian–Laird method. We only included genes that remained significant across the leave-one-dataset-out analyses. The reported p-value of each gene was adjusted by Benjamini–Hochberg FDR correction for multiple comparisons. We choose adjusted p-values <0.05 as our threshold for significance. The differentially expressed genes identified as significant in the discovery datasets were then validated in the validation dataset. Performance of the model to distinguish severe asthmatic subjects versus healthy controls was evaluated by area under the curve (AUC) analysis.
Pathway and network analyses
The differentially expressed gene lists from the meta-analysis and our bronchial thermoplasty dataset were compared. Pathway and network analyses were performed using MetaCore (Clarivate Analytics, Philadelphia, PA, USA). The significant common pathways were defined as pathways with FDR adjusted p-values <0.05. The network was built from the overlapping genes using the algorithm of shortest paths with two maximal steps.
Results
Table 1 summarises the age and sex of the study population in our bronchial thermoplasty dataset and the data source of the studies included in the meta-analysis. There were 12 subjects with a total of 25 bronchial thermoplasty sessions in our bronchial thermoplasty dataset. 50% of the subjects were males. The mean±sd age was 55±13 years. Among the 12 subjects, eight subjects were ex-smokers and four subjects were never-smokers. Also, two subjects had peripheral eosinophilia with counts >500 μL−1. Seven subjects received omalizumab at the time of the thermoplasty. The mean±sd exhaled nitric oxide fraction (FENO) was 26±19 ppb. One subject had FENO >50 ppb and six subjects had FENO <25 ppb. 11 of the 12 subjects completed all three bronchial thermoplasty sessions. Five subjects had three full sets of gene expression data, three subjects had two sets of gene expression data and the remaining four subjects had one set of gene expression data. The CONSORT flow diagram is shown in figure 2. In total, 133 severe asthma cases and 107 healthy controls were included in the meta-analysis datasets. From the meta-analysis datasets, 28 severe asthma cases and 42 healthy controls were reserved for the validation cohort.
In our longitudinal bronchial thermoplasty dataset, 116 genes demonstrated significant change in expression levels in response to bronchial thermoplasty sessions (FDR adjusted p-value <0.001) using a variance component score test (supplementary table S1).
For the meta-analysis, we identified four available public datasets with a total of 240 subjects that met our inclusion criteria. Among them, GSE43696 [19], GSE63142 [20] and GSE89809 [21] were used for the discovery dataset, while GSE64913 [22] was used for the validation dataset. 421 genes were differentially expressed in severe asthmatic subjects compared with healthy controls (FDR adjusted p-value <0.05). From the model obtained using the validation dataset, we determined that the 421 genes we identified correctly distinguished severe asthmatic subjects from healthy controls with an AUC value of 0.85 (95% CI 0.75–0.95) (supplementary table S2).
Comparing the longitudinal bronchial thermoplasty and meta-analysis gene lists, there were eight overlapping genes (genes associated with severe asthma and changed by bronchial thermoplasty): AMIGO2, CBX7, NR3C2, SETBP1, SHANK2, SNTB1, STXBP1 and ZNF853 (figure 3). All these overlapping genes were downregulated in severe asthmatic subjects compared with healthy controls in the differential expression meta-analysis. Two identified genes did not overlap between datasets, but were members of the same Rho GTPase family: RAC3 was differentially expressed in the meta-analysis and RHOBTB2 was differentially expressed in our bronchial thermoplasty dataset (table 2).
The differentially expressed, overlapping genes we identified from the meta-analysis and our longitudinal bronchial thermoplasty dataset were then subjected to downstream pathway and network analysis; 10 significant common pathways were identified as shown in table 3. Multiple pathways were found to be associated with neurophysiological processes and the cystic fibrosis transmembrane conductance regulator (CFTR). An additional pathway associated with gene expression changes from bronchial thermoplasty alone (SREBP1 and PDK2) is associated with a development role of interleukin (IL)-8 during angiogenesis (adjusted p-value 0.014). A network was built by using the shortest path with a maximum of two steps from the eight overlapping genes and Rho GTPase (figure 4). The key hub genes in this network were Rho GTPase family members RAC2 and RHOBTB2. Other key genes included AR (connects Rho GTPase and AMIGO2), PKD2 (connects Rho GTPase and STXBP1) and SGK1 (connects NR3C2 and SHANK2), which were differentially expressed in the meta-analysis (table 2).
Discussion
To the best of our knowledge, this is the first study exploring the potential mechanisms of bronchial thermoplasty using a transcriptomic approach to gene expression profiles in bronchial epithelial cells. Our method is robust and novel by combining the meta-analysis of public repositories and time-course longitudinal analysis of our patient cohort to identify the severe asthma-associated genes that were expressed differentially after serial bronchial thermoplasty treatments. Our methods utilise a longitudinal study design to test the trajectory of gene expression rather than simply comparing discrete time-points. By comparing our gene list with the results from the meta-analysis to identify overlapping genes, we focus on the possible disease-modifying mechanisms of bronchial thermoplasty rather than simply the overall pathways changed by bronchial thermoplasty. Our findings suggest that neurophysiological processes, CFTR activity and Rho GTPase are altered during bronchial thermoplasty sessions and represent possible important mechanisms in severe asthma.
Our study used airway epithelial cell samples instead of peripheral blood, which we suggest is a more robust representation of the changes that are occurring in the airway. We chose to sample airway cells that are untreated by the thermal energy of bronchial thermoplasty. The fact that we observed longitudinal gene expression changes over the course of bronchial thermoplasty sessions strongly suggests that we are observing fundamental changes to the entire lung that are induced by this therapy rather than a localised response to injury and repair. This finding is consistent with previous observations that bronchial thermoplasty influences important changes in adjacent untreated airways [10–12], which are likely to be, or reflect, the biological changes responsible for positive outcomes in severe asthma patients.
For example, NR3C2, identified in our study, encodes for a mineralocorticoid receptor and is bound by cortisol with high affinity. Single nucleotide polymorphisms residing in NR3C2 are found to be associated with asthma-related cytokines and chemokines, including interferon-γ, IL-13 and eotaxin [23]. Reduction in NR3C2 after bronchial thermoplasty is a potential mechanism for the observed decreased concentration of eosinophils in bronchoalveolar lavage of treated individuals [11]. Similarly, STXBP1 encodes for a syntaxin-binding protein and plays an important role in the release of neurotransmitters via regulation of syntaxin. This gene is involved in multiple neurophysiological processes as shown by our pathway analysis. Alteration of this gene after bronchial thermoplasty is consistent with the findings from previous studies that bronchial thermoplasty decreases the number of submucosal and airway smooth muscle-associated nerve fibres, as well as epithelial neuroendocrine cells [13]. Interestingly, both STXBP1 and SHANK2 have been found to independently regulate CFTR activity [24, 25]. Syntaxin 1A has been reported to regulate ion channels, possibly through direct inhibition of CFTR with direct protein–protein interaction [24]. SHANK2 physically and functionally is linked with cyclic nucleotide phosphodiesterase 4, an important target for novel therapies in asthma given its role in bronchoconstriction [25]. Although CFTR activity has been reported to be associated with chronic pulmonary diseases [26], its role in asthma is unclear. However, salt and water balance are critical components of airway calibre, and our findings suggest a novel pathway for future study.
Our network analysis demonstrates that the Rho GTPase family is a key hub connecting multiple gene networks. The Rho family of GTPases is a family of small signalling G proteins. The Rho GTPase pathway was found to impact calcium sensitisation of smooth muscle and regulate the contraction of the bronchial smooth muscle [27, 28]. Its downstream target, Rho kinase, has been investigated as a therapeutic target in the treatment of asthma [29, 30]. Based on our findings, bronchial thermoplasty may modify the function of members in the Rho GTPase family and therefore impact bronchial constriction in addition to reducing smooth muscle area.
As we obtained airway brush samples prior to each treatment session, our study is limited by the absence of an airway cell brush sample after subjects completed the final of the three bronchial thermoplasty sessions. The three samples we included for analysis were baseline (before any bronchial thermoplasty treatment), after the first session of bronchial thermoplasty (just prior to the second treatment) and after the second session of bronchial thermoplasty (just prior to the third treatment). Therefore, we are able to detect the processes altered from baseline through two treatment sessions, but not after completion of all three bronchial thermoplasty sessions.
Another limitation is the lack of robust clinical characterisation of the patients undergoing bronchial thermoplasty beyond that presented in table 1. Most of our patients were referrals from outside institutions. Therefore, we have no established baseline characteristics, particularly each patient's severe asthma phenotype (e.g. uncontrolled airway inflammation versus airway hyperresponsiveness). For similar reasons, we are unable to confirm a link between gene expression changes and related clinical outcomes in patients undergoing bronchial thermoplasty. Future studies linking changes in gene expression with specific clinical outcomes (i.e. changes in airway hyperresponsiveness) or specific severe asthma phenotypes (neutrophilic versus eosinophilic inflammation) are needed to better characterise those patients that will most likely benefit from bronchial thermoplasty.
Utilising a transcriptomic approach, we have identified potential mechanisms to explain the effectiveness of bronchial thermoplasty. We provided evidence that bronchial thermoplasty impacts neurophysiological processes, which is highly consistent with our current understanding of severe asthma pathobiology. We additionally identified novel targets, including CFTR and Rho GTPase family members. Future studies will need to focus on the relationship between these identified pathways and clinical response to therapy, and if differential expression of these pathways in severe asthmatic subjects identifies those patients most likely to respond to bronchial thermoplasty as a therapeutic option.
Supplementary material
Supplementary Material
Please note: supplementary material is not edited by the Editorial Office, and is uploaded as it has been supplied by the author.
Supplementary tables 1 and 2 00123-2018_supp_tables
Acknowledgements
We would like to thank Brian Morrissey, Michael Schivo, Christian Sebat, Nick Stollenwerk, Amir Zeki and Jill Steinbacher (UC Davis, Sacramento, CA, USA; B.M., M.S. and N.S. are also affiliated with VA Northern California Health Care System, Mather, CA, USA) for their help to collect airway specimens, and staff from the UC Davis Genomics Core for RNA sequencing.
Footnotes
This article has supplementary material available from openres.ersjournals.com.
Author Contributions: S-Y. Liao, A.L. Linderholm, N.J. Kenyon and R.W. Harper wrote the manuscript. S-Y. Liao analysed data. A.L. Linderholm, N.J. Kenyon and R.W. Harper designed the study and wrote the grant. A.L. Linderholm performed sample handling and RNA extraction. K.Y. Yoneda collected clinical data and specimens. N.J. Kenyon and R.W. Harper jointly supervised the research.
Conflict of interest: S-Y. Liao has no conflicts of interest.
Conflict of interest: A.L. Linderholm has no conflicts of interest.
Conflict of interest: K.Y. Yoneda has no conflicts of interest.
Conflict of interest: N.J. Kenyon has no conflicts of interest.
Conflict of interest: R.W. Harper has no conflicts of interest.
Support statement: This study was supported by an NIH Shared Instrumentation Grant 1S10OD010786-01 and a UC Davis Signature Research in Genomics Program Award. Funding information for this article has been deposited with the Crossref Funder Registry.
- Received July 31, 2018.
- Accepted December 19, 2018.
- Copyright ©ERS 2019
This article is open access and distributed under the terms of the Creative Commons Attribution Non-Commercial Licence 4.0.