Targeted lung denervation modulates the mucosal epithelial transcriptome in COPD

To the Editor:

Targeted lung denervation (TLD) is a novel treatment for COPD with evidence of a potential beneficial effect on exacerbation frequency [1]. It is a minimally invasive procedure involving the ablation of the airway nerve trunks of the vagus nerve that travel parallel to and outside the main bronchi, thereby reducing neuronal acetylcholine release. It has been postulated that TLD may share a similar mechanism of action with the inhaled long-acting muscarinic antagonists [2].

In a double-blind, randomised, sham-controlled study (AIRFLOW-2; clinicaltrials.gov identifier NCT02058459), a significant reduction (p<0.001) in respiratory adverse events, which included COPD exacerbations, at 3–6.5 months post-procedure was reported in the TLD arm. Furthermore, over 12.5 months of follow-up, the risk of a severe COPD exacerbation requiring hospitalisation, as assessed via time-to-first event analysis, was lower for the treatment group (p=0.039) [1].

We conducted a substudy to evaluate the post-treatment airway mucosal transcriptome using next-generation RNA sequencing (seq) of mucosal brush samples. The aim of this study was to explore gene expression changes between the post-treatment TLD and sham-control patients and provide hypotheses for future investigation into the mechanisms of TLD.

Airway mucosal brush samples collected at the 3-month follow-up visit in both sham-control group and treatment group were used for this study. For each patient, three brushes were collected from the right lower lobe, which constituted one sample. Samples were processed using methods as described previously [3, 4]. Data analysis was performed on a short-read dataset obtained using Illumina next-generation sequencing technology. RNA-seq was conducted using the Illumina NovaSeq 6000 sequencer by GenomeScan (www.genomescan.nl/). The procedure included data quality control, adapter trimming, alignment of short reads and feature counting. Library preparation was checked by calculating ribosomal (and globin) content. Checks for possible sample and barcode contamination were performed and a set of standard quality metrics for the raw dataset was determined using quality control tools (FstQC v0.34 and FastQA). Prior to alignment, the reads were trimmed for adapter sequences using Trimmomatic v0.30. To align the reads of each sample, the human reference GRCh37.75 was used. 25 samples (13 and 12 in the TLD and sham-control arms, respectively) passed quality control and were used in the analysis.

Differential gene analysis was performed using DESeq2 package v1.24.0 on the R platform v3.6.0. The full gene set was filtered for low read counts by excluding all genes with an average fragment per kilobase million (FPKM) of <1. A t-test was then performed for the filtered read counts for each gene (TLD versus sham-control); genes were ranked by nominal p-value and the corresponding q-values were calculated to correct for multiple testing. Unfortunately, no genes met the 0.25 false discovery rate (FDR) criterion for transcriptome-wide significance. As the main objective of this initial analysis was to generate hypotheses for future investigation, genes with a nominal p-value of <0.05 were included for the hypergeometric distribution overrepresentation analysis.

Gene set enrichment analysis uses the whole unfiltered gene set and ranks it according to expression levels. It then uses a database of gene sets to find biologically related pathways that are significantly expressed in the up- or downregulated ends of the ranked gene list. We used the Reactome database (https://reactome.org/), of which four gene sets were significantly downregulated in the TLD arm. The top two downregulated gene sets were related to acetylcholine, the main neurotransmitter involved in the parasympathetic nervous system: “Reactome highly calcium permeable postsynaptic nicotinic acetylcholine receptors” (p=0.008, q=0.047) and “Reactome acetylcholine binding and downstream events” (p=0.03, q=0.04).

In order to perform the HGD analysis, the FPKM >1 genes were separated into upregulated and downregulated genes. Of the upregulated genes, 40 had a nominal p<0.05, but no Reactome gene sets were found to be overrepresented on HGD analysis. Of the downregulated genes, 991 had a nominal p<0.05 and several immunity-related gene sets were overrepresented. On closer inspection, these immunity-related gene sets shared a number of common genes. These 35 genes were fed into the StringDB website (https://string-db.org/), which performs cluster analysis based on known protein–protein interactions. A large cluster of genes was noted to relate to the ubiquitin-proteasome system, which is the cellular machinery used to dispose of misfolded proteins (figure 1).

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 1

a) Mean average plot showing only modest log fold changes; high variance across the range of expression levels; and a lack of any significantly differentially expressed genes. The low number of biological replicates (i.e. participants sampled) may account for this. b) Principal component (PC) analysis plot showing horizontal clustering of the sham cases to the left with PC1 accounting for over half of the total variance within the sequencing read count dataset. This is suggestive of a significant difference in gene expression between the two groups. c) StringDB output for HGD-overrepresented downregulated genes. The large cluster of genes on the left relate to the ubiquitin-proteasome system.

In the present study, whole-transcriptome sequencing of mucosal brush samples has shown reduced expression of acetylcholine-related genes after TLD. This suggests reductions in mucosal acetylcholine pathways, which may have been indirectly induced by TLD itself. These acetylcholine-related gene sets refer to nicotinic receptors, which are known to be expressed in airway epithelial cells and may even correlate with the development of airway obstruction in COPD [5, 6]. No genes were significantly differentially expressed after correction for multiple testing. However, when genes with nominal p<0.05 were analysed with HGD and StringDB, we identified a cluster of genes common to several overrepresented immunity-related gene sets that relate to the ubiquitin-proteasome system (UPS).

The UPS removes denatured, misfolded, damaged or improperly translated proteins from cells. It is a highly complex pathway that is involved in many important cellular processes, including the regulation of immune and inflammatory responses as well as the cellular response to stress [7]. The UPS is known to regulate CD8⁺ T-cell response to viral infections, activate the NF-кB pathway and manage oxidative stress [8]. It is noted that all three processes are relevant to the pathophysiology of COPD and are increased during periods of exacerbation [9, 10]. Given that COPD exacerbations were reduced in the TLD arm of the AIRFLOW-2 study, we postulate that the cluster of downregulated UPS-related genes is a surrogate marker of reduced airway inflammation and oxidative stress, suggestive of a potential anti-inflammatory mechanism of TLD. Baseline data (table 1) show that there was a tendency toward more severe airways obstruction and worse health-related quality of life in the treatment group. One would expect airways inflammation to be worse in this group with upregulation of immunity-related gene expression. However, our results show the converse, and would support the hypothesis that TLD has suppressed inflammatory pathways. As mentioned previously, there were no significantly differentially expressed genes after the FDR correction was applied. This may be due to low expression levels of genes potentially modifiable by TLD that would require a higher number of biological replicates to attain statistical significance. Correcting for multiple testing is an important step in analysing gene expression data in order to minimise the number of false positive results. By using nominal p-values (<0.05) instead of FDR q-values as our threshold for differential expression, the risk of Type 1 error is amplified. However, we feel this is an acceptable adjustment to make given the aim of this study is to merely generate hypotheses for future TLD studies. Another major limitation was the lack of baseline data and therefore no between-group analyses. This makes it more difficult to interpret the data as a true treatment effect and leaves more room for confounding factors.

Cellular material in bronchial brush samples is unlikely to be purely epithelial in origin. For example, there is likely to be a significant inflammatory cell component within these samples. Bulk RNA-seq methods, such as that used in this study, are unable to resolve these differences in cell origin. Therefore, it must be noted that the gene expression results presented here are likely to have been influenced, in part, by nonepithelial cells. This would act to dilute any signal induced by TLD, which is more likely to modulate gene expression in the epithelium than in luminal inflammatory cells. This issue may be mitigated in future studies by using single-cell RNA-seq, a powerful tool that can assess gene expression at the level of individual cells instead of a global average across all cells in any given sample.

In conclusion, we observed trends in reduced acetylcholine-related gene expression after targeted lung denervation, potentially serving as novel indirect evidence that to some degree, denervation has occurred in these patients. Trends were also seen in downregulation of genes related to immunity and inflammation, specifically the UPS. This may reflect reduced airway inflammation and oxidative stress after TLD, and perhaps explain the mechanism behind its effect on reducing exacerbations in COPD.