A total of 27 pairs of tissue (tumor and normal adjacent healthy skin) were collected from the Department of Plastic Surgery, St. Josef Hospital, Catholic Clinics of the Ruhr Peninsula, Essen, Germany. While excising the BCC tissues with cold steel under local anesthesia, 4-mm punch biopsies were taken from the center of the tumor and from nonlesional epithelial skin (as normal, intraindividual controls). These samples were immediately placed in RNAlater (Qiagen, Hilden, Germany) and stored at −80°C. Tissue homogenization was performed with stainless steel beads of 5 mm (Qiagen) and TissueLyser LT (Qiagen). DNA was extracted with an AllPrep DNA/RNA/miRNA Universal Kit (Qiagen) according to the manufacturer’s protocol. All samples were quantified using a NanoDrop One (Thermo Scientific, Waltham, USA) and Qubit fluorometer 3.0 (Invitrogen) (Qubit dsDNA HS Assay (Life Technologies, Carlsbad, USA)), and DNA size and quality were tested using gel electrophoresis.

The library was prepared with 200 ng of high-quality DNA using the SureSelectXT Library Prep Kit (Agilent). A SureSelectXT Human All Exon V6 kit (Agilent) was used for exome capture. Sequencing was performed on an Illumina NovaSeq 6000 (San Diego, USA), generating 2x 100 bp paired-end reads. Library preparation, exome enrichment, and sequencing were performed at CeGaT, Tuebingen, Germany. Demultiplexing of the sequencing reads was performed with Illumina bcl2fastq (2.19). Adapters were trimmed with Skewer (version 0.2.2) (34). The Phred score was given with Illumina standard Phred encoding (offset +33). For each sample, two FASTQ files corresponding to forward and reverse reads were obtained. Next steps were done by us on the Poznan Supercomputing and Networking Center (PSNC) Eagle supercomputer. Paired-end reads were aligned to hg38 using BWA. PCR duplicates were marked and removed with the Picard package. Indel realignments with known sites and base quality score recalibration were performed with GATK version 4.1.2.0. SAM to BAM conversion was done using SAMtools. Somatic single-nucleotide variants were called with MuTect2 (version 4.1.0.0. with the use of the tumor-normal mode). Additionally, to avoid false-positive somatic mutations, we performed filtering for germline variants present in the gnomAD database (version 2.1.1). We also generated and flagged variants with a panel of normals (PoN) comprising variants representing commonly occurring sequencing noise that may mimic low allele-fraction somatic variants. We also added information about the localization of mutations in gene subregions (CDS, 5’UTR, 3’UTR, or introns) by use of an in-house Python script. From the list of somatic mutations, we additionally removed those that did not fulfill the following criteria: (i) at least five alternative allele-supporting reads in a tumor sample; (ii) frequency of alternative allele-supporting reads in a tumor sample of at least 0.05; and (iii) frequency of alternative allele-supporting reads in the tumor sample at least 5× higher than that in the corresponding normal sample.

To analyze mutational signatures, we used the web application Mutational Signatures in Cancer [MuSiCa; http://bioinfo.ciberehd.org/GPtoCRC/en/tools.html (35)], allowing the visualization of the somatic mutational profile of each analyzed sample and estimation of the contribution values of the predefined mutational signatures [(36); Catalogue Of Somatic Mutations In Cancer, COSMIC 2020]. Samples BCC14 and BCC21 were excluded from the signature analysis due to an insufficient number of mutations.

We defined genomic positions mutated in at least 3 (>10%) samples as hotspots. Mutations occurring in directly adjacent nucleotides were merged into one hotspot.

We defined genes with nonsynonymous mutations in a coding region in at least 5 samples, with mutations in a 5’UTR, in at least 4 samples, with mutations in a 3’UTR in at least 4 samples, and with mutations in introns (up to 40 nt from exon/intron boundaries) in at least 5 samples as frequently mutated. From the analysis, we excluded genes known to be commonly hypermutated with passenger mutations as a result of the increased background mutation rate but not related to cancer, listed in (37). To distinguish synonymous from nonsynonymous mutations, we used the SnpEff - genetic variant annotation and functional effect prediction toolbox (38), available on the Subio platform (Subio, Inc., Kagoshima, Japan, http://www.subio.jp). We also considered splice-site mutations located in introns up to +/-2 nt from exons as coding region mutations.

OncodriveFML (39) was run using the CADD score (hg38, version 1.6). The signature method was set as a complement, the statistical method was set to “amean”, and indels were included in the analysis using a max method (max_consecutive was set to 7 as default).

To identify chromosome arm-level and focal regions that were significantly amplified or deleted, we used GISTIC2 (40) with the following parameters: threshold for copy number amplifications and deletions, 0.2; confidence level to calculate the region containing a driver, 0.9; broad-level analysis; and the arm peel method to reduce noise.

To validate CNAs involving chromosome 9, i.e., chr9p duplications/amplifications (affecting JAK2, PDL1/CD274, and PDL2/CD273) and chr9q deletions (affecting PTCH1), we designed and generated an MLPA assay covering the entire chromosome 9. In total, the assay consisted of 20 probes, including (i) 7 probes distributed over the chr9p (n=5) and chr9q (n=2) arms, 2 probes located in or in close proximity to JAK2, PDL2, PDL1, and PTCH1 (in total 8 gene-specific probes), and 5 control probes (located on different chromosomes outside of chromosome 9 and regions of known cancer-related genes). The sequences and detailed characteristics of all probes as well as their exact positions are shown in Table S2 .

The MLPA probes and the probe-set layout were designed according to a previously proposed and well-validated strategy (41, 42). Shortly, each probe was composed of two half-probes of equal size, and the total probe length ranged from 93 to 172 nt. The target sequences for the probes were selected to avoid common SNPs, repeat elements, and sequences of extremely high or low GC content. The MLPA probes were synthesized by IDT (Skokie, IL, USA). The MLPA reactions were run according to the manufacturer’s general recommendations (MRC-Holland, Amsterdam, the Netherlands). All reagents except the probe mixes were purchased from MRC-Holland (http://www.mlpa.com). The products of the MLPA reaction were subsequently diluted 20x in HiDi formamide containing GS Liz600, which was used as a DNA sizing standard, and separated via capillary electrophoresis (POP7 polymer) in an ABI Prism 3130XL apparatus (Applied Biosystems, Carlsbad, CA, USA). The obtained electropherograms were analyzed using GeneMarker software v2.4.0 (SoftGenetics, State College, PA, USA). For each individual sample, the signal intensity of each probe was divided by the geometric average signal intensity of the control probes to normalize the run-to-run signal variation, and then the normalized signal of each probe in cancer samples was divided by the corresponding signal in the corresponding normal samples and multiplied by 2. The final MLPA result of each sample is presented on a bar-plot, in which the bars show the relative copy number value of the subsequent probes.

To compare the mutations recurring in BCC with mutations in other cancers, we used WES-generated somatic mutation datasets of 10,369 samples representing 33 cancer types generated and deposited in the TCGA repository (http://cancergenome.nih.gov). The full names and abbreviations of all TCGA cancer types are shown in Table S3 . Somatic mutations were identified against matched normal samples with the use of the standard TCGA pipeline (including the Mutect2, Muse, Varscan, and SomaticSnipper algorithms). We extracted somatic mutation calls (with PASS annotation only) localized in the annotated exons of BAD, DHODH, CHCHD2, FLG, and FLG2 (exon sequences were extended by 2 nt to enable identification of intronic splice-site mutations). The extraction was performed as described in our earlier study (43) with a set of in-house Python scripts available at (https://github.com/martynaut/mirnaome_somatic_mutations).

All mutations were annotated according to HGVS nomenclature (at the transcript and protein levels), and the effects of mutations were defined using the Ensembl Variant Effect Predictor (VEP) tool. For visualization of mutations on gene maps, we used ProteinPaint from St. Jude Children’s Research Hospital – PeCan Data Portal (44). The protein domains visualized on gene maps were positioned according to UniProt data (45). The comutation plot showing frequently mutated genes was created with the use of the Python library CoMut (46).

Target predictions were performed with the TargetScan Custom (release 5.2) web tool (47). The secondary RNA structures were predicted using mfold software (48) with default parameters. RNA sequence/structure functional motifs and transcription factor binding sites were analyzed with the RegRNA 2.0 (49) and MotifMap (50) web tools.

Specific statistical tests are indicated in the text, and a p-value <0.05 was considered significant. If necessary, p-values were corrected for multiple tests with the Benjamini-Hochberg procedure.

We performed WES on 27 paired tumor and corresponding intraindividual control skin DNA samples isolated from 22 nodular and 5 superficial BCC subtypes and corresponding healthy skin tissue. The average coverage of the targeted regions was 183x (185x in normal and 180x in tumor samples), ranging in different samples from 134x to 232x. In total, we identified 84,571 cancer-sample-specific somatic mutations ( Table S4 ), of which 42,380 (50.1%) were located in protein-coding (coding) regions, and the remaining 42,191 (49.9%) were located in noncoding regions ( Table 1 and Figure 1A ). The noncoding regions included (i) 5’UTRs, (ii) ~100 bp fragments of 3’UTRs adjacent to coding sequences (3’UTRs), (iii) exon-adjacent ~100 bp fragments of introns (introns), and (iv) sequences other than those classified above (i-iii), mostly intergenic sequences located upstream and downstream of the first and last gene exons (intergenic regions) (51). The average coverage of the mutated positions was 169x and was slightly higher in coding (195x) than in noncoding regions (142x), whereas the average fraction of reads mapping to alternative alleles was 0.35 (0.33 in coding and 0.40 in noncoding regions). The average mutation rate calculated based on the coding regions was 52.8 mutations/Mbp (ranging from 0.1 to 287.5), which, although slightly lower than that observed before in BCC (15, 52), is still higher than that in any other tested cancer type. Although somewhat counterintuitive, the lower mutation burden in our study than in other BCC studies (15, 52) may result from the much higher sequencing coverage in our study, which gave us much higher statistical power to filter out the fraction of false-positive mutations. The lower mutation burden in our study may also be explained by the identification in our cohort of two samples with an extremely low mutational burden (<0.2 mutations/Mbp). Most of the identified mutations were single-nucleotide substitutions (79,960 (94.5%), predominantly C>T transitions), followed by double substitutions (3,128 (3.7%), predominantly CC>TT transitions) and short (<4 nt) indels (1.483 (1.8%)) ( Table 1 and Figure 1B ). The higher frequency of indels in noncoding regions most likely results from the excess of low complexity sequences, which cause polymerase slippage.

To estimate the fraction of false-positive mutations, we resequenced (with Sanger sequencing) 52 mutations representing different types of alterations, including 39 substitutions and 13 indels ( Table S5 ). The analysis confirmed 51/52 of the mutations, indicating a very low (2%) fraction of false-positive mutations. The fraction may be even lower, as the only unconfirmed mutation (double substitution CC>TT in MYCN) was present in a low fraction of reads (7%), which is generally beyond the sensitivity of Sanger sequencing.

In the next step, we analyzed sample-specific mutational signatures to recognize the mutational processes playing a role in the mutagenesis of the analyzed BCC samples. Shortly, a mutational signature is a frequency pattern for different types of mutations (taking into account direct nucleotide context, -1 and +1 position) characteristic of particular cancer or cancer type. The pattern may reflect a main mutagenic process or a type of DNA repair deficiency that is specific to a given cancer. Originally based on analysis of single nucleotide variants, 30 distinctive mutational signatures were recognized in pancancer (36) but subsequently, the number of specific cancer signatures has been extended taking into account also other types of variants (53). The analysis showed that most of the samples were predominantly associated with signature 7 (average signature contribution (SC) = 0.7) and to a lesser extent with signature 11 (average SC = 0.2) ( Figures 1C, D ). Both signatures consist predominantly of C>T substitutions but differ in the sequence context of the substitutions. Signature 7 is associated with UV irradiation exposure and commonly occurs in melanoma and head and neck cancer. A hallmark of signature 7 is the frequent occurrence of double CC>TT substitutions resulting from UV radiation-induced pyrimidine dimers. Signature 11 was previously found in melanoma and glioblastoma multiforme, often in patients treated with the alkylating agent temozolomide, which is also used in BCC therapy. Only one sample (BCC22) showed a stronger association with signature 11 (SC = 0.6) than signature 7 (SC = 0.3). None of the analyzed samples showed an association with signatures 1, 2, 5, and 13, which are frequent in most cancer types. This may indicate that the deamination of 5-methylcytosine (5meC) predominantly induced by AID/APOBEC cytidine deaminases (attributed to the abovementioned signatures) does not play a role in the pathogenesis of BCC.

The comparison of the nodular and superficial BCC samples showed no substantial difference in terms of mutation burden or mutation types, with the exception of the contribution to mutational signature 7, which was higher for the nodular than superficial samples ( Figure 1E ), consistent with the higher UV radiation exposure of nodular BCCs.

As recurrent mutations may be indicators of the cancer-related function of the mutated genes, we first looked for hotspots defined as genomic positions mutated in at least 3 samples (>10% of the cohort). In total, we identified 43 hotspots, including 23 hotspots in coding and 20 hotspots in noncoding regions (8 in 5’UTRs, 1 in 3’UTRs, and 11 in introns) ( Table S6 ). Of the coding hotspots, 16 resulted in missense mutations, and 7 were synonymous substitutions. As the majority of synonymous mutations result from randomly occurring neutral alterations, we did not analyze the synonymous hotspot further. Although it has to be noted that the functionality of individual synonymous mutations cannot be unequivocally ruled out (51, 54, 55). For example, 315 (~2.1%) of the detected in our study synonymous mutations were predicted to be exonic splice-site mutations. Also, synonymous mutations located inside exons may affect different regulatory elements including exonic splicing enhancers and silencers (55). As shown in Table S6 , some of the hotspots were located in genes annotated in the COSMIC Cancer Gene Census (CGC) database and/or in genes playing a role in cancer or skin function.

Next, we looked at the overall frequency of mutations in the genes, separately analyzing mutations in coding regions, 5’UTRs, 3’UTRs, and introns (defined in Materials and Methods; listed in Table S7 ). Although they were not considered frequently mutated, in this section, we also report genes with any mutations in a coding region if they were detected in a pathway of a recurrently mutated gene. In the analysis of frequently mutated regions, we focused mostly on genes functionally related to cancer (annotated with CGC and a manual literature search) and genes playing a role in skin function.

To further investigate the mutations/mutated genes, we used OncodriveFML, which allows the prediction of the cancer driver potential of both coding and noncoding regions/genes based on functional mutation (FM) bias (39). As shown in Figures 5A–C and Table S8 , we identified 14 potential cancer driver genes based on mutations in coding regions (CDS-drivers), a disproportionately high number of 36 potential cancer driver genes based on mutations in 5’UTRs (5’UTR-drivers), and 7 potential cancer driver genes based on mutations in 3’UTRs (3’UTR-drivers). No potential cancer driver gene was identified based on the mutations in introns.

In addition to 4 CDS-drivers (PTCH1, TP53, TGFB1I1, and CARD6) also identified as frequently mutated, it is worth noting RORA, recently shown to play an important role in restraining allergic skin inflammation (103). Other interesting genes were PRDM9 and ZNF281, both of which play a role in DNA repair and have been shown to be responsible for frequent mutations in cancer (104, 105). None of these genes were previously implicated or identified as frequently mutated in BCC.

Among the 5’UTR-drivers, 6 were also identified as frequently mutated: DHODH, CHCHD2, and SPHK2 (described above), as well as POLR2M, NPC1, and NELL2. Additionally, it is worth noting IKBKB (mutated in 3 samples but not reported before as mutated in BCC) shown to act as a tumor suppressor in nonmelanoma skin cancers and noncancerous skin lesions; it was also shown that deletions of the gene lead to skin inflammation, hair follicle disruption, hyperplasia, and SCC development (106–109).

Among 3’UTR-drivers, two genes (mentioned above), i.e., BAD (the most significant 3’UTR-driver) and SMIM27 were also identified as frequently mutated. Additionally, it is worth mentioning the transcription factor gene POU3F2 (mutated in 3 samples), that plays a role in the invasiveness and metastasis of melanoma, and is controlled by miR-211 (110, 111) and miR-107 (112). Although the mutations were not located in the predicted miR-107 and miR-211 binding sites, they may affect the structure of the 3’UTR and thus indirectly change accessibility to these or other miRNA targets.

As somatic CNAs have not been extensively studied in BCC, in the next step, we performed analysis of both chromosome arm-level and focal CNAs [with GISTIC2 (40)]. At the chromosome arm level, we detected a significant recurring deletion of chr9q (q=1.4x10-6; occurring in 9 samples), involving PTCH1 ( Figures 4 , 6 ), and a significant recurring amplification of chr9p (q=0.05; occurring in 5 samples), involving a region with CD274 (also known as PDL1, encoding PD-L1), CD273 (also known as PDL2, encoding PD-L2), and JAK2 ( Figures 4 , 6 ). Although the loss of chr9q has been frequently observed in BCC (reported as loss-of-heterozygosity of PTCH1), gain of chr9p has been reported only in one case of rare metastatic BCC (113). To validate the chromosome 9 CNAs, we developed a multiplex ligation-dependent probe amplification (MLPA) assay with probes covering the entire chromosome 9 but especially focusing on the region containing PTCH1 (chr9q22.32) and the region harboring PDL1, PDL2, and JAK2 (chr9p24.1) ( Figure 6 ). The MLPA analysis confirmed CNAs in all tested samples as detected by GISTIC2, and examples are shown in Figure 6 .

CNA analysis also showed 54 regions of significant focal deletions, including 27 regions containing skin/cancer-related genes, and 56 significant amplifications, including 20 encompassing skin/cancer-related genes ( Figure 6 and Table S9 ). The elements involved in the most significant focal deletions were CDK11A (chr1p36.33; q=2.4x10-5; occurring in 6 samples), whose loss induces skin carcinogenesis (114); the LCE cluster (chr1q21.3; q=2.4x10-6; occurring in 4 samples), including genes such as LCE2 and LCE3, which play a role in maintaining skin barrier function and whose deletion has been associated with psoriasis (115); and the HLA-D cluster (HLA-DP, -DQ, and -DR, chr6p21.32; q=2x10-4; occurring in 3 samples), encoding components of major histocompatibility complex (MHC) class II molecules, whose increased expression has been associated with increased cancer immunogenicity and better prognosis in BCC, SCC and melanoma (116–122). The skin/cancer-related genes in the most significant focally amplified regions worth mentioning are STIM2 (chr4p15.2; q=0.16; occurring in 2 samples) (123), KLRB1/CD161 (chr12p13.31; q=0.007; occurring in 2 samples) (124, 125), and SPTLC3 (chr20p12.1 q=0.23; occurring in 2 samples) (126).