UVB induces cutaneous squamous cell carcinoma progression by de novo ID4 methylation via methylation regulating enzymes.

BACKGROUND: Little is known about whether UVB can directly influence epigenetic regulatory pathways to induce cutaneous squamous cell carcinoma (CSCC). This study aimed to identify epigenetic-regulated signalling pathways through global methylation and gene expression profiling and to elucidate their function in CSCC development.
METHODS: Global DNA methylation profiling by reduced representation bisulfite sequencing (RRBS) and genome-wide gene expression analysis by RNA sequencing (RNA-seq) in eight pairs of matched CSCC and adjacent normal skin tissues were used to investigate the potential candidate gene(s). Clinical samples, animal models, cell lines, and UVB irradiation were applied to validate the mechanism and function of the genes of interest.
FINDINGS: We identified the downregulation of the TGF-β/BMP-SMAD-ID4 signalling pathway in CSCC and increased methylation of inhibitor of DNA binding/differentiation 4 (ID4). In normal human and mouse skin tissues and cutaneous cell lines, UVB exposure induced ID4 DNA methylation, upregulated DNMT1 and downregulated ten-eleven translocation (TETs). Similarly, we detected the upregulation of DNMT1 and downregulation of TETs accompanying ID4 DNA methylation in CSCC tissues. Silencing of DNMT1 and overexpression of TET1 and TET2 in A431 and Colo16 cells led to increased ID4 expression. Finally, we showed that overexpression of ID4 reduced cell proliferation, migration, and invasion, and increased apoptosis in CSCC cell lines and reduced tumourigenesis in mouse models.
INTERPRETATION: The results indicate that ID4 is downregulated by UVB irradiation via DNA methylation. ID4 acts as a tumour suppressor gene in CSCC development. FUNDING: CAMS Innovation Fund for Medical Sciences (CIFMS) (2016-I2M-3-021, 2017-I2M-1-017), the Natural Science Foundation of Jiangsu Province (BK20191136), and the Fundamental Research Funds for the Central Universities (3332019104).

Evidence before this study

The literature indexed in PubMed already shows that cutaneous squamous cell carcinoma (CSCC) arises from the malignant transformation of keratinocytes in epidermis and skin adnexa. DNA methylation plays a crucial role in gene expression regulation, X chromosome inactivation, and tumourigenesis. Deregulation of ID proteins has been implicated in tumorigenesis. Epigenetic silencing of ID4 is detected in many types of cancer, including gastric, colorectal, breast, and lung cancers, where silencing of ID4 leads to uncontrolled growth and differentiation of cancer cells via the TGF-β/BMP/SMAD pathway.

Added value of this study

The results indicate that ID4 is downregulated by UVB irradiation via DNA methylation. ID4 acts as a tumour suppressor gene in CSCC development.

Implications for all the available evidence

Since abnormal methylation patterns are involved in the pathogenesis of CSCC, ID4 is a potential novel target for CSCC diagnosis and treatment.

Alt-text: Unlabelled box

Introduction

Cutaneous squamous cell carcinoma (CSCC) arises from the malignant transformation of keratinocytes in epidermis and skin adnexa. It is the second most common type of non-melanoma skin cancer in the United States, accounting for 20% of all skin cancers and being responsible for up to 9000 deaths each year [1,2]. Exposure to ultraviolet B rays (UVB), fair skin, and immunosuppression are the main risk factors for developing CSCC [1,2]. The incidence of CSCC has been increasing annually for the past three decades and is expected to continue to increase due to the aging population [1,2]. Understanding CSCC tumourigenesis and progression is critical for novel prevention and treatment strategies.

DNA methylation plays a crucial role in gene expression regulation, X chromosome inactivation, and tumourigenesis [3], including in CSCC [4]. In normal cells, normal methylation patterns ensure the proper regulation of gene transcription [5]. It is well accepted that aberrant DNA methylation occurs before malignant transformation [6,7]. Generalized hypomethylation and localized hypermethylation may lead to disrupted metabolism, mitotic recombination, copy number variations, chromosomal rearrangement, silencing of anti-oncogenes, and overexpression of proto-oncogenes [3]. DNA methylation changes have been extensively studied in many types of cancer [8], [9], [10], and DNA methyltransferases (DNMTs) play crucial roles [3]. It has been suggested that DNA methylation can be used as a cancer biomarker and as potential therapeutic targets in breast, gastric, liver, and lung cancers [11], [12], [13], [14]. Only a handful of studies investigated gene-specific DNA methylation changes in CSCC [15], [16], [17], [18], [19]. Currently, there are no studies on global methylation profiling in CSCC, correlating DNA methylation with gene expression, or mechanistic studies on candidate methylated genes and their association with clinical characteristics of CSCC.

Inhibitor of DNA binding/differentiation 4 (ID4) is one of the four ID proteins identified in mammalian cells and acts as a dominant-negative regulator of helix-loop-helix transcription factors [20]. It is a downstream mediator of the TGF-β/BMP/SMAD signalling pathway and regulates the growth and differentiation of embryonic tissues [21]. Deregulation of ID proteins has been implicated in tumourigenesis. Epigenetic silencing of ID4 is detected in many types of cancer, including gastric, colorectal, breast, and lung cancers, where silencing of ID4 leads to uncontrolled growth and differentiation of cancer cells via the TGF-β/BMP/SMAD pathway [20,[22], [23], [24]]. In estrogens receptor-positive breast cancer, gastric cancer, colorectal cancer, and leukaemia, ID4 is a tumour suppressor gene and is often inactivated by hypermethylation, leading to cancer growth [20,22,25,26]. On the other hand, ID4 is an oncogene in high-grade serous ovarian and liver cancers [27,28]. Currently, there are no studies investigating the role of ID4 in CSCC.

In this study, we first performed global methylation profiling by reduced representation bisulfite sequencing (RRBS) as well as transcription profiling by RNA-seq in eight pairs of CSCC and adjacent normal tissues to investigate the role of epigenetic changes in CSCC development and progression. Subsequently, we focused on DNA methylation regulation of one candidate gene, ID4, and its function in CSCC tumourigenesis.

Materials and methods

Ethics statement

This study was approved by the Medical Ethics Committee of Institute of Dermatology of the Chinese Academy of Medical Science and Peking Union Medical College (Institute of Dermatology, CAMS & PUMC). The study was performed according to the principles of the Declaration of Helsinki. All individuals participating in the study signed informed written consents. The animal experiments were approved by the Animal Ethics Committee of Institute of Dermatology, CAMS & PUMC and followed the Institutional Animal Care and Use guidelines.

Clinical skin specimens and cell lines

Patients admitted at Institute of Dermatology, CAMS & PUMC in 2015–2018 with the diagnosis of CSCC were approached for the study, without restrictions on age, sex, and ethnicity. A total of 202 pairs of fresh CSCC and adjacent normal tissues were collected and stored in liquid nitrogen. In addition, normal skin samples from distant non-sun exposed areas were collected from 60 CSCC patients, and skin samples from sun-exposed face areas were collected from 60 healthy patients (for neval removal) and stored in liquid nitrogen.

The human CSCC cell lines A431 and Colo16 and the immortalized keratinocytes HaCaT were purchased from the China Centre for Type Culture Collection. The cells were grown in DMEM media supplemented with 10% bovine serum at 37 °C with 5% CO2.

RRBS library construction and DMCs identification

About 1 µg of purified genomic DNA was digested with the methylation-insensitive enzyme MspI at 37 °C overnight. Subsequently, fragment ends were filled in and ligated with methylated paired-end Illumina indexing adapters. Adapters-ligated fragments were purified from excess free adapters before being subjected to bisulfite conversion using the EZ DNA Methylation Gold Kit (D5005, Zymo, USA). The bisulfite-converted DNA was PCR-amplified, analysed for fragment size using a bioanalyser, and quantified using both Qubit fluorometer with Quant-iT dsDNA HS Assay kit (Invitrogen) and Real-time quantity PCR (qPCR) with the KAPA library quantification kit (Cat. KK4824, Kapa Biosystems). An equal amount of amplified library from each sample was pooled and sequenced from both ends using an Illumina Hiseq 2500 Platform (Illumina, San Diego, CA, USA).

Quality control for the initial sequencing data was performed using the SOAPnuke software. Subsequent analysis was performed using BSMAP (version 2.74) and the February 2009 (GrCh37/hg19) build of the human genome. The percentage of methylation for each cytosine-phosphate-guanine (CpG) site was calculated. All CpG sites with at least 5 × sequencing coverage were used for the following analyses. Using the Student's t-test, differentially methylated CpG (DMC) sites between CSCC and adjacent normal tissues were identified (p < 0.05), and marked on the annotated human genome map. Finally, Gene Ontology (GO) and KEGG pathway enrichment analyses were performed using DAVID.

RNA-seq library construction and differentially expressed gene (DEG) identification

Total RNA was extracted from CSCC and adjacent normal tissues using Trizol (Invitrogen, Carlsbad, CA, USA). cDNA libraries were constructed using the NEBNext Ultra RNA Library Prep Kit for Illumina (New England Biolabs, USA). Library fragment size was measured using Bioanalyser. Libraries were quantified using both Qubit fluorometer with the Quant-iT dsDNA HS Assay kit (Invitrogen) and qPCR with the KAPA library quantification kit (Cat. KK4824, Kapa Biosystems). An equal amount of library cDNA from each sample was pooled and sequenced from both ends using the Illumina Hiseq 2500 Platform (Illumina, San Diego, CA, USA).

Initial sequencing data quality control was performed using the SOAPnuke software. All subsequent analysis was performed using Tophat (v2.0.13), and the reads were mapped to the February 2009 (GrCh37/hg19) build of the human genome. Transcripts were constructed using Cufflinks (v2.2.1) and merged to the annotated human genome map using Cuff merge in order to obtain the highest quality transcripts. RPKM analysis was used to determine gene expression, and the edgeR package (v.3.4.2) in the R software was used to identify DEGs between CSCC and adjacent normal tissues (gene expression fold change >2, p < 0.05). Finally, GO and KEGG pathway enrichment analyses were performed using DAVID.

DNA methylation analysis by MassArray

MethPrimer was used to design PCR primers to amplify bisulfite-converted BMP4, BMP8B, BMPR2, SMAD9, and ID4 genes (Table 1). Each amplicon was located in the promoter region and covered the majority of the CpG sites. Genomic DNA from CSCC and adjacent normal tissues was extracted using the QIAamp DNA Mini Kit (Qiagen). Bisulfite conversion was performed using the EZ DNA Methylation-Gold Kit (ZYMO, USA). PCR amplification was performed using the following program: 94 °C for 4 min; 45 cycles of 94 °C for 20 s, 56 °C for 30 s, and 72 °C for 1 min; 72 °C for 3 min; 4 °C for hold. After PCR amplification, free nucleotides were removed by alkaline phosphatase. In vitro transcription was performed using T7 RNA polymerase (Agena, USA) followed by base-specific enzymatic reaction. A nanodispenser was used to transfer the reaction mixture to a 384 SpectroCHIP for mass spectrometry analysis.

Table 1

Sequences of primers used in MassArray.

Species	Gene	Primer	Sequence (5 '→3 ')	Length	Tm (°C)
Human	BMP4	tag-FW	AGGAAGAGAGGGGTTTTTATTTTTAGAAAGGGAGG	10+25	60
		T7-RV	CAGTAATACGACTCACTATAGGGAGAAGGCTAAACTCCTAAACCCCCTCTACCTAT	31+25
	BMP8B	tag-FW	AGGAAGAGAGGGGTGTTTTAGAAGGGTTTTAGAGT	10+25	59
		T7-RV	CAGTAATACGACTCACTATAGGGAGAAGGCTAATCCCTACCCTACCCTACCC	31+21
	BMPR2	tag-FW	AGGAAGAGAGGGTTTTGTTTGTTTTTAGTTTGTGG	10+25	58
		T7-RV	CAGTAATACGACTCACTATAGGGAGAAGGCTTCAAAAAAATAATCTTTCCAATTCC	31+25
	SMAD9	tag-FW	AGGAAGAGAGAGAAAAGGTATTTGTTGTAGGGGTG	10+25	59
		T7-RV	CAGTAATACGACTCACTATAGGGAGAAGGCTTAACAATAAAATCCACATCCAACCT	31+25
	ID4	tag-FW	AGGAAGAGAGGGTTTGGAGTGGTTAGTTAATTAGG	10+25	58
		T7-RV	CAGTAATACGACTCACTATAGGGAGAAGGCTAAAAAACTACACATTCCATTCCATC	31+25
Mouse	ID4	tag-FW	AGGAAGAGAGGAGTGATTAGTTAATTAGGAGGATAGTG	10+28	56
		T7-RV	CAGTAATACGACTCACTATAGGGAGAAGGCTAAAAACCTAAAAACTAAACTCCCCC	31+25

Gene expression analysis by qPCR

The Primer 5.0 software was used to design mRNA-specific amplification primers for each target gene (Table 2). Total RNA was extracted from fresh skin tissues using the RNeasy Mini Kit (Qiagen, Germany), reverse transcribed into cDNA using the PrimeScriptTM RT mastermix (TaKaRa, Dalian, China), and PCR-amplified using the AceQ qPCR SYBR Green mastermix (Vazyme, Nanjing, China).

Table 2

Sequences of primers used in qPCR.

Species	Primers	Forward Sequence (5’ - 3’)	Reverse sequence (5’ - 3’)
Human	BMP4	F: ATGATTCCTGGTAACCGAATGC	R: CCCCGTCTCAGGTATCAAACT
	BMP8B	F: GGAGCCCCATTGGAAGGAG	R: CTCGGAGCGTCTGAAGATCC
	BMPR2	F: CGGCTGCTTCGCAGAATCA	R: TCTTGGGGATCTCCAATGTGAG
	SMAD9	F: CTAGGCTGGAAGCAAGGAGAT	R: GGGGAATCGTGACGCATTT
	ID4	F: TCCCGCCCAACAAGAAAGTC	R: CCAGGATGTAGTCGATAACGTG
	DNMT1	F: CCTAGCCCCAGGATTACAAGG	R: ACTCATCCGATTTGGCTCTTTC
	DNMT3A	F: CCGATGCTGGGGACAAGAAT	R: CCCGTCATCCACCAAGACAC
	DNMT3B	F: AGGGAAGACTCGATCCTCGTC	R: GTGTGTAGCTTAGCAGACTGG
	TET1	F: CATCAGTCAAGACTTTAAGCCCT	R: CGGGTGGTTTAGGTTCTGTTT
	TET2	F: CCAGACTATGTGCCTCAGAAATCC	R: GAAACGCAGGTAAGTGGGCTC
	TET3	F: CCCACGGTCGCCTCTATCC	R: CTGCGACATCCTTCTCAT
	β-actin	F: CCATCGTCCACCGCAAAT	R: GCTGTCACCTTCACCGTTCC
Mouse	ID4	F: CAGTGCGATATGAACGACTGC	R: GACTTTCTTGTTGGGCGGGAT
	DNMT1	F: ATCCTGTGAAAGAGAACCCTGT	R: CCGATGCGATAGGGCTCTG
	DNMT3A	F: CTGTCAGTCTGTCAACCTCAC	R: GTGGAAACCACCGAGAACAC
	DNMT3B	F: CGTTAATGGGAACTTCAGTGACC	R: CTGCGTGTAATTCAGAAGGCT
	TET1	F: ACACAGTGGTGCTAATGCAG	R: AGCATGAACGGGAGAATCGG
	TET2	F: AGAGAAGACAATCGAGAAGTCGG	R: CCTTCCGTACTCCCAAACTCAT
	TET3	F: TGCGATTGTGTCGAACAAATAGT	R: TCCATACCGATCCTCCATGAG
	β-actin	F: GTCCCTCACCCTCCCAAAAG	R: GCTGCCTCAACACCTCAACCC

UVB irradiation

An SS-04B ultraviolet phototherapy instrument (Sigma, Shanghai, China) was used to deliver UVB exposure to male C3H/HeN mice (150 mJ/cm2) as well as to the A431, Colo16, and HaCaT cell lines (10 mJ/cm2). Control mice were not irradiated, while irradiation was conducted for 4 days on half of the skin of the test mice, the remaining half being shaded.

Immunohistochemistry analysis

Tissue sections were deparaffinized and subjected to antigen retrieval. The endogenous peroxidase was blocked using 3% hydrogen peroxide. The slides were first blocked with normal goat serum at room temperature for 30 min to minimize nonspecific staining, then incubated overnight with primary antibodies against DNMT1 (1:100, Absin Bioscience Inc., Shanghai, China), DNMT3A (1:50, CST, USA), DNMT3B (1:50, Abcam, UK), TET1 (1:100, Absin Bioscience Inc., Shanghai, China), TET2 (1:100, Absin Bioscience Inc., Shanghai, China), and TET3 (1:100, Abcam, UK). Subsequently, the slides were incubated with HRP-labelled goat anti-rabbit/mouse secondary antibody at 37 °C for 20 min, counter-stained with hematoxylin, dehydrated, and stabilized with mounting medium.

Transfection

An ID4 expression lentivirus was constructed by Genechem Biomart (Shanghai, China) and used to infect the A431, Colo16, and HaCaT cells. To determine the action of DNMT1, TET1, and TET2 on ID4 gene in CSCC, we transfected DNMT1-siRNA (sense: 5′-GGAUGAGUCCAUCAAGGAATT-3′, antisense: 5′-UUCCUUGAUGGACUCAUCCTT-3′, 100 pmol/well) and negative control-siRNA (sense: 5′-UUCUCCGAACGUGUCACGUTT-3′, antisense: 5′-ACGUGACACGUU CGGAGAATT-3′, 100 pmol/well) into A431 and Colo16 cells in 6-well plates using Lipofectamine 2000 (Invitrogen, CA, USA), according to the manufacturer's protocol. The TET1 and TET2 expression lentiviruses were constructed by Genechem (Shanghai, China) and used to infect the A431 and Colo16 cells. The expression of the target genes was verified by western blot, as below.

Western blot

Total protein was extracted from infected cell lines using the RIPA lysis buffer (Beyotime, Jiangsu, China) containing 1% protease inhibitor cocktail (Sigma, USA). Protein concentration was measured using the BCA protein assay kit (Beyotime, Jiangsu, China). About 50 µg of protein from each cell line was loaded onto 10% SDS-PAGE and afterward transferred to Immun-Blot PVDF membranes (BioRad, USA). The PVDF membranes were incubated with rabbit antihuman ID4 antibody (1:1000, ab49261, Abcam, UK), DNMT1 antibody (1:1000, #5032, CST, USA), TET1 antibody (1:1000, ab191698, Abcam, UK), TET2 antibody (1:1000, #18950, CST, USA), and anti-β-actin antibody (1:1000, #4970, CST, USA) as an internal control. Images were developed using the 20 × LumiGLO® Reagent and 20 × peroxide (CST, USA).

Cell proliferation assay

For all the following experiments, the cells in each group had a good growth state and were in the logarithmic growth phase when inoculated. The differences in cell activity, apoptosis, migration, and invasion were detected 24 h later. Cells were plated on 96-well plates at 1 × 104 cells/well. After 24 h of incubation, 10 μL of Cell Counting Kit-8 (CCK-8) reagent (Dojindo Laboratories, Kumamoto, Japan) were added to each well and incubated at 37 °C for 4 h. Absorption at 450 nm was measured using a microplate reader.

Cell cycle analysis

Cells were washed and resuspended in PBS. Cold alcohol was gradually added to cell suspensions to reach a final alcohol concentration of 70%. Cells were fixed for at least 4 h at 4 °C, then incubated with 100 μL RNase (50 μg/mL) at 37 °C for 30 min. The reaction was stopped by adding 400 mL PI (50 μg/mL). The cell suspensions were kept at 4 °C in the dark for 30 min and then subjected to flow cytometry analysis.

Apoptosis analysis

The cells were digested with 0.25% pancreatic enzymes in the absence of EDTA. After the reaction was terminated, the supernatants were discarded, and the cells were resuspended in PBS. The cells were stained with 5 μL of AnnexinV-APC followed by 5 μL of 7-AAD according to the AnnexinV-APC/7-AAD Apoptosis Staining/Detection Kit protocol (eBioscience, San Diego, CA, USA). Cell suspensions were left at room temperature in the dark before flow cytometry analysis.

Wound healing assay

Cells were plated on 6-well culture plates at 1 × 106 cells/well. After the cells were attached to the bottom of the wells, a straight line was scratched in the cell layer. The wells were washed three times with PBS, and serum-free medium was added to each well. Photographs were taken at 0 and 24 h after the initial scratch. Cell migration rates were calculated.

Transwell invasion assay

About 600 μL of 10% FBS containing DMEM medium were added to the bottom wells of a 24-well plate, and an 8-μm transwell insert was added to each well, filled with the Matrigel Basement Membrane Matrix (#356234, BD Biosciences, Franklin Lake, NJ, USA). In the upper chambers, 200 μL of cell suspension in serum-free media were added. The plates were incubated for 24 h, and the upper transwell inserts were removed. The bottom wells were washed with PBS, then fixed with 70% alcohol for 1 h, and stained with 0.5% crystal violet for 20 min at room temperature.

Tumourigenesis in nude mice

BALB/c nude mice were randomly assigned into seven groups (A431, A431-NC, A431-ID4, Colo16, Colo16-NC, Colo16-ID4, and HaCaT), each group having five mice. After cleaning with alcohol, 1 × 107 cells of each cell line were injected into the right dorsal skin of each mouse. Tumour growth was monitored each day. The mice were sacrificed after 30 days by excessive anethesia. The tumours were resected and measured.

Statistical analysis

All analyses were carried out using SPSS 20.0 (IBM, Armonk, NY, USA). Continuous data were tested for normal distribution using the Kolmogorov–Smirnov test, are presented as means ± standard deviations, and were analysed using ANOVA and the LSD post hoc test. Categorical variables are presented as n (%) and were analysed using the chi-square test or Fisher's exact test, as appropriate.

Results

Methylome analysis

We performed RRBS [29] on eight pairs of CSCC and adjacent normal tissues. An average of 3.4 million CpGs with at least 5 × sequencing coverage was used in downstream methylation data analysis (Table S1). Principal component analysis of all samples using all commonly covered CpG sites showed that CSCC tissue samples were separated from adjacent normal tissues (Fig. S1a). In general, both CpG islands and shores of CpG islands were hypermethylated in CSCCs, while repeated sequences, especially LTR and satellites, were hypomethylated in CSCCs (Fig. 1a).

Fig. 1

Aberrant methylation and gene expression in cutaneous squamous cell carcinoma. (a) Violin plots of methylation of genic regions (upstream 2 kb from transcription start site, 5-UTR, CDS, intron, 3-UTR and downstream 2 kb from transcription end site) and regulatory elements (shores of CpG island, CpG island and repeat) in CSCCs (cyan) and normal tissues (purple). White points represent the median methylation level. (b) Relative enrichment of All-DMCs, Hyper-DMCs, and Hypo-DMCs within genic regions and regulatory elements. The enrichment was evaluated by the ratio of All-DMCs (gray), Hyper-DMCs (purple), or Hypo-DMCs (blue) counts divided by CpGs counts within these genomic elements out of all DMCs counts divided by all CpGs counts in human genome. (c) Volcano plot of all gene expressions. Differentially expressed genes are shown in red. (d) The percentage of DMCs and associated DEGs pairs showed a positive correlation (purple) or negative correlation (blue). Genes were grouped into upstream or gene bodies by the location of CpG island or repeat with DMCs. (e) Scatter plot of differential methylation and differential gene expression between CSCCs and normal tissues. These genes were grouped into upstream or gene body by the location of CpG island or repeat with DMCs and Pearson correlation coefficient was calculated in each group. (f) KEGG pathway and disease enrichment analysis of DEGs, which were negatively correlated with DMCs in CpG islands, or positively correlated with DMCs in repeat, located in promoter or gene body. The x-axis shows the q-value from hypergeometric test adjusted by the multiple test adjustment. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

A total of 50,184 DMCs were generated, and 63.22% of them were hypermethylated in CSCC (Table S2). Using these DMCs, CSCC tissues were clustered together by unsupervised clustering analysis (Fig. S1b). DMCs were enriched at satellite regions, but depleted at simple repeat and low complexity regions; hypermethylated DMCs were more frequent in genic regions and CpG islands, while hypomethylated DMCs were more frequent in repetitive sequences (Fig. 1b). Our results suggested that hypermethylation of CpG islands and hypomethylation of repetitive sequences were characteristics of CSCC.

Transcriptome analysis

We performed a transcriptome analysis on the same eight pairs of CSCC and adjacent normal tissues using RNA-seq. We identified 1053 DEGs, and 88.3% of them were downregulated in CSCC tissues (Fig. 1c and Table S3).

Integrative analysis of methylome and transcriptome of CSCC

Next, we correlated the gene expression levels of DEGs and DNA methylation levels of DMCs. We found that over 73% of DEGs showed negative correlation between gene expression and CpG island DMCs in the promoter region (77.15%) or in the gene body region (73.54%), while over 57% of DEGs showed positive correlation between gene expression and repetitive sequence DMCs in the gene body region (73.54%) or in the upstream region (56.94%) (Fig. 1d and e) (Table S4). Our results suggest that aberrant DNA methylation was correlated with gene expression changes in CSCCs.

Identification of aberrantly methylated genes in CSCC

We performed the KEGG pathway and disease category analyses for DEGs, showing correlation with DMCs. The top ten enriched KEGG pathways and disease categories could be grouped into three classes: (i) known cancer-related pathways; (ii) mental disease or nervous system pathway; and (iii) endocrine and metabolism (Fig. 1f).

To narrow down candidate methylated genes in CSCC, we first identified 339 DEGs that contained DMCs in the promoter region and showed a negative correlation between gene expression and DNA methylation; 23 genes showed both DNA methylation and gene expression differences between CSCC and adjacent normal tissues. Of these 23 genes, two (ID4 and UCHL1) were not only associated with malignancy but could also be targeted by the demethylating agent azacitidine [30,31] (Fig. 2a) (Table S5). The expression of ID4 and UCHL1 was significantly downregulated in CSCC tissues and correlated with increased levels of promoter methylation (Fig. 2b and c).

Fig. 2

Selection of ID4 for functional validation. (a) Overlap of gene with 10% methylation difference in promoter (blue), DEGs with differential TPM over 1 (purple). and azacitidine target genes (green). The number of genes associated with the tumour is shown in the ellipse (yellow). (b) The methylation level and mRNA expression level of 12 azacitidine target genes. The p-value of Wilcoxon test between CSCCs (purple) and normal tissues (yellow) was shown at the top. (c) Methylation level of CpGs in promoter of ID4. The DMCs are labelled by ‘*’, and the transcription is shown by black rectangle and arrowhead. (d) ID4 in TGF-β/BMP-SMAD-ID signal pathway. (e) Methylation of BMP4, BMP8B, BMPR2, SMAD9, and ID4 in CSCCs and normal tissues and p-value of differential analysis of each gene was indicated on the right. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Downregulation of TGF-β/BMP-SMAD-ID signalling pathway in CSCC

Because ID4 is the downstream effector of the TGF-β/BMP-SMAD-ID signalling pathway (Fig. 2d), we performed KEGG cluster profiler analysis on DEGs from RNA-seq data and showed that expression of five genes (BMP4, BMP8B, BMPR2, SMAD9, and ID4) in the TGF-β/BMP-SMAD-ID signal pathway was significantly lower in CSCC tissues than in adjacent normal tissues (Table S3), although only DNA methylation of ID4 was significantly higher in CSCC tissues than in adjacent normal tissues, based on RRBS data (Fig. 2c).

Confirmation of downregulation of the TGF-β/BMP-SMAD-ID pathway in CSCC using independent expanded clinical samples

We performed promoter DNA methylation analysis by bisulfite sequencing and mRNA expression by qPCR of BMP4, BMP8B, BMPR2, and SMAD9 on 30 pairs of CSCC and adjacent normal tissues, and of ID4 on 202 pairs of CSCC and adjacent normal tissues. Our results confirmed the expression downregulation of BMP4, BMP8B, BMPR2, SMAD9, and ID4 in CSCC tissues, and promoter methylation of ID4 (Fig. 2e and Fig. S2a). Finally, we showed that ID4 methylation levels were negatively correlated with its gene expression in CSCC tissues (p < 0.001, r = −0.455, Spearman correlation) (Fig. S2b).

Promoter DNA methylation of ID4 was regulated by UVB

Because UVB is the major carcinogen for CSCC, we tested whether UVB played a role in ID4 promoter DNA methylation. We collected sun-exposed normal skin in the head and neck region (n = 60) from normal patients and distal non-exposed normal skin from CSCC patients (n = 60). All samples were used for DNMTs and TETs expression detection, and 20 samples in each group were used for ID4 methylation and expression detection. ID4 DNA methylation was significantly elevated in sun-exposed normal skin tissues compared with non-exposed normal skin tissues (p < 0.001, t-test), while ID4 gene expression was significantly lower in sun-exposed skin tissues (p < 0.05). Expression of DNMT1 (p < 0.01, t-test) was significantly higher in sun-exposed skin tissues than in non-exposed skin tissues, while expression of DNMT3A, TET1, TET2, and TET3 was significantly lower in sun-exposed skin tissues than in non-exposed skin tissues (p < 0.001 for DNMT3A, TET1 and TET3; p < 0.01 for TET2, t-test) (Fig. S3).

Next, we determined ID4 expression and methylation changes after UVB irradiation in mouse dorsal skin (Fig. 3a). In both autologous (exposed vs. non-exposed skins in the same mouse) and allogeneic mouse (exposed vs. non-exposed skins from different mice) models, after 4 days of UVB exposure, ID4 expression was significantly downregulated in exposed skins compared to non-exposed skins, and ID4 DNA methylation was significantly higher in exposed skins (p < 0.05, t-test) (Fig. 3b). Expression of DNMT1 and DNMT3A was significantly upregulated in exposed skins (p < 0.05 for each, t-test), while expression of DNMT3B, TET1, TET2, and TET3 was significantly downregulated in exposed skins (p < 0.01 for DNMT3B and TET3; p < 0.05 for TET1 and TET2, t-test) (Fig. 3c).

Fig. 3

UVB irradiation induces alteration of gene expression and methylation in mouse dorsal skin and CSCC cell lines. Analyses of ID4 DNA methylation (b), and ID4, DNMTs, and TETs mRNA expression of the shaved dorsal skin of male C3H/HeN mice after being exposed to UVB radiation (150 mJ/cm2) (a) in both autologous and allogeneic mouse models (c). Immunohistochemistry analysis of DNMTs and TETs proteins (d). Histopathological changes of epidermis were observed by HE staining (e) in UVB-exposed skin, non-UVB-exposed skin, and non-UVB-exposed mouse skin. Analyses of ID4 DNA methylation, and ID4, DNMTs and TETs mRNA expression of A431 (f), Colo16 (g), and HaCaT (h) cell lines after being exposed to UVB radiation (10 mJ/cm2). Data are shown as means ± SD (n = 3). The statistical significance was assessed by the Student's t-test (*P < 0.05, **P < 0.01, ***P < 0.001).

Immunohistochemistry analysis of mouse skin tissues further confirmed the mRNA analysis. DNMT1, DNMT3A, and DNMT3B proteins were not detected in non-exposed skin tissues but were detected in the basal layer of epidermis of exposed skin tissues. Diffuse staining of TET1, TET2, and TET3 proteins was detected in the epidermis of non-exposed skin tissues. The staining of these proteins was reduced in the basal layer of epidermis of exposed skin tissues (Fig. 3d). HE staining showed that UVB caused a certain degree of pathological changes in mouse skin, such as epidermal hyperplasia, disordered cell arrangement, and partial cell morphology abnormalities (Fig. 3e).

Finally, we determined ID4 expression and DNA methylation changes in three cell lines (A431, Colo16, and HaCaT) after UVB irradiation. We showed that in all three cell lines, ID4 DNA methylation was increased, and ID4 expression was decreased 4–8 or 12 h after UVB irradiation (Fig. 3f–h). At the same time, we observed the upregulation of DNMT1 and downregulation of DNMT3A, DNMT3B, TET1, TET2, and TET3 in all three cell lines except the upregulation of TET1 in HaCaT (Fig. 3f–h).

Correlation of DNMTs and TETs expression with ID4 expression and methylation in paired CSCC and adjacent normal tissues

RNA-seq data on eight pairs of CSCC and adjacent normal tissues showed that the expression of DNMT1 was higher in CSCC than in normal tissues, and expression of TET2 was lower in CSCC than in normal tissues.

Next, we determined the mRNA expression of DNMTs and TETs in 40 additional pairs of CSCC and adjacent normal tissues by qPCR. The expression of DNMT1 was significantly higher in CSCC than in adjacent normal tissues, while the expression of TET1, TET2, and TET3 was significantly lower in CSCC than in adjacent normal tissues (Fig. 4a).

Fig. 4

Altered expression of DNMTs and TETs in CSCCs compared with adjacent normal tissues. (a) The relative expression levels of DNMTs and TETs in CSCC and adjacent normal tissues. Data are shown as means ± SD. The statistical significance was assessed by the Student's t-test. (b) The Pearson's correlation between expression level (2−△Ct) or methylation level of ID4 and expression levels (2−△△Ct) of DNMT1, DNMT3A, DNMT3B, TET1, TET2, and TET3, which were measured by MSP and qPCR in the same patient sets.

We correlated DNMTs and TETs expression with ID4 expression and methylation in these 40 CSCC tissues. DNMT1 expression positively correlated with ID4 methylation level (p < 0.05, r = 0.327, Pearson correlation) and negatively correlated with ID4 mRNA expression (p < 0.01, r = −0.484, Pearson correlation); DNMT3A/3B expression was not correlated with ID4 methylation and mRNA expression; TET1 and TET2 expression was negatively correlated with ID4 methylation (p < 0.01, r = −0.410 for TET1; p < 0.01, r = −0.427 for TET2, Pearson correlation) and positively correlated with ID4 mRNA expression (p < 0.05, r = 0.378 for TET1; p < 0.05, r = 0.318 for TET2, Pearson correlation); TET3 expression was not correlated with ID4 methylation and mRNA expression (Fig. 4b). Our data suggested that hypermethylation of ID4 in CSCC tissues was likely due to the upregulation of DNA methyltransferase activity and downregulation of DNA demethylation enzymes.

Finally, to determine whether the ID4 protein was regulated by DNMT1, TET1, and TET2, we transfected DNMT1-siRNA, TET1, and TET2 expression lentivirus into A431 and Colo16 cells. As shown in Supplementary Fig. 4, the expression of DNMT1 was decreased by the DNMT1-siRNA, which resulted in an increased expression of ID4 after 72 h of transfection (p < 0.05, t-test) (Fig. S4a). The expression of TET1 (Fig. S4b) and TET2 (Fig. S4c) were increased by their respective lentivirus vectors, both resulting in an increased expression of ID4 (p < 0.05, t-test).

Overexpression of ID4 reduced proliferation, migration, and invasion, but increased apoptosis in CSCC

In order to study the function of ID4 in CSCC progression, we constructed stable ID4 overexpressing A431, Colo16, and HaCaT cell lines using the lentivirus method (Fig. S5). We showed that high ID4 expression significantly reduced cell proliferation in A431 (p < 0.01, t-test) and Colo16 cell lines (p < 0.001, t-test) (Fig. 5a). High ID4 expression significantly reduced the percentage of S-phase cells in all three cell lines (p < 0.05, t-test) (Fig. 5b). High ID4 expression significantly increased the apoptosis level in A431 and Colo16 cell lines (p < 0.05, t-test) (Fig. 5c). High ID4 expression significantly reduced cell migration in A431 (p < 0.05, t-test), Colo16 (p < 0.05, t-test), and HaCaT (p < 0.01, t-test) cell lines by the wound healing assay (Fig. 6a). High ID4 expression also significantly reduced cell invasion in A431 (p < 0.001) and Colo16 (p < 0.001, t-test) cell lines by the Transwell assay (Fig. 6b).

Fig. 5

Cell proliferation, cell cycle, and apoptosis analysis of ID4-overexpressing cell lines. (a) Cell proliferation was assessed using the CCK-8 assay in ID4-overexpressing cell lines. (b) Cell-cycle distribution of ID4-overexpressing cell lines. The percentage of cells in each phase of the cell cycle was indicated for each panel. (c) Apoptosis analysis of ID4 overexpressed cell lines. Cells were double-stained with annexin V and propidium iodide (PI). Data are shown as means ± SD (n = 3). The statistical significance was assessed using the Student's t-test (**P < 0.01, ***P < 0.001).

Fig. 6

Wound repair, transwell, and xenograft assays of ID4. (a) Wound repair assays were monitored at 24 h in ID4 overexpressed cell lines. Data are shown as means ± SD (n = 3). (b) Transwell assays were used to evaluate the involvement of ID4 for invasion in ID4-overexpressing A431 and Colo16 cell lines. (c) Overexpression of ID4 suppressed xenograft tumour growth. Data are shown as means ± SD of the tumour volumes (n = 5). The statistical significance was assessed using the Student's t-test (*P < 0.05, **P < 0.01, ***P < 0.001).

Finally, we determined the effects of UVB irradiation in these high ID4 expression cell lines. In the high ID4 expression A431 cell line, ID4 methylation (p < 0.05, t-test) was significantly increased after 8 h of UVB exposure, and ID4 expression (p < 0.05, t-test) was significantly decreased. In the high ID4 expression Colo16 cell line, ID4 methylation (p < 0.05, t-test) was already significantly increased after 4 h of UVB exposure, but ID4 expression (p < 0.05, t-test) was significantly decreased only after 8 h of UVB exposure. Similarly, in the high ID4 expression HaCaT cell line, ID4 methylation (p < 0.05, t-test) was already significantly increased after 4 h of UVB exposure, but ID4 expression (p < 0.05, t-test) was significantly decreased only after 8 h of UVB exposure (Fig. S6).

Overexpression of ID4 reduced tumourigenesis in mouse CSCC models

To further investigate the role of ID4 in vivo, we determined the effect of ID4 expression on CSCC tumour development in nude mice. High ID4 expression in the A431 cell line delayed tumour development from 4 to 8 days, and high ID4 expression in the Colo16 cell line delayed tumour development from 15 to 21 days after transplantation. Furthermore, tumour size was significantly smaller in the high ID4 expression A431 (p < 0.001, t-test) and Colo16 (p < 0.01, t-test) cell lines 30 days after transplantation (data not shown). Our results indicated that high ID4 expression significantly decreased tumourigenesis of A431 and Colo16 in the nude mouse CSCC models (Fig. 6c).

Discussion

In the present study, we performed an integrative analysis of global DNA methylation and gene expression on matched CSCC and adjacent normal skin tissues. We observed the downregulation of the TGF-β/BMP-SMAD-ID4 signalling pathway in CSCC and increased methylation of ID4. We showed that ID4 methylation could be directly induced by UVB exposure in normal skin tissues, and increased ID4 methylation is associated with the upregulation of DNMT1 and downregulation of TETs. Silencing of DNMT1 led to increased ID4 expression, and overexpression of TET1 and TET2 also led to increased ID4 expression. Functionally, we showed that overexpression of ID4 reduced cell proliferation, migration, invasion, and increased apoptosis in CSCC cell lines. Overexpression of ID4 suppressed tumourigenesis in mouse skin cancer models. Our data indicate that UVB exposure could directly downregulate ID4 expression via DNA methylation to initiate cutaneous tumourigenesis.

In the preliminary stage of the study, UVB exposure in normal human and mouse skin tissues and cutaneous cell lines induced ID4 DNA methylation, upregulated DNMT1 and downregulated TETs. Similarly, we detected the upregulation of DNMT1 and downregulation of TETs accompanying ID4 DNA methylation in CSCC tissues. Therefore, those preliminary results could suggest that UVB regulates the methylation and expression of the ID4 gene by affecting the expression of DNMT1, TET1, and TET2. To determine whether the ID4 protein was regulated by DNMT1, TET1, and TET2, we transfected DNMT1-siRNA, TET1, and TET2 expression lentivirus into A431, and Colo16 cells. The results showed that the expression of DNMT1 was decreased by the DNMT1-siRNA, which resulted in an increased expression of ID4 at 72 h after transfection. The expression of TET1 and TET2 was increased by their respective lentivirus vectors, both resulting in an increased expression of ID4.

Compared to whole-genome DNA methylation analysis, RRBS is a cost-effective high-throughput global DNA methylation analysis method at a single-nucleotide level [32,33]. RRBS has been extensively used to characterise cancer methylome [34], [35], [36], [37]. To our knowledge, no previous studies reported the global methylation profiles in CSCC. Our RRBS results confirmed that similar to other cancer types, CSCC was characterized by hypermethylation of CpG islands and hypomethylation of repetitive sequences in the genome. To further characterize DNA methylation changes in CSCC, we performed integrative analysis of RRBS and RNA-seq on the same tissue samples and showed that the majority of DNA methylation changes in CSCC were inversely correlated with gene expression changes. Of 23 candidate genes that were hypermethylated and downregulated in CSCC, two of them (ID4 and UCHL1) were also targets for the demethylating agent azacitidine. Because ID4 is the downstream effector of TGF-β/BMP-SMAD-ID signalling pathway and the expression of five genes (BMP4, BMP8B, BMPR2, SMAD9, and ID4) in this pathway were all downregulated in CSCC, we focused our subsequent analyses on ID4.

ID4 is a dominant-negative regulator of the basic helix-loop-helix family of transcription factors. Promoter hypermethylation-mediated inactivation of ID4 plays an important role in the development of solid tumours. In breast cancer, ID4 is hypermethylated and downregulated in ER+ tumours, and is associated with poor survival [20,38,39]. In myeloid malignancies, ID4 acts as a tumour suppressor gene, and its inactivation by DNA methylation is associated with poor treatment outcome and overall survival [40]. In non-Hodgkin lymphoma, ID4 methylation inhibits apoptosis and promotes cell proliferation [30]. In prostate cancer, ID4 downregulation is by DNA methylation [41]. Our result indicates that ID4 was downregulated by DNA methylation in CSCC tissues compared with adjacent normal tissues. Our data were consistent with a previous report on ID4 methylation associated with parakeratosis [42].

Our data indicate that UVB irradiation could directly lead to ID4 methylation through the upregulation of DNMT1 and downregulation of TET1 and TET2. Accordingly, a silencing of DNMT1, which is involved in methylation [3,43], led to an upregulation of ID4, while overexpression of TET1 and TET2, which are involved in demethylation [44], also led to an increased ID4 expression. Although UVB irradiation is the major risk factor for CSCC, little is known whether UVB irradiation could influence the epigenome during CSCC progression. An early study using a microarray method on human keratinocytes showed that UVB irradiation does not directly induce detectable changes in DNA methylation [45]. On the other hand, recent studies using sequencing-based methods in mice skin cancer models reported distinct DNA hypermethylation patterns in UVB-exposed epidermal skin and UVB-induced skin tumours [46,47]. In both human keratinocytes and mouse skin cancer models, we showed that UVB exposure could directly lead to ID4 methylation. We showed that UVB exposure could significantly increase the expression of DNMT1 and DNMT3A, but not DNMT3B. In human clinical samples, only DNMT1 expression was significantly increased in CSCC tissues compared to adjacent normal tissues. Interestingly, we consistently detected significant downregulation of TET1, TET2, and TET3 in human CSCC tissues as well as UVB-exposed human CSCC cell lines and UVB-exposed normal mouse skin. TET dioxygenase-mediated oxidation of 5mC represents an active DNA demethylation process [48]. Our data suggested that TETs might coordinate with DNMT1 to regulate DNA methylation in skin cells. Interestingly, Wang et al. [49] showed that UVB irradiation upregulated all three TETs in the HaCaT cell line. The discrepancy might be due to differences in experimental procedures. Our data suggested that TETs might play a dominant role in regulating DNA methylation in skin cells, but future studies are needed to decipher the precise underlying mechanism.

The present study suggests that ID4 is a tumour suppressor gene in CSCC carcinogenesis. ID4 overexpression reduced cell proliferation, induced apoptosis, reduced migration, and invasion. In the mouse skin cancer model, ID4 overexpression reduced tumourigenesis. Future studies should investigate the signalling pathways regulated by ID4 during cutaneous carcinogenesis. Taken together, the results suggest that ID4 might be a target for the treatment or prevention of CSCC. Intra-tumoural delivery of functional ID4 proteins has been shown to inhibit prostate cancer growth [50], supporting the present study. Target strategy that could normalize the epigenetics could also be used for the treatment of CSCC [51,52].

Nevertheless, this study has shortcomings. First, only the methylation and mRNA expression levels of the ID4 gene in skin tissue specimens were detected, and the protein expression was not detected. Examining the protein expression would strengthen the present study, but most of the skin specimens used in this study were from the head and face. Such tumours are usually found early, and they are normally small. Moreover, unlike other types of tumours, the surgical treatment of facial tumours must consider both the efficacy and cosmetic requirements, so the removed surrounding non-cancerous/normal tissue will be as small as possible. Therefore, the amount of available material left after the histopathological diagnosis is often small, including cancerous and non-cancerous/normal tissues, limiting the analyses that are possible. In the present study, most specimens that could be used for the study were spent after extracting the DNA and RNA, and no tissue was left for protein analysis. Second, we used DNMT1-siRNA instead of lentiviral technology to silence the expression of DNMT1 in CSCC cells. We previously silenced DNMT1 alone in CSCC cells using lentivirus expressing DNMT1-shRNA, but the silencing of DNMT1 is lethal to the cells, and subsequent experiments could not be carried out [43,53]. Indeed, DNMT1 is critical to the function of epidermal stem cells. The enrichment of the DNMT1 protein was found in undifferentiated epidermal stem cells, maintaining proliferation tolerance and inhibiting differentiation. In the epidermal model tissues of DNMT1i-HaCaT immunodeficient mice, the absence of DNMT1 resulted in the early differentiation of basal cells and severe proliferation defects until the tissues disappeared [43]. In the human colorectal cancer cell line HCT116, the complete inactivation of DNMT1 leads to severe mitotic defects and cell death during mitosis or after cessation of tetraploid G1 state [53]. Third, only CSCC cells were used, and no analysis was performed on the regulation of ID4 gene by changing the expression of DNMT1, TET1, and TET2 in mouse skin. Such experiments will strengthen future studies. At present, 5-aza-dC is commonly used as a DNMT inhibitor. There have been no relevant studies using 5-aza-dC in animal model skin, and no compounds that regulate the expression of a specific methylation-related protease. Therefore, it is only possible to construct knockout and transgenic animals by knocking out DNMT1 or overexpressing TET1/2 in mice alone. While the lack of such experiments must be mentioned as a limitation of the study, we believe that the findings of this study are of interest. We are planning of conducting such experiments using knockout and transgenic animals in the future.

In conclusion, our data suggest that ID4 is downregulated by UVB irradiation via DNA methylation. ID4 acts as a tumour suppressor gene in CSCC development. ID4 is a potential novel target for CSCC diagnosis and treatment.

Authors’ contributions

JMJ designed and supervised the study. LLM and LFJ organized and carried out the whole experimental process and drafted the manuscript. XYD, SJF, and LQ carried out gene sequencing and the data analysis and figure/table presentation. YXY, FF, WQ, BWB, W Yan, ZK, and W Yi were responsible for sample recruitment. All authors contributed to the final version of the paper.

Declaration of Competing Interest

The authors declare that they have no conflicts of interest.