Cell-Type-Specific Chromatin States Differentially Prime Squamous Cell Carcinoma Tumor-Initiating Cells for Epithelial to Mesenchymal Transition.

INTRODUCTION
EMT is associated with cancer metastasis, tumor sternness, and resistance to
therapy (Mani et al., 2008; Nieto et al., 2016; Yang
et al., 2004). While the cancer cell of origin has been suggested to
control tumor heterogeneity, no study has demonstrated so far that the cancer cell
of origin controls EMT (Nieto et al., 2016).
Depending on the cancer cell of origin (multipotent and unipotent stem cells,
progenitors, and differentiated cells) initially targeted by oncogenic hits,
different tumor phenotypes may arise, differing by their differentiation,
aggressiveness, and EMT features.
The skin epidermis is an ideal model to assess whether the cancer cell of
origin controls EMT, as it is composed of spatially distinct cell lineages including
the interfollicular epidermis (IFE), the hair follicle (HF), and its associated
sebaceous glands, as well as the infundibulum that connects the HF to the IFE (Blanpain and Fuchs, 2014) (Figure 1A). During homeostasis, each of these distinct
epidermal lineages is self-sustained by its own pool of resident stem cells (SCs)
that can be genetically targeted by specific CreER mice (Blanpain and Fuchs, 2014), allowing the conditional
expression of oncogenes or deletion of tumor suppressor genes in different epidermal
lineages and the assessment of their capacity to induce tumor formation (Blanpain, 2013). In studying the cellular origin
of skin SCCs, the second most frequent skin cancer in humans, it has been previously
demonstrated that oncogenic KRas expression combined with p53 deletion in IFE cells
as well as in HF lineages leads to the development of different types of invasive
SCCs, sometimes associated with EMT features, demonstrating that different epidermal
lineages including the IFE and the HF were competent to induce skin SCCs (Lapouge et al., 2011; White et al., 2011). However, it remains unclear to what
extent the cellular origin of skin SCCs influences EMT in these tumors.
Here, we used genetically engineered mouse models coupled with lineage
tracing to assess whether the same oncogenic hits in different cell lineages of the
skin epidermis influence EMT. Interestingly, HF-derived tumors are much more prone
to undergo EMT as compared to IFE-derived tumors. Chromatin and transcriptional
profiling of these two different epidermal populations during tumorigenesis revealed
that the epigenetic and transcriptional landscapes of the cancer cell of origin
primed oncogene-targeted cells to develop into either well-differentiated SCCs or
more invasive tumors characterized by EMT, underscoring the importance of the cancer
cell of origin in controlling EMT.

RESULTS

The Cancer Cell of Origin Controls EMT in Skin SCC
To determine whether the cancer cell of origin controls EMT in skin
tumors, we assessed the tumor phenotypes following KRasG12D
expression and p53 deletion either in the IFE using K14CreER mice
(K14CreER/KRasG12D/p53fl/fl/Rosa-YFP), in which
low-dose tamoxifen (TAM) administration preferentially targets the IFE and the
infundibulum (Lapouge et al., 2011) or in
HF SCs and their progeny using Lgr5CreER mice
(Lgr5CreER/KRasG12D/p53fl/fl/Rosa-YFP) (Lapouge et al., 2012) Figures 1A–1C and S1A–S1D). These two CreER targeted exclusively epidermal
cells and not mesenchymal cells in the skin (Figures S1A–S1D). The
Rosa-YFP reporter gene was included in these different mouse models to track on
tissue sections and isolate by fluorescence-activated cell sorting (FACS)
YFP+ tumor cells, including those that may have lost the
expression of markers of their cell of origin (Figure 1B). As previously described (Lapouge et al., 2011), the kinetic of tumor appearance (around
6–9 weeks) and initial tumor growth rate were similar in K14-and
Lgr5-derived tumors. However, Lgr5CreER mice developed more tumors than K14CreER
mice (Figures
S2E–S2G). Histological analyses revealed that K14-derived
tumors were mostly well-differentiated SCCs, containing numerous keratin pearls.
In contrast, tumors arising from Lgr5CreER/KRasG12D/p53fl/fl/Rosa-YFP
mice were composed of distinct phenotypes (Figure
1D). The most frequent tumor phenotype consisted of mixed tumor
containing YFP+ tumor epithelial cells (TECs) and
YFP+ tumor mesenchymal-like cells (TMCs), resembling
human carcinosarcoma for which the epithelial or the mesenchymal origin remains
an open question (Paniz-Mondolfi et al.,
2015). The other tumors included well-differentiated SCCs resembling
the SCCs arising in K14CreER mice and tumors that were completely mesenchymal
and called spindle cell carcinoma (Figure
1D). Well-differentiated SCCs, regardless of whether they originate
from K14CreER or Lgr5CreER cells, expressed classical epidermal markers
including K14, Epcam, and E-cadherin, and no YFP+ TECs
expressed the mesenchymal marker Vimentin, which was only expressed by
YFP-negative cancer-associated fibroblasts that composed the tumor stroma (Figures 1E, 1F, and S1H). In Lgr5-derived mixed SCCs, the well-differentiated part of the
tumor expressed epithelial markers K14, Epcam and E-cadherin (Figures 1E and S1H). However, many cells located
in the underlying mesenchymal part of the tumor were YFP+
TMCs that had completely lost the expression of all classical epithelial markers
and expressed Vimentin (Figures 1E, 1F, and S1H). Likewise, in Lgr5-derived
mesenchymal tumors, YFP+ TMCs and mesenchymal stromal cells
were morphologically indistinguishable and were all K14, Epcam, and E-cadherin
negative (Figures 1E, 1F, and S1H). FACS analysis of Epcam
expression showed that the vast majority of YFP+ TECs in
well-differentiated SCCs, whether K14 or Lgr5-derived, expressed high levels of
Epcam (Figure 1G). Epcam expression was
always bimodal in mixed tumors with two distinct YFP+
Epcam+ TEC and YFP+
Epcam− TMC populations, while in mesenchymal tumors all
YFP+ TMCs were Epcam− (Figure 1G). FACS quantification of the proportion of
YFP+ Epcam+ TECs and
YFP+ Epcam− TMCs in a large number of
tumors showed that the majority of K14-derived tumors (75%, n =
63) were well-differentiated SCCs composed of Epcam+ cells
(Figures 1H and 1I). In contrast, the vast majority of Lgr5-derived
tumors (92%) were composed of cells that had partially (mixed tumors)
(65%, n = 224) or completely (mesenchymal tumors) (27%,
n = 94) undergone EMT (Figures 1H
and 1I). A similar degree of EMT was also
observed upon KrasG12D and p53 recombination using two other
HF-specific inducible CREs (K15CrePR and K19CreER) (Lapouge et al., 2011) (Figure S1I), further demonstrating
that the mesenchymal part of the tumor arises from the HF lineages. RT-PCR and
PCR analysis of the KrasG12D and p53 floxed alleles showed that KRas
was expressed at similar level in TECs and TMCs and that both alleles were
equally recombined in the different types of SCCs, showing that tumor
heterogeneity did not arise from different levels of oncogene expression or
recombination (Figures
S1J–S1L). To assess whether Lgr5-derived mixed tumors arise
from a clonal event, we induced mice with a lower dose of TAM (1 mg), which
strongly decreased the number of tumors per mouse (three versus 20) and thereby
the chance that two tumor-initiating clones fused together into a single tumor
mass. Despite the reduction of the number of tumors, these mice presented a
similar proportion of mixed tumors (Figures S2A–S2D). Using
multicolor confetti reporter mice, tumors labeled with one of the four colors
contained both Epcam+ and Epcam− labeled cells
expressing the same fluorescent protein (Figures S2E and S2F), further
supporting that a single tumor-initiating cell gave rise to a mixed tumor.
Altogether, these data demonstrate that KrasG12D expression and p53
deletion in HF SCs and their progeny induce SCCs that undergo spontaneous EMT
and indicate that the cellular origin of cancer dictates EMT in skin tumors.
We then assessed whether EMT is a progressive process that arises during
tumor progression or whether EMT occurs at the early step of tumorigenesis, by
analyzing the skin of the mice before the appearance of macroscopically visible
tumors. Small microscopic Lgr5-derived tumors were already composed of
K14+ and K14-negative cells (Figures S2G and S2H). In addition,
the proportion of Epcam+ cells in Lgr5-derived SCCs did not
correlate with their size (Figure S2I), further demonstrating that EMT can occur during the
early step of tumorigenesis, as it has been recently shown in pancreatic tumors
(Rhim et al., 2012).

Intrinsic Priming of HF to Undergo EMT Promoting Clonogenic and Metastatic
Potential
To define whether in the absence of their natural environment and
surrounding stromal cells, tumor cells are still biased to undergo EMT depending
on their cellular origin, we assessed the potential of Epcam+
TECs to undergo EMT in transplantation assays. To avoid the influence of the
tumor microenvironment, we transplanted freshly FACS purified tumor cells with
Matrigel subcutaneously into immunodeficient mice and assessed their ability to
reform secondary tumors that recapitulate the histology and heterogeneity of the
primary tumors (Beck and Blanpain, 2013)
(Figure 2A). Limiting dilution analyses
of transplanted Epcam+ TECs from K14 and Lgr5-derived tumors showed that
both TECs presented similar tumor propagating frequency irrespective of their
cell of origin (Figure 2B). Transplantation
of a single tumor Epcam+ TEC or their progeny gave rise to
mixed tumors containing TECs and TMCs, showing their ability to recapitulate the
histology of the tumor of origin at the clonal level (Figures S3A–S3L). The
histology and the proportion of Epcam+ TECs and
Epcam− TMCs from K14 and Lgr5 tumors were very similar to
the histology of the primary tumors, with K14CreER mice giving rise
preferentially to well-differentiated tumors containing high proportions of
Epcam+ TECs, whereas Lgr5-derived tumors were mostly
mixed or mesenchymal (Figures
2C–2G). These data
indicate that the cancer cell of origin influences the intrinsic ability of
tumor-initiating cells to undergo EMT.
Overexpression of EMT transcription factors (TFs) in cancer cell lines
promotes tumor renewal capacities in vitro and their ability to form secondary
tumors upon transplantation into immunodeficient mice (Puisieux et al., 2014), which suggested that EMT
promotes tumor stemness. To determine whether naturally occurring EMT in primary
skin tumors is associated with tumor stemness, we assessed TECs and TMCs
clonogenic potential in vitro and in vivo. The in vitro colony forming
efficiency and total cell output were higher in Epcam− TMCs
as compared to Epcam+ TECs (Figures S3M–S3T), showing
that spontaneous EMT in SCCs is associated with increased in vitro
clonogenicity. Transplantation experiments showed that Epcam−
TMCs were more efficient at forming secondary tumors (Figure 2B), although these secondary tumors were
exclusively composed of mesenchymal cells, showing that EMT cells, although
presenting enhanced renewal cannot revert back to an epithelial state following
subcutaneous transplantation (Figures 2F
and 2G).
To assess whether EMT increases the metastatic potential of primary skin
tumors, we injected intravenously TECs derived from K14 tumors and TECs and TMCs
coming from Lgr5 tumors and assessed the number of metastasis 4 weeks later. The
number of metastasis was much higher in mice injected with TMCs compared to TECs
(Figures 2H and 2I), demonstrating the greater metastatic potential of
skin tumor cells that have undergone EMT.
Altogether these data show that HF-derived tumors are intrinsically
primed to undergo EMT and EMT is associated with increased clonogenicity, tumor
propagation, and stemness as well as a higher capacity to give rise to
metastasis, demonstrating the functional importance of the cancer cell of origin
in primary skin SCCs.

Transcriptional Priming of EMT Genes in HF Lineages
To define the mechanisms that promote EMT in primary skin tumors, we
first established the transcriptional signature of TECs and TMCs in vivo. To
this end, we performed microarray analysis of FACS-isolated
YFP+ Epcam+ and
Epcam− tumor cells arising from IFE and HF lineages. Gene
ontology of the genome-wide transcriptional analysis confirmed by qRT-PCR and
immunostaining showed that K14 and Lgr5 derived Epcam+ TECs
were enriched for transcripts regulating the epithelial state, including
epidermal markers such as K14, K6, K1, and many TFs known to promote epithelial
adhesion, differentiation, and stratification such as p63
(Mills et al., 1999; Yang et al., 1999), Ovol (Lee et al., 2014; Watanabe et al., 2014), Grhl (Cieply et al., 2013; Hopkin et al., 2012; Xiang et al., 2012), Cebpa (Lopez et al., 2009), and Klf5 (Kenchegowda et al., 2011) TFs,
Esrp splicing factors that promote epithelial
differentiation and inhibit EMT (Warzecha et
al., 2010) as well as different cancer SC markers of
well-differentiated SCCs such as Sox2, Vegfa, and
CD133 (Beck et al.,
2011; Boumahdi et al., 2014)
(Figures 3A, 3B, and S4A).
In contrast, TMCs were strongly enriched in secreted proteins including
components of extracellular matrix (ECM), cell adhesion, ligands, and inhibitors
of developmental signaling pathways, regulators of angiogenesis, and various
types of mesenchymal collagens (Figure 3C).
RT-PCR and immunostaining confirmed the preferential expression of classical EMT
markers (Nieto et al., 2016) such as
EMT-related TFs Snai1, Zeb1/2, Twist1, and
Prrx1, EMT-related cytoskeleton (Vimentin)
and adhesion molecules (Cdh2/N-cadherin, Cadherin 11), and
specific ECM genes such as Col3a1 and Fn1
(Figures 3D and S4B).
To assess whether the gene expression in the cancer cell of origin
influences the type of tumors that isformed following oncogenic transformation,
we performed transcriptional profiling of normal IFE/infundibulum cells
(Lgr5-GFP negative, α6HiCD34−) and
Lgr5-GFP+ HF cells and compared the genes differentially
expressed between normal IFE and HF cells with the genes preferentially
expressed in TECs and TMCs they give rise to. Interestingly, 29% of the
of genes upregulated in TECs (414/1442) were already upregulated at the mRNA
level in IFE compared to HF (Figure 3E).
Likewise, 27% of genes upregulated in Epcam− versus
Epcam+ tumor cells (315/1148) were already upregulated at
the mRNA level in HF versus IFE (Figure
3F). Gene set enrichment analysis (GSEA), a bioinformatic approach that
takes into account all genes and their level of expression, showed a very strong
enrichment (normalized enrichment score [NES] >5) of IFE
genes within the TEC signature and of HF genes within TMC signature (Figures S4C and S4D),
further supporting the notion that transcriptional priming influences the
ability of the cancer cell of origin to undergo EMT during tumorigenesis.
The genes transcriptionally primed in IFE lineage and associated with
well-differentiated tumors comprised well-known transcription regulators of
epidermal differentiation such as Grhl1/3, Cebpa, Klf5, Ovol1,
and Pou3f1 (Cieply et al.,
2013; Hopkin et al., 2012;
Kenchegowda et al., 2011; Lee et al., 2014; Lopez et al., 2009; Watanabe et al., 2014; Xiang et al.,
2012), as well as genes associated with epithelial adhesion and
epidermal differentiation such as Klk5, 8,10, 11, Sprr1a, or
Tgm3 (Figure 3G). In
contrast, the genes primed in HF lineages and that are associated with EMT
included well-known HF markers such as Ltbp2, Grem1, Flstl1, S100A4,
Nfatc1, Tbx1, Tcf4, Tcf7l1, and
Ctgf (Chen et al.,
2012; Horsley et al., 2008;
Morris et al., 2004; Nguyen et al., 2009; Tumbar et al., 2004) and mesenchymal genes such as Col5a2,
Col6a1, Fn1, and MMP11 (Figure 3H).

Chromatin Landscape Remodeling during Skin Tumorigenesis
We next defined more globally the changes in the chromatin landscape
that occur during tumorigenesis and EMT and assessed whether and how the cancer
cell of origin is epigeneti-cally primed to give rise to different tumor
phenotypes. To this end, we performed Assay forTansposase-Accessible Chromatin
sequencing (ATAC-seq), a technique that allows the mapping of the open chromatin
regions with extremely high definition and with a very low amount of cells
(Buenrostro et al., 2013), on
FACS-isolated HF and IFE and their respective tumor cell populations, allowing
to define the open chromatin regions and the TFs associated with the remodeling
of these chromatin regions during tumorigenesis and EMT (Figure 4 and Table S1A). We first defined the
chromatin remodeling associated with tumorigenesis, by assessing the chromatin
regions (ATAC-seq peaks) that are changed by more than 3-fold, a stringent
threshold, between the cell of origin and the respective tumor populations. We
identified more open chromatin regions in 477 genes that are upregulated during
tumorigenesis across all the different tumor cell populations analyzed (K14
TECs, Lgr5 TECs, Lgr5 TMCs) (Figure 4A;
Table S1A), which
represent the common epigenetic and transcriptional changes associated with
tumor initiation irrespective of the cancer cell of origin and EMT. These genes
included ligands of the EGFR pathway (Areg, Ereg) and TFs known
to relay the activation of the Ras/MAPK pathway (Fos/Fosb,
Fosl1, Nfe2l2, Ets1) and other TFs promoting tumor
sternness and invasion (Twist1, Hmga2, Prrx1) (Beck et al., 2015; Copley et al., 2013; Ocahaet al., 2012)(Figures 4B and 4C). Interestingly, the chromatin regions within some of the EMT genes,
such as Snai1 and Zeb2, were already opened in
Epcam+ tumor cells despite the lack of protein expression
(Figures 4D and 4E), suggesting that the EMT program is epigenetically
primed in TECs. Motif enrichment analysis of the chromatin regions that opened
during tumorigenesis revealed a strong enrichment for the binding site of TFs
such as Jun/AP1 (65%), Ets1 (37%), Runx (29%), Nf-kb
(22%), and TEAD (25%) (Figure
4F).

Epigenetic Priming of the Cancer Cell of Origin to Undergo EMT
To define the chromatin remodeling that occurred during EMT, we assessed
the chromatin regions that were upregulated or downregulated by more than 3-fold
between TECs and TMCs from Lgr5-derived tumors (Figures 5A and 5B). GSEA
analysis showed that the opening of the chromatin was mainly associated with
gene activation while chromatin closing was associated with gene repression
during EMT (Figures S5A and
S5B). To get further insights into the gene regulatory network (GRN)
that controls EMT in primary skin tumors, we searched the TF motifs enriched in
the chromatin regions that are opened or closed during EMT. The motifs with the
highest statistical significance upregulated during EMT and positively
associated with gene expression in TMCs were Jun/AP1 (42%), NF1
(45%), Ets1 (10%), bHLH TFs (20%–45%),
Nfatc (27%), and Smad2 (37%) (Figure 5C). TFs binding these motifs (e.g., Runx1/2, Nfatc1/2,
Twist1/2, Tcf4) were upregulated in Epcam− TMCs (Figure S5C; Table S1B). The same core
of TFs including Jun/AP1 (54%), NF1 (48%), Ets1 (29%)
was also enriched in the more open chromatin regions in
Epcam+ TECs (Figure
5D). In addition to this core of TFs, a different set of specific TF
motifs including p63 (34%), Grhl (43%), Klf4/5 (18%),
Cepba (20%), orSox2 (25%) was highly enriched in the more open
chromatin regions in Epcam+ TECs (Figures 5D and S5D; Table S1C). While Zeb1 motif was
not particularly enriched in chromatin regions of the genes that were
upregulated in TMCs, Zeb1 motif was strongly enriched in chromatin regions that
become closed during EMT, supporting the notion that Zeb1 mainly acts as a
transcriptional repressor (Eger et al.,
2005). Smad2 motif was also strongly enriched in chromatin peaks that
are closed in TMCs (42% of the peaks), suggesting that transforming
growth factor β (TGF-β)/Smad2 axis acts both as transcriptional
activator and repressor during EMT (Figure
5D).
To assess whether the epigenetic landscape of the cancer cell of origin
promotes or restricts gene expression and EMT in skin tumors, we compared the
chromatin landscape of the cancer cell of origin to the tumor cells they
derived. Interestingly, we found that 139 genes transcriptionally upregulated
during EMT already presented more open chromatin regions in HF, suggesting that
these genes are epigenetically primed in HF to facilitate EMT during
tumorigenesis (Figures 5E, S5E, and S5F). These genes
including Tcf7l1, Ltbpl, and Zeb1 were
enriched for NF1 binding site (Figures 5F
and 5G), some of which share the same
active enhancers in HFs and TMCs (Figure S5G). About half
(45%) of the epigenetically primed genes were already upregulated in HF
cells (Figure S5F;
Table S1D). These
data indicate the chromatin landscape of HF primes these cells to undergo EMT
during tumorigenesis. We next analyzed whether the chromatin landscape of the
IFE primes these cells to give rise to well-differentiated tumors (Figure S5H).
Interestingly, 253 genes transcriptionally upregulated in
Epcam+ TECs presented more open chromatin regions in IFE
(Figures 5H and 5I), among which 112 genes were already
transcriptionally upregulated in IFE (Figure S5I; Table S1E). The common open
chromatin regions between IFE cells and Epcam+ TECs included
the regulatory regions of the key TFs promoting epidermal differentiation such
as p63, Klf5, Grhl1/3, or Cepba (Figures 5I and S5J). These commonly upregulated
chromatin regions between IFE and Epcam+ TECs were highly
enriched for Jun/AP1 (58%), p63 (43%), and Klf (22%)
binding sites (Figure 5J), suggesting that
p63 and Klf TFs contribute to the epigenetic priming of IFE cells to give rise
to well-differentiated tumors upon oncogenic transformation.

Activator and Repressor Functions of Epithelial and Mesenchymal TFs
To functionally challenge the bioinformatic predictions of the GRNs that
positively and negatively control EMT in primary skin SCCs (Figures 6A, 6B,
S6, and S7), we performed short
hairpin RNA(shRNA) knockdown (KD) of TFs predicted to control the GRNs
associated with EMT in vivo and assessed the impact of their KD on the
regulation of epithelial and mesenchymal gene expression. FACS-isolated
Epcam+ and Epcam− cells were cultured
on irradiated 3T3 feeder layers, and the tumor cells were infected with
lentiviruses expressing shRNA against p63, Klf5, Twist1, Zeb1,
or Snai1, the most statistically significant TF motifs found in
the open or close chromatin regions during EMT (Figures 6C and 6D).
Consistent with the high enrichment of Klf5 and p63 motifs in the open
chromatin regions of the genes associated with TEC state, shRNA against
p63 or Klf5 in TECs resulted in the
downregulation of key epidermal TFs (Grhl2, Cepba, and
p63) and epithelial markers (Ecadh, K5)
that presented p63 and Klf5 motifs in their open chromatin regions, supporting
the notion that p63 and Klf5 act as transcriptional activators of gene
expression (Figure 6C). The EMT genes
upregulated following p63 and Klf5 KD such as
Zeb were not enriched in p63 and Klf5 motifs in their open chromatin regions,
suggesting that this effect was indirect (Figure
6C). In contrast to the transcriptional activator function of the key
epithelial TFs, the classical EMT TFs Zeb1, Snai1, and Twist1 mainly act as
transcriptional repressors, as suggested by the presence of their binding sites
in the chromatin regions that are closed during EMT, and the upregulation of key
epidermal TFs such as p63, Cepba, and Grhl2
and epithelial genes such as Ecadh and K5
following their shRNA KD in vitro (Figure
6D). In addition to repressing gene expression, these EMT TFs could
also positively regulate the expression of some mesenchymal genes (e.g., Cdh11
by Zeb1). Altogether, the consequences of shRNA KD of several key TFs predicted
to regulate the expression of genes during EMT are consistent with the
bioinformatic prediction and support the notion that these key epithelial TFs
act mainly as transcriptional activators, whereas the classical EMT TFs act
mainly as transcriptional repressors.
TGF-β-induced EMT is the most commonly used method to promote
EMT in cancer cells (Moustakas and Heldin,
2014). To further challenge the bioinformatic prediction of the GRN
associated with EMT, we assessed how p-Smad2 activation mediated by
TGF-β signaling regulates the expression of the genes containing Smad2
binding sites in the chromatin regions that are remodeled during EMT. RNA-seq
performed 48 hr following the addition of TGF-β to
Epcam+ Lgr5-derived TECs showed that 20% of the
upregulated genes and 30% of the downregulated genes that contained
chromatin regions with Smad2 motifs remodeled during EMT were differentially
regulated by TGF-β treatment (Figures
6E and 6F), showing that
TGF-β/Smad2 axis activates and represses gene expression during EMT as
predicted by our GRN models.

ΔNp63 Primes the Cancer Cell of Origin toward Well-Differentiated
Tumors
Our bioinformatic analysis suggests that p63 regulates the epigenetic
and transcriptional priming of IFE cancer cell of origin toward
well-differentiated SCCs. Immunostaining for p63 in the different SCC subtypes
showed that p63 was only expressed in TECs and not in TMCs, consistent with this
possibility (Figure 7A). To test this
hypothesis, we generated
Lgr5CreER/KRasG12D/p53fl/fl/Rosa-ΔNp63-IRES-GFP
mice (Figure 7B), allowing constitutive
expression of ΔNp63 together with KRasG12D activation and p53 deletion
in HF lineages to assess whether the sustained expression of p63 restricts EMT
in HF-derived tumors. Whereas the number of tumors that formed in homozygous
ΔNp63-IRES-GFP mice was unchanged (Figure
7C), the proportion of well-differentiated SCCs was strongly
increased in Lgr5-derived tumors that expressed ΔNp63 (Figures 7D and 7E).
To further define the molecular mechanisms by which ΔNp63
promotes the well-differentiated tumor phenotypes, we performed transcriptional
profiling of FACS-isolated GFP+ Epcam−
tumor cells from Lgr5CreER/KRasG12D/p53fl/fl/Rosa
ΔNp63-IRES-GFP-induced mice. We found that the sustained expression of
ΔNp63 induced the upregulation of 566 genes in Epcam−
cells, among which 329 genes (58%) belong the TEC signature (p <
10−16) (Figure 7F;
Table S1F). GSEA
analysis further demonstrated the major enrichment of genes of the TEC signature
in Epcam− cells overexpressing ΔNp63 (Figure 7G). Many of these genes contained p63 motifs
in the open chromatin regions of Epcam+ TECs and IFE cells,
supporting the notion that p63 directly regulates their expression during
tumorigenesis (Figure 7H).
To functionally assess whether p63 restricts EMT induced by extrinsic
cues, we assessed the ability of TGF-β to induce EMT in K14 and
Lgr5-derived TECs overexpressing ΔNp63. K14-derived TECs were much more
resistant to TGF-β-induced EMT as compared to Lgr5-derived TECs (Figure 7I), further demonstrating the
intrinsic priming of Lgr5-derived TECs toward EMT as compared to IFE-derived
TECs. Interestingly, Lgr5 TECs overexpressing p63 were also much more resistant
to TGF-β-induced EMT as compared to Lgr5 TECs that did not overexpress
p63 (Figure 7I), showing that p63 opposes
TGF-β-induced EMT.
These data demonstrate that ΔNp63 acts functionally and
molecularly as a master regulator of the TEC fate in vivo and primes the IFE
cancer cell of origin toward well-differentiated SCCs.

DISCUSSION
Here, using lineage-tracing experiments allowing the expression of the same
oncogenic hits in different epidermal lineages, we showed for the first time that
cancer cell of origin controls EMT in skin SCCs. The transplantation of
FACS-isolated TECs in the absence of their underlying microenvironment lead to the
same biases to undergo EMT, suggesting that intrinsic factors promote EMT in
oncogene-targeted HF cells. Likewise, HF-derived TECs are much more sensitive to
TGF-β-induced EMT, further demonstrating that intrinsic factors and higher
propensity to respond to TGF-β signaling promote EMT in HF-derived
tumors.
Our genome-wide transcriptional analysis of the cancer cell of origin and
their corresponding tumor cells during skin tumorigenesis allows defining the
molecular determinants that regulate EMT in primary tumors in vivo. In addition to
these novel in vivo EMT transcriptional signatures, we unravel for the first time
the changes in the chromatin landscape occurring during cancer initiation and EMT in
primary tumors in vivo that provide novel insights into the mechanisms by which the
cancer cell of origin regulates EMT during tumorigenesis in vivo (Figure S7).
Bioinformatic analysis of the chromatin remodeling during tumorigenesis
reveals a core of TFs promotes gene expression during tumorigenesis independently of
the cell of origin and EMT. The most enriched motif (60%) in these common
open chromatin regions of tumor cells is API, which corresponds to the binding site
of Jun/Fos TFs. Transcriptional analysis revealed that several members of the
Jun/Fos family were upregulated in tumor cells as compared to their cell of origin,
including Fosll, Fosb, Fos, and Junb, consistent
with the well-known functions of Jun/Fos family members in relaying
transcriptional regulation downstream of oncogenic Ras/MAPK activation (Eferl and Wagner, 2003). Motifs for Ets,
another family of TFs known to relay oncogenic Ras/MAPK signaling in skin tumors
(Dittmer, 2015; Yang et al., 2015), and for Runxl, Tead, and Nfkb TFs,
which regulate skin tumorigenesis in vivo (Hoi et
al., 2010; Zanconato et al., 2015;
Zhu et al., 2009), were also highly
enriched in the open chromatin of all tumor cells.
In addition to these core TFs, lineage-specific TFs regulate the specificity
of the tumor phenotypes and EMT in skin tumors. P63, Klf, Grhl, and
Cepba are upregulated in TECs, and their binding sites are
enriched in the TEC open chromatin regions, suggesting they regulate the GRN that
promote squamous differentiation in oncogene-targeted IFE cells. Klf5 and p63 are
the most significant TFs associated with the priming of IFE toward
well-differentiated SCCs. KLF5 is mutated and acts as a genetic
driver in lung carcinomas, including lung SCCs (Campbell et al., 2016). Super-enhancer amplification associated with
increased KLF5 expression has been recently found in head and neck
SCCs (Zhang et al., 2016), and Klf5 prevents
apoptosis induced by TGF-β signaling (David
et al., 2016). P63 is the master regulator of epidermal stratification,
promotes SC renewal in stratified epithelia and is expressed in different human
carcinomas (Melino et al., 2015). Our
functional in vitro and in vivo gain- and loss-of-function experiments showed that
Klf5 and p63 act as master regulators of the epithelial state and prime the IFE
cells into well-differentiated tumors upon oncogenic Ras expression and cooperated
with the core TFs such as AP1 and Ets, to positively regulate gene expression in
TECs. By promoting the expression of micro-RNAs such as miR-200, which target key
EMT TFs (Wellner et al., 2009), p63 and Klf5
also participate to the down-regulation of EMT genes such as Zeb1 (Zhang et al., 2013; Zhao
etal., 2016).
During EMT, the core TFs cooperate with EMT-specific TFs including bHLH,
Runx, Nfat, and Smad2 TFs, several of which are well known to be preferentially
expressed in normal HF SCs and to regulate their function (Blanpain et al., 2004; Morris et al., 2004; Tumbar et al.,
2004). Our functional data using shRNA-mediated KD of the well-known
canonical EMT TFs confirm that, in contrast to the positive gene regulation mediated
by the key epithelial TFs, these EMT TFs mainly act as transcriptional repressors.
Our data indicate that the TGF-β/Smad2 axis rapidly activates and represses
a large number of genes associated with EMT that present Smad2 binding sites in
their open or close chromatin regions in good agreement with the well-known pro-EMT
role of TGF-β (Siegel and Massague,
2003) and the recent data showing that TGF-β-responsive cells
mediated cell invasion in skin SCCs (Paniz-Mondolfi
et al., 2015).
The remarkable and unexpected similarity in the chromatin and
transcriptional landscape between HF lineage and EMT cells suggests that many EMT
genes are primed in the HF lineages and facilitate the development of EMT in
HF-derived tumors. The EMT primed genes in HF lineages consist of TFs associated
with HF stemness and differentiation (Runx1, Nfatc1, Tcf7l1, Tbx1) (Chen et al., 2012; Horsley et al., 2008; Morris et al.,
2004; Nguyen et al., 2006; Tumbar et al., 2004), and a great number of
secreted molecules promoting TGF-β signaling (Ltbp1, 3) or inhibiting BMP
signaling (Grem1, Flstl1), ECM proteins (Col, Postn, Lox), leading to the autocrine
or paracrine formation of an EMT prone niche in the HF lineages that promotes EMT in
HF oncogene-targeted cells. The transcriptional and epigenetic priming of the cancer
cell of origin to undergo EMT does not exclude that part of the changes in gene
expression and epigenetic landscape can be also the consequence of regulatory
signals arising from the microenvironment. Future studies will be needed to identify
the importance and nature of extrinsic signals released in a paracrine or autocrine
manner by the tumor cells and their stromal cells that regulate the epigenetic
landscape of tumor cells and EMT.
In conclusion, our study demonstrates the functional importance of the
cancer cell of origin in regulating EMT and uncovers the molecular mechanisms by
which the cancer cell of origin promotes or restricts EMT in primary skin
tumors.

STAR★METHODS
Detailed methods are provided in the online version of this paper and
include the following:

KEY RESOURCES TABLE

CONTACT FOR REAGENT AND RESOURCE SHARING
Further information and requests for reagents may be directed to, and
will be fulfilled by, the Lead Contact, Pr. Cédric Blanpain
(Cedric.blanpain@ulb.ac.be).

EXPERIMENTAL MODEL AND SUBJECT DETAILS

Mice
Rosa26-YFP(Srinivasetal.,
2001), K14CreER (Vasioukhin
etal., 1999), Lgr5CreER (Barker et
al., 2007), KRasLSL-G12D (Tuveson et al., 2004) and p53fl/fl
(Jonkers et al., 2001) mice have
been imported from the NCI mouse repository and the Jackson Laboratories.
Wim Declercq (Ghent University, Belgium) generated the
Rosa26-ΔNp63-IRES-GFP. NOD/SCID/Il2Rγ null mice were
purchased from Charles River.
All mice used in this study were composed of males and females with
mixed genetic background. Mouse colonies were maintained in a certified
animal facility in accordance with the European guidelines and with approved
ethical protocol (#483N).

Primary cell culture
FACS isolated tumor YFP+EpCam+ or Epcam−
cells were plated on γ-irradiated 3T3 feeder cells in 6-well plates.
Cells were cultured in MEM medium supplemented with 10% FBS, 0.4
mg/ml hydrocortisone, 10 ng/ml EGF, 2×10−9M T3,
1% penicillin/streptomycin, 2mM L-glutamine. The feeders were
removed using PBS/EDTA (1mM). Cells are incubated at 37°C with
20% O2, and 5% CO2.

METHOD DETAILS
No randomization or blinding was performed in this study. Sample-size
and statistical methods are indicated in the quantification and statistical
analysis paragraph.

KRasG12D p53fl/fl induced skin tumors
Tamoxifen (TAM) was diluted at 25 mg/ml in sunflower oil (Sigma).
2,5 mg TAM was administered intraperitoneally (IP) to
K14CreER/KRasLSL-G12D/p53fl/fl/Rosa-YFP+/+
mice while
Lgr5CreER/KRasLSL-G12D/p53fl/fl/Rosa-YFP+/+
and
Lgr5CreER/KRasLSL-G12D/p53fl/fl/Rosa26-ΔNp63-IRES-GFP
mice were treated with 4 daily IP doses of 2,5 mg of TAM at P28 as
previously described (Lapouge et al.,
2012; Lapouge et al.,
2011). 1 mg of TAM was given to reduce the number of tumors in
Lgr5CreER model.

Monitoring of tumor growth
Tumor appearance and size were detected by daily observation and
palpation. Mice were euthanized when tumor size reached 1cm3 or
when mice presented signs of distress. Skin tumors were measured using a
precision calliper allowing to discriminate size modifications >
0,1mm. Tumor volumes were measured the first day of appearance of the tumor
and then, every week until the death of the animal with the formula V
= π × [d2 × D] /
6, where d is the minor tumor axis and D is the major tumor axis.

Antibodies
The following primary antibodies were used: anti-K14 (polyclonal
rabbit, 1:1000, Thermo Fisher Scientific), anti-GFP (chicken, 1:1000,
Abcam), anti-SOX2 (rabbit, 1:100, Abcam), anti-Vimentin (Rabbit, 1:200,
Abcam), Anti-ECadherin (rat, clone DECMA-1, 1:1000, eBioscience), anti-p63
(polyclonal rabbit, 1:200, Santa Cruz), Anti-Zeb1 (polyclonal rabbit, 1:300,
Bethyl), Anti-Zeb2 (polyclonal rabbit, 1:200, Sigma), anti-EpCam (rabbit
polyclonal, 1:200, Abcam). The following secondary antibodies were used:
anti-rabbit, anti-rat, anti-chicken, conjugated to AlexaFluor488 (1:400,
Molecular Probes), to rhodamine Red-X or to Cy5 (1:400, Jackson
ImmunoResearch).

Histology and immunostaining
For the staining on frozen sections, tissues were pre-fixed in
4% paraformaldehyde during 2 hr at room temperature, then washed in
PBS, incubated overnight in 30% sucrose at 4°C, and embedded
in OCT (Tissue Tek) for cryopreservation. Samples were sectioned at 5 mm
sections using CM3050S cryostat (Leica Microsystems GmbH). Nonspecific
antibody binding was blocked with 5% horse serum, 1% BSA,
and 0.2% Triton X-100 during 1 hr. Primary antibodies were incubated
overnight at 4°C in blocking buffer. Sections were rinsed in PBS and
incubated with secondary antibodies during 1 hr at room temperature. Nuclei
were stained with Hoechst (4 mM). Slides were mounted using Glycergel (Dako)
supplemented with 2.5% DABCO (Sigma-Aldrich).
For the staining on paraffin sections (for Sox2, Zeb1, Zeb2 and YFP
antibodies), 4 mm paraffin sections were deparaffinized and rehydrated.
Antigen unmasking was performed in citrate buffer (pH 6) at 98°C
during 20 min using the PT module. Endogenous peroxydase was blocked using
3% H2O2 (Merck) in methanol (VWR) during 20
min at room temperature. Endogenous avidin and biotin were blocked using the
Endogenous Blocking kit (Invitrogen) during 20 min at room temperature.
Primary antibodies were incubated overnight at 4°C. Anti-rabbit
biotinylated secondary antibodies were used, as well as Standard ABC kit,
and ImmPACT DAB (Vector Laboratories) for the detection of HRP activity.
Slides were mounted using SafeMount (Labonord).

Image acquisition
Imaging was performed on a Zeiss Axio Imager M1 (Thornwood)
fluorescence microscope with a Zeiss Axiocam MR3 camera and a Zeiss Axiocam
MRC5 camera for bright-field microscopy using Axiovision release 4.6
software. Brightness, contrast, and picture size were adjusted using
Photoshop CS6 (Adobe).

FACS Isolation of TECs and TMCs
Tumors were dissected, minced and digested in collagenase I (Sigma)
during 2 hr at 37°C on a rocking plate. Collagenase I activity was
blocked by the addition of EDTA (5 mM) and then the cells were rinsed in PBS
supplemented with 2% FBS. Before the staining, cells were blocked
during 20 min at room temperature in PBS supplemented with 30% FBS.
Cell suspensions were filtered through a 70 mm cell strainers (BD) then
through a 40 mm cell strainer to ensure the elimination of cell debris and
clumps of cells. Immunostaining was performed using PE-conjugated anti-CD45
(clone 30F11,1:100, eBioscience), PE-conjugated anti CD31 (clone MEC13.3;
1:100, BD PharMingen), and APC-Cy7-conjugated anti-Epcam (clone G8.8; 1:100,
Biolegend), during 30 min at 4°C on a rocking plate. Living tumor
cells were selected by forward scatter, side scatter, doublets
discrimination and by Hoechst dye exclusion. EpCam+ and EpCam- tumor
cells were selected based on the expression of YFP and the exclusion of
CD45, CD31 (Lin-). Fluorescence-activated cell sorting analysis was
performed using FACSAria and FACSDiva software (BD Bioscience). Sorted cells
were collected either in culture medium for in vivo transplantation
experiments or into lysis buffer for RNA extraction.

Tumor transplantation assays
The different FACS isolated populations of tumor cells (YFP+
Epcam+ from K14CreER and YFP+ EpCam+ from Lgr5CreER)
were collected in 4°C medium. Cells at different dilutions (1000 /
100 / 10 cells) were resuspended in 50 ml of Matrigel (50 ml, E1270, 970
mg/ml; Sigma) and injected subcutaneously to NOD/SCID/Il2Rg null mice
(Charles River, France). Triplicate injections per mouse were performed.
Secondary tumors were detected by palpation every week and their size
monitored until tumor reached 1cm3 or when mice presented signs
of distress, and the mice were sacrificed.

In vitro TGF-β treatment
FACS isolated tumor YFP+EpCam+ cells were plated on
g-irradiated 3T3 feeder cells in 6-well plates. For stimulation experiments,
media were supplemented with recombinant mouse TGF-β1 and
TGF-β2 (10ng/ml) (catalog number #76 66-MB and
#7346-B2 respectively, R&D Systems; resuspended with 4mM
HCl, 0,1% BSA).

Virus production, infection and selection
Stable knockdown cell lines were generated using lentiviral
pLKO/PuroR vectors (Sigma) after puromycin selection (2,5 mg/mL for TECs and
10 mg/mL for TMCs). Knockdown was confirmed by qRT-PCR. Three different
shRNA were used at the same time to target the same gene. The list of all
the shRNA used is listed in Table S2.
For virus production, 5×106 HEK293T cells were
seeded into 10 cm dishes and transfected with the vector of interest and
appropriate packaging plasmids psPax2 and pMD2.G (#12260 and
#12259 respectively, Addgene). Medium was changed 24 hr later and
next, supernatants were collected at 48 hr, and passed through a 0.45
μm filter. TECs of TMCs were plated in 6-well plate cells and
incubated with 40 μl/ml viruses when they reach 50% of
confluence, in the presence of polybrene (5 μg/ml). Medium was
changed 24 hr later and cells were selected with Puromycin for at least1
week.

RNA and DNA extraction, real-time PCR
RNA extraction from FACS isolated cells was performed using the
RNeasy micro kit (QIAGEN) according to the manufacturer’s
recommendations with DNAase treatment. After nanodrop RNA quantification,
the first strand cDNA was synthesized, using Superscript II (Invitrogen) and
random hexamers (Roche) in 50μl final volume. Control of genomic
contaminations was measured for each sample by performing the same procedure
with or without reverse transcriptase. Quantitative PCR assays were
performed using 1 ng of cDNA as template, SYBRGreen mix (Applied Bioscience)
and an Light Cycler ® 96 (Roche) real-time PCR system. TBP
housekeeping gene was used for normalization. Primers were designed using
Roche Universal ProbeLibrary Assay Design Center (https://lifescience.roche.com/webapp/wcs/stores/servlet/CategoryDisplay?tab=Assay+Design+Center&identifier=Universal+Probe+Library&langId=−1)
and are presented in Table
S2. Quantitative PCR Analysis was performed using Light Cycler
® 96 (Roche) and the DDCt method with TBP as a reference.

Microarray analysis
Total RNA was analyzed using Mouse whole genome 430 2.0 array from
Affymetrix at the AROS Applied Biotechnology A/S microarray facility (Aros,
Denmark) and Mouse whole genome 430 PM at IRB Functional Genomics Core
(Barcelona, Spain). 4 different biological replicates of FACS isolated
YFP+ EpCam+ cells from K14 tumors, 4 biological replicates
of FACS isolated YFP+ EpCam+ cells from Lrg5 mixed SCCs, 4
biological replicates of FACS isolated YFP+ EpCam− cells,
from Lgr5 mixed SCCs, 2 biological replicates of FACS isolated YFP+
EpCam+ cells from Lgr5 ΔNp63 and 2 biological replicates of
FACS isolated YFP+ EpCam- cells from Lgr5 ΔNp63 were
analyzed. Cells from Interfollicular Epidermis and infundibulum (Lgr5-GFP
negative, α6HCD34-) and Lgr5-GFP+ Hair Follicle cells have
been FACS sorted before KRasG12D p53cKO recombinaison and 2 biological
replicates of each case were analyzed.
All the results were normalized using the fRMA normalization using
R-bioconductor package fRMA with standard parameters. Cross platform
normalization was further performed to eliminate the batch effect using
ComBat unsupervised clustering using the Surrogate Variable Analysis, and
heatmap was generated using gplots (http://CRAN.R-project.org/package=gplots), all in
R-bioconductor.

ATAC-seq
Assay for transposase accessible chromatin (ATAC) followed by
sequencing was performed as following: 100000 sorted cells were collected in
1mL of PBS+3%FBS at 4°C. Cells were centrifuged,
then cell pellets were resuspended in 50 μL of lysis buffer (Tris
HCl 10mM, NaCl 10mM, MgCl2 3mM, Igepal 0,1%) and centrifuged (500 g)
for 25 min at 4°C. Supernatant was discarded and nuclei were
resuspended in 50 mL of reaction buffer (Tn5 transposase 2,5 μL, TD
buffer 22,5 μL and 25 μL H2O – Nextera DNA sample
preparation kit, Illumina). The reaction was performed for 30 min at
37°C and then blocked by addition of 5uL of clean up buffer (NaCl
900mM, EDTA 300mM). DNA was purified using the MinElute purification kit
(QIAGEN).

Library preparation and sequencing
DNA libraries were PCR amplified (Nextera DNA Sample Preparation
Kit, Illumina), and size selected for 200 to 800 bp (BluePippin, Sage
Sciences), following manufacturers’ protocols.

Alignment and Peak calling
More than 50 000 000 reads were mapped to mouse genomic DNA in each
condition. ATAC-seq reads (single-end or paired-end) were aligned to mouse
genome (NCBI37/mm10) using Bowtie2 (version 2.2.3) (Langmead and Salzberg, 2012) using option of
“–local” for single-end and
“–local–very–sensitive–local–dovetail–dovetail
-X 1000” for paired-end. Mitochondrial reads were excluded from
downstream analysis and duplicate reads were removed by Picard tools
(http://broadinstitute.github.io/picard/). Alignment data
tracks were visualized by Integrative Genomics Viewer (IGV) (Robinson et al., 2011).
Peak calling was performed on each individual sample by HOMER (Heinz et al., 2010) with parameters
setting of “-L 0 -C 3 -size 1000 -minDist 1000 -tbp 3 -o
auto.” Peaks from different ATAC-seq samples were merged for
downstream analysis.

Differential Peak analysis
Pairwise comparisons of ATAC peaks between two conditions were
performed by R package DESeq2 (Love et al.,
2014), with reads count of each peak calculated by HTSeq-count
(Anders et al., 2015).
Significance is defined as adjusted p value smaller than 0.001 and fold
change more than 3. Peaks were assigned to the nearest Refseq annotated
genes with 500kb range.

Motif analysis
De novo motif search was performed using program
offindMotifsGenome.pl in the HOMER package (Heinz et al., 2010) with parameters setting of “-size
−250,250 -S 15 -len 6,8,10,12,16.” Incidences of specific
motif was examined by the program of annotate-Peaks.pl in the HOMER package
with size parameter “-size 500.”

GSEA analysis
GSEA analysis was performed using ranked fold change values
(Epcam− over Epcam+) of ATAC peaks for the displayed
dataset. For the upregulated or downregulated genes, the highest fold change
peak was selected to represent the gene and enrichment score was calculated
following GSEA documentation.

RNA sequencing
RNA quality was checked by Bioanalyzer (Agilent). For RNA extracted
from TECs or TMCs, indexed cDNA libraries were obtained using the Ovation
Solo RNA-Seq Systems (NuGen) following manufacturer recommendations. The
multiplexed libraries (11 pM) were loaded and sequences were produced using
a HiSeq PE Cluster Kit v4 and TruSeq SBS Kit v3-HS (250 cycles) on a HiSeq
1500 (Illumina). Approximately 8 million of paired-end reads per sample were
mapped against the mouse reference genome (GRCm38.p4/mm10) using STAR
software to generate read alignments for each sample. Annotations
Mus_musculus.GRCm38.84.gtf were obtained from ftp.Ensembl.org. After
transcripts assembling, gene level counts were obtained using HTSeq and
normalized to 20 millions of aligned reads. Fold of changes (FC) were
computed on these values between the conditions.

QUANTIFICATION AND STATISTICAL ANALYSIS

Estimation of tumor propagating cell frequency
The frequency of tumor propagating cells was calculated using the
extreme limiting dilution analysis (ELDA) online software as previously
described (http://bioinf.wehi.edu.au/software/elda/) (Hu and Smyth, 2009). The statistical p value was
obtained using a Chi-square test.

Statistics
Statistical and graphical data analyses were performed using Prism 5
(Graphpad) software. All experiments shown were replicated at least twice.
All data in histograms represent mean ± SEM. Data of tumor
propagating cells frequencies (Figure
2B) represent percentage with 95% confidence interval.
Data were tested for normality using either D’Agostino and Pearson
omnibus normality test or Kolmogorov-Smrinov test (with
Dallal-Wilkinson-Lilliefor P value). Statistical significance was calculated
by Mann-Whitney test when sample size was small (Figures 1C, 1H, 2G, 5H, 7D, S1C, S1D, and S1G),
by Wilcoxon matched paired signed ranktest for non-parametric paired data
(Figures S1J and
S4Q–S4T), by Fisher’s exact test for analysis of
proportions (Figure 1I) and by Log-Rank
test for ranked observations (Figure S1E). Chi square test
for analysis of proportions when n is large (Figure 2B) using the Graphpad Prism 6 software, considering p
< 0.05 as statistically significant. All tests are two-sided.

DATA AND SOFTWARE AVAILABILITY
The accession number for the microarray data, RNA sequencing, and ATAC
sequencing reported in this paper is GEO: GSE88989.

ADDITIONAL RESOURCES
Annotations of Mus_musculus.GRCm38.84.gtf were obtained from ftp.Ensembl.org.

Supplementary Material
Supplemental Information