Intrinsic Molecular Characteristics of the Immune Subtypes of Salivary Mucoepidermoid Carcinoma

doi:10.21203/rs.3.rs-869310/v1

Download PDF

Research

Intrinsic Molecular Characteristics of the Immune Subtypes of Salivary Mucoepidermoid Carcinoma

https://doi.org/10.21203/rs.3.rs-869310/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background: Characterising the tumor microenvironment (TME) and immune landscape of cancer has been a promising step towards discovering new therapeutic biomarkers and guiding precision medicine; however, its application in salivary mucoepidermoid carcinoma (MEC) has been sparse. Here, we conducted a comprehensive transcriptomic study to understand the properties of the TME and immune profiles of MEC.

Methods: Molecular features in heterogeneous immunophenotypic subgroups of MEC and their intrinsic characteristics were determined by applying bioinformatic and immunoinformatic analyses on 20 matched primary MEC RNA-seq data.

Results: In this study, distinct two immunophenotypic subgroups, hot and cold MECs, were uncovered with their distinct molecular features, and potential immune-oncologic therapeutic options for each subgroup were suggested. In search for immunophenotype defining molecular features, tumor mutational burden, CRTC1-MAML2 fusion status, and its fusion neoantigen were not discriminable factors. However, we demonstrate that a significant inverse correlation between lipid metabolism activity and immunogenic state, and lipid metabolic regulators, such as MLXIPL and FASN, which are associated suppression of immune activity, were under-expressed in the immune-hot subgroup, contributing significant role in high immunity of immune active subgroup.

Conclusions: Our study has shown heterogeneous immunophenotypic MEC subgroups with their distinctive molecular characteristics and provided potential treatment options tailored to the immune context, which yields, for the first time, new insights into TME of MEC and may help the next step to studying this uncharted cancer.

Cancer Biology

Salivary mucoepidermoid carcinoma

tumor immune microenvironment

immune therapy

biomarker

Mucoepidermoid carcinoma (MEC) is the most common malignant neoplasm of the salivary gland, accounting for approximately 30–35% of all malignancies. MEC arises from heterogeneous cell types, including epidermoid, mucus-producing cells, and their intermediate cells (1). Patients with MEC generally show a relatively good prognosis, with a 5-year survival rate of 88.6% (1). However, poor prognosis factors such as age, vascular invasion, cytologic atypia, and high mitotic frequency, affect treatment efficiency; for instance, about 25% of the high-grade MECs eventually recur, limiting the 5-year survival rate to 50–56% (1). Unlike many other major cancers that have benefitted from genomics-based biomarker development, patient stratification, and precision medicine, the rarity and lack of underlying molecular/genomic characterisation hinder the use of such advances in MEC. Only a few driver mutations are known for MEC, including the t(11;19) translocation, which involves fusion of the transcriptional coactivator genes MAML2 and CRTC1 (2), copy number variations (CNVs) in DCC, SMAD4, and GALR1 (3) and additional somatic mutations in TP53 and POU6F2 (4). Except for the CRTC1-MAML2 fusion, which is known to have a favourable prognosis (2), no major genetic aberrations have been associated with or utilised for clinical use. As the discovery of novel high-frequency mutations is unlikely, alternative approaches with possible clinical implications are required for better treatment of MEC.

Recent advances in the profiling of immune cells and the tumor microenvironment (TME) warrant new therapeutic opportunities for immune-oncology, such as immune checkpoint blockades. Continuing efforts to understand the interplay between tumors and their surrounding components have led to the discovery of novel subtypes, uncovering of new biomarkers for improved prognosis, and optimisation of clinical management (5–7). In salivary gland cancer, the immune microenvironment and neoantigen profiles have been analysed for three subtypes: adenoid cystic carcinoma, myoepithelial carcinoma, and salivary duct carcinoma (8). However, no such study has been performed on MEC.

In the present study, we characterised the tumor immune microenvironment of MEC using whole transcriptome analysis of 20 tumor samples. Based on the level of immune and stromal cell infiltration, we profiled the TME and immune landscape to reveal immunophenotypically distinct subgroups. Deeper analysis of gene expression and protein interaction networks identified subgroup-associated intrinsic molecular signatures of tumors and potential immunotherapeutic biomarkers. We expect that our study will provide a basis for precision immuno-oncology of MEC.

Sample acquisition and preparation

Twenty patients who were diagnosed with MEC and treated at the Yonsei Head and Neck Cancer Center between March 2012 and April 2018, were enrolled in the study. Fresh frozen tissues from primary MECs and their matched normal tissues were retrieved and collected during the operation, after approval by the Institutional Review Board at Severance Hospital, Yonsei University College of Medicine (IRB 255-001). Informed consents were obtained from all patients. The inclusion criteria were as follows: 1) patients aged > 18 years; 2) patients with biopsy-proven MEC with a minimum diameter of 1 cm; and 3) patients with no previous history of chemotherapy or radiotherapy. All patients were treated surgically by complete resection, according to the National Comprehensive Cancer Network (NCCN) guidelines. The general clinical characteristics of the patients are presented in Table 1.

Total RNA sequencing

Total RNA was extracted using the RNeasy mini kit (Qiagen, Hilden, Germany) according to the manufacturer's recommendations. RNA quality was evaluated using the Agilent RNA 6000 Nano Kit and Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA). Total RNA (100 ng) was used to prepare sequencing libraries utilizing the TruSeq stranded total RNA sample preparation kit (Illumina, CA, USA), which combines RiboZero rRNA depletion with a strand-specific method similar to the dUDP method. The quality and quantity of cDNA libraries were evaluated using Agilent 2100 BioAnalyzer (Agilent) and the KAPA library quantification kit (Kapa Biosystems, MA, USA) according to the manufacturer’s library quantification protocol. Following cluster amplification of denatured templates, paired-end (2 × 100 bp) sequencing was performed using the Illumina HiSeq platform.

Data processing

FastQC (version 0.11.7) was used to determine read quality. No further read trimming and filtering were performed. Reads obtained from RNA sequencing were mapped against the UCSC hg38 reference genome (FASTA) and hg38 gene annotation (GTF) using STAR (version 2.6.0c) (9) under two default passes. The resulting alignment files (BAM) were used to quantify the read count per gene, using HTSeq (version 0.10.0) (10). All count data were normalised using the DESeq2 software package (version 1.24.0) (11) to normalise the raw counts into variance-stabilizing transformation (VST).

Calculation of TME

For constructing concise TME genes specific to microenvironment-related cells, gene signatures from CIBERSORT (12) and MCP-counter (13) were used. All 22 gene signature sets (LM22) of CIBERSORT were utilised; however, as they do not cover fibroblast and endothelial cell gene signatures, these were imported from the MCP-counter. Subsequently, using the gene set variation analysis (GSVA) software package (version 1.32.0) (14) single-sample gene set enrichment analysis (ssGSEA) function, a ssGSEA was performed on normalised expression data to measure the abundance of each cell subset in each patient.

Tumor subgrouping with TME estimation

ssGSEA was used to measure the abundance of immune and microenvironmental cells in the sample, using the ssGSEA function of the GSVA software package (version 1.32.0) (14). In total, 24 cell types were investigated by combining two gene signatures: the LM22 signature of CIBERSORT (12) for 22 major immune cell types, and MCP-counter (13) to supplement two microenvironmental cell types: fibroblasts and endothelial cells. Based on the ssGSEA scores, hierarchical clustering was performed using Euclidean distance and complete linkage methods utilizing the ComplexHeatmap software package (version 2.0.0) (15). The two largest clusters were used to define two immunologically distinct subgroups, which were later designated as immune-hot and immune-cold subgroups.

Robustness of the subgroups was tested by constructing clusters with an independent TME score from ESTIMATE (version 2.0.0) (16). To produce ESTIMATE scores, T-cell infiltration score (TIS), immune infiltration score (IIS), cytolytic score (CYT), and IFN-γ score were calculated as follows (17): 1) TIS was defined as the mean of the standardized values for CD8 + T, central memory T, and effector memory T cells, and Th1, Th2, Th17, and Treg cells; 2) IIS of a sample was defined as the mean of the standardized values for macrophages, DC subsets (total, plasmacytoid, immature, and activated), B cells, cytotoxic cells, eosinophils, mast cells, neutrophils, NK cell subsets (total, CD56bright, and CD56dim), and all T-cell subsets excluding T-gamma delta and T-follicular helper cells; 3) CYT score was calculated from the geometric mean of the expression of granzyme A (GZMA) and perforin (PRF1), which was demonstrated in the study by Rooney et al. (18); and, 4) IFN-γ signalling was scored with ssGSEA using the IFN-γ gene set defined in REACTOME (http://www. broadinstitute.org/ gsea/msigdb/cards/EACTOME_INTERFERON_GAMMA_ SIGNALING).

Comparative analysis of subgroups

Differential gene expression analysis between the immune-hot and cold subgroups was determined using the following significance criteria: 1) adjusted P < 0.05, and 2) Log₂(fold change) ≥ 1, using the internally implemented function of the DESeq2 software package (version 1.24.0) (11). Gene ontology (GO) and pathway enrichment analyses were performed using the clusterProfiler software (version 3.12.0) (19). For subtype tumor intrinsic analysis, differentially expressed gene (DEG) analysis was performed using the following significance criteria: 1) adjusted P < 0.01, and 2) Log₂(fold change) ≥ 3. For the gene set enrichment analysis (GSEA), genes were ranked by their log fold-exchange values, then GSEAPreanked with default settings was used (20) and the 50 hallmark genes from MSIgDB collection were utilised (21). The enrichment score (ES) denotes which gene set is overrepresented at the top or bottom of the ranked gene list.

Variant calling for mutation burden and neoantigen analysis

The STAR-aligned reads were processed for data clean-up following the best practice for RNA-seq short variant discovery (https://software.broadinstitute.org/gatk/best-practices). To perform variant calling, the pile-up of the alignment (mpileup) was created for each normal/tumor pair using Samtools (version 1.3.1) (22) and used for VarScan 2.4.3 (23). The output variant calls (VCF) were converted into the mutation annotation format (MAF) using vcf2maf v1.6.16 (https://github.com/mskcc/vcf2maf) and annotated with the variant effect predictor (VEP) (24). Tumor mutational burden (TMB) was calculated as the number of non-synonymous mutations. To predict neoantigens, vaxrank (version 1.0.2) (25) was used for somatic mutations and fusion (26) on the CRTC1-MAML2 fusion sites. HLA types of the patients were inferred using OptiType (version 1.3.2) (27).

Analysis of isoform switching

HISAT2 (version) (28) was used to quantify the transcript-level expression of the UCSC hg38 reference. The mapped transcripts were assembled using Stringtie (28) and used for downstream analysis. Alternative splicing and isoform changes in normal to tumor progression were examined using IsoformSwitchAnalyzeR (version 1.8.0) (29) at default setting and per the recommended procedures in the manual. To assess the protein-coding potential of transcripts, CPC2 (http://cpc2.cbi.pku.edu.cn/batch.php) (web version 2.0) was used. Protein domains were predicted using HMMER (https://www.ebi.ac.uk/Tools/hmmer/ search/hmmscan) (web-version 2.41.1), and PFAM was used as a profile database. Signal peptides were predicted using SignalP-5.0 (http://www.cbs.dtu.dk/services/SignalP) (web-version 5.0) with eukaryote and short-output options.

Network-based functional analysis of genes

The significant candidate genes were identified using the STRING database (30). A combined score of ≥ 0.4 was used as the cut off for a significant interaction. We developed a network of isoform-switched genes and PPIs. To visualise the protein interaction network and analyse hub proteins, the CytoHubba application (31) in Cytoscape (v.3.7.2) (32) was utilised to implement seven calculation methods: degree, bottleneck, closeness, betweenness, eccentricity, EPC, and MCC.

Immunological subgrouping of MECs

We initially profiled the tumor microenvironment in 20 patients with MEC from Yonsei Head and Neck Cancer Center (Table 1) using ssGSEA with respect to 24 immune and stromal cell types (Fig. 1A) (see Materials and methods for analysis procedures and Supplementary Figure S1 for overall workflow). Hierarchical clustering of the 20 MECs with the inferred cell abundance revealed two subgroups: eight tumors with high immune cell infiltration (referred to as immune-hot) and 12 with low infiltration (immune-cold). The infiltration patterns were concordant among most immune cell types, including myeloid progenitor cells, NK/T lineage cells, and stromal cells, with only a slight variation in B and mast cells. In addition, we confirmed that these clusters were consistently reproduced in an independent analysis with a different measurement algorithm (ESTIMATE (16), see Materials and methods), showing clear immunological discrepancies between the two subgroups (Fig. 1A and Supplementary Table S1) (Wilcoxon rank-sum test, P = 4.5 × 10^− 4 for immune score and P = 6.0 × 10^− 4 for stromal score, see Materials and methods). These results represent immunological heterogeneity and the existence of subgroups of MECs.

Further analyses confirmed distinct immune activities between the subgroups. The expression levels of CD8A and CD8B were significantly higher in the immune-hot group (Fig. 1B, Wilcoxon rank-sum test P = 1.1 × 10^− 3 for CD8A and P = 3.0 × 10^− 3 for CD8B) than in the immune-cold group, demonstrating distinct activities of cytotoxic T-cells. Similarly, scores for TIS (i.e. the average value of central memory, effector memory CD4⁺ and CD8⁺ T cells, and Th1, Th2, Th17, and Treg cells), and overall IIS (17) were high in immune-hot MECs (Wilcoxon rank-sum test P = 3.2 × 10^− 5 for TIS and P = 6.4 × 10^− 5 for IIS). Moreover, we noted elevation in the CYT score (P = 4.1 × 10^− 3) (18) and IFN-γ score (P = 3.0 × 10^− 4) (33) (see Materials and methods) for the immune-hot subgroups (Fig. 1C). After observing high levels of immune infiltration and immunity in the immune-hot MECs, we further characterised the immunophenotypes of the two MEC subgroups, focusing on the correlation between immune-cell infiltration and cytotoxic functions. IFN-γ, one of the final products of antitumor T-cell activity, was strongly positively correlated with the TIS (Pearson correlation coefficient r = 0.76, P = 1.0 × 10^− 4) and overall IIS (r = 0.84, P = 3.1 × 10^− 6; Supplementary Figure S2A). Together, these results suggested that immune cell infiltration and activity level-based MEC stratification exhibited two robust immunophenotypic subgroups.

Subgroup-associated immunophenotypic characteristics

We then investigated the phenotypic characteristics that confer immunological differences between the immune-hot and -cold subgroups. We initially noted that the TMB and neoantigen burden did not differ between the subgroups (Wilcoxon rank-sum test, P = 5.7 × 10^− 1 and P = 6.8 × 10^− 1 for immune-hot and -cold subgroups, respectively). In addition, presence of the CRTC1-MAML2 fusion, which is a known marker for favourable prognosis and elevated immunologic activity (34), was not associated with the immune subgroups (6/8 and 6/12 in immune-hot and -cold, respectively, Fisher's exact test, P = 3.7 × 10^− 1). Furthermore, no potential immunogenic neoantigens were predicted by the CRTC1-MAML2 fusion event (see Materials and methods). Based on these results, we presumed that the immunological differences between subgroups resulted from the machinery or regulatory mechanisms in the immunity cycle, rather than genome-level aberrations.

To profile the activities of subgroup-associated immune-regulatory mechanisms, expression levels of 78 immunomodulator (IM) genes (7) were used to infer the status of seven immune regulatory categories: antigen presentation, cell adhesion, immune checkpoint receptor, ligand, co-inhibitor, co-stimulator, and others (Fig. 2). We found that most of the IM genes and the categories, particularly cell adhesion (SELP, ITGB2, and ICAM1), antigen presentation (HLA Class I/II), co-inhibitor (PD-L1, PD-L2, BTN3A1, and BTN3A2), co-stimulator (CD80 and CD28), and ligand (CXCL9, CXCL10, IL2, IL10, OX40L, CD40LG, and IFNG), were consistently elevated in the immune-high group (all, P < 5.0 × 10^− 2), which implied that the immunologic differences may have risen at the early stages of the immunity cycle, such as the antigen-presenting machinery (APM). Of note, level of the APM (33) was significantly low in the immune-cold subgroup (P = 7.3 × 10^− 3) compared to that in the immune-hot subgroup (Supplementary Figure S2A), and was positively correlated with the TIS (r = 0.56, P = 9.5 × 10^− 3) and IIS (r = 0.67, P = 1.1 × 10^− 3) (Supplementary Figure S2B). These results indicate that the higher cytolytic potential in the immune-hot MECs could be a result of the higher level of APM expression with a cytolytic feed-forward loop along with an increased amount of immune cell infiltration, thereby leading to elevated IFN-γ levels and anti-cancer immune activity.

Transcriptomic landscape and hub genes of the immune subgroups

To illustrate the biological properties of immune-hot and immune-cold MECs, we conducted a functional analysis in the biological contexts based on gene expression profiles. We identified 1518 DEGs between the two subgroups, 1076 of which were upregulated in the immune-hot and 442 in the immune-cold MECs (see Materials and methods) (full list available in Supplementary Table S3). Gene ontology (GO)-based enrichment analysis (see Materials and methods) uncovered a global overrepresentation of innate and adaptive immune-related terms in the immune-hot subgroup, including T-cell activation, regulation of leukocyte activation, and regulation of lymphocyte activation (all, adjusted P-value < 1.0 × 10^− 2) (Supplementary Fig. 3A top) (Supplementary Tables S3 and S4). Likewise, cancer hallmark-based GSE analysis (see Materials and methods) revealed enrichment of inflammatory response and interferon-gamma response in the immune-hot group, confirming increased antitumor immune activity (Supplementary Fig. 3B). For the immune-cold group, ion transport-related terms, epidermis development, cornification, and chloride transport from GO analysis (all, adjusted P-value < 1.0 × 10^− 2) (Supplementary Fig. 3A, bottom), and metabolism-related cancer hallmark, such as fatty acid metabolism and adipogenesis, were enriched (Supplementary Fig. 3C).

Next, to improve understanding of factors that are shaping heterogeneity of immune state in MEC immune subgroups, we sought to assess tumor-intrinsic traits of them by comparing with matched adjacent normal tissues. First, we identified 893 DEGs with tumor-specific expression (see Materials and methods). Of these, 467 DEGs were over- or under-expressed from immune-hot MEC tumor and 326 from immune-cold MEC (full list available in Supplementary Tables S5-S8). As the top-ranked upregulated genes observed only in immune-hot MEC, CXCL13, which is a B cell chemoattractant and an important factor in TLS formation and initiation of secondary lymphoid organ development (35) (adjusted P-value = 1.55 X 10^–10 and log2 fold change 5.89) and DUOX2, responsible for the innate immune response (36) (adjusted P-value = 2.74 X 10^− 9 and log2 fold change 5.94), were observed, whereas MUC6 (mucin 6) was top ranked in immune-cold MEC (adjusted P-value = 1.34 X 10^− 6 and log2 fold change 5.47). POSTN1 (periostin), which affects tumor microenvironment remodelling during tumor progression, was one of the top ranking genes overexpressed in both immune subgroups (log2 fold change > 5). Then, also, to integrate biological impact of isoform-level changes in forming polar tumor microenvironment, we measured the isoform switches as mentioned in Kristoffer Vitting-Seerup et al. (29) and we identified 173 and 85 genes that underwent isoform switching in tumorigenesis in the immune-hot and -cold MECs, respectively. Among the 640 genes in immune-hot MECs, we identified T-cell regulatory genes (ADA, IL4R, and NFKBIZ) and lipid metabolism regulatory genes (MLXIPL, LEP, and FASN). From 411 genes from immune-cold MECs, we observed collagen and bone forming related genes (MMP13, COL10A1, and ITGB6) (see Materials and methods; full list available in Supplementary Tables S9-S10). These gene- and isoform-level intrinsic features provide a broad understanding of the distinctive intrinsic features that reflect the heterogeneous immunophenotypic state of each immune MEC subgroup.

To better study the functional significance of subgroup-specific genes at the system level, we constructed Protein-protein interaction (PPI) network by the union of DEGs and significant isoform-switching genes (Supplementary Table S11), which formed 610 nodes and 1429 edges, and 403 nodes and 896 edges for immune-hot and -cold subgroups, respectively (Supplementary Table S12). We then selected hub genes, important nodes with several interaction partners, with seven different topological analysis algorithms provided by cytoHubba (31), which appeared more than twice in the union of all top 20 genes from each algorithm, and resulted in 35 genes (15 upregulated and 20 downregulated) and 31 genes (2 upregulated and 29 downregulated; Supplementary Table S13) for immune-hot and -cold subgroups, respectively, and then searched for their functionality.

In the immune-hot MECs, the genes that are related to chemokine regulation and ECM (POSTN, COL11A1, MMP13, LUM, VCAN, and BGN) and lipid metabolism regulation (SCD, LPL, MLXIPL, PCK1, and DGAT2) were selected as hub genes (Fig. 3A). Overexpression of the ECM-related genes and lipid metabolism-related genes in cancers is associated with cancer progression, immune suppression, and tumor aggressiveness (37, 38). Here, as expected, ECM-related genes were over-expressed in the immune-hot MECs; however, lipid metabolic-regulator genes were under-expressed in the immune-hot tumor MEC samples, and this particular feature, downregulated lipid metabolism, was also captured in above GSEA analysis. Based on the aforementioned findings and inverse correlation between signature scores of fatty acid metabolism and adipogenesis obtained from cancer hallmark genesets with cytotoxic activity (Spearman correlation coefficient r = -0.49, P = 3.1 × 10^− 2 and r = -0.55, P = 1.4 × 10^− 2, respectively; Fig. 3B), we surmised that under-expressed lipid metabolic regulator genes in immune-hot MEC samples could be novel molecular characteristics and may be linked to agonistic anti-tumor immune activity in this group, as their overexpression is often related to a decrease in immune activity (38, 39). In the immune-cold MEC samples, all selected genes, except POSTN and RUNX2, were downregulated and are known to function as chemokines and their receptors (CXCR1, CXCR2, CCL2, CCR3, and CCR10) and in leukocyte migration and adhesion (ELANE, LEP, and SELL) (Fig. 3C).

Additionally, POSTN and RUNX2 were recurrently overexpressed in tumor samples from both MEC immune subgroups. POSTN participates in ECM structure formation and organisation (40, 41), and RUNX2 is a key transcription factor for osteoblast differentiation (42). At this point, we do not know their exact roles in MECs; however, based on our observations and previous studies, we conjecture that they could play a critical role in MEC ECM formation. In this analysis, we identified unique and common candidate molecular characteristics that were predicted to play a vital role in immune regulation and intrinsic tumor characteristics of two immunophenotypically distinct MECs. Further studies are warranted to better understand their role in MECs and their potential use as targets for novel classification and therapeutic implications.

Molecular features associated with immuno-oncology therapeutics

Based on the subgrouping and immunological analyses, we explored potential subgroup-specific molecular features that can be employed in the application of current immunotherapeutic to MECs. First, we examined the efficacy of immune checkpoint inhibitors (ICIs). The analysis of mRNA expression identified significant elevation of PD-1 and PD-L1 levels in the immune-hot MECs (P = 9.6 × 10^− 3 and P = 3.1 × 10^− 2) (Fig. 4A). In addition, co-inhibitor genes, including SLAMF7 (P = 1.0 × 10^− 2), BTN3A1 (P = 1.0 × 10^− 4), and BTN3A2 (P = 9.0 × 10^− 3), were highly expressed in the immune-hot subgroup. Extended analysis of a collection of immuno-oncological therapy biomarkers (43) with the immune-hot subgroup-specific DEGs identified elevated levels of additional immune checkpoint molecules (TIGIT, 4-IBB, TIM-3, PD-L2, and CTLA-4), T-cell targeted immunomodulators (ICOS, IL10, and OX40), adaptive immune systems (MICB), and cell adhesion molecules (ICAM1) in the immune-hot MECs (Fig. 4B and Supplementary Figure S4A). These results suggest that at least a subgroup of immune-hot MECs underwent T-cell exhaustion by checkpoint molecules and could be sensitive to ICIs. In contrast, in the immune-cold MECs, carbonic anhydrase IX (CA9), a surrogate marker of hypoxia-related response(44), and the Wnt signalling signal transduction regulators Dickkopf-1 (DKK1), which functions to confer tumor growth and metastasis(45), Ring Finger Protein 43 (RNF43), cytokine signalling regulator RANK (TNFRSF11A), and cell tight junction-associated claudin 3 (CLDN3) were highly activated (Fig. 4B and Supplementary Figure S4B). This suggests that different or combined therapeutic strategies are needed as the subgroups based on immune activity show distinct immunosuppressive TME characteristics.

In recent years, remarkable success of immunotherapy in cancer treatment has triggered tumor immune microenvironment research in other cancer type (33). With such efforts, understanding of how tumors interact with their surrounding environment has been broadened and utilised in prognosis assessment, treatment planning, and clinical outcome predictions. However, such efforts have been neglected with regard to MECs. Here, we made the first attempt to discover the molecular features that could be used in precision medicine and potential treatment biomarkers of MEC by landscaping its tumor immune microenvironment and analysing the intrinsic characteristics of its immunophenotypic subgroups.

In this study, we demonstrated two immunophenotypically heterogeneous MECs based on the estimated immune and stromal cell-infiltration abundance levels in tumor microenvironment. For immune-hot MEC, elevated immune activity (TIS and IIS score), expression level of the APM score, chemoattractants required for trafficking of T cells to tumors (CXCR9 and CXCR10), cytokines (TNFSF9, IL2, and INFG), and co-stimulators (CD80 and CD28) were observed. Concomitantly, increased expression levels of immune co-inhibitors, such as PD-L1, PD-L2, BTN3A1, and BTN3A2, and immune checkpoint genes (ICOS, OX40, TIM-3, and TIGIT) were observed as well. In the context of the cancer-immune cycle, active tumor-antigen presenting has led to enhanced immune-cell trafficking and T-cell infiltration, resulting in successful anti-tumor immunity, and such dynamic pro-immune TMEs were captured in this group. In contrast, immune-cold MEC was impoverished in immune activation signatures and enriched with genes that are associated with ECM signatures, such as epidermis development, cornification, skin development, and epidermal cell differentiation. Based on these findings, along with downregulated APM expression, we could infer that ineffective antigen-presenting mechanisms and over-expressed physical barriers prevent immune cell infiltration and promote immune evasion (37). Additionally, based on the distinct tumor immune microenvironment profiles observed in the two MEC subgroups, we believe PD-1/PD-L1 blocking immune checkpoint molecules could render therapeutic benefits to immune-hot MEC patients, and that other potential immuno-oncological therapy targets could also be utilised in treating immune-hot and -cold MEC. We hope our findings will aid in future immuno-oncologic therapeutic treatments for immune-hot and -cold MECs and provide wider therapeutic treatment options for MEC patients.

As tumor cells require more energy than normal cells and are the elementary characteristics of every transforming cell, lipid accumulation serves as an energy source and provides cell membrane structural integrity (46). Additionally, high activation of the lipid metabolic state in cancer is characterised by immune suppressive activity in the TME via enhancement of myeloid-derived suppressor cell activity and suppression of dendritic cells and CD1D (38, 39, 47). By conducting a multi-step transcriptomic analysis, we discovered the intrinsic molecular characteristics of the two MEC immune subgroups. It delineated that in immune-hot MEC, lipid metabolic regulator genes, such as MLXIPL (also known as ChREBP), a key regulator of fatty acid synthesis and lipid metabolism, fatty acid synthase (FASN) and stearoyl-CoA desaturase (SCD), were significantly downregulated. As they are known to dysregulate the TME at various levels and promote cancer growth levels (48), but based on our results, we believe that such findings may have contributed to forming elevated immune activity in immune-hot MEC, and further investigation is needed to ascertain the role of dysregulated lipid metabolism in immune-elevated MEC.

The CRTC1-MAML2 fusion has been known to be a good prognostic indicator in MEC (49) and increased tumor-infiltrating immune activity also results in a favourable prognosis. Therefore, we initially anticipated that the CRTC1-MAML2 fusion would be more significant in immune-hot MEC than in -cold MEC; however, there was no significant association between fusion gene status and MEC immune subgroups (Fisher’s exact test, P = 0.37). Based on this result, we must carefully consider a few reasons for this. First, although a higher proportion of fusion-positive patients was present in the immune-hot MEC subgroup (6/8 [75%] vs. 6/12 [50%]), our study sample size was not large enough to show a statistically significant difference between the two groups. Second, in this study, we performed CRTC1-MAML2 neoantigen analysis, but CRTC1-MAML2 fusion did not produce peptides with immunogenicity. This suggests that the CRTC1-MAML2 fusion is not directly related to the immunogenic role or the immunophenotype of MEC. Thus, although we cannot conclude the immunogenic role of CRTC1-MAML2 in MEC at this point, we might be able to do so with a larger cohort.

In this study, we discovered common molecular features that were overexpressed immune-hot and -cold subgroups compared to those in normal MEC immune subgroup samples, namely, POSTN and RUNX2. These are known to be overexpressed in bone development, maturation, and remodelling; epithelial-mesenchymal transition (EMT); ECM structuring and remodelling; and inflammatory processes. Their role as key regulators in such biological processes in other tumors has been well documented (50). We speculate that the roles of POSTN and RUNX2 in MEC are similar to those in previously studied tumors. However, further research is required to comprehensively understand their functional role and potential crosstalk in relation to the tumor immune microenvironment of MEC. It is worth mentioning that attempts have been made to target POSTN for new therapeutic options in other cancers, and we hope that the therapeutic value of POSTN in MEC will also be tested in the future.

In this study, we aimed to demonstrate the heterogeneity of the tumor immune microenvironment of MEC and the intrinsic characteristics of each MEC subgroup and provide potential immunotherapeutic targets. Significantly downregulated lipid metabolic regulators and immune-related genes were among the most distinguishing characteristics of the immune-hot MEC and -cold MEC. We discovered potential immunotherapeutic targets that could be used in immune-inflamed MEC. One limitation of this study is that we were not able to perform prognostic analysis with our study samples due to limited sample size and short follow up period. Although this study suggests promising immune-oncology targets for each immune-based subgroup, this study has limitation in that it cannot derive meaningful results from prognostic analysis due to insufficient sample size and follow-up period. Thus, in future studies, a larger cohort size and longer follow-up period would provide us with more comprehensive knowledge and information on the role of our findings in MEC.

In this study, we demonstrate lipid metabolic activity may play significant role in immunogenic state of MEC, but tumor mutational burden, CRTC1-MAML2 fusion status, and CRTC1-MAML2 fusion neoantigen. Moreover, potential immune-oncology therapeutic targets were proposed for future translational research purposes. Therefore, our findings may help the next step to studying this uncharted cancer.

CNV: Copy number variations

CYT: cytolytic score

DEG: Differentially expressed gene

ECM: Extracellular matrix

EMT: Epithelial-mesenchymal transition

GO: Gene ontology

GTF: Gene transfer format

GSVA: Gene set variation analysis

IM: Immunomodulator

IIS: Immune infiltration score

MAF: Mutation annotation format

MEC: Mucoepidermoid carcinoma

NCCN: National Comprehensive Cancer Network

PPI: Protein-protein interaction

ssGSEA: Single-sample gene set enrichment analysis

TIS: T-cell infiltration score

TMB: Tumor mutational burden

TME: Tumor microenvironment

VEP: Variant effect predictor

VST: variance-stabilizing transformation

Ethics approval and consent to participate

Fresh frozen tissues from primary MECs and their matched normal tissues were retrieved and collected during the operation, after approval by the Institutional Review Board at Severance Hospital, Yonsei University College of Medicine (IRB 255-001).

Consent for publication

Informed consents were obtained from all patients.

Availability of data and materials

The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

Competing interests

The authors declare that they have no competing interests

Funding

This research was supported by a grant of the Korea Health Technology R&D Project through the Korea Health Industry Development Institute (KHIDI), funded by the Ministry of Health and Welfare, Republic of Korea (grant number HI14C1324)

Authors' contributions

HK and MS designed the study, analysed, and interpreted the data. HK, MS, and BP visualized the data. YK and DK collected sample data. DK and SK supervised the study. HK, MS and SK wrote and revised manuscript. All authors read and approved the final manuscript.

Acknowledgements

Not applicable.

Janet-Ofelia GC, Rafael MV, Guillermo GA, Carlos-Enrique CV, Jose-Martin RM, Henry GM, et al. Mucoepidermoid Carcinoma of the Salivary Glands: Survival and Prognostic Factors. J Maxillofac Oral Surg. 2017;16(4):431–7.
Perez-de-Oliveira ME, Wagner VP, Araujo ALD, Martins MD, Santos-Silva AR, Bingle L, et al. Prognostic value of CRTC1-MAML2 translocation in salivary mucoepidermoid carcinoma: Systematic review and meta-analysis. J Oral Pathol Med. 2019.
Jee KJ, Persson M, Heikinheimo K, Passador-Santos F, Aro K, Knuutila S, et al. Genomic profiles and CRTC1-MAML2 fusion distinguish different subtypes of mucoepidermoid carcinoma. Mod Pathol. 2013;26(2):213–22.
Kang H, Tan M, Bishop JA, Jones S, Sausen M, Ha PK, et al. Whole-Exome Sequencing of Salivary Gland Mucoepidermoid Carcinoma. Clin Cancer Res. 2017;23(1):283–8.
Gnjatic S, Bronte V, Brunet LR, Butler MO, Disis ML, Galon J, et al. Identifying baseline immune-related biomarkers to predict clinical outcome of immunotherapy. J Immunother Cancer. 2017;5:44.
Mlecnik B, Bindea G, Angell HK, Maby P, Angelova M, Tougeron D, et al. Integrative Analyses of Colorectal Cancer Show Immunoscore Is a Stronger Predictor of Patient Survival Than Microsatellite Instability. Immunity. 2016;44(3):698–711.
Thorsson V, Gibbs DL, Brown SD, Wolf D, Bortone DS, Ou Yang TH, et al. The Immune Landscape of Cancer. Immunity. 2019;51(2):411–2.
Linxweiler M, Kuo F, Katabi N, Lee M, Nadeem Z, Dalin MG, et al. The immune microenvironment and neoantigen landscape of aggressive salivary gland carcinomas differ by subtype. Clin Cancer Res. 2020.
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21.
Anders S, Pyl PT, Huber W. HTSeq–a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31(2):166–9.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.
Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12(5):453–7.
Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 2016;17(1):218.
Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.
Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics. 2016;32(18):2847–9.
Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.
Senbabaoglu Y, Gejman RS, Winer AG, Liu M, Van Allen EM, de Velasco G, et al. Tumor immune microenvironment characterization in clear cell renal cell carcinoma identifies prognostic and immunotherapeutically relevant messenger RNA signatures. Genome Biol. 2016;17(1):231.
Rooney MS, Shukla SA, Wu CJ, Getz G, Hacohen N. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell. 2015;160(1–2):48–61.
Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7.
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50.
Liberzon A, Birger C, Thorvaldsdottir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 2015;1(6):417–25.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.
Koboldt DC, Chen K, Wylie T, Larson DE, McLellan MD, Mardis ER, et al. VarScan: variant detection in massively parallel sequencing of individual and pooled samples. Bioinformatics. 2009;25(17):2283–5.
McLaren W, Gil L, Hunt SE, Riat HS, Ritchie GR, Thormann A, et al. The Ensembl Variant Effect Predictor. Genome Biol. 2016;17(1):122.
Rubinsteyn A, Hodes I, Kodysh J, Hammerbacher J. Vaxrank: A computational tool for designing personalized cancer vaccines. bioRxiv. 2018:142919.
Wei Z, Zhou C, Zhang Z, Guan M, Zhang C, Liu Z, et al. The Landscape of Tumor Fusion Neoantigens: A Pan-Cancer Analysis. iScience. 2019;21:249–60.
Szolek A, Schubert B, Mohr C, Sturm M, Feldhahn M, Kohlbacher O. OptiType: precision HLA typing from next-generation sequencing data. Bioinformatics. 2014;30(23):3310–6.
Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc. 2016;11(9):1650–67.
Vitting-Seerup K, Sandelin A. IsoformSwitchAnalyzeR: analysis of changes in genome-wide patterns of alternative splicing and its functional consequences. Bioinformatics. 2019;35(21):4469–71.
Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019;47(D1):D607-D13.
Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol. 2014;8(Suppl 4):11.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.
Mandal R, Senbabaoglu Y, Desrichard A, Havel JJ, Dalin MG, Riaz N, et al. The head and neck cancer immune landscape and its immunotherapeutic implications. JCI Insight. 2016;1(17):e89829.
Birkeland AC, Foltin SK, Michmerhuizen NL, Hoesli RC, Rosko AJ, Byrd S, et al. Correlation of Crtc1/3-Maml2 fusion status, grade and survival in mucoepidermoid carcinoma. Oral Oncol. 2017;68:5–8.
Yang M, Lu J, Zhang G, Wang Y, He M, Xu Q, et al. CXCL13 shapes immunoactive tumor microenvironment and enhances the efficacy of PD-1 checkpoint blockade in high-grade serous ovarian cancer. J Immunother Cancer. 2021;9(1).
Cho SY, Kim S, Son MJ, Kim G, Singh P, Kim HN, et al. Dual oxidase 1 and NADPH oxidase 2 exert favorable effects in cervical cancer patients by activating immune response. BMC Cancer. 2019;19(1):1078.
Lu P, Weaver VM, Werb Z. The extracellular matrix: a dynamic niche in cancer progression. J Cell Biol. 2012;196(4):395–406.
Yan D, Adeshakin AO, Xu M, Afolabi LO, Zhang G, Chen YH, et al. Lipid Metabolic Pathways Confer the Immunosuppressive Function of Myeloid-Derived Suppressor Cells in Tumor. Front Immunol. 2019;10:1399.
Jiang L, Fang X, Wang H, Li D, Wang X. Ovarian Cancer-Intrinsic Fatty Acid Synthase Prevents Anti-tumor Immunity by Disrupting Tumor-Infiltrating Dendritic Cells. Front Immunol. 2018;9:2927.
Gonzalez-Gonzalez L, Alonso J. Periostin. A Matricellular Protein With Multiple Functions in Cancer Development and Progression. Front Oncol. 2018;8:225.
Bai Y, Nakamura M, Zhou G, Li Y, Liu Z, Ozaki T, et al. Novel isoforms of periostin expressed in the human thyroid. Jpn Clin Med. 2010;1:13–20.
Kawane T, Qin X, Jiang Q, Miyazaki T, Komori H, Yoshida CA, et al. Runx2 is required for the proliferation of osteoblast progenitors and induces proliferation by regulating Fgfr2 and Fgfr3. Sci Rep. 2018;8(1):13551.
Eddy JA, Thorsson V, Lamb AE, Gibbs DL, Heimann C, Yu JX, et al. CRI iAtlas: an interactive portal for immuno-oncology research. F1000Res. 2020;9:1028.
McDonald PC, Chafe SC, Brown WS, Saberi S, Swayampakula M, Venkateswaran G, et al. Regulation of pH by Carbonic Anhydrase 9 Mediates Survival of Pancreatic Cancer Cells With Activated KRAS in Response to Hypoxia. Gastroenterology. 2019;157(3):823–37.
Haas MS, Kagey MH, Heath H, Schuerpf F, Rottman JB, Newman W. mDKN-01, a Novel Anti-DKK1 mAb, Enhances Innate Immune Responses in the Tumor Microenvironment. Mol Cancer Res. 2021;19(4):717–25.
Baenke F, Peck B, Miess H, Schulze A. Hooked on fat: the role of lipid synthesis in cancer metabolism and tumour development. Dis Model Mech. 2013;6(6):1353–63.
Tiwary S, Berzofsky JA, Terabe M. Altered Lipid Tumor Environment and Its Potential Effects on NKT Cell Function in Tumor Immunity. Front Immunol. 2019;10:2187.
Airley RE, McHugh P, Evans AR, Harris B, Winchester L, Buffa FM, et al. Role of carbohydrate response element-binding protein (ChREBP) in generating an aerobic metabolic phenotype and in breast cancer progression. Br J Cancer. 2014;110(3):715–23.
Adams A, Warner K, Nor JE. Salivary gland cancer stem cells. Oral Oncol. 2013;49(9):845–53.
Maurya R, Bhattacharya P, Dey R, Nakhasi HL. Leptin Functions in Infectious Diseases. Front Immunol. 2018;9:2741.

Table 1. The clinical and pathological characteristics of Mucoepidermoid Carcinoma patient.

Abbreviation: Par, Parotid gland; Sublingual gland, SG; Submandibular gland, SMG; Minor salivary gland, MSG

Download PDF

Version 1

posted

You are reading this latest preprint version

Intrinsic Molecular Characteristics of the Immune Subtypes of Salivary Mucoepidermoid Carcinoma

Status:

Version 1

Abstract

Figures

Background

Materials And Methods

Sample acquisition and preparation

Total RNA sequencing

Data processing

Calculation of TME

Tumor subgrouping with TME estimation

Comparative analysis of subgroups

Variant calling for mutation burden and neoantigen analysis

Analysis of isoform switching

Network-based functional analysis of genes

Results

Discussion

Conclusions

Abbreviations

Declarations

References

Tables

Supplementary Files

Status:

Version 1