Epigenetic regulation of HOXA2 expression affects tumor progression and predicts breast cancer patient survival

Abstract

Accumulating evidence suggests that genetic and epigenetic biomarkers hold potential for enhancing the early detection and monitoring of breast cancer (BC). Epigenetic alterations of the Homeobox A2 (HOXA2) gene have recently garnered significant attention in the clinical management of various malignancies. However, the precise role of HOXA2 in breast tumorigenesis has remained elusive. To address this point, we conducted high-throughput RNA sequencing and DNA methylation array studies on laser-microdissected human BC samples, paired with normal tissue samples. Additionally, we performed comprehensive in silico analyses using large public datasets: TCGA and METABRIC. The diagnostic performance of HOXA2 was calculated by means of receiver operator characteristic curves. Its prognostic significance was assessed through immunohistochemical studies and Kaplan-Meier Plotter database interrogation. Moreover, we explored the function of HOXA2 and its role in breast carcinogenesis through in silico, in vitro, and in vivo investigations. Our work revealed significant hypermethylation and downregulation of HOXA2 in human BC tissues. Low HOXA2 expression correlated with increased BC aggressiveness and unfavorable patient survival outcomes. Suppression of HOXA2 expression significantly heightened cell proliferation, migration, and invasion in BC cells, and promoted tumor growth in mice. Conversely, transgenic HOXA2 overexpression suppressed these cellular processes and promoted apoptosis of cancer cells. Interestingly, a strategy of pharmacological demethylation successfully restored HOXA2 expression in malignant cells, reducing their neoplastic characteristics. Bioinformatics analyses, corroborated by in vitro experimentations, unveiled a novel implication of HOXA2 in the lipid metabolism of BC. Specifically, depletion of HOXA2 leaded to a concomitantly decreased expression of PPARγ and its target CIDEC, a master regulator of lipid droplet (LD) accumulation, thereby resulting in reduced LD abundance in BC cells. In summary, our study identifies HOXA2 as a novel prognosis-relevant tumor suppressor in the mammary gland.

Fig. 1. Transcriptome and DNA methylome profiling of human breast cancer tissues.

Workflow of the experimental design to identify tumor suppressor genes in breast cancer (BC) via RNA-sequencing and DNA methylation assays using non-tumor/normal and malignant breast tissue specimens extracting data from our (discovery) cohort of breast patients (a), or from TCGA dataset (b). a, b were created with BioRender.com. Volcano plots of differentially expressed genes (c, e), and differentially methylated genes (d, f) (reported according to their log₂ fold change and p-value in non-tumor/normal versus tumor breast samples), from our cohort of BC individuals (c, d) and from TCGA dataset (e, f). Of note, p-value ≤ 0.05 and |log₂FC|≥ 1 for significant differentially expressed genes (c, e), and p-value ≤ 0.05 and |log₂FC|≥ 0.5 for significant differentially methylated genes (d, f), were used as cut-off. Significant and insignificant modulations are in blue and grey, respectively. See Supplementary Tables S1–S5 for details. FC fold change; TCGA The Cancer Genome Atlas.

Fig. 2. Reduced expression and promoter hypermethylation of HOXA2 in breast cancer.

Scatter plots displaying significant (p-value ≤ 0.05) log₂ fold change (FC) values of differentially expressed (x-axis, |log₂FC| ≥ 1) against differentially methylated (y-axis, |log₂FC| ≥ 0.5) genes in non-tumor/normal versus tumor BC tissues from our cohort of samples (a) and from TCGA dataset (b). c Venn diagram showing overlapping hypermethylated and down-expressed genes in our cohort and TCGA datasets. The genes shared between the two cohorts of samples are listed. Plots showing the correlation between the methylation and expression of HOXA2 in non-tumor/normal and tumor breast samples in the whole cohort (d, e) or more specifically in patients where matched tumor and adjacent non-tumor specimens were available (f), from our cohort (d) or TCGA dataset (e, f). g Violin plot exhibiting the comparison of the expression (Exp) and methylation (Meth) of HOXA2 using individual breast cancer tissue samples and their adjacent non-tumor counterparts (paired samples) from TCGA dataset. d–g Individual values are reported as transcript per million (TPM) for expression, and as beta values for methylation. Coefficients of correlations (R) and significance are displayed on the graphs. h–n Boxplots illustrating the individual level of expression (h, i, l, m) and methylation (j, k) of HOXA2 gene in non-tumor/normal and malignant breast tissues of our cohort of patients (h, j), and derived from TCGA (i, k) and METABRIC (l, m) datasets. Box and whisker plots illustrate the level of expression of HOXA2 as log2 (FPKM + 10^-4) and of HOXA2 methylation as the average of mean beta values (0 = unmethylated, 1 = methylated) for each breast tissue sample. Corresponding data of the statistical analyses are shown in Supplementary Tables S1 and S4–S5. The differential expression of HOXA2 was measured using two probe sets (Pr1 in (l); Pr2 in (m)) in METABRIC dataset. Bars represent mean values ± standard deviation. Data show median, quartiles, and individual values. n Boxplot showing differential levels of methylation within different significant CpG loci of HOXA2 in non-tumor/normal and breast tumor samples derived from our cohort of patients. Each CpG site (according to the Infinium 450k BeadChip probe classification = CpG island regions, CpG-island shores (N and S), and CpG-island shelves) is color-coded in the right bar above. Individual values are expressed as beta values. See Supplementary Table S11 for complementary data. *p-value ≤ 0.05, **p-value ≤ 0.01, ***p-value ≤ 0.001. CpG Cytosine phosphate Guanine; FPKM fragments per kilobase of exon per million fragments mapped; METABRIC Molecular Taxonomy of Breast Cancer International Consortium; Pr1 probe set 1; Pr2 Probe set 2; TCGA The Cancer Genome Atlas.

Fig. 3. Expression level of HOXA2 is a diagnostic and prognostic indicator in breast cancer.

Diagnostic value of HOXA2 expression level in our cohort of breast samples (a), TCGA dataset (b), and METABRIC dataset (c, d) using two probe sets for HOXA2 transcript (c, probe set 1; d, probe set 2) through ROC curve analysis. e Representative immunohistochemistry images showing the staining of HOXA2 protein in human breast tissues (40x). From left to right panel: (i) high expression (more than 20% of positive cells), (ii) medium expression (from 1% to 20% of positive cells) and (iii) low expression (less than 1% of positive cells) of HOXA2 protein. f Graphical representation of the different expression levels of HOXA2 in 96 BC tissues obtained according to histological (G, grading) and clinical (T, tumor status and N, lymph node involvement) parameters. See Supplementary Table S13 for complementary data. Kaplan-Meier analysis of g OS, h RFS, and i DMFS in BC patients expressing high (red curves) versus low (black curves) levels of expression of HOXA2. Survival curves were generated with Kaplan-Meier Plotter online database. See Supplementary Table S14 for complementary data. Affymetrix HOXA2 ID: 214457_at. AUC area under the curve; BC breast cancer; DMFS distant metastasis-free survival; HR hazard ratio; METABRIC Molecular Taxonomy of Breast Cancer International Consortium; N lymph-node status; OS overall survival; Pr1 probe set 1; Pr2 Probe set 2; RFS relapse free-survival; ROC receiver operating characteristic; T tumor status; TCGA the cancer genome atlas.

Fig. 4. HOXA2 functions as tumor suppressor in brs.

HOXA2 mRNA expression was measured by RT-qPCR in a hTERT-HME1 cells transfected with a pool of three HOXA2-specific siRNA sequences (siPool) or the control siRNA (siUNR), and in b MCF7 cells transfected with HOXA2 plasmid (pHOXA2) or the control vector (pCMV). GAPDH was used as endogenous control; a % of knockdown (92.5%) was calculated as = (1−ΔΔCt) x100. Cell proliferation was evaluated by MTT assay in c hTERT-HME1 and d MCF7 cells following HOXA2 silencing (c) or HOXA2 forced expression (d). Cell proliferation after the silencing (e) or overexpression (f) of HOXA2 in breast cell lines was measured by clonogenic assay. Cell migration was detected by transwell assay in absence (g) or presence (h) of HOXA2 in hTERT-HME1 or MCF7, respectively. Cell invasion after HOXA2 knockdown (i) or HOXA2 overexpression (j) was detected by transwell assay. e-j Left panels, representative images of colonies, and transwell inserts stained with crystal violet. Right panels, corresponding quantification of colony formation, and migration and invasion efficiency. Evaluation of cell cycle perturbation by flow cytometry in k synchronized hTERT-HME1 cells transfected with siUNR or HOXA2-siPool, and in l synchronized MCF7 cells transfected with pCMV or pHOXA2 vectors, stained with Hoechst 33342 (10 μM). m Effect of the overexpression of HOXA2 on cell apoptosis by flow cytometry in MCF7 cells transfected with pCMV or pHOXA2 plasmids in the presence or absence of the pan-caspase inhibitor Z-VAD-fmk subjected to a double staining with DAPI and DiOC₆(3) for the detection of dying (DiOC₆(3)^lowDAPI⁻) and dead (DAPI⁺) cells. n Quantification of caspase 3/7 activity in MCF7 cells transfected with pCMV or pHOXA2 in the presence or absence of the pan-caspase Z-VAD-fmk by luminescence. o Western blot analysis of the apoptosis-related protein PARP in MCF7 cells transfected with pCMV or pHOXA2. Antibodies for Flag or β-actin were employed as control of overexpression of HOXA2 or as loading control, respectively. p Bar graph representing the HOXA2 mRNA expression in MCF7 cells after treatment with 5 μM of azacytidine (AZA) at different time points. mRNA expression was measured by RT-qPCR, normalized using GAPDH as endogenous control and calculated according to the 2^−∆∆Ct method. q Bar graph showing the percentage of distribution of synchronized MCF7 cells in each phase of the cell cycle in control condition and after treatment for 72 h with AZA (5 μM). a–q Data represent means ± SD from one representative experiment. Statistical significance was assessed by Student’s t, or by one/two-way ANOVA. AZA azacytidine; RLU relative light units; RT-qPCR reverse-transcription quantitative real time PCR; siUNR unrelated siRNA.

Fig. 5. Depletion of HOXA2 enhances breast tumor growth.

Cell proliferation (a, c) and migration (b, d) indexes of hTERT-HME1 HOXA2 knockout (HOXA2^KO) or wild type (WT) cells (a, b), and of MCF7 cells transfected with pHOXA2 or pCMV (c, d), measured in real-time by xCELLigence^® system. a–d Graphs show mean ± SEM from one representative experiment. Statistical significance was assessed by Student’s t-test. e–h Workflow of in vivo experimentation showing that hTERT-HME1 WT and HOXA2^KO cells were injected subcutaneously into the 4th abdominal mammary fat pad of NOD/SCID mice, and the development of the tumors was monitored over time (e: created with BioRender.com) (e). Tumor growth (f), size (g), and weight (h) at endpoint (n = 9/10 per group, mean ± SEM). f For comparing tumor growth curves, tumor growth p-values were calculated by means of a linear mixed-effect model modeling on the https://kroemerlab.shinyapps.io/TumGrowth/ platform. g Four representative images per group of dissected tumors are shown. Scale bar 1 cm. HOXA2^KO, HOXA2 knockout cell line; WT wild type.

Fig. 6. HOXA2 abundance modulates the expression of CIDEC and PPARγ in lipid metabolism of BC.

(a) Bubble plots showing the significant positive correlation values (R∈[0;1], p-value ≤ 0.05) between the expression of the three selected genes (HOXA2, CIDEC and PPARγ) in breast non-tumor/normal and cancer tissues from our cohort (left graph) and from TCGA dataset (right graph). Bubble color (from light to dark blue) indicates the strength (R value) of the correlation. Of note, a strong correlation is in dark blue and corresponds to a R = 1, while a weak correlation is in light blue and corresponds to a R = 0. Bubble size indicates the significance (p-value) of the correlation; bigger is the size of the bubble, lower is the p-value. (b, c) The mRNA expression of CIDEC (b) and of PPARγ (c) was measured in hTERT-HME1 cells wild type (WT) or knockout (HOXA^KO) for HOXA2 by RT-qPCR. GAPDH was used as endogenous control. (d, e) Western blot analysis of the expression of CIDEC and PPARγ in WT and HOXA2^KO cells (d) and in WT and HOXA2^KO cells transfected with HOXA2 plasmid (pHOXA2) or the control vector (pCMV) (e). Antibodies for Flag (e) or β-actin (d, e) were employed as control of overexpression of HOXA2 or as loading control, respectively. Experimental workflow employed for lipid droplet tracing (f) (f; created with BioRender.com). Confocal microscopy images (scale bar equals 10 μm) (g) and quantification of the number (h) of intracellular lipid droplets (LDs) normalized to cell area and expressed as fold change (FC) to control in hTERT-HME1 WT and HOXA2^KO cells with (+pHOXA2) or without (−pHOXA2) restoration of HOXA2 expression. LDs were stained with BODIPY 493/503 (green), and nuclei were stained with Hoechst (blue). b–e, h Data are presented as mean of ± SD from one representative experiment and analyzed with Student’s t-test or ANOVA. FC fold change; hTERT_WT hTERT-HME1 wild type cells; hTERT_HOXA2^KO hTERT-HME1 HOXA2 knockout cells; RT-qPCR reverse-transcription quantitative real time PCR.