Harnessing transcriptional regulation of alternative end-joining to predict cancer treatment

Abstract

Alternative end-joining (alt-EJ) is an error-prone DNA repair pathway that cancer cells deficient in homologous recombination rely on, making them vulnerable to synthetic lethality via inhibition of poly(ADP-ribose) polymerase (PARP). Targeting alt-EJ effector DNA polymerase theta (POLθ), which synergizes with PARP inhibitors and can overcome resistance, is of significant preclinical and clinical interest. However, the transcriptional regulation of alt-EJ and its interactions with processes driving cancer progression remain poorly understood. Here, we show that alt-EJ is suppressed by hypoxia while positively associated with MYC (myelocytomatosis oncogene) transcriptional activity. Hypoxia reduces PARP1 and POLQ expression, decreases MYC binding at their promoters, and lowers PARylation and alt-EJ-mediated DNA repair in cancer cells. Tumors with HIF1A mutations overexpress the alt-EJ gene signature. Inhibition of hypoxia-inducible factor 1α or HIF1A expression depletion, combined with PARP or POLθ inhibition, synergistically reduces the colony-forming capacity of cancer cells. Deep learning reveals the anticorrelation between alt-EJ and hypoxia across regions in tumor images, and the predictions for these and MYC activity achieve area under the curve values between 0.70 and 0.86. These findings further highlight the critical role of hypoxia in modulating DNA repair and present a strategy for predicting and improving outcomes centered on targeting alt-EJ.

Graphical Abstract

Introduction

Homologous recombination (HR) is a high-fidelity molecular mechanism crucial for accurately repairing double-strand breaks (DSBs) during the S and G2 phases of the cell cycle, when a sister chromatid is available as a template [1–3]. In cancer cells with defective HR—such as those with pathogenic variants in breast cancer gene 1 (BRCA1) or 2 (BRCA2) [4–6]—there is increased reliance on another DNA repair mechanism known as alternative end-joining (alt-EJ), also referred to as microhomology-mediated end-joining (MMEJ) [7–13]. This dependence forms the basis of the synthetic lethality exploited by inhibition of poly(ADP-ribose) polymerase (PARP) [7, 8], a strategy that has revolutionized cancer care [14]. Thus, assessing HR status has made PARP inhibition (PARPi) a standard treatment strategy for managing various cancers, including those of the breast, pancreas, prostate, and ovary [15–21]. Targeting alt-EJ is being extended to combination therapies and complementary approaches, including inhibition of DNA polymerase theta (POLθ), a central effector of alt-EJ/MMEJ or theta-mediated end-joining [10, 13, 22–26]. Inhibition of POLθ is synthetically lethal with the loss of BRCA1 or BRCA2, can synergize with PARPi, and can also be applied to target cancer cells that developed resistance to PARPi [10, 23, 24]. However, our understanding of the transcriptional regulation of alt-EJ genes, particularly in the context of tumorigenic processes influencing DNA repair, remains limited [9, 27, 28].

The pleiotropic cytokine transforming growth factor β (TGFβ) can enhance the expression of HR components and pathway functionality [29–34]. Gene expression signatures of TGFβ and alt-EJ competency are frequently found to be anticorrelated in tumors, in a cancer cell-autonomous manner [35]. Inhibition of TGFβ signaling induces alt-EJ even in cancer cells that are HR-proficient [33]. TGFβ inhibition results in HR deficit because TGFβ regulates BRCA1 expression via miR182 and cells shift to alt-EJ [30, 33]. Notably, neither alt-EJ execution nor expression of alt-EJ genes (LIG1, PARP1, or POLQ) is miR-182-dependent, which mechanistically separates the effects of TGFβ on HR from those on alt-EJ [33, 35, 36]. Here, by analyzing functional gene modules in breast cancer, we demonstrate that TGFβ signaling suppresses alt-EJ via hypoxia-inducible factor 1α (HIF1α). Hypoxia downregulates the alt-EJ competency signature, including PARP1 and POLQ genes, by altering MYC (myelocytomatosis oncogene) binding at the corresponding promoters. Hypoxia is associated with reduced alt-EJ activity, while HIF1α inhibition increases PARylation. Furthermore, inhibition of HIF1α or depletion of HIF1A expression, in combination with PARPi or POLθi, acts synergistically to reduce the colony-forming ability of cancer cells, including HR-proficient models. This inverse relationship between alt-EJ and hypoxia is applied in cancer treatment prediction by using deep learning of tumor pathology images. The approach provides a framework for classifying tumors based on alt-EJ, hypoxia, and MYC activity, offering novel insights into alt-EJ regulation and its potential for enhancing precision in cancer treatment and for advancing effective combination therapies.

Materials and methods

Gene expression and molecular datasets

Clinical and gene expression [RNA sequencing FPKM (fragments per kilobase of transcript per million mapped reads) UQ (upper quartile-normalized)] data from The Cancer Genome Atlas (TCGA) studies [37] were obtained from the Genomic Data Commons (GDC) Data Portal (https://portal.gdc.cancer.gov). The cancer types are named using the corresponding TCGA study abbreviations (https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables/tcga-study-abbreviations). The pediatric and hematological cancer normalized gene-expression data (reads per kilobase per million mapped reads (RPKM)) were also downloaded from cBioPortal [38]. The HR-deficiency scores in TCGA breast cancer (BRCA) corresponded to the combination of loss of heterozygosity, telomeric allelic imbalance, and large-scale transitions, computed in the HRDsum score [39, 40]. The TCGA kidney renal clear cell carcinomas (KIRC) with driver mutations of VHL were identified from cBioPortal [38], based on OncoKB annotations [41]. The TCGA tumors with HIF1A mutations, including frameshift, missense, and nonsense mutations were also identified from cBioPortal [38]. The single-cell processed count matrices of cancer cell lines induced to epithelial–mesenchymal transition (EMT) were downloaded from the Gene Expression Omnibus (GEO) reference GSE147405 [42]. The PRO-seq (Precision Run-On with sequencing) data of HIF1A^−/− and wild-type (WT) HCT166 cells in normoxia and hypoxia were downloaded from the GEO reference GSE145567 [43]. The Cancer Cell Line Encyclopedia (CCLE) [44] gene expression data were downloaded from the DepMap [45] portal (https://depmap.org/portal/; version 22Q2). The drug sensitivity data were also obtained from the DepMap portal, including the half-maximal inhibitory concentration (IC₅₀) and area under the dose-response curve (AUCdrug) from the CCLE [46] (version 24Q2) and the IC₅₀ values from the Genomics of Drug Sensitivity in Cancer (GDSC) screens [47–49] (release 8.1). PARP1 protein levels, as measured by mass spectrometry in the CCLE and GDSC cell lines, were obtained from the corresponding publications [50, 51]. The optimal number of cell line clusters was computed using the elbow method [52] and samples were clustered using the kmean function in R’s stats package (version 4.3.1). The MYC protein expression values in cancer cell lines were downloaded from The Cancer Proteome Atlas portal (https://tcpaportal.org/tcpa/download.html; reverse phase protein array data, level 4) [53]. SPATA2 software [54] was used to analyze the spatial transcriptomic data of glioblastoma (GBM) [55] and gene signatures were computed using its addGeneSet function. Tissues with an average gene signature score <20% were excluded from the analysis.

Gene signature analysis

The TGFβ (genes n = 50) and alt-EJ (genes n = 36) signatures have previously been reported and functionally validated against mutational signatures and molecular markers [33, 35, 36]. The hypoxia signature consisted of a 15-gene set [56], which was the best-performing among a group of functionally related signatures [57]. The MYC signature (genes n = 355) was derived from a MYC-centered regulatory network independent of a pluripotency network [58]. The signature scores were computed using the single-sample GSEA algorithm calculated using GSVA software [59] (version 1.43.1).

Non-negative matrix factorization and machine learning

Non-negative matrix factorization (NMF) [60] was used for unsupervised data dimension reduction and pattern identification in the TCGA BRCA dataset [37, 61]. The original gene-by-tumor matrix was decomposed into the production of two low-rank non-negative matrices, a basic matrix W and a coefficient matrix H. A rank of matrix factorization, k, from 2 to 50, was examined to preserve the information of the input dataset. The reconstruction error was computed using the Frobenius norm of the deviation between the original matrix X and its approximation WH, expressed as error = ||X – WH||_F. The latent features (rows of matrix H) were fed into a neural network with a structure comprising an input layer of the NMF ranks, two hidden layers, each composed of 128 neurons using a rectified linear unit activation function [62], a layer with a 0.5 dropout rate to mitigate overfitting [63], and an output layer that predicted the HRDsum score. The Adam optimizer [64] was used for training, with a learning rate of 0.001 and batch size of 32. A validation set consisting of 20% of the data was used to monitor overfitting during training. Early stopping with a patience of 10 epochs was employed to prevent overtraining [65]. Spearman correlations were calculated for the training and test datasets to assess the relationship between the predicted and real HRDsum scores [39, 40]. The SHAP (Shapley additive explanations) values were computed as described [66] to assess the importance of each gene module in the model, identifying the most positively and negatively influential sets. The rank of gene weights in each module was evaluated for associations with curated gene sets using the GSEA algorithm (version 4.2.3) [67] with the classic metric, accounting for the uncertain biological impact of gene weights in NMF modules, including those with noncontributory values (i.e. equal to 0).

Multivariate regression analysis

Nineteen cancer types from TCGA were included in the analysis. The ratio of alt-EJ to TGFβ signature scores was calculated and used as the dependent variable in the regression model. The following independent variables were included: age at diagnosis, cancer type, gender, tumor grade, normalized expression level of HIF1A, mutation status of MYC (genomic amplification: yes/no), HRDsum score [39, 40], and gene expression signature scores for hypoxia [56], MYC [58], and cell cycle phases (G1-S, S, G2, G2-M, and M-G1) [68]. Analogous analyses were conducted using expression levels of the cell proliferation marker MKI67 instead of cell cycle phase scores. The interaction term between the hypoxia and MYC signatures was included. The model was computed using the lm function in R software (version 4.3.1).

Cell culture and molecular analyses

The 22Rv1, DU145, MCF10A, MDA-MB-231, PC3, and T-47D cells were obtained from the American Type Culture Collection. The OVCAR4 and A2780 cells were a gift from Thomas Hamilton (Fox Chase Cancer Center) and Josefa Giménez-Bonafé (University of Barcelona), respectively. The microarray gene expression data for the MCF10A cell line are available under the GEO reference GSE8240 [69]. PrimePCR^TM SYBR® Green Assays (Bio-Rad) were used to assess the expression of the alt-EJ genes PARP1 and POLQ in total RNA of MCF10A cells, and gene control was peptidylprolyl isomerase (PPIA). Supplementary Table S1 lists the primer sequences of all genes examined in this study and the TaqMan assays used for validation. Gene expression values were measured using real-time quantitative polymerase chain reaction (qPCR) assays and differences computed by the ΔΔCt method [70]. The antibodies used for western blot were as follows: anti-actin (clone AC-15; dilution 1:2500; Sigma–Aldrich Merck, catalog A5441); anti-HIF1α, (described elsewhere [71]; dilution 1:500); anti-PAR (clone AM80; dilution 1:1000; Merck Millipore, catalog Ab-1/AM80); anti-tubulin (clone B-5-1-2; dilution 1:2000; Sigma–Aldrich Merck, catalog T6079); and anti-vinculin (dilution 1:2000; Sigma–Aldrich Merck, catalog V9131).

Measurement of alt-EJ activity

U2OS GFP-reporter cells (a gift from Prof. J. Stark) were employed to measure alt-EJ repair of I-SceI-induced DSBs [72]. The cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM) medium supplemented with 10% fetal bovine serum (FBS) and standard antibiotics, in an atmosphere of 95% air and 5% CO₂ at 37°C, or in physiological hypoxia, as detailed below. For the assays, 50 000 cells/well were seeded in 24-well plates and transfected the next day with an I-SceI expression plasmid (pCBASceI; Addgene, catalog 26477) using Lipofectamine 2000 (Thermo Fisher Scientific) according to the manufacturer’s instructions. The transfection medium was replaced after 3 h with fresh DMEM 10% FBS and the cells were left for 24 h in normoxic or hypoxic conditions (21% versus 2% O₂; hypoxia using Whitley H35 Hypoxystation, Don Whitley Scientific). The wells for the hypoxic conditions were treated with medium plus dimethylsulfoxide (DMSO) or PX-478 (5 μM). The cells were trypsinized, washed with phosphate-buffered saline (PBS), fixed with 10% formaldehyde, and analyzed for GFP positivity using a MACSQuant Analyzer 16 Flow Cytometer. Singlet cells were gated by plotting forward scatter height against forward scatter area. A logarithmic plot of the gated cells was obtained using PE-A (yellow fluorescence, X-axis) and FITC-A (green fluorescence, Y-axis). Untransfected cells were used to detect autofluorescence.

Colony-forming cell assays and drug synergism

Colony-forming cell (CFC) assays were performed in a 12-well format using the following cell seeding numbers: 22Rv1, 500 cells; A2780, 800 cells; DU145, 200 cells; MDA-MB-231, 300 cells; OVCAR4, 1000 cells; PC3, 300 cells; and T47D, 1500 cells. For high-density seeding, a 20-fold increase in cell counts was applied. In the assays, we evaluated the inhibitory effects of monotherapy and specified drug combinations, including novobiocin (100 μM), olaparib (2 μM), PX-478 (5 μM), rucaparib (2 μM), and talazoparib (1 nM). Cells were treated with the indicated drugs every 24–48 h and cultured for 6–10 days. After treatment, cells were washed with PBS, stained with crystal violet solution, and thoroughly rinsed with tap water. Image acquisition was performed using QuPath software (version 0.5.0). The script for image analysis and Bliss synergy scoring [73] has been deposited on Zenodo and is available at https://doi.org/10.5281/zenodo.14885044. Two-tailed Student’s unpaired t-test was used to assess the significance of the difference between control and treatment groups, with values of P< .05 considered significant. The Bliss independence model [73] assumes that two given drugs produce their effects independently and that significant deviation from their multiplicative effect indicates antagonism or synergism. We implemented the Bliss synergy score from the SynergyFinder [74]. A bootstrapping resampling method was applied to estimate the 95% confidence interval (CI). A score of 10 was interpreted as a 10% response beyond expectation. Dose-response viability matrices were used to evaluate the combination of PX-478 (1.5–40 μM) and rucaparib (2.5–80 μM) over 72 h, with data analyzed using Combenefit software [75].

Lentiviral production and transduction, and gene depletion

HEK-293FT cells were obtained from Invitrogen (Thermo Fisher Scientific) and maintained in DMEM 10% FBS at 37°C, in a 5% CO₂ incubator and split at 70%–80% confluence. Lentiviral particles and transduction were produced following standard protocols as previously described [76], using the packing plasmid psPAX2 (Addgene, catalog 12260) and envelope plasmid pMD2.G (Addgene, catalog 12259). The MISSION^®-validated short-hairpin RNA (shRNA) against human HIF1A expression was reported [77] (TRCN0000003810; Sigma–Aldrich) and the control vector was scrambled-pLKO.1 (Addgene, catalog 136035). The primer sequences for gene expression quantification are depicted in Supplementary Table S1.

ChIP assays

Cross-linking of proteins with DNA, fragmentation, and preparation of soluble chromatin followed by immunoprecipitation with specific antibodies were performed as previously described [78]. Briefly, 5 × 10⁶ cells subjected or not to hypoxia, as indicated, were incubated with 1% (v/v) formaldehyde in pre-warmed PBS for 10 min at 37°C. Cells were then washed in cold PBS, harvested, and lysed to isolate nuclei in hypotonic buffer containing 5 mM PIPES at pH 8.0, 85 mM KCl, and 0.5% NP-40. Nuclei were then resuspended, lysed in a buffer containing 1% sodium dodecyl sulfate, 10 mM ethylenediaminetetraacetic acid (pH 8.0), and 50 mM Tris/HCl (pH 8.1), and sonicated in 15-ml tubes with a Bioruptor sonication device (seven cycles of 30 s ON, 30 s OFF) to yield chromatin sizes of 150–300 bp. Subsequently, 30 μg of DNA/sample was used for immunoprecipitation with 17 μl anti-MYC (N-262; Santa Cruz Biotechnology, catalog SC-764) and with the same amount of anti-rabbit IgG (Cell Signaling Technology). Immunoprecipitated and input DNAs were purified and subjected to qPCR analysis using the primers amplifying the MYC-bound regions found in PARP1 and POLQ genes (Supplementary Table S1 and Text S1). Chromatin immunoprecipitation (ChIP) was quantified by real-time PCR using Roche LightCycler^® 480 RT PCR. The fold enrichment of the target sequence in the immunoprecipitated (IP) relative to the input (Ref) fractions was calculated using the comparative Ct (the number of cycles required to reach a threshold concentration) method, expressed by the term 2^{−(Ct(IP)-Ct(Ref))}.

γH2AX immunofluorescence and quantification

PC3 and T47D cells were cultured on coverslips in six-well plates and treated for 24 h with DMSO (control), rucaparib (2 μM), PX-478 (25 μM), novobiocin (100 μM), or a combination of PX-478 with either rucaparib or novobiocin. Following treatment, cells were fixed in 4% paraformaldehyde for 15 min at room temperature, permeabilized with 0.5% Triton X-100 for 15 min, and blocked in a solution of 5% normal donkey serum, 1% bovine serum albumin, and 0.01% Triton X-100 for 1 h. Cells were incubated overnight at 4°C with an anti-phospho-Ser139 γH2AX antibody (Sigma–Aldrich, catalog 05-636) diluted 1:200 in blocking solution. After primary antibody incubation, cells were washed and stained with a Cy3-conjugated donkey anti-mouse IgG secondary antibody (Jackson ImmunoResearch, catalog 715-165-151) diluted 1:400 in blocking solution containing 1 μg/ml DAPI for nuclear staining. Secondary antibody incubation was performed at room temperature for 1 h. Coverslips were mounted onto microscope slides using Fluoromount-G with DAPI (Invitrogen, catalog 00-4959-52), and images were acquired using a Zeiss LSM 980 confocal microscope. For image analysis, a custom macro was applied to quantify the integrated intensity of phospho-Ser139 γH2AX staining for individual nuclei, as defined by DAPI staining.

Deep learning tumor image modeling

We implemented the weakly supervised method of clustering-constrained-attention multiple-instance learning (CLAM) [79] for whole-slide image (WSI) processing. Feature extraction utilized the universal self-supervised model for pathology (UNI) [80], which was not trained on the histology images of TCGA. The TCGA WSIs were downloaded from the GDC Data Portal (https://portal.gdc.cancer.gov) in scanned virtual slide (SVS) format. Since SVS images contain multiple resolutions, we classified the slides based on the highest magnification level, which, for all selected cancer types, was 20× or 40×. Following the CLAM pipeline, all slides were split into patches. Specifically, the 20× images were divided into patches of 256 × 256 pixels, while the 40× images were split into 512 × 512 pixels to match the 20× resolution when downsampling. The patches were resized to 224 × 224 pixels, which was the input size for UNI, and the features of each patch were extracted into a 1024-dimensional vector. The data were split in an 80–10–10 ratio for training, testing, and validation when each predictor was represented in at least 20% of the primary tumors of a TCGA study. Otherwise, a 60–20–20 ratio was used to ensure enough positive samples for both testing and validation. We applied a dropout rate of 0.5 for training, and learning rate of 3e−5, 10-fold cross-validation, a multibranched CLAM model type, and a limit of 10 epochs, as the models converged relatively quickly. Other parameters were left at their default settings. The same parameters were used to evaluate the algorithm’s performance. All predictions based on signature scores utilized the top quartile as the threshold to distinguish between high and low tumor groups. For HRDsum scores, a cutoff of 42 was applied to both BRCA and serous ovarian (OV) tumors. However, for prostate (PRAD), where very few samples exceeded a score of 42, a cutoff of 20 was used to represent the top 25% of all samples. For generating heatmaps with attention scores, we used the same dropout rate, an overlap of 0.5, and applied a blur effect to produce smooth visual representations. The attention scores, which represent the contribution of each patch within a WSI to predict a class, were extracted for each correctly predicted sample and subsequently Z-score normalized. The top 20% of patches, ranked by the alt-EJ score, were selected for correlation analysis with the hypoxia and MYC scores in tumors.

Results

A cancer cell state defined by alt-EJ–TGFβ anticorrelation

HR-deficient tumors rely on alt-EJ and frequently overexpress components of this pathway [10, 35, 81–84]. To better understand the cancer cell states in which alt-EJ is active, we developed a machine-learning model to predict HR deficiency (HRD). Using the NMF method [60, 85], we identified co-regulated gene modules within the TCGA BRCA dataset [37, 61], resulting in 24 gene modules (Fig. 1A). These modules were used to predict HRD based on the HRDsum score, a genomic metric quantifying HRD in cancer cells by integrating measures of loss of heterozygosity, telomeric allelic imbalance, and chromosomal state transitions [39]. The model demonstrated high predictive accuracy, achieving Spearman correlation coefficients of 0.90 and 0.81 between predicted and observed HRDsum scores in the training and test subsets, respectively (Fig. 1B).

Figure 1.

Cancer gene modules associated with HR status. (A) Graph showing the reconstruction error of the original TCGA BRCA tumor-by-gene matrix using consecutive sets (n) of NMF modules. The set with minimum error is indicated. (B) Scatter plots showing the Spearman correlation coefficient (r) and associated P between the predicted and observed HR-deficiency score (HRDsum) using the machine learning model. Left and right panels show the results of the training and test sets, respectively. (C) Average SHAP values of NMF gene modules, with modules 2 and 11 (highlighted with circles) showing the strongest positive and negative contributions to HRDsum prediction, respectively. (D) GSEA outputs of the alt-EJ (left panel) and TGFβ (right panel) signatures in modules 2 (top panel) and 11 (bottom panel). The GSEA normalized enrichment scores (NESs) and statistical significance (P) are indicated. (E) Histogram depicting the GSEA “Cancer Hallmarks” gene sets positively [false-discovery rate (FDR)-adjusted P < .05; NESs > 1.0; red bars] and negatively (FDR-adjusted P < .05; NESs < 1.0; green bars) associated with module 11. The MYC V1 and V2 gene sets correspond to curated MYC target sets from The Molecular Signatures Database.

To understand the contribution of each gene module to HRD status, we calculated average SHAP values, which quantify the contribution of each module to the prediction. This analysis revealed that module 2 had the strongest positive association with HRDsum, while module 11 exhibited the strongest negative association (Fig. 1C). Consistent with these findings, the gene expression signature associated with alt-EJ activity [33, 35, 36] was positively correlated with module 2 and negatively correlated with module 11 (Fig. 1D, left panels). In contrast, the signature of TGFβ signaling activity [33, 35, 36] was positively correlated with both modules (Fig. 1D, right panels). The positive correlation of TGFβ signaling activity with both modules suggests its pleiotropic function and/or an indirect role in the transcriptional regulation of alt-EJ, independent of small mothers against decapentaplegic (SMAD) factors. Prior studies have shown that BRCA1 loss promotes luminal-to-basal transformation and EMT in breast cancer [86–88], a process mediated by TGFβ receptor 2 [89].

Since module 11 exhibited the strongest negative association with HRDsum, we performed additional analyses to examine its associations with curated gene sets representing various health and disease states [67]. In addition to EMT, this module showed positive associations with stemness-related pathways, including Hedgehog, NOTCH, and WNT/β-catenin signaling, as well as hypoxia signaling (Fig. 1E and Supplementary Table S2). Conversely, it displayed negative associations with gene sets related to E2F and MYC targets (Fig. 1E), suggesting a potential link to cell cycle regulation [90, 91].

Hypoxia signaling is anticorrelated with alt-EJ

TGFβ signaling is known to promote the expression and stabilization of HIF1α (encoded by the HIF1A gene) [92–95]. Consistently, the alt-EJ signature was found to be negatively correlated with both HIF1A expression and a 15-gene hypoxia signature [56] across most cancer types analyzed (n = 25; Fig. 2A). At the single-cell level, transcriptional analysis of four cancer cell lines exposed to TGFβ-1 [42] revealed progressive upregulation of both the TGFβ and hypoxia signatures, along with downregulation of the alt-EJ signature (Fig. 2B), supporting the pan-cancer observation.

Figure 2.

Hypoxia signaling anticorrelates with alt-EJ. (A) Forest plots showing the correlation (PCC and 95% CI) between the alt-EJ signature and HIF1A or the hypoxia 15-gene signature across various cancer types (study acronyms and tumor frequencies are displayed). Nominally significant correlations (P < .5) are indicated in the inset. (B) Plots of the TGFβ (left panel), hypoxia (middle panel), and alt-EJ (right panel) signature scores over time (pseudotime) in single cancer cells (A549, lung; DU145, prostate; MCF7, breast; and OVCA420, ovarian) exposed to TGFβ-1. The slopes from the numerical differentiation of the signatures are shown at the bottom. (C) Scatter plot showing the alt-EJ–TGFβ correlation (PCC) in TCGA KIRC tumors stratified by VHL status, as indicated in the inset. (D) Violin plot showing the overexpression of the alt-EJ signature in HIF1A mutated (MUT) relative to WT TCGA UCEC tumors. The significance (P) of the two-tailed Mann–Whitney test is indicated. (E) Violin plot showing the overexpression of the alt-EJ signature in HIF1A mutated (MUT) relative to WT TCGA pan-cancer (n = 19 cancer types). The significance (P) of the two-tailed Mann–Whitney test is indicated. (F) GSEA outputs of the alt-EJ signature in HCT116 HIF1A WT (left panel) and HIF1A-deteled (right panel) cells exposed to hypoxia versus normoxia. The GSEA–NESs and P-values are indicated. (G) Left panel, normalized activity of alt-EJ as measured by the EJ2-GFP reporter in UOS cells grown in hypoxia relative to hypoxia and exposed to PX-478, as indicated in the inset (n = 4 independent assays). The significance (P) of the one-way ANOVA test with Tukey correction is indicated. There are also indicated the P-values of the two-tailed Student’s paired-samples t-test. Right bottom panel, representative western blot results of PARylation in the above conditions. (H) Top-left panel, representative western blot results of PARylation in the following four UOS cell assays: normoxia; normoxia plus exposure to PX-478; hypoxia; and hypoxia plus exposure to PX-478, as indicated in the inset. Bottom left panel, quantification of PARylation relative to loading control in three independent assays. The significance of the two-tailed Student’s paired-samples t-test relative to normoxia is indicated. Bottom right panel, quantification of alt-EJ as measured by the EJ2-GFP reporter in UOS cells grown in normoxia relative to normoxia and exposed to PX-478.

We then explored the relationship between alt-EJ and TGFβ signatures in TCGA KIRC stratified by the status of the von Hippel–Lindau (VHL) tumor suppressor, which targets HIF1α for proteasomal degradation [96]. While the overall KIRC dataset did not show a significant anticorrelation between alt-EJ and TGFβ [n = 534; Pearson’s correlation coefficient (PCC) = −0.08; P = .078], a significant anticorrelation was observed in VHL-mutant tumors (n = 173; PCC = −0.18; P = .012; Fig. 2C). Next, we compared the alt-EJ signature between tumors with frameshift, missense, or nonsense mutations in HIF1A and those without mutations. Although HIF1A mutations are relatively rare, a substantial number of predicted deleterious variants are identified in TCGA by the OncoKB mutation database [41, 97] in cBioPortal [38] (Supplementary Table S3). The uterine corpus endometrial carcinoma (UCEC) dataset had the largest number of HIF1A mutations (n = 25), with mutated tumors showing significant overexpression of the alt-EJ signature (Fig. 2D). This pattern was also observed across all examined TCGA cancer types, where tumors harboring HIF1A mutations (n = 76) exhibited significantly higher alt-EJ signature expression (Fig. 2E). To further investigate the role of HIF1α in alt-EJ suppression, we analyzed transcriptomic data from HCT116 colorectal cancer cells comparing HIF1A-deleted (HIF1A^−/−) to WT (HIF1A^+/+) cells under normoxic or hypoxic (1% O₂ for 90 min) conditions [43]. WT cells showed suppression of alt-EJ under hypoxia, while this effect was absent in HIF1A-deleted cells (Fig. 2F).

To assess the impact of hypoxia on alt-EJ pathway activity, we used human osteosarcoma U2OS cells with a GFP-based reporter for alt-EJ repair [72] and measured PARylation in cell extracts [98]. PARylation, also known as polyADP-ribosylation—the addition of ADP-ribose molecules to target proteins—is essential for DNA repair and serves as a functional indicator of alt-EJ repair [99–101]. Cells cultured under hypoxia (2% O₂ for 24 h) exhibited a trend toward reduced alt-EJ repair compared with normoxia, with an average suppression of 19% (unpaired two-tailed Student’s t-test, P = .061; n = 4 independent assays; n = 2–3 replicates/assay; Fig. 2G). Treatment with PX-478 (5 μM), an inhibitor of both constitutive and hypoxia-induced HIF1α expression [102–105], restored alt-EJ repair to an average level similar to normoxia, though outliers were noted (Fig. 2G). At the molecular level, hypoxia reduced PARylation, which appeared to be restored by PX-478 treatment (Fig. 2G). To further assess this observation, we quantified PARylation in three additional assays under the same conditions and including cells treated with PX-478 (5 μM) in normoxia. A significant reduction in PARylation was only observed in cells grown under hypoxia compared with those in normoxia (Fig. 2H).

Interplay between hypoxia and MYC influences alt-EJ

Machine-learning modeling predicted that MYC promotes alt-EJ gene expression in TCGA BRCA, as shown by negative GSEA–NES associations (green bars in Fig. 1E). A positive correlation between MYC and alt-EJ expression was observed in 11 of 25 (44%) cancer types analyzed (Fig. 3A, left panel). When a MYC-driven gene signature [58] was analyzed, this positive association was confirmed across all cancer types (Fig. 3A, right panel). Supporting this, 83% (30/36) of the genes in the alt-EJ signature were predicted to be MYC targets based on ChIP data from the ChIP-Atlas [106], and promoter sequences of alt-EJ genes were significantly enriched for MYC binding across various cell types (Fig. 3B).

Figure 3.

Interplay between hypoxia and MYC regulates alt-EJ. (A) Forest plots showing the correlation (PCC and 95% CI) between the alt-EJ signature and MYC expression (left panel) or the MYC-driven signature (right panel) across various cancer types (study acronyms and tumor frequencies are indicated). Nominally significant correlations (P < .05) are indicated in the inset. (B) Violin plot showing MYC-binding (fold change) at the promoter regions (±5 kb from transcription start site) of alt-EJ genes, across tissue (inset). The results correspond to the tissue with ≥1 cell line (dots) with a significant MYC-binding at the alt-EJ promoter set (FDR-adjusted Fisher’s exact test P < .05). The cell lines with the highest MYC-binding in each tissue type are indicated. (C) Graph showing the time-course trends of the expression of the alt-EJ, hypoxia, MYC, and TGFβ signatures in MCF10A cells exposed to TGFβ-1 (500 pg/ml) for 0 to 6 h. The significance (P) of the difference in the hypoxia and TGFβ signature slopes relative to alt-EJ (top inset) and MYC (bottom inset) signature slopes are indicated. (D) Left panel, downregulation of PARP1 and POLQ1 (quantitative reverse transcription polymerase chain reaction (qRT-PCR) assays with specific primers and SYBR) in MCF10A cells cultured in hypoxia, as indicated in the inset. The bars show the mean ± standard deviation (s.d.), and the two-tailed Student’s unpaired-samples t-test P-values are indicated (n = 3 independent assays; 3 replicates/assay). Middle panel, overexpression of GLUT1, OCT4, TGFβ1, and VIM in MCF10A cells cultured in hypoxia, as indicated in the inset. Right panel, western blot showing overexpression of HIF1α in MCF10A cells under hypoxia. The loading control, α-tubulin (TUBA), is also shown. (E) Left panel, downregulation of PARP1 and POLQ1 (TaqMan assays) in MCF10A cells cultured under hypoxia, as indicated in the inset. The bars show the mean ±  s.d., and the two-tailed Student’s unpaired-samples t-test P-values are indicated (n = 2 independent assays; 4 replicates/assay). Middle panel, results for genes involved in HR or NHEJ pathways. Right panel, confirmation of overexpression of GLUT1 and OCT4 in the cells under hypoxia. (F) Left panel, downregulation of PARP1 and POLQ1 (qRT-PCR assays with specific primers and SYBR) in MCF10A cells cultured in normoxia and exposed to CoCl₂, as indicated in the inset. The bars show the mean ±  s.d., and the two-tailed Student’s unpaired-samples t-test P-values are indicated. Right panel, confirmation of GLUT1 overexpression in the cells exposed to CoCl₂. (G) Anticorrelation between the hypoxia signature and PARP1 protein expression in the GDSC and CCLE datasets. The correlation (PCC) and corresponding significance (P) are indicated. (H) Reduced MYC binding at the promoters of the PARP1 and POLQ genes in MCF10A cells grown under hypoxia, as indicated in the inset. The results are shown as the MYC-binding fold change relative to normoxia, including control isotype immunoglobulin (IgG). The unpaired two-sided Student’s t-test P-values are indicated (n = 2 independent assays; 3 replicates/assay). (I) Violin plot showing the distributions of the alt-EJ and hypoxia signature scores (inset) in the tertiles of the ratio of the alt-EJ/TGFβ signature scores in TCGA cancer types (n = 19). Trends in the alt-EJ and hypoxia signature are denoted by dashed lines. The paired two-sided Student’s t-test P-values are indicated for the signature comparisons in each tertile of the alt-EJ/TGFβ ratio. Significance (P) of the ANOVA of the alt-EJ/TGFβ ratio and signature score terms and their interaction are also indicated (bottom inset).

To experimentally assess the influence of TGFβ/hypoxia on alt-EJ, we used MCF10A epithelial cells, which can undergo EMT and acquire stemness upon exposure to TGFβ-1 or hypoxia [107–109]. Transcriptomic analysis of MCF10A cells exposed to TGFβ-1 (500 pg/ml for 6 h) [69] revealed upregulation of the TGFβ and hypoxia signatures, along with downregulation of the alt-EJ and MYC signatures (Fig. 3C). Semi-quantitative gene expression analysis in MCF10A cells under hypoxia (2% O₂ for 8 h) showed significant downregulation of PARP1 and POLQ by approximately 30%, while hypoxia-induced genes exhibited the expected opposite trend (Fig. 3D). Analogous TaqMan-based assays showed similar downregulation of PARP1 and POLQ in hypoxia, though statistical significance was limited to POLQ suppression (Fig. 3E). In agreement with prior studies demonstrating hypoxia-induced suppression of HR and NHEJ components [110, 111], BRCA1, PRKDC, and RAD51 were also downregulated (Fig. 3E, middle panel). Further analysis of MCF10A cells treated with the hypoxia mimetic CoCl₂ (200 μM for 24 h) showed significant reductions in PARP1 (25%) and POLQ (70%) expression (Fig. 3F). In addition, assessing the association between the hypoxia signature and PARP1 protein expression measured by mass spectrometry in hundreds of cancer cell lines (CCLE and GDSC datasets [50, 51]; POLθ was not reported) showed significant anticorrelation (Fig. 3G). Akin to these suppressions, ChIP assays demonstrated hypoxia-induced reductions in MYC binding at the PARP1 and POLQ promoters in MCF10A cells grown under hypoxia (2% O₂ for 8 h) compared with normoxia (Fig. 3H).

The opposing roles of hypoxia and MYC in regulating alt-EJ were further assessed using multivariate regression analyses of the alt-EJ/TGFβ ratio across TCGA studies, accounting for cell cycle signatures [68] or MKI67 expression. Hypoxia and MYC signatures were negatively and positively associated, respectively, with the alt-EJ/TGFβ ratio (β ≤ −0.13; P = 4 × 10⁻¹³ and β ≥ 0.38; P < 2 × 10⁻¹⁶; Supplementary Tables S4 and S5). A significant interaction between hypoxia and MYC signatures (P_interaction = .02) further supported their opposing roles, as confirmed by analysis across tertiles of the alt-EJ/TGFβ ratio (Fig. 3I). Incorporating the HRDsum score into the regression analysis showed that the effects of hypoxia (negative) and MYC (positive) on alt-EJ were independent of HR status (Supplementary Tables S6 and S7). Similarly, correlations between alt-EJ and hypoxia/TGFβ (negative) or MYC (positive) signatures were consistent in TCGA BRCA and TCGA OV tumors, regardless of HRD or HR proficiency as defined by HRDsum (Supplementary Fig. S1). These findings highlight that hypoxia/TGFβ and MYC exert opposing effects on alt-EJ gene expression, independent of HR status.

Combined inhibition of HIF1α and PARP or polθ synergistically reduces CFC

Given that hypoxia signaling is often abnormally activated in cancer [112], we hypothesized that targeting this pathway could relieve suppression of alt-EJ and increase sensitivity to PARP or POLθ inhibitors. To test this, we analyzed the basal signature scores from hundreds of cancer cell lines and their response to PARPi (CCLE dataset [44]). Unsupervised clustering of the alt-EJ, hypoxia, MYC, and TGFβ signatures revealed two major clusters (A and B; Supplementary Fig. S2A). Cluster B, characterized by high alt-EJ and MYC signature expression, and low hypoxia and TGFβ signature expression, showed greater sensitivity (i.e. lower IC₅₀) to PARPi (Supplementary Fig. S2B, C). Similar patterns were observed with three-cell clusters (Supplementary Fig. S2D).

We then performed CFC assays on HR-proficient human breast (MDA-MB-231, T47D) and ovarian (A2780, OVCAR4) cancer cell lines, as well as HR-proficient (DU145, PC3) and HR-deficient (22Rv1, BRCA2-mutant) prostate cancer cell lines, exposed to compound or drug inhibition of HIF1α and/or PARP or POLθ. For these assays, we used PX-478, PARP inhibitors (olaparib, rucaparib, and talazoparib), or the POLθ inhibitor novobiocin [24]. In two-dimensional viability assays, the IC₅₀ values of these drugs were frequently ≥10 μM, consistent with large screening datasets [45, 113] (Supplementary Tables S8 and S9). Cancer cells frequently exhibit some degree of HIF1α activation, even in normoxic conditions [114]. Consistently, the hypoxia signature scores of the selected cell lines (except of DU145, which lacked CCLE gene expression data) grown under normoxia were, on average, intermediate between the scores of cancer cell lines with deleterious HIF1A and VHL variants (Supplementary Fig. S3).

The combination of PX-478 (5 μM; IC₅₀ values between 10 and 42 μM) with rucaparib (PARPi, 2 μM; IC₅₀ values between 20 and 302 μM) showed evidence of synergy (Bliss [73] scores ≥ 16) in five cell lines (excluding A2780 and 22Rv1; n = 2 independent assays; n = 3 replicates/assay; Fig. 4A). In parallel, dose-response matrices over 72 h showed synergistic effects in six cell lines (excluding OVCAR4; Fig. 4B). The synergy was also assessed by combining PX-478 (5 μM) with two other PARP inhibitors (olaparib, 2 μM, and talazoparib, 1 nM) in three cell lines (OVCAR4, PC3, and T47D). These CFC assays yielded Bliss synergy scores ≥ 14 in 6/7 of the settings (n = 2 independent assays; n = 3 replicates/assay; Supplementary Fig. S4). Synergistic effects were also observed in high-confluence settings across all cell lines tested (Supplementary Fig. S5).

Figure 4.

Synergy between HIF1α and PARPi in reducing cancer cell CFC. (A) Results of CFC assays in breast, ovarian, and prostate cancer cell lines exposed to vehicle (DMSO), rucaparib (2 μM), PX-478 (5 μM), and the combination of rucaparib and PX-478. The graphs show the mean ± s.d. The unpaired two-sided Student’s t-test P-values for the comparisons with the drug combination group are indicated. The Bliss synergy score (95% CI) is indicated at the bottom of each panel (Bliss score > 10 indicates synergism). (B) Drug dose matrix of PX-478 and rucaparib in cell viability assays. The color scale indicates synergy or antagonism, with the numbers in the cells denoting the percent difference in loss of viability compared with expected values assuming no synergy. (C) Left panel, western blot results of PARylation and HIF1α levels in DU145 cells exposed for 24 h to vehicle (DMSO), rucaparib (2 μM), PX-478 (25 μM), or the combination of the two drugs. The molecular weights, in kDa, are indicated. The results of the loading control (TUBA) are also shown. Right panel, quantification (arbitrary units, AU) of PARylation under previous conditions. The significance (P) of the two-tailed Mann–Whitney test is indicated (mean ±  s.d.; n = 4 independent assays). (D) Western blot analysis of PARylation and HIF1α levels in DU145 cells exposed for 24 h to vehicle (DMSO), olaparib (2 μM), and/or PX-478 (25 μM). PARylation and HIF1α levels are denoted in each condition as the ratio of the corresponding signal to the loading control (β-actin, ACTB), normalized to the basal condition without drug treatment (set as 1). (E) Left panel, western blot results of total PARylation in DU145 cells grown under normoxia or hypoxia. The results of the loading control (ACTB) are also shown. Right panel, quantification (AU) of PARylation under previous conditions. The significance (P) of the two-tailed Mann–Whitney test is indicated (mean ±  s.d.; n = 5 independent assays). (F) Graph showing the ChIP assays of MYC binding at the promoters of the PARP1 and POLQ genes in DU145 cells grown under normoxic or hypoxic conditions, as indicated in the inset. The results are shown as MYC-binding fold change relative to normoxia, including the control IgG. The unpaired two-sided Student’s t-test P-values are indicated (n = 2 independent assays; 3 replicates/assay).

The molecular alterations underlying the synergy between HIF1α and PARPi were assessed in DU145 cells. PX-478 25 μM (DU145 IC₅₀ = 27.5 μM) for 24 h reduced HIF1α expression and increased PARylation, while the combination with rucaparib strongly suppressed both HIF1α and PARylation (Fig. 4C). An analogous assay substituting rucaparib with olaparib showed similar changes in PARylation and HIF1α levels (Fig. 4D). In addition, DU145 cells exposed to hypoxia (2% O₂ for 8 h) showed a significant decrease in PARylation (Fig. 4E). Furthermore, similar to MCF10A cells, MYC binding at the promoters of PARP1 and POLQ was significantly reduced in DU145 cells grown under hypoxia (2% O₂ for 8 h) compared with those grown under normoxia (Fig. 4F). To determine whether the observed alteration in MYC binding was influenced by changes in MYC expression, we analyzed MYC expression in the seven cell lines under normoxic and hypoxic conditions. This analysis did not reveal a consistent alteration in MYC expression levels (Supplementary Fig. S6A). Furthermore, the expression of MYC and its gene product showed no clear correlation with the hypoxia signature score in the CCLE dataset (Supplementary Fig. S6B).

Next, similar CFC assays were performed using the POLθ inhibitor novobiocin [24]. This drug was applied at a concentration of 100 μM, selected based on its reported high IC₅₀ values and minimal effect on HR-proficient cells [24, 115–119] (IC₅₀ values > 269 μM in the seven cell lines analyzed; Supplementary Table S9). The combination of PX-478 (5 μM) with novobiocin showed synergistic effects (Bliss synergy scores ≥ 12) in 4/7 of the cell lines (n = 2 independent assays; n = 3 replicates/assay Fig. 5A). The combination of novobiocin or rucaparib with PX-478 was also assessed using acriflavine, which reduces HIF1α signaling [120, 121]. Bliss synergy scores ≥ 14 indicated synergism in 67% (4/6) of the tests conducted in three cell lines at relative high confluence (OVCAR4, PC3, and T47D; n = 3 independent assays; n = 3 replicates/assay; Supplementary Fig. S7A). In these assays, exposure to acriflavine increased PARylation while rucaparib reduced this activity (Supplementary Fig. S7B).

Figure 5.

Synergy between HIF1α inhibition or HIF1A expression depletion and POLθ or PARPi in reducing cancer cell CFC. (A) CFC assays for vehicle (DMSO), novobiocin (100 μM), PX-478 (5 μM), and the combination of novobiocin and PX-478. The graphs indicate mean ± s.d. The unpaired two-sided Student’s t-test P-values for the comparisons with the drug combination group are indicated. The Bliss synergy score (95% CI) is indicated at the bottom of each panel (Bliss score > 10 indicates synergism). (B) Left panels, graphs showing the integrated density of phospho-Ser139 γH2AX signal per nucleus: top, PC3 cells; bottom, T47D cells. Data are presented as ≥ 100 nuclei per condition, with the mean value indicated by a horizontal line. The significance of the one-way ANOVA test is shown for all conditions relative to DMSO and for drug combinations relative to PX-478 alone. Right panels, representative images of phospho-Ser139 γH2AX immunodetection under the same conditions. The smaller images show the 4′,6-diamidino-2-phenylindole (DAPI) merge for nuclear staining. (C) Left panels, CFC assays for vehicle (DMSO), rucaparib (2 μM), and novobiocin (100 μM) in T47D cells transduced with pLKO or shHIF1A. Middle panel, semi-quantitative gene expression analysis of HIF1A in the assayed T47D cells; the graphs show the mean ± s.d. (n = 2 independent assays; 3 replicates/assay) relative to pLKO, and the PPIA gene was used as the control. The unpaired two-sided Student’s t-test P-value is indicated. Right panel, the graph shows the mean ± s.d. The unpaired two-sided Student’s t-test P-values for the comparisons with the vehicle are indicated. The Bliss synergy score (95% CI) is indicated at the bottom of each panel (Bliss score > 10 indicates synergism).

To test the prediction that combined inhibition of HIF1α and PARP or POLθ leads to increased DNA damage, we measured phospho-Ser139 γH2AX levels in PC3 and T47D cells. Cells were treated with DMSO, rucaparib (2 μM), PX-478 (25 μM), novobiocin (100 μM), or combinations of rucaparib and PX-478, as well as novobiocin and PX-478, for 24 h. PX-478 treatment alone caused a substantial increase in phospho-Ser139 γH2AX (Fig. 5B), consistent with prior findings [122]. Compared with PX-478 alone, the combination with rucaparib led to a slight, nonsignificant increase in γH2AX signal; however, the combination with novobiocin resulted in a significantly greater increase in phospho-Ser139 γH2AX in both cell lines (Fig. 5B). The underlying basis of these differential effects is unclear but may reflect the distinct roles of PARP1 and POLθ in DSB repair and alt-EJ/MMEJ, and the evaluation of HR-proficient backgrounds.

To assess the effect of targeting HIF1A, we performed CFC assays using a shRNA directed against its expression (shHIF1A) [77] or a scrambled shRNA as a negative control. The combination of shHIF1A with rucaparib (2 μM) or novobiocin (100 μM) showed a stronger inhibitory effect (unpaired two-tailed t-test P < .05) than the single-target assays in the three cell lines analyzed, with evidence of synergy in 50% (3/6) of the assays (Bliss synergy scores ≥ 13; Fig. 5C and Supplementary Fig. S8). However, transduction with shHIF1A alone caused strong inhibition across the three cell lines and, unexpectedly, frequently reduced PARylation and PARP1 expression relative to pLKO (Supplementary Fig. S9). These findings indicate that targeting HIF1A expression may induce additional molecular alterations that influence cell viability and drug combination effects.

Deep-learning prediction of alt-EJ and hypoxia in tumor pathology images for precision cancer treatment

The observed influence of HIF1α on alt-EJ offers potential for improving predictions of targeted therapy outcomes. To explore this, we analyzed spatial transcriptomics data from a GBM dataset [55]. Among the tumors in this study with informative expression of both alt-EJ and hypoxia signatures, 63% (5/8) showed a significant negative correlation between these signatures (PCC ≤ −0.10; P < 1 × 10⁻⁸), while only one tumor exhibited a positive alt-EJ–hypoxia correlation (PCC = 0.12; Fig. 6A).

Figure 6.

Deep learning classification of alt-EJ and hypoxia status. (A) Left panels, representative images of the inferred alt-EJ and hypoxia signature scores in a GBM tissue. Right panel, scatter plot of the alt-EJ–hypoxia signature (standardized score) correlation in GBM tumors (n = 8). The corresponding PCCs are shown: five are significantly (P < .05) negative, two are not significant, and one is significantly positive. (B) Pathology image, and alt-EJ and hypoxia attention maps of a TCGA basal-like (PAM50-classified) breast tumor developed in a women carrier of a BRCA1 pathological variant (tumor ID: TCGA-AN-A0XU). (C) Volcano plot for PCCs between the alt-EJ and hypoxia or MYC signature attention scores in breast tumors correctly classified as “high-alt-EJ and high-hypoxia”, “high-alt-EJ and low-hypoxia”, or “high-alt-EJ and high-MYC”, as indicated in the inset. The tumors with “high-alt-EJ and high-hypoxia” tended to show alt-EJ-hypoxia anticorrelation.

Next, we applied a deep-learning method [79, 80] to predict alt-EJ and hypoxia status, as well as common gene drivers, using tumor pathology images from TCGA BRCA [61], OV [123], and PRAD [124] cancers. The method demonstrated strong performance for both signatures, with area under the receiver operating characteristic curve (AUCROC) values of ≥0.70 in all settings. The AUCROC for alt-EJ was 0.768 ± 0.059, 0.702 ± 0.086, and 0.830 ± 0.0638, and for hypoxia, the AUCROC was 0.732 ± 0.051, 0.719 ± 0.044, and 0.621 ± 0.070 for BRCA, OV, and PRAD, respectively (Supplementary Table S10 and Fig. S10A). The method also performed relatively well in predicting MYC signature activity, with an AUCROC of 0.859 ± 0.058, 0.698 ± 0.064, and 0.821 ± 0.069 for BRCA, OV, and PRAD, respectively (Supplementary Fig. S10B). Among other features analyzed, the best performance were achieved in predicting TP53 mutation and tumor subtype (PAM50 [61]) in BRCA: AUCROC of 0.828 ± 0.068 and 0.892 ± 0.022, respectively (Supplementary Fig. S10C and D).

The analysis of the spatial distribution of inferred alt-EJ activity and hypoxia signaling in tumor images further emphasized their anticorrelation. This relationship was particularly evident in a basal-like breast tumor from a patient with a germline pathogenic BRCA1 variant. While the tumor was predicted to have high levels of both alt-EJ and hypoxia, the attention scores—indicating the model’s focus on specific tumor regions during signature prediction—showed an inverse pattern across tumor patches (Fig. 6B). To further investigate this observation, we examined the correlation between the signature attention scores in breast tumors correctly categorized as either “high-alt-EJ and high-hypoxia” (high–high; n = 37) or “high-alt-EJ and low-hypoxia” (high–low; n = 44). Consistent with the HIF1α-mediated suppression of alt-EJ, tumors in the high–high category frequently showed anticorrelations between alt-EJ and hypoxia scores, whereas positive correlations were common in the high–low category (Fig. 6C). Additionally, we observed that alt-EJ and MYC scores were positively correlated in most tumors (Fig. 6C). The molecular insights gained from this study, combined with the deep-learning approach, have the potential to improve patient stratification for precision treatment and guide future research on synergistic drug combinations.

Discussion

Based on gene expression, molecular, and cellular assays conducted across various cancer models and conditions, our study demonstrates that HIF1α suppresses alt-EJ gene expression and pathway function. Consequently, depletion or inhibition of HIF1α enhances alt-EJ activity, revealing a mechanism by which sensitivity to PARP and POLθ inhibitors can be increased in cancer cells, regardless of HR status. Incorporating these findings with the ability to predict alt-EJ and HIF1α competency using tumor pathology images could improve patient stratification and enhance therapeutic outcomes with PARPi and POLθi. Our pan-cancer observations and mechanistic insights build on the association between resistance to neoadjuvant talazoparib in BRCA1-mutant breast cancer and the basal expression of hypoxia and stem cell-related signatures [125]. Furthermore, our evidence aligns with recent findings that EMT, stemness, and hypoxia-related transcriptional programs contribute to cancer cell adaptation to PARPi [126]. Collectively, our data provide a rationale for future clinical trials combining PARPi or POLθi with HIF1α-targeting compounds [127, 128], which are predicted to synergistically enhance efficacy even in cancers with functional HR-mediated DNA repair.

The opposing effects of HIF1α and MYC on the transcriptional regulation of alt-EJ reveal a functional interaction between these master regulators that may influence cancer cell states and vulnerabilities. HIF1α and MYC generally have contrasting roles in regulating cell cycle progression and proliferation in cancer cells [129–131]. HIF1α promotes cell cycle arrest by partially displacing MYC from the promoters of cell cycle inhibitor genes [132]. By restricting proliferation, HIF1α also competes with MYC for the promoters of genes involved in HR- and mismatch-mediated DNA damage repair [110, 133]. HIF1α further represses HR-associated genes through additional mechanisms [111]. However, complete cell cycle arrest is expected only under severe hypoxic conditions (approximately <0.1% O₂). Most cancer cells within a tumor reside in microenvironments of moderate or mild hypoxia, typically between 0.5% and 5% O₂, with average measurements in human tumors often ranging between 1% and 2% O₂ [134]_. This level of oxygenation can still support HIF1α-mediated survival signals while allowing cancer cell proliferation [135, 136]. Our assays, conducted under conditions of mild or physiological hypoxia (2% O₂), demonstrated reductions in PARP1 and POLQ expression, PARylation, and alt-EJ activity. Additionally, cancer cells often exhibit some degree of HIF1α function even under normoxia, which may be driven by the rewiring of oncogenic pathways and metabolism, and increased oxidative stress, among other factors [137–139]. Recognizing that HIF1α generally suppresses alt-EJ reveals a mechanism by which mutagenic DNA repair processes may be exploited to treat HIF1α-driven aggressive cancer phenotypes [140]. These phenotypes are often resistant to DNA-damaging therapies, potentially due to reduced cell proliferation and the modulation of DNA repair pathways [141–143].

The evidence aligns with the frequent amplification and overexpression of MYC in HR-deficient breast, ovarian, and prostate tumors [123, 144–146]. Additionally, some PARPi-sensitive but HR-proficient cancers are characterized by MYC overexpression [147–149], and this oncogene has been associated with alt-EJ-mediated mutagenesis in leukemias [150]. Furthermore, synthetic lethality between PARP and MYC inhibition, independent of BRCA1/BRCA2 status, has been identified in triple-negative breast cancer [151]. Increased sensitivity to PARPi in MYC-driven cancers may be due to heightened replication stress, resulting in stalled replication forks and DSBs, impaired DSB repair, and/or increased reliance on alternative DNA repair mechanisms [150, 152–156]. Our study expands on these concepts by illustrating the interplay between HIF1α and MYC in the transcriptional regulation of key alt-EJ components, ultimately influencing alt-EJ pathway activity. These findings may also support the evaluation of PARPi and POLθi in tumors exhibiting “high MYC” transcriptional activity, further facilitated by the implemented deep-learning approach. Tumors, or specific regions within tumors, that are, or become, resistant to PARPi or POLθi may exhibit basal activation or acquire activation of HIF1α and/or suppression of MYC. Anticipating these mechanisms could improve therapeutic outcome predictions and guide future combination therapies aimed at overcoming therapeutic resistance.

While our study provides valuable insights into how TGFβ controls transcriptional regulation of alt-EJ by inducing HIFα and suppressing MYC, it is not without limitations. Although consistent molecular associations were observed in tumor data analysis, the findings are primarily derived from in vitro and cell-based experiments conducted under specific conditions that do not fully capture the heterogeneity and complexity of tumors. Additionally, studies incorporating a diverse range of cancer cell types with various genetic alterations, alongside complementary assays, are necessary to establish the broader applicability of these results. While we identified significant drug synergism across several settings, definitive causality remains to be confirmed through further mechanistic studies and an expanded repertoire of precision-targeting approaches, including those involving HIF1α. It is also conceivable that other oncogenic signaling pathways suppress alt-EJ, such those increasing cellular stress, causing metabolic rewiring, and/or modulating of DNA damage repair choice. In this context, synergies have been reported with the combined inhibition of PARP and PI3K/mTOR or MEK signaling [157–163]. Addressing these limitations in our study will be critical for validating and extending the clinical relevance of the findings.

Acknowledgements

We would like to thank Oriol Casanovas and Josep Maria Piulats for useful discussions. Our results are partly based upon data generated by the TCGA Research Network (https://www.cancer.gov/tcga), and we are grateful to the TCGA consortia and coordinators for providing the data and clinical information used here. We are also grateful to Jeremy Stark for providing the U2OS cell lines, to Isabel Puig (VHIO) for providing the shHIF1A, to Laura Soucek and Sandra Martínez (VHIO) for their advice on MYC analysis, and to the investigators who publicly deposited the data used in this work. The anti-HIF1α antibody was a gift from Edurne Berra (CIC bioGUNE).

Author contributions: M.A.P. conceived the study. M.H.B.-H., F.V., A.A., and M.A.P. supervised the study, acquired funding, and wrote the manuscript. R.E. performed bioinformatic, machine and deep learning analyses, and CFC assays. R.E., M.A.P.-C., and I.G. performed bioinformatic analyses. F.M.-J., J.S.-A., S.F.-M., F.M., A.F., A.S., L.R.-A., E.R.-A., I.G., A.B., L.F.-L., A.M.-T., and R.K. performed cellular and molecular analyses. R.T.S. and G.P.V. performed ChIP assays. M.M.-I., X.C.W., E.C., J.A.M., M.L.-C., P.M., I.R., A.R., I.F., A.V., X.S., J.C., M.A., A.H.-R., V.S., K.M., R.H., and J.B. performed data interpretation and provided resources.

Supplementary data

Supplementary data is available at NAR Cancer online.

Conflict of interest

M.H.B.-H. and M.A.P. have filed a patent application, submitted by the University of California and IDIBELL, covering the use of TGFβ and alt-EJ signatures for predicting cancer response to genotoxic therapies.

Funding

This study was funded by the “GINKGO Apac del Berguedà” patient association; the Instituto de Salud Carlos III [co-funded by the European Regional Development Fund (ERDF), “A way to build Europe”] grants PI21/0136 and PI24/01327 to M.A.P., and PI19/00342 and PI23/00513 to A.A.; the Asociación Española Contra el Cáncer (AECC) grant LABAE21170 to A.A.; the Department of Defense PC210340 PCRP-IDA grant to A.A.; the Ministerio de Ciencia, Innovación y Universidades, Agencia Estatal de Investigación [co-funded by the European Regional Development Fund (ERDF), “A way to build Europe”] grants PID2020-117815RB-I00 and PID2023-150836OB-100 to F.V., and CNS2023-145615 to A.A.; and the Generalitat de Catalunya, Agència de Gestió d'Ajuts Universitaris i de Recerca (AGAUR) grants SGR 2017-449 and 2021-00184 to F.V. and M.A.P., and 2021-00895 to A.A. R.E. and M.A.P.-C. were supported by contracts from the Departament de Salut, Generalitat de Catalunya, PERIS-PFI SLT017-20-000076 and PERIS-Suport SLT017/20/000072, respectively, and A.S. was supported by Generalitat de Catalunya AGAUR fellowship 2022-FI-B01068. We also acknowledge the support of the Centres de Recerca de Catalunya (CERCA) Program to IDIBELL and IDIBGI. The open-access publication of this article was funded by the Instituto de Salud Carlos III grants PI24/01327 to M.A.P., PI23/00513 to A.A., and PID2023-150836OB-100 to F.V.

Data availability

The data are available in the article, figures, and additional files. The code for the machine-learning analysis, Bliss synergy scoring, and deep-learning models have been deposited in https://doi.org/10.5281/zenodo.14885044. The MCF10A microarray data were deposited in GEO reference GSE8240.

Author notes

These authors should be regarded as Joint Second Authors.