Poor Outcome in Postpartum Breast Cancer Patients Is Associated with Distinct Molecular and Immunologic Features

In young women, ages 40 years or younger, breast cancer is associated with more proliferative disease, poorer prognosis, and increased mortality, making it the leading cause of cancer-related death in this age group (1, 2). Variations in the definition of pregnancy-associated breast cancer have led to conflicting results regarding the association between the exact timing of a breast cancer diagnosis in these young women and their prognosis (3, 4). When differentiating breast cancers diagnosed in the postpartum period (PP-BC) from breast cancers diagnosed during (Pr-BC) or outside (nulliparous, NP-BC) the context of a pregnancy, we (5, 6) and others (7–12) found that patients with PP-BC diagnosed up to 5 to 10 years postpartum had an approximately twofold increased risk of metastasis and/or death. This higher metastatic rate and increased mortality could not be fully attributed to differences in tumor characteristics nor intrinsic biological subtypes (13, 14).

Although the exact etiology underlying the poor prognosis of PP-BC remains poorly understood, the postpartum remodeling of the mammary gland to its pre-pregnant state, also known as mammary gland involution, is thought to contribute to the increased risk for metastasis and death. Evidence from preclinical animal models point to the activation of wound-healing–like programs, with a characteristic initial inflammatory response followed by an immunosuppressive phase (4, 15, 16). In rodent models of mammary gland involution, enhanced tumor growth and increased motility and invasion of cancer cells have been observed (17, 18), which correlated with changes in the type and proportion of immune cells in the involuting tumor immune microenvironment (TIME; ref. 19). In healthy women without breast cancer, it is thought that the remodeling phase of mammary gland involution occurs in a narrow timeframe after delivery, yet a deregulated immune profile, with similarities to microenvironments present during wound healing and tumor progression, can last up to several months or even years (20, 21). In fact, leukocytes and macrophages have been shown to persist in the postpartum gland beyond 18 months, and distinct immune signatures have been found to persist up to 10 years after delivery (21, 22). In women with PP-BC, no reports exist on studies specifically investigating postweaning cancerous breast tissue (22, 23). Molecular investigations in PP-BC tissue, focusing on women with either triple-negative or ER+ breast cancer and not considering outcome data or the lactational status of the mammary gland, found an upregulation of inflammation-associated gene signatures when compared with NP control patients (22–24). In line with these findings, enrichment in (CD8⁺) T cells, natural killer cells, and macrophages have been reported in the TIME of PP-BC tissue (20, 25, 26). However, information about the lactational behavior of these patients at the time of their breast cancer diagnosis and data on their clinical outcome are lacking.

To investigate whether distinct, prognostic molecular features common to all breast cancer molecular subtypes exist in postpartum tumor tissue, we collected treatment-naïve breast tumor tissue from a subset of age, grade, and molecular-subtype–matched patients with PP-BC, Pr-BC, and NP-BC from a previous large-scale retrospective cohort study, thereby ensuring that all breast cancer molecular subtypes were represented (6). In addition, we gathered detailed data on lactation behavior of patients with PP-BC, to delineate the start of the postweaning period as a surrogate for involution onset. Using unbiased transcriptome analyses, IHC and multiplex immunofluorescence (mIF) analyses, we identified immune components in the TIME of postweaning PP-BC (PP-BC_PW) tissue that appear to be correlated to the particularly poor prognosis in these patients. In particular, PP-BC_PW was shown to comprise a heterogenous entity in which patients, having decreased gene expression levels of Ig and B-cell–related genes while showing an increased plasma B-cell infiltration in their tumor tissue, had the poorest prognosis. Although tumor tissue of patients with PP-BC_PW was also characterized by an increased presence of CD8⁺ T cells, this was not linked to prognostic differences.

A total of 230 archival formalin-fixed paraffin-embedded (FFPE) tumor tissue specimens of young women ages 25 to 40 years with primary invasive breast cancer diagnosed between January 2005 and December 2017 were collected. Patient-, therapy-, and tumor-related characteristics were collected. Patients that received neoadjuvant chemotherapy were excluded from FFPE biopsy selection (a consort diagram explaining the selection of patients and samples is given in Fig. 1; details on host- and tumor-related prognostic parameters are summarized in Supplementary Table S1). FFPE tissue sections from lumpectomy or mastectomy were collected from biopsies taken for diagnostic purposes, with the tumor zone subsequently marked by means of hematoxylin and eosin (H&E) scoring by an expert breast cancer pathologist. Depending on the correlation between a patient's cancer diagnosis and pregnancy and/or lactation history, patients were classified as (i) PP-BCPW, if diagnosed within 2 years postweaning; (ii) PP-BCDL, if diagnosed during lactation; (iii) Pr-BC, if diagnosed during pregnancy; or (iv) NP-BC, if nulliparous (Fig. 1). Patients were matched for age, grade, and surrogate molecular subtype and only included if they had follow-up for at least 20 months. ER, PR, and HER-2 status were evaluated using IHC according to ASCO/CAP guidelines (27, 28). Additional in situ hybridization techniques were used to confirm HER-2 gene amplification according to each participating center's guidelines. Surrogate molecular subtype was classified as luminal A-like (ER positive, HER-2 negative, grades 1 and 2), luminal B-like (ER positive, HER-2 negative, and grade 3), luminal HER-2 (ER positive, HER-2 positive, and any grade), HER-2-like (ER negative, HER-2 positive, and any grade), or triple-negative breast cancer (TNBC: ER negative, PR negative, HER-2 negative, and any grade). The study was approved by Ethics Committee Research UZ/KU Leuven (S/25470). All retrospective medical data/biospecimen studies at the respective institutes have been executed pursuant to local legislation and international standards. Patients consented to the use of their personal data and biospecimens in research. Hence, the procedures comply both with (inter)national legislative and ethical standards.

Total RNA from the tumor area was isolated using the High Pure FFPET RNA Isolation Kit (Roche Life Science). RNA quality (DV₂₀₀ values) and quantity were assessed on an Agilent 2100 Bioanalyzer using the Agilent RNA 6000 Nano Chip. Biopsies that were taken longer than 15 years ago were excluded to guarantee good RNA quality and quantity (DV₂₀₀ > 200 nucleotides). Using 100 ng of total RNA from each sample (DV₂₀₀ > 30%), libraries were prepared with the Illumina TruSeq RNA Access Library Prep Kit according to the manufacturer's protocol (Illumina). Barcoded libraries were 4-plex pooled before HiSeq4000 (Illumina) sequencing with 50-bp single-end reads at a depth of 30 million reads. FASTQ files were aligned and quantified using Salmon v0.14.2 (Salmon, RRID:SCR_017036) with the GRCh38 (release 94) reference transcriptome (29). Sequence quality and alignment metrics were assessed and collected via FastQC (v.0.11.9) and MultiQC (v1.8). Samples were excluded when they had a library size <5 × 10^E6 reads, when <75% of reads mapped to the reference transcriptome, or when they had >5% of overrepresented sequences. Technical replicated were combined using DESeq2. Each sample was subsequently assigned a PAM50 classification (R package genefu, v2.8.0; ref. 30). There were three samples with high incongruity between surrogate molecular subtype and mRNA receptor expression (normalized count > 10), which were removed from further analyses.

Read count data were processed by R packages tximport (v1.18.0) and analyzed by DESeq2 (DESeq, RRID:SCR_000154, v.1.30.1). Low count filtering was applied to include only genes with a cpm of 1 in at least n of samples, where n was the smallest group of replicates (n = 51 samples out of 211 samples; ref. 31). Principal component analyses (PCA) and Kruskal–Wallis rank-sum tests were performed to investigate the association of PAM50 and sequencing preparation batch on DEG (Supplementary Fig. S1). As we aimed to identify genes and pathways that are specifically affected in PP-BC_PW tissue, yet are common to all breast cancer subtypes, we controlled for the identified influence of PAM50 and library preparation prep in all differentially expressed genes (DEG) analyses by including them as covariates. Pairwise comparisons (PP-BC_PW vs. NP-BC, PP-BC_PW vs. Pr-BC, PP-BC_PW vs. PP-BC_DL) and aggregate comparisons (PP-BC_PW vs. rest, i.e., a combination of NP-BC, Pr-BC, and PP-BC_DL) were performed. Fold changes were shrunk using Bayesian estimates from apeglm (32). A cut-off of absolute log2-fold change of ≥0.5 and FDR ≤0.05 was applied to select DEGs. Gene set enrichment was performed using Flexgsea (v1.3) with 1,000 permutations, using the hallmark, GO-BP, canonical pathway, transcription factor, and chemical-genetic perturbation gene sets from MSigDB (v7.0; ref. 33). FlexGSEA: Flexible Gene Set Enrichment Analysis (v1.3; Bismeijer and Kim, 2019); available from https://doi.org/10.5281/zenodo.2616660). The weighted Kolmogorov–Smirnov statistic was used to calculate enrichment scores. The STRING proteomic database was used to further visualize networks and interactions of proteins (34).

To estimate the abundance of immune cells from bulk RNA-sequencing data, the web platform for CibersortX was used to calculate the fractions (sample size/population size) of 22 human immune cell subsets (LM22) defined in the CibersortX package for each bulk RNA-sequencing sample (35). CibersortX was run in relative mode with batch correction (B-mode) enabled and quantile normalization disabled, with 100 permutations. A fragments-per-kilobase million (FPKM) expression matrix produced via DESeq2 was used as input. Tumors were assigned a “high” or “low” score based on the CibersortX fractions for each immune cell subset above or on/below the median, respectively. Tcr Receptor Utilities for Solid Tissue (TRUST) was used to analyse B-cell receptor sequences by BCR CDR3s per kilo BCR reads (CPK) in each sample (36). Complete CDR3 sequence was defined as CDR3 annotated with both V and J genes. In addition, antibody isotype abundance of immunoglobulin (Ig) heavy chains (IgA, IgD, IgG, and IgM) was determined using DEG data by adding the expression of heavy chains from the same isotype.

For mIF, the following antibodies were included: anti-CD3 (Thermo Fisher Scientific, Catalog No. RM-9107-S1, RRID:AB_149924, clone SP7, 1/400 dilution, 1 hour at RT), anti-FoxP3 (Abcam, Catalog No. ab20034clone 236A/E7, 1/50 dilution, 2 hours at RT), anti-CD20 (Dako, Catalog No. M0755, clone L26, 1/500 dilution, 1 hour at RT), anti-CD27 (Abcam, Catalog No. Ab131254, clone EPR8569, 1/500 dilution, 1 hour at RT), anti-CD8 (Dako, Catalog No. M7103, clone C8/144B, 1/100 dilution, 1 hour at RT), anti-CD138 (VWR, Catalog No. VWRKILM3825-C1, clone B-A38, 1/100 dilution, 1 hour at RT), and anti-PanCK (Abcam, Catalog No. Ab27988, clone AE1/AE3, 1/100 dilution, 2 hours at RT). mIF staining was performed on a Ventana Discovery Ultra automated stainer (Ventana Medical Systems), using the Opal Polaris 7-Color Manual IHC Kit (50-Slide Kit, Akoya Biosciences, Catalog No. NEL861001KT) on a subset of 47 PP-BC_PW, 11 Pr-BC, and 60 NP-BC patients with remaining FFPE tissue specimens. Slides were imaged using the Vectra Polaris automated imaging system (Akoya Biosciences). Scans were made with the MOTiF protocol. Using the InForm software version 2.5.0, the MOTIF images were unmixed into eight channels: DAPI, OPAL480, OPAL520, OPAL570, OPAL620, OPAL690, OPAL780, and Auto Fluorescence and exported to a multilayered TIFF file, which were fused with HALO software (version 3.2, Indica Labs). Two pathologists independently marked the tumor regions of interest. A differentiation was made for the intratumoral region (interacting immune cells with the tumor) and the peritumoral region (area surrounding the tumor parenchyma). To prevent variation, consecutive slides were superimposed using the image registration tool with synchronized navigation. ROI were annotated using the brush and flood annotation tools. The random forest classifier is used to distinguish between tumor and stroma within the ROI. The Indica Labs Highplex FL v4.0.2 analysis algorithm was used for analysis using AI nuclei segmentation. All annotation layers were analyzed, and both the summary data and cell object data were exported in comma-separated value files using the export manager in HALO. Human tonsil FFPE tissues with and without primary antibody were used as positive and negative controls.

To improve marker specificity, cells that appeared positive for both CD3 and CD20, CD8 and CD20, or FoxP3 and CD20, were assumed to be positive only for the marker with the highest intensity. CD138 staining was found to show specific expression on epithelial cells and in the extracellular matrix and was removed from further analyses. Every cell was assigned a single cell type based on marker positivity. For the MPIF26 panel, the decreasing priority is FoxP3⁺, CD3⁺FoxP3⁺, CD20⁺, PanCK⁺. For the MPIF27 panel, the decreasing priority is CD3⁺CD8⁺, CD3⁺CD8⁻, CD20⁺, PanCK⁺. Spatial statistics were calculated within the tumor area drawn by a pathologist. Patterns of cell distribution were analyzed using the following spatial point statistics: density (number of positive cells per mm² in tumor area), clustering between immune cells within a cell type, clustering between immune cells of different cell types and clustering of immune cells towards PanCK⁺ cells (37). Clustering was calculated with the L and L-cross measures at 30 μm, which can be understood as the average density of cells around a cell per square mm, minus the average L-measure of random Poisson patterns with the same number of cells and observation window (38). Only immune cell distribution was randomized, to avoid detecting clustering of PanCK⁺ cells. Tumors were assigned a “high” or “low” score based on the density and clustering scores above or on/below the median, respectively.

Tumor-associated plasma cells (TAPC) and tumor-infiltrating lymphocytes (TIL) were scored on H&E whole tumor sections, based on published methods (39–41). A TAPC score 0 was given when no plasma cells were present; a score of 1 referred to the presence of at least one plasma cell; a score of 2 referred to the presence of at least one cluster of five plasma cells; and a score of 3 referred to more than two fusing clusters of plasma cells. Tumors were then classified as “low” or “high” when the TAPC score was respectively 0 to 1 or 2 to 3. In a subsequent single-plex IHC experiment, whole tissue slides were incubated with an antibody specific for CD38 (clone E7Z8C, Catalog No. 51000S, Cell Signaling Technology, 1/40,000 dilution, 1 hour at RT), as a potential marker for plasma cells. The assessment of the slides was done with Slide Score. The CD38 score was expressed as a percentage of positive cells of the total number of stromal cells both inside and outside the tumor area (42, 43). Normal glandular tissue was not considered. Tumors were then assigned a “high” or “low” score when the percentage of the intratumoral areas occupied by cells labelled for CD38 was above or on/below the median, respectively. CD38⁺ cells were further used as a marker for plasma B cells, supported by the positive association with TAPC score, although we acknowledge its lack of specificity. The cut-off median percentages used were compatible with accepted clinical pathologic practices. Pathologists were blinded to the clinicopathologic information.

Differences in the distributions of the categorical characteristics were compared using χ² analyses. For pairwise comparisons, an overall P value, OR, and 95% confidence intervals (CI) were determined. Continuous variables were compared by means of one-way ANOVA, (pairwise) Wilcoxon rank-sum tests, and/or Kruskal–Wallis tests. Correlation between ranked parameters was investigated based on Kendall or Spearman correlations. Survival rates were calculated from the time of diagnosis to the time of death of any cause for overall survival (OS) and to the date of first systemic metastasis for distant recurrence-free survival (DRS). OS and DRS univariate and multivariate analyses were performed using Cox regression on TMM-log2 normalized counts and Kaplan–Meier curves (44). TMM normalization was chosen over size factor normalization to minimize outlier contribution. In the multivariate models including all patients, we adjusted for stage (accounting for tumor size and lymph node infiltration), surrogate molecular subtype, and therapy (Supplementary Table S1; Fig. 1). All statistical analyses were performed using R version 4.0.3.

Code for all analyses is publicly available at https://github.com/nki-ccb/ppbc/. Metadata, raw RNA-sequencing data, and spatial image data contain pseudonymized information and will be made available after completion of a data transfer agreement. Requests should be directed to F.A. (frederic.amant@uzleuven.be).

We collected tumor tissue from a subset (n = 230) of young women with primary invasive breast cancer, who were either diagnosed within 2 years postpartum (PP-BC, n = 91), during pregnancy (Pr-BC, n = 55), or never-been-pregnant (nulliparous, NP-BC, n = 84) and who were previously included in a large retrospective study evaluating clinical data (published in EJC by Lefrère and colleagues 2021 (ref. 6; Fig. 1; Supplementary Fig. S1). For all patients, clinicopathologic data were available and for PP-BC patients, breastfeeding history was used differentiate breast cancers identified during lactation (PP-BC_DL, n = 21) from those diagnosed postweaning (PP-BC_PW, n = 70). We delineated the postweaning period (never breastfed or after cessation of lactation) in patients with PP-BC as a surrogate for the mammary gland involution window, which has been indicated in preclinical animal models as being important in the stimulation of tumor growth and invasion of cancer cells in the mammary gland tumor microenvironment (TME). We first evaluated whether the poor cancer outcome, as previously identified for a large clinical PP-BC_PW patient cohort, was still identifiable for the patients in this selected group (6). Although we did not find statistically significant differences, patients with PP-BC_PW tended to have a two- to threefold increased risk for reduced survival rates (HR, 2.3–3.1-fold) and to develop metastatic disease more often (HR, 1.6–2.4-fold) compared with patients with PP-BC_DL, Pr-BC, and NP-BC, independent of differences in surrogate molecular subtype, tumor stage, and therapy characteristics [including surgery, (neo)adjuvant chemotherapy, radiotherapy, hormone therapy, and anti HER2-therapy; Supplementary Table S1 and Supplementary Fig. S2] and concurring with our previous publication (6).

To elucidate which biological mechanisms are associated with this poorer prognosis, we performed whole transcriptome sequencing analyses on tumor biopsy specimens to search for gene expression signatures that are specific to PP-BC_PW while being common to all breast cancer subtypes. As exploratory principal component analyses indicated that PAM50 subtype and sequencing batch contributed strongly to the variance in gene expression (Fig. 1; Supplementary Fig. S1), we controlled for these parameters in all downstream analyses. We compared the transcriptome in breast tumor tissue from patients with PP-BC_PW with the transcriptomes obtained from breast tumor samples from patients with Pr-BC, PP-BC_DL, and NP-BC in a series of pairwise and groupwise comparisons (Supplementary Fig. S3). Differential expression revealed the presence of 47 genes to be differentially expressed in PP-BC_PW tumor tissue when compared with all other patient groups, of which 35 were significantly upregulated (Fig. 2A; Supplementary Table S2). Upregulated genes included those involved in (mucin) glycoprotein biosynthesis processes (MUC2), cancer-testis antigens (POTEC and POTED), and breast cancer-related genes such as estrogen receptor 2 (ESR2, low expressed). Remarkably, there was a dominant cluster of Ig genes, consisting of almost half of the upregulated DEGs in PP-BC_PW tumor samples (Fig. 2B). This cluster, further referred to as “Ig gene signature,” consisted of 18 genes of the Ig superfamily and was expressed to a greater extent in PP-BC_PW tissue (Fig. 2C). STRING database analysis of the 47 DEGs confirmed the presence of a biological cluster of genes involved in the humoral immune response (Fig. 2D). Pathway enrichment analysis using Fisher exact testing indicated that 13 of the top 20 differentially expressed pathways in PP-BC_PW were related to the humoral immune response (Supplementary Fig. S4). Gene set enrichment analyses using gene ranking from DESeq2 with sample permutation (implemented via Flexgsea) did not further identify specific immunologic pathways differentially regulated in PP-BC_PW tissue. Transcription factor enrichment analyses did not identify significant differences in transcriptional regulators associated with PP-BC_PW compared with PP-BC_DL, Pr-BC, and NP-BC (lowest adjusted P value is FDR 0.56). Additional PAM50 subgroup analyses confirmed that PP-BC_PW-specific Ig gene signature upregulation was not driven by a specific PAM50 subtype (Supplementary Fig. S5). When applying K-means clustering to further identify patterns within the gene expression data of the 47 DEGs, the presence of three clusters was identified, which differed in their expression levels of the Ig gene signature (Fig. 2B). We found significantly more PP-BC_PW tumor samples with intermediate to high (81.8%) Ig gene signature expression levels compared with PP-BC_DL (44.4%), Pr-BC (58.7%), and NP-BC tissue (51.3%; P < 0.01; Fig. 2E). Interestingly, the genes belonging to the Ig gene signature were not differentially expressed in PP-BC_PW subgroups with less or greater than 6 months’ time interval between the start of the weaning period and their breast cancer diagnosis nor between PP-BC_PW subgroups with less than 1 month breastfeeding versus more than 1 month of breastfeeding prior to their diagnosis (Supplementary Fig. S6). These data indicate that the overexpression of this Ig gene signature did not differ based on the time period between the start of weaning nor on the duration of prior breastfeeding. When comparing the set of 47 PP-BC_PW-specific DEGs with published gene expression signatures from normal involuting mammary gland tissue, only two genes were found to be shared (i.e., CXCL14 and CELSR2; refs. 45, 46). Furthermore, hierarchical clustering and Spearman correlation analyses showed that the expression levels of well-known differentially expressed milk genes in our signature (i.e., LALBA, CSN1S1, CSN3, and CSN1S2AP) poorly correlated with the overexpression of genes from the Ig gene signature in PP-BC_PW tumor tissues, suggesting that the increased Ig gene signature in patients with PP-BC_PW is not driven by lingering milk production in normal postweaning mammary gland tissue (R² = 0.046; R = 0.22; P = 0.08; Supplementary Fig. S7). In conclusion, we observed an upregulation of the Ig gene signature in tumor tissue of patients with PP-BC_PW when compared with patients with PP-BC_DL, Pr-BC, and NP-BC, which cannot be fully attributed to normal lactational or other postweaning processes.

We next assessed whether the general overexpression of the Ig gene signature in PP-BC_PW tumor tissue was associated with patient outcome. We applied a cut-off of at least 15 patients per group for comparison analyses (ANOVA, Kruskal–Wallis, Wilcoxon tests) and 40 patients per group for inclusion in Cox regression analyses, to retain a statistical power of at least 70% (Fig. 1). Although no statistical significance was determined for patient groups below the cut-off value, we visualized the data from the comparison analyses to get additional insight in the distribution of the data. Patients with PP-BC_PW with low Ig gene expression levels (18.2%) were found to have a significantly increased risk for death (P = 0.04) and a nonsignificant increased risk to develop metastases (P = 0.078) compared with Pr-BC (41.2%) and NP-BC (48.7%) patients with low Ig gene expression levels; this was not the case among patients with PP-BC_PW with intermediate/high Ig gene expression levels (Fig. 2F; Supplementary Fig. S8). No similar association was found between Ig gene signature expression levels and prognosis, based on PAM50 status (Supplementary Fig. S9), further indicating that PAM50 is not associated with prognostic differences based on Ig gene signature expression levels. In parallel, the RNA transcriptome data were used to investigate via univariate Cox regression whether expression of specific genes was related to differential survival or metastatic rates in patients with PP-BC_PW versus the other patient groups. No genes were significantly associated with differences in survival rates among the groups. Overall, 95 genes were observed to be significantly correlated to metastasis (FDR < 0.05). Of the 59 genes that were significantly associated with an improved outcome, 20 were upregulated Ig and B-cell related genes (Supplementary Fig. S10A; Supplementary Table S3). From these 20, there were seven genes uniquely correlated with significant increased metastatic and/or death rates when expressed at low levels, specifically in PP-BC_PW, which included CD27, CD40LG, BAFF-R, SH2D1A, IGLV1–40, IGHV3–64D, and IGKV3–7 (Supplementary Fig. S10B). Together, these data suggest that a subgroup of patients with PP-BC_PW, that is, those with low Ig and B-cell–related gene expression in their breast tumor tissue, have a particular increased risk for death and metastatic disease compared with patients with Pr-BC and NP-BC as well as to patients with PP-BC_PW with high Ig gene signature expression levels.

We next aimed to examine whether the altered Ig gene signature expression levels were correlated with changes in (plasma) B-cell levels in breast tumor tissue of patients with PP-BC_PW and whether the presence of these cells also correlated with changes in survival rates. Exploratory CIBERSORTx analyses of the transcriptome data identified no differences in estimated fractions of naïve B cells, yet we found significantly decreased proportions of memory B cells (Wilcoxon rank-sum, P = 0.01) and significantly increased proportions of plasma B cells (Wilcoxon, P < 0.001) in PP-BC_PW samples compared with those of patients with NP-BC, but not to patients with PP-BC_DL and Pr-BC (Supplementary Figs. S11A and S11B). This was not correlated to differential outcomes when comparing patients with PP-BC_PW to NP-BC or Pr-BC (Supplementary Figs. S11C and S11D).

In parallel, we aimed to evaluate these findings at the cellular level, by means of mIF, IHC analysis, and microscopic evaluation of tumor biopsy slides. We used the mIF data on coexpression of CD20, CD3, and CD8 to identify—as non-plasma B cells (further referred to as CD20⁺ cells; Fig. 3A). Although CD20⁺ cell densities were not increased (Wilcoxon, P = 0.10; Fig. 3B), CD20⁺/CD20⁺ clustering (normalized cellular colocalization within a fixed radius of 0.30 mm) was significantly decreased (Wilcoxon, P = 0.03; Fig. 3D) in PP-BC_PW compared with NP-BC tissue, yet this difference was not associated with differences in prognosis (Fig. 3C–E; Supplementary Figs. S12A and S12B). Furthermore, no significant differences were found in clustering of CD20⁺ cells with tumor cells when comparing all patient groups (CD20⁺/PanCK⁺; Supplementary Fig. S13). Coexpression of CD20 and CD27 was used to exploratory examine differences in naïve (CD20⁺CD27⁻) and memory (CD20⁺CD27⁺) B-cell populations. Contrary to CibersortX findings, tissue from patients with PP-BC_PW did not display altered levels of CD20⁺CD27⁻ nor CD20⁺CD27⁺ cells compared with the other patient groups (Supplementary Fig. S14).

To evaluate the abundance of plasma B cells within the TIME, we resorted to assessment of tumor-associated plasma cells levels (TAPC) counting and single-plex IHC staining for CD38 in FFPE slides (Fig. 3F; Supplementary Table S4). TAPC can be readily identified in H&E tissue. They were graded based on their presence: “0” if no plasma cells were found, “1” for a few scattered plasma cells, “2” for a small cluster of five plasma cells, and “3” for a confluence of clusters of plasma cells present in the stroma surrounding the clustered tumor cells (Supplementary Fig. S15). PP-BC_PW tumor samples were found to have significantly higher TAPC counts (Kendall correlation P = 0.02; Wilcoxon P = 0.01) compared with both Pr-BC and NP-BC samples (Fig. 3G; Supplementary Fig. S16). In patients with PP-BC_PW, these high TAPC counts were found to associate with an increased risk for death (P = 0.03) and metastasis (P = 0.03) as opposed to patients with Pr-BC and NP-BC, in whom high counts were positively associated with improved prognosis (Fig. 3H; Supplementary Fig. S13C). In addition to TAPC assessment, intratumoral CD38⁺ cell densities were evaluated to further examine the presence of the plasma B-cell population within the TIME (Fig. 3I; Supplementary Table S4). CD38⁺ cell densities were significantly increased in PP-BC_PW tumor tissue compared with NP-BC tissue (P < 0.01; Fig. 3J). Again, patients with PP-BC_PW with a high intratumoral CD38⁺ density had a significantly increased risk for both death (P = 0.02) and metastasis (P = 0.01) compared with patients with NP-BC. No such outcome differences were observed among patients with low CD38⁺ densities (Fig. 3K; Supplementary Fig. S13D). Although CD38 is expressed on multiple immune cell types, it has the potential to be used as a marker for plasma B cells due to its high expression on this cell type (42). As there was a positive correlation between TAPC count and CD38⁺ intratumoral densities in the overall population (R² = 0.11; R = 0.33; P < 0.01; Supplementary Fig. S17), CD38⁺ can be interpreted as a plasma B-cell marker in our population, although caution remains warranted. Together these data indicate that, opposed to patients with Pr-BC and NP-BC, in patients with PP-BC_PW, an increased infiltration of TAPCs and CD38⁺ cells in the TIME was correlated to a significantly elevated risk for metastasis and death.

Given the low Ig gene signature expression levels and increased presence of plasma B cells in tumor tissue of PP-BC_PW, both parameters being associated with a poor prognosis in this patient group, we wondered whether changes in the antibody repertoire could explain this apparent contradiction. Further exploiting the RNA sequencing data, patients with PP-BC_PW were found to have significantly increased IgA and IgD gene expression levels but decreased IgG gene expression levels compared with patients with NP-BC (Wilcoxon, P values <0.01; Fig. 4A). IgD, but not IgA and IgG, gene expression levels were also significantly different in tumor tissue of PP-BC_PW compared with patients with Pr-BC (Wilcoxon, P < 0.01; Fig. 4A). IgM gene expression levels did not differ between patient groups. As IgA is the main Ig found in breast milk, we compared IgA gene expression levels with differentially expressed milk genes in our signature (i.e., LALBA, CSN1S1, CSN3, and CSN1S2AP). We found a positive correlation in all patient groups, with 19% to 35% of the variation explained through a correlation between IgA and the milk signature in PP-BC_DL, Pr-BC, and PP-BC_PW respectively, but only 9.8% to be explained within NP-BC. This indicates that high IgA expression is, at least partially, related to the lactating and/or postweaning processes in patients with PP-BC_PW (R² = 0.35; R = 0.59; P < 0.01; Supplementary Fig. S18). Increased IgG and decreased IgA gene expression levels tended to be correlated to an increased risk for death in patients with PP-BC_PW, although the difference did not reach statistical significance (P_IgG = 0.057; P_IgA = 0.052, Fig. 4B). No significant differences related to metastasis were observed. Also, no differences in neither IgD nor IgM gene expression levels were related to patient outcome (Fig. 4B; Supplementary Fig. S19). Only within the patient with PP-BC_PW group, low IgA and high IgG gene expression levels were significantly associated with high TAPC counts (R = −0.23, P = 0.02 and R = 0.36, P < 0.01, respectively; Fig. 4C and D; Supplementary Figs. S20–S21). Likewise, a comparable negative correlation for IgA (R = −0.24; P = 0.1) and a positive correlation for IgG (R = 0.25; P = 0.09) gene expression levels with CD38 intratumoral densities in PP-BC_PW was found, although these results were nonsignificant (Supplementary Fig. S22). Within PP-BC_PW, we observed a significant negative correlation between IgA and IgG gene expression levels (R² = 0.40; R = −0.63; P < 0.01; Supplementary Fig. S23), underscoring the inverse effects of IgA and IgG expression within this patient group. Finally, we could not find a significant difference in clonotype diversity between PP-BC_PW, PP-BC_DL, Pr-BC, and NP-BC tumor tissue as estimated with TRUST (Supplementary Fig. S24A; ref. 36), indicating no differences in immunoglobulin clonal expansion between patient groups nor any association with prognosis (Supplementary Fig. S24B).

In general, these data indicate that PP-BC_PW tumor tissue is characterized with increased IgA and decreased IgG gene expression levels. These data point to the existence of a heterogeneous PP-BC_PW population with an upregulated expression of the IgG gene and a decreased expression of IgA gene, correlated to increased TAPC counts, the latter being associated with a particularly poor prognosis.

Next, we investigated whether the differential association between plasma B-cell presence and prognosis in patients with PP-BC_PW compared with the other patient groups could be related to differences in interactions with other immune components. Histopathologic analyses of tumor biopsy slides did not identify any significant difference in the general TIL percentages between the patient groups (Fig. 5A). TIL percentages were found to positively correlate with CD4⁺ naïve T cells (R = 0.29; P = 1.5 × 10⁻⁰⁹), Ig gene signature levels (R = 0.3; P = 5.5 × 10⁻⁰⁸), and TAPC counts (R = 0.28; P = 4.2 × 10⁻⁰⁷) in all study groups (Supplementary Fig. S25). Exploratory CIBERSORTx analyses of the remaining immune cell subsets from the transcriptome data showed a significantly increased fraction of CD8⁺ cytotoxic T cells (Kruskal–Wallis, P = 0.02) in PP-BC_PW compared with PP-BC_DL, Pr-BC, and NP-BC tissues (Fig. 5B). Increased CD8⁺ fractions seemed to associate with a poorer prognosis in patients with PP-BC_PW, although this difference in OS/DRS was not statistically significant (log-rank, P_OS = 0.13 and P_DRS = 0.11; Fig. 5B; Supplementary Fig. S26A). No significant differences were observed between NP and postpartum patient groups among CIBERSORTx estimated fractions of natural killer, monocytes, macrophages, mast, nor dendritic immune cell subsets (Supplementary Fig. S27). This finding was further validated using our previously described mIF panels in which we stained for CD3⁺CD8⁺ T cells, further referred to as CD8⁺ T cells, and CD3⁺CD8⁻ T cells, further referred to as CD8⁻ T cells. mIF staining confirmed a significantly increased density of CD8⁺ T cells (number of cells/total tissue area; Wilcoxon, P = 0.005; Fig. 5C), yet decreased clustering of CD8⁺/CD8⁺ T cells (normalized cellular colocalization within a fixed radius of 0.30 μm; Wilcoxon, P = 0.02, Fig. 5D) in the TIME of patients with PP-BC_PW compared with what was seen in patients with NP-BC. Furthermore, compared with NP-BC tissue, PP-BC_PW tissue was characterized by decreased clustering between CD8⁺/CD8⁻ T cells (Wilcoxon, P = 0.01; Fig. 5E); between CD20⁺ B cells/CD8⁺ T cells (Wilcoxon, P = 0.01; Fig. 5F) and between CD20⁺ B cells/CD8⁻ T cells (Wilcoxon, P = 0.001; Fig. 5G). Although, in general, high clustering between these immune cells was associated with an increased risk for death and metastasis in especially patients with PP-BC_PW, these differences were not statistically significant (Fig. 5; Supplementary Fig. S26). No significant differences in clustering between CD8⁺ T cells/PanCK⁺ tumor cells nor CD8⁻ T cells/PanCK⁺ tumor cells were observed between patient groups (Supplementary Figs. S28A and S28B). Although clustering of CD8⁻/CD8⁻ T cells was decreased in PP-BC_PW compared with NP-BC, also no overall significant differences between patients’ groups were observed (Supplementary Fig. S28C). Together, these data demonstrate that the observed differences in CD8⁺ T-cell densities or distances to other immune cells in the TIME of patients with PP-BC_PW did not correlate with differential prognosis when comparing patient groups.

Here, we performed a comprehensive and large-scale analysis of PP-BC tissue, obtained from patients diagnosed with different molecular breast cancer subtypes, to specifically investigate if there is a distinct molecular signature in the postweaning mammary gland, as a surrogate for the postpartum involuting microenvironment, that is associated with the typical poor prognosis of this patient group. Interestingly, we identified decreased Ig gene signature expression levels and an increased tumor infiltration by plasma B cells in PP-BC_PW tumor tissue, being correlated with the poorest outcomes. This is the first time that an association with plasma B cells and breast cancer prognosis has been revealed in the context of PP-BC. These data may imply that plasma B cells could adopt protumoral roles in this PP-BC setting. At present, the factors defining the functional role of plasma B cells in cancer are not fully understood (42, 47, 48). Although there is much evidence regarding the significance of tumor-infiltrating T cells in breast cancer, little is known about the role of B cells and plasma B cells in breast cancer, and existing data are often conflicting (49–53). Local residing plasma B cells, typically occupying chronic inflammation sites in the TIME, have been found to produce high titers of tumor-specific antibodies that can drive antitumorigenic effects (54, 55), by promoting antibody-dependent cytotoxicity and phagocytosis (56), complement activation and enhancing antigen presentation by dendritic cells (57–60). Although antitumor effects of B cells are mainly attributed to the IgG isotype, the IgA isotype is more often associated with an immunosuppressive microenvironment in which B cells promote expansion of regulatory T cells and reduce cytolytic killing activity of T cells (55, 61). In our patient with Pr-BC and NP-BC groups, we observed such an association between high plasma B-cell densities, an upregulation of the Ig gene signature, increased IgG gene expression, decreased IgA gene expression levels, and improved prognosis. Although we found significant increased yet decreased gene expression levels of respectively IgA and IgG in PP-BC_PW compared with Pr-BC and NP-BC, no prognostic differences based on Ig-heavy subtype were observed between our patient cohorts. We did identify patients with PP-BC_PW with high plasma B-cell infiltration, previously correlated with poorer outcomes, to be skewed towards low IgA and high IgG gene expression levels. It has been shown that, dependent on the composition of the TIME, the phenotypes of B cells present and the antibodies they produce, certain types of tumor-infiltrating (plasma) B cells can exert protumor effects (55). In the PP-BC_PW setting, it might be possible that plasma B cells produce IgG antibodies that do not efficiently elicit T-cell responses, as has been reported for other tumor types (55, 62, 63). Another hypothesis that may explain the apparent contradictory correlation between the poor prognosis of PP-BC_PW, elevated plasma B-cell frequencies, and increased IgG gene expression levels in PP-BC_PW TIME and is that IgG gene expression is not B-cell derived but coming from cancer cells (47, 64, 65). Studies have shown that cancer-derived Ig share identical basic structures with B-cell–derived Ig but can have profound protumorigenic effects via different mechanisms, including promotion of tumor immune escape (66), inducing inflammation (67), and downregulation of the cytotoxic and natural killer cells (68). In contrast to B cells using the Ig transcription factor Oct-2, cancer cells tend to use Oct-1. Although no difference in DEG for Oct-1 nor Oct-2 were observed between study groups, high expression levels of the Oct-1 gene were significantly related with poor prognosis in our PP-BC_PW cohort (P = 0.01). Future IHC investigations on the presence and localization of immunoglobulin (subtypes) in PP-BC_PW tissue could help elucidating the cellular source of the upregulated IgG genes.

Apart from changes in the B-cell compartment, we also found an increased infiltration of CD8⁺ cytotoxic T cells in PP-BC_PW tumor tissue. This is in line with previous findings in biopsy tissue taken from patients with PP-BC, but for whom the diagnosis was not specifically delineated to the postweaning period or included only one specific breast cancer molecular subtype (20, 24–26). Yet, using patient outcome data, we were able to show that neither changes in cytotoxic CD8⁺ T-cell frequencies nor differences in colocalization with other lymphocytes in the TIME were correlated to the poor prognosis of patients with PP-BC_PW.

As mentioned, access to breast tumor tissue along with information on the patients’ lactational behavior and cancer outcome is a unique strength of this study, which allowed us to delineate the postweaning window and investigate the correlation between molecular alterations in postweaning breast cancer tissue and the patients' prognosis. Similar findings in all surrogate molecular subtypes suggests that the role of the plasma B-cell compartment in the prognosis of PP-BC_PW is a general phenomenon. As for some analyses, we lacked a sufficient number of patients in the Pr-BC and PP-BC_DL patient groups to reach significant power, future studies should focus on repeating similar analyses with higher patient numbers. Changes in Ig-heavy gene expression should be validated at the protein level. Despite initial inclusion of several Ig-heavy antibodies in our mIF panels, we were not able to optimize these panels accordingly due to a specific binding. Also, the immune milieu profiling by IHC and mIF was focused on a small subset of markers, with TAPC scoring and CD38⁺ staining being the sole markers to identify plasma cells. Although we observed a good correlation between TAPC and CD38 positivity, CD38 is not a specific marker for plasma B cells, and results should be interpreted with caution. Future research, where possible, focusing on prospective collection of fresh samples, larger sample sizes and single-cell (spatial) approaches that focus on individual B-cell subsets, using CD19, CD21, CD23, CD24, CD80, CD86, and TLR markers, can support further characterization of the full B-cell repertoire and their exact role in the PP-BC_PW TIME.

We identified PP-BC as a heterogenous entity and further distinguished a subgroup of patients with PP-BC_PW, characterized by changes in Ig gene expression and plasma B-cell levels, being associated with a particular poor outcome. Although we also found significant differences in presence and clustering of/with CD8⁺ T cells in the breast tumor tissue of patients with PP-BC_PW, these parameters had no prognostic value. The molecular and immunologic differences identified in this research when comparing tumor tissue from patients with PP-BC_PW with that from control patients may serve as an important starting point for additional biomarker research. In particular, further in-depth characterization of the role of B cells in the intratumoral context in PP-BC, may pave the way for personalized management strategies for these patients.

F.A. Peccatori reports personal fees from Roche Diagnostics, Ipsen, and Merck outside the submitted work. S. Loibl reports grants and other support from AZ, AbbVie, Amgen, DSI, Gileadalth, Celgene/BMS, Novartis, and Pfizer; grants from Molecular He; other support from Seagen, Sanofi, Relay, Olema, Eirgenix, Merck kG, Lilly, GSK, Pierre Fabre, Esai, MSD, and Incyte outside the submitted work; in addition, S. Loibl also has a patent for VM Scope with royalties paid, a patent for EP14153692.0 pending, a patent for EP21152186.9 pending, and a patent for EP15702464.7 pending. K.J. Jerzak reports personal fees from Amgen, AstraZeneca, Apo Biologix, Eli Lilly, Esai, Genomic Health, Gilead Sciences, Knight Therapeutics, Merck, Myriad Genetics Inc, Pfizer, Roche, Seagen, Novartis, and Viatris; and grants from AstraZeneca, Eli Lilly, and Seagen outside the submitted work. K.E. de Visser reports grants from Roche/Genentech, Dutch Cancer Society (KWF), Dutch Research Council (NWO), Oncode Institute, and ERC outside the submitted work; and other support from Macomics. L.F. Wessels reports grants from BMS outside the submitted work. No disclosures were reported by the other authors.

H. Lefrère: Data curation, software, formal analysis, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. K. Moore: Data curation, software, formal analysis, validation, visualization, writing–original draft, writing–review and editing. G. Floris: Funding acquisition, validation, investigation, methodology, writing–review and editing. J. Sanders: Investigation, methodology, writing–review and editing. I.M. Seignette: Investigation, writing–review and editing. T. Bismeijer: Software, formal analysis, visualization, writing–review and editing. D. Peters: Investigation, writing–review and editing. A. Broeks: Methodology, writing–review and editing. E. Hooijberg: Funding acquisition, validation, methodology, writing–review and editing. K. Van Calsteren: Resources, writing–review and editing. P. Neven: Resources, writing–review and editing. E. Warner: Resources, writing–review and editing. F.A. Peccatori: Resources, writing–review and editing. S. Loibl: Resources, writing–review and editing. C. Maggen: Resources, writing–review and editing. S.N. Han: Resources, writing–review and editing. K.J. Jerzak: Resources, writing– review and editing. D. Annibali: Writing—review and editing. D. Lambrechts: Conceptualization, software, supervision, funding acquisition, methodology, writing–review and editing. K.E. de Visser: Funding acquisition, validation, methodology, writing–review and editing. L. Wessels: Software, formal analysis, supervision, funding acquisition, methodology, writing–original draft, writing–review and editing. L. Lenaerts: Conceptualization, supervision, funding acquisition, methodology, writing–original draft, project administration, writing–review and editing. F. Amant: Conceptualization, supervision, funding acquisition, methodology, writing–original draft, project administration, writing–review and editing.

This work was supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme “[647047 to F. Amant].” This work was supported by the KWF kankerbestrijding, the Dutch Cancer Society “[11132 to F. Amant].” This work was supported by the Fonds Wetenschappelijk Onderzoek, the Flemish Research Foundation or “FWO” “[G0A9219N to F. Amant].” Prof. F. Amant is senior investigator for the FWO. This work was supported by the Kom op tegen Kanker (Stand up to Cancer), the Flemish Cancer Society “[3M150537 to H. Lefrère].” Prof. Giuseppe Floris is recipient of a postdoctoral mandate from KOOR in UZ-Leuven. G. Floris is recipient of a postdoctoral mandate from the Klinische onderzoeks- en opleidingsraad (KOOR) of the University Hospitals Leuven. We thank all participating patients and all (para-)medical staff involved in registering cases in the INCIP database (see www.cancerinpregnancy.org). We would like to acknowledge the NKI-AVL Core Facility Molecular Pathology & Biobanking (CFMPB) for supplying NKI-AVL Biobank material and lab support. We would like to thank Michel Faas for his lab and image analysis support during his student rotation period. We acknowledge institutional funding from the KWF to the Netherlands Cancer Institute.