Molecular characterization of pregnancy-associated breast cancer and insights on timing from GEICAM-EMBARCAM study

Introduction

Breast cancer (BC) is one of the most common malignancies in pregnancy, accounting for approximately 1 in 3000 pregnancies, and is the leading cause of cancer-related death among pregnant and breastfeeding women^1,2. Pregnancy-associated breast cancer (PABC) is mainly defined as a BC diagnosed during gestation or within the first years after childbirth, although the definition may vary depending on the length of this postpartum period^3,4. This malignancy in the prenatal and postnatal period is a challenging clinical situation to diagnose and manage, with the need to consider both the mother and the potential risks to the fetus. PABC is estimated to represent 0.2–3.8% of all BC cases, with the rate rising to 7% for women under 45 years of age^3,5. However, an increase in PABC cases has been reported in recent years and is expected to continue. This is likely because women tend to delay childbearing, which increases the risk of BC by about 5% per year of delay^5,6. Interestingly, the incidence of PABC is significantly higher in the postpartum period than during pregnancy, which may be due to pregnancy-related changes in the breast and diagnostic challenges that make the tumor difficult to detect during this period³.

It often affects younger women with a median age of 33 years at diagnosis, and there is significant concern regarding the particularly aggressive behavior of the disease. PABC tumors are typically negative for both estrogen/progesterone receptors (ER/PR) and the human epidermal growth factor receptor 2 (HER2), which results in a higher prevalence of triple-negative breast cancer (TNBC). Furthermore, they often exhibit less favorable tumor characteristics such as advanced TNM stages, high grade, larger tumor size, and a higher rate of lymph node involvement compared to non-pregnancy-associated breast cancer (non-PABC)^7,8,9. Overall, PABC is associated with a poorer prognosis compared to non-PABC, previously attributed to late diagnosis and advanced disease stage at presentation¹⁰. However, recent meta-analyses have reported that PABC exhibits a lower survival rate and a higher propensity for metastasis compared to non-PABC, even after adjusting for stage at diagnosis, particularly among those diagnosed shortly after birth^11,12,13. This finding suggests a potential influence of the mammary microenvironment during pregnancy and the reproductive history on the biology and prognosis of PABC¹⁴.

This hypothesis is supported by the fact that pregnancy triggers two potentially opposite effects on BC risk, both cancer-protective and cancer-promoting¹⁵. This dual effect is strongly influenced by the age at first pregnancy, with a more pronounced long-term protective effect for young mothers. Nevertheless, all mothers face a transient increase in the BC risk after a full-term pregnancy, but especially those who have their first child after the age of 35^16,17. Therefore, the period surrounding pregnancy may provide an opportunity for BC development due to the hormonal and immunological changes that occur during pregnancy, as well as the effects of breastfeeding and breast involution in the postpartum period, which could modify the tumor microenvironment and promote neoplastic formation and/or progression of pre-existing disease^18,19.

Although the definitive biological mechanisms driving PABC remain unclear, several hypotheses have been proposed about factors that could work together to contribute to PABC development. The unique hormonal milieu during pregnancy, marked by high levels of pregnancy hormones, may promote breast cell proliferation and potentially initiate oncogenic transformation^6,14. Furthermore, the cellular immunosuppression and immune tolerance, which prevent immune rejection of the fetus, may also allow mammary tumor cells to evade the immune system and proliferate without detection²⁰. Breast involution is a major physiological change where the mammary gland returns to a quiescent state resembling pre-pregnancy. This process involves apoptosis and clearance of the secretory mammary epithelium, creating a postpartum involuting microenvironment characterized by an initial inflammatory response followed by immunosuppression, which is known to be potentially pro-oncogenic^21,22. Indeed, several studies suggest that lack or short duration of lactation (<12 months) may confer an increased risk of developing aggressive tumors due to the accumulation of less differentiated cells that survive apoptosis during breast involution²³.

Given the complex nature of PABC, further research into the molecular mechanisms involved in its pathogenesis is essential to improve our understanding and management. Current diagnostic methods have limitations in the successful detection of BC in pregnant or breastfeeding women, likewise, treatment options are currently limited to existing guidelines for young women with BC without specific targeted therapies identified for the PABC population²⁴. Hence, it is crucial to uncover the molecular pathways and drivers of this neoplasm to explain its aggressive biological characteristics, and potential relationship with pregnancy, as well as to identify new biomarkers with preventive and clinical implications for PABC patients. However, the molecular and biological basis of PABC remains incomplete, as few studies have focused on the molecular characterization of PABC^25,26,27, and our knowledge of the genomic and immune profile of PABC tumors is still limited.

Therefore, this study aims to explore the molecular landscape of PABC by using molecular profiling approaches in order to identify key pathways and genes associated with PABC onset, considering the timing of diagnosis. Using the NanoString technology, we performed comparative gene expression analyses on formalin-fixed paraffin-embedded (FFPE) samples from young patients with BC diagnosed during gestation (PABC-GS), in the first year postpartum (PABC-PP), or not related to pregnancy (non-PABC) (Supplementary Fig. 1). In addition, we used STRING analysis to show the potential interaction network between the proteins encoded by the most significantly expressed genes in each patient group.

Results

Patients and tumor characteristics

Table 1 summarizes the main characteristics of the patients and tumors. A total of 59 BC patients were included in the present study, of which 55.9% (n = 33) were classified as PABC and 44.1% (n = 26) as non-PABC. Within the PABC group, 15 patients (45.5%) were diagnosed during gestation (PABC-GS subgroup), and the remaining 18 patients (54.5%) were diagnosed during the first year postpartum (PABC-PP subgroup). Notably, the proportion of patients with a family history of breast/ovarian cancer was significantly higher in the PABC than in the non-PABC group (57.6% vs 15.4%, p = 0.012).

Table 1 Clinicopathological characteristics of the study cohort

Full size table

According to the immunohistochemistry evaluation, the triple-negative subtype is prominent in PABC, with 48.5% of TN tumors, followed by 36.6% HR⁺/HER2⁻ and 15.5% of HER2⁺ tumors (p = 0.0595, Table 1). Regarding histological type, tumor size, histology grade, lymph node invasion, and metastasis, no statistical differences were observed between groups. But it is noteworthy that over half of PABC tumors (54.5%) were classified as high-grade (G3). Data on tumor-infiltrating lymphocytes (TILs) were only available for PABC patients, 51.5% of whom exhibited elevated TIL levels (above 10%). Notably, when comparing the subgroups, we found that 61.1% of PABC-PP tumors had high TIL scores, against 40.0% of PABC-GS tumors.

Basal-like is the predominant subtype in PABC patients

Regarding the classification of intrinsic subtypes by the PAM50 gene expression signature (Supplementary Table 1), we noticed that basal-like subtype was significantly more frequent in PABC patients compared to non-PABC group (48.48% vs 15.38%, p = 0.012) (Fig. 1a). While the luminal A subtype was the least common subtype in PABC group, it was the most common in the non-PABC group (12.12% vs 34.62%, p = 0.0578) (Fig. 1a). Upon analysis of the diagnosis timing, both PABC-GS and PABC-PP subgroups demonstrated similar proportions of basal-like tumors (46.67% and 50%, respectively). Notably, this subtype was also the most prevalent within each subgroup (Fig. 1b).

Molecular characterization of pregnancy-associated breast cancer and insights on timing from GEICAM-EMBARCAM study — **Fig. 1: Intrinsic subtypes classification of samples.**

Principal component analysis (PCA) score plot revealed distinct clusters separating most of the PABC and non-PABC samples (Fig. 1c). In addition, this analysis demonstrated a clear separation of basal-like samples from other subtypes, with the majority belonging to the PABC group (Fig. 1c).

PABC patients have a distinctive gene expression pattern

To identify differentially expressed genes (DEGs) in the PABC group, we performed DEG analysis using a threshold of absolute fold change (FC) ≥ 1.5 and a false discovery rate (FDR) of <0.05. 73 DEGs were identified in PABC, including 32 genes that were significantly upregulated and 41 genes that were significantly downregulated (Fig. 2a). To examine the differences in expression between both groups (PABC vs non-PABC group), we conducted unsupervised hierarchical clustering analysis using the 73 DEGs previously identified, revealing a predominant clustering pattern for both groups (Fig. 2b).

**Fig. 2: DEGs between PABC and non-PABC.**

The top ten most significant DEGs between both groups (FDR < 0.01) were genes involved in several biological processes including DNA damage repair, proliferation or adhesion, and migration (Supplementary Table 2). Remarkably, four of the top 5 upregulated genes in PABC samples (DEPDC1, CCNA2, CDKN3, and FAM83D), exhibiting nearly a two-fold increase compared to non-PABC samples, play a role in the DNA damage repair pathway (Supplementary Table 2). Except for the PSAT1 gene, which is related to the triple-negative biology and whose expression in PABC is almost 4-fold higher than in non-PABC (Supplementary Table 2). Furthermore, there was significant up-regulation of these genes in both PABC subgroups (GS and PP subgroups) compared to the non-PABC group. Likewise, downregulated genes also exhibited a similar and lower expression level in both PABC subgroups compared to non-PABC (FDR < 0.05, Fig. 2c).

Furthermore, after adjusting DEG analysis for intrinsic subtype, more than 60% of the genes initially identified remained significant (p-value < 0.05), with CDKN3 prominently highlighted (Supplementary Fig. 2a). The DNA damage repair pathway was the most enriched signature even after this adjustment in pathway scoring analysis (PSA) (Supplementary Fig. 2b). Indeed when we re-evaluated the key DNA repair-related genes (DEPDC1, CCNA2, CDKN3, and FAM83D) in the unadjusted DEG analysis within each molecular subtype group (LumA, LumB, HER2+, Basal-like), we noted that these genes showed consistently higher expression in PABC samples across most subtypes, providing evidence that DNA damage repair pathway appears to be enriched in PABC independently of intrinsic subtype (Supplementary Fig. 2c).

Given the distinct transcriptomic signature for PABC tumors, we determine whether these 73 DEGs represent a coordinated change in tumor biology. Search tool for the retrieval of interacting genes (STRING) clustering analysis generates a network and predicts the interactions between the proteins encoded by this set of genes, providing an overall view of the multiple cellular functions in which these related genes participate. Thus, the functional enrichment analysis revealed a dominant biological cluster of 29 genes involved in the regulation of the cell cycle (FDR < 2.9e-12), with over half of these genes closely clustered together (Fig. 2d) Notably, the top 5 genes with higher interconnectivity within the constructed PPI network were MYC, FOXM1, CCNA2, CDC20, and CDKN3; all of the above hub genes were upregulated in PABC, except PTEN, which was the only one downregulated (Fig. 2d).

To get an overview of the biological processes in which these DEGs were involved, we performed a PSA on PABC samples compared with non-PABC. We observed enrichment in pathways associated with DNA damage repair, proliferation, and molecular subtype signatures (Fig. 3a). When examining the differences between samples within PABC subgroups (PABC-GS and PABC-PP), PABC-PP samples consistently showed slightly higher scores for DNA repair and proliferation signatures compared to PABC-GS samples (Fig. 3b).

Considering the potential significance of these signatures in PABC tumors, we sought to identify these biological pathways’ key genes. We observed that a substantial number of DEGs in PABC were part of these pathways, with 22 genes belonging to the DNA repair pathway (18 upregulated and 4 downregulated) and 20 genes to the proliferation pathway (11 upregulated and 9 downregulated) (Fig. 4a–d). Notably, this includes relevant proliferation-related genes such as the MYC oncogene and the tumor suppressor gene PTEN, as well as the oncogenic transcription factor FOXM1 involved in DNA repair (Supplementary Table 3). Regarding the subtypes signature, 9 DEGs were found in PABC (8 upregulated and 1 downregulated) (Fig. 4e, f), most of them were basal-like-related genes according to PAM50 subtyping (CDC20, KIF2C, MYC, PTTG1, TYMS, and UBE2C) (Supplementary Table 3).

**Fig. 4: DEGs for the most relevant signatures in PABC.**

Additionally, a heat map was represented for each signature using these DEGs in order to analyze the gene expression pattern by a group of samples, and particularly the DNA damage repair heat map showed a clear and strong overexpression pattern in PABC for most of those genes (Fig. 4b).

Specific differential genes for PABC-GS and PABC-PP subgroups

As PABC patients appear to have different biological features depending on whether the BC occurs during or after pregnancy^4,11, we focused on identifying genes that might show differential expression (DE) either during gestation (PABC-GS) or in the immediate postpartum period (PABC-PP).

We found 50 DEGs specific for PABC-GS, 20 DEGs specific for PABC-PP, and ten genes common to both subgroups (Fig. 5a–c). In the heatmaps performed with the DEGs, it is evident that PABC-GS has a larger number of significant genes than PABC-PP (60 vs 30 genes), however, clustering by sample type is still evident in both heatmaps (Fig. 5d, e). Among the most characteristic genes of each subgroup with the highest degree of significant DE, AGTR1 was highly downregulated in PABC-GS while SMAD5 was significantly downregulated in PABC-PP (Fig. 5b, c). The matched 10 DEGs to both subgroups, FGF10, CCL21, MAP2K4, and TCEAL1 (downregulated genes), as well as CDKN3, DEPDC1 and FAM83D (upregulated genes), were identified as distinctive genes for all PABC patients regardless of the time of BC onset (Fig. 5a).

**Fig. 5: Specific DEGs for each subgroup of PABC vs non-PABC tumors.**

PABC-PP tumors have an enrichment of immune infiltration signature

We performed DE analysis to identify significant differences in gene expression patterns between PABC-PP vs PABC-GS subgroups. This analysis identified 71 DEGs (with an absolute FC ≥ 1.5 and a p-value < 0.05), with 38 upregulated and 33 downregulated genes in the PABC-PP subgroup (Fig. 6a). The clustering heat map did not show a clear separation of the PABC subgroups using these DEGs. Nonetheless, we did observe a clustering among several samples according to their subgroup and a pattern of overexpression of certain genes in PABC-PP is evident (Fig. 6b). Top10 DEGs between the two subgroups (p-value < 0.01) are listed in the Supplementary Table 4 with box plots showing the level of expression of these genes in each sample by subgroup (Fig. 6c). Among the most upregulated genes, AGTR1 stands out, whose expression was highly increased in some PABC-PP samples compared to PABC-GS samples. In contrast, CDH2 and CLDN1, which play a role in adhesion and migration, exhibit consistent downregulation in PABC-PP samples compared to PABC-GS (Supplementary Table 4).

**Fig. 6: DEGs between PABC-PP and PABC-GS.**

STRING clustering analysis was also performed but with the identified 71 DEGs between the subgroups (PABC-PP vs PABC-GS) to generate a gene interaction network for a global view of the interactions and functions among these genes. The functional enrichment analysis generated a statistically significant biological cluster of 28 genes associated with the immune response (Fig. 6d, FDR < 3.21e-06). It is noteworthy that more of the 40% of DEGs (30/71) between subgroups are related to this process, of which 25 are upregulated in PABC-PP. Notably, of the top 5 genes with higher interconnectivity in the constructed PPI network, CD44, STAT1, CCR2, and FOXP3 were identified as hub genes upregulated in PABC and involved in the immune response (Fig. 6d).

We then conducted a PSA to provide a thorough understanding of the enriched biological processes between subgroups. Although not as pronounced, the PABC-PP samples showed a high score for the immune infiltration signature, while the PABC-GS samples displayed a high score for the epithelial-mesenchymal transition (EMT) and adhesion-migration pathway (Fig. 7a). Additionally, as immune infiltration was found to score highest in the PABC-PP subgroup, we examined the volcano plot of the DE analysis from PABC-PP vs PABC-GS, highlighting the genes involved in this signature and we also conducted a heatmap of these genes. As a result, although there was no clear clustering by sample subgroup in the heatmap, a tendency for up-regulation of immune infiltration-related genes in PABC-PP samples was observed in the volcano plot (Fig. 7b). Of the 32 genes in the immune infiltration signature, 28 were upregulated in the PABC-PP subgroup. Notably, the immune checkpoint gene PDCD1/PD-1 was identified as one of the most significant genes in this set (|FC| ≥ 1.5 and p-value < 0.05) (Fig. 7b). In summary, the PSA revealed that the immune infiltration signature was enriched in the PABC-PP subset. This was consistent with the upregulation of individual genes observed in the DE analyses and the immune response cluster identified in the STRING network.

**Fig. 7: Pathway and cell-type profiling analysis between subgroups.**

Finally, given the potential importance of the immune infiltrate in PABC-PP, we conducted a cell-type profiling analysis (CPA) to examine transcriptomic changes of specific immune cell types by measuring the abundance of selected marker gene transcripts across subgroups. Despite limited statistical significance, the quantification of cell populations based on gene expression profiles revealed an increased score of B cells, neutrophils, regulatory T lymphocytes, and macrophages in PABC-PP compared to PABC-GS (Fig. 7c).

Discussion

Despite significant advances in early detection and understanding of the molecular basis of BC, deficits remain for those that occur during or shortly after pregnancy. In this study, we analyzed the clinicopathological data and the specific gene expression profile of PABC, as well as the distinctive patterns in relation to the time of BC diagnosis. The results revealed molecular features that are unique to PABC, especially in tumors diagnosed shortly after childbirth, highlighting the need to differentiate between two PABC subsets based on diagnosis timing.

Our study showed a higher prevalence of the basal-like subtype in PABC compared to non-PABC. Moreover, PSA revealed a strong association between PABC tumors and BC molecular subtyping by gold-standard PAM50 signature. Basal-like tumors are often classified as TNBC because most are typically negative for ER, PR, and HER2; and both subtypes are associated with aggressive pathologic features and fewer therapeutic options^28,29. Certainly, our research demonstrated a high level of agreement between the triple-negative and basal-like subtypes assessed by IHC and the PAM50 assay, with 48.5% of PABC tumors being classified as triple-negative and 48.48% as basal-like. Other studies also confirmed this finding among PABC while the luminal A subtype was uncommon^8,30. TNBC subtype, which is reported to account for 50–60% of all PABC cases, has been associated with several molecular alterations^18,31. In our work, Phosphoserine aminotransferase 1 (PSAT1) was identified as one of the most significantly up-regulated genes in PABC with an expression nearly 4-fold higher than in non-PABC. Indeed, a recent study has linked the high expression of PSAT1 to the TNBC clinical grade, suggesting that multiple pro-tumourigenic pathways may enhance PSAT1 expression, thereby promoting the migration and invasion potential of TNBC³².

Some investigations have reported that the protective effect of pregnancy against BC varies by molecular subtypes. In fact, parity and early age at first birth have been associated with a marked reduction in the luminal BC risk, whereas neither parity nor age at first birth protects against TNBC^33,34,35. Elevated hormone levels during pregnancy are thought to decrease ER expression levels, which may contribute to suppressing ER-positive tumors and increase the development of ER-negative tumors in PABC^35,36. This could explain the low frequency of luminal PABC tumors observed in our study. Currently, as women tend to delay pregnancy, breast tissue is more likely to have accumulated cells with mutations and abnormal cells with the potential for malignant transformation. This risk is further compounded by a shorter duration of breastfeeding, which may lead to an undifferentiated progenitor cell population with a higher potential for carcinogenesis^37,38,39. As a result, there is a greater risk of developing an undifferentiated tumor phenotype, such as TNBC/basal-like BC. The higher frequency of these tumors in PABC is particularly relevant as it may contribute to their poorer overall prognosis and impact clinical management. Thus, it should be considered when evaluating optimal treatment strategies for these patients.

Our gene expression analysis supports the notion that pregnancy strongly influences BC transcriptome profile. Notably, our results are consistent with previous studies suggesting that PABC may be a distinctive and challenging entity of BC, exhibiting unique molecular signatures that differ from other BC in young women^25,26,40. In this regard, of the 73 DEGs identified in PABC, STRING functional analysis revealed a statistically significant cluster of 29 genes involved in cell-cycle regulation, in agreement with a previous study that observed an enrichment of cell cycle-related genes in PABC epithelia and stroma²⁶. Furthermore, pathway analysis revealed that PABC tumors were characterized by an enrichment of genes associated with relevant biological processes such as DNA damage repair and cell proliferation. In fact, four of the five most significantly up-regulated genes in PABC are related to DNA repair (DEPDC1, CCNA2, CDKN3, and FAM83D); highlighting cyclin A2 (CCNA2) which is known to be involved in cell cycle control and whose dysregulation appears to be closely associated with chromosomal instability and tumor proliferation in several cancers; particularly there is evidence of overexpression in the TNBC⁴¹. Nonetheless, we have evidence that these key genes and the DNA damage repair signature appear to be enriched in PABC independently of the intrinsic subtype, suggesting that they are distinctive features of PABC rather than merely a reflection of TNBC characteristics.

Of note, PABC also shows a dysregulated expression of other relevant cancer-related genes involved in cell cycle, proliferation, and DNA repair. These include the oncogenes MYC and FOXM1, as well as the tumor suppressor gene PTEN, which play key roles in breast progression and tumorigenesis^42,43. It would be interesting to further investigate the role of these genes in PABC, as well as to identify them as potential biomarkers for diagnosis and/or treatment. Taken as a whole, these results suggest that PABC tumors may exhibit increased cell cycle activity, heightened proliferative capacity, and aberrant DNA damage repair, which could lead to faster tumor growth, and a higher risk of metastasis. These findings could provide a biological explanation for the more aggressive behavior of PABC compared to non-PABC and could potentially lead to new treatment options for PABC patients targeting DNA repair or cell cycle checkpoints. Furthermore, we observed that a higher proportion of PABC patients had a positive family history of breast/ovarian cancer. This suggests that a heritable component may be present in PABC^44,45,46. However, this study lacks enough germline mutation data. Further research into germline mutations in well-known cancer genes, such as BRCA1/2, and other potential genes is required to unravel their role and contribution to PABC. This is necessary to understand how they may influence clinical and molecular characteristics, potentially leading to differences in tumor biology and treatment response in PABC patients.

Finally, we independently analyzed the gene expression profile of PABC patients diagnosed during pregnancy (PABC-GS) and those diagnosed within 12 months postpartum (PABC-PP), given their distinct biological features that may impact the BC transcriptomic profile⁴. Thus, our results could indicate that a BC arising in the postpartum period is associated with an increase of immune infiltration genes, unlike BC which develops during gestation⁴⁷. It is also noteworthy that PABC-PP exhibits increased expression of potentially relevant cancer targets involved in the immune response, such as the programmed cell death protein 1 (PDCD1/PD-1). Weaning-induced breast involution provides an active immune milieu characterized by massive epithelial cell death, stromal remodeling, and infiltration of immune cells with immunosuppressive features which has been demonstrated to promote BC development and metastasis^48,49. In fact, postpartum involution shares many attributes with wound healing, including up-regulation of genes involved in immune cell infiltration, which creates a pro-tumorigenic immune microenvironment that has been shown to persist beyond the timeframe of postpartum mammary gland remodeling^21,22,48. Additionally, our analysis between the PABC-PP and PABC-GS subgroups revealed altered expression of genes related to adhesion-migration and EMT, which may reflect changes in the mammary stroma and tissue remodeling due to the involution proces^50,51. Therefore, our findings indirectly support the hypothesis that postpartum breast involution could impact the tumor immune milieu and enhance tumor growth and metastasis, facilitating the dissemination of pre-existing tumor cells in an immune-evasive microenvironment, and thus contributing to the pronounced worse outcomes associated with a postpartum BC diagnosis^6,52,53.

Furthermore, our results suggest an immune phenotype in PABC-PP characterized by an increased abundance of several cell types, including neutrophils, regulatory T lymphocytes, macrophages, and B cells. Consistent with these findings, when comparing patients´ characteristics, we observed a greater number of tumors from the PABC-PP subgroup with high levels of TILs (>10%) compared to PABC-GS. It has been reported that tumor-infiltrating immune cells play a dual role in the onset and progression of cancer: they can protect against tumor progression by killing neoplastic cells, but they can also shape tumor immunogenicity and establish an immunosuppressive environment that facilitates the dissemination of tumor cells⁵⁴. Evidence of macrophage, neutrophil, and lymphocyte activation has been found in the involuting mammary gland⁵⁵. Regarding the specific role of immune cell populations present in mammary tumor microenvironments, regulatory T cells (Tregs) and M2-type tumor-associated macrophages (M2-TAMs) are considered to exert an immunosuppressive action that can drive tumor growth⁵⁴. Tregs are a distinct specialized subpopulation of T cells characterized by the expression of the nuclear transcription factor Forkhead box P3 (FoxP3), which appears to have a major role in disrupting the immune control of cancer^56,57. Several studies have shown that along with M2-TAMs, they promote tumor growth by suppressing the immune system through the secretion of cytokines such as IL-10 and TGF-β, leading to a poor prognosis for BC patients⁵⁸. In fact, an increased influx of Tregs FoxP3+ and macrophages with T cell suppressive function have already been observed in the involuting mammary glands of rodents and BC patients^48,51. In the case of B-cells, both pro- and anti-tumorigenic features were attributed, as BC cells can induce a regulatory phenotype in B-cells which support tumor progression by inhibiting Th1 mechanisms and enhancing Treg responses⁵⁹. Our work found that B-cells were the main cell type highlighted in PABC-PP, which is in line with an interesting recent study in which an increased plasma B-cell infiltration in post-weaning PABC tissue was correlated with the poorest outcomes⁶⁰.

Taken together, these findings collectively suggest that pregnancy-related immunological changes may exert a substantial influence on PABC, and given the pivotal role of the immunosuppressive cell influx in postpartum BC, this underscores the need to further investigate their tumor immune microenvironment and its impact on BC development and progression. An additional conserved mechanism of immune suppression between pregnancy, mammary gland involution and malignancy is altered immune checkpoint molecules such as PD-1 and PD-L1. It has recently been reported that in cases of postpartum BC, there is an increase in PD-1 expression on T cells and that anti-PD-1 treatment is effective in reversing involution-associated tumor growth⁶¹. This knowledge could help to explore the potential benefits of using immunotherapy for these patients, such as targeting tumor-associated macrophages or immune checkpoint pathways like the PD-1/PD-L1 axis.

Despite the lack of statistical power between subgroups, most likely due to the small sample size caused by the difficulty in recruiting PABC patients, these preliminary insights are highly relevant and warrant additional studies in a larger cohort to validate these findings and improve our understanding of BC biology during different stages of pregnancy, lactation, involution and (early and distant) postpartum. Of note, we need further confirmation in a large cohort of patients, in order to draw robust conclusions about the broader molecular landscape of BC subtypes. Additionally, longer follow-up and investigation of parity impact, the specific timing of diagnosis during pregnancy, and the duration of breastfeeding, are required to better elucidate the aggressive and recurrence risk associated with PABC.

Currently, against the traditional definition of PABC, there is a growing recognition that a BC diagnosed during gestation should be considered a separate and distinct clinical entity from a BC discovered in the postpartum period (which recent data suggest can extend to 5–10 years after birth)⁴. Our research supports this distinction, but also demonstrates that they are molecularly distinct entities, where differences at gene expression level are evident even in tumors detected in the immediate postpartum period (one year after childbirth), compared to tumors diagnosed during pregnancy.

In conclusion, our study provides relevant findings on PABC which displays a unique differential gene expression pattern linked to precise biological pathways and a network of protein interactions that could potentially contribute significantly to PABC aggressiveness. Furthermore, it became clear that it is crucial to distinguish between BC that develops specifically during gestation and that which arises during the early postpartum period, justifying the need for future studies to investigate their disparities and the possible clinical implications.

Methods

Study design: patients and sample collection

Building on the research of the GEICAM/2012-03 study⁶², the Registry Study of Pregnancy and Breast Cancer (EMBARCAM GEICAM/2017-07 study) (ClinicalTrials.gov identifier: NCT04603820) is a multicenter, observational, ambispective study enrolling BC patients, including those diagnosed during the gestation and up to one year postpartum.

The present analysis involved an age-matched subset cohort of PABC patients from both the GEICAM/2017-07 study and GEICAM/2012-03 studies, and non-PABC patients from the EpiGEICAM study (control cohort)^63,64. The inclusion criteria for PABC patients were: female sex, a new BC diagnosis during pregnancy or within one year postpartum, and age ≤ 45 years at diagnosis. Non-PABC patients were: female sex, a new BC diagnosis outside of pregnancy or one year postpartum, and aged ≤ 45 years. A total of 66 archival FFPE tumor tissue specimens of these BC patients diagnosed between 1992 and 2015 were collected. H&E slides corresponding to the FFPE tumor tissue samples of the patients were reviewed by an expert BC pathologist to confirm the diagnosis and determine the tumor surface area. Clinicopathological data for all included patients were also collected. Informed consent was obtained from each patient for the collection of clinicopathological and biological data following the ethical guidelines of the Declaration of Helsinki and good clinical practice guidelines.

The study cohort initially consisted of 66 samples from BC patients, but seven were excluded because the samples did not pass the quality analysis. This resulted in a total of 59 patients with samples available for gene expression analysis. The patients were classified according to parity as (i) PABC, if diagnosed during pregnancy or within 12 months of delivery; (ii) non-PABC. The PABC group was further divided into two subgroups according to the specific time at BC diagnosis: during gestation (PABC-GS) or in the first year postpartum (PABC-PP).

RNA isolation and nCounter gene expression profiling assay

Based on the size of the patient’s breast tumor surface, a total of 3–5 tissue sections, 10 µm thick, were obtained for RNA isolation using the RNeasy FFPE Isolation Kit (Qiagen, Hilden, Germany), following the manufacturer’s protocol. RNA yield and purity were assessed using NanoDrop™ 2000/Spectrophotometer 2000c (Thermo Fisher Scientific, Waltham, MA, USA) and Invitrogen™ Qubit® 3.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). RNA quality control was performed using an Agilent 4200 TapeStation system (Agilent Technologies, Santa Clara, CA, USA).

Tumor RNA sample gene expression was analyzed using the nCounter® Breast Cancer 360^™ V2 Panel (BC360^™ NanoString Technologies, Seattle, WA, USA). This gene expression platform covers 23 key BC pathways and processes, including the BC intrinsic subtypes (Luminal A, Luminal B, HER2-Enriched, Basal-like) by PAM50® signature⁶⁵. This panel consists of 758 gene-specific probe pairs for the targets, 18 housekeeping genes for normalization, 6 exogenous positive controls, and 8 exogenous negative controls. In a single reaction, probes were hybridized in solution with 150–500 ng of total RNA overnight at 65 °C following the instructions and kits provided by NanoString. Digital images were processed within the nCounter Digital Analyzer and gene counts from all the 66 samples were obtained for bioinformatics analysis. A quality control (QC) of raw data was performed using nSolver™ 4.0 Analysis Software (NanoString Technologies). Most samples (89.4% of the original cohort) satisfied quality criteria metrics by NanoString, and only seven samples were excluded due to analytical failure.

Raw transcriptome data from the 59 samples were background threshold-corrected using negative probes with a mean minus two standard deviations. Then, gene expression values were subjected to housekeeping gene normalization by calculating the geometric mean using the geNorm algorithm. Normalized data were log2-transformed for the following analyses conducted using the nCounter Advanced Analysis 2.0 and ROSALIND® platform (NanoString Technologies).

Differential gene expression analysis

Data were subjected to unsupervised hierarchical clustering. Expression values were normalized using z-scores and visualized in heat maps. Additionally, PCA was used to identify expression patterns between the different study groups. DE analysis was performed to identify the most significant DEGs between groups and subgroups. FC and p-values were calculated using the log-linear negative binomial model (fast method in the nCounter® Advanced Analysis)⁶⁶. p-values were adjusted using the Benjamini-Hochberg method to estimate the FDR. The clustering of genes for the final heatmap of DEGs was performed using the Partitioning Around Medoids method.

Pathway and immune cell population analysis

For PSA, each sample gene expression profile was condensed into a small set of pathway scores using the first principal component, which represented the average expression change (log2 scale) for the associated genes⁶⁷. We then performed a gene set analysis (GSA), which overlays DE data for sets of genes grouped by biological function, to identify the DEGs of the enriched pathways. DGE and PSA adjusted for intrinsic subtype as a confounding factor were also analyzed. We also performed a CPA, which uses the expression levels of cell type-specific marker genes to measure the relative abundance of immune cell populations among different groups of samples⁶⁸.

Protein–protein interaction (PPI) network construction

The interaction network of the protein products from the DEGs identified in the DE analysis of the PABC group and subgroups was explored by bioinformatic analysis using the STRING (version 11.5) database⁶⁹. We set the strength of the interactions at >0.4 as the threshold for analysis. The PPI network represented the genes as nodes and the interactions as edges. Nodes with a high degree of connectivity were called hub genes. We analyzed the functional annotations of this set of significant genes provided by STRING.

Statistical analysis

Statistical analysis of the clinicopathological data and PAM50 subtypes distribution between groups was performed using GraphPad Prism 8 (version 3.5). Quantitative variables were presented using descriptive statistics while qualitative variables were presented as numbers and percentages. The comparison between groups with qualitative data was performed by using Chi-square and Fisher’s exact test, whereas quantitative data were analyzed by the Mann-Whitney test. p-values < 0.05 were considered statistically significant.

Gene expression data were analyzed using the nCounter Advanced Analysis Software 2.0 and ROSALIND® platform (NanoString Technologies). In the DE analysis, p-values were corrected for multiple testing using the Benjamini-Hochberg FDR adjustment⁷⁰. To select statistically significant DEGs, a threshold of absolute FC ≥ 1.5 and an FDR-adjusted p-value < 0.05 were applied, except for DE subgroup analyses with fewer samples where an unadjusted p-value < 0.05 was accepted.