Abstract
Purpose: Lymphocytic infiltration of tumors predicts improved survival in patients with breast cancer. Previous studies have suggested that this survival benefit is confined predominantly to the basal-like subtype. Immune infiltration in ovarian tumors is also associated with improved prognosis. Currently, it is unclear what aspects of the immune response mediate this improved outcome.
Experimental Design: Using The Cancer Genome Atlas mRNA-seq data and a large microarray dataset, we evaluated adaptive immune gene expression by genomic subtype in breast and ovarian cancer. To investigate B-cells observed to be prognostic within specific subtypes, we developed methods to analyze B-cell population diversity and degree of somatic hypermutation (SHM) from B-cell receptor (BCR) sequences in mRNA-seq data.
Results: Improved metastasis-free/progression-free survival was correlated with B-cell gene expression signatures, which were restricted mainly to the basal-like and HER2-enriched breast cancer subtypes and the immunoreactive ovarian cancer subtype. Consistent with a restricted epitope-driven response, a subset of basal-like and HER2-enriched breast tumors and immunoreactive ovarian tumors showed high expression of a low-diversity population of BCR gene segments. More BCR segments showed improved prognosis with increased expression in basal-like breast tumors and immunoreactive ovarian tumors compared with other subtypes. Basal-like and HER2-enriched tumors exhibited more BCR sequence variants in regions consistent with SHM.
Conclusion: Taken together, these data suggest the presence of a productive and potentially restricted antitumor B-cell response in basal-like breast and immunoreactive ovarian cancers. Immunomodulatory therapies that support B-cell responses may be a promising therapeutic approach to targeting these B-cell infiltrated tumors. Clin Cancer Res; 20(14); 3818–29. ©2014 AACR.
This article is featured in Highlights of This Issue, p. 3627
Immunomodulatory therapies such as those targeting CTLA4 and PDL1 (CD274) have been shown to be effective in a number of tumor types. These treatments primarily target the T-cell component of the adaptive immune system. In this work, we provide further evidence that tumor-infiltrating B-cells are also important in antitumor immunity in basal-like breast and immunoreactive ovarian cancer, suggesting that immunomodulatory therapy may be effective in these tumor types. By identifying basal-like breast cancer as an important setting for immunomodulatory treatment, we provide a rationale for the use of targeted therapies in many clinically triple-negative breast cancers, in which no targeted therapies currently exist.
Introduction
The role of tumor-infiltrating lymphocytes (TIL) in breast cancer is not fully understood, although multiple studies have shown an association between the presence of TILs and an improved prognosis (1–5). TILs in breast tumors are predominantly cytotoxic (CD8+) T-cells (6, 7), and the proportion of CD8+ T-cells may be prognostic (4, 5, 8). In contrast, TILs of the regulatory T-cell phenotype (CD4+CD25+FOXP3+ Tregs) are associated with poorer outcomes in breast cancer (9, 10). The role of B-cell TILs in human breast cancer is not as clear as that of T-cell TILs. Using gene expression profiling, our group and others have showed that gene signatures representing B-cells, plasmablasts, plasma cells, and immunoglobulin predicted favorable clinical outcome in estrogen receptor-positive (ER+) and estrogen receptor-negative (ER−) breast tumors (11–15). In this article, these are referred to as B-cell signatures; although plasmablasts and plasma cells are known to infiltrate some breast tumors, we use the term “B-cell TIL” here to refer to any TIL in the B-cell lineage. The presence of B-cell TILs as assessed by immunohistochemistry (IHC) has also been shown to be an independent prognostic feature in breast cancer (16). Studies of small numbers of breast tumors have shown the B-cell response in these tumors to be clonally expanded, with evidence of having undergone class switching and somatic hypermutation (SHM; refs. 17–22). This strongly suggested that in some breast tumors there may be a clonally restricted, antigen-directed B-cell antitumor response. Several studies have identified auto-antibodies in patients with breast cancer, including antibodies against improperly processed β-actin in some medullary breast cancers, although the association between such auto-antibodies and patient survival is unclear (18, 21, 23). Together, these findings provide evidence that B-cell TILs may be important in affecting breast cancer biology and progression.
Human breast cancer is a heterogeneous disease, with individual tumors varying according to morphology, natural history, and response to therapy. Gene expression analyses have identified at least five distinct genomic subtypes of breast cancer: luminal A, luminal B, HER2-enriched, basal-like, and claudin-low, as well as a normal-like group (24–28). The prognostic value of both T- and B-cell TILs may be restricted to a subset of highly immune-infiltrated breast tumors (14). Basal-like breast tumors, in particular, seem to have beneficial TILs (5, 15). Multiple groups have identified signatures of lymphocyte-related gene expression that are overrepresented in basal-like breast tumors and predict better survival (14, 15); in contrast, luminal A breast tumors show low levels of lymphocytic infiltrate (5).
Comprehensive genomic profiling of multiple tumor types in The Cancer Genome Atlas (TCGA) has shown that there is a strong similarity between basal-like breast cancer and serous ovarian cancer (24). These two tumor types exhibit a similar mutational spectrum and share many of the same driver events (i.e., TP53 loss, RB1 loss, c-MYC gain, etc.). Like basal-like breast cancer, many ovarian tumors are rich in TILs. Analysis of TCGA serous ovarian cancer gene expression identified four genomic subtypes: mesenchymal, proliferative, differentiated, and immunoreactive (29). The immunoreactive subtype, in particular, showed high expression of T-cell chemokine ligands and lymphocyte-related genes. Furthermore, a number of studies have shown that the presence of T- and B-cell TILs is a positive prognostic feature in ovarian cancer (30–33). As in breast cancer, the precise role of B-cell TILs is less understood than that of T-cell TILs. These data suggest that, like basal-like breast cancer, serous ovarian cancer may be a likely candidate for identifying a productive antitumor T-cell and/or B-cell TIL response.
If there is an effective, subtype-specific antitumor response mediated by B-cell TILs, this presents the possibility of subtype-specific immunogenic epitopes that could promote development of a subtype-specific antibody response. Although some studies have identified antigen-directed TIL clones in breast tumors (18–22), currently the degree to which TILs are antigen-directed is unknown. The development of a mature B-cell response following antigen stimulation depends on a number of processes that occur during the germinal center reaction, including clonal expansion and antibody class switching (34, 35). While in the germinal center, B cells undergo SHM, whereby mutations at BCR loci are introduced to enhance B-cell receptor (BCR) affinity. Mutations occur preferentially in “hot spot” nucleotide positions, particularly within the antigen-binding complementarity-determining regions (CDR), and they favor replacements and transitions (36). In a tumor antigen–driven response, the TIL population is expected to be enriched for one or more “dominant clones” exhibiting BCR characteristics consistent with SHM. Here, using a novel approach to characterize B-cell responses from short read mRNA-seq data, we demonstrate subtype-specific enrichment of B-cell gene segments in basal-like and HER2-enriched breast cancer and in immunoreactive ovarian cancer. We show evidence of clonal restriction of the B-cell response in these three tumor types, and mutation patterns consistent with SHM in basal-like and HER2-enriched breast cancer. These findings suggest an important role for the endogenous B-cell response specifically in these tumor subtypes.
Materials and Methods
Datasets
The breast cancer dataset used for all analyses except the survival analysis of gene expression signatures was the TCGA dataset of 819 mRNA-seq samples, comprising 728 breast tumors and 91 normal breast samples (see TCGA data portal at https://tcga-data.nci.nih.gov/tcga/, CGHub at https://cghub.ucsc.edu/). This 728 sample set is an extension of the 480 tumors previously profiled by microarray (24), but these 480 plus 350 new samples have all been assayed by mRNA-seq using Illumina 2 × 50-bp sequencing as described by TCGA in an evaluation of lung squamous carcinoma samples (37). Gene expression values were represented as RSEM (RNA-seq by expectancy maximization) data normalized within sample to the upper quartile of total reads as previously described (37). These data and further details about data processing are available at the TCGA data portal under the V2_MapSpliceRSEM workflow (https://tcga-data.nci.nih.gov/tcgafiles/ftp_auth/distro_ftpusers/anonymous/tumor/brca/cgcc/unc.edu/illuminahiseq_rnaseqv2/rnaseqv2/unc.edu_BRCA.IlluminaHiSeq_RNASeqV2.mage-tab.1.6.0/DESCRIPTION.txt). Genomic subtype was assigned within the set of 728 mRNA-seq samples using the PAM50 assay (38). The training set of breast samples used in the PAM50 assay is 50% clinically ER+; therefore, the mRNA-seq data were normalized to reflect the training set (https://genome.unc.edu/pubsup/breastGEO/Guide%20to%20Intrinsic%20Subtyping%209-6-10.pdf). On the basis of the clinical data taken from the TCGA data portal on September, 2012 (https://tcga-data.nci.nih.gov/tcga/), of the 728 samples, 157 were ER−, 535 were ER+, two were ER indeterminate, 29 did not have ER status assays performed, and five did not have available data, indicating that 77% of the mRNA-seq samples were ER+. To normalize the data similar to the PAM50 training dataset, in which 50% of samples are ER+, all 157 ER− samples were selected, as well as 157 randomly selected ER+ samples. The median gene expression for the PAM50 intrinsic gene list was calculated on the basis of this subset of samples. To perform platform correction for mRNA-seq, these median values were then subtracted from all 728 samples before running the PAM50 assay as previously described (https://genome.unc.edu/pubsup/breastGEO/PAM50.zip; ref. 38). Because of the short median follow-up time (17 months) of the TCGA dataset, survival analyses of gene expression signatures were carried out using a microarray-based gene-expression dataset of 855 breast tumors with published intrinsic subtype calls (140 basal-like, 90 claudin-low, 144 HER2-enriched, 243 luminal A, 162 luminal B, and 76 normal-like) and clinical data (combined data from the following datasets: GSE2034, GSE12276, GSE2603, and the NKI295; ref. 39). Survival analyses of BCR segment expression, however, used the TCGA mRNA-seq dataset. For all analyses of ovarian cancer, we used the TCGA serous ovarian cancer mRNA-seq dataset, which, like the breast cancer dataset, represents new mRNA-seq data, again using Illumina 2 × 50-bp sequencing, on a subset of the 500 cases from the TCGA ovarian project (29). This mRNA-seq dataset consists of 266 tumors with follow-up data (https://tcga-data.nci.nih.gov/tcga/, https://cghub.ucsc.edu/). The TCGA barcodes and genomic subtypes of all breast and ovarian mRNA-seq samples are included as Supplementary Data S1. Unique sample IDs for downloading TCGA mRNA-seq data from CGHub are included as Supplementary Data S3.
Gene expression signatures and survival analyses
Immune gene expression signatures were established using unsupervised hierarchical clustering of mRNA-seq expression data for 728 breast tumor samples (see Supplementary Data). Gene dendrogram nodes corresponding to genes characteristically expressed in specific immune cell types were identified and validated through DAVID functional annotation clustering and Ingenuity Systems Analysis (www.ingenuity.com; refs. 40, 41). Gene lists for all five signatures are included in Supplementary Data S2. Additional lymphocyte gene signatures were obtained from published studies: IGG_Cluster (11), B_Cell (42), and B_Cell_60 gene (12) are B-cell signatures, and T_Cell (42), CD8 (42), LCK (43), and TNBC_T-Cell (15) are T-cell signatures, with the CD8 signature specifically representing CD8+ T-cells.
Survival analyses were performed by Kaplan–Meier analysis and log-rank testing, and HRs were derived from the Cox proportional hazards model. For analysis of the prognostic value of BCR segment expression, samples were divided into high- and low-expression groups of equal number for Kaplan–Meier analysis and log-rank testing. To evaluate the prognostic value of gene expression signatures, the Cox proportional hazards model was used with each signature tested as a continuous variable. Multivariate Cox proportional hazards models were used to test the prognostic value of individual gene expression signatures when combined with the other clinical and genomic variables.
To provide a control for the number of prognostic BCR gene segments, 353 (the number of BCR gene segments tested) random genes were selected and P values were calculated for the association of each of these genes with overall survival (breast) or progression-free survival (PFS; ovarian) as in the BCR gene segments. The number of significant (P < 0.05) P values was calculated from this set and 95% confidence intervals were calculated through bootstrap resampling.
BCR diversity
The method for estimating sequence diversity of a BCR gene segment for an individual sample/tumor using paired-end mRNA-seq data is outlined in Supplementary Fig. S3. Read pairs mapping to the EntrezGene genomic location (www.ncbi.nlm.nih.gov/entrez/query.fcgi?db=gene; ref. 44) of a given BCR gene segment were identified (mapped by MapSplice; ref. 45). The sequence of these read pairs was compared with the hg19 reference genome to identify nonreference bases. The genomic position and nucleotide identity of all nonreference bases were identified for each read pair. Each observed pattern of nonreference bases was then assigned a score representing the number of read pairs containing exactly that pattern of nonreference bases. This set of observed patterns and their corresponding count was used to calculate the effective number of species, which is a diversity function isomorphic to Shannon entropy, as described by Jost (46).
De novo assembly and SHM analysis
De novo assembly of BCR variable (V) gene segments from paired-end mRNA-seq reads was performed using the Assembly-based Re-Aligner (ABRA) algorithm (Mose and colleagues; unpublished data). To generate ABRA contigs, unmapped reads and reads mapping to a BCR variable region of interest were first split into overlapping K-mers, where K = 31. K-mers that were composed exclusively of nonambiguous bases with the quality score >20 were assembled into a de Bruijn graph. K-mers with less than 100 observations were then pruned from the graph. The graph was then traversed to produce all possible contigs. This set of contigs was used for SHM analyses.
To determine whether the sequence of a BCR variable gene segment was consistent with SHM, the reference sequence for that gene segment was first established by Smith–Waterman alignment to each IMGT (the international ImMunoGeneTics information system http://www.imgt.org founder and director: Marie-Paule Lefranc, Montpellier, France; refs. 47–53) reference allele sharing the same gene segment family and selection of the closest match. The segment sequence was compared with the reference sequence, and mutated and nonmutated bases were counted within SHM hotspots (WRCY, RGYW, WA, and TW sequences) and non-hotspot regions (54). Mutated and nonmutated (i.e., reference and nonreference) bases were counted again within CDRs. Then, χ2 testing was used to determine whether the distribution of mutated bases was consistent with the mutation pattern expected in SHM; χ2 testing was conducted separately for the whole V segment sequence and for CDRs.
The method for producing an overall estimate of the degree of SHM for a BCR variable segment in mRNA-seq data is outlined in Supplementary Fig. S4. ABRA contigs were first constructed for that V gene segment and quantitated by Burrows-Wheeler alignment (BWA) paired-end alignment (55) of all unmapped reads and reads mapping to that segment to the set of ABRA contigs. Multiple mappings were allowed, as long as both read pairs mapped to the same contig. For each contig, the number of SHM hotspot mutations, SHM hotspot nonmutated bases, non-hotspot mutations, and non-hotspot nonmutated bases were counted for the whole segment and within CDRs. This information was weighted by the BWA alignment score for that contig. These weighted values were summed across all contigs and χ2 testing was used to determine whether the mutation pattern across the whole segment or within CDRs indicated the presence of SHM. For this analysis, only segments with a number of mapped reads greater than or equal to 0.04 × the length in bp of the segment (approximating even coverage of depth 2) were considered for sequence analysis in a given sample.
Results
B-cell gene expression signatures are prognostic in breast and ovarian cancer
Increased expression of B-cell gene signatures has been shown to be favorably prognostic in breast cancer (11–13). To explore the role of B cells and other lymphocyte cell types in the different intrinsic subtypes of breast cancer, immune cell–associated genomic signatures were newly derived from unsupervised clustering of mRNA-seq data from 728 TCGA breast cancer samples (Fig. 1A). Gene dendrogram nodes containing characteristic lymphocyte genes were selected as potential gene signatures. The identities of these signatures were confirmed through functional annotation analysis and gene pathway–based analysis (40, 41). These and other previously published lymphocyte gene signatures were used to confirm the prognostic value of TILs on a genomic level, and to assess whether this benefit is isolated to one or more intrinsic subtypes. We first evaluated a gene expression microarray dataset of 855 breast tumors, using a univariate Cox proportional hazards model, for prognostic value by subtype of lymphocyte gene signatures (39). For the newly derived B_Cell_cluster signature derived from unsupervised clustering here, the IGG_Cluster previously developed by our group (11), and three B-cell gene signatures generated by others (12, 14, 15), overall expression in breast tumors was greater in the basal-like and HER2-enriched subtypes (Fig. 2A and C). Similar to previous work (12, 16, 56), high expression was associated with better metastasis-free survival in basal-like and HER2-enriched tumors with greatest difference in hazards in the basal-like subtype (Table 1).
. | . | IGG_Clustera . | TNBC_B-Cellb . | B_Cella . | B_Cell_60 genec . | B_Cell_clusterd . | T_Cella . | TNBC_T-cellb . | LCKa . | CD8a . | T_Cell_clusterd . | CD8_clusterd . | |||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
. | n . | HR . | P . | HR . | P . | HR . | P . | HR . | P . | HR . | P . | HR . | P . | HR . | P . | HR . | P . | HR . | P . | HR . | P . | HR . | P . |
Breast | |||||||||||||||||||||||
All breast | 855 | 0.847 | 6.61E−04 | 0.901 | 1.56E−03 | 0.585 | 1.52E−02 | 0.888 | 4.79E−04 | 0.782 | 9.38E−04 | 0.562 | 5.31E−03 | 0.901 | 3.83E−02 | 0.888 | 3.37E−02 | NS | 0.858 | 2.27E−02 | 0.865 | 1.81E−02 | |
Luminal A | 243 | NS | NS | NS | NS | NS | NS | NS | NS | NS | NS | NS | |||||||||||
Luminal B | 162 | NS | NS | NS | NS | NS | 0.304 | 1.59E−02 | NS | NS | NS | NS | NS | ||||||||||
HER2-enriched | 144 | 0.755 | 4.67E−03 | 0.827 | 1.12E−02 | 0.323 | 9.16E−03 | 0.791 | 2.43E−03 | 0.641 | 4.34E−03 | 0.328 | 6.01E−03 | 0.745 | 4.74E−03 | 0.725 | 6.39E−03 | 0.658 | 3.29E−02 | 0.687 | 5.03E−03 | 0.705 | 4.15E−03 |
Basal-like | 140 | 0.599 | 1.20E−04 | 0.686 | 4.24E−05 | 0.17 | 2.39E−03 | 0.682 | 6.05E−05 | 0.491 | 3.56E−04 | 0.126 | 2.44E−04 | 0.644 | 1.00E−03 | 0.619 | 1.34E−03 | 0.472 | 1.53E−03 | 0.496 | 1.29E−04 | 0.548 | 2.35E−04 |
Claudin-low | 90 | NS | NS | NS | NS | NS | NS | NS | NS | NS | NS | NS | |||||||||||
Ovarian | |||||||||||||||||||||||
All ovarian | 266 | 0.845 | 7.40E−03 | NS | NS | 0.857 | 4.09E−03 | 0.856 | 5.63E−03 | NS | NS | NS | NS | NS | NS | ||||||||
Immunoreactive | 77 | 0.675 | 7.43E−03 | NS | NS | 0.729 | 1.68E−02 | 0.725 | 1.22E−02 | 0.177 | 2.25E−02 | NS | NS | NS | 0.686 | 2.99E−02 | 0.707 | 4.28E−02 | |||||
Mesenchymal | 42 | 0.673 | 3.32E−02 | NS | NS | 0.714 | 4.31E−02 | NS | NS | NS | NS | NS | NS | NS | |||||||||
Differentiated | 65 | NS | NS | NS | NS | NS | NS | NS | NS | NS | NS | NS | |||||||||||
Proliferative | 78 | NS | NS | NS | NS | NS | NS | NS | NS | NS | NS | NS |
. | . | IGG_Clustera . | TNBC_B-Cellb . | B_Cella . | B_Cell_60 genec . | B_Cell_clusterd . | T_Cella . | TNBC_T-cellb . | LCKa . | CD8a . | T_Cell_clusterd . | CD8_clusterd . | |||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
. | n . | HR . | P . | HR . | P . | HR . | P . | HR . | P . | HR . | P . | HR . | P . | HR . | P . | HR . | P . | HR . | P . | HR . | P . | HR . | P . |
Breast | |||||||||||||||||||||||
All breast | 855 | 0.847 | 6.61E−04 | 0.901 | 1.56E−03 | 0.585 | 1.52E−02 | 0.888 | 4.79E−04 | 0.782 | 9.38E−04 | 0.562 | 5.31E−03 | 0.901 | 3.83E−02 | 0.888 | 3.37E−02 | NS | 0.858 | 2.27E−02 | 0.865 | 1.81E−02 | |
Luminal A | 243 | NS | NS | NS | NS | NS | NS | NS | NS | NS | NS | NS | |||||||||||
Luminal B | 162 | NS | NS | NS | NS | NS | 0.304 | 1.59E−02 | NS | NS | NS | NS | NS | ||||||||||
HER2-enriched | 144 | 0.755 | 4.67E−03 | 0.827 | 1.12E−02 | 0.323 | 9.16E−03 | 0.791 | 2.43E−03 | 0.641 | 4.34E−03 | 0.328 | 6.01E−03 | 0.745 | 4.74E−03 | 0.725 | 6.39E−03 | 0.658 | 3.29E−02 | 0.687 | 5.03E−03 | 0.705 | 4.15E−03 |
Basal-like | 140 | 0.599 | 1.20E−04 | 0.686 | 4.24E−05 | 0.17 | 2.39E−03 | 0.682 | 6.05E−05 | 0.491 | 3.56E−04 | 0.126 | 2.44E−04 | 0.644 | 1.00E−03 | 0.619 | 1.34E−03 | 0.472 | 1.53E−03 | 0.496 | 1.29E−04 | 0.548 | 2.35E−04 |
Claudin-low | 90 | NS | NS | NS | NS | NS | NS | NS | NS | NS | NS | NS | |||||||||||
Ovarian | |||||||||||||||||||||||
All ovarian | 266 | 0.845 | 7.40E−03 | NS | NS | 0.857 | 4.09E−03 | 0.856 | 5.63E−03 | NS | NS | NS | NS | NS | NS | ||||||||
Immunoreactive | 77 | 0.675 | 7.43E−03 | NS | NS | 0.729 | 1.68E−02 | 0.725 | 1.22E−02 | 0.177 | 2.25E−02 | NS | NS | NS | 0.686 | 2.99E−02 | 0.707 | 4.28E−02 | |||||
Mesenchymal | 42 | 0.673 | 3.32E−02 | NS | NS | 0.714 | 4.31E−02 | NS | NS | NS | NS | NS | NS | NS | |||||||||
Differentiated | 65 | NS | NS | NS | NS | NS | NS | NS | NS | NS | NS | NS | |||||||||||
Proliferative | 78 | NS | NS | NS | NS | NS | NS | NS | NS | NS | NS | NS |
We next performed a similar analysis on TCGA ovarian cancer data. Overall B-cell gene signature expression was increased in immunoreactive ovarian tumors (Fig. 2E and F). Several B-cell gene signatures were prognostic for PFS in the immunoreactive ovarian tumor subtype, which was not true for the other subtypes (Table 1). T-cell signatures (14, 15, 43) were also evaluated, and showed a similar pattern of expression and prognostic value (Fig. 2B and D). In multivariate survival analysis of individual immune signatures with other clinical and genomic features in breast cancer, most B-cell and T-cell signatures remained significantly prognostic (Supplementary Table S1).
For both breast and ovarian tumors, B-cell gene signature expression strongly correlated with other lymphocyte gene signatures, including those representing T-cells and macrophages (Supplementary Fig. S1). Likelihood ratio testing was performed to assess the independence of immune gene signatures as predictors of survival. Conditioning on clinical variables (node status, age, and hormone receptor status for breast; stage, grade, and age for ovarian) and adding either B-cell or T-cell gene signatures to the model, only one gene signature was needed to significantly improve the predictive power of the model (data not shown). In accordance with the high degree of correlation between immune cell signatures, adding further signatures for the same cell type did not improve the model. Each ordering of immune signatures was tested to ensure this finding was not specific to particular signatures. Including both B-cell and T-cell signatures in the model, in breast cancer only one signature significantly improved the model; in ovarian, one B-cell signature and one T-cell signature each significantly improved the model. Together, these analyses indicate an improved outcome for patients with specific subtypes of breast and ovarian cancer. This correlated with the presence of B-cells, plasmablasts, and/or plasma cells in the tumor microenvironment, which suggests that a productive endogenous B/plasma cell response may be present in the tumor microenvironment.
Specific BCR gene segment expression is prognostic in basal-like breast cancer
Next, we wished to determine whether the B-cell gene signature found in patients with basal-like breast cancer was consistent with an antigen-specific response. Other groups have shown clonal expansion and SHM in breast B-cell TILs, suggesting an antigen-directed response in those samples (18, 20, 22). Actively responding antigen-specific B-cell populations are characterized by clonal expansion; thus, we expect B-cell clonal expansion in patients where an effective, antigen-directed antitumor response is occurring. Because the clonal diversity of a B-cell population can be inferred by the diversity of the BCRs they express, there should be a prognostic benefit in samples with increased expression of specific BCR gene segments (i.e., immunoglobulin heavy chain and light chain variable, joining, diversity, and constant region segments). It has been shown that the BCR protein from breast cancer TILs is mainly produced by plasma cells, not B cells (13). Here, we will continue to use the term “B-cell” to refer to the heterogeneous group of BCR-producing cells in the B-cell lineage.
We first calculated the expression levels of all 353 BCR gene segments available in the IMGT database across the breast and ovarian tumor data. Breast HER2-enriched and basal-like tumors, as well as ovarian immunoreactive and mesenchymal tumors, showed high expression widely across BCR segments (Fig. 3). We analyzed prognosis by expression of each individual BCR segment, and then compared this with an identical number of randomly selected genes using a bootstrap procedure to assess the significance of this finding. In basal-like breast tumor subtypes, we identified a significantly greater number of BCR segments predictive of overall survival than expected by random sampling (Fig. 4A and B). No other breast cancer subtype demonstrated a greater number of prognostic BCR segments than expected by chance (Fig. 4A). Similarly, in patients with ovarian cancer significantly more BCR segments were predictive of PFS in immunoreactive tumors than in the other subtypes (Fig. 4C and D). The mesenchymal and differentiated ovarian subtypes also showed significantly more prognostic BCR segments than expected by chance. In breast tumors, this finding cannot be explained solely by increased overall expression of BCR gene segments in the basal-like subtype, as the highest expression of BCR segments was found in the HER2-enriched subtype (Fig. 3). Prognostic segments were discovered in multiple gene families, including IgHV, IgHJ, IgHC, IgKV, IgKJ, IgKC, IgLV, IgLJ, and IgLC. IgKC, which has been previously identified as prognostic in several solid tumor types, including breast cancer (13), predicted PFS in ovarian cancer but did not attain significance in the breast dataset. Individual representative plots of overall survival or PFS by high versus low expression of representative prognostic segments are shown in Supplementary Fig. S2.
Because B-cells undergo SHM following antigen stimulation in the germinal center reaction, reads mapping to each germline BCR gene segment are expected to contain many corresponding single-nucleotide variations. Each group of mapped reads corresponding to a BCR gene segment would then exist as a population, the diversity of which should inversely relate to the degree of clonal expansion in the tumor infiltrate. We used the Shannon entropy index normalized as the effective number of species as a measure of diversity (46). For this analysis, we calculated the diversity per patient of each BCR gene segment; a description of this procedure is shown as Supplementary Fig. S3. Basal-like, HER2-enriched, and luminal B breast tumors, and immunoreactive ovarian tumors, include a subset of tumors with high expression of low-diversity segments (Fig. 5). This finding is consistent with the presence of a clonally expanded B-cell population within those tumor subtypes that is absent in other subtypes.
Analysis of SHM patterns in mRNA-seq data
SHM in BCR gene segments is characterized by mutations that favor defined local sequence regional “hotspots” and CDRs, due to bias in the enzymatic activity that facilitates the mutation process (36). To evaluate the degree of SHM represented in our data, we made use of the novel de novo assembly algorithm ABRA to assemble unique contigs from reads that map to each BCR variable (V) segment locus, followed by analysis of the contigs for the presence or absence of SHM. These contigs allowed us to analyze SHM mutation patterns across a V segment or its CDRs, rather than interrogating each mRNA-seq read pair separately. An overview of this method is given as Supplementary Fig. S4.
We applied our method of analyzing SHM in mRNA-seq data to the TCGA breast and ovarian datasets. For the top 10 most highly expressed BCR V gene segments in breast or ovarian tumors in our datasets, the basal-like and HER2-enriched breast subtypes were enriched for tumors with V gene segments consistent with SHM (Supplementary Fig. S5). Immunoreactive ovarian tumors showed a high proportion of segments with mutation patterns suggestive of SHM, but it was not significantly higher than the proportion observed in other ovarian subtypes. The presence of SHM sequence characteristics from TILs is suggestive of the presence of antigen-experienced B-cells, potentially against tumor antigens, in the tumor microenvironment.
Discussion
We define here four characteristics of an active, antigen-driven, antitumor B-cell response that can be identified from mRNA-seq data, namely: (i) increased expression and prognostic value of B-cell gene signatures, (ii) increased expression and prognostic value of BCR gene segments, (iii) decreased diversity of highly expressed BCR gene segments, and (iv) mutation patterns consistent with BCR SHM. All four conditions were found in basal-like breast cancers, and three of these conditions were found for immunoreactive ovarian tumors and HER2-enriched breast tumors. These findings support the hypothesis that a productive B-cell–driven endogenous antitumor response may be generated in many basal-like breast and immunoreactive ovarian carcinomas. To our knowledge this represents the first inference of BCR repertoire characteristics from mRNA-seq data.
Investigations into the anticancer adaptive immune response have largely been focused on T-cells. Accordingly, current cancer immunotherapy is directed at modifying the T-cell immune response through modulating targets like CTLA4 and PDL1. In this work, we show that the presence of tumor-infiltrating B-cells correlated with overall and PFS suggesting that B-cells play an important role in antitumor immunity. We do show that the expression of B-cell genes was highly correlated with the expression of T-cell genes. By further demonstrating that in specific breast and ovarian cancer subtypes B-cell TILs are clonally expanded and enriched for SHM, we provide evidence that B-cell TILs are not merely a surrogate marker for an antitumor T-cell response. Although it is technically possible that previously expanded B-cell clones may be trafficked to the tumor independent of their antigen-binding capability, previous studies showing clonal evolution within breast tumors make this unlikely (19–22), as does the association with specific tumor subtypes. Tumor antigen-directed B-cell responses, which we suggest are present in many basal-like breast and immunoreactive ovarian tumors, may provide a novel way to clinically target these tumor types. Interestingly PD1 (PDCD1), which is expressed on activated T- and B-cells, is currently a very promising target for immunotherapy. Previous work has shown that blocking the interaction of PD1 with PDL1 and PDL2 (PDCD1LG2) enhances the activation, proliferation, and cytokine production of human B-cells in the presence of Toll-like receptor stimulation (57). As immunotherapy advances in breast cancer, it will be important to evaluate B-cell TILs to investigate whether and how anticancer immunotherapies may modulate the B-cell compartment.
As more immunomodulatory treatments become available for cancer therapy, one critical issue is the identification of the specific patients who may benefit from such therapy. This work highlights the subtype association of clonally restricted B-cell responses. Previous studies in ovarian cancer have been mixed as to the importance of B-cell TILs, perhaps because of the heterogeneity of the B-cell response across the subtypes of ovarian cancer. Milne and colleagues highlighted the high-grade serous histologic subtype as being selectively associated with TILs predictive of disease-specific survival (32). Here, we further identify the immunoreactive genomic subtype of serous ovarian cancer as containing prognostic TILs. Among basal-like breast tumors, it is interesting to note that patients with high B-cell infiltration as assessed by gene signatures were also significantly younger than other patients with basal-like breast cancer (data not shown), corroborating previous work highlighting this group (58).
Other investigations of B-cell TILs in breast cancer have found the survival benefit associated with B-cell gene expression to be dependent on proliferation (12, 13, 56). However, we do not see this association in our data. Although the basal-like and HER2-enriched subtypes are both highly proliferative, we observed no survival benefit for B-cell TILs in luminal B tumors, which are also characterized by high proliferation. Furthermore, likelihood ratio testing conditioning on clinical variables and genomic subtype (data not shown) demonstrated that proliferation did not significantly increase the predictive ability of the model in breast cancer.
This work again underscores the similarity between basal-like breast and ovarian tumors. Previous genomic studies have established that serous ovarian tumors resemble basal-like breast cancer in terms of their mutational profiles and DNA copy-number changes (24). Here, we show that this similarity extends to the immune component of the tumor microenvironment. In terms of the immune response, basal-like breast cancer bears more similarity to ovarian immunoreactive cancer than to the luminal A breast cancer subtype. This adds further weight to the notion that the therapeutic approach to basal-like breast and ovarian cancer could be similar.
The claudin-low subtype of breast cancer is also known to have abundant TILs (26), which we confirm by expression of TIL gene signatures (Fig. 2A and B). However, unlike the basal-like and HER2-enriched subtypes, TIL abundance was not associated with a survival benefit within claudin-low tumors (Table 1). This could potentially be due to different immunosuppressive mechanisms within the tumor microenvironment, or it is possible that claudin-low tumors elicit a nonspecific inflammatory response in contrast with other high-immune breast tumors. We were unable to assess the BCR sequence diversity of claudin-low TILs as very few (fewer than 10) claudin-low tumors have been identified within the TCGA dataset. If TILs within claudin-low tumors are not productive or antigen-directed, misclassification of these tumors may limit the effectiveness of immunomodulatory treatments within triple-negative breast cancers.
There are several standard approaches to analyzing the adaptive immune response present in the tumor microenvironment and tumor-draining lymph node. IHC and immunofluorescence can be performed, although specific antibodies often require frozen tissue; similarly, flow cytometry can be performed on frozen tissue. One of the benefits of the approach described here is the potential to perform the analysis on formalin-fixed paraffin-embedded tissue as methods for mRNA-seq from these samples have been demonstrated (59), and will continue to improve and become standardized. As this is available on a substantially greater number of patients compared with frozen tissue, this approach could allow for a much larger group of patients to be analyzed for the presence of adaptive immune signatures. Indeed, there is a large and growing number of human tissue samples with available mRNA-seq data; through the methods described here, these samples may be analyzed for antigen-directed B-cell responses.
Given our data that the presence of B-cells in the tumor microenvironment in patients with basal-like and HER2-enriched breast cancers and immunoreactive ovarian cancer is predictive of outcome, the role that endogenous B-cells play at these sites of tumor growth is a critical question. Plasma cells could generate antitumor antibodies that could be important in the early control of the growth of breast cancer cells, but which may ultimately become lost during tumor progression. Alternatively, B-cells may function as antigen-presenting cells to activate tumor-specific T-cells, which in turn may be inhibited via immunosuppressive mechanisms such as the PD1/PDL1 interaction. Future work is needed to determine the mechanism by which B-cells affect tumor growth in these different molecular subtypes of cancer, and whether and how this could be harnessed to improve endogenous antitumor immune responses.
The most difficult breast tumors to treat clinically are often of the “triple-negative” class, defined as such by the lack of cell-surface expression of ER, progesterone receptor, and HER2 (60). The majority (60%–80%) of triple-negative tumors are basal-like (61) and, thus, the basal-like subtype represents a critical target for the development of novel therapeutics. The presence of BCR characteristics associated with overall survival and consistent with a productive antitumor endogenous B-cell response suggests that methods to enhance or induce antitumor B-cells in patients with basal-like breast cancer may be clinically efficacious.
Disclosure of Potential Conflicts of Interest
C.M. Perou is an employee of and has ownership interests (including patents) in Bioclassifier LLC and University Genomics. No potential conflicts of interest were disclosed by the other authors.
Authors' Contributions
Conception and design: M.D. Iglesia, B.G. Vincent, C.M. Perou, J.S. Serody
Development of methodology: M.D. Iglesia, B.G. Vincent, J.S. Parker, C.M. Perou
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): B.G. Vincent, K.A. Hoadley, L.A. Carey, C.M. Perou
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): M.D. Iglesia, B.G. Vincent, J.S. Parker, L.A. Carey, C.M. Perou, J.S. Serody
Writing, review, and/or revision of the manuscript: M.D. Iglesia, B.G. Vincent, K.A. Hoadley, L.A. Carey, C.M. Perou, J.S. Serody
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): C.M. Perou
Study supervision: C.M. Perou, J.S. Serody
Acknowledgments
The authors thank the TCGA Network for the mRNA-seq data for the breast and ovarian samples.
Grant Support
This work was supported by grants from the NCI Breast SPORE program grant P50-CA58223-09A1, The Cancer Genome Atlas (U24-CA143848-05), RO1-CA-138255, F31-CA177188-01, the Breast Cancer Research Foundation, and the University Cancer Research Fund.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.