Molecular subtyping reveals immune alterations associated with progression of bronchial premalignant lesions

Jennifer E Beane, Sarah A Mazzilli, Joshua D Campbell, Grant Duclos, Kostyantyn Krysan, Christopher Moy, Catalina Perdomo, Michael Schaffer, Gang Liu, Sherry Zhang, Hanqiao Liu, Jessica Vick, Samjot S Dhillon, Suso J Platero, Steven M Dubinett, Christopher Stevenson, Mary E Reid, Marc E Lenburg, Avrum E Spira, Jennifer E Beane, Sarah A Mazzilli, Joshua D Campbell, Grant Duclos, Kostyantyn Krysan, Christopher Moy, Catalina Perdomo, Michael Schaffer, Gang Liu, Sherry Zhang, Hanqiao Liu, Jessica Vick, Samjot S Dhillon, Suso J Platero, Steven M Dubinett, Christopher Stevenson, Mary E Reid, Marc E Lenburg, Avrum E Spira

Abstract

Bronchial premalignant lesions (PMLs) are precursors of lung squamous cell carcinoma, but have variable outcome, and we lack tools to identify and treat PMLs at risk for progression to cancer. Here we report the identification of four molecular subtypes of PMLs with distinct differences in epithelial and immune processes based on RNA-Seq profiling of endobronchial biopsies from high-risk smokers. The Proliferative subtype is enriched with bronchial dysplasia and exhibits up-regulation of metabolic and cell cycle pathways. A Proliferative subtype-associated gene signature identifies subjects with Proliferative PMLs from normal-appearing uninvolved large airway brushings with high specificity. In progressive/persistent Proliferative lesions expression of interferon signaling and antigen processing/presentation pathways decrease and immunofluorescence indicates a depletion of innate and adaptive immune cells compared with regressive lesions. Molecular biomarkers measured in PMLs or the uninvolved airway can enhance histopathological grading and suggest immunoprevention strategies for intercepting the progression of PMLs to lung cancer.

Conflict of interest statement

JEB, SAM, JDC, CP, JV, MER, MEL, and AES received commercial research grants from Janssen Pharmaceuticals. MEL is a consultant to Veracyte. AES, MS, CM, and CS are employed by Johnson and Johnson. SJP was employed by Johnson and Johnson during the study and is now employed by Covance. GD, KK, GL, SZ, HL, SSD, and SMD declare no competing interests.

Figures

Fig. 1
Fig. 1
Endobronchial biopsies divide into four distinct molecular subtypes that correlate with clinical and molecular phenotypes. a Genes (n = 3936) organized into nine gene co-expression modules were used to discover four molecular subtypes (Proliferative, Inflammatory, Secretory, and Normal-like) across the 190 DC biopsies using consensus clustering. The heatmap shows semi-supervised hierarchal clustering of z-score-normalized gene expression across the 3936 genes and 190 DC biopsies. The top color bar represents the four molecular subtypes: Proliferative (n = 52 samples), Inflammatory (n = 37 samples), Secretory (n = 61 samples), and Normal-like (n = 40 samples). To the left of the heatmap, barplots for each module show the mean module GSVA score for each subtype. To the right of the heatmap, a summary of enriched biological pathways is listed for each module. b Bubbleplots showing significant associations (p < 0.01, two-sided Fisher’s exact test) between the molecular subtypes and genomic smoking status, biopsy histological grade, and the predicted LUSC tumor molecular subtypes. The columns represent the four molecular subtypes (Proliferative, Inflammatory, Secretory, and Normal-like) and the diameter of the circle is proportional to the number of samples within each subtype that have the row phenotype. c Boxplot of MKI67 expression values in biopsies with normal or hyperplasia histology (n = 8, 16, 26, 18 in Proliferative, Inflammatory, Secretory, and Normal-like subtypes, respectively). The MKI67 expression levels of the Proliferative subtype are significantly greater than non-Proliferative subtype samples (FDR = 3.4e-10, linear model). d Boxplot of expression values of MKI67 in biopsies with dysplastic histology (n = 33, 11, 19, 9 in Proliferative, Inflammatory, Secretory, and Normal-like subtypes, respectively). The MKI67 expression levels of the Proliferative subtype are significantly greater than non-Proliferative subtype samples (FDR = 3.1e-8). e Immunofluorescent staining demonstrating the increased MKI67 and KRT5 staining and reduced TUB1A1 staining in the Proliferative subtype. The representative samples shown for the Proliferative and Inflammatory subtypes have dysplasia histology, whereas the samples shown for the Secretory and Normal-like subtypes have normal histology (Magnification ×  200). In the boxplots, the upper and lower hinges correspond to the first and third quartile, center line represents the median, and whiskers extend from the hinge to the largest or smallest value at most 1.5 times the distance between the quartiles. Source data are provided as a Source Data file
Fig. 2
Fig. 2
Phenotypic associations with the molecular subtypes are confirmed in an independent sample set. a The 190 DC biopsies and the 3936 genes were used to build a 22-gene nearest centroid molecular subtype classifier. The heatmap shows semi-supervised hierarchal clustering of z-score normalized gene expression across the 22 classifier genes and 190 DC biopsies training samples. b The 22-gene nearest centroid molecular subtype classifier was used to predict the molecular subtypes of the 105 VC biopsies. The heatmap shows semi-supervised hierarchal clustering of z-score normalized gene expression across 22 genes and 105 VC. The rows of the heatmap give the gene name and module membership, and the column color bar shows molecular subtype membership. c Bubbleplots showing significant associations (p < 0.01 by two-sided Fisher’s exact test) between the VC molecular subtypes and smoking status, biopsy histological grade, and the predicted LUSC tumor molecular subtypes. The columns represent the four molecular subtypes (Proliferative, Inflammatory, Secretory, and Normal-like) and the radius of the circle is proportional to the number of samples within each subtype that have the row phenotype. Source data are provided as a Source Data file
Fig. 3
Fig. 3
Normal-appearing bronchial brushes predict the presence of proliferative lesions. a The DC (left) and VC (right) cohorts, showing the number of brushes (y axis) classified as proliferative orange) that have at least one biopsy (y axis) classified as proliferative at the time the brush was sampled. Brushes/biopsies classified as not proliferative are turquoise. b Boxplots of GSVA scores for modules 4, 5, 6, and 7 (y axis) across all brushes (n = 86 in DC and n = 48 in VC) and biopsies (n = 190 in DC and n = 105 in VC) from each cohort classified as Proliferative or not Proliferative (x axis). The red asterisk indicates significant differences between the Proliferative subtype versus all other samples (FDR < 0.05, linear model). In the boxplots, the upper and lower hinges correspond to the first and third quartile, center line represents the median, and whiskers extend from the hinge to the largest or smallest value at most 1.5 times the distance between the quartiles. Source data are provided as a Source Data file
Fig. 4
Fig. 4
Immune alterations are associated with lesion outcome in the Proliferative subtype. Boxplots of Module 9 GSVA scores across DC a and VC biopsies b within the Proliferative subtype. There is a significant difference between the progressive/persistent versus regressive biopsies (p = 0.002 (DC) and p = 0.03 (VC), linear models). c Top: heatmap of z-score-normalized gene expression across the 112 genes in Module 9 in the DC biopsies (left) and the VC biopsies (right). Each heatmap is supervised by Module 9 GSVA scores. Top color bars indicate the histological grade of the biopsies and their progression status. Bottom: heatmap of xCell results indicating the relative abundance of immune cell types across the DC biopsies (left) and the VC biopsies (right). Immune cell types displayed are significantly associated with lesion progression/persistence (FDR < 0.05 in both the DC and VC, linear model). d Representative histology where the dashed yellow line denotes the separation of epithelium and stromal compartments. Top panels: a progressive severe dysplasia has reduced presence of immune cells demonstrated by the marked reduction in expression of M2 macrophages (CD68/163 staining, double-positive cells indicated by the yellow arrows) and CD8 T cells (sample corresponds to *P in c). Bottom panels: a regressive moderate dysplasia has increased presence of immune cells including M2 macrophages (CD68/163 staining double-positive cells indicated by the yellow arrows) and CD8 T cells (samples correspond to *R in c). e Boxplots of the percentages of CD68 and CD163, CD68, CD163, CD4, and CD8 positively stained cells between progressive/persistent and regressive biopsies (p < 0.001, linear model, for all comparisons). The x axis labels indicate the number of regions (R) enumerated across (P) subjects for each stain and outcome group depicted in the boxplot. Biopsies were included in the analysis if their clinical outcome was concordant with the Module 9 score. In the boxplots, the upper and lower hinges correspond to the first and third quartile, center line represents the median, and whiskers extend from the hinge to the largest or smallest value at most 1.5 times the distance between the quartiles. Source data are provided as a Source Data file

References

    1. Campbell JD, et al. The case for a Pre-Cancer Genome Atlas (PCGA) Cancer Prev. Res. (Philos.) 2016;9:119–124. doi: 10.1158/1940-6207.CAPR-16-0024.
    1. Beane J, Campbell JD, Lel J, Vick J, Spira A. Genomic approaches to accelerate cancer interception. Lancet Oncol. 2017;18:e494–e502. doi: 10.1016/S1470-2045(17)30373-X.
    1. Auerbach O, Stout AP, Hammond EC, Garfinkel L. Changes in bronchial epithelium in relation to cigarette smoking and in relation to lung cancer. N. Engl. J. Med. 1961;265:253–267. doi: 10.1056/NEJM196108102650601.
    1. van Boerdonk RAA, et al. Close surveillance with long-term follow-up of subjects with preinvasive endobronchial lesions. Am. J. Respir. Crit. Care Med. 2015;192:1483–1489. doi: 10.1164/rccm.201504-0822OC.
    1. Merrick DT, et al. Persistence of bronchial dysplasia is associated with development of invasive squamous cell carcinoma. Cancer Prev. Res. 2016;9:96–104. doi: 10.1158/1940-6207.CAPR-15-0305.
    1. Ishizumi T, McWilliams A, MacAulay C, Gazdar A, Lam S. Natural history of bronchial preinvasive lesions. Cancer Metastas-. Rev. 2010;29:5–14. doi: 10.1007/s10555-010-9214-7.
    1. Beane J, et al. Reversible and permanent effects of tobacco smoke exposure on airway epithelial gene expression. Genome Biol. 2007;8:R201. doi: 10.1186/gb-2007-8-9-r201.
    1. Wang L, Wang S, Li W. RSeQC: quality control of RNA-seq experiments. Bioinformatics. 2012;28:2184–2185. doi: 10.1093/bioinformatics/bts356.
    1. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559. doi: 10.1186/1471-2105-9-559.
    1. Campbell JD, et al. Distinct patterns of somatic genome alterations in lung adenocarcinomas and squamous cell carcinomas. Nat. Genet. 2016;48:607–616. doi: 10.1038/ng.3564.
    1. Cancer Genome Atlas Research Network. Comprehensive genomic characterization of squamous cell lung cancers. Nature. 2012;489:519–525. doi: 10.1038/nature11404.
    1. Wilkerson MD, et al. Lung squamous cell carcinoma mRNA expression subtypes are reproducible, clinically important, and correspond to normal cell types. Clin. Cancer Res. 2010;16:4864–4875. doi: 10.1158/1078-0432.CCR-10-0199.
    1. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7. doi: 10.1186/1471-2105-14-7.
    1. Beane Jennifer, Mazzilli Sarah A., Tassinari Anna M., Liu Gang, Zhang Xiaohui, Liu Hanqiao, Buncio Anne Dy, Dhillon Samjot S., Platero Suso J., Lenburg Marc E., Reid Mary E., Lam Stephen, Spira Avrum E. Detecting the Presence and Progression of Premalignant Lung Lesions via Airway Gene Expression. Clinical Cancer Research. 2017;23(17):5091–5100. doi: 10.1158/1078-0432.CCR-16-2540.
    1. Bindea G, et al. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity. 2013;39:782–795. doi: 10.1016/j.immuni.2013.10.003.
    1. Aran D, Hu Z, Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. 2017;18:220. doi: 10.1186/s13059-017-1349-1.
    1. Merrick DT, et al. Analysis of c-ErbB1/epidermal growth factor receptor and c-ErbB2/HER-2 expression in bronchial dysplasia: evaluation of potential targets for chemoprevention of lung cancer. Clin. Cancer Res. 2006;12:2281–2288. doi: 10.1158/1078-0432.CCR-05-2291.
    1. Jeanmart M, et al. Value of immunohistochemical markers in preinvasive bronchial lesions in risk assessment of lung cancer. Clin. Cancer Res. 2003;9:2195–2203.
    1. Lantuejoul S, et al. Telomere maintenance and DNA damage responses during lung carcinogenesis. Clin. Cancer Res. 2010;16:2979–2988. doi: 10.1158/1078-0432.CCR-10-0142.
    1. Wistuba II, et al. Sequential molecular abnormalities are involved in the multistage development of squamous cell lung carcinoma. Oncogene. 1999;18:643–650. doi: 10.1038/sj.onc.1202349.
    1. Wistuba II, et al. High resolution chromosome 3p allelotyping of human lung cancer and preneoplastic/preinvasive bronchial epithelium reveals multiple, discontinuous sites of 3p allele loss and three regions of frequent breakpoints. Cancer Res. 2000;60:1949–1960.
    1. Nakachi I, et al. Application of SNP microarrays to the genome-wide analysis of chromosomal instability in premalignant airway lesions. Cancer Prev. Res. (Philos.) 2014;7:255–265. doi: 10.1158/1940-6207.CAPR-12-0485.
    1. Massion PP, et al. Recurrent genomic gains in preinvasive lesions as a biomarker of risk for lung cancer. PLoS ONE. 2009;4:e5611. doi: 10.1371/journal.pone.0005611.
    1. Gower AC, Spira A, Lenburg ME. Discovering biological connections between experimental conditions based on common patterns of differential gene expression. BMC Bioinformatics. 2011;12:381. doi: 10.1186/1471-2105-12-381.
    1. Rahman SMJ, et al. The airway epithelium undergoes metabolic reprogramming in individuals at high risk for lung cancer. JCI Insight. 2016;1:e88814. doi: 10.1172/jci.insight.88814.
    1. Keith RL, et al. Oral iloprost improves endobronchial dysplasia in former smokers. Cancer Prev. Res. 2011;4:793–802. doi: 10.1158/1940-6207.CAPR-11-0057.
    1. Lam S, et al. A randomized phase IIb trial of myo-inositol in smokers with bronchial dysplasia. Cancer Prev. Res. 2016;9:906–914. doi: 10.1158/1940-6207.CAPR-15-0254.
    1. Ridker PM, et al. Effect of interleukin-1β inhibition with canakinumab on incident lung cancer in patients with atherosclerosis: exploratory results from a randomised, double-blind, placebo-controlled trial. Lancet. 2017;390:1833–1842. doi: 10.1016/S0140-6736(17)32247-X.
    1. Spira A, et al. Effects of cigarette smoke on the human airway epithelial cell transcriptome. Proc. Natl. Acad. Sci. USA. 2004;101:10143–10148. doi: 10.1073/pnas.0401422101.
    1. Steiling K, et al. A dynamic bronchial airway gene expression signature of chronic obstructive pulmonary disease and lung function impairment. Am. J. Respir. Crit. Care Med. 2013;187:933–942. doi: 10.1164/rccm.201208-1449OC.
    1. Spira A, et al. Airway epithelial gene expression in the diagnostic evaluation of smokers with suspect lung cancer. Nat. Med. 2007;13:361–366. doi: 10.1038/nm1556.
    1. Whitney DH, et al. Derivation of a bronchial genomic classifier for lung cancer in a prospective study of patients undergoing diagnostic bronchoscopy. BMC Med. Genomics. 2015;8:18. doi: 10.1186/s12920-015-0091-3.
    1. Silvestri Gerard A., Vachani Anil, Whitney Duncan, Elashoff Michael, Porta Smith Kate, Ferguson J. Scott, Parsons Ed, Mitra Nandita, Brody Jerome, Lenburg Marc E., Spira Avrum. A Bronchial Genomic Classifier for the Diagnostic Evaluation of Lung Cancer. New England Journal of Medicine. 2015;373(3):243–251. doi: 10.1056/NEJMoa1504601.
    1. Beane J, et al. A prediction model for lung cancer diagnosis that integrates genomic and clinical features. Cancer Prev. Res. (Philos.) 2008;1:56–64. doi: 10.1158/1940-6207.CAPR-08-0011.
    1. Wang S, et al. Expression of CD163, interleukin-10, and interferon-gamma in oral squamous cell carcinoma: mutual relationships and prognostic implications. Eur. J. Oral Sci. 2014;122:202–209. doi: 10.1111/eos.12131.
    1. Gettinger S, et al. Impaired HLA class I antigen processing and presentation as a mechanism of acquired resistance to immune checkpoint inhibitors in lung cancer. Cancer Discov. 2017;7:1420–1435. doi: 10.1158/-17-0593.
    1. Pereira C, et al. Genomic profiling of patient-derived xenografts for lung cancer identifies B2M inactivation impairingimmunorecognition. Clin. Cancer Res. 2017;23:3203–3213. doi: 10.1158/1078-0432.CCR-16-1946-T.
    1. Gao J, et al. Loss of IFN-γ pathway genes in tumor cells as a mechanism of resistance to anti-CTLA-4 therapy. Cell. 2016;167:397–404.e9. doi: 10.1016/j.cell.2016.08.069.
    1. Kotsakis, A. et al. Prognostic value of circulating regulatory T cell subsets in untreated non-small cell lung cancer patients, Sci. Rep.6, 39247 (2016).
    1. Wu Si-Pei, Liao Ri-Qiang, Tu Hai-Yan, Wang Wen-Jun, Dong Zhong-Yi, Huang Shu-Mei, Guo Wei-Bang, Gou Lan-Ying, Sun Hui-Wen, Zhang Qi, Xie Zhi, Yan Li-Xu, Su Jian, Yang Jin-Ji, Zhong Wen-Zhao, Zhang Xu-Chao, Wu Yi-Long. Stromal PD-L1–Positive Regulatory T cells and PD-1–Positive CD8-Positive T cells Define the Response of Different Subsets of Non–Small Cell Lung Cancer to PD-1/PD-L1 Blockade Immunotherapy. Journal of Thoracic Oncology. 2018;13(4):521–532. doi: 10.1016/j.jtho.2017.11.132.
    1. International Agency for Research on Cancer. Who classification of tumours of the lung, pleura, thymus and heart. (World Health Organization, 2015).
    1. Dobin A, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21. doi: 10.1093/bioinformatics/bts635.
    1. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323. doi: 10.1186/1471-2105-12-323.
    1. Pedersen BS, Quinlan AR. Who’s who? Detecting and resolving sample anomalies in human DNA sequencing studies with peddy. Am. J. Hum. Genet. 2017;100:406–413. doi: 10.1016/j.ajhg.2017.01.017.
    1. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–140. doi: 10.1093/bioinformatics/btp616.
    1. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010;26:1572–1573. doi: 10.1093/bioinformatics/btq170.
    1. Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. 2012;28:882–883. doi: 10.1093/bioinformatics/bts034.
    1. Chen EY, et al. Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics. 2013;14:128. doi: 10.1186/1471-2105-14-128.
    1. Subramanian A, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. USA. 2005;102:15545–15550. doi: 10.1073/pnas.0506580102.
    1. Liberzon A, et al. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 2015;1:417–425. doi: 10.1016/j.cels.2015.12.004.
    1. Ritchie ME, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47. doi: 10.1093/nar/gkv007.
    1. Law CW, Chen Y, Shi W, Smyth GK. Voom: precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biol. 2014;15:R29. doi: 10.1186/gb-2014-15-2-r29.
    1. Yoshihara K, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 2013;4:2612. doi: 10.1038/ncomms3612.
    1. Dvorak, A., Tilley, A. E., Shaykhiev, R., Wang, R. & Crystal, R. G. Do airway epithelium air-liquid cultures represent the in vivo airway epithelium transcriptome? Am. J. Respir. Cell Mol. Biol. 44, 465–473 (2011).

Source: PubMed

3
Subscribe