Genomic and Transcriptomic Landscape of an Oral Squamous Cell Carcinoma Mouse Model for Immunotherapy

Head and neck cancer is among the most prevalent cancer types worldwide, with an annual incidence of over 500,000 (1). Patients with head and neck cancer are usually treated with surgery, radiotherapy, chemotherapy, or a combination of these treatments, and the corresponding 5-year survival rate is 63% (2). The FDA-approved the use of immune checkpoint inhibitors (ICI) in 2016, and thus expanded treatment options for recurrent and metastatic head and neck squamous cell carcinoma (HNSCC). ICIs targeting programmed death-1 (PD-1) have been reported to exhibit notable efficacy in HNSCC (3–5). However, only a small proportion of patients (<20%) benefit from single-agent immunotherapy. Therefore, developing other treatment options to improve the effectiveness of ICIs in most patients with HNSCC is imperative (5).

Head and neck cancer development represents a sequential pathogenic process that involves a progression from inflammation to invasive cancer (6). Studies examining the genetic heterogeneity of head and neck cancer in humans have primarily focused on genetic lesions in flank tumors. Recent whole-genome studies have demonstrated that most tumors have inactivating mutations in TP53, and mutations in NOTCH1, FAT1, and PIK3CA have also been frequently observed (7, 8).

4-nitroquinoline-1-oxide (4NQO) is an ideal carcinogen for the development of an experimental oral carcinogenesis model (9, 10), and the pathologic processes induced by 4NQO are similar to those observed in human oral squamous cell carcinoma (OSCC). In this study, we used an autologous mouse model of OSCC induced by the oral administration of 4NQO. Cell exposure to 4NQO resulted in DNA alterations, a process similar to tobacco-induced mutagenesis (11), and DNA alterations in Trp53, Notch1, Fat1, and Pik3ca were also observed in mouse tumors. In addition, the 4NQO-induced mouse model, which exhibits a fully functioning immune system, has been extensively used in immunotherapy studies (11–13). However, these studies have focused either on the DNA alterations that occur at the transition from oral premalignant lesions to OSCC or on the immune landscape of lymph nodes using flow cytometry analysis. In the current study, we focused on the primary 4NQO-induced tumors on tongues and covered all developmental stages. The genomic and RNA profiles generated in this study offer a comprehensive picture of the genomic alterations and immune infiltration during tumor development and can be used for testing drug responses to ICI-based therapies.

CD40 plays a secondary and stimulatory role in the activation of multiple immune cells (14). Agonistic CD40-targeting antibodies have been shown to inhibit tumor progression in both the preclinical and clinical settings using several tumor models, including a head and neck cancer model (15–17). The immune modulatory effect of these antibodies is based on the broad pro-inflammatory effects of CD40, its ligand on dendritic cells (DC) and macrophages, and downstream B-cell and T-cell activation (18, 19). Single-agent anti-CD40 therapy has been demonstrated to have moderate antitumor activity, and the most common adverse events of this treatment include cytokine release syndrome, anemia, lymphopenia, and liver function abnormalities (20–22). This study additionally examined the effectiveness of agonistic anti-CD40 therapy in the 4NQO-induced model and propose that alterations in the M1/M2 macrophage ratio correlate with the antitumor activity of the therapy. The findings of this study provide critical insights into the effectiveness of ICI-based treatment strategies in inhibiting the occurrence of HNSCC.

Six-week-old male C57BL/6J mice were purchased from the National Laboratory Animal Center (Taipei, Taiwan). The mice were housed under specific pathogen-free conditions at the Laboratory Animal Center of the College of Medicine, National Taiwan University. Domestication continued for 1 week before the start of each experiment. All animal experiments in this study were conducted in accordance with approved protocols [Institutional Animal Care and Use Committee (IACUC), Approval No: 20200055] and guidelines specified by the IACUC of the College of Medicine, National Taiwan University, and conformed to the criteria outlined in the Guide for the Care and Use of Laboratory Animals.

A batch of 6-week-old male C57BL/6J mice were administered with regular water containing or not containing 100 μg/mL of 4NQO (Sigma-Aldrich) for 10 weeks and subsequently fed with regular water for 19 weeks. The time at which mice were first administered 4NQO was week 0. A total of 40 4NQO-treated mice and 4 normal mice were euthanized during the study period, with 8 to 12 mice being euthanized at each of the following weeks: 0, 10, 16, 22, and 29 (Supplementary Fig. S1A). For anti-CD40 treatment, 4NQO-fed mice (n = 30) were randomly divided into 4 groups and were treated at week 3 post-4NQO feeding: normal (negative control, receiving regular water), control (isotype antibody), anti–PD-1, and agonistic anti-CD40. Mice were treated a total of 5 times (once every 3 days) via intraperitoneal injection and sacrificed at week 17. Antibodies used included isotype control antibody (InVivoPlus rat IgG2a isotype control, anti-trinitrophenol; Bio X Cell), rat anti–PD-1 (anti-CD279, InVivoPlus, Bio X Cell), and agonistic anti-CD40 (clone FGK4.5, InVivoMab, Bio X Cell). Harvested tongues were cut in half with scissors, and one tongue section was first fixed with 14% formalin (Choneye Pure Chemicals) and stained with hematoxylin and eosin (H&E) for histopathologic examination, and the other section was frozen. Not all 4NQO-bearing mice were used in the whole-exome sequencing (WES) and RNA sequencing (RNA-seq). On the basis of the pathologic stage, part of the frozen tongues was selected and used for subsequent extraction of DNA/RNA for WES and bulk RNA-seq (Supplementary Fig. S1A).

The lesion regions (yellow half circle, Supplementary Fig. S1A) on the frozen tongues from the selected 4NQO-induced OSCC mice were dissected. Genomic DNA and total RNA were extracted simultaneously from the same sample by using the AllPrep DNA/RNA Mini Kit (catalog no. 80204, Qiagen) according to the manufacturer's protocol. 20 mg of harvested tumors were disassociated in RLT Plus buffer (kit reagent) with 1% of 2-mercaptoethanol (Sigma-Aldrich), and the tumor was homogenized by the soft tissue homogenizing CK14 tube (P000933-LYSK0-A, Bertin Corp) through repeated centrifugation (6,000 g, 20 sec, 3 times repeat, Tomy MX-301 highspeed centrifuge). The supernatant in the column was transferred to the DNA-containing column, and RNA remained in the collection tube. For DNA, AW1 and AW2 buffer were added in order (both kit reagents) following high-speed centrifugation (10,000 rpm). DNA was purified by adding of EB buffer (kit reagent). For RNA, 70% ethanol (Sigma-Aldrich), RW1 buffer (kit reagent), and RPE buffer (kit reagent) in order were used for washing. RNA was resuspended in RNase-free water (Qiagen). The concentration of DNA and RNA was quantified using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific).

Genomic DNA was isolated from the tongues of selected 4NQO-treated mice with three pathologic stages (n = 19; Supplementary Table S2) and the tongue of a normal mouse (n = 1) by using the AllPrep DNA/RNA Mini Kit (Qiagen). The WES of all DNA samples was performed by the Genomic BioSci & Tech Co. Briefly, exome capture libraries were prepared using the SureSelectXT Mouse All Exon Kit (Agilent Tech) by the manufacturer's instructions; the prepared libraries in 1 μg of DNA input were subsequently sequenced using the Illumina NovaSeq 6000 system operated under a 150-base paired-end read configuration. Adaptor sequences were trimmed, and low-quality reads were removed using Cutadapt (v3.5) with parameters –minimum-length = 50 –trim-n –max-n 5. The reads were then aligned to the mouse genome assembly (GRCm38.p6) by using BWA (v0.7.17). Polymerase chain reaction duplicates were marked using Picard (v2.25.4), and base quality score recalibration was performed using GATK (v4.2.3.0). Single-nucleotide variants (SNV) and small indels were called using Mutect2 implemented in GATK (v4.2.3.0) and annotated using ANNOVAR (v20210819). The COSMIC mutational signatures (23) determined for each sample were fitted using the R package deconstructSigs (24).

Total RNA of the tumors from the tongues of selected 4NQO-treated mice (n = 28; Supplementary Table S2) and the tongues of normal mice (n = 3) was simultaneously isolated with genomic DNA extraction by using the AllPrep DNA/RNA Mini Kit (Qiagen) as described above. The RNA integrity number (RIN) of each sample ranged between 8.4 and 9.9 (RNA 6000 Nano kit, Agilent 2100 Bioanalyzer, Agilent Tech). The RNA-seq of total RNA samples was assessed by the Genomic BioSci & Tech Co. Briefly, 2 μg of RNA input per sample in TruSeq stranded mRNA mouse libraries were prepared using the TruSeq stranded mRNA sample preparation kit and subsequently sequenced using the Illumina NovaSeq 6000 system operated under a 150-base paired-end read configuration. Adaptor sequences were trimmed using Cutadapt (v3.5) with parameters –minimum-length = 50 –trim-n –max-n 5. Clean reads were aligned to the mouse genome assembly (GRCm38.p6) by using STAR (v 2.7.8a) with two-pass mode. A gene count matrix was generated using featureCounts (v2.0.1) based on the gene annotation of GENCODE vM25. The trimmed mean of M values (TMM) normalization method was used to adjust for differences in library size between the samples. The immune cell composition of each sample was deconvoluted using the CIBERSROTx (25) with the mouse immune cell matrix (26). The sample-level gene set score was calculated using the single-sample gene set enrichment analysis (GSEA) implemented in the R package gsva for gene sets obtained from MSigDB (v7.4).

The trajectory of tumor progression was considered to be as follows: normal (NOR) → hyperplasia (HYP) → mild dysplasia (DYS) → moderate to severe dysplasia (CIS) → invasive carcinoma (ICA). For each gene, differential expression analysis was performed for each pairwise comparison (i.e., normal vs. HYP; HYP vs. DYS; DYS vs. CIS; and CIS vs. ICA) along the trajectory by using the R package limma. An adjusted P value of < 0.05 was considered to indicate differential expression. If a significant increase in expression was noted in a specific pairwise comparison, the change in expression for that comparison was scored as 1; if a significant decrease in expression was noted in a specific pairwise comparison, the change in expression for that comparison was scored as −1. If no significant change was noted, a score of 0 was assigned. The relative gene expression was determined by examining all genes at a baseline of 0. Subsequently, these scores were added across pairwise transitions. Genes with a score of 0 at all stages were removed, and the remaining genes were included in gene clusters based on their transition patterns. An overrepresentation analysis was performed for each gene cluster using gene sets obtained from MSigDB (v7.4) using the R package clusterProfile.

The tumor regions of interest (ROI) were circled on 10-μm-thick section slides of formalin-fixed, paraffin-embedded (FFPE; normal n = 1, vehicle n = 5, agonistic CD40 n = 3) and total RNA was extracted for NanoString assay by Cold Spring BioTech Co. NanoString raw data were processed using the R package NanoStringQCPro. The library size factor for each sample was calculated using the TMM method based on the expression of housekeeping genes (including Rpl19, Abcf1, G6pdx, Oaz1, Eef1g, Ppia, Eif2b4, Polr1b, Edc3, Hprt, Tubb5, Sf3a3, Polr2a, Tbp, Nubp1, Hdac3, Alas1, Sap130, Gusb, and Sdha), which was used to normalize gene counts across all samples. Preranked GSEA was performed using gene sets from MSigDB (v7.4) in the R package clusterProfiler. Genes were ranked on the basis of the score: −log₁₀ (P value) × sign (log₂FC); the corresponding P value and log₂FC were derived from the assessment of differential expression between samples with and without anti-CD40 treatment by using the R package limma.

H&E staining was performed at the IACUC of the College of Medicine, National Taiwan University (Taipei, Taiwan). FFPE tumor tissues were processed for IHC staining with the following antibodies: anti-keratin 17 (1:200, Cell Signaling Technology) and anti-ki67 (1:500, Proteintech). Antigen retrieval was performed by immersing the tissue sections in 0.01 mol/L citric acid buffer solution (pH 6.0), followed by autoclaving at 121°C for 10 minutes to enhance detection. After antigen retrieval and a specific cooling-off period, endogenous peroxidase was inactivated through the addition of 3% hydrogen peroxide to PBS for 10 minutes at room temperature. The antibodies were incubated with the tissue sections at 4°C overnight in a humidified chamber. After incubation with the primary antibodies, alkaline phosphatase–based IHC staining was performed using a detection kit (EnVision Detection System, Peroxidase/DAB, Rabbit/Mouse, Agilent Tech, Santa Clara) in accordance with the manufacturer's instructions. The whole specimen section slides were examined and acquired by the 20× lens TissueFAXS High-Speed Multispectral Tissue Cytometer (TissueGnostics). Images analysis of ki67 was processed using StrataQuest software (TissueGnostics). Briefly, cell numbers were counted by recognizing the nuclear signal first, and subsequently the intensity of the target marker was normalized to the negative controls (without primary antibody). All the selected ROIs (0.5 × 0.5 mm) on slides were processed by the same analysis setting. IHC stain counts shown on the dot plot mean total counts of every independent region (ten regions per treatment group). All statistical tests were performed in GraphPad Prism version 8.0 (GraphPad software).

Multiplex immunofluorescence staining was performed using the Opal 7-Color Automation IHC Kit (NEL861001KT, Akoya Biosciences) in accordance with the manufacturer's protocol. FFPE tissues with and without primary antibody treatment were used as negative controls in each staining process. After deparaffinization, rehydration, and endogenous peroxidase blocking, the sections were washed with Tris-buffered saline with 0.1% of Tween 20 (TBST; chemicals are all from Sigma-Aldrich) and blocked with a 1x Antibodydiluent/blocking buffer (ARD1001EA, Akoya Biosciences) for 10 minutes; subsequently, samples were incubated with the primary antibody overnight. Each of the sections was then incubated with the Opal anti-Ms + Rb horseradish peroxidase (ARH1001EA, Akoya Biosciences) for 10 minutes, followed by incubation with an Opal fluorophore (Opal520, Opal540, Opal570, Opal620, Opal650, or Opal690) for 10 minutes. Regarding the Opal dyes used to label the corresponding primary antibodies, Opal620 was used to label rabbit anti-CD4 (Abcam), Opal570 was used to label rabbit anti-CD8 (Bioss Antibodies Inc), Opal650 was used to label rabbit anti-MHC II (Bioss), Opal520 was used to label rabbit anti-CD19 (HistoSure, Synaptic Systems GmbH), Opal520 was used to label rabbit anti-CD206 (Bioss), and Opal690 was used to label rabbit anti-F4/80 (HistoSure). Bound primary and secondary antibodies were removed through repetitive heat-mediated antibody-stripping treatments. After each section was washed with cold water and TBST, the staining and antibody removal processes were repeated using a different Opal fluorophore. Finally, after staining with the sixth Opal fluorophore, the tissue specimens were stained with 1X spectral 4′,6-diamidino-2-phenylindole (DAPI; FP1490, Akoya Biosciences) for 5 minutes and mounted in ProLong Diamond Antifade Mountant (Thermo Fisher Scientific). Images were acquired using the Vectra Polaris automated quantitative pathology image system (Akoya Biosciences) in accordance with the manufacturer's instructions. Data analysis was performed using inForm tissue analysis software version 2.4.6 (Akoya Biosciences). Each antibody staining was combined using spectral library information (Opal 7-Color Kit) to associate each fluorochrome component with an individual staining component. All cell populations (blank, normal tissue, and tumor cells) from each panel were characterized and quantified using the cell segmentation and phenotype cell tool implemented in inForm Tissue analysis software was performed under a pathologist's supervision. The number of each marker was quantified and normalized by the area of each FFPE section slide (density/mm²). Subsequently, quantified expression data of each immune marker were imported in GraphPad Prism software for drawing dot plots and examining the statistical difference.

The raw data of RNA-seq and NanoString generated during this study are available via Gene Expression Omnibus with accession numbers GSE229289 and GSE217407. The raw data from WES are available via European Nucleotide Archive with accession number PRJEB61274.

To evaluate genetic alterations and gene expression during OSCC development, we used a carcinogen-induced mouse model and collected tongues at different stages (Fig. 1A; Supplementary Fig. S1A). The tongue of the 4NQO-fed mice exhibited a rough surface at week 10. Subsequently, visible white patches and surface toughness began to appear at week 16 and worsened at weeks 22 and 29 (Fig. 1B). All the tissue irregularities were classified as tongue lesions and progressed in a time-dependent manner (Fig. 1B). Tongue samples were collected from the 4NQO-treated mice and examined at different time points (weeks 10, 16, 22, and 29) for pathologic stage detection through histologic staining (Fig. 1C). H&E stains revealed five pathologic stages: NOR, HYP, DYS, CIS, and ICA (Fig. 1C; Supplementary Fig. S1B). The NOR stage was characterized by a well-oriented stratified epithelium with palisading basal cells. In the HYP stage, there was an increase in the number of cells. The dysplasia stages showed losses of epithelial cell polarity, with nuclear pleomorphism and increased keratinization (Fig. 1C). The ICA stage demonstrated squamous cell carcinoma with stromal invasion (Fig. 1C). The occurrence of carcinogenesis associated with the 4NQO treatment period (Fig. 1D). The observed histologic changes, Ki67 overexpression (Fig. 1E), and increased keratin 17 expression indicated that OSCC development was a linear process. To investigate the carcinogenesis process, we collected tongue samples from the mice and analyzed them through WES and RNA-seq.

We investigated genetic changes in the 4NQO-induced OSCC tongues compared with the normal tissues using WES. We identified 27,234 SNVs, 221 small indels, and 155 doublet base substitutions in 19 OSCC samples. Overall, the mutational burden was higher in the CIS and ICA stages than in the HYP and DYS stages (Fig. 2A). All categories of somatic single-base substitutions were identified in the 4NQO-induced neoplasms, and the C > A mutation dominated the spectra (Fig. 2B and C); this finding is consistent with that obtained for a previous 4NQO-induced oral cancer mouse model (11). We conducted a principal component analysis (PCA) of the mutational spectra, and the analysis revealed the trajectory of the pathologic stages (Fig. 2D), indicating the relative abundance of different mutations with the progression of lesions. The abundance of the C > A mutation increased significantly with the progression of lesions (Supplementary Fig. S2A).

Mutational signature analysis of individual tumors can reveal the etiology of mutational processes. We determined the mutational signatures of each sample based on COSMIC mutational signatures (version 2). The compositions of the mutational signatures in all samples were similar (average cosine similarity > 0.86), except for that in one HYP sample. However, a slight difference was observed between precancerous (HYP and DYS stages) and cancerous lesions (CIS and ICA stages; Fig. 2E; Supplementary Fig. S2B). Most of the signatures comprised COSMIC Signatures 4, 6, 14, 18, 20, 24, and 29 (Fig. 2E; Supplementary Fig. S2C). Regarding the etiologies of the signatures, Signature 4, which constituted the largest component in the 4NQO-induced lesions, associates with tobacco smoking in HNSCC (27). Moreover, Signatures 6 and 20 associate with defective DNA repair (28), and Signature 24 associates with exposure to aflatoxin in OSCC (29). However, the etiologies of Signatures 14 and 18 remain unknown, but these signatures are commonly identified in human HNSCC (30). We examined the contributions of these signatures across the pathologic stages (Supplementary Fig. S2D) and determined that the contribution of Signature 4 increased with the progression of lesions (Fig. 2F). However, Signature 6 negatively correlated with the pathologic stages (Fig. 2F).

Finally, we analyzed all genes that were mutated in our samples. A total of 9,611 genes were mutated in at least one of the samples (Supplementary Table S1). Several of the mutated genes identified in this study are frequently identified in human HNSCC, including TP53, CSMD3, NOTCH1, and SYNE1 (Fig. 2G). We compared the frequency of mutated genes in our model with that in another 4NQO-induced mouse OSCC model (31). We noted a significant correlation between the two models (Supplementary Fig. S2E), demonstrating the robustness of the 4NQO-induced OSCC model. These results suggest that the 4NQO-induced OSCC model recapitulates human HNSCC. Studies examining changes in RNA expression in each 4NQO-induced OSCC stage may provide insights into the molecular mechanism underlying the carcinogenesis of HNSCC.

We analyzed the gene expressions and performed a PCA of the whole transcriptome of all OSCC samples and found samples clustered on the basis of the pathologic stages. The lesion samples formed a gradient across the principal component (PC) space, organizing a cancer progression trajectory (Fig. 3A). We investigated genes whose expression changed along the progression trajectory (NOR → HYP → DYS → CIS → ICA) and clustered genes by the relative expression change between stages. A total of 70 expression patterns that exhibited at least one transition change were identified, and 79.1% of the genes were captured in the first 12 clusters (Fig. 3B; Supplementary Fig. S3). We conducted an overrepresentation analysis on the genes in each cluster by using hallmark gene sets obtained from MSigDB. The gene clusters correlated with various biological processes and were dependent on the pathologic stage (Fig. 3C). Genes in cluster 2, involved in the epithelial–mesenchymal transition, exhibited an inflection point between DYS and CIS, suggesting that genes in cluster 2 are involved in tumor invasiveness. In addition, the scores of gene sets related to tumor invasiveness and stemness increased in the CIS and ICA stages (Fig. 3D). A critical transitional point of tumor progression might be triggered between DYS and CIS stages. Although clusters 8 and 9 contained genes involved in responses to inflammation and IFNγ, they exhibited distinct expression patterns: the expression of cluster 8 increased suddenly in the HYP stage, whereas that of cluster 9 increased continually (Fig. 3B). Furthermore, we examined the IFNγ-response signatures and found that the expression of Ifng was high only in the HYP stage; its expression was suppressed after the DYS stage (Fig. 3D and E, middle and left, respectively). However, the expression of immune response–related genes was increased in the CIS and ICA stages (Fig. 3D, bottom). Although an increase in inflammatory responses was noted, the expression of genes related to immune suppression and T-cell exhaustion was also high in invasive tumors, including Lag3 and Pdcd1 (Fig. 3E; Supplementary Fig. S4). According to the stage-dependent expression of genes related to IFNγ responses, inflammatory responses, and immune suppression, an abnormal balance of immune responses may occur during the carcinogenesis process.

The dynamics of immune cells play a critical role in tumor progression. Hence, we examined the immune microenvironment of individual samples by using CIBERSORTx (Fig. 4A). PCA was conducted to examine immune cell composition, revealing that samples were aligned along the trajectory of tumor progression (Fig. 4B). The number of tumor-infiltrating lymphocytes (TIL), quantified from CIBERSORTx, increased during tumor progression (Fig. 4C). Similarly, the gene expression of Cd45 was higher in invasive lesions than in the HYP and DYS stages (Fig. 4C). We subsequently examined immune cell abundance across the pathologic stages. The number of immune cells differed significantly across the pathologic stages (P < 0.05, ANOVA; Fig. 4D and E; Supplementary Fig. S5). The numbers of DCs and CD8⁺ T cells increased with the advancement of the pathologic stage (Fig. 4E). However, M1 macrophages started to increase in the HYP stage but decreased after the DYS stage, and CD4⁺ T cells decreased in the carcinoma stage (Fig. 4E). Although the number of CD8⁺ T cells increased, these cells might be dysfunctional or exhausted owing to the high expression of exhaustion markers in the CIS stage, including Pdcd1, Lag3, and Havcr2 (Supplementary Fig. S4). In our 4NQO-induced OSCC samples, the number of M1 macrophages decreased, whereas those of M0 macrophages and M2 macrophages increased (Fig. 4F). Although the abundance of M2 macrophages did not differ significantly across the pathologic stages (Fig. 4D), the expression of the M2 signature increased (Fig. 4F).

Opal multiplex immunofluorescence staining using immune system–based markers serve as a robust tool for immune profiling. It enables the simultaneous detection of multiple markers in the same tissue section in FFPE tumor tissues. We used 5-μm-thick FFPE sections and detected signals for B cells, CD4⁺ and CD8⁺ cells, total macrophages, M1 macrophages, and M2 macrophages. M2 macrophages substantially accumulated, indicating their leading role in the OSCC tumor microenvironment (TME; Fig. 5A and B). The expression of the M1 macrophage marker was relatively weak but increased in the HYP stage. The numbers of B cells, CD4⁺ T cells, and CD8⁺ T cells increased in the carcinoma stages (Fig. 5C). The expression pattern of immune markers observed in multiplex immunofluorescence staining (Fig. 5A; Supplementary Fig. S6) was similar to that observed in our RNA-seq analysis (Fig. 4D and E). We further analyzed the inForm-processed opal data to examine the distance between M1 or M2 macrophages and CD8⁺ T cells in the HYP and CIS stages. The dot plot revealed no significant differences between the HYP and CIS stages regarding the distances between either the M1 macrophage and CD8⁺ T cells or the M2 macrophages and CD8⁺ T cells (Supplementary Fig. S7).

We expected that the 4NQO-induced OSCC mouse model would serve as a reliable screening platform for preclinical testing in HNSCC. Accordingly, we examined whether the effect of ICI treatment could be evaluated using the 4NQO-induced tumor model. We tested the effects of two ICIs, namely anti–PD-1 and agonistic anti-CD40 (Fig. 6A). We treated the 4NQO-fed mice with an isotype-control antibody or agonistic anti-CD40 or anti–PD-1 5 times (once every 3 other day) from week 3 after 4NQO exposure. We examined changes in the morphology of the tongue lesions in accordance with the aforementioned definition of the morphologic changes (Fig. 1D). At week 17 after 4NQO treatment, the control tongue samples exhibited several lesions and severe keratinized tissues growing outward from the surface. Some over-differentiated tumors were larger in the anti–PD-1–treated tongues than in the control tongues (Fig. 6B). The lesions were significantly smaller in the agonistic anti-CD40 treatment group than in the control and anti–PD-1 groups (Fig. 6B). Tongues exhibiting a whitish and papillary appearance were subjected to H&E staining for histopathologic analysis (Fig. 6C). The mucosa was thicker in the agonistic anti-CD40 treatment group than in the healthy group, but the structures remained normal (Fig. 6C). Slight neoplasia or hyperplasia was observed in some parts of the tongues. Overall, 4NQO-induced cancerous lesions were reduced in the agonistic anti-CD40 treatment group relative to those in the vehicle and anti–PD-1 treatment groups.

We examined immune filtration in the 4NQO-induced primary tumors. We observed the accumulation of B cells, CD4⁺ T cells, or CD8⁺ T cells in the agonistic anti-CD40–treated tumors (Fig. 6D and E; Supplementary Fig. S8A). Abundant M2 macrophages were observed around the mucosa and muscle regions in the control tumors (Fig. 6D), and the abundance of M2 macrophages decreased in the agonistic anti-CD40–treated primary tumors (Fig. 6D). The number of M1 macrophages was relatively low in the control tumors and increased significantly in the agonistic anti-CD40–treated tumors (Fig. 6F). The number of M2 macrophages was significantly decreased in the agonistic anti-CD40–treated tumors (Fig. 6F). Therefore, the ratio of M1/M2 macrophages increased in the agonistic anti-CD40–treated tumors (Fig. 6F).

To investigate functional changes in the immune microenvironment after agonistic anti-CD40 treatment, we profiled tumors by using the NanoString Pan-cancer Immunology Panel. The GSEA results for these tumors revealed that genes and signatures related to antitumor immunity were highly expressed in the agonistic anti-CD40–treated tumors (Fig. 6G). B cells are the major target cells of agonistic anti-CD40, and we found an increase in the expression of genes related to B-cell activation and humoral immunity in the agonistic anti-CD40–treated tumors relative to that in the control tumors (Fig. 6G, top). In addition, agonistic anti-CD40 treatment induced the expression of genes involved in macrophage activation and increased the expression of signatures related to the M1/M2 ratio (Fig. 6G; Supplementary Fig. S8B, bottom). Signatures related to complement activity were highly expressed in the agonistic anti-CD40–treated primary tumors than in the control tumors (Supplementary Fig. S8B).

We identified 41 differentially expressed genes by comparing tumors that received and did not receive agonistic anti-CD40 treatment (Fig. 6H). Some of these genes were related to IFNγ responses, inflammatory responses, and macrophage activation. The data indicate that agonistic anti-CD40 signaling may enhance immune responses by modulating the expression of these selected genes.

Because human and mouse models are genetically heterogenous, our OSCC mouse model can reflect part of human OSCC TME. In this study, changes in immune cell composition in each of the pathologic stages in the primary TME were elucidated. When 4NQO was initially administered, the morphology of mucosae remained normal, the M1 macrophages were activated, and the antitumor immune cells defended against the external carcinogen. With persistent exposure to the carcinogen, M2 macrophages, CD8⁺ T cells, and other immune cells were activated, potentially influencing the outcome of OSCC. In the early stages of carcinogenesis following exposure to 4NQO, the agonistic anti-CD40 treatment was effective in inhibiting tumor growth and inducing specific alterations in the compositions of immune cells, resulting in the creation of an environment hostile toward tumor development. Our OSCC tumor model with a clear immune signature fingerprint and specific genomic mutations can thus be used to examine the effects of a variety of ICIs or combination treatments.

Several studies have analyzed the genomic mutation profiles in 4NQO-induced mouse models, revealing high heterogeneity between human and mouse samples (11, 32). Wang and colleagues used a syngeneic animal model to describe several differences in DNA mutation and immune landscapes between two different OSCC-derived cell lines. Our WES data from mouse tumors revealed several DNA mutations that may not be found in large human databases. Consequently, immune responses occurring in the early pathologic stages, including inflammation, immune cell recruitment, and immune effector activation, may be ignored. Low ryanodine receptor 2 (Ryr2) expression associates with poor prognosis in thyroid carcinoma (33). Somatic mutations in Ryr2 have been frequently observed in head and neck cancer (34). Fat3, a tumor suppressor, was mutated from the dysplasia to the carcinoma stages. Fat3 mutation was previously reported in various solid tumors and might lead to immune dysfunction (35). Our WES data indicated increased inflammatory activation in the TME. The recruitment or expansion of inflammatory immune cells can lead to their activation, polarization, and differentiation to rapidly defend and kill tumor cells in the TME (36, 37).

We observed that activation of some immune cells in the primary TME was dependent on pathologic stage. The composition of the immune signature changed dynamically, with a transition point occurring between the HYP and DYS stages, supporting the idea that the immune signature varies at different stages and correlates with the progression of the oral premalignant lesions into OSCC (38). Most immune cells exhibited stage-dependent activation. Macrophages regulate innate and adaptive immune responses by recruiting other leukocytes such as CD8⁺ T lymphocytes and natural killer cells. Whether tumor-associated macrophages (TAM) in the primary TME exhibit protumor or antitumor function depends on the ratio of M1 and M2 macrophages. From our data, pro-inflammatory M1 macrophages were activated in the HYP stage but rapidly decreased in the next stages. M0 macrophages and anti-inflammatory M2 macrophages were activated from the dysplasia to the carcinoma stages. M1 macrophages produce antitumor chemokines, inducible nitric oxide synthase, interleukin-1β, and TNFα to rapidly kill tumors (39, 40). Previous studies have reported that in patients with HNSCC, M2 macrophages are abundant, and a high infiltration of CD68⁺ TAMs associates with a poor prognosis and shorter disease-free survival time (41, 42). Moreover, other studies have revealed that an increase in CD68⁺ TAMs positively correlates with a tumor size and metastasis (43, 44). The depletion of DUSP1 inhibits tumor growth in a subcutaneous HNSCC-bearing mice tumor model (45). Ex vivo, an increase in the pro-inflammatory cytokine IL1β is observed in bone marrow–derived macrophages with Dusp1^−/− genetic background.

The efficacy of agonistic anti-CD40 treatment may result from its ability to modulate the immune composition. Agonistic anti-CD40 treatment activated antitumor activity by increasing the number of M1 macrophages, suggesting that agonistic anti-CD40 may alter macrophage polarization in the TME. Our GSEA data revealed that complement and humoral immunity, lymphocytes, and B cells were involved in the effect of agonistic anti-CD40. Complement-associated immunity and humoral immunity are related to the activation of B cells, suggesting that B cell–mediated immune responses are sensitive to agonistic anti-CD40 treatment. Furthermore, we examined changes in the expression of specific genes owing to agonistic anti-CD40 treatment and noted changes in the expression of 14 genes. Signaling via Toll-like receptors (TLR) can alter the TME, and TLR1, TLR2, TLR7, and TLR8 are reported to be associated with macrophage activity (46). Survival analysis reveals that TLR7 expression significantly affects survival, and lower TLR7 expression associates with a poorer prognosis in The Cancer Genome Atlas database (47). Increased TLR7 expression leads to an increase in the M1/M2 ratio in TAMs (48). In addition, CD40 activation induces macrophage polarization toward the M1 phenotype (49). These findings suggest that agonistic anti-CD40 treatment upregulates TLR7 expression, thereby delaying tumor progression. Other methods of modulating macrophage polarization and activation, such as targeting M1 or M2 macrophages specifically, are also worthy of further exploration.

ICIs are shown to exhibit antitumor activity in a 4NQO-induced mouse model that focused on established tumors (49). We first observed immune signature changes during early carcinogenesis. We noticed an increase in the number of CD8⁺ T cells and M1 macrophages in primary tumors using RNA-seq and multiplex immunofluorescence staining. On the basis of our preliminary drug treatment tests, including the use of anti–PD-1 and agonistic anti-CD40, we observed that drug treatment at the early stages of carcinogenesis can effectively block the development of OSCC and prevent the formation of oral premalignant lesions in mice. In the current study, in contrast to the effect of anti–PD-1, CD40/CD40 ligand mediate the immune response by modulating specific immune cells. Preclinical studies have demonstrated that CD40-mediated gene activation can exert an antineoplastic effect in several tumor models. One study reveals that anti-CD40 treatment causes increases in the numbers of CD8⁺ T cells and M1 macrophages, whereas anti–PD-1 and anti–CTLA-4 therapies lead to depletion of regulatory T cells (49). Another study indicates that CD40 activation induces the polarization of macrophages from the anti-inflammatory M2 phenotype into the pro-inflammatory M1 functional phenotype, resulting in enhanced antitumor efficacy (50). Further, a clinical study reveals that in patients with advanced solid tumors, treatment with bispecific anti-FAPxCD40 MP0317 (NCT05098405) reduces the M1/M2 macrophage ratio in vitro (51).

In the current study, composition changes in immune cells in the TME at various pathologic stages were described. In the early stages of carcinogenesis following 4NQO exposure, the mucosae's structure remains unaffected, triggering activation of M1 macrophages and induction of antitumor responses. With prolonged exposure to the carcinogen, M2 macrophages, CD8⁺ T cells, and other immune cells are activated, potentially impacting the development of OSCC. In conclusion, our OSCC tumor model, with a clear immune signature fingerprint and specific genomic mutations, can be used to examine the effects of a variety of ICIs or combination treatments.

No disclosures were reported.

Y.-M. Lee: Conceptualization, data curation, methodology, writing–original draft, writing–review and editing. C.-L. Hsu: Conceptualization, data curation, methodology, writing–original draft, writing–review and editing. Y.-H. Chen: Data curation, formal analysis, validation, methodology, writing–review and editing. D.-L. Ou: Supervision, visualization, writing–review and editing. C. Hsu: Supervision, visualization, writing–review and editing. C.-T. Tan: Conceptualization, supervision, funding acquisition, visualization, project administration, writing–review and editing.

The work was supported by the High-throughput Genomics & Dig Data Analysis Core, Department of Medical Research, National Taiwan University Hospital for technical support and National Center for High-performance Computing (NCHC) for providing computational and storage resources. We thank the staff of the Second Core Lab, Department of Medical Research, National Taiwan University Hospital for microscope technical support in IHC staining. This manuscript was edited by Wallace Academic Editing.

This work was supported by the National Taiwan University Hospital under grant MS251; the National Health Research Institutes of Taiwan and National Taiwan University Hospital Hsin-Chu Branch under grant NHRI-111-B06; the Hsin-Chu Science Park under grant B11102; and the Ministry of Science and Technology of Taiwan under grant MOST 109–2314-B-002–057-MY3.

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.