Abstract
We examined gene expression, germline variant, and somatic mutation features associated with pathologic response to neoadjuvant durvalumab plus chemotherapy in basal-like triple-negative breast cancer (bTNBC).
Germline and somatic whole-exome DNA and RNA sequencing, programmed death ligand 1 (PD-L1) IHC, and stromal tumor-infiltrating lymphocyte scoring were performed on 57 patients. We validated our results using 162 patients from the GeparNuevo randomized trial.
Gene set enrichment analysis showed that pathways involved in immunity (adaptive, humoral, innate), JAK–STAT signaling, cancer drivers, cell cycle, apoptosis, and DNA repair were enriched in cases with pathologic complete response (pCR), whereas epithelial–mesenchymal transition, extracellular matrix, and TGFβ pathways were enriched in cases with residual disease (RD). Immune-rich bTNBC with RD was enriched in CCL-3, -4, -5, -8, -23, CXCL-1, -3, -6, -10, and IL1, -23, -27, -34, and had higher expression of macrophage markers compared with immune-rich cancers with pCR that were enriched in IFNγ, IL2, -12, -21, chemokines CXCL-9, -13, CXCR5, and activated T- and B-cell markers (GZMB, CD79A). In the validation cohort, an immune-rich five-gene signature showed higher expression in pCR cases in the durvalumab arm (P = 0.040) but not in the placebo arm (P = 0.923) or in immune-poor cancers. Independent of immune markers, tumor mutation burden was higher, and PI3K, DNA damage repair, MAPK, and WNT/β-catenin signaling pathways were enriched in germline and somatic mutations in cases with pCR.
The TGFβ pathway is associated with immune-poor phenotype and RD in bTNBC. Among immune-rich bTNBC RD, macrophage/neutrophil chemoattractants dominate the cytokine milieu, and IFNγ and activated B cells and T cells dominate immune-rich cancers with pCR.
We found that high tumor mutation burden and immune-rich microenvironment are independently associated with pathologic complete response (pCR) to anti–programmed death ligand 1 therapy plus chemotherapy in basal-like triple-negative breast cancer (bTNBC), whereas lack of pCR and immune-poor phenotype are associated with higher expression of TGFβ pathway and epithelial/mesenchymal markers. Immune-rich bTNBC with residual disease are characterized by higher expression of CCL-3, -4, -5, -8, -23, CXCL-1, -3, -6, -10, IL1, -23, -27, -34, and more abundant in macrophage markers. Immune-rich bTNBC with pCR are characterized by activated T- and B-cell markers and expression of IFNγ, IL2, -12, -21, CD79A, and GZMB. No mutation in single genes was associated with response, but cancers with pCR had significantly more germline variants and somatic mutations in the PI3K, DNA damage repair, MAPK, and WNT/β-catenin pathways that could affect chemotherapy sensitivity.
Introduction
Multiple randomized trials demonstrated increased pathological complete response (pCR; ypT0is/ypN0) rates when an anti–programmed cell death protein 1 (PD-1; pembrolizumab) or anti–programmed death ligand 1 (PD-L1; atezolizumab, durvalumab) antibody is included with standard-of-care neoadjuvant chemotherapy in triple-negative breast cancer (TNBC; refs. 1–5). Patients who achieve pCR have excellent long-term survival, and two of these randomized trials also reported statistically significant improvement in recurrence-free survival with immunotherapy (2, 6). The pCR rates after combined neoadjuvant immunotherapy + chemotherapy range between 44% and 65% depending on the type of chemotherapy regimen indicating that many patients continue to have residual disease (RD) after therapy.
Several studies examined molecular predictors of pCR in neoadjuvant immunotherapy + chemotherapy trials and demonstrated that cancers with higher levels of tumor-infiltrating lymphocytes (TIL) and greater expression of PD-L1 protein on immune cells have higher pCR rates compared with cancers with lesser immune infiltration (7). The presence of TILs and PD-L1 expression are associated with the expression of a broad range of immune gene expression signatures that also predict pCR (8, 9). Tumor mutation burden (TMB) recently emerged as an independent predictor for pCR (10). However, all of these markers are predictive of pCR with or without immunotherapy (11–13), no validated molecular markers exist that could identify patients who selectively benefit from inclusion of immunotherapy with their neoadjuvant chemotherapy. We recently reported that MHC class II protein expression on tumor cells may identify cancers that are selectively benefitting from neoadjuvant immunotherapy, but this observation will require independent validation (14).
The goal of this study was to comprehensively characterize molecular features of basal-like TNBC (bTNBC) that achieved pCR after neoadjuvant anti–PD-L1 therapy plus chemotherapy compared with cases with RD in the overall study population and in the immune-rich subset. We focused on bTNBC to minimize molecular heterogeneity in our relatively small sample set and because of the clinical differences between bTNBC and other less frequent TNBC molecular subtypes (15, 16). We performed whole-exome and whole-transcriptome RNA sequencing (RNA-seq) along with histologic assessment of pretreatment needle biopsies collected during a single arm phase II clinical trial to identify candidate markers of response (7). We assessed the association between our candidate response markers on the chemotherapy alone and chemotherapy plus durvalumab arms of the GeparNuevo randomized trial.
Materials and Methods
Patient population and biospecimens
Pretreatment core needle biopsies for research were obtained from patients with stage I to III TNBC who enrolled in a single arm neoadjuvant clinical trial (NCT02489448) and received durvalumab concurrent with weekly nab-paclitaxel x 12 followed by durvalumab plus dose dense doxorubicin/cyclophosphamide x 4 treatments. Primary efficacy results were previously published (7). Sixty female patients were enrolled in the trial, 2 patients were not evaluable for pathologic response, and 1 patient withdrew consent; therefore, the biomarker population includes 57 patients (pCR n = 26, RD n = 31; Supplementary Fig. S1). All patients provided written informed consent for research on their donated tissues, including germline DNA sequencing. Ethical approval was obtained from the Yale Human Investigations Committee (HIC; Yale University, New Haven, CT; HIC no. 1409014537). The validation data included targeted mRNA sequencing results of 2,559 transcripts generated from pretreatment biopsies of 162 patients enrolled in the GeparNuevo trial (NCT02685059) who received durvalumab or placebo every plus nab-paclitaxel x 12 weeks, followed by durvalumab or placebo plus epirubicin/cyclophosphamide x 4 treatments (2). The GeparNuevo protocol was approved by the respective ethics committee, institutional review board, and national competent authority. All studies were conducted in accordance with the Declaration of Helsinki.
Isolation of RNA and DNA
For the Yale cohort, RNA and DNA were extracted from one biopsy collected in RNAlater™ (Qiagen) and stored at -80˚C. After homogenization with the TissueLyser II bead-milling system (Qiagen), DNA was isolated using the AllPrep DNA/RNA/miRNA Universal Kit, and the flow-through RNA was extracted with RNeasy Plus Kit (Qiagen) following the manufacturer's instructions. The quality and concentration of isolated DNA and RNA were tested on the Agilent 2100 Bioanalyzer system. DNA and RNA-seq were performed at the Yale Center for Genome Analysis.
RNA-seq and data processing
Paired-end sequencing of 100-bp fragments of total RNA for a targeted depth of 50 million reads was performed using the Illumina NovaSeq platform. Quality was assessed using FastqQC v 0.11.5 (17), adapter sequences were trimmed with Trimmomatic v0.36 (RRID:SCR_011848; ref. 18). Sequencing reads were aligned to human genome, hg38, with STAR v2.5.3a (RRID:SCR_004463; ref. 19) using two-pass mode and default parameters; alignment quality and strandedness was checked using RSeQC v2.6.4 (RRID:SCR_005275; ref. 20). Gene expression was quantified using RSEM v1.3.0 (RRID:SCR_013027) (21) and ENSEMBL release 91 (RRID:SCR_002344) was used to annotate reads to human genes. One specimen was excluded due to poor RNA quality. Molecular subtyping was performed using the AIMS, SCMOD2, and PAM50 methods. All methods identified the same 6 cases as non–basal-like and outlier analysis using principal component analysis (PCA) and uniform manifold approximation and projection (UMAP) of transcriptomic data confirmed these cases as distinct from the remaining samples. (Supplementary Figs. S1 and S2). All six nonbasal cases had RD. In order to work with a molecularly homogeneous bTNBC set we excluded these 6 cases from further analysis, resulting in 50 bTNBC cases (pCR n = 25, RD n = 25).
Gene expression differences between pCR and RD samples were determined using “DESeq2” R package (RRID:SCR_000154; ref. 22). To adjust for variable tissue composition from case to case, we added a previously published stromal score that quantifies tumor stromal content calculated using the ESTIMATE R package (23). Differentially expressed genes were defined as log-fold change > 1, the Benjamini–Hochberg method was used to adjust for multiple comparisons and adjusted P < 0.05 was considered significant. Gene set enrichment analysis (GSEA) was implemented using the fgsea R package (24). To quantify biological and immune processes, the NanoString Hallmarks of Cancer and Biological Pathways and Processes gene sets and a collection of previously published immune gene signatures (5, 8) were used (Supplementary Tables S1–S2). Immune signature-based classification into immune-high versus immune-low status was performed using the median values of the Tumor Inflammation Signature (TIS), GeparSixto, NHI 5-gene, STAT1 gene signatures, and expression of IFNγ single gene, respectively. Gene signature expressions were compared between pCR and RD groups using the Mann–Whitney test. Multivariate association between gene signatures expressions and pCR were assessed using logistic regression adjusted for age (continuous variable), tumor size (T1 vs. ≥ T2), nodal status (N0 vs. N1–N3), and stromal score (from ESTIMATE algorithm). Benjamini–Hochberg corrected P < 0.05 was considered significant.
For the GeparNuevo validation cohort, formalin-fixed, paraffin-embedded (FFPE) tissues were processed using an HTG EdgeSeq instrument (HTG Molecular Inc.) with the Oncology Biomarker Panel according to the manufacturer's instructions as previously described (2, 10). Fisher exact test and Pearson χ2 were used to evaluate categorical variables [pCR status; stromal TILs (sTIL; high, low); ref. 10]. Univariate and multivariate logistic regression models with adjustments for clinical covariates (as previously described) were used for assessment of predictive value of genes for pCR (10).
Whole-exome sequencing
Genomic DNA (1 μg) from tumor biopsies and matched peripheral blood buffy coats of 57 patients were sheared to a mean fragment length of 140 bp and exomes were captured using the NimbleGen SeqCap EZ v2 kit. The resulting library was sequenced on an Illumina HiSeq 4000 instrument in paired-end 75-cycles mode to achieve an average target sequencing depth of 232x for tumor samples and 207x for matched normal samples. Reads were filtered by Illumina CASAVA 1.8.2 software, trimmed at the 30 end using FASTX v0.0.13 (RRID:SCR_005534), and aligned to the human reference genome (GRCh38) by Burrows-Wheeler Aligner v0.7.15a (RRID:SCR_010910; ref. 25). PCR duplicates were removed with MarkDuplicates (Picard v 2.17.11, http://broadinstitute.github.io/picard/, RRID:SCR_006525) algorithm. Indelrealigner and RealignerTargetCreator kits of GATK (v3.4; ref. 26) were used to align indel regions. Mutect (v.1.1.4; RRID:SCR_000559; ref. 27) was used to identify somatic single-nucleotide variants (SNV). We used IndelGenotyper (36.3336) of GATK (v3.4; RRID:SCR_001876) for somatic indel calling. We applied the HaplotypeCaller algorithm of GATK (26) to call high quality germline variants with default parameters. To control for the false positive rate of germline variants calling, we filtered low-quality variants with the following criteria: DP < 4, QD < 2.0, FS > 60.0, MQ < 35.0, MQRankSum < −12.5, and ReadPosRankSum < −8.0. TMB was calculated as the total number of exonic somatic mutation divided by total length of exome capture probes (34 MB).
The functional impact of germline missense variants was predicted using MetaSVM (28) and annotation from the ClinVar database (RRID:SCR_006169; ref. 29). We considered a missense variant as high functional impact if it was classified as deleterious by MetaSVM or Pathogenic/Likely-Pathogenic in ClinVar. Loss-of-function (LoF) variants including frameshift indels, stop gain, and stop loss variants were also considered as high functional impact, as well as variants annotated as high-confidence loss of function in gnomAD (30). We used 723 cancer census genes from the Catalogue Of Somatic Mutations In Cancer database (COSMIC; release v94, May 28th, 2021; RRID:SCR_002260; ref. 31) to generate the oncoplots for both germline and somatic mutations. Associations between pCR and gene- or pathway-level germline variants and somatic mutations were assessed using logistic regression. For the gene-level analysis, only genes affected in at least five out of 57 cases were considered for the association test.
To assess mutations at the pathway level, we collected 107 canonical biological pathways from the NanoString Hallmarks, NanoString Metabolic, and MSigDB Pathway databases (Supplementary Table S2) and considered a pathway mutated if it had ≥1 member gene with mutation. To assess significance of pathway level mutation, we first calculated ORs of the response category (pCR or RD) versus the gene or pathway status (mutated versus wild-type) using logistic regression, next we randomly permuted the pCR or RD labels for 1,000 iterations and the OR for each gene or pathway was recalculated. The proportion of random permutations showing an OR greater than the OR of the unperturbed data was defined as the P value.
PD-L1 IHC and sTILs assessment
sTILs were assessed on FFPE hematoxylin and eosin–stained 4-μm sections. The slides were digitally scanned and independently scored by two pathologists. The sTILs score was calculated as the area occupied by mononuclear inflammatory cells over the total intratumoral stromal area (32). Immune-high cancers were defined as sTILs ≥ 30% (33). PD-L1 protein expression was assessed with chromogenic IHC using the VENTANA PD-L1 (SP263) Assay following the manufacturer's instructions. PD-L1 positivity was defined as ≥1% tumor and/or immune cells staining positive (32,34). Mann–Whitney test was used for sTILs scores and the Fisher exact test was use for PD-L1 IHC positivity to determine if there were significant (P < 0.05) differences between pCR and RD.
Data and material availability
All data associated with this study are presented in this paper or Supplemental Materials. The whole exome and transcriptomic data from the Yale clinical trial (NCT02489448) are deposited in National Center for Biotechnology Information (NCBI) database of Genotypes and Phenotypes (dbGaP) under bioproject number PRJNA558949. To access the GeparNuevo (NCT02685059) dataset please refer to https://gbg.de/en/research/trafo.php.
Results
Differentially expressed genes and enriched pathways between bTNBC with pCR and RD
One hundred forty-three and 66 genes were significantly overexpressed in cancers that achieved pCR and RD, respectively (Fig. 1A; Supplementary Table S3). GSEA showed that adaptive immunity (P < 0.001), cancer driver genes (P < 0.01), cell cycle and apoptosis (P < 0.05), DNA repair (P < 0.05), humoral immunity (P < 0.001), innate immunity (P < 0.001), and JAK–STAT pathways (P < 0.001) were enriched in patients with pCR (Fig. 1B; Supplementary Table S4). Epithelial–mesenchymal transition (EMT; P < 0.05), extracellular matrix (ECM; P < 0.01), and TGFβ (P < 0.05) pathways were enriched in patients with RD (Fig. 1B; Supplementary Table S4). The leading-edge genes from the enriched pathways showed that genes that regulate T-cell and B-cell activities drove the pathway enrichment in pCR, and genes that impacted macrophages, fibroblasts, and cancer cell response to cytotoxic therapy (i.e., decreased DNA repair machinery) drove pathway enrichment in RD (Fig. 1C; Supplementary Table S4). These findings confirm that high levels of immune gene expression are characteristics of highly chemotherapy-sensitive cancers. However, a subset of immune-rich bTNBC fail to achieve pCR and what drives this difference remains unknown.
Differentially expressed genes and pathways between immune-rich TNBC with pCR and RD
We also performed differential gene and pathway expression analysis between cases with pCR versus RD restricted to immune-rich cancers only. To assess how the results might depend on the definition of immune richness, we applied two histology based methods including sTILs ≥ 30% and PD-L1 IHC positivity, and five immune gene signatures dichotomized at the median (TIS > 0.1249, GeparSixto gene signature > 0.26, NHI 5-gene score > 0.0014, STAT1 score > 0.07695, and IFNγ gene expression ≥ 0.20) to define immune-rich status (Supplementary Tables S1–S2 and S5; refs. 5, 8). The ranking of differentially expressed genes varied substantially depending on how immune-rich status was defined (Fig. 2,A,–G; Supplementary Tables S6–S13). However, GSEA revealed highly consistent differences at pathway level between cancers with RD versus pCR, regardless of how immune-rich status was defined (Fig. 3A; Supplementary Tables S14–S20). The pathways that were significantly enriched in RD despite high immune infiltration included inflammation (P < 0.05) and innate immunity (P < 0.05), the TGFβ pathway and EMT were also consistently enriched but failed to reach statistical significance. In cancers with pCR, adaptive immunity (P < 0.05) and cancer driver gene pathways (P < 0.01) were significantly enriched; several other pathways including DNA repair, cell cycle, apoptosis, chromatin modifications, PIK3A, and RAS were also consistently enriched but failed to reach statistical significance. Examination of the leading-edge genes from the significantly enriched pathways revealed a substantially different cytokine/chemokine milieu with RD compared with pCR (Fig. 3B; Supplementary Tables S14–S20). In cancers with RD, the dominant chemokines were CCL-3, -4, -5, -8, -23, CXCL-1, -3, -6, -10, and cytokines were IL1A/B, -23A, -27, -34. The chemokines CCL-3, -5 and CXCL-1, -6 are major chemoattractants for tumor-associated macrophages (TAM) and neutrophils that can exerting protumorigenic effects (35). IL34 promotes macrophage differentiation and IL1 is the classical macrophage-derived proinflammatory cytokine. Indeed, the leading-edge genes of the innate immunity pathway enriched in RD included all the hallmarks of a strong macrophage presence; high expression of CSF1, CSF1R, CD14, scavenger receptor MARCO, and Toll-like receptors (TLR) 1, -2, -3, -4, -5, -6.
In immune-rich cancers with pCR, the dominant cytokines were IFNγ, IL2, -12A/B, -21, and chemokines CXCL-9 and CXCL-13 and its receptor CXCR5. IFNγ and IL2 are the quintessential immune growth factors that play critical roles in activating and sustaining T-cell response. IL12 induces Th cell differentiation and increases the cytotoxic activity of T cells; it also inhibits TAMs and myeloid-derived suppressor cells. IL21 regulates differentiation of B cells into plasma cells and increases cytotoxicity of T cells. CXCL-9 is a chemoattractant for activated T-cells, and CXCL-13 is a chemoattractant for B cells. Consistent with this highly immune activating cytokine milieu, the leading-edge genes also included many T-cell (CD3, CD5, CD6, CD7, CD40LG) and B-cell (MS4A1, CD19, CD38, CD22, CD37, CD79A) markers, human leukocyte antigen class II (HLA-D) molecules that present antigens to T cells, and granzymes that mediate cytotoxicity.
When we compared immune-poor cancers with pCR versus RD we identified different sets of differentially expressed genes, with less than 40% overlap with the genes associated with response in immune-rich cancers, suggesting that different processes are involved in determining response or resistance depending on the immune microenvironment. Cancers with pCR were enriched in the adaptive (P < 0.01), humoral (P < 0.05), and innate immunity pathways (P < 0.05) despite belonging to the overall immune-poor subset. Patients with RD were enriched in the angiogenesis (P < 0.01), ECM (P < 0.05), and RAS pathways (P < 0.05; Supplementary Tables S21–S34; Supplementary Fig. S4).
The GeparNuevo randomized trial used an essentially identical durvalumab plus chemotherapy arm as our study and therefore represents an ideal validation cohort that also provides an opportunity to test the immunotherapy-specific predictive role of our response associated genes. We tested if the leading-edge genes that distinguished immune-rich cancers with pCR from those with RD in the Yale cohort were also differentially expressed between pCR and RD in immune-rich cancers from the GeparNuevo trial. Only 36 of the leading-edge genes associated with pCR in immune-rich bTNBC had expression data available from the GeparNuevo samples. Supplementary Table S35 lists the gene level validation results. Most importantly, we observed that from our gene list IFNG and IL21 were significantly positively associated with pCR in the chemotherapy alone arm and CXCL9, CXCL13, CD79A, and cytotoxins GZMA and GZMB were positively associated with pCR only in the durvalumab arm. Chemokines CXCL1 and CXCL3 were positively associated with RD in chemotherapy alone arm whereas CSF1, TLR3, CCL5, CXCL10, and CCL4 were associated with RD in durvalumab arm only. An immune-rich pCR signature created from the mean value of the scaled expression of IFNG, IL2, IL21, CD79A, and GZMB that individually showed a weak association with pCR (P < 0.2) in GeparNuevo, showed significantly higher expression in cases with pCR in the durvalumab arm (P = 0.040) but not in the placebo arm (P = 0.923) or in immune-poor cancers irrespective of treatment (Fig. 4).
Germline and somatic mutation landscape associated with response
We identified 206 protein coding genes with high functional impact germline variants affecting at least 5 (out of 57) patients. Among genes affected by germline variants, there were no significant differences in variant frequency by pathologic response after adjusting for multiple comparison (Fig. 5A; Supplementary Tables S36). Eight patients had germline BRCA1/2 mutation (5 pCR, 3 RD; P = 0.2). Somatic mutations affected 3,422 distinct genes, among these 1342 were mutated in only one cancer (Fig. 5B; Supplementary Tables S37). The most frequently mutated gene was TP53 (24 pCR, 22 RD). There was no statistically significant difference in somatic mutation frequencies by pathologic response for any gene after adjustment for multiple comparison. Next, we assessed associations between response and pathway level mutations separately for somatic mutations and high functional impact germline variants. Four pathways were significantly enriched in high functional impact germline variants including PI3K, DNA damage repair, MAPK, and WNT/β-catenin signaling pathways (P < 0.05; Table 1; Supplementary Table S38). These pathways were more frequently affected in cases with pCR. Somatic mutations were enriched in 22 pathways (P < 0.05) including the same four pathways identified in the germline analysis (Table 1; Supplementary Table S39). Higher TMB was significantly associated with pCR and was independent of immune gene signature expression (Table 2; Supplementary Table S40).
Pathway . | Source . | Mutation . | pCR (n = 26) . | RD (n = 31) . | Mutation rate in pCR . | Mutation rate in RD . | Permutation P value . | OR . | OR lower than 95% CI . |
---|---|---|---|---|---|---|---|---|---|
PI3K | Nanostring metabolic | Germline | 16 | 9 | 0.6154 | 0.2903 | 0.0082 | 3.9111 | 2.2232 |
DNA damage repair | Nanostring metabolic | Germline | 10 | 5 | 0.3846 | 0.1613 | 0.0307 | 3.2500 | 1.7254 |
MAPK | Nanostring metabolic | Germline | 8 | 4 | 0.3077 | 0.1290 | 0.0435 | 3.0000 | 1.5141 |
WNT_BETA_CATENIN_SIGNALING | MSigDB Hallmark | Germline | 16 | 13 | 0.6154 | 0.4194 | 0.0493 | 2.2154 | 1.2870 |
INTERFERON_ALPHA_RESPONSE | MSigDB Hallmark | Somatic | 12 | 5 | 0.4615 | 0.1613 | 0.0048 | 4.4571 | 2.3808 |
MAPK | Nanostring metabolic | Somatic | 25 | 23 | 0.9615 | 0.7419 | 0.0076 | 8.6957 | 2.8966 |
Myc | Nanostring metabolic | Somatic | 24 | 22 | 0.9231 | 0.7097 | 0.0089 | 4.9091 | 2.1290 |
DNA damage repair | Nanostring metabolic | Somatic | 25 | 23 | 0.9615 | 0.7419 | 0.0131 | 8.6957 | 2.8966 |
Transcriptional regulation | Nanostring metabolic | Somatic | 25 | 23 | 0.9615 | 0.7419 | 0.0175 | 8.6957 | 2.8966 |
Wnt_pathway | Nanostring hallmarks | Somatic | 24 | 23 | 0.9231 | 0.7419 | 0.0184 | 4.1739 | 1.7971 |
Cell cycle | Nanostring metabolic | Somatic | 24 | 23 | 0.9231 | 0.7419 | 0.0184 | 4.1739 | 1.7971 |
Transcriptional_misregulation | Nanostring hallmarks | Somatic | 24 | 23 | 0.9231 | 0.7419 | 0.0217 | 4.1739 | 1.7971 |
WNT_BETA_CATENIN_SIGNALING | MSigDB Hallmark | Somatic | 24 | 23 | 0.9231 | 0.7419 | 0.0217 | 4.1739 | 1.7971 |
MYC_TARGETS_V2 | MSigDB Hallmark | Somatic | 7 | 3 | 0.2692 | 0.0968 | 0.0250 | 3.4386 | 1.6221 |
DNA_REPAIR | MSigDB Hallmark | Somatic | 25 | 25 | 0.9615 | 0.8065 | 0.0268 | 6.0000 | 1.9645 |
HYPOXIA | MSigDB Hallmark | Somatic | 16 | 12 | 0.6154 | 0.3871 | 0.0287 | 2.5333 | 1.4670 |
INTERFERON_GAMMA_RESPONSE | MSigDB Hallmark | Somatic | 15 | 10 | 0.5769 | 0.3226 | 0.0313 | 2.8636 | 1.6481 |
UV_RESPONSE_UP | MSigDB Hallmark | Somatic | 12 | 8 | 0.4615 | 0.2581 | 0.0320 | 2.4643 | 1.3956 |
APICAL_SURFACE | MSigDB Hallmark | Somatic | 8 | 4 | 0.3077 | 0.1290 | 0.0333 | 3.0000 | 1.5141 |
Cell_cycle_and_apoptosis | Nanostring hallmarks | Somatic | 25 | 24 | 0.9615 | 0.7742 | 0.0344 | 7.2917 | 2.4113 |
Cytokine & chemokine signaling | Nanostring metabolic | Somatic | 25 | 24 | 0.9615 | 0.7742 | 0.0344 | 7.2917 | 2.4113 |
E2F_TARGETS | MSigDB Hallmark | Somatic | 25 | 24 | 0.9615 | 0.7742 | 0.0344 | 7.2917 | 2.4113 |
PI3K | Nanostring metabolic | Somatic | 25 | 24 | 0.9615 | 0.7742 | 0.0360 | 7.2917 | 2.4113 |
Pathway . | Source . | Mutation . | pCR (n = 26) . | RD (n = 31) . | Mutation rate in pCR . | Mutation rate in RD . | Permutation P value . | OR . | OR lower than 95% CI . |
---|---|---|---|---|---|---|---|---|---|
PI3K | Nanostring metabolic | Germline | 16 | 9 | 0.6154 | 0.2903 | 0.0082 | 3.9111 | 2.2232 |
DNA damage repair | Nanostring metabolic | Germline | 10 | 5 | 0.3846 | 0.1613 | 0.0307 | 3.2500 | 1.7254 |
MAPK | Nanostring metabolic | Germline | 8 | 4 | 0.3077 | 0.1290 | 0.0435 | 3.0000 | 1.5141 |
WNT_BETA_CATENIN_SIGNALING | MSigDB Hallmark | Germline | 16 | 13 | 0.6154 | 0.4194 | 0.0493 | 2.2154 | 1.2870 |
INTERFERON_ALPHA_RESPONSE | MSigDB Hallmark | Somatic | 12 | 5 | 0.4615 | 0.1613 | 0.0048 | 4.4571 | 2.3808 |
MAPK | Nanostring metabolic | Somatic | 25 | 23 | 0.9615 | 0.7419 | 0.0076 | 8.6957 | 2.8966 |
Myc | Nanostring metabolic | Somatic | 24 | 22 | 0.9231 | 0.7097 | 0.0089 | 4.9091 | 2.1290 |
DNA damage repair | Nanostring metabolic | Somatic | 25 | 23 | 0.9615 | 0.7419 | 0.0131 | 8.6957 | 2.8966 |
Transcriptional regulation | Nanostring metabolic | Somatic | 25 | 23 | 0.9615 | 0.7419 | 0.0175 | 8.6957 | 2.8966 |
Wnt_pathway | Nanostring hallmarks | Somatic | 24 | 23 | 0.9231 | 0.7419 | 0.0184 | 4.1739 | 1.7971 |
Cell cycle | Nanostring metabolic | Somatic | 24 | 23 | 0.9231 | 0.7419 | 0.0184 | 4.1739 | 1.7971 |
Transcriptional_misregulation | Nanostring hallmarks | Somatic | 24 | 23 | 0.9231 | 0.7419 | 0.0217 | 4.1739 | 1.7971 |
WNT_BETA_CATENIN_SIGNALING | MSigDB Hallmark | Somatic | 24 | 23 | 0.9231 | 0.7419 | 0.0217 | 4.1739 | 1.7971 |
MYC_TARGETS_V2 | MSigDB Hallmark | Somatic | 7 | 3 | 0.2692 | 0.0968 | 0.0250 | 3.4386 | 1.6221 |
DNA_REPAIR | MSigDB Hallmark | Somatic | 25 | 25 | 0.9615 | 0.8065 | 0.0268 | 6.0000 | 1.9645 |
HYPOXIA | MSigDB Hallmark | Somatic | 16 | 12 | 0.6154 | 0.3871 | 0.0287 | 2.5333 | 1.4670 |
INTERFERON_GAMMA_RESPONSE | MSigDB Hallmark | Somatic | 15 | 10 | 0.5769 | 0.3226 | 0.0313 | 2.8636 | 1.6481 |
UV_RESPONSE_UP | MSigDB Hallmark | Somatic | 12 | 8 | 0.4615 | 0.2581 | 0.0320 | 2.4643 | 1.3956 |
APICAL_SURFACE | MSigDB Hallmark | Somatic | 8 | 4 | 0.3077 | 0.1290 | 0.0333 | 3.0000 | 1.5141 |
Cell_cycle_and_apoptosis | Nanostring hallmarks | Somatic | 25 | 24 | 0.9615 | 0.7742 | 0.0344 | 7.2917 | 2.4113 |
Cytokine & chemokine signaling | Nanostring metabolic | Somatic | 25 | 24 | 0.9615 | 0.7742 | 0.0344 | 7.2917 | 2.4113 |
E2F_TARGETS | MSigDB Hallmark | Somatic | 25 | 24 | 0.9615 | 0.7742 | 0.0344 | 7.2917 | 2.4113 |
PI3K | Nanostring metabolic | Somatic | 25 | 24 | 0.9615 | 0.7742 | 0.0360 | 7.2917 | 2.4113 |
Note: Pathways with permutation P < 0.05 are shown.
Abbreviation: CI, confidence interval.
Modela . | Multivariate . |
---|---|
pCR–TMB | OR 1.62 (1.09–2.61) |
P = 0.0279 | |
pCR–GeparSixto | OR 2.86 (1.37–6.80) |
P = 0.0091 | |
pCR–TIS | OR 3.06 (1.44–7.56) |
P = 0.0073 | |
pCR–TMB+GeparSixto | OR 1.83 (1.16–3.29) |
P = 0.0213 | |
pCR–TMB+TIS | OR 1.81 (1.15–3.28) |
P = 0.0249 |
Modela . | Multivariate . |
---|---|
pCR–TMB | OR 1.62 (1.09–2.61) |
P = 0.0279 | |
pCR–GeparSixto | OR 2.86 (1.37–6.80) |
P = 0.0091 | |
pCR–TIS | OR 3.06 (1.44–7.56) |
P = 0.0073 | |
pCR–TMB+GeparSixto | OR 1.83 (1.16–3.29) |
P = 0.0213 | |
pCR–TMB+TIS | OR 1.81 (1.15–3.28) |
P = 0.0249 |
aModel + age + T size + N status.
Discussion
Low levels of TILs and low expression of a broad range of immune genes were associated with lack of pCR after chemotherapy plus durvalumab therapy. The lower immune infiltration was accompanied by higher expression of TGFβ and mesenchymal features of the cancer. TGFβ is an important negative regulator of cellular immunity and has been implicated in immune evasion and resistance to PD-L1 blockade in multiple cancer models (36–38). These observations suggest that targeting TGFβ could remove a barrier to immune infiltration and create a more immune competent tumor microenvironment in immune-cold cancers.
While immune-rich cancers more frequently achieve pCR, a substantial minority continues to have viable residual cancer at surgery. We, therefore, examined transcriptional differences between immune-rich bTNBC that had pCR versus those with RD. Immune-rich cancers that achieved pCR were characterized by activated T cells and B cells, high expression of immunoglobulins, granzymes, granulysin, and HLA class II antigens. The dominant cytokines in the tumor microenvironment were IFNγ, IL2, -12, -21, CXCL-9, -13, and CXCR5. These are classic chemoattractants and activators of T cells, B cells, and mediators of adaptive immunity (39). In contrast, immune-rich bTNBCs with RD were enriched in genes associated with myeloid/macrophage activity including monocyte chemoattractants CCL5 and CXCL10, ILs that these cells secrete (IL1, -23, -27, -34), and classic TLRs (TLR-1, -2, -3, -4, -5, -6) that provide pathogen recognition and subsequent activation of innate immunity. RD samples were also enriched in genes inhibiting the complement system including CD46 and CD55. An impaired complement system can hinder both antibody-dependent and -independent cell death and diminish macrophage-mediated phagocytosis. These results suggest that some immune-rich cancers have a possibly dysfunctional innate, rather than adaptive immune response to the cancer that makes these cancers less responsive to cytotoxic and PD-L1–directed therapies. Altering the cytokine environment by targeting IL1 (40) and the macrophage monocyte lineage might alter the balance between a dysfunctional innate and more effective adaptive immune response in otherwise immune-rich TNBC.
In a single-arm anti–PD-L1 plus chemotherapy trial it is not possible to determine which response marker, if any, is selectively predictive of benefit from the combination versus individual components. Multiple studies have demonstrated that higher immune infiltration that can be captured by a large number of different immune gene signatures due to their highly correlated coexpression is associated with higher pCR rate to chemotherapy with or without immune checkpoint inhibitors. Our results indicated that there are significant differences in the cytokine and immune milieu of immune-rich cancers that achieved pCR with durvalumab plus chemotherapy versus those that did not. We therefore tested if the leading-edge genes that distinguished immune-rich cancers with pCR from those with RD in the Yale cohort were also overrepresented in cancers with pCR in immune-rich cancers from the GeparNuevo trial, and if any of the genes could predict benefit selectively from durvalumab. We could only perform partial validation due to many missing genes in the GeparNuevo data, but reassuringly several genes showed a similar trend as seen in the Yale cohort and several cytokines showed a differential predictive role by treatment arm. We created a five-gene signature from genes individually weekly associated with pCR that showed significantly higher expression in immune-rich TNBC with pCR in the durvalumab arm (P = 0.040) but not in the placebo arm (P = 0.923) or in immune-poor cancers irrespective of treatment.
Our study has limitations, as we could only partially validate our observations in the similar GeparNuevo trial due to missing information on many candidate genes. The GeparNuevo data was generated on a different RNA-seq platform and represent targeted sequencing with a different dynamic range than whole transcriptome RNA-seq. These differences decrease the power and accuracy of our validation attempt. We also recognize that the candidate genes were identified in the Yale cohort but the five-gene durvalumab plus chemotherapy predictive signature itself was selected from the GeparNuevo data, and therefore further validation on independent data will be required. Due to colinear expression of many immune genes, it is entirely possible that other genes could also provide the same, or even better, response discriminating function. However, the combined analysis of these two trials suggests that there are immunologic differences between immune-rich TNBC that achieve pCR and those that do not. These differences can inspire new therapeutic strategies and may hold the key for developing new biomarkers for treatment selection.
We also examined associations between pathologic response and germline variants in coding genes and somatic mutations. We found no statistically significant differences in somatic mutation or germline variant frequencies by pathologic response for any gene. However, TMB was significantly higher in cancers with pCR. When mutations were mapped to biological pathways, we found that cancers with pCR had significantly more germline variants and somatic mutations in the PI3K, DNA damage repair, MAPK, and WNT/β-catenin pathways. There were no statistically significantly more frequently mutated pathways in cases with RD. The more frequent mutations in cancer relevant signaling pathways and DNA repair genes might lead to a more immunogenic cancer, however we found no positive correlation between TMB and immune gene expression, similar to an earlier study (41). It is more likely that mutations in these cancer relevant pathways directly lead to increased chemotherapy sensitivity due to DNA damage repair deficiency (42–44).
In conclusion, genes in the TGFβ pathway are associated with immune-attenuated phenotype and lack of pCR. Among immune-rich cancers that fail to achieve pCR, macrophage/neutrophil and innate immunity related chemoattractants dominate the cytokine milieu, whereas in cancers with pCR, IFNγ, and activated B and T cells and adoptive immunity-related markers dominate the tumor microenvironment. Inhibitors of complement cascade blockers, TGFβ inhibitors, and modulators of TAMs may improve immunotherapy efficacy in basal-like TNBC.
Authors' Disclosures
K.R.M. Blenman reports scientific advisory board membership at CDI Labs. T. Karn reports a patent for EP18209672 pending. A. Silber reports personal fees from AstraZeneca during the conduct of the study. C. Denkert reports personal fees from MSD Oncology, Daiichi Sankyo, Molecular Health, AstraZeneca, Merck, Roche, and Eli Lilly and Company; and grants from Myriad and German Breast Group outside the submitted work; in addition, C. Denkert has a patent for WO2020109570A1 pending. B.V. Sinn reports a patent for PCT/EP2019/083124 pending. M. Rozenblit reports grants from American Society of Clinical Oncology (ASCO) Young Investigator Award outside the submitted work. D.L. Rimm reports grants and personal fees from Amgen, AstraZeneca, Cepheid, Konica - Minolta, Eli Lilly and Company, and NextCure; personal fees from Cell Signaling Technology, Danaher, Fluidigm, GSK, Merck, Monopteros, NanoString, Odonate, Paige.AI, Regeneron, Roche, Sanofi-Aventis, and Ventana; and personal fees from Verily outside the submitted work; in addition, D.L. Rimm has a patent for Rarecyte with royalties paid. S. Loibl reports grants and other support from Abbvie, AstraZeneca, and Celgene; other support from Amgen, Bayer HealthCare, BMS, Eirgenix, GSK, Eli Lilly and Company, Merck, Pierre Fabre, Prime/Medscape, and Samsung; grants, nonfinancial support, and other support from Daiichi Sankyo, Gilead, Novartis, Pfizer, and Roche; and non-financial support and other support from Puma and Seagen outside the submitted work; in addition, S. Loibl has a patent for EP14153692.0 pending, a patent for EP21152186.9 pending, a patent for EP15702464.7 issued, a patent for EP19808852.8 pending, and a patent for Digital Ki67 Evaluator issued and with royalties paid. L. Pusztai reports grants and personal fees from AstraZeneca during the conduct of the study; personal fees from Pfizer, Merck, Novartis, Bristol-Myers Squibb, Genentech, Seagen, Syndax, Personalis, and Natera; and other support from Bristol-Myers Squibb, Pfizer, Seagen, and Merck outside the submitted work. No disclosures were reported by the other authors.
Authors' Contributions
K.R.M Blenman: Conceptualization, data curation, formal analysis, supervision, visualization, methodology, writing–original draft, writing–review and editing. M. Marczyk: Formal analysis, visualization. T. Karn: Validation, writing–review and editing. T. Qing: Formal analysis, visualization, writing–review and editing. X. Li: Formal analysis, visualization, writing–original draft, writing–review and editing. V. Gunasekharan: Formal analysis, writing–review and editing. V. Yaghoobi: Methodology, writing–review and editing. Y. Bai: Methodology, writing–review and editing. E.Y. Ibrahim: Data curation, writing–review and editing. T. Park: Data curation, writing–review and editing. A. Silber: Data curation, writing–review and editing. D.M. Wolf: Formal analysis, writing–review and editing. E. Reisenbichler: Methodology, writing–review and editing. C. Denkert: Validation, writing–review and editing. B.V. Sinn: Validation, writing–review and editing. M. Rozenblit: Data curation, writing–review and editing. J. Foldi: Data curation, writing–review and editing. D.L. Rimm: Methodology, writing–review and editing. S. Loibl: Validation, writing–review and editing. L. Pusztai: Conceptualization, resources, data curation, supervision, funding acquisition, investigation, writing–original draft, project administration, writing–review and editing.
Acknowledgments
This work was supported by an NCI grant (R01CA219647, to L. Pusztai), a Susan Komen Foundation Leadership Award (SAC160076, to L. Pusztai), investigator awards from the Breast Cancer Research Foundation (AWDR11559, to L. Pusztai and D. Rimm), and research grant M82 from H.W. & J. Hector-Foundation (to T. Karn).