Comprehensive Analysis of the Immunogenomics of Triple-Negative Breast Cancer Brain Metastases From LCCC1419.

Background: Triple negative breast cancer (TNBC) is an aggressive variant of breast cancer that lacks the expression of estrogen and progesterone receptors (ER and PR) and HER2. Nearly 50% of patients with advanced TNBC will develop brain metastases (BrM), commonly with progressive extracranial disease. Immunotherapy has shown promise in the treatment of advanced TNBC; however, the immune contexture of BrM remains largely unknown. We conducted a comprehensive analysis of TNBC BrM and matched primary tumors to characterize the genomic and immune landscape of TNBC BrM to inform the development of immunotherapy strategies in this aggressive disease.
Methods: Whole-exome sequencing (WES) and RNA sequencing were conducted on formalin-fixed, paraffin-embedded samples of BrM and primary tumors of patients with clinical TNBC (n = 25, n = 9 matched pairs) from the LCCC1419 biobank at UNC-Chapel Hill. Matched blood was analyzed by DNA sequencing as a comparison for tumor WES for the identification of somatic variants. A comprehensive genomics assessment, including mutational and copy number alteration analyses, neoantigen prediction, and transcriptomic analysis of the tumor immune microenvironment were performed.
Results: Primary and BrM tissues were confirmed as TNBC (23/25 primaries, 16/17 BrM) by immunohistochemistry and of the basal intrinsic subtype (13/15 primaries and 16/19 BrM) by PAM50. Compared to primary tumors, BrM demonstrated a higher tumor mutational burden. TP53 was the most frequently mutated gene and was altered in 50% of the samples. Neoantigen prediction showed elevated cancer testis antigen- and endogenous retrovirus-derived MHC class I-binding peptides in both primary tumors and BrM and predicted that single-nucleotide variant (SNV)-derived peptides were significantly higher in BrM. BrM demonstrated a reduced immune gene signature expression, although a signature associated with fibroblast-associated wound healing was elevated in BrM. Metrics of T and B cell receptor diversity were also reduced in BrM. Conclusions: BrM harbored higher mutational burden and SNV-derived neoantigen expression along with reduced immune gene signature expression relative to primary TNBC. Immune signatures correlated with improved survival, including T cell signatures. Further research will expand these findings to other breast cancer subtypes in the same biobank. Exploration of immunomodulatory approaches including vaccine applications and immune checkpoint inhibition to enhance anti-tumor immunity in TNBC BrM is warranted.

Introduction

Triple-negative breast cancer (TNBC) lacks the expression of hormone receptors estrogen (ER) and progesterone (PR) as well as human epidermal growth factor receptor 2 (HER2). TNBC is also the most aggressive subtype of breast cancer, with a predilection for brain metastases; up to 50% of patients with metastatic TNBC will develop brain metastases (BrM) during their disease course (1). Patients with TNBC BrM face a poor prognosis with a median survival following diagnosis of less than 6 months (2). Despite progress in the treatment of ER+ and HER2+ breast cancer BrM with the advent of brain-penetrant targeted therapies, the prognosis for TNBC BrM remains largely unchanged over the past decade (3). Thus, studies to better understand the biology of TNBC BrM to identify new therapeutic targets are needed.

Previous studies have conducted sequencing of primary and metastatic tissues, including BrM, from melanoma (4), lung cancer (5), breast cancer (6–8), and multiple solid tumor types (9). A seminal work by Brastianos et al. demonstrated that some solid tumors (including BrM) undergo branched evolution during the metastatic process. These studies have led to a growing appreciation that BrM can be biologically different from not just their primary tumors but also extracranial metastases, including differential acquisition or loss of targetable alterations. Breast cancer brain metastasis (BCBrM) have demonstrated mutations and/or copy number alterations in clinically targetable genes and pathways such as HER2 (6, 7, 9), BRAF (8), PI3K/Akt (9), CDK (6), ATM (8), and CRYAB (10) not seen in primary tumors. BrM can also be metabolically different from primaries, with increased oxidative phosphorylation (4). Preclinical studies have demonstrated that these targets can alter the brain metastatic potential and/or growth of BrM in breast cancer models (4, 11–13). These findings have led to the first genomically guided clinical trial in BrM to match alterations present in BrM to an appropriate brain-penetrant inhibitor (NCT03994796).

While these and other prior studies have made significant progress in the genomic characterization of BrM in recent years, no studies have yet focused exclusively on TNBC BrM specifically, and few studies have looked at comprehensive RNA and DNA sequencing-derived features of the tumor immune microenvironment. Using the LCCC1419 Biobank of metastatic breast cancer samples, we have collected and analyzed BrM, matched primaries, and normal tissue from 25 patients with TNBC through both whole-exome sequencing (WES) and mRNA sequencing. We report the somatic mutational landscape of TNBC BrM compared to primary tumors and implement a comprehensive neoantigen prediction pipeline to elucidate potentially immunogenic peptides arising from traditional and alternative neoantigen sources. Utilizing mRNA gene expression, we evaluated the tumor immune microenvironment of TNBC BrM relative to primary tumors and correlate these features with overall survival. To our knowledge, this study represents the largest evaluation of TNBC BrM through WES and RNA-seq and is the first to analyze gene expression and immunogenomics in addition to the mutational landscape.

Materials and Methods

Patient Consent and Tissue Collection

Archival formalin-fixed, paraffin embedded (FFPE) tumor tissues were obtained from patients with clinically determined triple-negative breast cancer (TNBC) based on either a primary or metastatic site, with known metastasis to the brain. The patients consented to participation in either the UNC Health Registry (UNC IRB 09-0605), opened on 04/16/2009 and consented between 11/2014 and 06/2016, or to a clinically annotated biobank study at the University of North Carolina at Chapel Hill under an Institutional Review Board (IRB)-approved protocol (LCCC1419) which opened on 10/31/2014 and consented from 11/2014 to 11/2018. Brain metastases tissues were available from n = 19 patients (n = 19 with RNA, n = 17 with DNA), while primary tumor tissue was available from n = 16 patients (n = 15 with RNA, n = 13 with DNA). Matched whole blood samples were available for n = 22 patients (DNA). Of these cases, n = 9 included matched RNA primary TNBC and BrM tissue pairs from the same patient, and n = 6 had matched DNA triplet samples (primary, BrM, and blood), all n = 6 of which also had RNA for both primary and BrM samples.

DNA Whole-Exome Sequencing and Variant Calling

FFPE tumor/tissue blocks and normal fresh frozen blood samples were collected and inventoried through honest brokers in accordance with IRB standards through the UNC Health Registry. UNC patient samples were inventoried through UNC Surgical Pathology Core (SP), while patient samples from outside UNC were inventoried through UNC Tissue Procurement Facility. Tissue blocks were sectioned by UNC Translational Pathology Laboratory (TPL). Twenty-two sections were made at a time: two 5-μm sections at the beginning and end of sectioning for pathologist review and 20 10-μm sections for DNA/RNA isolation on glass slides. TPL pathologists reviewed the samples, circling areas with more than 50% tumor mass for DNA/RNA isolation. When required, the process was repeated to collect more DNA/RNA. UNC BioSpecimen Processing Facility performed all DNA/RNA isolations. DNA from tumor-enriched cores were extracted using the Maxwell 16 FFPE Tissue LEV DNA Purification Kit. DNA was exome-captured and amplified with Agilent SureSelect XT (G9611B) and capture library (5190-8881). WES was performed using the Illumina HiSeq 2500 or NextSeq 500 platform with multiple samples per lane using 2 × 100 paired-end chemistry. RNA was isolated from the same sections with QIAGEN AllPrep FFPE, and libraries were prepared with Illumina TruSeq Stranded with RiboZero Gold (RS-122-2301). mRNA-Seq libraries were run at 2 samples per lane on an Illumina HiSeq2500 sequencer in high-output mode using 2 × 50 paired-end chemistry.

WES was performed on FFPE tumor tissue, with peripheral blood mononuclear cells serving as a matched normal. Library preparation was performed with the SureSelect XT Human All Exon V6 + UTR kit (Agilent, Santa Clara, CA, USA), and pooled samples were sequenced on the HiSeq2500 platform (Illumina). The resulting somatic and germline WES sequencing files were aligned to Hg38 using bwa (v0.7.17), sorted, and indexed, and duplicates were marked using biobambam2 (v2.0.87). BAMs were re-aligned with Abra2 (v2.22), followed by somatic and germline variant detection with Strelka2 (v2.9.10), Cadabra (from Abra2 v2.22), and Mutect2 (GATK v4.1.4.0). The capture of exonic sequences was verified using the Picard (v2.21.1) CollectHsMetrics tool, and the quality of sequencing data was verified using FastQC (v0.11.8) and the Picard suite’s CollectAlignmentSummaryMetrics, CollectInsertSizeMetrics, QualityScoreDistribution, and MeanQualityByCycle tools. Variants with matched normals were filtered by the following criteria: protein-coding mutations only, Cadabra indel quality >10.5, Mutect2 indel quality >6.8, or (single-nucleotide variant) SNV quality >9.2, Strelka2 indel quality >15.2 or SNV quality >19.7. Additionally, Cadabra indels with quality <35 required a supporting high-quality call from either Strelka2 or Mutect2, and Strelka2 calls with SomaticEVS <20 similarly required a matching call from either Mutect2 or Cadabra. Variants for tumor-only samples were detected by Mutect2 and filtered to retain the protein-coding mutations. The remaining variants required at least 5 supporting reads and a minimum read depth of 40 or 10 supporting reads and a minimum read depth of 80 if MAF <5%. The variants with a MAF >5% in normal tissue were dropped, as were the variants appearing at rates above 1% in any subpopulation in either GnomAD or 1000 Genomes databases. To counter FFPE artifacts, C>T and G>A substitutions required a minimum MAF of 10%. Tumor mutational burden (TMB) was calculated from small indels and substitutions identified by WES and divided by the megabases adequately covered by sequencing reads. Whole-exome sequencing data for both tumor and germline was available for 28 samples (representing 22 unique patients) at baseline. Tumor-only whole-exome sequencing data—without matched normal—was available for a further 2 samples (representing 2 patients), for a total of 30 samples across 24 patients with WES data. Oncoplots were created in R using maftools v2.10.0, with variant genes limited to those implicated in breast cancer from the COSMIC Cancer Gene Census Tier 1 list (https://cancer.sanger.ac.uk/cosmic).

Copy Number Variation and Subclonal Heterogeneity Analysis

LCCC1419 patients with matched DNA normal and DNA tumor samples (n = 22) were processed through trim galore v. 0.6.2. The resulting trimmed FASTQs were aligned to the human reference genome FASTA file Homo_sapiens_assembly38.fasta (from the GATK hg38 file bundle) with the bwa mem command from BWA v. 0.7.17 with default parameters. The resulting SAM files were sorted, converted to BAM files, and indexed using SAMtools v. 1.9. The matched samples from each patient were then processed through the Sequenza v. 3.0.0 workflow (gc_wiggle with window size of 50 basepairs, bam2seqz, seqz_binning, sequenza.extract, sequenza.fit, and sequenza.results all with default parameters). The resulting patient-level segment files were then modified to conform to the format required by GISTIC. This conversion was performed with standard BASH scripting and included a log2(x) - 1 transformation of Sequenza’s raw depth ratio value for GISTIC. The resulting modified segment files were then concatenated across samples and ran through GISTIC v. 2.0.23 using the hg38.UCSC.add_miR.160920.refgene.mat file (from docker://shixiangwang/gistic:1.2) for the -refgene parameter. CNVKit was run using the “batch” mode for primary tumors and BrM separately (with their matched normal samples, respectively). All normal samples within each group were pooled together to generate a pan-sample normal control. Agilent’s SureSelect Human All Exon V6+UTRs probed bed file was provided for the –targets parameter, a standard hg38 refFlat.txt file was provided for annotations, and a k50 umap mappability BED file was provided for the –access parameter. Outputs from Sequenza were paired with each patient’s corresponding MuTect2 somatic variant calls to create PyClone-VI input files using Sequenza copy number information and MuTect2’s variant read support information. Files were created for each primary tumor and BrM sample separately, and paired input files were created for patients that had both primary tumor and BrM data available. The resulting files were run through PyClone-VI using default parameters.

Neoantigen Prediction Using Genomics Data

Tumor antigens were predicted from a comprehensive set of genomic sources (single-nucleotide variations, insertions/deletions, gene fusions, alternative splice variants, cancer testis antigens, overexpressed self-antigens, and viral and endogenous retroviral antigens) using methods developed by our group and implementations of methods developed by others (14–22). Briefly, whole-exome sequencing was used to identify tumor-specific genetic variants (single-nucleotide variations, insertions/deletions, gene fusions), and RNA sequencing was used to confirm the expression of these variants. RNA sequencing data alone were used to evaluate for expressed alternative splice variants, viral and endogenous retroviral antigens. Cancer testis antigens/overexpressed self-antigens were evaluated using RNAseq data, but WES data was used to incorporate germline variants. RNA sequencing data was also used to infer tumor MHC haplotypes via HLAProfiler, the most accurate tool for MHC haplotype inference (23). Peptide fragments generated in silico are evaluated for predicted binding affinity to tumor MHC alleles using NetMHCpan-4.1 (24). Peptides with predicted binding affinity <500 nM were considered positive binders (e.g., potential tumor antigens), while peptides with predicted binding affinity <50 nM considered strong binders and more likely to be tumor antigens (25).

RNA-Seq Data Processing

RNA-Seq Paired FASTQs were run through trim galore v. 0.6.2 using –paired parameter. STAR v. 2.7.0f was used to index the reference genome Homo_sapiens.GRCh38.dna_sm.primary assembly.fa from GATK and to map trimmed reads to reference (using parameters –quantMode TranscriptomeSAM – outSAMtype BAM SortedByCoordinate –sjdbGTFfile Homo_sapiens.GRCh38.100.gtf). Gffread v. 0.11.7 was used to create a transcriptome reference using the reference genome and the gtf file Homo_sapiens.GRCh38.100.gtf. The “toTranscriptome” alignments from STAR were used with Salmon v. 1.1.0 using “salmon quant -l a”. Sample quality was assessed using MultiQC v1.9 with Picard CollectRnaSeqMetrics, and samples with less than 20 M coding reads were excluded as this threshold has been found to approximate the minimal sequencing depth to achieve equivalent detection to microarrays (26). Counts were log2-transformed and the upper quartile normalized for further downstream analysis. Some patients’ tumors were sequenced multiple times (technical replicates), and in such cases, gene-level expression values were averaged across technical replicates.

Intrinsic Subtype Analysis

Intrinsic subtype analysis was performed according to the methods described in Picornell et al (27). The R package heatmaply (28) was used for heat map visualization with hierarchical clustering based on average linkage.

Immune Gene Signature Expression Analysis

Thirty-two immune gene signatures were chosen to reflect the diversity of tumor-infiltrating immune cell populations and to minimize redundancy (refer to associated.gmt file in the ). The binfotron R package (29) was used to compute the differential gene expression [along with the DESeq2 dependency (30)], produce a volcano plot, and calculate immune signature metagene scores (median log2 expression values) for downstream analysis. ssGSEA was performed using the R packages GSVA and GSEABase (31, 32). CIBERSORTx immune cell fraction imputation using the LM22 matrix was also performed (33).

TCR/BCR Repertoire Analysis

Immune chain inference was performed on RNA-Seq samples via MiXCR 2.1.2 for TCR chains (34) and V’DJer 0.12 (35) for BCR chains. The consensus BCR contigs from V’DJer were quantified using Salmon 0.13.11 (36). Repertoire diversity was calculated using a model-based approach, which improves estimations of diversity in part by minimizing known sources of estimate bias (37).

Survival Analyses

Multivariable Cox proportional hazards modeling was performed using the survival ( 38 ) R package. The model included cancer stage, age at primary tumor diagnosis, and race. Time from initial diagnosis, from diagnosis of any metastasis regardless of anatomical location, or from diagnosis of BrM to an event was interrogated. Hazard ratios and 95% confidence intervals were returned for gene signature covariates and visualized using the forestplot ( 39) R package.

Accession Numbers and Data Sharing

Sample information for RNA-seq and DNA-seq fastQ runs, including the clinical information, were uploaded to the NCBI’s dbGaP repository (accession no. phs002457.v1.p1) and SRA.

Results

Patient Cohort Characteristics

Tissues and blood from patients with clinical TNBC (n = 25) were included in this analysis, including BrM tissues (n = 19), primary breast tumors (n = 17), and whole blood samples (n = 22). The specimen numbers by tissue type and analysis (IHC, RNA, or DNA) are outlined in . Patient demographics are included in , with individual clinical–pathological characteristics and specimen availability presented in . The majority of patients were Caucasian (n = 17, 68%), with African American, Asian, and other ethnicities represented [n = 6 (24%), n = 1 (4%), and n = 1 (4%), respectively]. One male was included in the cohort. Median age at breast cancer diagnosis was 46.7 years (range, 29–70.9), while median age at BrM diagnosis was 51.5 years (range, 31.7–72). The majority (n = 12, 48%) of patients were initially diagnosed with stage II disease prior to recurrence; a minority (n = 2, 8%) were diagnosed with de novo stage IV TNBC. In addition to BrM, other sites of disease included the liver (n = 8, 32%), bone (n = 15, 60%), lung (n = 16, 64%), and non-local lymph nodes (n = 18, 72%). Ten patients (40%) were initially diagnosed with a solitary BrM, while 7 patients were diagnosed with 5 or greater BrM (28%). BrM was supratentorial in n = 23 (92%) patients and infratentorial in n = 13 (52%) patients. The median progression-free survival (e.g., time from primary TNBC diagnosis to the diagnosis of any metastasis) was 1.8 years (range, 0–19.5). The median OS from primary TNBC diagnosis was 3.7 years (range, 0.9–19.8), while the median OS from BrM diagnosis was 1.2 years (range, 0 – 8.9).

Table 1

Relevant demographic, subtype, and clinical diagnostic information for the LCCC1419 TNBC cohort.

Characteristics	Number (%)
Demographic information
Means of enrollment (n = 25)
LCCC 1419 consent Health Registry consent Waiver of consent	6 (24%)16 (64%)3 (12%)
Race (n = 25)
African American Asian Caucasian Other	6 (24%)1 (4%)17 (68%)1 (4%)
Ethnicity (n = 25)
Hispanic Not Hispanic Unknown	2 (8%)22 (88%)1 (4%)
Sex (n = 25)
Female Male	24 (96%)1 (4%)
Smoking status (n = 25)
Never smoker Current smoker Former smoker	14 (56%)4 (1%)7 (28%)
Subtype information
Subtypes by primary resection (n = 25)
Luminal A (ER/PR+, HER2-) Luminal B (ER/PR+, HER2+) HER2 (ER-, PR-, HER2+) TNBC (ER-, PR-, HER2-) Mixed (two primaries tested with different results) Unknown	1 (4%)0 (0%)0 (0%)23 (92%)0 (0%)1 (4%)
Subtypes by CNS resection (n = 21)
Luminal A (ER/PR+, HER2-) Luminal B (ER/PR+, HER2+) HER2 (ER-, PR-, HER2+) TNBC (ER-, PR-, HER2-) Radiation necrosis Unknown	1 (5%)0 (0%)0 (0%)15 (71%)1 (5%)4 (19%)

Intrinsic Subtype Classification of Primary TNBC and Brain Metastases

Intrinsic molecular subtype analysis using RNAseq data illustrated that the majority of clinically determined TNBC samples were of the basal subtype ( ). Of the BrM (n = 19), 16 were classified as basal (84%), with normal-like (n = 2) and HER2-enriched (n = 1) comprising a small fraction of the cohort. Similarly, primary tumors (n = 15) were predominantly classified as basal (n = 13, 87%), with the remaining tumors being normal-like (n = 2). Notably, the 4 samples that were called normal-like by PAM50 analysis had a basal subtype as the second highest identity probability. There were 2 cases with discordant receptor classification between primary tumor and BrM by immunohistochemistry (n = 1 ER+/PR+/HER2- Luminal A primary converted to a TNBC BrM, and n = 1 TNBC primary converted to an ER+/PR+/HER2- Luminal A BrM) ( ; ). Despite the potential subtype switching between primary and BrM, these samples were included in the downstream analyses.

Mutational, Somatic Copy Number Alteration, and Subclonal Analyses of Primary TNBC Tumors and Brain Metastases

First, we analyzed the tumor mutational burden (TMB) of primary tumors (n = 13) relative to BrM (n = 17) using WES. On average, BrM harbored a greater mutational load than primary tumors (median 3.33 vs. 1.78 mutations/Mb, respectively, p < 0.05; ). Upon analysis of matched primary-BrM WES pairs (n = 6), however, there was no significant difference between tumor location and TMB (median 2.80 vs. 1.88 mutations/Mb, respectively, p = 0.69; ). An analysis of shared mutations within matched pairs revealed varying degrees of mutational conservation between anatomical locations ( ). The degree of variant sharing between matched pairs ( ) was generally greater than the degree of mutations shared between primaries and BrM globally ( ), highlighting interpatient tumor heterogeneity and mutational divergence. We assessed whether a survival benefit was conferred by increasing TMB, as TMB has been considered a proxy for neoantigen burden (40, 41). There was no significant association between TMB and survival (p = 0.07) from the time of primary TNBC diagnosis in the context of a model that included standard clinicopathological features (age at diagnosis of primary tumor, stage, and race) ( ). Next, we examined the mutational spectrum of genes with known associations to breast cancer development (42). We found that these genes were altered in 70% of combined primary and BrM samples, with TP53 being the most commonly mutated gene (mutated in 50% of samples, n = 15), in accordance with its known relevance to TNBC (43) ( ). The next most frequently altered genes were MAP3K13 and PIK3CA, which were mutated in 13% (n = 4) and 10% (n = 3) of samples, respectively; all other genes were mutated less frequently, occurring in 7% or less (n ≤ 2) of samples ( ).

Figure 1

Mutational analysis of the LCCC1419 cohort. (A) Increased tumor mutational burden (TMB) was observed for BrM (n = 17) relative to primary (n = 13); *p <0.05, Wilcoxon rank-sum test. (B) No difference was observed between TMB for matched pairs (n = 6). (C) Circos plot showing the total number of variants and the proportion of variants shared between matched primary/BrM (n = 6). (D) Oncoprint displaying the mutational spectrum of cancer-associated genes with known etiology to breast cancer (COSMIC Tier 1, Sanger Institute). Representation of the somatic copy number alterations in (E) primary TNBC and (F) BrM as determined by CNVkit/GISTIC 2.0 analysis. Significant amplicons (Amp) or deleted (Del) regions are annotated (q <0.25, green line). The 1q21.2 amplification and the 20p13 deletion, which are shared genomic features of primary and BrM TNBC, are highlighted by colored boxes. Relevant oncogenes and tumor suppressors are annotated on the plots, and stars indicate genomic regions where copy number alterations are known to contribute to breast cancer/aggressive basal breast cancer.

We next analyzed the recurrent copy number alteration patterns in primary and BrM samples using Sequenza/GISTIC 2.0 (44, 45). The primary TNBC samples harbored 3 significant amplicons and 1 deleted region, whereas BrM was more profoundly altered with 15 significant amplicons and 12 regions of deletion (q < 0.25; ). At this level of genomic resolution, only 2 regions were commonly altered between primary and BrM (11p13 amplicon and 13q11 deletion). In breast cancer, these two sites are previously known to be amplified or deleted, respectively (46, 47). Interestingly, a number of amplicons/deleted regions identified in this cohort are known to be associated with breast cancer/aggressive basal breast cancer, such as gains of 1q, 8p11-12, 8q, 12p13, 13q34, 17q, and 19q and deletions of 3p, 4p16.3, 8p, 11p15, 17p, and 19p13 (48–58). Documented breast cancer oncogenes (NOTCH2, ENSA, PIK3CA, CD44, WT1, BCL2L2, AKT2, and TFF3) and tumor suppressors (BRCA2 and PRKCDBP) reside or are in close proximity to some of these significantly amplified/deleted genomic regions, and these alterations potentially contribute to TNBC progression and BrM development.

Since somatic copy number alteration (SCNA) detection tools are prone to high false positive rates as well as issues with precision and accuracy (59, 60), we also performed SCNA calling using CNVkit (61), which combines all normal samples into a pooled reference to increase performance. In contrast to Sequenza, this method identified a greater number of SCNA in primary TNBC relative to BrM, showing that, as a collective group, primary TNBC samples harbored 27 significant amplicons and 17 deleted regions, whereas BrM had 8 significant amplicons and 4 regions of deletion (q < 0.25; ). Despite notable differences between the two methods, there was a concordance in the results as well, with corroboration of 8p11.22 and 11p13 amplification in primary TNBC and validation of 14q11.2 and 19q13.2 amplification and 4p16.3 deletion in BrM ( ; ). The CNVkit SCNA analysis also highlights the potential importance of oncogenes (e.g., RCP, CD44, WT1, BCL2L2, and AKT2) and tumor suppressors (e.g., PRKCDBP) to TNBC etiology and metastatic progression, as significant amplicons/deleted regions harbor these genes.

Finally, we examined the ploidy, tumor purity, and subclonal makeup of tumors in this cohort. No differences in cellular ploidy were noted between primary TNBC and BrM (median ploidy of 3.2 and 3.25, respectively; ). Similarly, no significant differences in tumor purity were observed between groups (median purity of 0.57 and 0.76, respectively; ). An analysis of subclonal tumoral architecture using Pyclone-VI (62) showed no difference between the number of clones per tumor in each location, where a median of 3.5 clones per tumor in primary TNBC (range, 1–5 clones/tumor) and a median of 4 clones per tumor in BrM (range, 3–6 clones/tumor) were observed ( ). A subclonal assessment in patients with matched primary/BrM pairs showed that some pairs had similar clonal constituency between anatomical sites (e.g., patients L-01-054, L-01-085, L-02-119, and L-03-016), whereas other pairs showed signs of divergent clonal evolution (e.g., L-02-120 and L-03-011) ( ). Interestingly, 5 of 6 matched pairs the dominant subclone harbored the highest mutational burden, which was reflected in the analysis of unmatched tumors as well (not shown), suggesting that increased mutational load may endow these subclones with a selective growth advantage.

Tumor Antigen Landscape

We next performed a comprehensive analysis of the neoantigen landscape in this cohort. We queried a range of neoantigen sources, including single-nucleotide variants (SNVs), insertion/deletion events (InDels), splice variants, structural fusion events, cancer testis antigens (CTAs)/self-antigens, endogenous retroviruses (ERVs), and viral sources excluding ERVs (63) ( ). The predominant antigen sources in both primary TNBC and BrM were CTAs/self-antigens and ERVs ( ). Upon comparison of the number of predicted neoantigen-derived peptides, there were significantly more SNV-derived MHC class I-binding peptides in BrM as compared to primary TNBC (p = 0.005), with no differences seen between groups with respect to other neoantigen sources ( ). This analysis together shows that TNBC harbors a diverse set of potentially therapeutically actionable neoantigen-derived peptides.

Figure 2

Tumor antigen sources among the LCCC1419 patients. A comprehensive bioinformatics prediction pipeline that exports (A) the number of neoantigen-derived MHC class I-binding peptides (Kd < 500 nM) broken down by antigen source was employed. Some patients’ tumors did not have associated tumor or normal whole-exome sequencing data, and as such, antigen sources that require DNA sequencing data (single-nucleotide variants, InDels, cancer testis antigens, or fusion events) are not able to be queried in these cases (denoted by †). (B) Distribution of the number of neoantigen-derived MHC class I-binding peptides (Kd < 500 nM) broken down by antigen source, corresponding to (A); **p <0.01 (Wilcoxon rank-sum test).

Comparison of Immune Gene Signatures Between Primary TNBC and BrM

Immune gene signatures (IGS) representing multiple components of the immune system, including B cells, T cells, natural killer cells, and innate immune cells along with immune cell phenotype frequencies, were evaluated between primary TNBC tumors (n = 15) and BrM (n = 19) using RNA-Seq ( ; ; ). The majority of IGS across each of these categories were lower in TNBC BrM compared to primary TNBC. A gene signature associated with fibroblast-associated wound healing [Chang_Serum_Response_Up (64)] was significantly higher in BrM relative to primary tumors (q < 0.05). RNAseq expression data from primary tumors and BrM were also assessed using CIBERSORTx (65) to determine relative frequencies of 22 immune cell subtypes (LM22) to tumor composition. In this analysis, naïve B cells and M1 macrophages were lower in BrM compared to primary tumors, while eosinophils and neutrophils were higher in BrM tissues (q < 0.05) ( ). The expression of the 20-gene immunologic constant of rejection signature [ICR (66)], which is representative of Th1-mediated immunity, cytotoxic function, and tissue-specific destruction (e.g., GVHD, autoimmunity, and allograft rejection), was also significantly reduced in BrM relative to primary TNBC ( ). Additionally, the blood transcriptional modules reported by Rinchai et al. (67) were queried against our dataset and showed a significant reduction of B and T cell modules in BrM relative to primary TNBC ( ), concordant with the abovementioned data. These results are together consistent with an overall immune-excluded brain tumor microenvironment (TME) in the context of TNBC BrM.

Figure 3

Immune gene signature metagene analysis showed an overall immune cell deficit in BrM relative to primary triple-negative breast cancer. The colored bars above the plot indicate both the immune cell category assigned to the respective signatures and whether the signatures were increased/decreased in the primary tumor relative to BrM. Wilcoxon rank-sum test was performed on Z-transformed signature scores to determine the statistical significance after false discovery rate correction. Significance codes: † q < 0.1, *q < 0.05.

T and B Cell Repertoire Analysis

We used RNA-Seq data from primary TNBC tumors (n = 15) and BrM (n = 19) to perform T cell and B cell repertoire (TCR/BCR) profiling. Relative to primary TNBC, TNBC BrM had lower read counts of T cell receptor alpha and beta (TRA, p < 0.001 and TRB, p < 0.01), with BCR heavy chain and light chain abundance showing trending but non-significant differences ( ). This result is in accordance with RNA-seq data that showed less T cell and B cell abundance in the primary samples relative to BrM ( ). Repertoire diversity was indexed as modeled Shannon entropy (37), which is a diversity index that accounts for both the richness of the sample (e.g., the number of unique TCR/BCR sequences) and relative species abundance (evenness) (68, 69). Thus, a large Shannon entropy score reflects a more diverse distribution of TCR/BCR sequences. The modeled Shannon entropy (TCR/BCR diversity) was lower for BrM compared to primary tumors (TRA, p < 0.01 and TRB, p < 0.05) ( ). A comparison of matched BrM and primary TNBC pairs only, however, did not show a reduction of TCR/BCR read counts and modeled Shannon entropy ( ).

Figure 4

T cell and B cell repertoire analysis revealed adaptive immune cell deficit in BrM relative to primary triple-negative breast cancer. The distribution of read counts and modeled Shannon entropy for all samples is displayed in (A, B), respectively. The same information is displayed, respectively, in (C, D) for matched pairs (note that the number of matched pairs varies due to the presence/absence of relevant reads). Wilcoxon rank-sum test was performed on log10-transformed (read counts) or raw (modeled Shannon entropy) values to determine the statistical significance. Significance codes: † p < 0.1; *p < 0.05; **p < 0.01; ***p < 0.001.

Differential Gene Expression and Pathway Analysis Support an Immune Cell Deficit in TNBC BrM Relative to Primary TNBC Tumors

Gene expression was evaluated by utilizing RNA-Seq data between the primary tumor (n = 15) and BrM (n = 19) tissues. In total, there were 1,669 differentially expressed genes (DEGs) between these 2 groups, with 935 genes upregulated and 734 genes downregulated in the BrM tissues compared to primary tumors (q ≤ 0.1; ; ). Gene ontology (GO) analysis of DEGs revealed a significant enrichment of immune-related terms in the primary TNBC tumors compared to the BrM (particularly terms reflecting adaptive immune system involvement), whereas GO terms associated with the nervous system were significantly higher in the BrM relative to the primary TNBC tumors ( ). Canonical pathway analysis (Ingenuity Pathway Analysis, IPA) of DEGs in primary tumors versus BrM illustrated a similar preponderance of immune signaling-related pathways as well as nervous system-related pathways associated with DEGs in BrM relative to primary tumors ( ). Upstream regulator analysis (IPA) further demonstrated an association of immune-related signaling activity with DEGs in primary tumors relative to BrM (e.g., IFNG, NFKB, CD3, CSF2, and IL-1β) and an association of potential oncogenic drivers (e.g., TCF7L2, mTOR, and SH3TC2) with regulation of BrM DEGs ( ; ).

Figure 5

Differential gene expression analysis and Gene Set Enrichment Analysis further support the immune cell deficit in BrM relative to primary triple-negative breast cancer (TNBC). (A) Volcano plot displaying differentially expressed genes (DEGs) in primary TNBC relative to BrM, with the legend showing color-coded levels of significance. (B) Top Gene Ontology terms associated with DEGs. DEGs with LFC > |1| and q < 0.1 (equating to 468 genes for primary_vs_BrM and 463 genes for BrM_vs_primary) were subjected to PANTHER overrepresentation test (dotted line represents q = 0.1). (C) Canonical pathway analysis (Ingenuity Pathway Analysis, IPA) of DEGs. The pathways displayed were significant at q <0.05 (dotted line) and were associated with a significant z-score (z >|2|) which indicates associative activity. (D) Upstream regulator analysis (IPA) displaying top regulators (z >|2|, q < 0.1) identified to be associated with an active or inhibited state in primary versus BrM TNBC [see for the full list of DEGs (q < 0.1) and upstream regulators (z >|2|, q <0.1).

Adaptive Immune Cell Signatures Are Associated With Improved Survival for Patients With TNBC BrM

We examined the survival association of standard clinicopathological variables (age at diagnosis of primary tumor, stage, and race) with different time metrics to event: (1) time from diagnosis of primary TNBC to death, (2) time from diagnosis of any metastatic disease to death, and (3) time from diagnosis of BrM to death. Of these variables, only older age was significantly associated with poor survival using each of these time metrics ( ), which was similar to other recent reports (70, 71). Next, survival associations relative to IGS expression were evaluated using multivariable CoxPH models in both primary TNBC and BrM. The IGS features in primary TNBC tumors which were associated with improved survival following metastatic diagnosis included T cell, B cell, and dendritic cell (DC) signatures ( ). Interestingly, a fibroblast serum response/wound healing signature (64) (“Chang_Serum_Response_Up”) was associated with a significantly poorer survival (p = 0.025) in BrM ( ) after a diagnosis of metastasis.

Discussion

In this study, we examined the genomic and transcriptomic landscape of TNBC BrM and primary tumors to further the understanding of TNBC BrM etiology and the tumor immune microenvironment. Despite recent progress in the treatment of ER+ and HER2+ BCBrM with newer brain-penetrant, targeted therapies, the treatment options for TNBC BrM remain largely restricted to chemotherapy and local therapy due to lack of known targets. A growing appreciation for the role of immunotherapy in the treatment of TNBC highlights the need to better understand the immune context of BrM as we consider incorporation of immunotherapy into the care of our patients (72).

Through whole-exome sequencing, we report that BrM, as a group, exhibited a greater TMB than primary tumors, though this observation was not recapitulated in matched tissue pairs. We suspect that this was due to underpowering of our study to assess matched pairs (n = 6) in the context of TNBC tumor biological heterogeneity (73). An analysis of shared variants showed that matched primary TNBC and BrM samples were more alike than inter-patient primaries and inter-patient BrM, showing that TNBC is a heterogeneous disease with potentially non-redundant mechanisms of tumorigenesis. A subclonal analysis of matched pairs also showed that some patients displayed patterns of divergent evolution between primary tumors and BrM. The mutational spectrum of genes with known causality to breast cancer tumorigenesis was also queried. TP53 was commonly mutated in this cohort (50%), supporting its causal role in the development of TNBC, while other genes such as MAP3K13 and PIK3CA were mutated at a lower frequency.

Copy number variation analysis revealed common and unique genomic alteration events between primary TNBC and BrM. 11p13 was commonly amplified in both primary tumors and BrM. This genomic location harbors the CD44 gene, which is used to discern breast CSCs, although it has been shown that it is not likely a driver of amplification of this region in basal breast cancer (74). WT1, which also resides at 11p13, has been shown to promote a mesenchymal phenotype in breast cancer cells as well as to elicit resistance to taxane therapy (47). Regions of 1q were also commonly amplified, which supports a known role for this genomic location in breast cancer development (48). The amplification and increased expression of ENSA (1q21.3) have recently been shown to drive TNBC progression via positive regulation of cholesterol biosynthesis (58). 13q11 was deleted in both primary and BrM TNBC, and this site is proximal to BRCA2 (located on 13q13.1). Whether or not the loss of 13q11 has any BRCA2-regulatory functionality is unknown, although deletions in 13q and 14q are common in BRCA2-mutated breast cancers (75). Common deletion of 20p13 was also observed. While this location is known to be deleted in colon cancer (76), its association with breast cancer has yet to be explored. Other common deletions identified included 1p36 and 11p15.5. In ductal breast carcinoma, the 1p36 deletion is associated with grade, ERBB2 loss, and loss of BCL2 expression (57) and is known to be a common feature underlying breast cancer development and the carcinogenesis of various cancer types (77). BrM-specific deletion at 11p15.5 (region harboring PRKCDBP) was also observed, and the chromosomal loss of this region is associated with BCBrM, with PRKCDBP identified as a putative tumor suppressor (53).

In primary TNBC, notable amplicons were associated with both arms of chromosomes 1 and 8. These locations are associated with breast cancer cytogenetics and pathology (48, 78, 79) and harbor genes (NOTCH2 and RCP, respectively) associated with breast cancer etiology (49, 80). There were also several notable alterations specific to BrM. Regions harboring the oncogenes BCL2L2 (14q11.2), AKT2 (19q13.2), and TFF3 (21q22.3) were amplified in BrM. BCL2L2 is an anti-apoptotic protein that has an oncogenic role in many solid tumor types, and it has been found to contribute to breast cancer progression through its upregulation via hypermethylation of the negative-regulatory miR-129-2 (81). While PIK3CA was only mutated in 10% of evaluated samples in this study, AKT2 upregulation via genomic amplification may have a significant impact on TNBC BrM progression. Dysregulation of the PI3K/AKT/mTOR axis is a common feature of TNBC (82), and this pathway represents a promising target in this disease context. TFF3 is also associated with breast cancer metastasis, where its expression predicts poor survival (83), and it is also associated with residual invasive disease following neoadjuvant chemotherapy in breast carcinoma (84). Interestingly, TFF3 was also found to be significantly upregulated in T cell-cold tumors of diverse tissue types, and it was in the top percentile of genes differentially expressed in T cell-cold versus T cell-hot breast cancers (85), which suggests its potential as an immunotherapy target. These mutational and copy number analyses together highlight potential causative genomic alterations contributing to TNBC progression and BrM.

A systematic evaluation of the neoantigen landscape in LCCC1419 was undertaken here. Using a suite of bioinformatics prediction software, we analyzed tumor-associated antigens (e.g., CTAs/self-antigens), traditional tumor-specific antigens (TSAs; e.g., SNVs), and alternative TSAs [e.g., derived from splice variants, chromosomal structural variants, InDels, ERVs, and other viral antigens (63)]. We found that both primary TNBC and BrM harbored substantial numbers of high-affinity MHC class I-binding peptides derived from CTAs and ERVs relative to other antigen sources. CTAs are known to be associated with aggressive hormone-negative breast cancers and poor survival; however, they have also been associated with robust immunogenicity in some contexts (86). ERVs, which are evolutionary remnants of viral insertional mutagenesis, are also potentially powerful immunogens (18). Although ERV transcriptional regulation is often epigenetically silenced in normal cells, tumor cell-specific derepression is known to occur and is associated with a response to immune checkpoint blockade (ICB) in multiple cancer types (18, 87, 88). As such, antigens derived from CTAs and ERVs may be invaluable immunotherapeutic targets for vaccine strategies targeting TNBC and BrM lesions. Relative to primary TNBC, we also found an elevated SNV mutational load associated with BrM. This augmented TSA burden in BrM also represents a potential vulnerability to be targeted by combination immunotherapeutic approaches, including neoantigen vaccine strategies.

A comprehensive analysis of transcriptomic data derived from this cohort was performed to further understand the difference between the tumor immune microenvironment of primary and BrM TNBC. We found that BrM lesions harbored significantly less immune infiltrate than primary tumors. This is not surprising, as the brain has historically been considered an immunologically protected organ (89). A recent study with RNA array data in BCBrM, agnostic to subtype of BC, has similarly reported reduced immune scores in BCBrM relative to primary tumors (90). The general dearth of immune involvement in the BrM spanned both adaptive (T and B cell) and innate (DC, eosinophils, and mast cells) immune populations, indicative of a broad immune deficit relative to primary tumors and again similar to recent reports (90). Interestingly, BrM displayed an elevated expression level of genes involved in a serum-induced fibroblast wound healing response (64). This finding may suggest that, relative to primary TNBC, BrM lesions are more reliant on aberrant wound healing properties, requiring increased levels of stromal involvement for growth and maintenance, as seminally put forth by Dvorak (91). We also observed that BrM had significantly decreased TCR (TRA/TRB) abundance and diversity as compared to primary tumors, and this association was verging on significance for certain immunoglobulin classes. These metrics are important, as increased TCR abundance and diversity have been associated with a response to ICB in multiple solid tumor types (92). DEGs between primary and BrM TNBC also reflected a BrM-specific immune deficit. Gene Ontology and canonical pathway analysis showed that genes that exhibited relatively lower expression levels in BrM versus primary tumors were enriched for terms related primarily to an adaptive immune response. An upstream regulator analysis further supported these findings, with IFNG being the putative regulator with the highest significance. In BrM, this upstream regulator analysis further demonstrated the importance of mTOR signaling but also showed that TCF7L2 and SH3TC2 may be important players in BrM development. TCF7L2 variants have been found to be associated with breast cancer incidence (93, 94). Additionally, this gene is a positive regulator of Wnt signaling, regulates the MYC oncogene, represses the cell cycle inhibitors CDKN2C/CDKN2D, and is a transcriptional driver of various oncogenes, contributing to the progression of colon cancer and other cancer types (95).

We performed survival analyses examining the prognostic potential of IGS in the context of standard clinicopathological features. Signatures representing levels of T cells, B cells, NK cells, and DC cells in the primary tumors were associated with improved survival. The collective association of these IGS with favorable survival is likely indicative of the high degree of expression correlation structure (and thus co-infiltration levels of the associated immune cell types) observed in this cohort ( ). Moreover, the favorable association between T cell, B cell, and DC IGS and survival may indicate that patients with higher anti-tumor immune infiltrate in their primary TNBC may have a higher propensity to develop long-lasting immunological memory that functions to stave off metastatic spread. Similarly, levels of signatures reflective of gamma-delta T cells and ICB responsiveness [e.g., Vincent_IPRES_Responder signature (96)] in BrM were associated with improved survival from the time of BrM diagnosis, indicating that elevated immune involvement in the brain TME may be beneficial to patient survival. Conversely, the aforementioned fibroblast wound healing signature [Chang_Serum_Response_Up (64)] was associated with poor survival in BrM, indicative of a deleterious quality of this signature and the underlying biology that it represents.

While this study represents the largest series focused on TNBC BrM to date, to our knowledge, it is mainly limited by low power, particularly regarding matched pairs (with only n = 6 matched WES and n = 9 matched RNA-seq pairs). An additional limitation is the inability to corroborate adaptive immune receptor repertoire inference with amplicon sequencing, which was precluded due to inadequate specimen nucleic acid abundance. Future work will expand these, and additional analyses to additional TNBC samples, as well as to other BC subtypes in the LCCC1419 biobank, including HER2+ and ER/PR+ BCBrM, to enable a comparison of BrM across the spectrum of BC. Utilization of in vivo murine models for testing the relevance of these findings, including the assessment of vaccine strategies and ICB as potential therapeutic approaches for TNBC BrM, is warranted.

In summary, we report the genomic characterization of BrM compared to primary tumors from TNBC patients, including some matched pairs, with a focus on the immune landscape. Utilizing both WES and RNA-seq analytical pipelines, we demonstrated that BrM exhibited increased TMB and SNV mutational load, reduced immune gene signature expression and TCR receptor abundance/diversity metrics, and increased expression of a wound healing signature. A prediction of elevated levels of CTA- and ERV-specific neoantigen peptides was confirmed in both anatomical locations, supporting the continued development of vaccine and immune checkpoint inhibition approaches in TNBC. IGS, including T cell-related immune signatures in primary and BrM TNBC, correlated with improved survival in this patient cohort. We expect that these results and the data reported herein will be valuable in understanding TNBC BrM biology going forward and provide further rationale for the application of immunotherapeutic approaches in this disease.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The name of the repository and accession number can be found below:

https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs002457.v1.p1.

Ethics Statement

The studies involving human participants were reviewed and approved by the University of North Carolina at Chapel Hill Institutional Review Board. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

Conception: CA, BV, MM, and LCa. Design: ER, AVS, MS, MM, BV, SC, CA, and LiC. Acquisition/analysis: MS, LuC, LiC, AG, AW, ER, SV, SC, MW, JP, MM, AVS, BV, CA, and AH. Interpretation: ER, AVS, MS, SC, JP, MM, MW, CA, and BV. Manuscript preparation/editing: all authors. All authors contributed to the article and approved the submitted version.

Funding

This study received funding from Susan G. Komen Career Catalyst Award (BV), V Foundation for Cancer Research Translational Grant (BV and CA), AACR grant (CA), and Translating Duke Health (CA).

Conflict of Interest

JP is an inventor on patent applications for the Breast PAM50 assay. BV and JP report equity and consulting fees for GeneCentric Therapeutics. CA receives research funding from PUMA, Lilly, Merck, Seattle Genetics, Nektar, Tesaro, G1-Therapeutics, ZION, Novartis, and Pfizer; compensated consultant role: Genentech (1/2019-), Eisai (1/2019-), IPSEN (2/2019-), Seattle Genetics (11/15/2019-11/15/2020); Astra Zeneca (3/2020-6/2020), Novartis (5/2020-5/2022), Immunomedics (10/1/2020-9/22/2021), Elucida (9/2020); and royalties from UpToDate and Jones and Bartlett. LC has received institutional research funding from Syndax, Novartis, NanoString Technologies, AbbVie, Seattle Genetics, and Veracyte. An immediate family member has a royalty-sharing agreement and investorship interest in licensed IP to startup company Falcon Therapeutics that is designing neural stem-cell-based therapy for glioblastoma multiforme. She has uncompensated relationships with Sanofi, Novartis, G1 Therapeutics, Genentech/Roche, GlaxoSmithKline, AstraZeneca/Daiichi Sankyo, Aptitude Health, Exact Sciences and Eisai.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.