Local mutational diversity drives intratumoral immune heterogeneity in non-small cell lung cancer

Abstract

Combining whole exome sequencing, transcriptome profiling, and T cell repertoire analysis, we investigate the spatial features of surgically-removed biopsies from multiple loci in tumor masses of 15 patients with non-small cell lung cancer (NSCLC). This revealed that the immune microenvironment has high spatial heterogeneity such that intratumoral regional variation is as large as inter-personal variation. While the local total mutational burden (TMB) is associated with local T-cell clonal expansion, local anti-tumor cytotoxicity does not directly correlate with neoantigen abundance. Together, these findings caution against that immunological signatures can be predicted solely from TMB or microenvironmental analysis from a single locus biopsy.

Fig. 1

Total mutational burden correlates with T-cell clonality. a Schematic of sampling strategy and experimental workflow. Tissues from multiple loci within the whole tumor were resected and subjected to high-throughput sequencing. Non-synonymous mutations, HLA typing, predicted neoantigens, transcriptomic profiling, and T-cell repertoire were analyzed to characterize TME heterogeneity. b–d Heat maps depicting the inter-population and intratumoral distribution of non-synonymous mutations, predicted neoantigens (with binding strength <500 nM), and dominant T-cell clones (frequencies ≥0.5%) in all sequenced subjects; presence (blue) or absence (white) is indicated for every tumor focus. Samples were grouped according to individual patients. e, f Scatterplot showing correlation between total mutation load and expanded properties of the T-cell repertoire. T-cell clonality (e) and Simpson diversity index (f) were used to depict the T-cell repertoire composition. Enrichment of highly expanded clones results in higher values for clonality and Simpson diversity. R coefficient of Pearson's correlation. Shaded area, 95% confident interval for the correlation. g–i Correlation between mutation load (in log2 scale) with expression of interferon-gamma, granzyme-A, and cytolytic activity (measured as the geometric mean of granzyme-A with perforin 1) in log2 of transcript per kilobase million (TPM)

Fig. 2

Machine learning classifies TME into hot and cold immunophenotypes. a Visualization based on two-dimensional coordinates from multidimensional scaling (MDS) of proximity matrix from the input variables in NSCLC. Color indicates cytolytic activity (product of PRF1 and GZMA) for each sampling site. To categorize the samples by an unsupervised method, the Gaussian expectation maximization algorithm was employed to perform categorization under Gaussian mixture models. The contour shows the estimated probability density for the two categories. Left, cold area; right, hot area. b, c Expression profiling of prognostic genes in the POPLAR study for hot and cold areas. Expression values are transformed TPM format and in log2 scale. Horizontal bar in boxplot, median value. Statistics based on two-tailed Mann–Whitney U-test. d Total mutational burden in the categorization. e Samples from each individual are labeled by color. Contours guiding immunologic categorization are shown as in a. f Upper panel, pie chart showing the proportion of hot vs cold area samples in all sequenced tissues. Lower panel, proportion of patients harboring multiple distinct immunological statuses simultaneously. Homogeneous, patients with one kind of immunological status; Heterogeneous, patients with both immunological statuses. g Correlation among immunological status, TMB, and PD-L1 expression. Upper panel, the proportion of immunological statuses found in each individual. Middle and lower panels, TMB and PD-L1 expression for each sample. For each individual, the samples in the middle and lower panels were labeled in the same colors. h Heatmap showing the expression of feasible immunotherapy-predictive immune genes for 44 NSCLC samples. Expression values for each gene are normalized into z-scores

Fig. 3

Correlation of TMB and immune cell infiltration heterogeneity. a Single-sample gene set enrichment analysis identifying the relative infiltration of immune cell populations for 44 NSCLC tumor samples with available RNA-sequencing data. The relative infiltration of each cell type is normalized into a z-score. b Validation by flow cytometry. Samples were collected separately from the sequenced batch. c Correlation between infiltration of cell types executing anti-tumor immunity (ActCD4, ActCD8, TcmCD4, TcmCD8, TemCD4, TemCD8, Th1, Th17, ActDC, CD56briNK, NK, NKT) and cell types executing pro-tumor, immune suppressive functions (Treg, Th2, CD56dimNK, imDC, TAM, MDSC, Neutrophil, and pDC). R coefficient of Pearson's correlation. The shaded area represents 95% confident interval. d Scatterplot showing correlation between heterogeneities of total mutation burden and immune cell infiltration. An average pairwise correlation coefficient was calculated to quantify the divergence of immune cell infiltration, and presented as 1-value for clear visualization. Higher values, heterogeneous immune cell infiltration; lower values, homogeneous immune cell infiltration. The variation in genomic mutations was determined by either intratumoral heterogeneity (ITH) or coefficient of variance (CV). Higher values, heterogeneous; lower values, homogeneous mutation pattern. R coefficient of Pearson's correlation. e Phylogenetic trees generated by a parsimony ratchet approach based on the distribution of all detected mutations are shown for patients with heterogeneous immunological status; trunk and branch lengths are proportional to the number of non-synonymous mutations acquired. Red label, samples in hot area; blue label, samples in cold area