Airway transcriptomic profiling after bronchial thermoplasty.

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.

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.

FIGURE 1

Schematic of the analysis workflow. DEG: differentially expressed genes; BT: bronchial thermoplasty.

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.

TABLE 1

Characteristics of study subjects

	Study
	BT dataset	GSE43696	GSE63142	GSE89809	GSE64913
Use	Longitudinal gene expression change	Meta-analysis (discovery)	Meta-analysis (discovery)	Meta-analysis (discovery)	Meta-analysis (validation)
PubMed ID [ref.]		24518246 [19]	25338189 [20]	28933920 [21]	28045928 [22]
Data type	RNA sequencing	Microarray	Microarray	Microarray	Microarray
Severe asthmatic subjects n	12	38	56	11	28
Males in severe asthmatic subjects %	50.00	26.32	NP	45.45	57.14
Mean±sd age in severe asthmatic subjects years	55.17±13.07	44.40±10.06	NP	53.00±12.04	40.57±13.96
Median (interquartile range) FEV1 % pred#	59 (54–68)	56 (40–71)	56±20 (mean±sem)	63 (52–75)	64 (26–123) (mean (range))
Maintenance oral corticosteroids# %	41.67	NP	NP	28.57	23.53
Healthy controls n	0	20	27	18	42

BT: bronchial thermoplasty; FEV1: forced expiratory volume in 1 s; NP: not provided in the data. #: value extracted from the original paper, not the published RNA microarray data.

FIGURE 2

CONSORT flow diagram. BT: bronchial thermoplasty.

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).

FIGURE 3

Number of unique and common genes in each analysis. BT: bronchial thermoplasty.

TABLE 2

Differentially expressed genes from meta-analyses and the clinical cohort of patients undergoing bronchial thermoplasty (BT)

Gene symbol	Description	Differential expression meta-analysis		Genes changed over BT sessions: FDR adjusted p-value¶
		Effect size#	FDR adjusted p-value¶
Overlapping genes
SNTB1	β1-syntrophin	−1.07	1.55×10−17	<0.001
CBX7	Chromobox protein homologue 7	−0.82	1.13×10−8	<0.001
ZNF853	Zinc finger protein 853	−0.72	1.56×10−7	<0.001
NR3C2	Mineralocorticoid receptor	−0.97	1.09×10−6	<0.001
STXBP1	Syntaxin-binding protein 1	−0.70	1.54×10−6	<0.001
SETBP1	SET-binding protein	−0.64	3.47×10−6	<0.001
AMIGO2	Amphoterin-induced protein 2	−0.92	3.84×10−6	<0.001
SHANK2	SH3 and multiple ankyrin repeat domains protein 2	−0.46	1.36×10−4	<0.001
Rho GTPase
RAC3	Ras-related C3 botulinum toxin substrate 3	0.85	2.40×10−5	0.247
RHOBTB2	Rho-related BTB domain-containing protein 2	0.17	2.64×10−1	<0.001
Unique genes in the network
AR	Androgen receptor	−0.81	3.23×10−7	0.090
PKD2	Protein kinase D2	−0.83	7.47×10−9	0.157
SGK1	Serum/glucocorticoid regulated kinase 1	0.92	4.80×10−6	1.000

FDR: false discovery rate. #: standardised mean difference (log2 scale); ¶: p-value adjusted by the Benjamini–Hochberg procedure.

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).

TABLE 3

Significant pathways (false discovery rate (FDR) adjusted p-value <0.05) associated with overlapping genes between differentially expressed genes from meta-analyses and genes that changed over serial bronchial thermoplasty sessions

MetaCore pathway name	Genes n	p-value	FDR adjusted p-value#	Genes from data
Regulation of CFTR activity	62	3.601×10−4	3.601×10−3	STXBP1, SHANK2
Cytoskeleton remodelling neurofilaments	25	1.221×10−2	2.438×10−2	STXBP1
Neurophysiological process role of CDK5 in pre-synaptic signalling	28	1.367×10−2	2.438×10−2	STXBP1
Neurophysiological process GABA-B receptor signalling in pre-synaptic nerve terminals	28	1.367×10−2	2.438×10−2	STXBP1
Muscle contraction nNOS signalling in skeletal muscle	28	1.367×10−2	2.438×10−2	SNTB1
Transport aldosterone-mediated regulation of ENaC sodium transport	30	1.464×10−2	2.438×10−2	NR3C2
Neurophysiological process Ephrin-B receptors in dendritic spine morphogenesis and synaptogenesis	35	1.706×10−2	2.438×10−2	SHANK2
wtCFTR and ΔF508 traffic/membrane expression	47	2.286×10−2	2.807×10−2	STXBP1
Neurophysiological process synaptic vesicle fusion and recycling in nerve terminals	52	2.527×10−2	2.807×10−2	STXBP1
Neurophysiological process constitutive and activity-dependent synaptic AMPA receptor delivery	59	2.863×10−2	2.863×10−2	SHANK2

CFTR: cystic fibrosis transmembrane conductance regulator; CDK: cyclin-dependent kinase; GABA: γ-aminobutyric acid; nNOS: neuronal nitric oxide synthase; ENaC: epithelial sodium channel; wt: wild-type; AMPA: α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid. #: p-value adjusted by the Benjamini–Hochberg procedure.

FIGURE 4

Network built from overlapping genes and Rho GTPase (SNTB1 and ZNF853 were not found to interact with other genes in the network), modified from the MetaCore network. Circled: overlapping genes; green arrow: activation; red arrow: inhibition; grey arrow: unspecified.

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.

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