Endometrial tumorigenesis involves epigenetic plasticity demarcating non-coding somatic mutations and 3D-genome alterations

Abstract

Background: The incidence and mortality of endometrial cancer (EC) is on the rise. Eighty-five percent of ECs depend on estrogen receptor alpha (ERα) for proliferation, but little is known about its transcriptional regulation in these tumors.

Results: We generate epigenomics, transcriptomics, and Hi-C datastreams in healthy and tumor endometrial tissues, identifying robust ERα reprogramming and profound alterations in 3D genome organization that lead to a gain of tumor-specific enhancer activity during EC development. Integration with endometrial cancer risk single-nucleotide polymorphisms and whole-genome sequencing data from primary tumors and metastatic samples reveals a striking enrichment of risk variants and non-coding somatic mutations at tumor-enriched ERα sites. Through machine learning-based predictions and interaction proteomics analyses, we identify an enhancer mutation which alters 3D genome conformation, impairing recruitment of the transcriptional repressor EHMT2/G9a/KMT1C, thereby alleviating transcriptional repression of ESR1 in EC.

Conclusions: In summary, we identify a complex genomic-epigenomic interplay in EC development and progression, altering 3D genome organization to enhance expression of the critical driver ERα.

Keywords: 3D genome organization; Endometrial cancer; Epigenetic plasticity in tumor development; Estrogen receptor; Gene regulation; Non-coding somatic mutations.

Fig. 1

ERα cistrome changes upon tumorigenesis. a Schematic workflow of the multi-omics approach applied in this work. Boxplot depicting the distribution of peak number detected by ERα (b) or H3K27ac (c) ChIP-seq in healthy (blue) or tumor (orange) endometrial tissues. d Venn diagram of the overlap between consensus ERα binding sites detected in healthy (blue) and tumor (orange) tissues. Consensus peaks of each group have been defined as peaks present in at least 75% of the samples belonging to that specific group. e MA plot of differential binding analyses performed on the consensus ERα ChIP-seq peaks (healthy vs tumor). Differential binding sites (FDR ≤ 0.05 and |log₂(FoldChange)| ≥ 1) are highlighted in pink. Tornado plot of the ERα (f) or H3K27ac (g) ChIP-seq signal at tumor-depleted (blue, upper blocks) or tumor-enriched (orange, lower blocks) ERα consensus peaks for all the 4 healthy (blue, left heatmaps) and 5 tumor (orange, right heatmaps) endometrial tissues. On the top of each heatmap is plotted the average density signal for each region group. h Stacked bar plot depicting the genomic localization frequency of tumor-depleted and tumor-enriched ERα ChIP-seq consensus peaks. i Stacked bar plot displaying the distance to associated gene TSS frequency distribution of tumor-depleted and tumor-enriched ERα ChIP-seq consensus peaks. j Wordcloud of the top-50 enriched motifs at tumor-depleted (blues) and tumor-enriched (oranges) ERα consensus peaks. Size and color intensity are proportional to the −log₁₀(E-value). k The upper part of the heatmap shows the individual tumor-depleted (blues) and tumor-enriched (oranges) ERα consensus peaks colored by signal intensity. On the lower part, each black bar represents an overlapping peak of ChIP-seq data of several targets publicly available for the Ishikawa endometrial cancer cell line. l Heatmap of the percentage of tumor-depleted or tumor-enriched ERα consensus peaks overlapping with each ChIP-seq target in Ishikawa cells from (k). Ranking is performed by descending number of overlaps in tumor-depleted peaks

Fig. 2

3D genome landscape is remodeled during tumorigenesis. a Schematic of the Hi-C library preparation starting from fresh frozen tissues. Ten 30-µm-thick slices of flash-frozen tissue are cross-linked 25 min in 2% formaldehyde. Then, electronically homogenized tissues are filtered using a 75-µm cell strainer and subjected to Hi-C library preparation. Hi-C libraries are then sequenced and resulting reads are analyzed by the snHiC [24] pipeline. b On the left, 40-kb-resolution matrices of the average Hi-C score in tumor (n = 3) and healthy endometrial tissues (n = 3) at chromosome 6. On the right, differences of the scores shown on the left side (tumor − healthy). c Relative Hi-C contact probability (RCP) as function of the distance for each individual sample (40-kb resolution). d Violin plot of the distribution of the short-range (<2 Mb) over long-rage (>2 Mb) Hi-C contacts ratio in each individual sample (40-kb resolution). For each comparison, the Wilcoxon’s test P value is indicated. e Distribution of the Hi-C loop distance (10-kb resolution) for loops overlapping with tumor-depleted, tumor-enriched, or shared ERα consensus peaks in combined healthy (blue) or tumor (orange) endometrial tissues. f PE-SCAn results in healthy and tumor tissues at common, tumor-depleted, and tumor-enriched ERα binding sites. For each panel, top and middle rows depict the 3D and 2D, respectively, representation of the Hi-C contact frequency distribution in healthy (left) and tumor (right) tissues; lower row shows the difference of Hi-C contact frequencies between tumor (orange) and healthy (blue) tissue scores. g Compartment polarization ratio (100-kb resolution), defined as (AA + BB)/(AB + BA), for each individual sample. h Top: saddle plot of A/B compartments interactions as computed in g for each individual sample. Bottom: difference of saddle-score with the reference H.005.A2 (healthy tissue), where orange indicates a higher score in tumor samples while purple indicates a higher score in healthy tissues. i Sankey plot depicting the proportion of compartment state transition (100-kb resolution) of healthy compartments (left) upon tumorigenesis (right). j Upset plot showing the overlaps between different compartment transition states (100-kb resolution) in healthy compared to tumor as in g (A-to-A, B-to-B, A-to-B, B-to-A) and different categories of ERα consensus peaks. k Stacked bar plot of the proportion of compartment transition status as in g for individual genomic bins overlapping with only one of the different categories of ERα consensus peaks

Fig. 3

Tumor-enriched ERα binding sites represent non-coding regions target of somatic mutation in metastatic endometrial cancer tumor. Genomic localization (a), mutational signature (b), and counts (c) distribution of the somatic mutations detected by WGS in the primary endometrial cancer samples form a selection of samples from the TCGA-UCEC cohort (n = 41). d Stacked bar blot depicting the fraction of ERα peaks overlapping with somatic mutations identified in the primary cohort. e Variant set enrichment (VSE) analysis depicting the enrichment of differentially bound ERα sites over endometrial cancer risk loci identified by genome-wide association study (GWAS, P < 1 × 10⁻⁵) [27]. Density plot represents distribution for the Z-score from matched controls defining the null distribution. Blue (tumor-depleted) or orange (tumor-enriched) vertical lines represent observed enrichments, with tumor-enriched enrichment being statistically significant (P_adj = 0.0060). Gray vertical lines define 0, 25, 75, and 100 percentiles of the distribution, while dotted black line indicates the median. f Scheme of the number metastases and metastatic site, for all metastatic samples used for the WGS analyses. Genomic localization (g) and mutational signature (h) distribution of the somatic mutations detected by WGS in the metastatic samples described in f. i In the outer plot, it is shown the number of different types of somatic variants detected in each metastatic sample. In the inner plot, stacked bar plot of the most-frequently protein sequence mutated genes and relative number and type of somatic variants identified in metastatic samples. j Stacked bar blot depicting the fraction of ERα peaks overlapping with somatic mutations identified in the metastatic cohort. k MA plot of differential expression analyses in Ishikawa cell lines upon 8 h of 10 nM β-estradiol (E2) stimulation. Black dots indicate genes which promoter has been linked, by H3K27ac Hi-C analyses, to a tumor-depleted or tumor-enriched ERα-bound regulatory element bearing somatic mutations in metastatic samples (n = 311). Differentially expressed genes (|Fold Change expression| ≥ 1.5 and P_adj < 0.05) upon estradiol stimulation have been labeled and highlighted by a red dot (n = 12). l Ishikawa endometrial cancer cells were treated or not for 6 or 12 h with 10 nM estradiol. Whole-cell extracts were analyzed by immunoblotting with antibodies against ERα and GAPDH (loading control). m Progression-free Kaplan-Meier curve of endometrial cancer patients (TCGA data) divided into two groups using the median of ESR1 expression (FPKM) as cutoff

Fig. 4

Short-range chromatin contacts at the ESR1 locus are stronger in tumors. a ChIP-seq genomic tracks for ERα (upper block) and H3K27ac (lower block) in healthy (blue) and tumor (orange) endometrial primary tissues at the ESR1 locus. ESR1 Enhancer 1 and Enhancer 2 are indicated, as well as the mutations found by WGS analyses in metastatic samples. b Representation of the averaged 40-kb resolution Hi-C matrix at the ESR1 locus for the 3 healthy (top) and 3 tumor (middle) tissues or the score difference (bottom) tumor − healthy, where orange indicates higher scores in the tumors while purple higher scores in healthy. The matrices scores are dived by the sum of the matrix. Black lines indicate the topologically associated domains (TADs) identified. c 4C-seq genomic tracks at the ESR1 locus using the ESR1 TSS as view-point (VP) for 2 healthy (top, blue) and 3 tumor (middle, orange) endometrial tissues, and the average difference of score tumor − healthy (bottom). With green arcs are depicted loops detected by H3K27ac HiChIP in Ishikawa endometrial cancer cells. The ribbon around the 4C-seq signal lines indicates the standard error mean (SEM) among biological replicates. d In the top row, observed/expected matrices of sequence-based machine learning prediction in a ±500-kb window surrounding the ESR1 locus. In order from left to right can be found: wild-type sequence, point mutation in a genomic desert (negative control), deletion of the full Enhancer 1 sequence (positive control), introduction of SNVs found by WGS analyses of metastatic endometrial cancer samples. On the bottom row is showed the difference of observed/expected score over the wild-type sequence, where orange indicates higher scores in the altered sequence and purple a higher score in the wild-type one

Fig. 5

EHMT2/G9a is a negative regulator of ESR1 expression in endometrial cancer. a Schematic workflow used to perform DNA-oligo protein pull-downs. Biotin-conjugated wild-type or mutated DNA-oligos are immobilized on streptavidin magnetic beads and mixed with Ishikawa nuclear lysates. Captured proteins are then dimethyl labeled and analyzed by mass spectrometry. b A DNA oligo with the sequence of the ±25 bp surrounding the ERα binding site in the ESR1 Enhancer 1, in wild-type or chr6:152,002,679–TCT-to-T form, was used to perform DNA-oligo protein pull-downs in Ishikawa cells as described in a. The scatter plot shows the log₂ ratios of all identified and quantified proteins in both experiments plotted against each other. Proteins significantly enriched at the wild-type sequence are highlighted in red, and proteins significantly enriched at the mutant sequence are highlighted in blue. c Gene expression correlation heatmap for all corresponding 24 differentially bound proteins identified in b in 589 endometrial cancer patients (TCGA). Dashed lines indicate the separation between positive and negative correlation scores. Genes are ranked by correlation with ESR1 gene expression. d Progression-free Kaplan-Meier curve of endometrial cancer patients (TCGA data) divided into two groups using the median of EHMT2/G9a expression (FPKM) as cutoff. e EHMT2/G9a and IgG (negative control) ChIP followed by western blot in Ishikawa endometrial cancer cells. Antibodies against EHMT2/G9a and ERα have been used. For EHMT2/G9a, two different exposure images were used for input and IP, as indicated by the vertical dashed bar. f EHMT2/G9a ChIP in Ishikawa cells stimulated for 6 h with 10 nM β-estradiol. Bar plot shows percentage of enrichment over the input (% Input) at the ESR1 Enhancer 1, ESR1 promoter, and CDK12 promoter (negative control) analyzed by quantitative PCR (qPCR). Mean of 2 independent experiments is shown. g Ishikawa endometrial cancer cells were stimulated for 72 h with 10 nM β-estradiol in combination or not with 100 nM ICI-182,780 (fulvestrant, negative control). In all conditions, cells where incubated either with a non-targeting (NT) siRNA or with an siRNA against EHMT2/G9a. Then, whole-cell extracts were analyzed by immunoblotting with antibodies against ERα, EHMT2/G9a, and GAPDH (loading control). h Differential ERα peak centered EHMT2/G9a average ChIP-seq density profiles in Ishikawa cells upon 6 h treatment with 10 nM β-estradiol (purple) or DMSO (blue, control). Ribbon indicates SEM. Paired Wilcoxon test was performed on the average score of the highlighted area for each individual region analyzed; P values are indicated. i Schematic representation of the ESR1 gene regulation working model. Upon tumorigenesis, ERα is re-located to an endometrial cancer specific ESR1 enhancer (Fig. 1f–g). ERα interacts with EHMT2/G9a (e) and binds both enhancer and promoter regions of ESR1, in a hormone-dependent fashion (f). In a subset of metastatic endometrial cancers, a somatic mutation is acquired at the tumor-specific ESR1 enhancer (Fig. 4a). In vitro analyses probing this mutation revealed a loss of EHMT2/G9a DNA binding capacity (b), and EHMT2/G9a knockdown in endometrial cancer cell lines leads to an increase of ESR1 expression levels (g)