Dnmt3a deficiency in the skin causes focal, canonical DNA hypomethylation and a cellular proliferation phenotype.

Mice with homozygous Dnmt3a deficiency develop normally but die within 3 wk of birth ( 28 ). To study the effects of Dnmt3a deficiency in the skin, we therefore used transgenic mice expressing Cre recombinase under the control of the Keratin 14 ( Krt14 ) promoter, which exhibits a well-established pattern of expression that affects the epidermal skin and its appendages ( 29 ). With appropriate crosses, we obtained mice that were wild type (WT) for Dnmt3a ( Krt14 - Cre − × Dnmt3a flox/flox , hereafter “ Dnmt3a WT ”), heterozygous ( Krt14 - Cre + × Dnmt3a flox/+ , hereafter “ Dnmt3a HET ”), or homozygous for the floxed allele ( Krt14 - Cre + × Dnmt3a flox/flox , hereafter “ Dnmt3a KO ”). The floxing efficiency of the Dnmt3a flox/flox alleles in epidermal cells with Krt14 -Cre was >90%, resulting in a complete loss of Dnmt3a protein expression ( SI Appendix , Fig. S1 A ). Whole genome bisulfite sequencing (WGBS) was used to assess DNA methylation of isolated murine epidermal cells. Dnmt3a KO mice exhibited a small but significant decrease in global DNA methylation in all annotated regions examined, including CpG islands, island shores, island shelves, gene bodies, promoters, enhancers, and transcriptional start sites ( Fig. 1 A , Top ). Loss of one allele of Dnmt3a ( Dnmt3a HET ) resulted in a statistically significant change in methylation from the Dnmt3a WT in all regions except gene bodies. On a global scale, the overall decrease in CpG methylation was very small, similar to what has been observed in hematopoietic tissues ( Fig. 1 B , Top ) ( 18 , 30 ). Dnmt3a deficiency results in focal, canonical CpG methylation loss in murine epidermal cells. ( A ) WGBS of isolated murine epidermal cells, comparing independent samples from 11 Dnmt3a WT vs. 5 Dnmt3a KO vs. 3 Dnmt3a HET mice. Dot plots represent average methylation values over annotated regions for all CpGs ( Top ) and for DMRs within each annotated region ( Bottom ). Significance was determined by Student’s t test with denotation of “ns” for not significant and levels of significance as * adjusted P < 0.05, ** P < 0.01, *** P < 0.001. ( B ) Methylation density plots for the indicated genotypes across all genomic CpGs ( Top ) and at DMRs ( Bottom ). ( C ) Targeted deficiency of Dnmt3a ( Krt14 -Cre + × Dnmt3a fl/fl , Dnmt3a KO ) results in 23,983 DMRs in epidermal cells compared to wild-type control, Dnmt3a WT . Heatmap representation of WGBS data demonstrates DMRs defined by differences between Dnmt3a WT vs. Dnmt3a KO epidermal cells. Heterozygosity for Dnmt3a ( Krt14 -Cre + × Dnmt3a fl/+ , Dnmt3a HET ) results in an intermediate phenotype. Each column represents an independent mouse. ( D ) Annotated regions associated with DMRs. *Of 16,315 DMRs involving gene bodies, 3,612 unique genes are involved. ( E ) Integrated genome viewer (IGV) view of the G0s2 gene CpG methylation of Dnmt3a WT versus Dnmt3a KO epidermal cells. Normalized methylated CpG reads are represented in bars (scale 0 to 1.0 for each row). Differentially methylated regions are indicated in the bottom track and highlighted in yellow. We identified differentially methylated regions (DMRs) between Dnmt3a WT and Dnmt3a KO samples using previously described approaches ( 18 , 30 ). DMRs were defined as having greater than 10 CpGs with a mean methylation difference of more than 0.2 (i.e., 20%) between Dnmt3a WT and Dnmt3a KO skin, with a false discovery rate (FDR) of less than 0.05. Adjacent DMRs within 50 base pairs were merged. Based on these criteria, we identified 23,983 DMRs with an average size of 0.95 kilobases (Kb), encompassing 23.53 megabases (Mb) of DNA, which represents 0.87% of the genome ( Dataset S1 ). Virtually all DMRs were hypomethylated (23,981/23,983; 99.99%), and the degree of hypomethylation did not vary based on their location in annotated genomic regions ( Fig. 1 A , Bottom ). Skin from Dnmt3a HET mice demonstrated a more subtle methylation phenotype at these DMRs ( Fig. 1 B , Bottom ). The DMRs in Dnmt3a-deficient skin are represented by a heatmap, where each column represents an independent mouse, and each row, the mean methylation value of an individual DMR for all samples ( Fig. 1 C ). There is remarkable consistency (i.e., canonicality) among the methylation phenotypes at each DMR among independent mice with the same genotypes. The methylation values of the Dnmt3a HET samples are plotted passively and exhibit an intermediate level of methylation at many DMRs ( Fig. 1 C ), demonstrating that haploinsufficiency for Dnmt3a also has a methylation phenotype in the skin ( Fig. 1 A and B ). Dnmt3a-dependent DMRs occur most commonly within gene bodies; 16.4% of all annotated protein-coding genes were found to have at least one DMR. DNA methylation within most CpG islands is low in Dnmt3a WT skin, as expected. A small fraction of these islands (142/16,923; 0.89%) exhibit significantly higher methylation than total CpG islands in Dnmt3a WT skin, and these rare islands generally become hypomethylated in Dnmt3a KO skin ( Fig. 1 D and A , Bottom ). These DMRs usually represent intragenic CpG islands, of which ∼65% are methylated in human embryonic stem cells ( 31 ). Indeed, 129/142 (91%) Dnmt3a KO DMRs associated with CpG islands are located within gene bodies. The shores and shelves flanking CpG islands exhibited greater proportions of DMRs (3.7% and 8.85%, respectively), while the methylation of transcriptional start sites (TSSs) was minimally affected by Dnmt3a deficiency (0.86%, Fig. 1 D ). The methylation pattern near the G0s2 gene is shown as an example of a Dnmt3a KO DMR that occurs in a gene body ( Fig. 1 E ). This gene was found to be up-regulated in Dnmt3a KO skin tumors ( 11 ). Dnmt3a deficiency does not result in any gross abnormalities in unperturbed skin ( 11 ). Further, flow cytometric evaluation of cellular populations within the skin revealed minimal differences between Dnmt3a WT and Dnmt3a KO ( SI Appendix , Fig. S1 B ). To interrogate the transcriptomic consequences of Dnmt3a deficiency, we employed single-cell RNA sequencing to analyze adult murine epidermis in telogen (the resting phase of the hair cycle). We analyzed two mice from each genotype ( Dnmt3a WT and Dnmt3a KO ) at 9 to 10 wk of age; 10,482 cells from Dnmt3a WT and 15,403 cells from Dnmt3a KO met quality standards. Unbiased clustering based on expression signatures defined subpopulations of cells in the skin ( Fig. 2 A ) that were remarkably similar among all four samples ( Fig. 2 B ). Single-cell RNA sequencing reveals that Dnmt3a deficiency causes discrete gene expression changes in murine keratinocytes. ( A ) UMAP representation of single-cell RNA sequencing data from whole epidermis of two pairs of Krt14 -Cre − × Dnmt3a fl/fl ( Dnmt3a WT ; 10,783 cells) and Krt14 -Cre + × Dnmt3a fl/fl mice ( Dnmt3a KO ;16,324 cells) with unbiased graph-based clustering demonstrating known skin populations represented by lineage defining genes ( Right ), including HFB, uHF, bIFE, sbIFE, and sebaceous gland (Seb), Krt18-expressing cells (Krt18), proliferative cells (Pro), melanocytes (Mel), T cells (T-cell), and Langerhan’s cells (Lang). ( B ) Comparison of keratinocyte lineage cells in the Dnmt3a WT versus Dnmt3a KO mice, excluding inflammatory cells and melanocytes (excluded clusters 6, 8, 9, and 10). UMAP projections for each pair of samples with the indicated genotypes are shown. ( C ) No major population shifts are observed in epidermal lineages affected by Krt14 -Cre mediated Dnmt3a deficiency, with the exception that the cells dominated by a proliferative signature (Pro, cluster M) are increased (6.1% in Dnmt3a KO versus 1.8% in Dnmt3a WT ), as indicated by the bar graphs of keratinocyte lineage subclusters. Subclusters of the hair follicle bulge (A and B), upper hair follicle (C and D), basal interfollicular epidermis (E–G), suprabasal interfollicular epidermis (H–K), and sebocytes (L) are relatively unchanged. ( D ) Canonical gene expression changes were noted in single-cell RNA sequencing comparing keratinocyte lineage cells in Dnmt3a WT versus Dnmt3a KO epidermis ( n = 2 mice per genotype). ( E ) There are 174 differentially expressed genes between Dnmt3a WT and Dnmt3a KO keratinocyte lineage cells with 137 up-regulated in the Dnmt3a KO and 37 down-regulated genes with fold change >2 and FDR <0.01. Unique populations segregated into distinct clusters with graph-based clustering and are represented by uniform manifold approximation and projection (UMAP) dimension reduction plots shown in Fig. 2 A ( 32 ). In total, 19 clusters were identified by unbiased methods (Louvain method graph-based clustering, Dataset S3 ), which were subsequently categorized into 10 populations based on well-established, experimentally validated markers representing known populations in the epidermis ( Fig. 2 A , Right ). The hair follicle bulge (HFB) contains follicular stem cells and is marked by Cd34 expression ( 33 ), while the upper hair follicle (uHF) is marked by beta defensin ( Defb6 ) expression ( 34 , 35 ). Subdivided into the isthmus and infundibulum, the upper hair follicle bridges the follicle with the interfollicular epidermis and sebaceous gland, which is marked by the expression of Mgst1 ( 34 ). The basal interfollicular epithelium (bIFE) contains slow cycling, proliferative cells that act as an epidermal stem cell compartment ( 36 ), and is marked by “basal” keratin gene expression ( Krt5 and Krt14 ), as well as Thbs1 ( 37 ). The suprabasal interfollicular epithelium (sbIFE) is postmitotic and marked by the expression of the genes that encode keratins Krt1 and Krt10 and an early differentiation marker metallothionein Mt4 ( 34 ). Melanocytes were marked by expression of the genes encoding tyrosinase ( Tyr ) and dopachrome tautomerase ( Dct ). Inflammatory cells, including T cells and Cd207-expressing Langerhans cells, comprised 1.3% and 1.1% of total cells in the samples, respectively. Finally, one population was dominated by expression of several genes associated with cellular proliferation (“Pro”), including genes encoding cyclins, cyclin-dependent kinases, and later cell cycle components like Ube2c , which encodes part of the anaphase promoting complex (see Fig. 4 below); this suggests that cells committed to mitosis have a signature dominated by components necessary for execution of the cell cycle ( 38 ). The Krt14 -Cre conditional knockout of Dnmt3a targets the epidermal keratinocyte lineages; we therefore excluded skin-infiltrating hematopoietic cells from subsequent analyses. Keratinocyte lineage UMAP representations of single-cell data demonstrate a high degree of similarity among samples, and unbiased clustering demonstrated similar distributions between populations ( Fig. 2 B , clusters A–M). Small, but statistically significant differences in population sizes between Dnmt3a WT and Dnmt3a KO epidermal keratinocytes were noted in the follicular bulge cluster 2 and upper hair follicle 1 (clusters B and C), as well as basal and suprabasal interfollicular epidermis (clusters E and I); however, the most dramatic difference was in the proliferative (Pro) cluster (cluster M), which is expanded 3.4-fold in Dnmt3a KO skin (1.8% vs. 6.1%, P = 2.2*10 −16 , Fig. 2 C ). In addition to this population shift, there were canonical gene expression changes associated with Dnmt3a deficiency ( Fig. 2 D ). In total, there were 174 differentially expressed genes (DEGs) exhibiting a fold change of 2 or greater with a FDR of 0.01 or less. Most of these genes (137) were more highly expressed in the Dnmt3a KO cells, while 37 were expressed at lower levels ( Fig. 2 D and E ). Of the 23,983 total DMRs identified by WGBS, 16,772 were within 10 kb of a gene, representing 5,573 unique genes. Of these, 2,181 genes (39%) were differentially expressed, as defined by an FDR of 1, while 47 of these genes had an expression change of greater than twofold (38 up-regulated, and 9 down-regulated in Dnmt3a KO cells, Fig. 3 A ); however, the magnitude of methylation change was not predictive of expression changes. Further, some genes (e.g., Krt7 and Meis2 ) exhibited dramatic but opposing expression changes, despite associated DMRs that were hypomethylated to a similar extent ( Fig. 3 B ). The location of DMRs within annotated regions (e.g., enhancers, CpG islands, etc.) did not predict the magnitude or direction of expression changes ( SI Appendix , Fig. S2 ). Differentially methylated regions are associated with differentially expressed genes. ( A ) Heatmap from single-cell RNA sequencing represents the 47 differentially expressed genes between Dnmt3a WT and Dnmt3a KO keratinocyte lineage cells that are associated with a Dnmt3a DMR and exhibit an expression fold change >2 with an FDR <0.01. Grayscale bars represent cells from independent samples, two from each genotype. ( B ) Relationship between gene expression and methylation changes between Dnmt3a KO and Dnmt3a WT keratinocytes. There are 5,968 differentially expressed genes with an FDR 2 denoted as red dots. Gray dots represent fold change between 1 and 2. ( C ) IGV views of differentially methylated regions associated with Meis2, which has annotated DMRs in CpG islands, shelves, shores, and gene body. Gene expression of Meis2 changes by −5.77-fold (FDR 5.7*10 −249 ). Dnmt3a KO expresses Meis2 in 172/15,403 (1.11%) of cells compared to Dnmt3a WT 1,106/10,482 (10.6%). ( D ) Differentially methylated regions associated with Hoxc8 , which has DMRs in shores and gene body. Hoxc8 expression changes by −6.72-fold (FDR 2.51*10 −300 ) and is expressed in 234/15,403 (1.52%) Dnmt3a KO cells compared to 1,375/10,482 (13.1%) Dnmt3a WT . ( E ) Differentially methylated regions associated with Krt7, which has DMRs in shores, shelves, promoter, and gene body. Expression is changed by 7.19-fold in Dnmt3a KO (FDR 9.2*10 −182 ), and cells expressing increase to Dnmt3a KO 2,551/15,403 (16.6%) from Dnmt3a WT 573/10,482 (5.47%). ( F ) Differentially methylated regions associated with Irx1 , which has DMRs in CpG islands, shores, promoter, and gene body. Expression is changed by 2.02-fold (FDR 1.16*10 −109 ) in Dnmt3a KO keratinocyte lineage cells, and cells expressing increases from 2,873/15,403 (18.7%) Dnmt3a KO cells to 885/10,482 (8.44%) Dnmt3a WT . ( G ) Simple Western demonstrates protein level expression changes in Irx1 and Krt7 from epidermal cells isolated from the indicated Dnmt3a WT and Dnmt3a KO littermates (bracketed). Actin serves as a loading control. Meis2 is a homeobox transcription factor involved in developmental processes and regulation of tumor proliferation ( 39 ), and Hoxc8 modulates Wnt signaling in the hair follicle niche to control hair follicle regeneration ( 40 ). Both are DMR-associated genes and are expressed at significantly lower levels in most Dnmt3a KO epidermal cells ( Fig. 3 C and D and SI Appendix , Fig. S3 A ), suggesting that the normal expression of these two genes may be linked to DNA methylation marks added by Dnmt3a. In contrast, some genes are differentially expressed only in subpopulations of skin cells. For example, Keratin 7 ( Krt7 ), which is frequently overexpressed in aggressive subtypes of cutaneous squamous cell carcinoma, and associated with poor survival in patients with esophageal squamous cell carcinoma ( 41 , 42 ), is predominantly expressed in the upper hair follicle in Dnmt3a WT mice; in Dnmt3a KO mice, its expression increases in the upper hair follicle as well as the basal and suprabasal interfollicular epidermis compartments ( Fig. 3 E and SI Appendix , Fig. S3 B ). Similarly, Irx1 is primarily expressed in the hair follicle bulge and upper hair follicle cells of Dnmt3a WT mice, but its expression increases in all compartments in Dnmt3a KO mice ( Fig. 3 F and SI Appendix , Fig. S3 B ), a finding that was validated by Western blotting ( Fig. 3 G ). Decreased methylation in the IRX1 promoter is associated with increased expression and metastatic capability in osteosarcoma ( 43 ). The most significant population shift between Dnmt3a WT and Dnmt3a KO keratinocyte lineage cells occurred in a cluster that was dominated by the expression of genes associated with cell division ( Figs. 2 B and C and 4 A ). This population (cluster M) was defined by the expression of late interphase cyclins (e.g., Ccna2 and Ccnb2 ), as well as later cell cycle genes like Cenpf ; these genes are rarely expressed in the cells of any other clusters ( Fig. 4 B ). The size of this “proliferative cluster” expanded from 1.8% of total epidermal cells in Dnmt3a WT mice to 6.1% in Dnmt3a KO mice (3.4-fold, Fig. 4 A ), similar to the increase in the fraction of cells expressing Mki67 (from 2.8 to 10.2%, 3.6-fold, Fig. 4 C ). The expression of cell cycle genes in cluster M was equally high in Dnmt3a WT and Dnmt3a KO cells, as expected, since all cells in this cluster are likely to be actively dividing, regardless of their genotype ( Fig. 4 D ); however, Dnmt3a KO cells in clusters A–L demonstrated increased expression levels of cell cycle genes, including Mki67 (although the expression levels of these genes were far lower than that of proliferating cells), as well as the essential cyclin-dependent kinase Cdk1 , which is a differentially expressed gene associated with a DMR ( Fig. 3 A ). To validate the proliferative phenotype with an orthogonal method, we injected 5-bromo-2′-deoxyuridine (BrdU) intraperitoneally into young (8 to 9 wk old) Dnmt3a WT and Dnmt3a KO mice, and generated single-cell suspensions of epidermal cells 24 h later for flow cytometric analysis of BrdU uptake. As predicted by the single-cell RNA-sequencing (scRNA-seq) data, the fraction of actively cycling, Cd45-negative epidermal cells was significantly increased in the Dnmt3a KO mice ( Fig. 4 E ). Epidermal subpopulations likewise demonstrated significantly increased BrdU incorporation in Dnmt3a KO basal interfollicular epidermis cells (and trended upwards in the upper hair follicle); however, no change was detected in the follicular bulge stem cell compartment, or the postmitotic suprabasal interfollicular epidermis, neither of which proliferates in the telogen state ( SI Appendix , Fig. S4 A and B ). Dnmt3a-deficient keratinocytes exhibit a proliferative phenotype. ( A ) UMAP of Dnmt3a WT and Dnmt3a KO single-cell RNA sequencing data representing a pair of samples for each indicated genotype, highlighting the population (in red) identified by unbiased clustering as dominated by expression of genes involved in cell cycling (cluster M). ( B ) Proliferative gene expression in the proliferative cluster (cluster M) versus the aggregate of the remaining keratinocyte linage clusters (clusters A–L). These data represent the aggregate of Dnmt3a WT and Dnmt3a KO cells. Each dot is an individual cell. Expression values are represented by least-squares (LS) means. ( C ) The number of cells expressing the proliferative marker gene Mki67 increases from 298/10,485 (2.8%) cells in Dnmt3a WT to 1,509/14,815 (10.2%) cells in the Dnmt3a KO . ( D ) Single-cell RNA-sequencing expression data of multiple cell cycle relevant genes, comparing the proliferative cells (cluster M) versus the rest (cluster A–L). *FDR <0.01 and ns, not significant. Expression values are represented by LS means. ( E ) BrdU uptake in vivo in keratinocytes assessed by flow cytometry (Cd45-negative, single-cell epidermal suspension) comparing Dnmt3a WT ( n = 7) and Dnmt3a KO ( n = 5) mice injected 24 h earlier with BrdU. BrdU uptake is significantly higher in Dnmt3a KO epidermal keratinocytes (5.6% BrdU+ in Dnmt3a WT vs. 7.5% BrdU+ in Dnmt3a KO , P = 0.023). Patients with TBRS have germline mutations in DNMT3A and exhibit several common clinical features, including overgrowth, characteristic facies, and neuropsychiatric disorders ( 44 , 45 ), as well as an increased risk of developing AML and perhaps other malignancies. We previously demonstrated that the peripheral blood cells from a nonleukemic patient with a heterozygous germline DNMT3A R882H allele exhibited a focal, canonical hypomethylation phenotype similar to that of AML samples initiated by the same mutation, suggesting that hypomethylation precedes transformation ( 30 ). To determine whether this mutation causes a methylation phenotype in human skin, we analyzed DNA methylation in two skin samples from the same TBRS patient (at ages 9 and 14) and his unaffected brother (at ages 13 and 17). Global DNA methylation levels were minimally altered (similar to Dnmt3a KO skin). Nearly all of the 1,234 DMRs identified (1,221/1,234; 98.9%) were hypomethylated ( Fig. 5 B ); the smaller number of DMRs (compared to Dnmt3a KO mouse skin) probably reflects the fact that the R882H mutation is dominant negative (80% reduced activity) rather than null, and also the reduced power to detect DMRs with only two samples from each genotype. Plotting the mean methylation values of each DMR between both samples from the two brothers demonstrated the canonical nature of most of these DMRs ( Fig. 5 C ). These DMRs were most likely to occur in gene bodies, similar to that seen in Dnmt3a KO mice ( Figs. 1 D and 5 E ). A germline, dominant-negative DNMT3A mutation causes focal DNA methylation alterations in human skin. ( A ) Whole genome bisulfite sequencing analysis for whole skin samples comparing two samples from a TBRS patient and his brother, comparing CpG methylation globally ( Top ) or at DMRs ( Bottom ). Level of significance * adjusted P < 0.05 and ** P < 0.01. ( B ) Density plots demonstrate global CpG methylation ( Top ) and CpG methylation at DMRs ( Bottom ) in the TBRS patient and his unaffected sibling. ( C ) Heatmap illustrates 1,234 differentially methylated regions in the skin of the TBRS patient, compared with his sibling. ( D ) IGV visualization of CpG methylation at the IRX3/Irx3 and FZD1/Fzd1 loci. Each row represents an independent sample with Dnmt3a WT vs. Dnmt3a KO represented in blue and red ( Top ) and sibling vs. TBRS patient in purple and yellow ( Bottom ). Gray boxes highlight areas of differential methylation. ( E ) Gene bodies are the most frequently involved annotated region by TBRS DMRs. *There are 963 DMRs that reside in annotated gene bodies, representing 637 unique genes. ( F ) There are 23,983 DMRs in the Dnmt3a KO , 16,315 of them are within 10 kb of 3,612 unique genes. Of the 1,234 TBRS DMRs, 1,086 are located within 10 kb of genes in the human genome. Of the 1,086 TBRS-DMR associated genes, 707 have a mouse homolog, with 432 genes overlapping DMR-associated genes in Dnmt3a KO mice and the TBRS patient. At many loci, there was a striking similarity in DNA methylation patterns between the DMRs of homologous genes in mice and humans. For example, DMRs near the homeobox gene IRX3 in human TBRS skin are similar to DMRs near Irx3 in mouse skin ( Fig. 5 D , Left ). Another gene important for development, FZD1 , shares a similar methylation phenotype between the TBRS skin and mouse Dnmt3a KO skin ( Fig. 5 D , Right ). The Wnt signaling pathway is important both for skin and bone development, and multiple genes of the Wnt pathway also exhibited similar Dnmt3a-dependent methylation patterns between mouse and human skin ( SI Appendix , Fig. S5 ). Of the 1,034 DMR-associated genes identified in the TBRS patient, 432 overlapped with Dnmt3a KO DMR-associated genes identified in mouse skin ( Fig. 5 F ).

The purpose of this study was to define the baseline epigenetic characteristics of Dnmt3a-deficient skin. Using WGBS, we defined the focal, canonical DNA methylation phenotype of Dnmt3a-deficient mouse epidermal cells and of a human patient with a germline DNMT3A R882H mutation. The small changes in the global methylation state of epidermal cells reflects the restricted places in the genome where Dnmt3a acts, and the small sizes of these regions (on average, less than 1 Kb). There were 1,234 DMRs identified in the skin of the TBRS patient, which affected 3.19% of annotated gene bodies. Similarly, Dnmt3a-dependent methylation in murine epidermis occurred most often in gene bodies—more than in promoters or CpG island shores, which were the next most common regions affected. However, the magnitude and location of these differentially methylated regions were not predictive of gene expression changes detected by scRNA-seq. Overall, there were 5,968 DEGs in Dnmt3a KO skin cells, but the number of DEGs that were close to a DMR (i.e., within 10 Kb) represented a subset of these genes (2,181/5,968; 36.5%). Only a small subset of these exhibited an expression fold change of >2 (2.15% of DMR associated DEGs, 0.8% of all DEGs). Like multiple previous studies ( 18 , 30 , 46 ), the relationship of DMRs to gene expression patterns was unpredictable at a global level, again demonstrating that DNA methylation is but one factor in determining the expression of a specific gene in a specific cellular context. Prior studies demonstrated accurate remethylation of Dnmt3a-dependent DMRs in Dnmt3a KO bone marrow with induced expression of wild-type DNMT3A, as well as correction of a subset of gene expression changes ( 18 ). However, it is not yet clear whether remethylation was directly responsible for altering the expression of these differentially expressed genes. Aging and sun exposure have been demonstrated to cause DNA methylation changes, which were particularly apparent comparing UV-exposed skin in older patients to UV-protected skin in younger patients. Large, predominantly hypomethylated blocks were found to cover 99 Mb of the genome in older, UV-exposed skin ( 5 ). Although the cause of this phenotype is unknown, a number of possibilities exist, including a progressive decline in expression (or function) of the DNA methyltransferases, or increased activity of the demethylases (i.e., TET gene family members). In mice, Dnmt1 and Dnmt3a are expressed in the basal layer of keratinocytes during development; however, in adults, the expression of Dnmt1 is restricted to the hair follicle, and Dnmt3a expression is restricted to the interfollicular epidermis ( 8 , 11 ). Notably, Dnmt3a expression decreases in this compartment in adult mice ( 11 ), though it is not yet clear whether these findings contribute to the increased risk of skin cancers in elderly patients. To begin to address this question, Rinaldi et al. demonstrated that epidermal Dnmt3a deficiency resulted in decreased tumor latency and increased tumor numbers in mouse skin treated with DMBA/TPA, suggesting that the threshold for oncogenic transformation is lower in mice with Dnmt3a deficiency ( 11 ). In the same study, mice with Dnmt3b deficiency had no phenotype, and compound deficiency of Dnmt3a and Dnmt3b (double knockout/DKO) resulted in tumor latency that was indistinguishable from that of Dnmt3a deficiency; however, a broader spectrum of tumors arose in the DKO mice, including basal cell carcinomas, as well as more aggressive tumors (including squamous cell and spindle cell carcinomas). These tumors developed a reliance on Pparg ; we did not detect up-regulation of Pparg in the premalignant state, suggesting that dysregulation may occur during or after transformation. In contrast, certain gene expression changes (i.e., increased expression of G0s2 ) persisted from the Dnmt3a-deficient preneoplastic state to the resulting tumors ( 11 ) ( Fig. 1 E and Dataset S2 ). Together, these data strongly suggested that Dnmt3a deficiency creates a premalignant state for keratinocyte carcinomas in mice, with some gene expression changes that persist in tumors. We observed that Dnmt3b deficiency causes only 20 DMRs in the skin, and that Dnmt3a and Dnmt3b compound deficiency (DKO) adds little to the Dnmt3a KO methylation phenotype ( SI Appendix , Fig. S6 ). Therefore, in mouse skin, Dnmt3a is the dominant de novo DNA methyltransferase, and the loss of its function in aging skin could likewise be relevant for human skin cancer pathogenesis. Although Dnmt3a deficiency results in minimal skin population changes, it is associated with an expanded number of cells with a proliferative signature in the Dnmt3a KO epidermis. Despite this, the skin of young mice is not measurably thickened, and spontaneous skin cancers do not occur, suggesting that robust compensatory mechanisms are active. Multiple critical cell cycle genes, including Ccna2 or Cdk1 , are up-regulated in both proliferating and nonproliferating cells as well; conversely, the cyclin-dependent kinase inhibitor Cdkn1a is down-regulated in nonproliferating Dnmt3a-deficient cells ( SI Appendix , Fig. S7 A ). The expression of Mki67 is tightly regulated and absent in quiescent (G0) cells and more highly expressed in G2 through M phase cells. Cells that are more frequently in the cell cycle express higher levels of Mki67 than quiescent cells ( 47 ). Indeed, we found that Mki67 expression is elevated even in the nonproliferating cells of Dnmt3a KO mice, which may reflect the proliferative history of these cells. Together, these results suggest that Dnmt3a deficiency may prime epidermal keratinocytes to enter the cell cycle, or remove antiproliferative “brakes” from mitogenic signals. While proliferative signaling was implicated, Cdk1 was the only direct proliferative driver associated with a DMR and expressed with a fold change >2. Enhanced proliferation has been observed in the setting of DNMT3A deficiency in other premalignant and malignant conditions. Decreased DNMT3A expression has been associated with nonmalignant adenomyosis, where ectopic endometrial cells invade the myometrium. These cells exhibit increased proliferation with DNMT3A knockdown, and decreased proliferation with overexpression ( 48 ). Additionally, a model of Dnmt3a deficiency in KRAS G12D driven lung tumors demonstrated that at 24 wk after induction, Dnmt3a KO tumors were six times larger than their Dnmt3a WT counterparts and exhibited a markedly higher proliferative index ( 49 ). Together, these findings suggest that epigenetic changes caused by Dnmt3a deficiency can result in a hyperproliferative state in several cellular contexts. The specific genes responsible for initiating this phenotype are not yet clear, but several good candidates (i.e., Cdk1 , Hras , E2f1 , and Cdkn1a ) are defined in this study ( SI Appendix , Fig. S7 B – D ). Proliferative priming may therefore represent one mechanism that is relevant for the premalignant state associated with DNMT3A haploinsufficiency or deficiency in multiple tissues. The resulting expanded proliferative capacity may exist in equilibrium until cooperating mutations occur, which appear to be essential for transformation in many model systems. In addition to clonal hematopoiesis, premalignant mutation-driven clonal expansion has now been shown to occur in many other systems, including the skin ( 2 ), liver ( 50 ), esophagus ( 51 , 52 ), and bronchial epithelium ( 53 ). Mutationally driven, expanded clones usually do not progress to frank malignancies, suggesting that additional epigenetic or genetic alterations are required to cause transformation. In clonal hematopoiesis, the expansion of hematopoietic stem/progenitor cells is usually caused by mutations in epigenetic regulatory genes, most commonly in DNMT3A ( 23 , 24 , 26 , 54 ). These clonally expanded cells are at increased risk of transformation to myelodysplastic syndromes or AML, likely by altering the epigenetic “fitness” of these cells for transformation ( 55 – 57 ). Similar epigenetic fitness alterations may modify tumor susceptibility in skin and other organs by a variety of mechanisms. We identified similarities between the methylation phenotypes of mice and humans with decreased DNMT3A activity: DMRs in the epidermis of Dnmt3a KO mice shared considerable gene-level similarity with a TBRS patient with a germline DNMT3A R882H mutation. In fact, 432 out of the 1,086 gene-associated DMRs identified in the skin of this patient overlapped with gene-associated DMRs found in Dnmt3a KO skin. Several of these overlapping genes have known roles in development, including IRX3/Irx3 ; this locus shares a high degree of evolutionary conservation in noncoding sequences between human, mouse, and zebrafish, suggesting the presence and conservation of both functional and regulatory elements ( 58 ). Additionally, important pathways in skin development (like the WNT pathway) have multiple members with similar DNMT3A/Dnmt3a -dependent methylation patterns ( SI Appendix , Fig. S5 ). Recent descriptions of AMLs and central nervous system (CNS) tumors in young patients with germline loss-of-function DNMT3A mutations strongly suggest that DNMT3A acts as a tumor suppressor in both hematopoietic and CNS cells ( 59 ). Although the skin phenotype and skin cancer susceptibility of children with TBRS has not yet been described, our data suggest that these patients may be at elevated risk, and that appropriate primary prevention measures should be taken to protect these patients from sustained UV light exposure. Finally, these data also suggest a role for DNMT3A in keratinocyte carcinomas arising in aging and sun-exposed skin; additional studies to define this role are in progress.

Samples were obtained as part of a study approved by the Human Research Protection Office at Washington University School of Medicine (approval #201011766), which explicitly allows for potentially identifying genomic studies, including whole genome sequencing. Written informed consent was obtained from the parents of the patient with TBRS (and his unaffected brother) for whole genome bisulfite sequencing studies, in accordance with the institutional review board approved protocol at Washington University School of Medicine. All mouse work was performed in accordance with institutional guidelines and approved by the Animal Studies Committee at Washington University. For all studies, we harvested single-cell suspensions of dorsal epidermis in second telogen (7 to 10 wk of age), or in older mice, when skin is >90% telogen based on pale or pink skin color. Details for mouse strain and epidermal single-cell suspensions can be found in SI Appendix , Supplementary Materials and Methods . WGBS was performed as described previously ( 17 ) with the exception that DNA was generated from epidermal cells. Alignment and DMR calling details are in SI Appendix , Supplementary Materials and Methods . Single-cell suspensions of murine epidermis were freshly prepared and libraries generated using the 10× Chromium Controller system. Data analysis was performed using the standard tools available in the Partek Flow software suite (Partek, Inc.). Full details of the analysis are in SI Appendix , Supplementary Materials and Methods . Freshly isolated single-cell suspensions of murine epidermis were lysed in radioimmunoprecipitation assay (RIPA) buffer or suspended in staining buffer for further analysis by flow cytometry. Lysis and staining conditions (as well as antibody information) are described in SI Appendix , Supplementary Materials and Methods .