Multi-omics analysis reveals distinct non-reversion mechanisms of PARPi resistance in BRCA1- versus BRCA2-deficient mammary tumors

Abstract

BRCA1 and BRCA2 both function in DNA double-strand break repair by homologous recombination (HR). Due to their HR defect, BRCA1/2-deficient cancers are sensitive to poly(ADP-ribose) polymerase inhibitors (PARPis), but they eventually acquire resistance. Preclinical studies yielded several PARPi resistance mechanisms that do not involve BRCA1/2 reactivation, but their relevance in the clinic remains elusive. To investigate which BRCA1/2-independent mechanisms drive spontaneous resistance in vivo, we combine molecular profiling with functional analysis of HR of matched PARPi-naive and PARPi-resistant mouse mammary tumors harboring large intragenic deletions that prevent reactivation of BRCA1/2. We observe restoration of HR in 62% of PARPi-resistant BRCA1-deficient tumors but none in the PARPi-resistant BRCA2-deficient tumors. Moreover, we find that 53BP1 loss is the prevalent resistance mechanism in HR-proficient BRCA1-deficient tumors, whereas resistance in BRCA2-deficient tumors is mainly induced by PARG loss. Furthermore, combined multi-omics analysis identifies additional genes and pathways potentially involved in modulating PARPi response.

Keywords: BRCA1; BRCA2; CP: Cancer; PARP inhibitor; breast cancer; homologous recombination; multi-omics; therapy resistance.

Graphical abstract

Figure 1

HR restoration drives PARPi resistance in BRCA1-deficient tumors (A) Outline of the generation of matched PARPi-naive and PARPi-resistant KB1P(M) and KB2P tumors and of the experimental approach. (B) Schematic representation of the RAD51 IRIF formation assay. Cryopreserved PARPi-naive and PARPi-resistant tumors were orthotopically transplanted into syngeneic recipient mice, and upon outgrowth to 500 mm3, DNA damage was inflicted by locally applied ionizing radiation (IR) at a dose of 15 Gy. 2 h post-irradiation, tumors were isolated, and fixed tissues were used for RAD51 immunofluorescence imaging. (C and D) Quantification (C) and representative images (D) of the RAD51 IRIFs for the different matched KB1P(M) and KB2P tumor pairs and control KP tumors; IR, irradiated; NIR, non-irradiated; scale bar, 10 μm; data in (C) represented as percentages of positive cells (≥5 foci/nucleus) per imaged area (single data point, typically 100–200 cells/area). n = 5 per imaged area. Data are represented as mean ± standard deviation (SD); ^∗∗∗∗p < 0.0001, ^∗∗∗p < 0.001, ^∗∗p < 0.01 (two-tailed Mann-Whitney U test, followed by Benjamini-Hochberg [BH] correction). (E) Pie charts showing the outcome of the RAD51 IRIF assay in PARPi-resistant KB1P(M) and KB2P tumor cohorts; percentages and numbers of individual tumors analyzed are indicated; n, total number of individual tumors analyzed from the indicated models. p = 0.0001 (two-tailed Fisher’s exact test).

Figure 2

Alterations in previously reported PARPi resistance factors (A) Heatmap (left) of gene expression changes between matched resistant versus naive KB1P(M) and KB2P tumors for the previously reported PARPi resistance genes. Genes for which loss or gain of function have been reported to drive resistance are indicated in blue or red, respectively. In the heatmap, genomic alterations are marked by different symbols (-: copy-number deletion; +: stop-gained SNV; †: frameshift SNV; #: loss-of-function SV). The resistant tumors with transcriptional alterations (lower or higher than 2-fold compared with matched naive tumors for each loss- or gain-of-function driver) are marked by thicker borderlines. The top panel of the heatmap indicates RAD51 IRIF status, which is a proxy for HR status (positive: blue; negative: red; tumors with both positive and negative areas: blue and red; tumors for which RAD51 IRIF was not determined even though it was expected to be negative: pink). The bottom panel of the heatmap indicates the tumors with alterations in at least one known gene. Resistance mechanisms associated with each gene are categorized by different color bars (light green: HR restoration; purple: restoration of fork stability; light blue: PARP signaling). Frequencies for dysregulation of each gene (by either genomic or transcriptional alterations) are shown in the bar plots next to the heatmap (right). The genes preferentially altered in specific tumor types were assessed by the Fisher’s exact test (^∗p < 0.05). (B–D) Scatterplots comparing the alteration frequency of each PARPi resistance factor in the different resistant tumor types. The size of the circle is proportional to the sum of the alteration frequency of the two resistant tumor types compared, and circles are colored if statistically significant (Fisher’s exact test, p < 0.05). The color of the circle indicates the resistance mechanism associated with each resistance factor as mentioned in (A). (E) DIDS outlier scores computed from gene expression data for known resistance factors. Red (positive score) and blue (negative score) indicate upregulation and downregulation of each factor in a subset of resistant tumors compared with naive tumors. The genes with significant DIDS scores are marked with an asterisk (permutation-based exact test, p < 0.05). (F–H) Dot plots of Mad2l2 gene expression in KB1P(M) tumors (F) and Ezh2 (G) and Rif1 (H) gene expressions in KB2P tumors where significant DIDS scores were detected. (I) RAD51 and 53BP1 IRIF status in KB1P(M) PARPi-resistant tumors measured by in situ IRIF assay. (J) Dot plots of Trp53bp1 gene expression in KB1P(M) tumors with 53BP1 IRIF status.

Figure 3

HR-deficient PARPi-resistant BRCA1-KO tumors show increased immune cell infiltration (A) Radar chart showing pathways enriched by upregulated and downregulated genes in KB1P(M) and KB2P resistant tumors compared with naive tumors based on gene expression data. Gene sets from MSigDB Hallmark and KEGG were used for these enrichment analyses. The scale of the axis is represented by false discovery rate (FDR; Fisher’s exact test followed by BH correction). The dotted line indicates an adjusted p value <0.25. (B) Co-functionality network constructed by STRING for the immune-related genes that are significantly upregulated in RAD51-negative KB1P(M) resistant tumors compared with naive tumors. The genes in the network were annotated by MSigDB Hallmark and KEGG and colored depending on the annotated pathways. (C) Quantification and representative images of IHC analysis of markers for different immune cells including leukocytes (CD45), T cells (CD3), B cells (B220), macrophages (F4/80), PD1, and PD-L1 in KB1P(M) tumors. Fisher’s exact test was performed to compare the number of samples with IHC staining levels below or above 5% and 10% between naive and RAD51-positive or RAD51-negative KB1P(M) resistant tumors. Asterisks denote the statistically significant enrichment based on the sample grouping by IHC staining level 5% or 10%.

Figure 4

Multi-omics analysis identifies potential PARPi resistance factors/pathways (A) Schematic overview of the analysis to identify resistance-specific genomic and transcriptional alterations from WE-seq, LCWG-seq, and RNA-seq data. For each PARPi-naive and -resistant RAD51-positive and RAD51-negative KB1P(M) and KB2P tumor, (1) deleterious SNVs, indels, SVs, and focal amplifications/deletions (resistance-specific genomic alterations) and (2) DEG sets identified by either limma or DIDS analysis (resistance-specific transcriptional alterations) were selected. DEG sets were extended to capture the genes with both homogeneous (by limma) and heterogeneous behavior (by DIDS) between naive and resistant tumors. (B) Venn diagrams showing the overlaps between the genes having resistance-specific genomic alterations (SNVs, indels, SVs, and copy-number focal gains and losses) and transcriptional alterations (DEGs) in each resistant tumor type. (C) Heatmap of pathways significantly enriched by the genes with resistance-specific genomic and transcriptional alterations. Enrichment, -log₁₀(FDR) computed by Fisher’s exact test followed by BH correction, by upregulated and downregulated genes are represented by red and green in the heatmap, respectively. Gene sets from MSigDB Hallmark and KEGG were used for these enrichment analyses. (D) Pathway enrichment map for genes with resistance-specific genomic and transcriptional alterations constructed by EnrichmentMap. Node size represents the size (number of genes) of each gene set, and edges represent mutual overlaps between the gene sets (minimum similarity score 0.3). Node and border colors represent enrichment by the genes with resistance-specific transcriptional and genomic alterations, respectively. (E) Venn diagrams showing the overlaps of the genes having resistance-specific alterations with either genomic or transcriptional alterations across the three resistant tumor types. (F) List of DDR-associated genes with resistance-specific alterations that were identified in multiple resistant tumor types.