Metastases and tumor recurrence have a major prognostic impact in head and neck squamous cell carcinoma (HNSCC); however, cellular models that comprehensively characterize metastatic and recurrent HNSCC are lacking. To this end, we obtained genomic, transcriptomic, and copy number profiles of the UM-SCC cell line panel, encompassing patient-matched metastatic and recurrent cells. UM-SCC cells recapitulate the most prevalent genomic alterations described in HNSCC, featuring common TP53, PI3K, NOTCH, and Hippo pathway mutations. This analysis identified a novel F977Y kinase domain PIK3CA mutation exclusively present in a recurrent cell line (UM-SCC14B), potentially conferring resistance to PI3K inhibitors. Small proline-rich protein 2A (SPRR2A), a protein involved in epithelial homeostasis and invasion, was one of the most consistently downregulated transcripts in metastatic and recurrent UM-SCC cells. Assessment of SPRR2A protein expression in a clinical cohort of patients with HNSCC confirmed common SPRR2A downregulation in primary tumors (61.9% of cases) and lymph node metastases (31.3%), but not in normal tissue. High expression of SPRR2A in lymph node metastases was, along with nonoropharyngeal location of the primary tumor, an independent prognostic factor for regional disease recurrence after surgery and radiotherapy (HR 2.81; 95% CI, 1.16–6.79; P = 0.02). These results suggest that SPRR2A plays a dual role in invasion and therapeutic resistance in HNSCC, respectively through its downregulation and overexpression.
The current study reveals translationally relevant mechanisms underlying metastasis and recurrence in HNSCC and represents an adjuvant tool for preclinical research in this disease setting. Underlining its discovery potential this approach identified a PIK3CA-resistant mutation as well as SPRR2A as possible theragnostic markers.
Development of regional or distant metastases and recurrence after definitive therapy are arguably the most outstanding issues in the management of patients with head and neck squamous cell carcinoma (HNSCC). Presence of lymph node metastases has a major prognostic impact, as it reduces survival of HNSCC patients roughly by half (1). Moreover, lymph node metastases do not necessarily respond to given therapeutic approaches in the same way as their matched primaries do, an observation likely underlying different intrinsic biological mechanisms (2). Such a heterogeneous behavior poses an obvious challenge to achieve simultaneous local (primary tumor) and regional (lymph node metastases) disease control. Indeed, while most early-stage (with the exception of oral squamous cell carcinoma) and some advanced-stage primaries can be usually managed by radiotherapy (with or without chemotherapy/cetuximab), presence of lymph node metastases often requires combined approaches, encompassing neck dissection and irradiation (3). Distant metastases develop in approximately 10%–30% of HNSCC patients after definitive therapy and carries a meager prognosis (approximately 5% overall 5-year survival; ref. 4).
With respect to therapeutic resistance, recurrences after definitive therapy in HNSCC occur in 30%–50% of the cases, depending on primary tumor anatomic location, stage at presentation, and human papillomavirus (HPV) status among other factors (5). Management options for treatment failures after multimodal therapy, including extensive salvage surgery, reirradiation, or palliation, are limited and often of limited efficacy. Moreover, such approaches carry almost systematically high rates of severe toxicity and/or functional impairment of breathing, swallowing, airway protection, and speech production (6, 7).
The inescapable corollary of these clinical considerations is the pressing need to develop and implement more effective therapeutic strategies for metastatic and recurrent HNSCC.
The ability to tailor individual management approaches and to implement them requires thorough characterization of the molecular alterations in HNSCC and extensive preclinical validation. In this sense, concerted next-generation sequencing (NGS) approaches in HNSCC have been published over the last few years, providing an unprecedented amount of information on the most frequent alterations and mechanisms underlying tumorigenesis (8–12). With respect to preclinical validation, cancer cell lines and their related in vivo models represent major research tools in preclinical oncology when it comes to addressing issues related to tumor biology and target discovery. Comprehensive genomic characterization is often unavailable for commonly used cell lines, a limitation that hinders the generalization of results and accounts at least to a certain extent for the discordant responses to novel therapies seen in preclinical versus clinical settings (13).
Recognizing the crucial nature of this issue, Barretina and colleagues (14) established genomic profiles of a large number of human cancer cell lines. In their study, the authors explored correlations between genomic background and responses to particular therapeutic agents, emphasizing the relevance of genetically driven management approaches. Several other studies provide similar datasets for cell lines from different tumor entities, and databases such as the catalogue of somatic mutations in cancer (COSMIC) systematically gather results from these studies to facilitate access to the research community (www.cancer.sanger.ac.uk/cosmic). With respect to HNSCC, Martin and colleagues (15) obtained full exome and transcriptome profiles of 21 HNSCC cell lines, showing that most typical genomic alterations found in HNSCC patients were recapitulated in these cell lines. While this study and others cover the spectrum of alterations in primary tumors, patient-derived models to address research questions pertaining to mechanisms of invasion and therapeutic resistance have so far not been comprehensively characterized from a genomic perspective.
Paralleling the sequencing approaches performed in primary HNSCCs, two recent studies investigated the genomic landscape of metastatic and recurrent HNSCC (16, 17). Hedberg and colleagues (16) performed whole-exome sequencing from 13 HNSCC lymph node metastases and 10 recurrent tumors and their matched primaries, while Morris and colleagues (17) carried out a targeted mutational approach in 410 cancer-related genes and copy number profiles of several head and neck tumors, including 53 HNSCCs.
Here we provide comprehensive genomic profiles of a HNSCC cell line panel derived from primaries and their patient-matched metastases/recurrences (18). Our results show that the UM-SCC panel is representative of HNSCC from a molecular perspective, as the most common aberrations found in HNSCC are recapitulated in the cell lines. Moreover, the current datasets underline the discovery potential of the UM-SCC panel and represent a unique toolset for the preclinical study of invasion and treatment resistance mechanisms in HNSCC.
Materials and Methods
Patients' tissue and clinical data were collected after approval by the Bernese Cantonal Ethical Committee (Protocol 050/14). The UM-SCC cell lines were obtained from surgically excised tissues from patients treated at the University of Michigan (Ann Arbor, MI) who gave written informed consent for the use of their tissues for research studies, including the development of permanent cell lines.
Cell lines and reagents
UM-SCC cell lines were provided by Professor T. Carey (University of Michigan, Ann Arbor, MI). Culture, maintenance, and cell line identity confirmation were performed as described previously (18). For group analysis, cell lines were divided in invasive/metastatic (UM-SCC-10B, -14C, -17B and -22B) and recurrent (UM-SCC-11B, -14B, -74B and -81B), and compared with their matched primaries (“A” series).
PI3K was inhibited with pictilisib (GDC-0941, AbMole BioScience). This drug was dissolved in DMSO and stored at −20°C.
For chromosome analysis, cells were grown as monolayers on coverslips (35 mm dish, 22 × 22 mm Glass, MatTek corporation) and harvested in situ. Cultures in exponential growth were arrested at metaphase by adding 10 μg/mL colcemid (KaryoMax Colcemid, Gibco) for 20 minutes at 37°C. Cells were then harvested with hypotonic saline solution (60 mmol/L KCl), fixed with fresh fixative (3:1 methanol: acetic acid, Scharlau), dried in a drying chamber (25°C and 40% humidity, CDS-5, Thermotron), aged for 1 hour at 93°C on a hot plate (Präzitherm, Huber&Co AG), and stained using Giemsa (azur eosin methylene blue solution, Scharlau) and trypsin (Difco Trypsin 250, Dr. Grogg Chemie). Capturing and karyotyping was performed with the commercial software Genikon (Nikon).
Transcriptome and exome profiling
Nucleic acids were extracted following lysis of cell pellets in 600 μL RTL lysis buffer (AllPrep DNA/RNA/miRNA Universal Kit, Qiagen). The lysates were homogenized in CK14 Precellys homogenization tubes (Labgene Scientific) using the Minilys homogenizer (Bertin Technologies). DNA and RNA were purified from the homogenized lysates using a robotic workstation (Qiacube, Qiagen) following manufacturer's instructions and their quality assessed using the Fragment Analyzer (Advanced Analytical Technologies Inc.). RNA sequencing libraries were prepared using the Illumina TruSeq Stranded Total RNA reagents (catalog number RS-122-2201; Illumina) according to the protocol supplied by the manufacturer and using 750 ng of total RNA. Exome sequencing libraries were produced as follows: 3 μg of genomic DNA was fragmented into 150–200 base pair fragments using the Covaris Sample Preparation System (S-series). Samples were prepared using SureSelect library prep kit (Agilent). The resulting adaptor-tagged libraries were then hybridized with exon biotin-coated baits from the SureSelect V5 kit (Agilent) for 24 hours at 65°C. Hybrids were captured with streptavidin-coated beads, PCR amplified and indexed.
Cluster generation was performed with the RNA and DNA libraries using the Illumina HiSeq PE Cluster Kit v4 cBot reagents (catalog number PE-401-400) and sequenced on the Illumina HiSeq 2500 using HiSeq SBS Kit V4 reagents (catalog number FC-401-4002). Sequencing data were processed using the Illumina Pipeline Software version 1.84. Exome sequencing initial number of reads averaged 164 ± 41 (SD) million per sample. The rate of mapping to the hGRC37 genome was of 99.9% ± 0.1% (SD).
RNA sequencing initial number of reads averaged 69 ± 13 (SD) million per sample. Reads were first trimmed to remove polyA and Illumina TruSeq adapter sequences using cutadapt (19), then aligned to the human reference hGRC37 genome using the STAR aligner (20). Reads that uniquely mapped to the reference genome averaged 92.6% ± 1.6% (SD). The number of counts were summarized at the gene level using featureCounts (21). Gene expression values were computed from the RPKM (reads per kilobase per million) values produced by the rpkm function of the edgeR package (version 3.16.5) by adding a pseudocount of 1 and log2-transforming the results.
Somatic mutation and copy number variation calling
Exome sequences were processed according to the Genome Analysis Toolkit (GATK) best practices pipeline, which involves trimming, alignment using BWA-MEM, marking of duplicates, local realignment around indels, and base recalibration (Supplementary Methods; ref. 22). Single-nucleotide polymorphisms (SNP) calling was performed against a panel of eight normal tissue samples taken from the salivary gland of independent patients obtained within the frame of a different study, using the algorithms of MuTect (23) and VarScan2 (24) with default settings. The variants identified with both MuTect and Varscan2 were annotated using ANNOVAR (25). To specifically find mutations involved in recurrence and metastatic progression, the variant calling was also performed comparing each recurrent or metastatic cell line directly against its primary tumor counterpart, by using the same algorithms.
We determined known cancer-relevant genes contained in chromosomal regions that were amplified or lost, using the curated cancer-driver gene database derived from Rubio-Perez and colleagues (https://www.intogen.org/downloads; ref. 26). The GISTIC2 software was used to find copy number alterations recurring in multiple samples.
Somatic copy number variation between primary cell lines and metastatic/recurrent counterparts was estimated using the VEGAWES R package (27) and represented in Circos plots (28). Use of VEGAWES for CNV (copy number variation) calling was previously validated by establishing karyotypes of two randomly selected cell lines (UM-SCC-17A and UM-SCC-17B). Karyotypes and VEGAWES calling was performed and showed a strikingly high concordance (Supplementary Fig. S1).
Transcriptome and pathway-level analyses
Paired differential gene expression analyses were performed with the generalized linear modeling functions glmFit and glmLRT of the edgeR package (version 3.16.5) using integer read counts as input (29). Only genes whose expression was consistently up- or downregulated in all 4 cell lines of either group (metastases and/or recurrences) were kept. The false discovery rate (FDR) values computed by edgeR were used to select the differentially expressed genes. Genes with an FDR < 0.25 were kept. Genes/cell lines biclustering and heatmaps were constructed using correlation distance as metrics and the average linkage agglomeration algorithm and heatmaps (R package nclust version 1.9.0, http://bcf.isb-sib.ch/Resources.html).
Pathway analyses were performed with the single sample GSEA (ssGSEA) and the preranked GSEA methods available at the GenePattern web interface of the Broad Institute (http://software.broadinstitute.org/cancer/software/genepattern/). Gene directional likelihood ratios resulting from differential gene expression analyses were used as weights in the preranked GSEA method and pathways with an FDR < 0.25 were kept. The single samples GSEA signature score of these significant pathways are presented in heatmaps created with the pheatmap package. We used the Hallmark collection of the MSigDB portal (http://www.broadinstitute.org/gsea/msigdb). The network of 24 proteins SPRR2A interacts with was obtained from the STRING database (medium confidence, all interaction sources, interactions are weighted with the combined score reported by STRING).
TCGA data analysis
Level 3 TCGA RNA sequencing gene expression data was acquired using the RTCGA package. We used the RSEM normalized expression values and log2-transformed them adding a pseudocount of 1. The gene classifier that allows categorizing TCGA patients into HNSCC molecular subtypes was kindly provided by V. Walter, University of North Carolina at Chapel Hill, Chapel Hill, NC (30). HNSCC subtype scores are defined as the correlation of the sample gene expression with the subtype centroids. TCGA patients and UM-SCC lines were then classified as subtypes with the highest correlation score. To analyze correlation of SPRR2A RNA levels and survival, we included patients for which both transcriptomic data (FPKM-UQ) of primary tumors and survival/follow-up information was available to calculate overall survival and disease-free survival. A total of 498 patients out of the 806 available (as of February 27, 2018) fulfilled these criteria.
Real-time PCR and genomic locus sequencing
Total RNA was extracted from 80% confluent cell cultures using TRIzol lysis according to manufacturer's instructions (Roche). Reverse transcription of mRNA was performed using the Omniscript RT Kit (Qiagen). Quantitative PCR of SPRR2A was performed using a 7900HT fast real-time PCR system with a TaqMan assay (Applied Biosystems). The mean Ct was determined from triplicate experiments and mRNA levels normalized to those obtained for 28S ribosomal protein. Changes in expression were determined by calculations of ΔΔCt.
A 1 kbp region around the SNP of interest was amplified by PCR, purified on agarose gel, and sequenced by the Sanger capillary method. The resulting electropherograms were manually analyzed to confirm the presence of the SNP of interest.
Cell proliferation, clonogenic assays, and cell migration
Cell viability was determined using a resazurin sodium salt reduction assay (Sigma). Briefly, at the indicated experimental endpoints cells cultured in 96-well plates were supplemented with medium containing 44 μmol/L resazurin. Resazurin reduction was colorimetrically estimated 1 and 6 hours later (570/600 nm) using a Tecan Reader (Tecan Group Ltd.). The measurements obtained in every treatment condition were normalized to vehicle-treated controls. Results are representative of three independent experiments.
For clonogenic assays, cells were plated and left to attach overnight. Irradiation was delivered using a 137Cs research irradiator (MDS Nordion) with a dose rate of 0.75 Gy/minute. The surviving fraction was normalized to the plating efficiency (PE = colonies formed/cells plated). Radiosensitization was evaluated using the radiation enhancement ratio (RER). A RER significantly superior to 1 according to one-sample t test was deemed to indicate radiosensitization (31).
Cell migration was assessed with the Oris Cell Migration Assembly Kit (AMS Biotechnology) after 3 days of culture. Briefly, cells were plated in 96-well plates in presence of a stopper which creates a 2-mm diameter cell-free circular area in the middle of each well. After cells are attached, cell stoppers are removed and cells may freely migrate. Pictures were captured at baseline and after 72 hours using a Leica DC 300F microscope. The invaded area in each treatment condition was determined with the ImageJ software (imagej.nih.gov/ij/).
HNSCC patient cohort
Validations in patient-derived tissues were based on a retrospective cohort including HNSCC formalin-fixed paraffin-embedded (FFPE) specimens of patient-matched primary tumors, lymph node metastases, and normal reference tissue. All tumors were stage IV-A and IV-B according to the UICC 7th edition (2009). All patients were treated with curative intent in the Departments of Radiation Oncology and Otorhinolaryngology-Head and Neck Surgery, Inselspital Bern University Hospital (Bern, Switzerland) and received radiotherapy, either as primary therapeutic strategy or following surgery for primary tumors, with or without concomitant systemic therapy (either cisplatin/5-fluorouracil or cetuximab). All patients received neck dissections (Supplementary Table S1).
Construction of tissue microarrays and IHC
Tissue microarrays (TMAs) were constructed using a new-generation tissue arrayer. Briefly, FFPE blocks were retrieved from the archives of the Institute of Pathology, University of Bern (Bern, Switzerland). Punches with a diameter of 0.6 mm were transferred to receptor TMA blocks. IHC staining was performed with an automated system BOND RX (Leica Biosystems). Sections of 4 μm were deparaffinized and rehydrated in dewax solution (Leica Biosystems). Endogenous peroxide activity was blocked in H2O2 for 4 minutes. Samples were incubated with specific primary antibodies for 30 minutes at room temperature: SPRR2A (1:500; Novus Biologicals). Slides were scanned using the Pannoramic Midi digital slide scanner (3DHISTECH Ltd.). Staining intensity was scored using the following scores: 0, no detectable staining; 1, weak staining; 2, moderate to strong in up to 50% of tumor cells; and 3, strong staining in more than 50% of tumor cells. For subsequent ad hoc group comparisons, scores 0 and 1 were considered as “low staining intensity” and 2 and 3 as “high staining intensity.”
Descriptive and comparative analyses were performed using GraphPad v5.03 (GraphPad Software, Inc.). Student t test was performed for intergroup comparison, association between stainings and clinicopathologic features were evaluated using χ2 test or Fisher exact test. Univariate survival analysis was plotted according to the Kaplan–Meier method and compared with the log-rank test. Multivariate survival analysis was performed using Cox proportional hazards regression analysis, using established risk factors for decreased disease control or survival. All P values were two-sided and statistical significance was set at 0.05.
Genomic profiles of UM-SCC cell lines recapitulate the molecular heterogeneity and most common alterations in HNSCC
We employed whole-exome sequencing, copy number variation (CNV), and RNA sequencing to characterize the UM-SCC panel. This panel, which includes 15 cell lines derived from 7 different patients (3 females and 4 males), encompasses all anatomical HNSCC subsites (Fig. 1A and B; ref. 18).
We applied the molecular HNSCC subtype classifier initially reported by Chung and colleagues (32) and subsequently refined and used by a number of studies, including the TCGA (9, 30). After validation of the subtype classifier in the TCGA cohort, UM-SCC cell lines were classified as atypical (n = 9), mesenchymal (n = 5), and classical (n = 1; Fig. 1C; Supplementary Fig. S2A and S2B). Importantly, no cell line displayed high centroid correlation for any given molecular subtype, reflecting a high degree of molecular heterogeneity, which fully parallels observations in HNSCC (30). Equally relevant in terms of expression profiles, we found that most cell lines cluster with their patient-matched counterparts, with the exception of UM-SCC-11A/B (Supplementary Fig. S2C). UM-SCC-11B clustered better with UM-SCC-74A/B than with its own matched counterpart, an observation most likely because these three cell lines have a clearly higher mesenchymal centroid correlation and phenotype than UM-SCC-11A (Fig. 1B and C). In addition, gene expression correlation was not necessarily determined by primary anatomic location or specimen's origin (primary, metastatic, etc.), a finding congruent with previous studies (Supplementary Fig. S2D; refs. 9, 16, 17, 30).
Also consistent with previous NGS studies, the mutational profiles of the UM-SCC panel featured common alterations of TP53, members of the Notch, the PI3K and the Hippo pathway, and of the histone methyltransferases KMT2D and NSD1 (Fig. 1D; Table 1; Fig. 2A; Supplementary Document S1; refs. 8, 9, 12, 15, 33). Conversely, however, mutations in other frequently altered genes such as CDKN2A, CASP8, and TGFBR2 were not found in this panel (Fig. 2A). Alterations in TRAF3, commonly associated with HPV-positive HNSCC, were not encountered in this panel of exclusively HPV-negative cell lines (9).
|.||UM-SCC panel .||Public databases/Publications (%) .|
|Gene .||% .||TCGA (9) .||Stransky (12) .||Agrawal (8) .||Pickering (10) .|
|.||UM-SCC panel .||Public databases/Publications (%) .|
|Gene .||% .||TCGA (9) .||Stransky (12) .||Agrawal (8) .||Pickering (10) .|
NOTE: Data were retrieved from the cBioPortal for Cancer Genomics.
To explore recurrently mutated genes in an unbiased manner, we focused on genes whose mutation rate across cell lines was high and whose length was short to middle, as long genes are statistically more likely to accumulate passenger mutations (Fig. 1E). Among these were NOD2, a gene involved in immune responses whose mutations have been reported in several types of cancers (34). Another interesting finding was presence of ADAMTSL4, member of a gene family involved in cellular adhesion and angiogenesis (35). In addition, our unbiased approach revealed a number of genes known to be consistently false-positive findings in mutational studies (36).
Finally, we performed gene expression pathway analysis using single sample gene set enrichment analysis (ssGSEA), with the hallmark gene signature collections available from the MSigDB molecular signature database. We selected pathways whose correlation matrices displayed at least one low correlation coefficient between two cell lines (r ≤ 0.5), to identify variable activation patterns across cell lines (Fig. 1F). Frequently altered HNSCC pathways were equally included. Mirroring previously published data, gene expression pathway analysis displayed very heterogeneous degrees of reliance on specific pathways (Fig. 1F). Considering such patterns of expression is most relevant in preclinical research when considering genomic correlates of responses to given therapies and allows selection models to address specific research questions in a genomically guided manner (9, 14). Altogether, these observations show that the UM-SCC panel is genomically representative of HNSCC.
Landscapes of UM-SCC cell lines display differentially distributed mutations with potential therapeutic implications
Two recent sequencing studies demonstrated that metastatic and recurrent HNSCC specimens are mutationally similar to their patient-matched primaries (16, 17). Mirroring such findings, variant allele frequency (VAF) in the UM-SCC panel displayed high concordance between cell lines derived from matched pairs of primaries and recurrences/metastases (Supplementary Fig. S3A and S3B).
Even though mutations exclusively present in metastatic/recurrent tumors seem to be overall rare, such differential distribution is potentially relevant for implementation of therapeutic approaches. For instance, Hedberg and colleagues (16) found DDR2 mutations in some metastatic/recurrent specimens, absent in their patient-matched primaries. Ectopic expression of such DDR2 mutations in cell lines rendered them exquisitely sensitive to dasatinib (16), underlying possible implementation of genomically guided personalized approaches in metastatic/recurrent HNSCC.
In spite of the small number of samples analyzed in this study, we identified a novel kinase domain PIK3CA mutation in UM-SCC-14B (recurrence), absent in UM-SCC-14A (primary) and UM-SCC-14C (metastasis; Fig. 2A). The allelic frequency of this mutated variant was 25% (Fig. 2B). Presence of the PIK3CA T2930A mutation (resulting in a F977Y protein change) was confirmed by direct sequencing of the region of interest (Fig. 2C). Next, to assess the impact of this mutation on cell proliferation, UM-SCC-14A, UM-SCC-14B, and UM-SCC-14C cells were exposed for 72 hours to the pan-class I PI3K inhibitor pictilisib. Interestingly, pictilisib was less effective at impairing proliferation in UM-SCC-14B than in the other two matched cell lines (Fig. 2D).
Given that the PI3K pathway is one of the most altered in HNSCC and PIK3CA among the top mutated oncogenes, PI3K signaling is considered a potential therapeutic target in HNSCC (9, 33). Our observation suggests existence of drug-resistant PIK3CA mutations and strongly stresses the importance of mutationally guided pretherapeutic patient stratification. We consulted the TCGA datasets through the cBioPortal for Cancer Genomics (www.cbioportal.org) to assess presence of this mutation. Although this specific mutation was not found, as UM-SCC-14B was derived from a recurrent tumor following radiotherapy, it is possible that the F977Y mutation be treatment-induced.
Differentially distributed mutations in more than one cell line were only found in metastatic lines (Supplementary Table S2). More specifically, DCHS2 mutations were found in three metastatic cell lines. DCHS2 encodes for the dachsous homolog 2 protein, also known as cadherin-27 or protocadherin-23, and is related to Hippo pathway signaling (37). In addition, three different mutations in the dystrophin gene (DMD) were found in UM-SCC-14C and UM-SCC-22B. Nevertheless, as both DCHS2 and DMD are very large genes, the likelihood that identified mutations are mere passenger events is high, especially in the case of DMD, the largest human gene. Most relevant perhaps was the A1G variant of TLR8 (Toll-like receptor 8) in 2 metastatic cell lines (Supplementary Table S2). This variant has been reported to activate NFκB signaling, as well as modify cytokine secretion and metabolic patterns in a number of infectious diseases, resulting in more efficient protection against disease. Along the same lines, TLR8 is thought to potentially confer enhanced tumor invasiveness by stimulating an NFκB-mediated proinflammatory response (38). Our results warrant further characterization of the role of TLR8 in metastatic HNSCC.
Comparative gene-copy number variation profiles in UM-SCC cell lines
As next step, we performed somatic CNV calling to find whether specific chromosomal aberrations were recurrently present in metastatic or recurrent cell lines (Fig. 3A and B; Supplementary Document 2). Using GISTIC2, we found three chromosomal regions significantly amplified in metastatic cell lines (4q35.1, 9q34.3, and 14q24.1), and only one significantly deleted in recurrent cell lines (21q22.3; Supplementary Fig. S4). The amplified regions in metastatic cell lines contained genes encoding for established drivers such as IRF2, FAT1, NOTCH1, MAX, and MLH3. With the exception of FAT1, all genes were overexpressed at least in 3 of 4 metastatic cell lines (Fig. 3A). The third amplified region (14q24.1) in metastatic cell lines contains, among others, MYC associated factor X gene (MAX) and MLH3, two genes whose mutations have been involved in pheochromocytoma and colorectal cancer development, respectively (39, 40). The only significantly deleted region in recurrent cell lines (21q22.3) contains RUNX1, a tumor suppressor that potentially interacts with p53 (41).
Taken together, our findings suggest that, while differences in CNV between primary and matched metastatic/recurrent tumors seem to be overall discrete, genes involved in progression and therapeutic resistance may be often altered and most likely play relevant roles in these processes.
Differentially regulated pathway analysis in metastatic and recurrent cell lines highlights the relevance of epithelial-to-mesenchymal transition as a mechanism of therapeutic resistance
To characterize differentially regulated pathways in metastatic and recurrent cell lines, we performed pathway enrichment analysis using the gene likelihood ratios obtained by differential gene analysis as weight factors (Fig. 4A). Among the frequently enriched pathways in recurrent tumors, we identified several proinflammatory signatures (IFNα/γ response, IL6/JAK/STAT3 and IL2/STAT5 signaling, and TNFα signaling via NFκB). Alterations of proinflammatory signaling are frequent in HNSCC tumors and cell lines, and have been linked to potential therapeutic resistance (9, 15, 42). Another process involved in resistance to chemotherapeutic agents, xenobiotic metabolism, was equally upregulated in 3 of 4 recurrent cell lines. Metabolic reprogramming is indeed thought to confer a selection advantage to tumor cells upon genotoxic therapy (43).
Most prominently, epithelial-to-mesenchymal transition (EMT) was upregulated in all recurrent cell lines, while no consistently enriched pathways could be identified in metastatic cell lines (Fig. 4A). In line with the observed EMT upregulation, we found that HNSCC subtype score progression in recurrent but not metastatic cell lines is characterized by an increase in the mesenchymal score (Fig. 4B). To further assess the differences in EMT's enrichment in metastases and recurrent cell lines, we performed differential expression analysis and compared the results for genes in the EMT group versus all genes in metastatic and recurrent cell lines (Fig. 4C). In contrast to recurrent cell lines, EMT genes with P values < 10−3 were absent in metastatic cell lines. Focusing on 14 EMT genes with an FDR < 0.25 in the differential expression analysis, we found these genes to be well connected and consistently upregulated in recurrent cell lines (Fig. 4D).
The consistent upregulation of EMT in recurrent but not metastatic cell lines is congruent with recent findings in pancreatic and mammary carcinomas, supporting the notion that EMT may have a more important role in therapeutic resistance than in invasiveness (44–46).
Comparative expression profiles of UM-SCC cell lines unveil a potential role of SPRR2A as effector of therapeutic resistance and potential prognostic marker in regionally invasive HNSCC
The differential expression analysis results were further used to identify individual genes involved in metastasis and/or recurrence in the UM-SCC panel. The ten most differentially expressed genes in metastatic and/or recurrent UM-SCC cell lines include ECM components (SRGN), genes involved in extracellular matrix (ECM) remodeling (MMP2 and KLK5), or epidermal differentiation and homeostasis (PRRX1, ALPK2, FLG, and IVL; Fig. 5A). A full list of differentially regulated genes (FDR < 1) is provided in Supplementary Document S3.
Interestingly, small proline-rich protein 2A (SPRR2A) was the only gene differentially expressed in both metastatic (FDR = 0.089) and recurrent (FDR = 0.004) cell lines (Fig. 5A). Downregulation was confirmed in 3 of 4 cell lines tested by means of qPCR (Fig. 5B). SPRR2A is one of the 14 SPRR-family of proteins and acts as a keratinocyte cross-linking protein, ensuring tissue integrity when facing injury and oxidative stress (47).
To gain further mechanistic insights into the potential roles of SPRR2A in metastasis and recurrence, we used Cytoscape to visualize the network of established SPRR2A interactors according to the STRING database and their fold-change results from the transcriptomic analysis. Genes of the SPRR family were coordinately downregulated both in metastatic and recurrent cell lines when compared with their matched primaries (Fig. 5C and D). Other genes involved in epithelial homeostasis, such as POU2F2 and CNFN, were similarly downregulated.
The oncogenic role of SPRR2A is complex. Specht and colleagues (48) demonstrated that stable SPRR2A transfection of cholangiocarcinoma cell lines injected intrasplenically led to inability to form liver metastases, however increasing local tumor invasiveness. Conversely, cell lines with downregulated SPRR2A readily developed liver metastases. From a mechanistic perspective, SPRR2A binds SH3 domain–containing kinases, including Src, Yes1, and Abl (47). In our transcriptomic dataset, SPRR2A downregulation did not seem to have an impact at the transcriptional level on Src-related kinases (Fig. 5C and D). However, given that interactions between SPRR2A and Src-related kinases occur at the protein-protein level, it is anticipated that SPRR2A downregulation may have an impact on prosurvival signaling mediated by these kinases.
While a limited number of studies have addressed the role of SPRR2A in hepatobiliary malignancies, the relevance of SPRR2A in HNSCC has not been explored. Thus, we proceeded to assess expression of SPRR2A by qPCR in 51 patient-matched sets of HNSCC primary tumors and lymph node metastases, demonstrating significant downregulation of SPRR2A mRNA levels in metastatic specimens (Fig. 6A). Next, to explore the potential clinical relevance of SPRR2A, protein expression levels were assessed by IHC in 426 primary tumors from 181 patients, 530 lymph node metastases from 236 patients, and 295 normal tissues (mainly salivary gland tissue) from 237 patients. SPRR2A was highly expressed in 49.53% of primaries and 26.98% of lymph node metastases, while its expression was not detected in normal tissues (Fig. 6B).
Matched sets of primary tumor, lymph node metastasis, and normal tissue, as well as full clinical and survival data were available for 147 patients. Patients were classified according to SPRR2A expression (low: scores 0 and 1; vs. high: scores 2 and 3), in primary tumors and lymph node metastases (Supplementary Table S3). In this set of matched samples, SPRR2A was expressed according to 4 different patterns (Fig. 6C). In 44.9% of the cases, SPRR2A was highly expressed in primary tumors only, in 14.29% in lymph node metastases only, and in 17.01% in both primary and matched lymph node metastases. SPRR2A was not detected in 23.81% of the cases (Fig. 6C). SPRR2A expression was not associated with tumoral stage, site, or development of distant metastases. SPRR2A low expression in primary tumors was significantly associated with local recurrences (P < 0.01) in the univariate time-independent analysis. SPRR2A low expression in lymph node metastases was significantly more prevalent in patients older than 60 years old (P = 0.02), while SPRR2A high expression in lymph node metastases was more prevalent in patients that experienced regional recurrences (P = 0.04) and those who got concomitant systemic therapy (P = 0.04) (Supplementary Table S3).
In the univariate survival analysis, only high SPRR2A expression in lymph node metastases was significantly associated with reduced regional recurrence-free survival (RRFS; HR = 2.80, 95% CI: 1.16–6.79, P = 0.02; Fig. 6D; Supplementary Table S4). In the multivariate analysis, nonoropharyngeal location of the primary tumor and high SPRR2A expression in lymph node metastases were both independent predictors of poor regional recurrence-free survival (Supplementary Table S4). Our multivariate models could not identify significant predictors of overall survival or local recurrence-free survival. Surgical resection of the primary tumor and oropharyngeal primary site predicted improved distant metastasis-free survival (Supplementary Table S4).
We next performed similar analyses using the TCGA datasets available for 498 patients and found no correlation between SPRR2A RNA levels and survival (Supplementary Fig. S5). This observation is not surprising, as data is derived from primary tumors and protein levels are not available in the TCGA cohorts.
Because the clinical data seems to suggest that high expression of SPRR2A may contribute to radioresistance and low expression may facilitate metastasis, we performed two supplementary in vitro approaches. For this purpose, UM-SCC-10A/B and UM-SCC-22A/B cells were plated and allowed to migrate for 72 hours. The UM-SCC-10 pair migrated faster than UM-SCC-22 cells. Nevertheless, for both pairs of cell lines, the SPRR2A low expressing 10B and 22B migrated significantly faster than 10A and 22A, respectively (Fig. 6E). To assess radiosensitivity, cells were irradiated with increasing doses (0, 2, and 4 Gy), and left to proliferate for 7 days. The SPRR2A-high line UM-SCC-10A was significantly more resistant to irradiation than its matched SPRR2A-low cell line UM-SCC-10B. We did not find significant differences between UM-SCC-22A and 22B (Fig. 6F).
In this study, we provide the most comprehensive molecular characterization of a panel of cell lines derived from metastatic and recurrent HNSCC specimens to date with the aim to provide a suitable tool for preclinical research in the metastatic and recurrent setting of HNSCC. Molecular portraits of HNSCC, including comparisons between HPV-positive and HPV-negative tumors, as well as primaries versus lymph node metastases and recurrent tumors, have emerged in recent years (8, 9, 11, 12, 16, 17). Implementation of novel findings in clinical practice, however, requires thorough preclinical validation in well-characterized models (14).
Our current results show that the UM-SCC panel is representative of the most common and most characteristic molecular alterations found in HNSCC. Indeed, mutational profiling of UM-SCC cells revealed typical alterations affecting TP53, PI3K, and the NOTCH pathways. A major aspect of interest was the observation of mutations exclusively present in metastases or recurrences, and the subsequent possibility to study differential sensitivities to diverse management approaches in primaries versus their matched metastases and/or recurrences based on their mutational background. For instance, it was recently demonstrated that exclusive presence of DDR2 mutations in lymph node metastases resulted in increased sensitization of HNSCC preclinical models to dasatinib (16). In this study and hinting at a potential resistance mechanism, we identified a novel kinase domain PIK3CA mutation (F977Y) that impairs responses to PI3K inhibitors in cell lines (Fig. 2D). This mutation was found in a cell line derived from a recurrent HNSCC specimen and is possibly treatment induced. Thus, given the potential relevance of the PI3K pathway as a therapeutic target along with the high rate of alterations along this pathway, elucidating the implications of such alterations in terms of treatment responses is of utmost importance prior to implementation of PI3K pathway inhibitors in recurrent HNSCC.
Mutations of components of the Hippo pathway were found in all cell lines (Figs. 1D and 2A). Aberrant Hippo pathway activation has been reported in several types of human malignancies, including lung, colorectal, ovarian, and liver cancer (49). Most prevalently mutated were DCHS2 and FAT isoforms 1 to 4. FAT1 is one of the most commonly mutated genes in HNSCC and is thought to act upstream of the Hippo pathway mediating crosstalk with the Wnt pathway, thus contributing to acquisition of invasive features (9, 49). Nevertheless, the oncogenic significance of FAT1, DCHS1, and DCHS2 mutations is not fully clear, as these genes are large and high prevalence of mutations in these genes may be due, in part, to their size. Moreover, while FAT and DCHS are well established Hippo pathway modulators in Drosophila melanogaster, their precise role in mammalians are not fully clear (49). Besides mutations, we found hemicentin 1 (HMCN1), an increasingly recognized key upstream regulator of Hippo pathway core kinases MST1/MST2 and LATS1/LATS2, to be consistently upregulated in metastatic cell lines (Fig. 5A). Given the marked differences in Hippo pathway members' mutations in primary versus metastatic and recurrent HNSCCs, the impact of Hippo signaling in such settings warrant further investigation.
CNV profiles were highly preserved in metastatic and recurrent versus primary tumors, with only 3 significant peaks in metastatic and 1 in recurrent cell lines (Fig. 3A and B). Such discrete differences are consistent with previously findings (16). Nevertheless, altered regions contained relevant drivers or tumor suppressors and underline potentially relevant mechanisms in development of metastases or recurrences.
Expression data revealed that transcriptomic differences between primary tumors and their matched metastases were more marked than in the case of recurrences. Among the frequently upregulated pathways in recurrent cell lines were pathways involved in the inflammatory response and, most significantly, EMT (Fig. 4A). While EMT is considered a hallmark process of metastasizing tumors, emergent evidence links EMT to treatment resistance (45). Indeed, recent studies in mammary and pancreatic tumor models convincingly show that impairment of EMT by blockade of either ZEB1 and ZEB2 or Snail and Twist does not lead to decreased metastatic ability (44, 46). Importantly however, inability to undergo EMT resulted in increased sensitivity to cyclophosphamide and gemcitabine in the same models (44, 46).
When considering individual cell lines, EMT signatures featured most prominently in UM-SCC-74A, UM-SCC-74B, and UM-SCC-11B, all derived from patients previously managed with chemoradiotherapy and with phenotypical characteristics of EMT (Fig. 1A and B; ref. 18).
The heterogeneity in terms of pathway analysis seen in the UM-SCC panel reflects the heterogeneity found in actual tumors and needs to be considered when exploring novel therapeutic approaches (9, 14, 30, 32). Relevant to this point, the characterization provided in UM-SCC panel allows investigation of mechanisms and target validation taking into consideration the molecular background of such cell lines. An important implication of such approach is the possibility to unveil heterogeneity in treatment responses and potentially identifying predictive biomarkers.
Finally, the finding that most clearly illustrates the discovery and validation potential of the UM-SCC panel concerns the identification of SPRR2A as potential effector of therapeutic resistance and biomarker in HNSCC. In cholangiocarcinoma and liver cancer models SPRR2A is activated by STAT-3 both in biliary epithelial cells and cholangiocarcinoma, coorchestrating EMT through interaction with ZEB-1 (47, 50). In addition, SPRR2A reduces p53 acetylation by impairing p300-p53 interactions, leading to inhibition of p53 transcriptional targets (50). Moreover, SPRR2A acts as an SH3 domain ligand, hence promoting protection against oxidative stress and DNA damage (47). In our FFPE validation cohort, most primary tumors (61.91%) highly expressed SPRR2A, while only a minority of lymph node metastases did (31.30%), perhaps reflecting necessary SPRR2A downregulation in certain clones within the primary tumor to give rise to metastases (48). Importantly, however, our results emphasize an additional aspect on the role of SPRR2A in therapeutic resistance. Indeed, as high SPRR2A expression in lymph node metastases is an independent risk factor for development of regional recurrence after therapy, it could be hypothesized that after metastatic colonization SPRR2A expression may be restored in neoplastic cells. In this context, SPRR2A would contribute to the stress response upon radiotherapy/chemoradiation by virtue of its antioxidative stress ability and through protein–protein interactions with multiple prosurvival signaling molecules through its SH3 domain ligand activity (including Src, Yes1, and others) (47). In vitro migration and radiosensitization assays seem to suggest that indeed high SPRR2A expression may contribute to radioresistance, while low SPRR2A expression may stimulate invasive features.
Factors that determine SPRR2A reestablished expression in metastases are yet to be fully investigated and may provide novel mechanistic insights with potential therapeutic relevance, as SPRR2A cannot be directly targeted.
In addition, it would be of high interest to determine SPRR2A expression levels in recurrent lymph node metastases after radiotherapy and performing prospective validation of SPRR2A as a biomarker for risk stratification.
Limitations of this study include the absence of HPV-positive cell lines and limited amount of cell lines included. However, it should be noted that availability of metastatic/recurrent and primary tumor patient-matched HNSCC cell lines, especially HPV-positive ones, remains limited.
Disclosure of Potential Conflicts of Interest
O. Elicin is a consultant/advisory board member for AstraZeneca, Merck, and Merck Serono. T.E. Carey has ownership interest (including stock, patents, etc.) in licensed cell lines to Millipore EMD. No potential conflicts of interest were disclosed by the other authors.
Conception and design: L. Nisa, D. Barras, M. Medová, D.M. Aebersold, R. Giger, M. Delorenzi, Y. Zimmer
Development of methodology: L. Nisa, D. Barras, M. S. Dettmer, R. Giger, L. Ho
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): L. Nisa, M. Paliaková, J. Koch, B. Bojaxhiu, O. Eliçin, M. S. Dettmer, R. Giger, M.D. Caversaccio, T.E. Carey, T.A. McKee
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): L. Nisa, D. Barras, M. Medová, M. Medo, M. Paliaková, J. Koch, P. Angelino, R. Giger, M. Delorenzi, Y. Zimmer
Writing, review, and/or revision of the manuscript: L. Nisa, D. Barras, M. Medová, D.M. Aebersold, O. Eliçin, R. Giger, U. Borner, M.D. Caversaccio, T.E. Carey, T.A. McKee, M. Delorenzi, Y. Zimmer
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): L. Nisa, D. Barras, M. Medová, J. Koch, B. Bojaxhiu, O. Eliçin, M. S. Dettmer, R. Giger, U. Borner, M.D. Caversaccio, M. Delorenzi
Study supervision: L. Nisa, D. Barras, M. Medová, R. Giger, M. Delorenzi, Y. Zimmer
This work was supported by Swiss National Science Foundation (grant no. 31003A_156816, to Y. Zimmer) and by a Bernische Krebsliga grant (to M. Medová). The authors gratefully acknowledge the help of Sabina Gallati, PhD and Nijas Aliu, MSc, from the Division of Human Genetics, Department of Pediatrics (Inselspital, Bern, Switzerland) for establishing karyotypes.
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.