Mutant Lef1 controls Gata6 in sebaceous gland development and cancer

Abstract Mutations in Lef1 occur in human and mouse sebaceous gland (SG) tumors, but their contribution to carcinogenesis remains unclear. Since Gata6 controls lineage identity in SG, we investigated the link between these two transcription factors. Here, we show that Gata6 is a β‐catenin‐independent transcriptional target of mutant Lef1. During epidermal development, Gata6 is expressed in a subset of Sox9‐positive Lef1‐negative hair follicle progenitors that give rise to the upper SG. Overexpression of Gata6 by in utero lentiviral injection is sufficient to induce ectopic sebaceous gland elements. In mice overexpressing mutant Lef1, Gata6 ablation increases the total number of skin tumors yet decreases the proportion of SG tumors. The increased tumor burden correlates with impaired DNA mismatch repair and decreased expression of Mlh1 and Msh2 genes, defects frequently observed in human sebaceous neoplasia. Gata6 specifically marks human SG tumors and also defines tumors with elements of sebaceous differentiation, including a subset of basal cell carcinomas. Our findings reveal that Gata6 controls sebaceous gland development and cancer.


Introduction
At the core of the canonical Wnt signaling pathway is nuclear translocation of b-catenin and consequent transcription of downstream target genes through b-catenin binding to members of the lymphoid enhancer-binding factor/T-cell factor (Lef/Tcf) transcription factor family (Klaus & Birchmeier, 2008;Nusse & Clevers, 2017). The Wnt pathway plays a central role in stem cell maintenance and fate specification in mammalian epidermis (Watt & Collins, 2008;Lim & Nusse, 2012), controlling the balance between hair follicle (HF) and sebaceous gland (SG) differentiation. During embryonic life and during the postnatal hair cycles, activation of bcatenin triggers HF growth (Huelsken et al, 2001;Lowry et al, 2005;Donati et al, 2014). Ectopic HF can also be generated upon transient activation of epidermal b-catenin (Lo Celso et al, 2004;Silva-Vargas et al, 2005), in particular in Lrig1-positive and Lgr6-positive stem cell (SC) populations of the upper pilosebaceous unit (Kretzschmar et al, 2016). While Wnt signaling favors HF over SG fate, an N-terminally truncated form of Lef1 (DNLef1), unable to bind b-catenin, converts HF into keratinized epidermal cysts with ectopic sebocytes (Merrill et al, 2001;Niemann et al, 2002;Donati et al, 2017).
Deregulation of the Wnt/b-catenin pathway occurs in several skin cancers. Transgenic mice overexpressing a stabilized form of bcatenin develop HF tumors: pilomatricomas or trichofolliculomas (Gat et al, 1998;Lo Celso et al, 2004). In humans, stabilizing mutations in b-catenin are found in a majority of pilomatricomas (Chan et al, 1999) and pilomatrix carcinomas (Lazar et al, 2005).
While genetic deletion of b-catenin from the epidermis is not associated with tumor development (Huelsken et al, 2001;Malanchi et al, 2008), transgenic mice expressing DNLef1 under the control of the keratin 14 promoter (K14DNLef1) spontaneously develop skin tumors, most of which are sebaceous adenomas and sebaceomas (Niemann et al, 2002). In K14DNLef1 mice, DNLef1 decreases endogenous Lef1 expression and acts as a dominant negative inhibitor of b-catenin (Niemann et al, 2002). Therefore, Wnt/b-catenin targets are downregulated  and the hair follicle cycle is compromised (Niemann et al, 2002). The DNLef1 transgene acts as a tumor promoter in chemical carcinogenesis experiments. Tumors that develop upon DNLef1 expression exhibit sebaceous differentiation rather than the papillomas and squamous cell carcinomas characteristic of wild-type (WT) mice (Niemann et al, 2007). Consistent with these findings, mutations in the N-terminus of Lef1 that prevent b-catenin binding are found in approximately 30% of human benign sebaceous tumors (Takeda et al, 2006) and 20% of eyelid sebaceous carcinomas (Jayaraj et al, 2015).
We recently reported that the transcription factor Gata6 plays a role in sebaceous lineage determination and is highly upregulated in the junctional zone (JZ) of K14DNLef1 mice . Therefore, in the present work, we sought to explore the role of Gata6 in establishing the sebaceous lineage and to understand the link between Gata6 and Lef1 in the context of sebaceous gland tumors. Our findings demonstrate that sebaceous fate during cancer is controlled by Lef1 and Gata6.

Results
Gata6 is a DNLef1 target gene As reported previously (Merrill et al, 2001;Niemann et al, 2002;Petersson et al, 2011Petersson et al, , 2015, in adult K14DNLef1 mice, the HF undergoes conversion to multilayered, keratinized cysts with ectopic sebocytes (Fig 1A and B). Based on our earlier studies, we now believe that these cysts represent abnormal Gata6-positive sebaceous ducts (SD) while the associated mature sebocytes do not express Gata6 . To identify direct DNLef1 target genes, we performed chromatin immunoprecipitation and nextgeneration sequencing (ChIP-Seq) on keratinocytes isolated from K14DNLef1 epidermis ( Fig 1C). As expected, the intact DNLef1 DNA binding domain bound the same DNA core motif as Lef1, as, for example, on Vgll4 and Runx1 loci ( Fig 1C).
To understand the transcriptional role of DNLef1 in the epidermis, we compared DNLef1 direct target genes to the gene signatures of the three major epidermal compartments (interfollicular epidermis (IFE), HF, and SG) obtained upon micro-dissection of adult tail skin from WT mice (Fig 1D and E;Donati et al, 2017). We observed that DNLef1 mainly bound to repressed genes belonging to the SG signature ( Fig 1D and E), consistent with a known TLE/Grouchodependent repressive function of Lef1 in the absence of b-catenin (Niemann et al, 2002;Ramakrishnan et al, 2018).
Further analysis of the ChIP-Seq data revealed that Gata6 was a direct transcriptional target of DNLef1 (Fig 1C). We confirmed by RT-qPCR that Gata6 expression was increased in K14DNLef1 keratinocytes as compared to WT cells (Fig 1F). RNAi-mediated knockdown of Lef1 led to a striking Gata6 downregulation (Fig 1F), while DNLef1 overexpression in the SebE6E7 sebocyte cell line increased Gata6 expression (Fig 1G). These results are in agreement with the expansion of Gata6-positive cells observed in K14DNLef1 mice (Fig 1H;Donati et al, 2017).
To uncover the in vivo changes in gene expression linked to DNLef1 expression, we compared the gene expression profiles of flow-sorted total basal keratinocytes (basal, Itga6 + Cd34 À ) and bulge stem cells (HFSC, Itga6 + Cd34 + ) in WT and K14DNLef1 transgenic mice ( Fig EV1A). To detect early molecular events, we collected cells from 9.5-week-old mice, when the HFs are in the resting (telogen) phase of the hair growth cycle (Oh et al, 2016), and the DNLef1 phenotype is not yet fully apparent (Niemann et al, 2002). Differentially expressed genes (DEG) between the different keratinocyte populations were validated by RT-qPCR ( Fig EV1B). DNLef1 expression favored gene expression characteristic of the pilosebaceous unit while repressing the IFE gene signature (Fig 1I), and differently impacted several signaling pathways in the three compartments ( Fig EV1C). The correlation between bulge and SG signatures in K14DNLef1 mice significantly increased when we selected genes that were direct targets of DNLef1 (white box in Fig 1J), suggesting a direct role of DNLef1 in the formation of ectopic SG in the HF. In addition, Gene Ontology (GO) of DEG in WT as compared to K14DNLef1 cells showed enrichment for the terms "tissue morphogenesis" and "gland development" (Fig 1K).
A Hematoxylin and eosin (H&E) staining of WT and K14DNLef1 HF from back skin at 11.5 weeks (early anagen). White arrowheads: ectopic SG and epidermal cysts. B Tail epidermal whole mounts from WT and K14DNLef1 mice labeled with anti-Fabp5, Krt15, and counterstained with Dapi. C Average signal intensity of DNLef1 binding sites detected by a Lef1 antibody (red) as compared to input alone (gray) and negative control (Neg. Ctl., black). DNA motif analysis revealed Lef1 classical consensus sequence. Representative plots of ChIP-Seq reads aligned to the Vgll4, Runx1, and Gata6 loci. D Heat-map (Pearson's correlation) of differentially expressed genes (DEG) between different micro-dissected regions (IFE, HF, and SG) from WT mice (left panel) . DNLef1 Chip-Seq data show DNLef1 direct transcriptional targets in black (right panel). E Gene Set Enrichment Analysis comparing gene location of DNLef1 peaks and SG gene expression signature. F RT-qPCR analysis of Gata6 in WT and K14DNLef1 primary keratinocytes transfected with siRNA targeting DNLef1 or a scrambled sequence (scr.). Data are means AE SEM of three independent wells. G RT-qPCR analysis of Gata6 in human SebE6E7 sebocytes transfected for 48 h with a mock plasmid or a Lef1-or DN34Lef1-expressing plasmid. DN34Lef1 is the human ortholog of murine DN32Lef1, which is expressed in K14DNLef1 transgenic mice (Takeda et al, 2006). Data are means AE SEM of three independent experiments. H Dorsal skin sections of WT and K14DNLef1 mice stained for Krt14 and Gata6. I Heat-maps of DEG from  or sorted (Itga6 + Cd34 À , basal; Itga6 + Cd34 + , HFSC) keratinocytes in comparison with IFE, SG, and HF gene signatures ranked from high to low expression. J Heat-map depicting similarities, as Pearson's correlation coefficient, between DEG in keratinocytes from WT or K14DNLef1 mice and IFE vs SG vs HF transcriptome analysis (left panel). Additional correlation with DNLef1 Chip-Seq data is depicted in the right hand heat-map. K GO enrichment analysis of DNLef1 direct target genes among DEG of WT and K14DNLef1 mice.
Data information: (A, B, and H) Scale bars: 50 lm. (F, G) Statistical analyses were performed with an ordinary one-way ANOVA: (ns) not significant; *P < 0.05; ***P < 0.0005.  regulators of the SG (Chang et al, 2009;Niemann & Horsley, 2012). In parallel, DNLef1 directly downregulates genes such as Klf5, important for IFE identity (Ge et al, 2017), Igfbp3 that is expressed in SG (Dahlhoff et al, 2016), and Slco2a1. Loss-of-function mutations in Slco2a1 are associated with sebaceous hyperplasia (Guo et al, 2017). We conclude that DNLef1 has activating and repressive functions that converge in promoting a SD/SG phenotype in K14DNLef1 mice.

WT
Gata6 expression is not dependent on Lef1 or b-catenin during epidermal development Although Gata6 is a DNLef1 transcriptional target, DNLef1 did not trigger ectopic Gata6 expression in the IFE of K14DNLef1 mice (Fig 1H), suggesting that the upper pilosebaceous compartment is the only permissive niche for DNLef1-mediated Gata6 expression. During early HF/SG morphogenesis (from E15.5 to E18.5), Gata6 and Lef1 were not co-expressed in epidermal cells (Fig 2A), and at E18.5, Lef1 À/À epidermal Gata6 expression was indistinguishable from control epidermis ( Fig 2B). In embryonic WT skin, endogenous Gata6 appeared mainly in suprabasal cells located in the upper part of stage 4 HF (late hair peg; Fig 2A). While Gata6 and Lef1 did not co-localize in developing HF, Gata6 expression was detected in a subset of Sox9-positive cells in stage 4 and 5 HF (Fig 2A). In anagen HF of adult WT mice, endogenous Lef1 and Gata6 did not co-localize ( Fig EV2A). From birth, Gata6 was expressed, as in the adult epidermis, in the JZ and upper SG, but not in mature sebocytes ( Fig EV2B). Consistent with the restriction of Sox9 expression to bulge cells in adult HF (Nowak et al, 2008), there was little co-expression of Gata6 and Sox9 at P1 (Fig EV2B). Gata6-positive cells did not colocalize with another bulge marker, CD34. Cells expressing low levels of Gata6 co-expressed Lgr6 and Lrig1 in neonatal P1 mice, as in adult animals . No co-expression of Gata6 and Tcf3/4 was observed. Gata6 did not co-localize with Ki67 ( Fig EV2B), indicating that Gata6-positive cells were not proliferative, in agreement with the role of Gata6 in terminal differentiation of the SD lineage . During HF morphogenesis in human embryonic skin, Gata6 expression was initiated at a similar HF stage to the mouse and was located in the JZ, SD, and upper SG of more mature HF (Fig EV2C).
To test whether Gata6 expression was regulated by canonical Wnt signaling, we used inducible and constitutive epidermal bcatenin gain-of-function mouse models (Fig 2C and D). K14DNb-CateninER mice in which N-terminally truncated b-catenin is fused with the estrogen receptor ligand binding domain (ER) (line D2, 12 copies; line D4, 21 copies) served as the inducible model (Lo Celso et al, 2004). Upon application of tamoxifen (4OHT) for 1 or 6 days, HF entered anagen in both mouse lines. In the D4 line, this was accompanied by massive thickening of the existing HF and induction of ectopic follicles. However, we did not observe any ectopic expression of Gata6 ( Fig 2C).
For constitutive b-catenin activation, we generated K14Cre/bCat Flox(ex3)/+ mice in which stabilized b-catenin lacking the GSK3b phosphorylation sites becomes prematurely active in all basal keratinocytes during embryogenesis. In these mice, the IFE and SG adopt a HF fate (hair shaft; Kretzschmar et al, 2016). We observed induced ectopic Sox9 expression in K14Cre/bCat Flox(ex3)/+ epidermis. Nevertheless, premature and broad epidermal b-catenin activation, even in this developmental context, did not trigger Gata6 expression ( Fig 2D).
We conclude that although Gata6 is a mutant Lef1 target gene, Gata6 expression is neither Lef1 nor b-catenin-dependent during pilosebaceous unit morphogenesis.

The Gata6 lineage participates in sebaceous gland morphogenesis
We next investigated the role of Gata6 in SG morphogenesis. We did not observe any major SG abnormalities when Gata6 was deleted via the Krt5 promoter (cKO; Donati et al, 2017). However, an indepth analysis of tail epidermis did indicate that in the absence of Gata6 the proportion of hypotrophic SG was significantly higher than in WT and heterozygous mice, albeit constituting a minority of the total SG ( Fig 3A).
By overexpressing Gata6 in primary murine keratinocytes, we previously demonstrated that Gata6 induces a SD/SG transcriptional program upon differentiation in vitro . In addition to SD genes, Gata6 triggered the expression of genes associated with sebocyte differentiation, such as Pparg and Fasn, and downregulated genes associated with the androgen receptor gene signature that is a distinctive feature of the base of the SG . To overexpress Gata6 in developing epidermis, we performed in utero lentiviral infection (Beronja et al, 2013) with a Gata6-ires-GFP virus or a control virus expressing GFP alone (Figs 3B and C). This protocol is characterized by a low infection efficiency, with most transduced cells carrying only one transgene to avoid uncontrolled overexpression and off-target effects (Beronja et al, 2013). Gata6 overexpression in epidermal cells at E9.5 recapitulated aspects of the K14DNLef1 phenotype, in particular the formation of epithelial cysts with ectopic sebocytes (Fig 3B). Gata6 overexpression also resulted in the ectopic expression of the differentiating

of 19
The EMBO Journal 38: e100526 | 2019 ª 2019 The Authors sebocyte marker Fasn in the SG and in the HF but not in the IFE, suggesting that Gata6 is able to promote SG identity in a cell compartment-dependent manner ( Fig 3B). Gata6-induced ectopic sebocytes did not complete their maturation, as shown by the absence of LipidTOX staining ( Fig 3C). This is consistent with the absence of Gata6 staining in differentiated sebocytes in mouse skin (Fig EV2B;Donati et al, 2017). The origins of the SG during development are still a subject of debate . Primitive sebocytes can be detected in stage 5 HF (bulbous hair peg) as described by Paus et al (1999). While it is generally assumed that the entire gland is derived from a single lineage expressing Sox9 (Nowak et al, 2008;Frances & Niemann, 2012), lineage tracing via retrovirus-mediated LacZ transduction of dermabraded skin has shown that individual SG can be derived from more than one cell population (Ghazizadeh & Taichman, 2001). To examine this without having to damage the skin or rely on candidate SG markers, we generated WT:GFP chimeric mice ( Fig EV3A) as previously described (Arwert et al, 2010). Aggregation chimeras are powerful tools that have been extensively used to infer developmental mechanisms (Tam, 2003). We assessed GFP expression in SG from 5 chimeric mice. In labeled SG, GFP staining was frequently localized exclusively to the upper SG ( Fig EV3A). Therefore, we postulate that the SG is formed from at least two progenitors, one of which exclusively populates the upper SG and ducts.
To test whether Gata6 progeny were involved in morphogenesis of the upper SG, we performed lineage tracing by inducing recombination in Gata6EGFPCreERT2 (Gata6creER) × Rosa26-fl/STOP/ fl-tdTomato (ROSA-dTom) mice at E16.5 and E18.5 and analyzing whole mounts of tail and back skin 2 days and 17 days (P13) postlabeling (Figs 3D and EV3B and C). In line with our Gata6 expression analysis (Fig 2A), only cells in the upper part of stage 4-5 HFs were labeled 2 days after recombination ( Fig EV3B). HF labeling efficiency was higher at E18.5 due to the higher abundance of stage 4-5 HF, but Gata6 progeny were never observed in earlier HF stages (Fig EV3B). At P13, Gata6 progeny were specifically located in the upper SG and JZ in approximately 90% of labeled pilosebaceous units (Fig 3D). Gata6 progeny were only rarely found in the lower SG or in the lower part of the HF (Fig EV3C).
Taken together, these data suggest that the Gata6 lineage is responsible for generating the upper part of the SG, including the SD.

Gata6 specifies number and type of K14DNLef1 tumors
We next investigated the role of Gata6 in epidermal carcinogenesis. Unlike WT mice, K14DNLef1 mice develop tumors spontaneously, or after a single DMBA application (Niemann et al, 2007). In both contexts, K14DNLef1 tumors have elements of sebaceous differentiation. Gata6 was expressed in all K14DNLef1 tumors, whereas it was undetectable in DMBA/TPA-induced papillomas in WT mice ( Fig 4A). Ablation of Gata6 in K14DNLef1 mice (K14DNLef1:cKO) resulted in an increased rate of tumor formation following DMBA treatment, and in an increased number of tumors per mouse (Fig 4B), indicating that Gata6 acted as a tumor suppressor. In addition to increasing the overall number of tumors per mouse, loss of Gata6 led to an increase in papilloma-like tumors and a decrease in tumors with SG elements (Fig 4C). K14DNLef1 tumors strongly expressed the SD markers Plet1 and ATP6v1c2 , while WT and K14DNLef1:cKO tumors did not, indicating that SG differentiation was significantly reduced upon Gata6 loss ( Fig 4D  and Appendix Fig S1).

Gata6 acts as a tumor suppressor by controlling the DNA mismatch repair response
Muir-Torre syndrome is a rare genetic condition that predisposes patients to developing sebaceous skin tumors and visceral malignancies. This autosomal dominant variant of Lynch syndrome is caused by mutations in DNA mismatch repair (MMR) genes, resulting in microsatellite instability (Eisen & Michael, 2009a,b;John & Schwartz, 2016). Since loss of Gata6 in cultured mouse keratinocytes leads to DNA damage, triggering apoptosis (Wang et al, 2017), we investigated whether the tumor suppressive function of Gata6 might be due to an effect on MMR gene transcription.
Computational analysis of RNA-Seq data from Gata6-deficient keratinocytes (Wang et al, 2017) revealed that DNA metabolism ◀ Figure 3. Gata6 progenitors contribute to SG morphogenesis.
A Bright-field images of whole-mount WT and Gata6 cKO epidermis, showing abnormal SG upon loss of Gata6. Quantification of hypotrophic SG in WT, heterozygous (Het), and Gata6 cKO epidermal whole mounts. Data are boxplots with indication of means. Box limits are minimum and maximum values. An average of 165 HF per mouse (from 3 to 4 mice per genotype) was analyzed. *P < 0.05; **P < 0.005; unpaired Student's t-test. B Schematic representation of plasmid constructs used for in utero lentiviral infection. Whole mounts or sections of adult infected epidermis with empty vector (EV) or Gata6-ires-GFP (G6OE) lentivirus were stained for GFP and Fasn. In vivo overexpression of Gata6 leads to ectopic Fasn expression in the HF/SG unit. White dotted lines define SG, and yellow dotted lines define a cyst. Note that the cyst is mostly negative for Fasn in agreement with its SD-like phenotype. White arrows indicate GFP-positive infected cells. These cells are stained with Fasn only on G6OE expression. H&E-stained skin from G6OE mice shows a cyst with SG elements in the HF unit (black arrows). Staining for Gata6 (both endogenous and exogenous) shows that Gata6  genes (including a number of genes related to DNA repair and replication) were over-represented in downregulated genes (Fig 5A and B). When we intersected Gata6 ChIP-Seq data      in the genes downregulated on loss of Gata6 (Fig 5A and B). The specific MMR genes Msh3, Exo1, Pold1, Rfc3, and Pcna were downregulated in Gata6 cKO keratinocytes (Fig 5A), although they were not identified as direct transcriptional targets of Gata6 . Using ChIP-Atlas, an extensive database of publicly available ChIP-Seq experiments (Oki et al, 2018 and http://chip-atlas.org), we analyzed Gata6 ChIP-Seq experiments performed in human ES cell-derived mesendodermal cells and in human colon, gastric and pancreatic adenocarcinoma cell lines. We identified eight genes in the MMR pathway for which a Gata6-binding peak was found within 10 kb from the transcription start in at least 2 different datasets: Msh3, Rfc2, Pold4, Ssbp1, Pold3, Exo1, Rfc1, and Mlh3 ( Fig EV4). It is possible that these were not found in Gata6-overexpressing primary mouse keratinocytes  because the ChIP-Seq analyses were performed too soon after lentiviral infection for Gata6 to induce MMR genes. Nevertheless, the MMR genes identified as Gata6 direct transcriptional targets using ChIP-Atlas closely match the MMR genes that were downregulated in Gata6 cKO keratinocytes and confirm that the MMR pathway is regulated by Gata6.
To further study the functional impact of Gata6 on the MMR pathway, we used two methods that are key to diagnosing Muir-Torre syndrome (Eisen & Michael, 2009b) and Lynch syndrome (Giardiello et al, 2014): microsatellite instability (MSI) testing and immunohistochemistry for MMR proteins. We first evaluated MSI in skin tumors from K14DNLef1 mice and K14DNLef1:cKO mice. We analyzed the relative allele frequency of five microsatellite regions, as previously described (Woerner et al, 2015;Germano et al, 2017;Keysselt et al, 2017). Three markers (Bat64, AA003063, and L24372) showed a shift in allele distribution, two of which were significant (Fig 5C). This indicated a more unstable microsatellite phenotype in K14DNLef1:cKO mice.
The MMR pathway corrects errors within newly synthesized DNA strands during replication and mainly relies on MutSa and MutLa complexes formed by Msh2/Msh6 and Mlh1/Pms2, respectively (Jiricny, 2006). Mlh1 (Fig 5D and Appendix Fig S2A) and Msh2 (Fig 5E and Appendix Fig S2B) expression were significantly reduced in K14DNLef1:cKO as compared to K14DNLef1 tumors. The broad downregulation of MMR genes, in particular of Msh3, in Gata6-deleted skin is likely to explain the downregulation of Mlh1 and Msh2 (Figs 5D and E, and Appendix Fig S2) through destabilization of MMR protein complexes.
Altogether, these results indicate that the increased incidence of tumors in K14DNLef1:cKO mice could be due to a decrease in expression of Gata6-dependent MMR proteins.

Gata6 expression is a hallmark of human skin tumors with sebaceous differentiation
To examine whether the observations in mice were relevant to human skin tumors, we first confirmed that Gata6 was expressed in the same locations in human as in mouse adult skin. Indeed, Gata6 was expressed in the upper SG, SD, and JZ ( Fig 6A). We then analyzed Gata6 expression in a large panel of benign and malignant human skin samples (N = 73; Fig 6B and Appendix Fig S3). Eighteen out of 19 human sebaceous tumors were Gata6-positive, regardless of their malignancy grade. In addition, the majority of the tumors harboring elements of SG differentiation were positive for Gata6: sebaceous nevus, syringocystadenoma papilliferum (SCAP, commonly associated with sebaceous nevus), folliculosebaceous cystic hamartoma (FSCH), steatocystoma multiplex (genetic ◀ Figure 5. Gata6 controls the DNA mismatch repair response during sebaceous tumorigenesis. A Heat-map highlighting differentially expressed genes (DEG) between WT and Gata6 cKO keratinocytes from Wang et al (2017). Gata6 direct transcriptional targets from Chip-Seq data  are shown in dark gray. B KEGG pathway analysis performed on downregulated genes upon Gata6 loss. C MSI analysis performed on skin tumors from DMBA-treated K14DNLef1 and K14DNLef1:cKO mice (5 mice in each group). Representative microsatellite profiles of K14DNLef1 and K14DNLef1:cKO tumors. Peak heights were normalized to the highest peak in each microsatellite profile to obtain the relative frequency of each allele for five different markers. "0" indicates the position of the highest peak in bp. *P < 0.05; **P < 0.001; paired t-test. D Representative immunohistochemistry staining for Mlh1 performed on skin tumors from DMBA-treated K14DNLef1 and K14DNLef1:cKO mice. Four mice were included in each group. A technical control is displayed (without primary antibody incubation). Semi-quantitative analysis of Mlh1 staining and statistical analysis (chi-square test) were performed. E Representative immunohistochemistry staining for Msh2 and semi-quantitative analysis performed as in (D).
Data information: (D and E) Scale bar: 250 lm. ▸ Figure 6. Gata6 expression is conserved in human normal skin and found in tumors harboring sebaceous differentiation.
A Normal human scalp skin sections stained with Gata6 and pankeratin antibodies. Endogenous Gata6 expression was found in the JZ, upper SG, and SD. As in mouse skin, Gata6 was not expressed in sweat glands (SwG).         Fig S3). The percentage of Gata6-positive cells was typically in the same range within each tumor group, except in the case of sebaceous carcinomas, which exhibited more variability in Gata6 expression (Fig 6B and Appendix Fig S3). Two BCC subtypes could be distinguished on the basis of Gata6 expression. Gata6 was selectively expressed in BCC with cysts. We speculate that these BCC may arise from a different region of the HF compared to cyst-less BCC. In addition, although the number of samples of the benign skin disorder SCAP was limited, we observed Gata6 expression in almost 70% of cases (Fig 6B and Appendix Fig  S3). This is intriguing because SCAP is thought to derive from sweat glands (Yamamoto et al, 2002), yet normal sweat glands do not express Gata6 (Fig 6A).
In support of our findings, we reanalyzed published RNA-Seq data from human SCC, BCC, and sebaceous carcinomas (SebC; North et al, 2018). We found that Gata6 expression was significantly increased in SebC as compared to BCC and SCC (Fig 6C), which confirms our previous observations (Fig 6B and Appendix Fig S3).
Therefore, our results establish Gata6 as a key histological marker of human skin tumors that originate from, or have differentiated elements of, the SG.

Mutation and downregulation of Gata6 are features of sebaceous carcinomas with a high mutational burden
North et al (2018) have distinguished three subclasses of human SebC based on their mutational profile: the ocular and cutaneous pauci-mutational SebC (with a low prevalence of mutations); SebC with a MSI mutational signature (intermediate prevalence of mutations) and SebC with a UV mutational signature (highest somatic mutation burden). Within this dataset, we found GATA6 mutations in approximately 30% of SebC harboring a MSI or UV damage signature (Fig 7A). GATA6 missense mutations were mostly deleterious (Fig 7B). In addition, analysis of RNA-Seq data showed that Gata6 expression was significantly lower in UV-induced SebC than in other SebC (Fig 7C). In agreement with our observations in mice (Fig 4B  and C), UV-related SebC expressing a low level of Gata6 were less differentiated and more aggressive than the pauci-mutational and MSI-mutant tumors (North et al, 2018).
Despite the different mutational mechanisms associated with MSI-related and UV-induced SebC, Gata6-mutated tumors (Gata6mut) displayed a higher number of mutations/Mb (Fig 7D), of somatic single-nucleotide variants (SSNV) ( Fig 7E) and Indel ( Fig 7F) than Gata6 wild-type tumors (Gata6wt). In addition, Gata6mut tumors displayed a trend of downregulation in MMR genes when compared to Gata6wt tumors ( Fig 7G). The limited number of Gata6mut SebC samples did not allow us to test for significance. However, these results suggest that Gata6 affects DNA damage pathways in human sebaceous tumors as in mice.
Our data indicate a role of Gata6 in the physiopathology of sebaceous tumors, particularly in relation to DNA mismatch repair.

Discussion
Gata6 is widely expressed in the heart and in endoderm-derived tissues, including the lungs, liver, pancreas, stomach, and intestine (Maeda et al, 2005). Mice lacking Gata6 die during embryonic development as a result of endoderm defects (Morrisey et al, 1998). In contrast, Gata6 expression in the epidermis is limited to a population of SD/SG progenitors in the developing HF and to the JZ, upper SG, and part of the infundibulum in adult skin (Figs 2A and 6A, and EV2A-C; Donati et al, 2017).
Although Gata6 is strongly upregulated in the epidermis of K14DNLef1 mice (Fig 1H) and is a direct DNLef1 target gene (Fig 1C), we saw no evidence for co-expression of Gata6 and Lef1 in developing or adult epidermis and Gata6 was not induced upon b-catenin activation (Fig 2), even though Gata6 synergizes with or activates Wnt signaling in a number of contexts (Afouda et al, 2008;Zhang et al, 2008b;Whissell et al, 2014).
During skin morphogenesis, Gata6 was first expressed in the upper stage 4 HF (Fig 2A), consistent with a recent study that reported co-localization with Krt79 and Lrig1 (Mesler et al, 2017). The Gata6-positive population does not originate from Shh-positive progenitors (Mesler et al, 2017) and an alternative possibility is that it arises from Sox9 progenitors (Fig 2A). Regulation of Sox9 by Gata6 has been shown in several organs, including the pancreas (Carrasco et al, 2012) and cardiac valves (Gharibeh et al, 2018).
Gata6 gain and loss-of-function experiments revealed that Gata6 regulates JZ, SD, and upper SG differentiation in WT mice (Figs 3 and EV3B and C). In K14DNLef1 mice, a DNLef1-Gata6 ◀ Figure 7. Gata6 is mutated or downregulated in human sebaceous carcinomas with a high mutational burden.
ª 2019 The Authors The EMBO Journal 38: e100526 | 2019 transcriptional cascade triggers formation of multilayered cysts and ectopic SG (Fig 1A and B). Most of the cysts represent aberrant SD and JZ, since they express SD/JZ markers , rather than IFE as previously proposed (Niemann et al, 2002). Gata6 overexpression in utero led to formation of cysts, reminiscent of ducts, with ectopic SG features (Fig 3B and C). Our functional data together with our lineage tracing experiments (Figs 3D and EV3B and C) allow us to propose that the Gata6-positive cells in the developing HF give rise to the upper SG. The lower SG has a distinct origin that has yet to be identified . Gata6 is overexpressed in a variety of cancers, including pancreato-biliary cancers (Kwei et al, 2008), colon cancers (Shureiqi et al, 2007;Tsuji et al, 2014), esophageal adenocarcinomas (Lin et al, 2012), breast cancers (Song et al, 2015), and adrenal tumors (Vuorenoja et al, 2007). Conversely, it is lost in certain ovarian cancers (Cai et al, 2009) and may act as a tumor suppressor gene in lung cancer (Cheung et al, 2013) and astrocytoma (Kamnasaran et al, 2007). We observed a strong association between Gata6 and sebaceous carcinogenesis in mouse and human skin. Gata6 was expressed in sebaceous tumors of K14DNLef1 mice ( Fig 4A) and deletion of Gata6 reduced the proportion of tumors with sebaceous differentiation (Fig 4C). Furthermore, human skin tumors with sebaceous differentiation expressed Gata6 (Fig 6B and C). Thus, Gata6 could be a useful biomarker for diagnostic pathology of sebaceous tumors.
We observed that Gata6-positive cells were Ki67-negative in WT mouse back skin (Fig EV2B). However, we previously showed in K14DNLef1 mice that cell proliferation is a downstream effect of DNLef1 overexpression (Niemann et al, 2002). Ki67 staining is found in the outer root sheath and in some cells at the periphery of the dermal cysts. In addition, Ki67-positive cells are twice as frequent in the basal layer of K14DNLef1 as WT IFE (Niemann et al, 2002). Thus, Ki67 is mostly expressed in the Gata6-negative regions of K14DNLef1 epidermis. This suggests that Gata6 does not stimulate proliferation, consistent with previous observations . Therefore, two cell populations can be distinguished: actively proliferating K14 + /DNLef1 + /Gata6 À keratinocytes and infrequently proliferating K14 + /DNLef1 + /Gata6 + keratinocytes of the JZ and upper SG. K14DNLef1 mice develop sebaceous tumors at high frequency. This suggests that tumor initiation involves cooperation between these two cell populations. Our data indicate that Gata6 is likely to be responsible for the sebaceous differentiation observed in the tumors (Fig 4C) but also acts as a tumor suppressor (Fig 4B). In addition, we and others have confirmed the existence of proliferative cells in K14DNLef1 sebaceous tumors by showing Ki67 expression (Niemann et al, 2002), BrdU incorporation (Niemann et al, 2007), and isolating tumor-propagating cells that form secondary tumors in serial transplantation assays (Petersson et al, 2015). How the two cell populations cooperate to initiate tumors remains to be elucidated.
Sebaceous tumors are closely associated with MMR deficiency, for example, in the context of Muir-Torre syndrome (Eisen & Michael, 2009a,b;John & Schwartz, 2016). In this syndrome, mutations in the MMR genes Msh2 or, less frequently, Mlh1 and Msh6, predispose cells to DNA base errors (shown by the acquisition of MSI status; John & Schwartz, 2016). The estimated frequency of MSI is at least 60% of all sebaceous neoplasms, while only 3% of SG hyperplasia (Kruse et al, 2003;Jessup et al, 2016). Furthermore, Msh2-deficient mice develop sebaceous tumors (Reitmair et al, 1996).
We observed that Gata6 was required for expression of MMR genes. Microsatellite stability and expression of Msh2 and Mlh1 were reduced upon Gata6 knockout in K14DNLef1 mice (Figs 5C,D and E,and Appendix Fig S2). This would explain why tumor incidence was increased in K14DNLef1:cKO mice (Fig 4B). Paradoxically, while MSI-linked tumors develop and progress rapidly, their prognosis is often better. The increased mutation load resulting from MMR inactivation generates multiple neo-antigens and stimulates immune surveillance and cancer clearance (Germano et al, 2017). MMR-deficient tumors show an excellent response to pembrolizumab, an immunotherapy targeting the PD-1 receptor (Le et al, 2015). In addition to MMR regulation, Gata6 may also suppress tumor formation by reducing proliferation and increasing SD/SG differentiation .
Like K14DNLef1 mice, mice overexpressing DNLef1 under the control of the Krt15 promoter (which labels HF bulge cells) develop sebaceous tumors (Petersson et al, 2015). These tumors are more aggressive than those of K14DNLef1 mice and are associated with increased DNA damage (Petersson et al, 2015). These results could be explained by the lack of Gata6 expression in the bulge.
In humans, sebaceous tumors harbor a high level of expression of Gata6 (Fig 6B and C, and Appendix Fig S3). However, sebaceous carcinoma can also acquire Gata6 mutations correlating with the total mutational burden (Fig 7). Gata6mut tumors display reduced expression of several MMR genes (Fig 7G). In addition, a reduction in Gata6 expression is associated with less differentiated and more aggressive UV-related SebC, independent of GATA6 mutational status (Fig 7C).
The link between Gata6 and MMR that we have found could open up new therapeutic approaches to target sebaceous carcinogenesis. It is tempting to speculate that the MMR DNA repair pathway is necessary for maintenance of healthy SG, because lipid production is associated with the formation of reactive oxygen species (Bek-Thomsen et al, 2014;Ibrahim et al, 2014) and subsequent DNA damage.

Mice
All animal procedures were subject to local ethical approval and performed under a UK Government Home Office license (PPL 70/ 8474). K14DNLef1 (Niemann et al, 2002), Lef1 À/À (Van Genderen et al, 1994), K14DNb-CateninER (D2 and D4 lines;Lo Celso et al, 2004), K14Cre/bCat Flox(ex3)/+ (Zhang et al, 2008a), epidermal Gata6 conditional knockout (cKO) , Gata6EGFPCreERT2 ), Lgr6EGFPCreERT2 (Snippert et al, 2010, and Rosa26-fl/STOP/fl-tdTomato (Madisen et al, 2010) mice were previously described. The K14DNb-CateninER transgene was activated by one or six topical applications of 1.5 mg 4-hydroxytamoxifen (4OHT) (Sigma) . For lineage tracing experiments, pregnant females were injected intraperitoneally with a dose of 50 lg/g of tamoxifen (Sigma) at E16.5 and E18.5. Samples were collected at 2 or 17 days after recombination. Ultrasound-guided lentiviral in utero injection procedures were performed as previously described (Beronja et al, 2013). Lentiviral particles were produced by transfecting HEK293 cells with a Trans-Lentiviral Packaging System in combination with Precision LentiORFs control and Gata6 (Dharmacon). GFP chimeric mice were generated as previously described (Arwert et al, 2010). For skin carcinogenesis experiments, K14DNLef1 and K14DNLef1:cKO mice received a single 100 nmol dose of DMBA (7,12-dimethylbenz(a)anthracene; Niemann et al, 2007). Tumor incidence and burden were assessed once a week by two independent researchers. As a control for Appendix Fig S2A and B, a WT mouse was UVB-irradiated with a total dose of 500 mJ/cm 2 . No specific method for randomization, blinding, or estimation of sample size was used. Male and females were used. All efforts were made to minimize suffering of mice.

Human tissue
All human tissue samples were collected, diagnosed, and processed for research in accordance with the recommendations of the relevant local ethics committees in compliance with the UK Human Tissue Act and approved by the National Research Ethics Service (08/H0306/30), German Medical Council, and/or the Japanese Ministry of Health, Labor, and Welfare. Human embryonic and fetal tissues were obtained with appropriate ethical approval from the UK Human Developmental Biology Resource (www.hdbr. org). Informed consent was obtained from all subjects. The experiments conformed to the principles set out in the WMA Declaration of Helsinki and the Department of Health and Human Services Belmont Report.

ChIP-Seq, microarray, and computational analysis
ChIP-Seq was performed as previously described (Mulder et al, 2012). Briefly, crosslinked material corresponding to $ 10 7 K14DNLef1 cultured primary keratinocytes (collected from 9.5week-old mice) was incubated overnight with 10 lg of Lef1 (Cell Signaling) or Flag antibody (Sigma), and immune-precipitated DNA was sequenced. Base-calling, genome alignment, filtering against potential PCR duplications and peak calling were performed as previously described .
RNAs from FACS-purified Cd34 + Itga6 + and Cd34 À Itga6 + epidermal cells collected from WT and K14DNLef1 epidermis (9.5-weekold mice) were provided to the Paterson Institute Microarray Core Facility to perform gene expression profiling using the mouse Exon1.0ST Affymetrix platform. Differential expression analysis was carried out on normalized data as described previously .
Gene Set Enrichment Analysis was performed to compare expression signature genes in SG vs IFE and HF (ranked Z-score values of differentially expressed genes (DEG); Donati et al, 2017) with respect to a gene subset composed of the nearest genes to the DNLef1 peaks from Chip-Seq data.
Human Gata6 ChIP-Seq datasets were interrogated using the "Target Genes" module of ChIP-Atlas (Oki et al, 2018 and http:// chip-atlas.org). This module predicts genes directly regulated by a given protein, based on binding profiles of all public ChIP-Seq data for particular gene loci.

Cell culture and siRNA transfection
Primary mouse keratinocytes were isolated from dorsal skin as previously described (Jensen et al, 2010) and cultured on confluent irradiated 3T3-J2 fibroblast feeders in calcium-free FAD medium (DMEM: Ham's F12, 3:1, 1.8 × 10 À4 M adenine) supplemented with 10% fetal calf serum, hydrocortisone (0.5 lg ml À1 ), insulin (5 lg ml À1 ), cholera toxin (8.4 ng ml À1 ), and epidermal growth ª 2019 The Authors The EMBO Journal 38: e100526 | 2019 factor (10 ng ml À1 ). siRNAs for the negative control and mouse Lef1 (Ambion) were introduced into cells by nucleofection using the Amaxa 96-well shuttle system (Lonza) as previously described (Mulder et al, 2012). The SebE6E7 sebocyte line was obtained and cultured as described previously (Lo Celso et al, 2008). All cell stocks were routinely tested for mycoplasma contamination and were negative. SebE6E7 cells were transfected with a control plasmid or human full-length or DN34Lef1 plasmids (Takeda et al, 2006) using jetPRIME transfection reagent (Polyplus transfection) following the manufacturer's instructions.

RT-qPCR
For RT-qPCR, total RNA was isolated using RNeasy kits (Qiagen). cDNAs were generated using the SuperScript III Supermix (Invitrogen) or QuantiTect Reverse Transcription Kit (Qiagen) and analyzed using Power SYBR Green (Applied Biosystems) or SYBR Green Master Mix (Life Technologies) and custom-made primers.

Microsatellite instability analysis
Microsatellite instability in sebaceous tumors was determined as previously published (Woerner et al, 2015;Germano et al, 2017). DNA was extracted from paraffin-embedded tumors using a Relia-Prep FFPE gDNA Miniprep System Kit (Promega). Amplification was performed on 20 ng DNA with a Type-it Microsatellite PCR Kit (Qiagen). The cycling profile was as follows: 1 cycle at 95°C for 5 min, 95°C for 30 s, 56°C for 90 s, and 72°C for 30 s for a total of 28 cycles; final extension at 60°C for 30 min. The following labeled primers were used: mBat64, forward fluorescein-GCCCACACTCCTGAAA ACAGTCAT and reverse CCCTGGTGTGGCAACTTTAAGC; AC09 6777, forward JOE-TCCCTGTATAACCCTGGCTGACT and reverse GCAACCAGTTGTCCTGGCGTGGA; AA003063, forward Tamra-ACGTCAAAAATCAATGTTAGG and reverse CAGCAAGGGTCCC TGTCTTA; U12235, forward JOE-GCTCATCTTCGTTCCCTGTC and reverse CATTCGGTGGAAAGCTCTGA; and L24372, forward fluorescein-GGGAAGACTGCTTAGGGAAGA and reverse ATTTGGCTTT CAAGCATCCATA. PCR fragments were separated using Applied Biosystems Big-Dye Ver 3.1 chemistry on an Applied Biosystems model 3730 automated capillary DNA sequencer. Raw data were visualized with Geneious software and further analyzed as previously described (Keysselt et al, 2017).

Statistics and reproducibility
Statistical analyses were performed using GraphPad Prism software. No statistical method was used to predetermine sample size. Panels showing images are representative of at least two independent experiments as indicated in each panel of the figure legends.
Expanded View for this article is available online.