The Genomic Landscape of Actinic Keratosis

Actinic keratoses (AKs) are lesions of epidermal keratinocyte dysplasia and are precursors for invasive cutaneous squamous cell carcinoma (cSCC). Identifying the specific genomic alterations driving the progression from normal skin to skin with AK to skin with invasive cSCC is challenging because of the massive UVR-induced mutational burden characteristic at all stages of this progression. In this study, we report the largest AK whole-exome sequencing study to date and perform a mutational signature and candidate driver gene analysis on these lesions. We demonstrate in 37 AKs from both immunosuppressed and immunocompetent patients that there are significant similarities between AKs and cSCC in terms of mutational burden, copy number alterations, mutational signatures, and patterns of driver gene mutations. We identify 44 significantly mutated AK driver genes and confirm that these genes are similarly altered in cSCC. We identify azathioprine mutational signature in all AKs from patients exposed to the drug, providing further evidence for its role in keratinocyte carcinogenesis. cSCCs differ from AKs in having higher levels of intrasample heterogeneity. Alterations in signaling pathways also differ, with immune-related signaling and TGFβ signaling significantly more mutated in cSCC. Integrating our findings with independent gene expression datasets confirms that dysregulated TGFβ signaling may represent an important event in AK‒cSCC progression.


INTRODUCTION
Actinic keratoses (AKs) are dysplastic epidermal keratinocyte lesions resulting from chronic UVR exposure. It is generally accepted that AKs are cutaneous squamous cell carcinoma (cSCC) premalignancies (Siegel et al., 2016): at least twothirds of cSCC arise from AK, although fewer than 0.6% AK per year will progress to cSCC (Criscione et al., 2009).
Identifying genomic alterations driving AK-cSCC progression is challenging because of the high level of background mutations in keratinocytes associated with UV exposure. cSCCs have an average tumor mutational burden (TMB) of 50 mutations per megabase DNA, making them among the most mutated of all human cancers (South et al., 2014).
Epigenomic alterations and immunological factors (reviewed in Bottomley et al. [2019]) add further layers of complexity in characterizing this progression.
Several molecular genetic studies have examined AK and cSCC at the levels of gene expression, chromosomal instability, and mutations in known cancer genes. These studies generally show that AK and cSCC possess similar genetic alterations with some conflict on whether AKs are more or less genomically unstable (summarized in Supplementary  Table S1). However, most AK-cSCC genetic studies have used targeted techniques on fixed tissue, which have many limitations and biases. Only three previous studies have used whole-exome sequencing (WES) to study AK; these studies have included 10 samples and reported similar mutational spectra in AK and cSCC with frequent mutations in known cSCC tumor suppressor genes TP53, NOTCH1, NOTCH2, and FAT1 (Albibas et al., 2018;Chitsazzadeh et al., 2016;Rodríguez-Paredes et al., 2018) (Supplementary Table S2).
In this study, we present the largest WES study of AK conducted to date and include AK from both immunocompetent (IC) and immunosuppressed (IS) patients. We demonstrate that AKs are strikingly similar to cSCC at the genomic level with similar patterns of driver genes and copy number alterations. We also demonstrate mutational signature 32 in all samples from patients exposed to azathioprine, providing further evidence for its role in keratinocyte carcinogenesis at the precursor stage. By integrating our mutational pathway analysis with independent gene expression datasets, we have identified that the dysregulation of TGFb signaling may represent a critical gatekeeper pathway resulting in AK-cSCC progression.

Patients and samples
A total of 37 AKs from 37 patients (median age ¼ 70 years, range ¼ 48e86 years) were included (Supplementary Table S3): 23 AKs from IS patients and 14 AKs from IC patients. WES data from cSCC samples previously analyzed for 16 of these patients were also available for comparison (Inman et al., 2018).

Mutational burden and mutational signatures
WES targeted 334,378 exons from 20,965 genes and yielded a mean coverage of Â53 with 83% of targeted bases covered by Â20 and 94% of them covered by Â10. In total, 80,511 somatic mutations (range ¼ 12e8,326 per AK) were identified. Of these, 49,853 were nonsynonymous (range ¼ 9e4,993 per AK) with a median of 1,676 total and 1,071 nonsynonymous mutations per AK (Figure 1a and Supplementary Tables S4 and S5). This corresponds to a mean TMB of 43.5 mutations per megabase DNA, similar to the TMB of 50 mutations per megabase DNA we previously reported in cSCC (Inman et al., 2018).
AKs from IS patients were significantly more mutated than those from IC patients when comparing total nonsynonymous mutations (median ¼ 1,609 vs. 729 mutations, respectively; Wilcoxon test, P ¼ 0.03; Figure 1b). The results remained significant when controlled for per sample read depth (Supplementary Table S5). AKs from IS patients had a significantly higher median TMB than those from IC patients (55.9 and 22.3 mutations per megabase DNA, respectively; Wilcoxon test, P ¼ 0.03). There was no significant difference in nonsynonymous mutation rates from AK compared with those from cSCC (median ¼ 1,070 vs. 912; Wilcoxon test, Nine single-base substitution mutational signatures (Alexandrov et al., 2020(Alexandrov et al., , 2013 were identified (Figure 1c). A total of 33 AKs (89%) displayed characteristic UVR signatures (7a/b, 13e97%), mainly C > T transitions at dipyrimidine sites. A total of 20 AKs (54%) harbored variable levels (5e100%) of signature-5 mutations, and 35 AKs (95%) had low levels of signature-1 mutations (1e13%). The etiologies of signatures 5 and 1 are unknown, but both are commonly found at low levels in most cancers and are clock-like because the number of mutations correlates with age.
Signature 32 was detected in 22 AKs (59%), exclusively in patients exposed to azathioprine ( Figure 1c) (Fisher's exact test, P < 0.00001), reproducing our findings in cSCC. In contrast to cSCC, the prevalence was not significantly associated with the duration of azathioprine exposure (r ¼ 0.21, Supplementary Table S6).

Somatic copy number aberrations
Copy number aberrations (CNAs) and loss of heterozygosity (LOH) events were also examined (Figure 2 and  Supplementary Table S7). Across our AK cohort, a median of 9% (range ¼ 0e62%) of each AK genome was affected by CNAs (Supplementary Table S8), similar to cSCC (Wilcoxon test, P ¼ 0.68) and independent of immune status (Wilcoxon test, P ¼ 0.26). There was no significant correlation between mutational burden and the proportion of the genome affected by CNAs (r ¼ 0.25, P ¼ 0.13).
GISTIC (Mermel et al., 2011) analysis identified 19 significantly deleted regions (q < 0.25) encompassing 859 genes, including 27 known cancer genes (Catalogue Of Somatic Mutations In Cancer database [Tate et al., 2019]). 9p21.3 was the most significantly deleted segment in AK and harbors 10 cancer genes, including CDKN2A, JAK2, and MLLT3 (Supplementary Figure S1a and Supplementary Table S10). Nine regions were significantly gained (q < 0.25) encompassing 59 genes, none known to be cancer related (Supplementary Figure S1b and Supplementary Table S11). Analysis in IC and IS subgroups did not resolve any major differences because of the limited power within each group (Supplementary Table S12).
IS and IC subgroup analysis using our SMG pipeline identified 23 SMGs (detected by both OncodriveFM and MutsigCV, P < 0.05 cutoff was implemented because of smaller sample sizes) in the IS group and 10 SMGs in the IC group (Supplementary Tables S16 and S17). Onco-driveCLUST was unsuccessful for this analysis. TP53, NOTCH1, and NOTCH2 were the only SMGs common to IS, IC, and the combined cohort.

Comparison between AK and cSCC SMGs
We integrated somatic mutations and CNAs to produce a mutation OncoPrint for the 44 AK SMGs ( Figure 3b) and compared this with the OncoPrint for the 22 cSCC SMGs from our previous series. TP53, NOTCH1, and NOTCH2 were the only shared SMGs, and 6 of the total 63 SMGs showed differed frequencies between AK and cSCC (Supplementary Table S18). CCDC71L (chi-square test, P ¼ 0.004) and STRN (P ¼ 0.014) were altered at higher frequency in AK, whereas LCLAT1 (P ¼ 0.032), HERC6 (P ¼ 0.029), MAP3K9 (P ¼ 0.043), and TMEM51 (P ¼ 0.048) were altered at a higher frequency in cSCC. When adjusted for multiple comparisons, none of the 63 SMGs were altered significantly between AK and cSCC.
We also compared AK SMGs with well-differentiated and moderately and/or poorly differentiated group-specific SMGs (Inman et al., 2018). There were no significant differences between SMGs in AK and those in well-differentiated cSCC, but 10 moderately and/or poorly differentiated cSCC SMGs were significantly more altered in cSCC than in AK (PRB1, TMEM51, LRP1, POLH, ACVR2A, RPLP1, ZZEF1, VWF, ICAM1, and HECTD4) (Supplementary Table S18). The mutation distribution along the protein domains was similar for   Figure S2), with the exception of the PIK3CA hotspot in AK, as described earlier.

AK and cSCC SMGs (Supplementary
Genes previously implicated in cSCC, such as CDKN2A and HRAS, were found altered at similar rates between our AK and cSCC cohorts (54.1% and 45% altered for CDKN2A, respectively, and 16.2% and 22.5% altered for HRAS, respectively; Supplementary Table S18).
We next compared the cancer cell fractions of nonsynonymous mutations between AK and cSCC across the 63 SMGs to assess their clonality and to identify whether mutations in any genes are more clonally dominant in AK or cSCC using synonymous and untranslated region mutations as a negative control. Two genes (CACNA1C and KCNK5) showed significant evidence that nonsynonymous mutations in them were more clonally dominant in AK than in cSCC, suggesting that these genes may be under direct selection in the AK lineage (Supplementary Table S19 and Supplementary Figure S3).

Inference of clonal evolution and order of genetic changes
To assess the levels of intrasample heterogeneity, we performed clonality analyses using Expanding Ploidy and Allelefrequency on Nested Subpopulations (Andor et al., 2014), which estimates the number of clones and their proportions within each sample. A total of 31 AKs (84%) had sufficient numbers of somatic mutations (!200) for analysis. Variation in intrasample heterogeneity was identified, with a median of four (range ¼ 1e9) clones per AK ( Figure 4a) compared with a median of six in cSCC (Wilcoxon test, P < 0.01, Figure 4b).
For each of the 44 SMGs, we assessed the order of mutation acquisition inferred from the aggregate frequencies at which they were found to be clonal or subclonal. More than 70% of mutations in 22 SMGs (TP53, NOTCH1, FAT1, ADAMTS9, CCDC71, RB1, PBRM1, ADAM19, CHEK2, KIF13B, DAB2, ABI3BP, GPR161, XAB2, GPS1, IMPA1, USP34, GXYLT1, EPB41L2, PAX3, KCNK5, and INSIG2) were clonal, indicating that alterations in these genes are more likely to be early events in AK development. Compared with the 22 SMGs mentioned earlier, a relatively larger proportion of subclonal nonsynonymous mutations, implying later development, were found in the remaining 22 SMGs, although more clonal than subclonal mutations were still observed for most of them, except for PIK3CA, RAB22A, and KIF24 ( Figure 4c). For the three SMGs shared by both AK and cSCC, patterns of clonality were similar, with a larger proportion of clonal mutations observed in TP53 and NOTCH1 than in NOTCH2.
Seven AKs had more than one nonsynonymous mutation in TP53, and we estimated the timing of multiple mutations in these lesions. In four samples (AK05, AK33, AK34, and AK35), all TP53 mutations appeared to be clonal and thus likely early events on the basis of clonality analysis. Two samples (AK04 and AK32) showed a mixture of early and late (i.e., subclonal) events, and one sample (AK38) exhibited both TP53 mutations as late events (Supplementary Figure S4aed). These data suggest that although TP53 mutations are commonly observed as early events, additional mutations can occur later on in the development of AK.

Comparison of intrapatient AK with cSCC
A total of 16 AKs (43%) also had cSCC WES from the same individual (mostly from anatomically distinct sun-exposed fields) (Supplementary Table S20). There was a significant positive correlation in the mutational signature profiles of   Figure S5a and b). Seven AK-cSCC pairs had overlapping CNA segments involving a median of 20% of total CNA segments, most commonly present on chromosomes 3, 9, and 17. The direction of CNAs was largely concordant, although not statistically significant (chisquare test, P ¼ 0.06) (Supplementary Table S21).

Mutational pathway analysis
We compared significantly mutated pathways between AK and cSCC (Figure 5a and b). Many immune system signaling pathways were significantly more mutated in cSCC than in AK (TCR, Fc epsilon-RI, RIG-I-like receptor, and chemokine signaling). TGFb, adipocytokine, GnRH, and insulin signaling were also significantly more mutated in cSCC ( Figure 5a). Several metabolism pathways were differentially mutated: ether lipid, pyruvate, or a-linolenic among others were significantly more mutated in cSCC (Figure 5b).

Integration of genomic drivers, CNAs, and gene expression profiles
We integrated the existing gene expression profiles of normal skin, AK, and cSCC to investigate the potential tumorsuppressing or -promoting roles of AK SMGs through the pattern of their expression in disease progression. We used five independent datasets (described in the Materials and Methods section). The SMG hierarchies across the datasets were broadly concordant with two main clusters observed: one consisting of SMGs upregulated in normal skin and  AK29  AK26  AK33  AK04  AK40  AK17  AK14  AK15  AK02  AK34  AK21  AK35  AK23  AK32  AK43  AK12  AK10  AK31  AK38  AK41  AK05  AK06  AK24  AK44  AK45  AK09  AK25  AK30  AK42  AK37  AK03  AK08  AK39  AK28  AK22  AK07    progressively downregulated in AK and cSCC, and the other cluster demonstrating the reverse pattern ( Figure 6a). In AK versus normal skin, five SMGs were significantly downregulated in at least two datasets (KIF24, KCNK5, EPB41L, INSIG2, and ABI3BP), suggesting a tumor suppressor role ( Figure 6 and Supplementary Figure S6a). NOTCH1 was the only significantly upregulated SMG (Figure 6b). In cSCC versus AK, IMPA1 was the sole SMG to be significantly differentially expressed in at least two datasets, being upregulated in cSCC, suggesting a tumor promoter role (Figure 6c and Supplementary Figure S6b).
Because the TGFb signaling pathway was significantly more mutated in cSCC than in AK, we examined the expression patterns of TGFb pathway genes in normal skin, in skin with AK, and in skin with cSCC. The datasets show two clusters of progression-dependent dysregulation of TGFb signaling (Figure 6d). One cluster showed an upregulation of genes in the normal skin that become increasingly more downregulated, progressing from AK to cSCC, and the second cluster had the reverse pattern (Figure 6d).
We assessed the gene expression of significantly deleted and gained regions identified by our GISTIC analysis, which   revealed 55 genes in the deleted regions that were downregulated (in at least two datasets) and two genes in the gained regions that were upregulated (in at least two datasets) in AK compared with those in the normal control (Supplementary Figure S7aec).

DISCUSSION
This study, the largest AK genomic dataset to date, demonstrates that AKs and cSCC are strikingly similar at the genomic level in terms of TMB, patterns of driver genes, and CNAs. The AK TMB of 43.5 mutations per megabase DNA is higher than previously reported (Albibas et al., 2018;Chitsazzadeh et al., 2016), which is likely a consequence of having a majority of the AKs from IS individuals because we have demonstrated immunosuppression results in significantly higher mutation rates. The TMB from the AKs of IC patients (30.4) is similar to the 34.5 reported previously (Albibas et al., 2018). We detected the expected UV signature (7a/b) in most AKs. Furthermore, signature 32 was present in all AKs from patients exposed to azathioprine, further implicating this drug in the early stages of cSCC development (Inman et al., 2018). The significant positive correlation in the matched AK-cSCC mutation profiles also provides further evidence that the same underlying mutational processes are at play in both AK and cSCC even when they are from different anatomical sites.
We identified 44 AK SMGs, including many classical tumor suppressor genes (TP53, NOTCH1, NOTCH2, and FAT1), that are consistently mutated in cSCC, but importantly, these genes are also mutated at low levels and under strong positive selection in physiologically normal sun-exposed skin (Martincorena et al., 2015). However, CDKN2A is not mutated in the normal skin, and we previously postulated that it may have a gatekeeper role in cSCC (Inman et al., 2018). CDKN2A was not identified as an SMG, but the loss of 9p21.3-a CDKN2A gene locus-was identified as the most significantly deleted CNA, and there was no significant difference in the frequency of deletion in AK compared with that in cSCC (54% vs. 45%, respectively, chisquare test, P ¼ 0.43, Supplementary Table S18), suggesting that this deletion is an early event and may play an important role in AK pathogenesis. The oncogene PIK3CA was significantly altered in almost half of all AKs and cSCC, a higher frequency than seen in previous cSCC studies (Hafner et al., 2010;Janus et al., 2017;Pickering et al., 2014). We also identified a hotspot activating splice-site mutation in three AKs. Activating mutations in PIK3CA result in the activation of the phosphoinositide 3-kinases/protein kinase B/mTOR pathway, which is commonly seen in other organs' squamous cell carcinomas (Vivanco and Sawyers, 2002).
J Thomson et al. HMCN1 was identified as an SMG, being altered in 19 AKs (51%). Although not a confirmed cancer gene, its role warrants further investigation because it is frequently mutated in cSCC (Durinck et al., 2011;Pickering et al., 2014) and AKs (Albibas et al., 2018). It is postulated to be involved in cancer cell invasion and metastasis (Timpl et al., 2003) through its function as an extracellular matrix protein.
We show that AKs have the same levels of genomic instability as cSCC with many shared CNAs between them, notably loss of 9p, chromosome 13, and 5q. Apart from an early study that suggested that AKs have more LOH than cSCC (Rehman et al., 1996), our findings are contrary to those of more recent studies that found that AKs have relatively few CNAs (García-Díez et al., 2019;Hameetman et al., 2013). However, most studies agree on the most common sites of chromosomal instability, with 9p loss-the region where CDKN2A is located-ubiquitous among AK studies (Ashton et al., 2003;García-Díez et al., 2019;Kanellou et al., 2008;Rehman et al., 1996).
Interrogating published expression datasets against our SMG list found KIF24, KCNK5, EPB41L2, and ABI3BP downregulated in AK compared with those in the normal skin, suggesting tumor suppressor roles. ABI3BP was significantly downregulated in cSCC compared with that in the normal skin in one previous study (Prasad et al., 2014) and similarly downregulated in the esophagus with squamous cell carcinoma compared with that in the normal esophagus (Jiang et al., 2018). It is a tumor suppressor in several cancers where it functions to promote cellular senescence (Wakoh et al., 2009) and has a role in cell-substrate adhesion (Hodgkinson et al., 2013;Latini et al., 2008). KCNK5, a twopore domain potassium channel, is in the top 1% of underexpressed genes in melanoma and top 5% of underexpressed genes in breast, colorectal, renal, and liver cancers, and there is increasing interest in the role of potassium channels in cancer (Williams et al., 2013). NOTCH1 mRNA was significantly upregulated in AK compared with that in the normal skin (in two datasets), suggesting tumor promoter function, which contrasts with previous findings in cSCC where it is inactivated early in cSCC pathogenesis (South et al., 2014). We also observe early mutational inactivation of NOTCH1 in AK, and subsequent loss of expression may facilitate progression from AK to cSCC. NOTCH proteins have opposing tumor suppressor and promoter roles in different cancer types and contexts (Ranganathan et al., 2011), and this requires further investigation in AK to cSCC progression.
We have previously shown that the dysregulation of TGFb signaling through inactivation of its receptors in stem cells is an early driver event in cSCC pathogenesis and plays a likely tumor suppressor role (Cammareri et al., 2016). TGFb signaling pathway genes were significantly more mutated in cSCC and were also dysregulated in expression datasets. From the GSE45216 dataset (the largest dataset), TGFBR2 becomes progressively more underexpressed in the transition from normal skin through AK to cSCC, consistent with the previous studies that have demonstrated a key role for TGFBR2 in TGFb tumor-suppressive function (Han et al., 2005). Taken together, these findings further support the hypothesis that TGFb dysregulation is a critical step in AK-cSCC transition.
Epigenomic alterations may also play a part in driving AK progression, and this is supported by the significant differences we have observed in AK and cSCC expression profiles. Recent work landscaping the methylomes of AK and cSCC has directly addressed this hypothesis. A complex nonlinear evolution of distinct DNA-methylation patterns during the progression of AK to cSCC and metastasis has been demonstrated by one group (Hervás-Marín et al., 2019), but others failed to show any differences in AK and cSCC methylomes (Rodríguez-Paredes et al., 2018). Epigenomic studies to date are limited, and further research is needed.
A limitation of our study is that although we clinically diagnosed and then histologically confirmed AK, we did not specifically grade AK and analyze them according to grades, that is, AK-I, -II, and -III. There is a possibility that each of these AK grades might have a different molecular profile. However, AKs are frequently heterogeneous histologically and may include combinations of AK IeIII. Consequently, there is a possibility that small foci of AK III and/or cSCC in situ were included in these samples, and this may have affected the results. Laser capture microdissection of samples may have helped to minimize this risk.
In conclusion, our data demonstrate that AKs already possess the majority of genomic alterations present in cSCC. Significant molecular alterations, which we have uncovered and which may contribute to evolution from AK to cSCC, include alterations in key signaling pathways, particularly in TGFb and immune system signaling; mutations in specific genes, including ABI3BP and IMPA1; and differences in intrasample heterogeneity. These will be the focus of future investigations of the molecular pathogenesis of AK progression.

Collection of patient samples and clinical data
The 4 mm punch biopsies of AK diagnosed clinically and confirmed histologically were obtained from participants after written informed consent was obtained. Processing of AK and venous blood for matched germline DNA and clinical data collection was done as previously described (Inman et al., 2018). The study was approved by the East of Scotland Research Ethics Service (reference: 08/S1401/ 69) and the East London and City Health Authority Local Ethics Committee and conducted according to the Declaration of Helsinki Principles.

DNA extraction and genetic analysis
DNA extraction and genetic analysis were performed as previously described (Inman et al., 2018).
WES data processing, somatic variant calling, annotation, and visualization WES data were analyzed and annotated using our established pipeline (Inman et al., 2018). Maftools was used to summarize, visualize, and compare AK and cSCC Mutation Annotation Format files and to make lollipop plots of gene mutation distributions (Mayakonda et al., 2018). Mutation signatures across the exomes were identified on the basis of the nonnegative matrix factorization approach previously described (Alexandrov et al., 2020(Alexandrov et al., , 2013 using the Catalogue Of Somatic Mutations In Cancer mutational signatures.

Identification of CNA and LOH using WES data
Analyses for CNA and LOH events from WES data were based on a combinational approach previously described (Okosun et al., 2016;Tawana et al., 2015). Genes targeted by copy-gain, copy-loss, and copy-neutral LOH in each sample were identified. Thresholds for CNA calls were fully described in our previous study (Inman et al., 2018).

Tumor subpopulation identification and clonality analysis
Expanding Ploidy and Allele-frequency on Nested Subpopulations (Andor et al., 2014) was used to estimate clonal expansions and cellular frequency of each clonal and subclonal population as previously described (Inman et al., 2018). Tumor purity was also estimated on the basis of the cellular frequency of the dominant clone. All somatic variants were assigned to their nested clones.