Head and Neck Squamous Cell Carcinomas Are Characterized by a Stable Immune Signature Within the Primary Tumor Over Time and Space

Morphologic heterogeneity in cancer has long been appreciated by pathologists; more recent data show that tumor cell clones may undergo marked diversification over space and time as a result of mutational divergence on the one hand and selection by, for example, immunologic attack on the other hand (1–7). It has further become evident that immune attack appears to vary according to its location in the cancer tissue: when quantifying the density of tumor-infiltrating immune cells by microscopy, the spatial distribution of immune cells with respect to tumor cells needs to be taken into account (8). Given that immunotherapeutics only benefit a minority of patients (9), heterogeneity of immune cell distribution might contribute to treatment failure; and in that case small biopsy samples might not accurately represent the tumor microenvironment and immune status of a patient's tumor. In clinical practice, however, it is difficult to sample multiple tumor areas at any time other than at surgical resection.

Valuable prognostic information can be gained from simple enumeration of T-cell infiltrates in both human papillomavirus expressing [HPV(+)] and HPV independent [HPV(−)] head and neck squamous cell carcinoma (HNSCC; refs. 10, 11). Transcriptomic analysis of primary HNSCC using RNA-sequencing (RNA-Seq) allows further characterization of global molecular signatures (12) and detailed assessment of immune infiltrates (13). Such evaluation of gene expression provides unprecedented insight into biological processes that occur in tissue. Like morphologic assessments, molecular data are shaped by tumor-driving mutations, tissue heterogeneity, disease progression, or selective pressures from treatment (14–16). Technical issues may further affect interpretation: in breast cancer, significant differences were identified between diagnostic tissue core biopsies when compared with tumor excisions. Surgical ischemia was identified as a significant confounding factor in transcriptomic analysis (17).

Knowledge of variability is critical for the understanding and effective evaluation of immunotherapy interventions using checkpoint inhibitors (e.g., anti-CTLA4, anti-PD1, and anti-PDL-1 antibodies; refs. 18, 19). The mechanisms that underpin immunotherapeutic success as opposed to treatment failure remain poorly understood but must be linked to the particulars of the preexisting immune infiltrate (5, 20, 21). To make sense of treatment failure therefore, and of any changes detected in longitudinal assessments, understanding of natural variability is critical, as is the understanding of how results from different assays (e.g., histologic and transcriptomic assessment) correlate with each other.

Using high-resolution transcriptomic evaluation with RNA-Seq (22), we assessed variability in HNSCC within an individual cancer between multiple samples (replicates). In addition, we wondered whether replicates at diagnosis would be similar or different to replicates from a resection specimen of the same cancer in one person, taken about a month later. We were particularly interested in immunologic readouts and included both HPV(+) and HPV(−) HNSCC. Up to six tumor samples (replicates) per patient were collected from a total of 16 patients. Gene expression analysis using RNA-Seq revealed consistent hierarchical clustering and a high level of correlation between replicates in the same patient. Detailed analysis of immunologic gene signatures showed clustering of replicates by patient. Overall, our data show a stable transcriptional immunologic signature between multiple samples from the same patient and between diagnosis and surgical resection. These suggest that if applied to treatment study, the approach has the potential to allow definition of treatment consequences. Our data further suggest that immunologic heterogeneity within the primary tumor is not the cause for treatment failure in HNSCC.

Following LREC approval and written informed consent, multiregion tumor samples (replicates) were collected under general anesthesia and were snap frozen immediately at diagnostic biopsy and surgical resection. Pre-resection replicates were taken before vascular ligature to minimize hypoxic time. Spatial heterogeneity was assessed by collection of multiple replicates at either diagnostic sampling or at resection. Samples were taken at least 10-mm apart to maximize the chance of capturing spatial heterogeneity. The comparison between replicates collected at diagnostic sampling to those replicates from the same patient at resection were further used to assess change over time (temporal heterogeneity), in the recognition of the fact that this comparison may also be shaped by any spatial heterogeneity in any given cancer. Patient demographics, clinical details and number of replicates are shown in Supplementary Table S1. For RNA extraction, cryosections (10 μm) were cut and RNA isolated with the RNeasy Mini Kit (Qiagen Ltd.). From 16 patients, 44 samples were collected. Samples with transcript bias or RIN score <7 were removed (n = 5); this led to the exclusion of two patients who as a result only had a single sample but no replicate. The full analysis included 14 patients and n = 37 tumor replicates; from 6 patients 15 samples were collected at a single timepoint, 8 patients had replicates (n = 22) from two timepoints, both diagnosis and resection, with an average of 26 days between them (Supplementary Table S1). Clinical information for an additional cohort of patients also assessed between diagnostic biopsy and surgical resection is shown in Supplementary Table S1.

RNA quality was assessed using the Agilent 2100 Bioanalyzer generating an RNA integrity number (RIN; Agilent Technologies UK Ltd.). Poly-A mRNA was purified from total RNA (100 ng) using the Poly(A) Purist Mag Kit (Life Technologies Ltd.), according to the manufacturer's instructions. RNA was converted into a library for sequencing using the TruSeq stranded mRNA Sample Preparation Kit (Illumina Inc.) followed by hybridization to the flow cell for single-end (SE 35 bp) sequencing on the HiSeq 2000 (Illumina Inc.). Low-quality SE read data from FASTQ files were trimmed or removed. High-quality reads were mapped to the human genome (hg19) using TopHat (ref. 23; version 2.0.9) and, following the removal of multimapping reads, converted to gene specific read counts for annotated genes using HTSeq-count (version 0.5.4; ref. 24). A sample was considered positive for HPV if there were greater than 500 HPV assigned reads against the HPV reference genome (12) and confirmed with Illumina Pathseq (data not shown; ref. 25). Data are available at: ArrayExpress accession E-MTAB-4546.

Tumor grade and differentiation were recorded from formalin fixed, paraffin-embedded (FFPE) tissue. IHC staining for CD8a (anti-CD8a antibody; clone: C8/144B; DAKO) was performed on each frozen tumor tissue replicate for each case using a region immediately adjacent to that used for RNA-Seq analysis. An insufficient amount of material remained for cases #6, #11, and #13. An additional group of 9 patients were added in which FFPE from diagnostic material was compared with resection specimen again using anti-CD8a for IHC evaluation. Tumor-infiltrating lymphocytes (TIL) were quantified using an Olympus CKX41 with an average of 10 high-power (×400) fields; an average intratumoral score per high-power field (HPF) was calculated.

Raw counts from RNA-Seq were processed in Bioconductor package DESeq, variance was estimated and size factor normalized. Correlations and statistics were performed in Corrplot 0.73, in an R statistical environment (3.3). Qlucore Omics software (3.2) was used to visualize the normalized RNA-Seq data using Principle component analysis (PCA) and hierarchical clustering. These were performed on data using the setting: mean = 0, variance = 1 normalization in the Qlucore Omics software (3.2). Where the clustering is agglomerative, average linkage criterion was used to determine the distance between sets of samples and hierarchical clustering is based on Euclidean distances, where the Euclidean distance will be equal to sqrt(1 − r)⁁2, where r is the Pearson correlation coefficient. To identify the 4,000 most variable genes, variance was set at 0.2879. Briefly genes with a SD less than a specified variance cutoff from the maximal SD were removed using Qlucore Omics Explorer (3.2) as previously reported (22). Data were visualized in heat maps where each row represents normalized gene expression values for a given gene; each column represents the gene expression for a given tumor: red shading denotes greater gene expression; blue shading denotes lower gene expression. Hierarchical clustering (cluster method = average linkage; distance measure = Pearson correlation) of genes and tumors based on their expression profile is reflected in the dendrograms to the left and the top of the heatmap, respectively. A curated list of immune genes relating to immune cell markers, effector function, and exhaustion/regulatory function were derived from the gene ontology terms GO:0002250 adaptive immune response, GO:0002449 lymphocyte-mediated immunity and GO:0002456 T-cell–mediated immunity. As a measure of similarity and dissimilarity for intrapatient and interpatient samples, we quantified Euclidean distances for the 4,000 most variable genes (mean = 0, variance = 1 normalized) in the R statistical environment (3.3) between all replicates from the same patient to all other replicates, where the Euclidean distance is equal to the sqrt(1 − pearson)⁁2. This measures distances between samples in the hierarchical tree, with smaller distances between replicates indicating a closer relationship. The median distances were calculated for all replicates to all other replicates (interpatient) and compared to the median distances from intrapatient replicates, for example, case #14 median distance for intrapatient derived from n = 6 and compared to median distance to all other cases n = 31 (each replicate comparison is plotted). The Wilcoxon matched-pairs signed rank test was used to compare median distances between replicates from the same patient to those from different patients. Reads per kilobase per million mapped reads (RPKM) were used to display expression of genes between replicates in dot plot comparisons. Spearman correlation analysis between CD8A gene expression and CD8 T-Cell IHC count were performed; correlation analysis was also carried out on CD8 T-cell IHC count between diagnostic biopsy and surgical resection from FFPE samples in the second cohort.

The top 4,000 most variable genes across the cohort were displayed by variance filtering of the RNA-Seq data. These were hierarchically clustered and displayed in a heatmap (Fig. 1A). The 37 replicates from 14 patients evaluated at either diagnosis or resection cluster precisely by patient in 9 of 14 cases. The same genes were visualized with PCA in Fig. 1B showing grouping of replicates by patient, complementing the hierarchical clustering. In addition to this, the assessment of all genes (n = 18,979) and all cases was performed using hierarchical clustering and PCA (Supplementary Fig. S1A and S1B) again showing replicates grouping by case. To quantify the level of variation between replicates, the Euclidean distance was assessed for each replicate (Fig. 2A). The graph displays each case and the median intrapatient (circle) distance compared with the median interpatient distance (diamond). The length of the line connecting circles and diamonds is a measure of difference when replicates from one case are compared with all other samples. A long line indicates large differences (e.g., case #13), conversely a short line highlights cases where this difference is small. The minimum and maximum distance for each replicate is shown in Supplementary Fig. S2. The individual Euclidean distances are presented as a scatter graph in Fig. 2B, demonstrating that tumor samples from the same patient (intrapatient) are significantly closer to each other than to those from other patients (interpatient; P ≤ 0.0001).

Correlation analysis of the top 4,000 most variable genes was carried out across all samples (Fig. 3). The r values in the correlation matrix range from r = 1 (highly correlated, yellow) to r = 0.5 (less correlated, blue with values <0.5 also shown as blue). It shows that intrapatient tumor replicates are more correlated to each other than interpatient replicates. Case #14 shows a variable correlation between samples, but there is a higher correlation between its own replicates than to samples from other cases. The global median correlation coefficient for intrapatient replicates was r = 0.82 compared with an r = 0.63 for interpatient replicates. When cases are displayed in a correlation matrix, cases #13 and #9 have similar gene expression; this was also visible in the heatmap (Fig. 1A). Our data suggest that each patient's tumor has a distinct transcriptomic landscape. Correlation analysis of all genes across all samples with the inclusion of the less variable genes increases the correlation between samples (r scale = 0.9–1); again replicates from the same patient are most similar to each other (Supplementary Fig. S3).

The prognostic importance of T-cell infiltration in HNSCC is well documented (10–12). Gene expression for immune cell lineage markers, effector function, exhaustion, and regulatory genes were therefore represented in a heatmap (Fig. 4). Samples are ordered by case number and cases with replicates at a single timepoint were grouped, as were cases with replicates from two timepoints. We observe remarkable visual consistency in the immunologic signature between replicates from the same patient including the markers of immune effector function (IFNG/GZMA; ref. 5) and targets of cancer immunotherapy (PDL1/CTLA4; ref. 26). In addition to this, the RPKM of CD3E, GZMA, IFNG, CTLA4, and PDL1 (CD274) are displayed as dot plots (Supplementary Fig. S4A–S4E). Cases #1, #8, #5, #6, and #7 contain a consistent and low relative level of immune-gene expression, whereas cases #16, #9, and #13 have consistently high immune-gene expression. However, in some cases, replicates #2, #3, #6, #10, and one sample from case #14, the level of immune-gene expression is more variable, this is reflected in the hierarchical clustering of the tumor replicates for immune genes only (Supplementary Fig. S5).

We evaluated how the gene expression for CD8A compared with the assessment of CD8 cell counts by standard IHC. Using tumor tissue sections that were immediately adjacent to those used for RNA-Seq, we were able to assess 11 of 14 cases. Gene expression for CD8A is presented as dot plots for cases grouped by whether we had replicates from a single timepoint or two timepoints (Fig. 5A). CD8⁺ counts (average of 10 HPF) evaluated by IHC are shown in Fig. 5B. Visually, the cases with high CD8A transcript count appear to be equally rich in CD8 T cells. This can be quantified by assessing the Spearman correlation of the CD8A gene expression to CD8 T-cell IHC count with an r = 0.82 (Supplementary Fig. S6A). We asked if this degree of immunologic homogeneity assessed by IHC was unique to our cohort and evaluated the immune cell counts in an additional tumor cohort of paired samples (Fig. 6). The paired tumor replicates between diagnostic biopsy and surgical resection showed comparable levels of CD8⁺ T cells with a Spearman correlation of r = 0.95 (Supplementary Fig. S6B).

In solid cancers, morphologic heterogeneity is commonly observed and raises the concern of how representative of the whole cancer biopsy samples are. This is important clinically as biopsy samples are used to make treatment decisions and failure to appreciate heterogeneity could contribute to treatment failure. In addition, if tumor samples are used to gain detailed insights into molecular events resulting from therapeutic intervention(s), understanding transcriptomic heterogeneity becomes critical. Recent high-resolution studies have identified that beyond morphologic features, genetic features can vary within the geography of an individual's cancer (4, 6, 7). Therefore, we wished to quantify tumor and immunologic heterogeneity at high resolution using RNA sequencing in multiple tumor biopsy samples from individual patients. Our focus was to understand gene expression that reflects adaptive immune attack; as boosting this is the key for success of immunotherapy. We have previously shown that in HNSCC, TIL density is tightly linked to outcome, in both HPV(+) and HPV(−) cancers (11, 13). Here we used RNA-Seq to quantitate the immune signatures and global tumor gene expression profiles from tumor biopsies, separated in space (more than one sample at the same time) and time (samples taken at diagnosis and compared with samples from resection). Molecular immune cell quantification was then correlated to data generated by the current “gold standard,” IHC assessment of HNSCC patients, which is known to be important for outlook (10, 11, 13, 27).

The gene expression data generated using RNA-Seq for our patient cohort were assessed by PCA, hierarchical clustering, and correlation analysis of the top 4,000 most variable genes, as well as evaluation of specific immune-genes. Hierarchical clustering and Euclidean distance measures showed that multiple samples from the same patient were significantly more similar to each other than to those from other patients, this was also mirrored in the correlation analysis. Our data suggest a single sample is able to capture the key immune characteristics of the patient's primary HNSCC, irrespective of whether patient samples were taken at the same time or at different timepoints. Tumor biopsies clustered by patient and also displayed similar expression levels of the immune genes. In contrast, there is considerable variability between patients (11, 13, 28). We have previously identified that the global transcriptomic signature of T cells is similar between HPV(+) and HPV(−) tumors, but that HPV(−) cancers have a lower overall density of immune cell infiltrates (11, 13). We find here that the similarity of tumor replicates is maintained in both HPV(+) and HPV(−) HNSCC, but recognize that our series is limited by its size.

Overall our data are suggestive of a consistent immune signature in an individual tumor across both location and time, at least within the window between diagnosis and resection; additionally, our data demonstrate that tumors can be grouped by the immune-gene expression profiles. The evaluation of gene expression of key markers for immune attack (IFNG/GZMA; ref. 5) and targets of cancer immunotherapy (PDL1/CTLA4; ref. 26) further confirm a stable immune signature in the individual patient.

Reassuringly, gene expression of CD8A correlated strongly (r = 0.82) with the CD8 T-cell count assessed by manual counting of cells using IHC. Cell counts generated by IHC are also stable over time and space, this observation was confirmed in an independent cohort of paired samples from biopsy and resection specimens.

Nonetheless intrapatient tumor replicates are not a perfect match with each other and in some cases (cases #2, #3, #10, #11, #14, and #15) the tumor samples did not cluster perfectly by patient: In case #10, gene expression differences were observed at the immune-gene level which were also reflected in the IHC of CD8. Conversely, case #2 displayed a stable CD8 and immune cell marker signature but variability in the expression of effector function genes. Cases #2, #6, #10, and #14 were laryngeal tumors, which have been highlighted as having a higher level of mutational heterogeneity compared to other sites in HNSCC (29). We are currently evaluating whether the degree of transcriptomic heterogeneity links to the genomic heterogeneity identified by other groups (1, 30).

Our data support that in HNSCC the global tumor and adaptive immune signatures are stable across space and time between replicate samples from the same tumor. Our data have a number of implications. They suggest that immunologic heterogeneity is not likely to be a key reason for immunotherapy failure for the primary cancer in HNSCC; paired analysis of primary and metastatic disease is necessary to evaluate this question further for patients with later stage disease. Our data also imply that it is important to examine in other tumors whether a consistent immunologic “fingerprint” is present even in cancers with known genetic heterogeneity. Finally, our data support that from a background of transcriptomic stability in untreated patients, RNA sequencing may be useful for the detection of change resulting from treatment effects.

No potential conflicts of interest were disclosed.

Conception and design: O. Wood, J. Clarke, G.J. Thomas, P. Vijayanand, E. King, C.H. Ottensmeier

Development of methodology: O. Wood, J. Clarke, G.J. Thomas, P. Vijayanand, E. King, C.H. Ottensmeier

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): O. Wood, A.H. Mirza, P. Vijayanand, E. King, C.H. Ottensmeier

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): O. Wood, J. Clarke, J. Woo, C.H. Woelk, P. Vijayanand, E. King, C.H. Ottensmeier

Writing, review, and/or revision of the manuscript: O. Wood, J. Clarke, C.H. Woelk, G.J. Thomas, E. King, C.H. Ottensmeier

Study supervision: C.H. Ottensmeier

This work was funded by Cancer Research UK (grant number C491/A15951).

The authors would like to thank the patients and clinical trials assistants for providing and handling of the samples, respectively. The authors acknowledge the use of the IRIDIS High Performance Computing Facility, and the research histology team with their associated support services at the University of Southampton, in the completion of this work. Finally, the authors would like to thank the La Jolla Bioinformatics and sequencing core for processing of the transcriptome libraries.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.