TP53 mutations predict poor response to immunotherapy in patients with metastatic solid tumors

TP53 is the most commonly mutated gene across all cancer types. R175H mutation was considered structural mutation where the mutation causes misfolding of the protein and leads to a significant conformational alterations within p53's DNA binding domain. The aim of this study was to explain the reason why R175H worse the response to immunotherapy by analyzing tumor immune microenvironment through the expression of immune cells and PD‐1.


| BACKGROUND
TP53 is a transcription factor that regulates cellular proliferation, DNA repair, and apoptosis as a tumor suppressor. 1,2,3 The gene encoding this transcription factor, TP53, is the most commonly mutated gene across all cancer types. 4 Up to 50% of esophageal adenocarcinomas (ADCs), lung squamous cell carcinomas (SqCCs), pancreatic ADCs, and colorectal ADCs harbor TP53 mutations. [5][6][7][8][9] TP53 mutations most frequently occur in the DNA-binding domain located in amino acid residues 102 and 292. 3,10,11 Of the several mutational hotspots in TP53, eight hotspot mutations constitute 28% of all TP53 mutations despite missense TP53 mutations occurring in nearly 190 codons. 10,12 In colorectal cancer, approximately 37% of all TP53 mutations were comprised with missense mutations in p53 residues R175, R273, R248, and R282 which were four most commonly mutated. 1 In this study, we sought new loci which could be influential to survival period. In addition, R175H is known as the most frequent and significant hotspot in TP53 mutants. The recent studies have suggested that TP53 can undergo gainof-function (GoF) mutations despite its role as a tumor suppressor gene. 11,13 Previous research has suggested that ectopic expression of TP53 with R175H and R273H mutations enhances tumorigenic potential and increases the expression of drug resistance genes. 11,13,14 R175H mutation was considered structural mutation where the mutation causes misfolding of the protein and leads to a significant conformational alterations within p53's DNA binding domain. 1 In this study, we tried to explain the reason why R175H worse the response to immunotherapy by analyzing tumor immune microenvironment through the expression of immune cells and PD-1.
Furthermore, we aimed to analyze the pan-cancer prevalence and profile of mutations in the TP53 gene. We also investigated the effect of TP53 mutations at each amino acid locus and structural domain. We sought to determine the incidence of TP53 mutations and the association between each TP53 locus and survival outcome at a pan-cancer level using a targeted deep sequencing (TDS) panel. Next, we analyzed the impact of TP53 mutation on treatment response to chemotherapy and immunotherapy.
Discussion: Each TP53 mutations indicated differential treatment outcomes following chemotherapy or immunotherapy in patients with metastatic cancer.
Functional analysis including RNASeq suggested that TP53 mutation downregulated immune response.
Conclusion: Overall, we found each TP53 mutation to indicate different prognoses in patients with metastatic tumors undergoing chemotherapy and ICI treatment. Further validations, including a prospective cohort study or a functional study, would be particularly valuable in advancing the knowledge on this aspect and developing improved prognostic parameters.

| DNA extraction
Metastatic tumor regions were micro-dissected for most tumor tissues, except for the samples used in genomic DNA extraction. Genomic DNA was isolated from formalinfixed paraffin-embedded (FFPE) tissue fragments and purified using AllPrep DNA/RNA FFPE Kit (Qiagen). DNA concentration was measured using a Qubit dsDNA HS assay kit (Thermo Fisher Scientific), and 40 ng DNA was used as input for library preparation. DNA integrity number, which is a measure of DNA fragment size and consequently DNA quality, was determined using the Genomic DNA ScreenTape assay on an Agilent 2200 TapeStation system (Agilent Technologies).

| Library preparation and data analysis
Following the manufacturer's protocol, we prepared DNA library applying a hybrid capture-based TruSight Oncology 500 DNA/RNA NextSeq Kit. Using unique molecular identifiers (UMIs) in TruSight Oncology 500 (TSO 500), we determined the unique coverage of each position and lowered background noise. In the process of DNA library preparation, we can detect low variant allele frequencies (VAFs) with reducing errors, thereby guaranteeing high specificity.
We analyzed sequencing data for genomic alterations, including SNVs, CNVs, and fusions. SNVs and small indels which of variant allele frequency (VAF) were less than 2% were excluded. Average copy number more than 4 were considered as amplification. Only amplification was analyzed in the TSO 500-CNV analysis. RNA translocationsupporting reads of which were more than 4-12 were included in translocation. Annotation of data output from the TSO 500 pipeline (Illumina) were performed using the Ensembl Variant Effect Predictor (VEP) Annotation Engine. The variations were marked according to a 4-tier system suggested by the American Society of Clinical Oncology/College of American Pathologists. TMB and microsatellite instability (MSI) statuses was determined by TSO 500 pipeline (Illumina).

| PROVEAN and SIFT score prediction
To know whether an amino acid substitution or indel cause significant impacts on the biological function of a protein, we used a software tool named PROVEAN. 15 After files containing information about amino acid substitutions were uploaded, the files showing the prediction results were saved. This analysis also provided SIFT scores with PROVEAN scores. 16

| Sequencing of whole transcriptome
Total RNA concentration and quality were calculated using Quant-IT RiboGreen (Invitrogen). After samples were run on the TapeStation RNA ScreenTape (Agilent), we used 100 ng of total RNA for sequencing library construction with the TruSeq RNA Access Library Prep Kit (Illumina) following the manufacturer's protocol. After total RNA was fragmented into small pieces, they were copied into first-strand cDNA using SuperScript II reverse transcriptase (Invitrogen) and random primers, following by second-strand cDNA synthesis. The purified products were enriched using PCR to make the cDNA library. After normalization and pooling six into a single hybridization/capture reaction, pooled libraries were incubated with biotinylated oligos which were correspond to the coding regions of the genome. Using streptavidin-conjugated beads, targeted library molecules were collected by hybridized biotinylated oligo probes. The collected libraries were quantified using KAPA Library Quantification kits for Illumina Sequencing platforms following the qPCR Quantification Protocol Guide (KAPA BIOSYSTEMS, #KK4854). Indexed libraries were then subjected to an Illumina HiSeq2500 (Illumina, Inc.), and Macrogen Inc. (South Korea) performed paired-end sequencing.

| Preprocessing of wholetranscriptome sequencing data
Annotation of RNA sequence were performed with ENSEMBL (version 98) and align to the human reference genome (GRCh38) was conducted using STAR (version 2.6.1; ref. 17). We quantified in units of transcript per million (TPM) by applying RSEM (version 1.3.1; ref. 18) pipeline. TPM values less than 1 were considered as zero.

Whole-Transcriptome Sequencing
To explore the activity of signaling pathways, we used gene set variation analysis (GSVA 18 ) and GSEA program in default mode. Enrichment analysis was performed using DAVID (Database for Annotation, Visualization, and Integrated Discovery) (https://david.ncifc rf.gov/). The calculation of composition of immune cells was performed by xCell using immunedeconv R package.

| Statistical analysis
Descriptive statistics are reported as proportions and medians. Data are presented as number (%) for categorical variables. Correlations between clinical characteristics and TP53 status were analyzed by a t-test, Fisher's exact test, chi-square test, or two-way analysis of variance, as appropriate. OS was defined as the time from the first treatment to the date of death. PFS was defined as the time from treatment initiation to the date of disease progression or all-cause mortality. Kaplan-Meier estimates were used in the analysis of all time-to-event variables, and the 95% confidence interval (CI) was computed for the median time to an event. All statistical analyses were performed using R for Windows (version 4.1.2, https://cran.r-proje ct.org/bin/windo ws/ base/). RStudio desktop was used to draw all graphics (RStudio Team, Boston, MA, USA; https://www.rstud io.com/produ cts/rstud io/downl oad/). The following R packages were used for graphics: ggplot2 (3.3.5), Ggally
Before detailed analysis, we also categorized TP53 mutations into damaging and tolerated mutations according to SIFT predictions, 15,16 and further analysis surrounding the relationship between TP53 mutation and TMB, HRD, and PD-L1 (combined positive score: CPS) status was performed. The median values of TMB scores were 6.3, 6.3, and 4.7 in damaging mutation, tolerated mutation, and wild-type TP53 groups, respectively ( Figure 1D). The median number of mutations in HRD genes was 1.0 in all groups ( Figure 1E), whereas the median CPS scores were 1.3, 1.4, and 1.4, respectively, in each group ( Figure 1F).

| Prognostic value of each TP53 mutation locus
Details of mutation loci among cancer types are described in Figure 2A. Most TP53 mutations were missense mutations, whereas nonsense mutations were frequently observed in R196 (n = 24/27, 88.9%), R213 (n = 24/27, 88.9%), R306 (n = 18/18, 100%), and R342 (n = 15/18, 83.3%) loci (Figure 2A). To identifying which loci mutation was influential to survival period, we performed survival analysis with all loci presented in Figure 2A. As a result, the patients having mutation in R213, R342, C242, and R282 had shorter survival period than those having mutation in other loci. ( Figure S2). In addition, we assessed the prognostic value of the TP53 R175H mutation, which is a well-established marker of poor prognosis. We found that the OS and PFS of patients with R175H were shorter than those of other SNVs following both chemotherapy and ICI treatment ( Figure S2). Therefore, we evaluated the biological functional damage of mutation in each loci by calculating PROVEAN score. The PROVEAN score of each locus showed that R196 and R231 loci had lower scores (median PROVEAN score: −14.0 and 14.1, respectively) than the other loci (median: −5.8) (p < 0.001) ( Figure 2B). In addition, C242 and R342 loci had relatively lower scores (median: −10.7 and −9.5, respectively) than other loci (p < 0.001) ( Figure 2B). We further analyzed the relationship between PD-L1 status and TP53 mutation loci in 187 cancers treated with immune checkpoint inhibitors (ICIs) ( Figure 2C). In this analysis, PD-L1 positivity defined as CPS ≥1 was associated with TP53 R342 (41.2%) and R213 (40.7%) mutations, whereas PD-L1 positivity was inversely associated with C242 (15.0%) and R282 (13.2%) mutations (p < 0.001). There was not clear correlation between PD-L1 and loci, although each loci has various PD-L1 positivity depending on mutation site. These results suggested the loci of mutation could affect PD-L1 status. In the patients harboring TP53 mutations, the most frequently  co-occurring mutation was found in APC, followed by KRAS, HIST1H1C, LRP1B, and SPTA1 ( Figure 2D).

| The difference of gene expression profile as a factor contributing to survival in R175H mutation
It is known that there are some hotspots in TP53 mutations. Among them, although R175H mutations is most frequent mutation, both the correlation between R175H mutation and the result of immunotherapy and the changes of tumor immune microenvironment by R175H mutation have not yet been sufficiently reported. In an effort to elucidate potential transcriptional signature or molecular determinants that are associated with response to ICI therapy, we compared gene expression profiles for tumors having R175H mutation in TP53 with TP53 wildtype tumors. In total 64 CRC tumors, six had R175H in TP53 gene. The genes (n = 100) with the greatest difference in their expression between two groups were presented as the heat map ( Figure 3A).In fold enrichment analysis by GSEA, chemokine signaling pathway was reduced in tumors having R175H mutation compared to wild-type tumors (normalized enrichment score: −1.83, false discovery rate[FDR] q-value: 0.004) ( Figure 3B upper panel). In addition, antigen-presenting pathway (normalized enrichment score: −1.83, FDR q-value: 0.004) and immune network for IgA production pathway (normalized enrichment score:-1.87, FDR q-value: 0.003) were downregulated in R175H mutation patients ( Figure 3B middle and lower panel).
With 50 genes whose expression was the most lowered in R175H mutated tumors, we performed the analysis of fold enrichment pathways with DAVID bioinformatics tools ( Figure 3C). As a result, following pathways were downregulated in R175H-mutated tumors: adaptive immune response (p < 0.001), antigen processing and presentation of peptide antigen via MHC class I (p = 0.024), positive regulation of B-cell activation (p = 0.024), chemokine-mediated signaling pathway (p = 0.025), innate immune response (p < 0.001) ( Figure 3C). We also analyzed immune cell abundance base on expression of gene related to immunity using xCell. Especially, Plasma B cells (p = 0.032), granulocyte monocyte progenitor cells (p = 012), macrophages (p < 0.001), M1 macrophages (p = 0.017), naïve CD4(+) T cells (p = 0.042), CD8(+) T cells (p = 0.008), NK T cells (p = 0.005) were significantly lowered in R175H patients ( Figure 3D). According to previous article, PD-1 expression balance between CD8 + T cells and Treg cells predicts the clinical efficacy of PD-1 blockade therapies. Therefore, we checked PD-1 and PD-L1 expression by RNA-seq, showing the expression of two genes were significantly lowered in R175H mutant samples (p = 0.002 and p < 0.001) ( Figure 3E). The reduction in immune cells and the expression of PD-1 could explain some part of reason why R175H mutation patients had shorter survival period to immunotherapy.

TP53 mutation
Copy number variants (CNVs) associated with TP53 mutations were evaluated. Across copy number variations, only amplifications (copy number ≥4) were reported according to the our hospital reporting system (deletions were not reported). In patients with TP53 mutation, 69 genes were amplified in total. Among them, MYC (n = 132), RICTOR (n = 73), EGFR (n = 68), CCNE1 (n = 64), LAMP1 (n = 61), and BRCA2 (n = 54) were the most frequently amplified genes in cancers harboring TP53 mutations ( Figure 5A). Among these CNV events, CCND1, FGF4, and FGF19 in chromosome 11 were strongly related ( Figure 5B). According to survival analysis, tumors with TP53 mutation and co-existing copy number amplification of the three aforementioned genes conferred poor survival outcomes compared to those with only TP53 mutation following chemotherapy (p < 0.05) ( Figure 5C-E). However, there was no such significant difference pertinent to ICI treatment (data not shown).

| Features of good responses to chemotherapy according to machine learning
Finally, we tried to predict the response to chemotherapy using the neuralnet R package. A neural network is a computational system that creates predictions based on existing data. A neural network consists of input layer, hidden layer, and output layer. The model was designed to have five hidden layers with a threshold of 0.01 ( Figure S4). We tested the number of hidden layer from 1 to 10, and we chose 5 hidden layers because prediction rate was highest. We divided our data to 80% for training set and 20% for test set. We input (1) the CNV data of CCND1, FGF4, and FGF19; (2) TP53 hotspot mutations (C242, R282, R342, and R175); (3) deleterious mutations predicted by PROVEAN; and (4) TP53 NES or E4F1binding domain mutations and predicted the treatment response (CR/PR or SD/PD). Among the four factors, the TP53 mutation domain was the most weighted (0.35), followed by PROVEAN score (0.28), TP53 loci (0.24), and CNV (0.14) ( Figure S4). The accuracy of this model was 0.84 (95% CI: 0.8035, 0.8637); its sensitivity was 1.00, and its specificity was 0.59 ( Figure S4), showing this model's sensitivity is very high. According to this model, the most important factor determining the response to chemotherapy was the domain was involved in TP53 mutation ( Figure S4). Next, the response is determined by whether the mutation is predicted as deleterious or not as determined by the PROVEAN score. The third factor was the presence of the mutation locus in a hotspot. The final factor was co-occurrence with the amplification of FGF19, FGF4, or CCND1 ( Figure S4).

| DISCUSSION
TP53 is a well-known tumor suppressor gene, but its mechanisms of action are complex. TP53 is involved in genomic instability, angiogenesis, replicative cell immortality, tumor invasion, and tumor-promoting inflammation. [19][20][21][22][23] In this study, we investigated the prevalence and profile of TP53 mutations in patients with different types of cancers who underwent clinical NGS analysis. In addition, we analyzed genomic mutational factors affecting the survival of these patients. All genomic mutations were investigated after dividing them according to mutational locus (amino acid position), structural domains, and degree of functional protein damage. Lastly, we developed a machine learning model that could predict the response to chemotherapy/ immunotherapy and calculate the importance of each factor. According to our model, domain harboring mutation was the most important factor which could determine the response of chemotherapy/immunotherapy, following by PROVEAN score, loci where mutation occurred and co-occurring CNVs.
First, in terms of domain, structural domain analysis showed that mutations occurring in the NES or E4F1binding domain were critical for survival after both chemotherapy and immunotherapy. The NES domain exports  protein and the attached RNA into the cytoplasm through the nuclear pore. 24 Many studies have demonstrated the importance of the NES domain in proteins related to cancer. [25][26][27] For example, EGFR NES mutants with increased nuclear accumulation promote the aggressiveness of cancer cells, including TKI resistance. 28 Moreover, singleallele mutations in the EGFR NES domain are capable of driving lymphomagenesis, suggesting that the EGFR NES mutant possesses strong oncogenic activity. 28 Similar to this mechanism, TP53 proteins bearing mutations in the NES domain cannot export proteins to the cytoplasm, which accumulate in the nucleus and block the function of TP53. E4F1 is a key posttranslational regulator of TP53, which modulates its effector functions involved in alternative cell fates: growth arrest or apoptosis. E4F1 mediates mutually exclusive posttranslational modifications of TP53. 29 E4F1-dependent Ub-TP53 conjugates are associated with chromatin, and their stimulation coincides with the induction of a TP53-dependent transcriptional program specifically involved in cell cycle arrest and not apoptosis. 29 According to this study, mutations in the E4F1-binding domain of TP53 interfere with posttranslational modifications of TP53, suggesting loss of TP53 function.
In terms of loci, Specifically, missense mutations in the four most commonly mutated p53 residues (R175, R248, R273, and R282) comprise most of all TP53 mutations. Among these hotspot mutations, p53-R175H has the highest occurrence. Although losing the transactivating function of the wild-type p53 and prone to aggregation, p53-R175H gains oncogenic functions by interacting with many proteins. We found another important loci (R213, R342, C242, and R282) which can cause short survival period. Our study also suggested that the R175H mutation was associated with shorter PFS and OS than other mutations and wild-type TP53 following both chemotherapy and ICI treatment. Recent studies have suggested that TP53 can undergo GoF mutations despite its role as a tumor suppressor gene. 11,13 Ectopic expression of TP53 with the R175H mutation enhances tumorigenic potential and increases the expression of drug resistance genes. 11,13,14 The R175H mutation regulates doxorubicininduced apoptosis. 30,31 TP53 R175H mutations in lung cancer result in drug resistance against etoposide and cisplatin. 32 Functionally, tumors harboring TP53 GoF mutations have been found to progress more rapidly and confer poorer survival outcomes than tumors with deleterious TP53 mutations. 33,34 In this study, we tried to look for the reason why R175H mutation made the response to chemotherapy/immunotherapy worse by analyzing tumor microenvironment using RNA-seq. In RNA sequencing data, we found that the expression of TP53 was elevated in TP53 R175H compared to wild-type and GoF of TP53 could exert on tumor environment. By analysis GO term related differentially expressed gene between normal and R175H mutation patients, both innate and adaptive immune response pathways including cytokine-mediated signaling pathway were lowered in R175H-mutant tumors. Furthermore, some immune cells proportions were significantly lowered in R175H patients. In addition, PD-1 and PD-L1 expression were also significantly lower in R175H mutant tumors compared to others. According to previous article, PD-1 expression balance between CD8 + T cells and Treg cells predicts the clinical efficacy of PD-1 blockade therapies. The reduction of immune cells and the expression of PD-1 could explain some part of reason why the response to immunotherapy in patients harboring TP53 R175H mutation was worse than patients with wildtype TP53. Moreover, a recent study targeting the R175H locus of TP53 revealed a cancer treatment strategy including therapeutic antibodies targeting neoantigens derived from a common TP53 mutation. 35 We suggest that other loci of TP53 mutations might also be effective targets for T-cell therapy in the near future.
We also focused on copy number amplification in tumors harboring TP53 mutations. Interestingly, FGF4, FGF19, and CCND1 copy number amplifications were associated with each other. These three copy number amplifications with TP53 mutations conferred poor survival outcomes following chemotherapy. Interestingly, chromosome 11 allelotypes and TP53 have been reported to collaborate to induce chemical carcinogenesis. 36,37 Further functional research will confirm the interaction between TP53 and gene amplification in chromosome 11.
Our study had a limitation. In our database, TP53 R175H mutation had worse prognosis in patients treated with immune checkpoint inhibitor. However, we could not validate our result. We had tried to validate our finding for a predictive role of R175H mutation but public database including the cancer genome atlas did not have treatment information of immune checkpoint inhibitor. Moreover, we had tried to search the references for specific TP53 genetic alterations and immunotherapy, but there were few articles to analyze the relationship between TP53 genetic alteration and immunotherapy. Further functional study would be warranted to support our data. Overall, we found each TP53 mutation to indicate different prognoses in patients with metastatic tumors undergoing chemotherapy and ICI treatment. This information would be of great value for the determination of patient prognosis in the metastatic setting. The model showed high sensitivity and specificity. We expect to improve this model by collecting and applying more data in the future. Further validations, including a prospective cohort study or a functional study, would be particularly valuable in advancing the knowledge on this aspect and developing improved prognostic parameters.

AUTHOR CONTRIBUTIONS
Ji-Yeon Kim: Conceptualization (lead); data curation (equal); formal analysis (lead); funding acquisition (equal); writing -original draft (lead); writing -review and editing (equal). Jaeyun Jung: Data curation (lead); formal analysis (lead); methodology (lead); resources (equal); software (lead); writing -original draft (equal); writing -review and editing (equal). Kyoung-Mee Kim: Methodology (equal); resources (equal); supervision (equal). Jeeyun Lee: Project administration (lead); supervision (lead); writing -original draft (equal); writing -review and editing (lead). Young-Hyuck Im: Investigation (equal); project administration (equal); supervision (lead); writing -original draft (equal); writing -review and editing (equal). F I G U R E 5 Copy number variants in patients with TP53 mutations and their effect on overall survival. (A) Circled bar graph representing the genes in which amplification occurred in patients with TP53 mutations. Bar height indicates the number of CNVs in each gene. The black point in the middle of the bar shows the mean value of copy number in each gene. The genes are ordered and colored by chromosome number. (B) Correlation between the occurrence of TP53 mutation and the occurrence of each CNV. CNVs represented in the graph indicate the most frequent CNVs among all CNVs (frequency ≥20). Kaplan-Meier-estimated OS curve in patients treated with chemotherapy who have copy number amplifications in CCND1 (C), FGF19 (D), and FGF4 (E) genes (red line: group with CNVs in specific genes and SNVs in the TP53 gene; blue line: group with SNVs only in the TP53 gene).