Whole-exome Sequencing in Penile Squamous Cell Carcinoma Uncovers Novel Prognostic Categorization and Drug Targets Similar to Head and Neck Squamous Cell Carcinoma Free

Penile squamous cell carcinoma (PSCC) is a rare cancer in the United States with approximately 2,000 new cases and 400 deaths yearly (1). The incidence of PSCC is significantly higher in developing countries worldwide (2). Human papillomavirus (HPV) is the known etiology for approximately 30%–50% of PSCC cases globally (3–5). The current treatment approach for locally advanced PSCC, using neoadjuvant platinum-based chemotherapy followed by surgical consolidation offers long-term disease control for less than 30% of patients with locally advanced PSCC (6); unfortunately, the majority of the patients will have disease progression or relapse. Relapsed PSCC after first-line chemotherapy is associated with poor survival, with an estimated median overall survival (OS) of 5 months (7) and very limited treatment options (8). Thus, it is important to perform molecular investigations to identify recurrent alterations associated with PSCC and potential molecular targets for drug development in this rare disease with fatal outcomes (9, 10). The current evidence investigating the genomic landscape of PSCC has been limited by small sample size, limited HPV-positive (HPV⁺) cases, and heterogeneous molecular testing methods (11–16). Prior efforts have identified two different pathways for penile oncogenesis, one that is HPV related and a second that is not. In addition, only a few recurrent point mutations in known tumor suppressor genes have been reported in PSCCC, including TP53, CDKN2A, and PIK3CA (11–16). Therefore, the extent to which other genes are contributing to the disease, and thus, constitute additional targets for potential therapeutic exploration is not well characterized. Therefore, a pressing need exists for the identification of new therapies, based on an improved insight of the disease biology. Other studies have found strong evidence for molecular convergence between squamous cell carcinomas (SCC) from different anatomic sites (17), but this has not been explored yet in PSCC. In this study, we performed the largest whole-exome analysis to date of 34 PSCCs to explore the genetics of this rare cancer.

PSCC and matched normal penile tissues (n = 34 pairs) were obtained from the University of Texas MD Anderson Cancer Center (UT MDACC, Houston, TX) tissue bank of prospectively collected samples. The diagnosis of PSCC was made based on the review of an expert pathologist specializing in genitourinary cancers (P. Rao), who assessed and reported the tumor histologic features per the American Joint Committee on Cancer 8th edition (18). All specimens were obtained from patients who were evaluated at UT MDACC (Houston, TX), and following informed consent under protocols approved by the institutional review board. Paired PSCC and normal tissues were obtained snap frozen from the same patient at the time of surgery. Most tumors were collected prior to receiving chemotherapy (9/10), when indicated clinically. The percentage of malignant tissue was calculated by histologic examination, following hematoxylin and eosin (H&E) staining. All PSCCs analyzed comprised at least 50% malignant penile cells. The studies were conducted in accordance with the Declaration of Helsinki and informed written consent was obtained from each subject or each subject's guardian.

For IHC for p16^INK4a, all samples were tested, and staining was performed on a BenchMark Autostainer (Ventana Medical Systems) as described by the manufacturer's protocol using a prediluted mouse mAb (CINtec p16 Histology, clone E6H4, Ventana Medical Systems). Microscopic slide evaluation for p16^INK4a staining patterns was classified as follows: pattern 0, complete absence of p16INK4a staining in all the neoplastic cells; pattern 1, spotty and discontinuous individual staining in some of the neoplastic cells; pattern 2, a more extensive, albeit discontinuous staining pattern with small clusters of positive neoplastic cells; and pattern 3, entire and continuous cytoplasmic or nuclear staining in all neoplastic cells, except in hyperkeratotic or parakeratotic areas when present. On the basis of prior analyses, only pattern 3 was considered as positive for p16^INK4a overexpression (19, 20). The p16^INK4a staining patterns were confirmed by an experienced pathologist (P. Rao).

HPV genotyping for all samples was performed using Cobas HPV Assay (Roche Diagnostics) as described elsewhere (21). First, the formalin-fixed, paraffin-embedded (FFPE) tissue sections were prepared from each block (5 μm) with the last section stained with H&E and evaluated to determine the quality of the specimen. Two to three unstained tissue sections were used for specimen preparation. Briefly, the FFPE tissue sections were deparaffinized by xylene (×3) for 10 minutes each, followed by an ethanol wash (×3) to remove paraffin. Deparaffinized tissue pellets placed in 1.5 mL microcentrifuge tube were air dried and resuspended in 200 μL of Lysis Reagent (10% proteinase K in an ATL Buffer, Qiagen). The mixture was incubated in a 56°C water bath for 1 hour and diluted into 2 mL with 50% ethanol. After deparaffinization and specimen pretreatment, 1 mL solution was applied to Cobas 4800 system, following the manufacturer's instruction. Cobas HPV assay is an FDA-approved HPV testing assay for cervical cancer screening in Pap cytology specimens. The Cobas 4800 system is a highly automated HPV testing platform, including DNA extraction, a real-time PCR thermocycler, and an analyzer for the testing result reporting. In each testing run, positive and negative controls, as well as an internal control with β-globin amplification to determine the presence of DNA, were included. The testing results were reported as HPV 16 and 18 genotypes, as well as other collectively 12 high-risk HPV types (HPV 31, 33, 35, 39, 45, 51, 52, 56, 58, 59, 66, and 68). The FFPE specimen in Cobas 4800 system for HPV testing was validated in-house in our institution.

Genomic DNA was extracted from tumor and peripheral nonmalignant frozen tissue using the QIAmp DNA Tissue Kit (Qiagen) following the manufacturer's instructions. Genomic DNA was quantified by PicoGreen (Invitrogen) and quality was accessed using Genomic DNA Tape for the 2200 TapeStation (Agilent Technologies). Only samples with more than 200 ng of genomic DNA (optimal quantity threshold) and with DNA integrity number higher than 7 (optimal quality threshold) were considered for library preparation. Exome capture was performed with the SeqCap EZ Exome Probes v3.0 (Roche), which covers 64 megabases (Mb) of human genome (GRCh37). The captured libraries were sequenced on a HiSeq 4000 (Illumina Inc.) for 2 × 100 paired end reads, using Cycle Sequencing v3 Reagents (Illumina). Only sequenced samples with >80% of data with quality score ≥ Q30 (99.9% of base call accuracy) were used in the analysis.

The resulting raw output BCL files containing the sequence data were converted into FASTQ files using CASAVA 1.8.2. FASTQ files were aligned to the reference genome (human Hg19) using BWA (22). The aligned BAM files were subjected to mark duplication, realignment, and recalibration using Picard and GATK (23) before any downstream analyses. The mean number of total reads generated for the 34 tumor samples was 185.1 ± 26.6 million, with mean of 99.3% ± 0.1% of DNA mapping rate and mean coverage of 141.1× ± 20.5×. For the 34 paired normal tissue samples, the mean number of reads generated was 99.8 ± 14.8 million, with 99% ± 0.1% of mean DNA mapping rate and mean coverage of 70.2× ± 9.2×. Somatic mutation calls were performed using MuTect (24), and insertions and deletion (indel) calls using Pindel (24). Resulting data were combined into mutation annotation format (MAF) files, and then results were manually filtered on the basis of the following parameters: nucleotide changes were considered only when we observed a minimum of 20× coverage for tumor samples and 10× for normal samples, and with variant allele fraction (VAF) higher or equal to 0.02 for tumors and lower than or equal to 0.02 for normal samples for single-nucleotide variants (SNV). For INDELs, the minimum VAF considered was 0.05. All variants described as SNPs at 1000 Genome Project database (25), ESP6500 (26), and EXAC (27), with frequency reported as higher than 0.01, were excluded. All events occurring in regions with local repeats were also excluded.

PSCC filtered variants organized in an MAF file were analyzed using the R package, Maftools (28), which congregates multiple genomic analytic tools, as briefly described below. Estimation of hypermutated genomic regions was performed by assessing interevent distances in each chromosome. Genomic areas exhibiting six or more consecutive events with an average intermutation distance of 1,000 bp were considered hypermutated. PSCC driver genes were estimated by mutation positional analysis, which searches for mutation enrichment in specific loci, or hotspots, as described by Tamborero and colleagues (29). We also assessed the enrichment level of 10 oncogenic signaling pathways (RTK/RAS pathway, Wnt pathway, Notch pathway, Hippo pathway, PI3K pathway, cell-cycle pathway, Myc pathway, TGFβ pathway, p53 pathway, and Nrf2 pathway) in our PSCC samples on the basis of the curated set of genes described by Sanchez-Vega and colleagues (30). In our work, we considered all missense and truncated mutations occurring in genes described in their gene set.

Estimation of gene druggability was performed by accessing data from the DGIdb database that include gene–drug interaction information (31). Tumor variants were classified into one of the 96 nucleotide substitution classes that consider the mutant base and the immediate 3′ and 5′ nucleotides around it. The classified variants were then used for tumor mutation signature estimation by two methods. The first method aimed to identify APOBEC-related mutation enrichments. APOBEC-induced mutations tend to be characterized by an enrichment of TpCpW motif, in which W corresponds to an A or T. Using the approach proposed by Roberts and colleagues (32), the enrichment of C>T mutations associated with the TpCpW motif was compared with all C>T events and other background events to calculate an enrichment score. This enrichment score was then evaluated by Fisher exact test, and then, samples were categorized into APOBEC-enriched and non-APOBEC–enriched groups. The second method of signature estimation focused on the identification of PSCC mutational patterns and association of these patterns with the version 3 of the mutation signatures [part of the v89 Catalogue of Somatic Mutations in Cancer (COSMIC), released in May 2019] from COSMIC database (cancer.sanger.ac.uk; ref. 33). Initial step involves the estimation of the number of different mutation patterns among the PSCC samples by nonnegative matrix factorization analysis and cophenetic correlation. Then, the nucleotide substitution classes are decomposed into the estimated number of PSCC signatures and finally compared with COSMIC mutation signatures and the best match is identified by cosine similarity calculation.

Aiming to compare PSCC genomic data with other known squamous carcinomas, we retrieved MAF files from The Cancer Genome Atlas (TCGA) studies with carcinoma cohorts (34) using TCGA mutations R package, which included the following sites: bladder (dx.doi.org/10.7908/C1MW2GGF), cervix (dx.doi.org/10.7908/C1MG7NV6), esophagus (dx.doi.org/10.7908/C1BV7FZC), head and neck (dx.doi.org/10.7908/C18C9VM5), and lung (dx.doi.org/10.7908/C1×34WXV). Because these studies included several histologic variants of carcinomas, we further filtered these MAF files by including only those with squamous histology, as described by Campbell and colleagues (35). The list of included samples is described in the Supplementary Table S1.

MAF files of each cohort were analyzed by using the Maftools R package, in which we searched for differentially mutant genes between PSCC and squamous TCGA cohorts. Finally, we estimated APOBEC and COSMIC mutation signatures, as described previously, for all TCGA squamous tumors and PSCC samples together, by compiling mutation information in a single MAF file.

Associations between categorical variables were assessed by Fisher exact test, or χ² test, whenever appropriate. Associations between categorical and continuous variables were assessed by Mann–Whitney U test for pairwise comparison and by Kruskal–Wallis test followed by Tukey test for three or more categorical features. Correlations between continuous variables were assessed by Pearson correlation test. For survival analysis, we performed the log-rank test and univariate Cox regression. Analyses were performed using JMP 15.0 Software.

Collection and use of patient samples for this study were approved by the UT MDACC (Houston, TX) Institutional Review Board (PA15-0138).

Thirty-four untreated primary PSCC specimens and corresponding matching normal penile tissue samples were subjected to whole-exome sequencing (WES). Clinicopathologic characteristics of these patients are described in Table 1. The median age at diagnosis for this cohort was 61 years (range, 39–83 years). Twenty-one patients were Caucasians (62%), 11 were Hispanic (32%), and 2 (6%) were African-Americans. The cohort included 10 (29%) tumors that were HPV⁺ as assessed by Cobas HPV test, but only eight of those were p16 positive (p16⁺) by IHC staining. Most patients (n = 24, 70%) had usual SCC histology, while six (17%) tumors were basaloid, two (6%) were papillary not otherwise specified (NOS), and two (6%) were classified as verrucous carcinoma. Pathologic primary cancer staging associated with high metastatic potential risk was observed in (T2–T4) 24 patients (71%) with T2–T4 tumors, in 25 patients (74%) with grade 2–3 tumors, in 16 (47%) patients with lymphovascular invasion, and in 12 (46%) patients with perineural invasion. On the other hand, most patients (n = 18, 53%) had clinically N0 disease at the time of diagnosis, defined as no evidence of palpable lymph nodes based on physical exam and nodal imaging at diagnosis. Pathologically, 21 patients (68%) with enough follow-up were pN0. In this cohort, 10 (29.4%) patients received neoadjuvant chemotherapy, the majority (n = 9, 90%) received the combination chemotherapy regimen with cisplatin, ifosfamide, and paclitaxel, followed by consolidation inguinal and pelvic lymph node dissection in all 10 patients. The median follow-up time for the entire cohort was 32.2 months, during this period, 5 (14.7%) patients died [only 3 (8.8%) died from PSCC], 3 (8.8%) patients were alive with disease, and 26 (76.4%) remained alive without evidence of disease. Only 8 (24%) patients had disease progression or recurrence. A detailed description of the clinical and pathologic characteristics is presented in Supplementary Table S1.

WES was performed on 34 tumor–normal pairs. All samples were sequenced using a whole-exome capture system covering 64 Mb of the human genome (CRCh37). Sequencing depth of targeted regions was 141.1× ± 20.2× for tumor samples and 70.2× ± 9.1× for normal samples, with a mapping rate of 99.3% ± 0.1% and 99.1% ± 0.1%, respectively. Collectively, the samples presented 4,275 somatic exonic mutations, including 3,717 missense, 327 nonsense, 1,607 silent, 75 splice site, 95 frameshift indels, and 61 in-frame indels. The full set of nucleotide variants found in PSCC samples is described in the Supplementary Table S2. The 30 most frequent exonic mutations identified in PSCC are described in the oncoplot in Fig. 1A.

Patient's tumor mutation burden (TMB) was determined by dividing the total number of SNVs by the amount of the genome covered by whole exome sequencing (∼64 Mb). TMB was then represented by the number of SNVs per Mb. The mean TMB among patients with PSCC was 1.96 (ranging from 0.01 to 21.2). On the basis of Tukey outlier test, 2 patients (14T and 27T) were classified as hypermutators (Fig. 1B). Analysis of intermutation distances further confirmed this classification by identifying areas of hypermutated loci (six consecutive mutations with an average distance of less than or equal to 1,000 bp) only in these two samples (Fig. 1C).

On the basis of the manually curated gene set presented by Sanchez-Vega and colleagues (30), we observed that multiple genes associated with oncogenic pathways were highly mutated in PSCC samples (Supplementary Fig. S1). Notch pathway was the most frequently affected in PSCC samples, in which 70.6% of samples (n = 24) had mutations in at least one of the 17 Notch genes, including NOTCH1 (n = 15, 44.1%), EP300 (n = 5, 14.7%), and FBXW7 (n = 5, 14.7%; Fig. 1D). Other highly affected pathways were: Hippo pathway with 17 mutant genes among 20 (58.8%) patients, in which 12 (35.3%) exhibited FAT1 mutations; RTK-RAS with 22 mutant genes among 16 patients (47.1%), which included three (8.8%) cases with HRAS mutations and other three with NF1 mutations; p53 pathway with only three mutant genes among 15 patients (44.1%), but with TP53 mutations involving 12 (35.3%) of them; and PI3K, cell cycle, and Wnt pathways were mutated in 11 (32.3%) patients, and PIK3CA, CDKN2A (n = 10, 29.4%), and CHD8 (n = 3, 8.8%) mutations were the most frequent for each pathway, respectively. Nrf2 (n = 4, 11.7%), TGFβ (n = 3, 8.8%), and MYC (n = 3, 8.8%) were less enriched among PSCC samples, and no gene was mutated in more than two samples. TMB was significantly higher among patients with mutations in Wnt [mutation (Mut) = 4.15 ± 5.8 vs. wild-type (WT) = 0.9 ± 0.7; P = 0.0003], RTK-RAS (Mut = 3.16 ± 5.1 vs. WT = 0.90 ± 0.7; P = 0.010), Notch (Mut = 2.44 ± 0.8 vs. WT = 0.80 ± 0.3; P = 0.014), and p53 (Mut = 2.01 ± 1.6 vs. WT = 1.92 ± 4.7; P = 0.013) pathways.

Positional enrichment analysis indicated that mutations in PIK3CA and CDKN2A were significantly enriched in specific protein loci (FDR = 0.007 and FDR = 0.039, respectively), suggesting these may be potential driver mutations in PSCC (Supplementary Fig. S2). The PIK3CA mutations, p.E545K/G, E542K, and E453K, are all considered gain-of-function oncogenic mutations by OncoKB, while the CDKN2A mutations, p.R80X, p.E88X, p.G65fs, p.D108G, and p.R115fs, are classified as likely loss-of-function and likely oncogenic (36). The most frequent SNVs among patients with PSCC were characterized by C>T (mean of 48% of all SNVs/patient) changes, followed by C>G (16.5%), C>A (15.4%), T>C (12.1%), T>G (4.2%), and T>A (3.8%; Fig. 1A).

We estimated the mutational patterns for each patient based on the frequency of single-base substitutions in the context of immediate 5′ and 3′ bases, classifying 96 different mutation types. By performing cophenetic analysis, we identified two distinct mutation patterns among PSCC samples, which we called mutation pattern 1 (MP1) and 2 (MP2; Fig. 2A and B). Enrichment of each of these two mutation patterns for each individual sample is described in the Fig. 2C. MP1 was the most prevalent in nine PSCC samples, while MP2 was more prevalent in the other 25 cases. These mutation patterns were then compared with known cancer mutation signatures from COSMIC database (Mutation Signatures V3, May 2019). Cosine similarity scores indicated the closest matches between PSCC mutation patterns and COSMIC signatures (Fig. 2D). The MP1 showed the highest cosine similarity to the SBS2 (0.767) COSMIC signature, which has as a proposed etiology the activity of AID/APOBEC family of cytidine deaminases (C>T – TC context). MP1 also presented high cosine similarity with SBS13 (0.618) and SBS7a (0.644) COSMIC signatures. SBS13 is also associated with AID/APOBEC cytidine deaminases family, and both are related to local hypermutation phenomenon called kataegis (C>G – TC context). Interestingly, MP1 enrichment was positively correlated with increased TMB (CC = 0.71; P < 0.0001). By the other side, SBS7a signature is commonly found in tumors associated with sun exposure and are likely caused by UV radiation (C>T − TC/CC).

The MP2 had the highest cosine similarity with the signature SBS6 (0.836), which is associated with defective DNA mismatch repair system, microsatellite instability, and small indels, as well as the other enriched signature SBS15 (0.708). Another highly enriched signature in MP2 was the SBS1 (cosine similarity = 0.796), which is associated with deamination of 5-methylcytosine to thymine, especially in the NpCpG context.

Considering that the MP1 was highly associated with mutation signatures associated with APOBEC cytidine deaminases activity, we also assessed the enrichment of APOBEC-related mutagenesis pattern by using an algorithm described by Roberts and colleagues (2013), which considers the rate of TpCpW (in which W may be an adenine or a thymine) changes as a marker of APOBEC-related mutagenesis process (32). We estimated that 38% (n = 13) of PSCC samples had an enriched APOBEC-related mutation signature (Fig. 2C) and supports the existence of this subset of PSCC. These APOBEC-enriched samples had a significantly higher TMB (3.87 ± 5.43) compared with non-APOBEC–enriched samples (0.78 ± 0.55; P < 0.0001).

High p16 immunoexpression is a clinically approved surrogate marker of HPV infection for cervical and oropharyngeal tumors. Although HPV detection and p16 staining were associated in the majority of PSCC cases, in two samples (17T and 2T) HPV was detected by HPV PCR testing, but the tumors were p16 negative by IHC. For both cases, nonsense CDKN2A mutations were detected, which could explain the absence of p16 immunoexpression. However, the case 17T showed TP53 and CASP8 mutations, which are rare among HPV⁺ tumors from TCGA database. Whether these two cases were HPV driven is unknown.

HPV detection was significantly associated with mutations in the ARPP21, CMYA5, and RPGRIP1 genes (all mutant patients were HPV⁺; P = 0.004, P = 0.020, and P = 0.020, respectively). P53 pathway genes (TP53 and CHECK2) were mutated in only two (20%) HPV⁺ cases, while more than half (n = 13, 54.2%) of the HPV-negative cases had mutations in one of these genes (P = 0.128). No Nrf2 pathway genes were mutated in HPV⁺ PSCC.

p16 immunoexpression was significantly associated with CSPG4 (3/4 mutant patients were p16⁺; P = 0.032) and TP53 (all p16⁺ tumors were WT; P = 0.030). Although not significant (P = 0.072), no p16⁺ case presented CDKN2A mutations. P53 pathway genes affected 14 p16-negative tumors (53.8%), but only one (12.5%) p16-positive case (P = 0.053). No Nrf2 pathway mutations were detected among p16⁺ tumors.

Local recurrence was observed in 20% (n = 3) of patients with mutations in p53 pathway, while no local recurrence was detected in patients with no p53 pathway mutations (P = 0.076). In fact, all local recurrence cases occurred among TP53-mutant patients (P = 0.037). Although development of regional recurrence was rare (n = 3), 2 (66.7%) of these patients were TGFβ pathway mutants (P = 0.015) and all had Notch pathway mutations (P = 0.538).

Patients with mutations in the Notch pathway showed a trend toward worse OS (log-rank P = 0.150), and all deceased patients (n = 5) were mutant for this pathway (20.9% of all Notch pathway–mutant patients). On the other hand, patients with mutations in the PI3K pathway tended to have improved progression-free survival when compared with patients with no mutations in this pathway (log-rank P = 0.196). Progression occurred in 1 mutant patient (9.1%), while 7 (30.4%) patients without PI3K pathway mutations progressed.

Considering the mutational pattern scores, enrichment of MP1 was significantly associated with worse OS status (Fig. 2E; P = 0.011). In fact, univariate Cox regression showed that increased MP1 score was associated with lower OS [HR, 10.2 (1.13–92.9); Wald test P = 0.039]. Because MP1 score was highly enriched by mutational signatures associated with APOBEC cytidine deaminases activity, we also explored the impact of APOBEC scores in these samples. High APOBEC score was significantly associated with mortality (P = 0.009) and presence of lymph node metastasis (P = 0.024). Increased APOBEC score was significantly associated with poor OS by univariate Cox regression [HR, 2.41 (1.11–6.77); Wald test P = 0.042]. When samples were dichotomized in APOBEC-enriched and non-APOBEC, a trend toward worse survival among APOBEC-enriched samples was observed, although not statistically significant (Fig. 2F; P = 0.060; HR, 6.2; 95% confidence interval, 0.7–56).

Among the mutant genes found in PSCC samples, 11 have been described in the DGIdb database as targetable by known drugs (PIK3CA, 338 drugs; TP53, 193 drugs; CACNA1C, 44 drugs, CDKN2A, 31 drugs; NOTCH1, 26 drugs; FBXW7, 11 drugs; EP300, 10 drugs; MUC2, nine drugs; CASP8, six drugs; LAMA1, one drug; and TTN, one drug). Drug–gene interaction information is detailed in the Supplementary Table S3. From these drug–gene interactions, two of these have been the basis for the FDA drug authorization of a PI3Kα-specific inhibitor and CDK4/6 inhibitor to treat patients with advanced breast cancer (37, 38). In addition, deleterious mutations in DNA damage responsive (DDR) genes are frequently associated with response to PARP inhibitors and platinum chemotherapy (39–41). The DDR pathway is important in tumor biology, allowing cancer cells a mechanism to resist damage by chemotherapy and radiotherapy (42). BRCA1/2 are the most well-described genes in the pathway, but several others (ATM, ATR, PALB2, etc.) are involved with DDR and are mutated in many cancers. However, much remains unknown about their prevalence in PSCC, therefore, we report the prevalence of likely pathogenic variants in DDR genes. We identified nine genes (BRCA1/2, ARID1A, ATR, CHEK2, PARP1, FANCA, PALB2, and RAD51) as DDR-related on the basis of PubMed searches and the NCBI Gene and Biosystems Databases. A total of 26 pathogenic or likely pathogenic variants in DDR genes were identified in 8 of 34 (23.5%) patients. These variants were found in the following genes: ARID1AM 2 patients (5.9%), BRCA2 2 patients (5.9%), ATR 2 patients (5.9%), CHEK2 2 patients (5.9%), PARP1 2 patients (5.9%), FANCA 1 patient (2.9%), PALB2 1 patient (2.9%), and RAD51 1 patient (2.9%).

We compared the genomic alterations of PSCC with head and neck squamous cell carcinoma (HNSC; n = 501), cervical carcinoma (n = 241), esophageal carcinoma (n = 95), bladder carcinoma (n = 47), and lung squamous carcinoma (LUSC; n = 471) using scalable available whole-exome data (34). TMB of PSCC tumors was the lowest of all squamous tumors analyzed, being significantly lower than any other squamous carcinoma investigated (Fig. 3A). When comparing the TMB of PSCC tumors with other solid malignancies (including nonsquamous tumors from the abovementioned sites) and lymphomas, we found that PSCC has an intermediate TMB among TCGA-sequenced tumors (Supplementary Fig. S3).

Considering the 30 genes with highest average mutation frequency among squamous carcinomas from TCGA, PSCC cohort had a similar mutation frequency in 22 of them (73.3%; Supplementary Table S4). The other eight genes that differentially mutated between PSCC and TCGA squamous cohorts were TP53, TTN, CSMD3, RYR2, KMT2D, FAT1, CDKN2A, and NOTCH1 (Fig. 3B). In addition, PSCC showed mutations in KDM6A, CASP8, LAMA1, DNAH6, and MUC2 genes that are not frequently mutated in the majority of TCGA squamous tumors (Fig. 3B). Patients with PSCC had increased NOTCH1 mutation frequency compared with patients with bladder carcinoma, cervical carcinoma, esophageal carcinoma, and LUSC. Similarly, FAT1 mutants were more frequent in PSCC than in cervical carcinoma and esophageal carcinoma cohorts, CDKN2A was more frequent in PSCC than cervical carcinoma; DNAH6 in PSCC than HNSC, and CASP8 was more frequent in PSCC than esophageal carcinoma and LUSC. TP53, TTN, CSND3, RYR2, and KMT2D were significantly less frequent among PSCC samples. In addition, we compared the mutation frequency of DDR (BRCA1, BRCA2, ARID1A, ATR, CHEK2, PARP1, FANCA, PALB2, and RAD51) genes between PSCC and the squamous cohorts. PSCC cohort exhibited a high number of patients with mutant DDR genes (n = 8, 23.5%), similarly to LUSC (n = 116, 24.6%) and bladder carcinoma (n = 13, 27.7%). HNSC (n = 85, 17%), esophageal carcinoma (n = 13, 13.7%), and cervical carcinoma (n = 44, 18.3%) cohorts showed lower number of DDR-mutant patients, but no cohort was significantly different from PSCC cohort.

To assess the genomic similarity between PSCC and other squamous tumors, we evaluated the frequency of the 96 mutation types in PSCC compared with the squamous TCGA cohorts. By performing cophenetic analysis, we estimated three different mutation patterns among the six squamous cohorts, which we named Sq-MP1, Sq-MP2, and Sq-MP3. As described previously, these mutation patterns were compared with COSMIC signatures, which indicated that Sq-MP1 best match [cosine similarity (CS) = 0.755] was SBS2 (APOBEC cytidine deaminase – C>T), Sq-MP2 best match (CS = 0.879) was SBS6 (defective DNA mismatch repair), and Sq-MP3 best match (CS = 0.947) was SBS4 (exposure to tobacco smoking mutagens; Supplementary Fig. S4). PSCC samples were primarily enriched by Sq-MP2 (n = 23, 67.6%) followed by Sq-MP1 (n = 8, 23.5%) and Sq-MP3 (n = 3, 8%; Fig. 3C). Distribution of each signature among cohorts showed a significant difference between PSCC and LUSC (P < 0.0001), bladder carcinoma (P = 0.006), and cervical carcinoma (P < 0.0001) cohorts. However, no significant difference between PSCC and HNSC (P = 0.208) and esophageal carcinoma (P = 0.387) was observed, suggesting a similar mutation signature between these three cohorts. We further compared the background of mutations among these cohorts by assessing the APOBEC mutation signature estimation, as described previously (Fig. 3D). PSCC samples differed significantly from bladder carcinoma (P = 0.0002) and cervical carcinoma (P < 0.0001) cohorts that exhibited a high enrichment of APOBEC-related mutations. No differences were observed when PSCC was compare with LUSC, esophageal carcinoma, and HNSC cohorts.

Our study provides several novel insights into PSCC biology and highlights the similarities with HNSC, taking advantage of WES analysis, in a rare cancer with limited therapeutic choices (7, 8, 10, 43). The goal was to exploit this technology to gain potential insights into new therapeutic strategies in this rare and lethal cancer. To this extent, while our sample size was small, it represents the largest PSCC cohort reported using the WES approach. We have also analyzed the data in the context of the HPV status and genotype (11, 15, 16). Another limitation to our study is that even though our sample is reflective of the stage at presentation, but only a limited number had advanced stage and no other penile cancer histology beside squamous cell, which limits the generalization of our conclusions. This highlights the importance of ongoing international collaborative efforts aiming to develop the PSCC tumor tissue and blood banking infrastructure that would allow the field to overcome this limitation.

In summary, our systematic analyses of PSCCC genomics identified that the NOTCH signaling deregulation pathway is implicated in more than half of PSCC samples showing an enrichment of Notch pathway (n = 24, 70.6%) alterations. Others have shown enrichment of the Notch pathway at lower rates, this is related to the sequencing technique and analysis used that limited their investigation. This is of interest for future clinical investigation, as recent work showed similar findings with oncogenic mutations found in 67% of human HNSC cases that converge onto the NOTCH signaling pathway. Thus, NOTCH inactivation is a hallmark of both PSCC and HNSC (44). In the era of targeted therapy, a rare cancer like PSCC, could be included with other SCCs umbrella clinical trials targeting a specific common targetable alteration. The above findings provide the potential rationale for the ongoing and future clinical trials in patients with NOTCH1-mutated HNSC using PI3K/mTOR inhibitors, to potentially include patients with PSCC along with HNSC cohorts (44). It is known that the aberrant expression of the essential NOTCH coactivator of the Mastermind-like (MAML) family provides an alternative mechanism to activate NOTCH signaling in human cancers. Interestingly, the work by Necchi and colleagues identified that the MAML2 (HR, 10.411; P = 0.003) was associated with poor OS in patients with penile cancer that received cisplatin chemotherapy (45).

Our data also provide some evidence that PIK3CA and CDKN2A could be driver mutations in PSCC given positional enrichment analysis indicating that mutations in PIK3CA and CDKN2A were significantly enriched in specific protein loci. Future studies should include serial sampling from patients with PSCC or precursors lesions for functional analysis to confirm this hypothesis. Also, a significant proportion of patients (23.5%) was found to have pathogenic or likely pathogenic variants in DDR genes, which could indicate another potential clinical benefit from PARP inhibition alone or in combination with immune checkpoint inhibition strategies (46, 47). These findings could lay the foundation for precision trials using targeted drugs in patients with chemotherapy refractory PSCC.

This study builds on the prior molecular efforts in PSCC, specifically the work published by Ali and colleagues (16), Jacob and colleagues (13), and Feber and colleagues (14), and addresses some of these studies limitations. Our results differed as the sequencing approaches presented by Ali and colleagues and Jacob and colleagues were not comparable with our methodologic approach. In this study, we employed WES, which covers most of the exonic regions in the human genome providing coverage of more than 20,000 human genes. Ali and colleagues and Jacob and colleagues used a target sequencing technology called “comprehensive genomic profiling.” Despite its name, this technology covers alterations in approximately 300 human genes, which represents only 1.5% of the genome coverage we had in our study. In this way, target sequencing technologies are not appropriate to explore overall tumor mutation profile because only a selected fraction of the existent genes is being represented. Beyond covering less of the genome, some of the genes included in such target sequencing approaches are not even fully covered by the assay, underrepresenting the mutation frequency of certain genes.

On the other hand, Feber and colleagues performed whole-exome analysis, however, there are important technical differences between the studies to highlight (14). While the mean sequencing depth generated in our study was of 141×, in Feber and colleagues' study, the mean sequencing depth was only 60×. This difference in the amount of sequencing data generated can have an important impact on the final mutation profile results. In their study, the mean mutation rate per patient was 30, which is much lower than we what found, with a mean of 125.7 mutations per patient. Even adjusting the mutation rate per Mb sequenced, their mutation rate, 1.78 (0.72–7.5) mutations per Mb, was lower than what we found in our study, 1.96 (0.01–21.2). The differences can also be highlighted at the gene level. The most frequently mutated genes in their analysis, FAT1 and TP53, affected only four of the 27 sequenced samples, which represents only 14.8% of their cases. On the other hand, FAT1 and TP53, which were among most commonly mutated genes in our cohort, affected 35% of our samples. These differences indicate that the mutation frequency in their study is probably underestimated. Another possible effect of the low sequencing depth in their study is the presence of false positive findings. Two of their cases (7.4%) showed GPS1 (what they called CSN1) mutations. None of our PSCC cases showed such mutation, and the mutation frequency among the other TCGA squamous tumors was pretty low or absent: bladder carcinoma = 6.4% (3/47), cervical carcinoma = 1.2% (3/241), esophageal carcinoma = no mutations found, HNSC = 0.2% (1/501), and LUSC = 1.3% (6/471). Considering its low frequency of mutations in this gene in squamous tumors, we suspect that the true frequency of GPS1 mutations in PSCC is in the low single-digit range, and it will require a larger cohort of PSCC cases to be confident in the true mutation frequency.

Furthermore, by performing cophenetic analysis, we identified two distinct mutation patterns among PSCC samples. The first (MP1), exhibited the highest cosine similarity to the SBS2 and SBS13 COSMIC signatures, both associated with oncogenic activity of AID/APOBEC family of cytidine deaminases and related to local hypermutation phenomenon called kataegis (C>G – TC context). On the other hand, the second (MP2), was more common and had the highest cosine similarity to SBS6 COSMIC signature, associated with defective DNA mismatch repair system and microsatellite instability. Interestingly, MP1 enrichment was positively correlated with increased TMB (CC = 0.71; P < 0.0001) and correlated with significantly worse survival of patients with PSCC in comparison with those in the MP2 subset. This study is the first to show the existence of a subset of PSCC (38%) with enrichment of APOBEC-related mutation signatures, higher TMB, and worse OS in comparison with the non-APOBEC–enriched subset. These findings can provide the rationale for drug investigation, to assess whether patients with MP1 subset could clinically benefit from immunotherapy with check point inhibitors, while those with the MP2 subset could benefit from PARP inhibitors.

Finally, considering the limited genomic understanding of PSCC, we aimed to compare these tumors with other well-characterized squamous carcinomas using TCGA data. On our initial analysis, we identified that PSCC appears to be genomically distinct in ways from other squamous solid tumors, with a significantly lower TMB and having the highest enrichment of NOTCH1 mutations, but similar to HNSC. When we evaluated the 96 mutation types in PSCC compared with the squamous TCGA cohorts, by performing cophenetic analysis and comparison with COSMIC signatures, PSCC appeared to have a mutation signature that was significantly different from LUSC, bladder carcinoma, and cervical carcinoma cohorts, but similar to the mutation signature of HNSC and esophageal carcinoma. On this note, it is of interest that the first-line chemotherapy (utilized in the InPACT trial) regimen for metastatic penile cancer using cisplatin, paclitaxel, and ifosfamide has been shown to have objective responses in 50% of patients and was initially chosen because of its efficacy in HNSC (6, 48, 49). We further compared the APOBEC mutation signature estimation of PSCC samples, which differed significantly from bladder carcinoma and cervical carcinoma cohorts, but appeared similar to LUSC, esophageal carcinoma, and HNSC cohorts. These findings highlight that the genetic signature of PSCC is most similar to HNSC and, to a lesser degree, to esophageal carcinoma.

Future collaborative efforts examining larger cohorts, coupled with functional molecular analysis and taking benefit of recently available preclinical models (50) are needed for both validation and better characterization of these findings. This work is a first step toward providing the field with the rationale for new therapeutic strategies.

J. Chahoud reported grants from ASCO CCF during the conduct of the study and other from Pfizer outside the submitted work. C.R. Pickering reported grants from NIH outside the submitted work. C.A. Pettaway reported personal fees from Wolters Kluwer outside the submitted work. No disclosures were reported by the other authors.

J. Chahoud: Conceptualization, data curation, formal analysis, investigation, writing–original draft, writing–review and editing. F.O. Gleber-Netto: Conceptualization, data curation, software, formal analysis, methodology, writing–original draft, writing–review and editing. B.Z. McCormick: Data curation, writing–review and editing. P. Rao: Data curation, formal analysis, investigation, writing–review and editing. X. Lu: Formal analysis, methodology, writing–review and editing. M. Guo: Investigation, writing–review and editing. M.B. Morgan: Data curation, software, formal analysis, writing–review and editing. R.A. Chu: Data curation, software, formal analysis, writing–review and editing. M. Martinez-Ferrer: Software, visualization, writing–review and editing. A.K. Eterovic: Resources, software, writing–review and editing. C.R. Pickering: Data curation, software, formal analysis, supervision, methodology, writing–original draft, writing–review and editing. C.A. Pettaway: Conceptualization, data curation, software, formal analysis, supervision, funding acquisition, investigation, methodology, writing–original draft, writing–review and editing.

This work was supported by funding from the University of Texas MD Anderson Cancer Center Human Papillomavirus related cancers moonshot program and the Conquer Cancer Foundation (CCF). J. Chahoud was supported by the CCF American Cancer Society of Clinical Oncology young investigator award.

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