Oestrogen receptor phosphorylation profiles and in silico PAM50 subtyping reflect sexual dimorphism in breast cancer.

Introduction

Sexual dimorphism is increasingly recognised in biology and has been reported in cancer, with breast cancer (BC) showing the biggest dichotomy [1]. While BC is still rare in males, the numbers receiving a BC diagnosis have increased globally from 8,500 diagnoses in 1990 to 23,100 in 2017, with age‐standardised rates in 100,000 person years of 0.46 and 0.61, respectively [2]. There is evidence that increased obesity aligns with this increase [3]. That we are an ageing population may also contribute, particularly to reduced overall survival (OS), but this does not account for the lack of improvement in BC‐specific survival observed for men diagnosed with stage III or IV BC over the past 30 years [4].

A striking feature of male BC is its almost universal positivity for both oestrogen receptors (ER) ‐α and ‐β [5,6]. ERα, used to determine clinical outcome in female BC, is expressed in almost 90% of all male BCs in contrast to the 60–70% expression observed in females and confirmed in two of the largest studies on male BC [6,7]. Male BC is thus considered, somewhat counterintuitively, to be ER driven with the most recent ASCO guidelines advocating tamoxifen as first line endocrine therapy [8]. When bound to ligands like tamoxifen, ERα becomes phosphorylated. Many phosphorylation sites exist on ERα and S104, S118, S167, S282, S294, T311, and S559 are the best characterised in female BC [9]. The impact of ERα phosphorylation in male BC is unknown.

Reflective of its propensity for high ERα‐positivity, male BC is predominantly of Luminal phenotype [6,10]. This is a well‐recognised intrinsic subtype of female BC [11] and from this, PAM50, a 50‐gene subtype predictor, including ERα, PR (progesterone receptor), and HER2 (human epidermal growth factor receptor 2), was developed [12]. Marketed as Prosigna™, the commercial PAM50 assay can estimate the risk of future recurrence in hormone receptor‐positive female BC [13]. Further stratification of female BC into 10 subgroups has been achieved by the METABRIC group [14]. Similar attempts have been made to identify specific subtypes of male BC. Comparative genomic hybridisation (CGH) of tumour DNA from 56 cases applied to a BAC array identified male‐simple and male‐complex stratification. Cases in the smaller male‐simple cluster were associated with better prognosis: small size, reduced S‐phase and fewer genomic alterations while male‐complex aligned closely with the Luminal B subtype identified in females [15]. Two subgroups were identified in the same cohort of male BC patients using transcriptomics and an unsupervised clustering approach: Luminal M1 and M2 which were not represented by the intrinsic subgroups of female BC [16]. Combining the CGH and expression data derived from these two studies identifiedTAF4andCD164as the only two genes which were common between male and female BC, indicating that candidate driver genes of BC are not the same between sexes [17]. Two clusters were also observed through whole transcriptome analysis of male BCs with germline mutations in the most relevant BC susceptibility genes (BRCA1/2, PALB2,RAD50,RAD51D) with higher expression of immune response genes and high scores of gene‐expression signatures associated with aggressive phenotype and reduced survival [18]. Other groups have also shown sex‐specific differences with reduced frequency of 16q copy number in male compared to female BC [19], and less frequent promoter hypermethylation [20]. While informative, these studies were limited in terms of case numbers, a common issue when studying male BC.

The aim of this study was twofold. First, we employed an integrated bioinformatics approach using published bulk RNA‐seq datasets to agnostically investigate the intrinsic subtypes of male BC and their transcriptomic characteristics. Second, based on these findings, we used immunohistochemistry to examine the expression and prognostic value of ERα phosphorylated at S104, S118, S167, and S294 in male BC.

Materials and methods

A pipeline of the integrated bioinformatics approach is shown in Figure1with specific details presented below.

Experimental pipeline for cluster discovery. Following preprocessing, PAM50 gene signatures for each case were derived from individual datasets which were retained in the metadata. These datasets were merged, batch‐adjusted, and split randomly into discovery (n= 97) and validation (n= 98) sets. Cluster discovery and differential gene expression (DGE) analysis were performed separately on the datasets. GEX module signature score calculation and GSEA were performed on common differentially expressed genes in the discovery and validation sets. Created in BioRender. Chatterji, S. (2026).https://BioRender.com/ujqg65o

Dataset selection and pre‐processing

Datasets were identified from The Cancer Genome Atlas (TCGA) [21], National Centre for Biotechnology Information (NCBI) – Gene Expression Omnibus (GEO) [22], and ArrayExpress [23], and by contacting authors of previously published studies [18,24]. Five non‐RNA‐seq high throughput datasets were excluded [17,25,26,27,28] with only bulk RNA‐seq datasets sequenced on Illumina platforms included to minimise technical bias. Details of the datasets used are provided (supplementary material, TableS1). Sample identifiers for each patient included in datasets obtained from NCBI‐GEO and TCGA are given in supplementary material, TableS2. From this, 195 cases were suitable for analysis. For consistency, all were mapped to Entrez gene IDs using the biomaRt package [29,30] in R. In case of multiple mappings to the same gene ID, the averages of the raw counts were used. When data were available in the form of FASTQ files, alignment to the reference human genome (Ensembl GRCh38.p14) with STAR v2.7.11 [31] and raw counts were obtained using HTSeq v2.0.4 [32]. These steps were performed using Python 3.11. For visualisation, HUGO Gene Nomenclature Committee (HGNC) gene symbols were used. Manual pre‐processing of each dataset was done separately to account for the differences in the formats of the available processed data as detailed inSupplementary Materials and Methods. PAM50 intrinsic subtypes were determined by applying thegenefupackage [33] to each dataset separately. The resulting subtypes were retained in the clinical information associated with each patient for further analysis.

Integration, batch correction, and patient randomisation

Pre‐processed datasets were combined using the merge function in R. Only those genes that were common to all datasets were included for further analysis (n= 15,674). Since datasets were combined from multiple sources, it was necessary to minimise batch effects due to technical differences while retaining the biological signals. For batch correction, the ComBat‐seq method [34] available from the sva package in R [35] was applied to the filtered read counts. Batch correction was confirmed by principal component analysis (PCA) both before and after application of ComBat‐seq (supplementary material, FigureS1). The batch adjusted dataset was randomly sampled using the sample function in R into discovery (n= 97) and validation sets (n= 98). Subsequent analyses were performed independently on these datasets.

Unsupervised hierarchical clustering and differential gene expression (DGE) analysis

Generation of gene expression matrices and subsequent DGE analysis for the discovery and validation sets were performed with the DESeq2 package [36]. The variance stabilising transformation method was used to normalise gene expression matrices using the vst package [37]. Subgroups with common gene expression patterns were then estimated based on the 2000 most variable transcripts in each dataset using unsupervised hierarchical clustering following Ward's minimum variance method [38]. The top 2,000 variable genes were analysed with the assumption that these probably represented the major sources of variation in the data, in accordance with other BC studies [39,40]. Stability of the resulting clusters was evaluated by multiscale bootstrap resampling using pvclust [41], where the number of bootstrapped datasets was set to 10,000, as per package recommendation to minimise error. Veracity of the estimated clusters was confirmed using clusterpval [42].

Differential expression between the clusters was determined based on the log2 (fold change) ≤1.0 (downregulation) and ≥1.0 (upregulation), with an adjustedpvalue of <0.05. False discovery rate was controlled for using the Benjamini Hochberg method.

Gene expression (GEX) module signatures

GEX signatures associated with genes involved in proliferation (AURKA), apoptosis (CASP3), ERα signalling (ESR1), HER2 signalling (ERBB2), invasion and metastasis (PLAU), immune response (STAT1), and angiogenesis (VEGF) were calculated for each patient as follows:

$$\text{Module Score} = \frac{\Sigma_{i} w_{i} x_{i}}{\Sigma_{i} \mid w_{i} \mid}$$

For each gene 𝑖, its expression was denoted as 𝑥, while the value of 𝑤 was +1 or −1 depending on whether gene 𝑖 was up or downregulated in the female BC study that defined these signatures [43]. These module scores were calculated for both the discovery and validation sets, and the difference in the scores between clusters was determined using the Wilcoxon rank‐sum test.

Gene set enrichment analysis (GSEA)

GSEA was performed using Enrichr [44] to determine the biological relevance of differentially expressed genes (DEGs) between the clusters. This was applied to DEGs that were common to discovery and validation sets. Gene ontology (GO) enrichment analysis was performed based on biological process (BP), cellular component (CC), and molecular function (MF) [45]. Pathway analysis was performed by interrogating the Reactome Pathway database [46], the National Centre for Advancing Translational Sciences (NCATS) BioPlanet database [47], the Kyoto Encyclopaedia of Genes and Genomes (KEGG) [48], WikiPathways [49], and Molecular Signatures Database (MSigDb) Hallmarks [50]. Significance of the key terms from each database was determined using Benjamini Hochberg FDR‐adjustedpvalues of <0.05. Enrichment scores were calculated as the −log10 of the adjustedpvalue. The key terms were sorted according to decreasing enrichment scores and the top 10 terms from each database that reached statistical significance were reported for each cluster.

Immunohistochemistry (IHC)

A cohort of 506 male BCs was represented on tissue microarrays (TMAs) as described previously [28]. Three separate cores were sampled from different areas of each tumour case and placed into TMAs. Clinical characteristics are shown in supplementary material, TableS3. Ethical approval was obtained from the Leeds (West) Research Ethics Committee (06/Q125/156), Greater Glasgow Health Board (Network Approval: TR000269), and the NHS Grampian Tissue Bank Committee (Network Approval: TR000269). Four phosphorylated ER epitopes, which had been validated previously in female BC [9] were detected immunohistochemically (Leica BOND III Protocol F; BOND Refine Detection Kit). Details of antibodies used are provided in supplementary material, TableS4.

Following IHC, TMA slides were scanned (ZEISS AxioScan Z1 slide scanner 20×) generating whole slide images (WSIs). These were visualised using QuPath [51] open‐source image analysis software (v0.4.4510), which was also used for all subsequent image analysis steps, including stain vector normalisation, TMA dearraying, tissue detection, annotation of tumour and stromal regions, and stain quantification. The Allred system was used for nuclear staining and percentage of positive tumour cells for cytoplasmic staining. Scoring and assessment was done blinded to patient characteristics and outcome.

Statistics

All statistical analyses and data visualisation were performed in R 4.3.0. OS was defined as the time between diagnosis to death from any cause or the date of last follow‐up. Log‐rank tests and Cox proportional hazards models based on 5‐ and 10‐year OS were used for prognostic assessment. Correlation between analytical and clinicopathological variables was assessed using Fisher's exact test with the Bonferroni method to correct for multiple testing. Apvalue <0.05 after adjustment for multiple testing was considered significant.

Results

Unsupervised hierarchical clustering reveals two stable clusters in male BC

Two stable clusters were identified in the discovery and validation sets (supplementary material, FigureS2), termed Cluster C1 and C2. We tested the null hypothesis that there was no difference in gene expression levels between the two clusters of patients. For both our discovery and validation datasets, we identified numerous significant gene expression differences between Clusters C1 and C2 (p< 0.001). To rule out potential confounding from specific studies in our meta‐analysis, we examined whether cluster separation in either the discovery or validation data was driven by any of the source datasets by Fisher's exact test between the clusters and the source datasets. No such evidence was found as no significant correlation was seen between these with any one cluster in either the discovery or validation sets (p= 0.07;p= 0.269) (supplementary material, TableS5). Cluster C1 was smaller in both sets: discovery set (n= 14; approximately unbiased [AU] probability = 76%), validation set (n= 22, AU probability = 85%). Cluster C2 had 83 patients in the discovery set and 76 in the validation set (AU probabilities 77%, 86%, respectively). These checks confirmed data fidelity.

Differential gene expression (DGE) analysis

DGE was analysed between Clusters C1 and C2 in the discovery and validation sets. This revealed 1,168 DEGs in the discovery set (409 upregulated and 759 downregulated genes in Cluster C1, with the opposite pattern of regulatory state seen in Cluster C2). In the validation set, 1,294 genes were differentially expressed (367 upregulated and 926 downregulated genes in Cluster C1, with Cluster C2 exhibiting the opposite pattern of regulation; supplementary material, FigureS3). There were 603 common differentially expressed genes between the discovery and validation data sets, of which 230 were downregulated and 373 were upregulated in Cluster C1. No significant difference in the mean log2 (fold change) was observed between the discovery and validation sets for the common genes (p= 0.24). Expression patterns of the 604 common DEGs in the discovery and validation sets are shown in Figure2. Upregulation ofSLC39A6,NAT1, andPGR, which showed the highest differential expression in both discovery and validation sets, was observed in Cluster C1.

Expression of the 604 common DEGs in the discovery and validation sets and their relationship with PAM50. Heatmaps generated from the discovery (A) and validation (B) sets. Hormone receptor status and predicted PAM50 phenotype are reported for each case.SLC39A6,NAT1, andPGRwere upregulated in Cluster C1 of both the discovery and validation sets (black vertical line).

Survival was assessed in a subset of the cohort that had complete OS information. This was available for 108 of 195 cases analysed. While cases from Cluster C2 had a trend towards poorer survival, no significant differences were observed between cases from Clusters C1 and C2 at either 5‐ (p= 0.08) or 10‐year OS (p= 0.29). This is shown in supplementary material, FigureS4.

To resolve whether the two clusters of male BC aligned with previous female BC subtypes, datasets were assessed further for enrichment of PAM50 features. Their enrichment in each cluster was examined with Fisher's exact test. Cluster C1 was significantly correlated with the Luminal A subtype in both discovery and validation sets (adjustedp= 0.009 andp= 0.01, respectively; supplementary material, TableS6). Interestingly, when stratified by ERα and PR status, no enrichment was observed with the predicted PAM50 classification in either dataset (supplementary material, TableS7).

Gene set enrichment analysis (GSEA)

We applied GSEA to examine how predefined groups of genes were altered in male BC. GSEA scores and their associated genes are shown in Table1. GO and pathway terms enriched in Cluster C1 included genes involved in basement membrane (GO: CC), oxidative stress response (WikiPathways), arylesterase activity (GO: MF), KRAS signalling activation and late oestrogen response (both MSigDb Hallmarks). MAPK signalling upregulation was seen in multiple platforms (NCATS BioPlanet, KEGG, WikiPathways). In Cluster C2, genes involved in extracellular matrix (ECM) organisation and regulation were upregulated. Terms relating to this included ECM‐receptor interaction (NCATS BioPlanet, KEGG), ECM organisation (GO: Biological Process, Reactome Pathway), collagen‐containing ECM (GO: CC), IL‐1 regulation of ECM and TGF‐β regulation of ECM (both NCATS BioPlanet). Genes involved in the Hedgehog (NCATS BioPlanet) and Hippo signalling pathways (KEGG), oncostatin M (NCATS BioPlanet), regulation of cell population proliferation (GO: Biological Process), hormone activity, receptor ligand activity, cytokine activity (all GO: MF), epithelial to mesenchymal transition, downregulated genes due to KRAS signalling activation, early oestrogen response (all MSigDb Hallmarks), and differentiation of white and brown adipocytes (WikiPathways) were also enriched in Cluster C2. Significantly enriched GO terms and pathways that had high gene counts are shown in Table1.

Table: Enrichment terms and their associated genes, along with their regulation in Clusters 1 and 2, sorted in descending order of enrichment score

GEX module signature analysis

To identify biological pathways involved in male BC, GEX module scores were calculated for each patient. These were compared between Clusters C1 and C2 for both discovery and validation sets. In both datasets, Cluster C1 had significantly higher scores for the ERα signalling (ESR1signature) module (p= 7.3e‐05 for discovery set;p= 1.8e‐06 for validation set). On the other hand, significantly higher scores for proliferation (AURKAsignature) and invasion and metastasis (PLAUsignature) modules were seen for Cluster C2 in both discovery (p= 4.2e‐04 andp= 4.7e‐03, respectively) and validation sets (p= 1.5e‐04 andp= 4.7e‐06, respectively). Violin plots showing the above distributions are given in Figure3.

Distribution of GEX module scores that had significant differences between Clusters C1 and C2 in both discovery and validation sets. Violin plots showing the GEX modules that had significantly different scores between Clusters 1 and 2 in both the discovery and validation sets. These were (A and B) proliferation (AURKAsignature), (C and D) invasion and metastasis (PLAUsignature), and (E and F) ERα signalling (ESR1signature) modules.

No other GEX modules examined showed significant differences between Clusters C1 and C2 consistently in both datasets. The HER2 signalling module (ERBB2signature) had higher scores in Cluster C2 (p= 0.03), but only in the discovery set. Similarly, the angiogenesis module (VEGF signature) had higher scores in Cluster C2 (p= 0.03), but only in the validation set. No differences were seen in the apoptosis (CASP3 signature) and immune response (STAT1 signature) modules between the clusters in either dataset. Violin plots showing the above comparisons are shown in supplementary material, FigureS5.

Immunohistochemical expression and prognostic role of phosphorylated ERα

GSEA data showed upregulation of oestrogen response genes and MAPK signalling pathways. The MAPK pathway significantly influences ERα function, primarily through the phosphorylation of serine residues. Since the presence of ERα is always consistently higher in male than female BC, we reasoned that this may be due to ERα phosphorylation, a post‐translational modification that alters ERα activity and which has not been explored in male BC. While S104, S118, S167, and S294 all displayed nuclear and cytoplasmic staining (supplementary material, FiguresS7–S9), only nuclear immunoreactivity, determined by the Allred system, was considered as this is its typical cellular location. Staining trends were comparable when manual scoring (Allred) was compared to QuPath‐generated values.

Due to missing data, clinicopathological information available was specific to each biomarker. Correlation between available clinicopathological variables and expression of each phosphorylated ERα was tested and is shown in Table2. Low S104 expression was significantly associated with high tumour grade (p= 0.02), which lost significance after Bonferroni adjustment. No associations were seen for S118. Low nuclear expression of S167 was associated with PR and AR negativity (p= 0.02 and 0.002, respectively). High nuclear S294 correlated significantly with ERα and AR negativity (p= 0.03 and 0.04, respectively) and high tumour grade (p= 0.04), which was lost upon Bonferroni correction. No other significant associations were observed. Prognostic values were investigated using multivariate analysis (Table3) for 5‐ and 10‐year OS; however, no significant associations were found.

Table: Clinicopathological associations of ERα phosphorylated at S104, S118, S167, and S294 according to cell location

Table: Multivariate analysis to assess the impact of ER phosphorylated at S104, S118, S167, and S294 on overall survival at 5 and 10 years

Discussion

Molecular subtypes were first described in female BC some 25 years ago [11,52] and have been gradually refined over the years. Transcriptomic profile comparisons between female and male BC have identified sex‐specific pathways related to ECM organisation, metabolism and protein translation [26,28]. Incomplete overlap of hormone receptor signalling has also been reported [6,24,53]. More recently, 30 male and 54 female BCs were profiled in a 90‐gene expression PCR‐based assay which showed a distinct separation of sexes. Four genes (PI15,AZGP1,PRRX1,AGR2) were upregulated, and five (PGR,SFRP1,PLA2G2A,S100A2,CHI3L1) were downregulated in males, with the sex chromosomesRPS4Y1 and XISTup‐ and downregulated, respectively [54].

Attempts to stratify male BC using genomic approaches identified a two‐cluster profile by two independent research groups, with significant differences in prognosis and regulation of pathways relating to proliferation, invasion and metastasis observed [16,18]. Since the current reporting framework for all BC is based on female BC and as sex‐specific discrepancies, notably the almost universal ER‐positivity in male BC, have been previously reported [6,7],a prioristratification based solely on clinical features is challenging. Given the rarity of male BC, this warrants an agnostic approach. To address this, previously published bulk RNA‐seq datasets were combined, re‐split randomly into discovery and validation sets, with clustering and gene expression analysis performed independently on both. This approach also identified two clusters (C1 and C2) in both the discovery and validation sets, with Cluster C2 harbouring a more aggressive phenotype, in terms of proliferation, invasion and metastasis.

Cluster C2 was significantly enriched for GO and pathway terms related to ECM regulation and organisation, with collagen‐containing ECM showing highest enrichment. An immunohistochemical study of male BC showed that strong expression of collagen IV in the stroma predicted shorter disease‐free survival [55]. Deregulation of the collagen‐containing microenvironment has been reported in several cancers with p53 and JAK–STAT oncogenic pathways implicated [56,57]. In pancreatic cancer, mutatedKRASproto‐oncogenes interacted with the EMT‐regulator Snail, resulting in enhanced collagen production [58]. Crosstalk between collagen and cancer cells can also be facilitated through TGF‐β/Smad signalling [59]. Interestingly, Cluster C2 was significantly enriched for genes downregulated due to KRAS signalling activation and those involved in EMT and TGF‐β mediated ECM regulation. It could be speculated that regulation of collagen‐containing ECM in some male BCs is similar to findings in pancreatic cancer where collagen production increased due to interactions between EMT‐regulators andKRAS[58,60].

In the context of ECM, upregulation of genes involved in the Hippo signalling pathway in Cluster C2 is notable. This concurs with our previous findings that expression of TAZ and YAP, which are also Hippo transducers, along with their targets CTGF and AXL were predictors of poor OS in male BC [61,62]. In agreement with this, a trend for poorer survival was observed in Cluster C2 for both 5‐ and 10‐year timepoints, which may potentially be significant in a larger and more balanced cohort. Thus, deregulation of collagen in the ECM may be an important mediator of male BC cancer progression. Detailed interrogation of the collagen‐containing ECM using surrogate biomarkers such asCOL2A1, COL4A1, COL11A1, andLOX(all significantly upregulated in Cluster C2) alongside morphometric analysis could provide insight into potential functional roles. This is currently under investigation in our group.

We found enrichment ofESR1signature (ERα signalling module) in Cluster C1 and correlation with Luminal A phenotype. Intriguingly, no significant associations were found between hormone receptor status in male BC and the predicted PAM50 subtypes, which show strong correlation in FBC [13]. PAM50 has been applied in other male BC datasets, where a statistical association between immunohistochemical and PAM50 subtyping was seen [63,64]. However, questions have been raised over the predictive ability of PAM50 in male BC since over half that were classified as Luminal A by immunohistochemistry were grouped differently by PAM50 [64]. It is noteworthy from these previous studies [63,64] that PAM50 subtypes were derived using the Prosigna™ PAM50 assay [65], rather than usinggenefupredictions conducted in this work. This may question the concordance between computation and assay‐derived phenotypes. These challenges have been recognised in BC [66]. Lack of association between the PAM50 intrinsic subtypes and hormone receptor status from ourgenefupredictions may have arisen from the high imbalance in ERα/PR expression groups. However, this is unavoidable in male BC due to its almost universal ER‐positivity [6,7]. It may also reflect the modest numbers of cases available. Nevertheless, our data showed that PAM50 used to stratify female BC does not apply in males.

While qualitative examination of the relative expression of the PAM50 genes in the Clusters C1 and C2 did not reveal any discernible patterns, consistent upregulation ofPGR,NAT1, andSLC39A6was observed in Cluster C1.NAT1upregulation in Cluster C1 was consistent with previous findings [16,18]. Differential expression ofPGRdid not match the consistent PR positivity across the study cohort, suggesting a non‐linear relationship between gene and protein expression.

GO terms related to oestrogen response were enriched in both Clusters C1 and C2, albeit through different gene sets. Cluster C1 was significantly enriched for late response pathways to 17β‐oestradiol (E2), while Cluster C2 was enriched for early response to E2. A meta‐analysis of oestrogen response in MCF‐7 cells showed that early response genes were related to cell growth and proliferation, while late response genes were related to cellular assembly and organisation, cell cycle, DNA replication, recombination and repair [67] – this was in agreement with our GEX module score findings as Cluster C2 was highly enriched forAURKA(proliferation) module score and exhibited a trend for poorer OS. This meta‐analysis also reported high regulation of the AhR (Aryl hydrocarbon Receptor) signalling pathway for both early and late response. The AhR pathway has been reported to inhibit ERα signalling in rat studies including direct inhibition by AhR/ARNT (Aryl hydrocarbon nuclear translocator) heterodimerisation, proteasomal degradation of ERα and increased CYP1A1 and CYP1B1 expression, inhibiting E2 synthesis [68]. Notably, Cluster C1 had significant enrichment for the GO term arylesterase activity. This may suggest increased inhibition of the ERα pathway in this cluster, with increased expression of the hormone receptor as a compensatory response.

We reasoned that despite overwhelming ERα positivity in male BC, the signalling pathway itself may be deregulated. Noting that our GSEA data showed upregulation of oestrogen response genes and MAPK signalling pathways in both clusters we hypothesised that this might impact on ER‐phosphorylation. Phosphorylation is fundamental in regulating steroid hormone receptor activity. ERα is phosphorylated at multiple sites, governed by kinase activity and phosphorylation at S104, S106, and S118 contributes to ERα activity [69]. ERα phosphorylation at S118 characterises an intact oestrogen‐dependent signalling pathway in BC and is associated with a better clinical outcome in female patients treated with tamoxifen [70,71]. Also in female patients, ERα phosphorylation at S167 predicts significantly longer survival, particularly after endocrine relapse [72,73]. While ER phosphorylated at S104, S118, S167 was detected in male BC, this did not impact on survival, contrasting the findings in female BC reported above. The association of low S104 expression that we observed with high tumour grade in male BC has been reported in female BC but was weaker [74] Alongside our deep learning work which showed that attention‐based machine‐learning trained on WSIs of female BC could accurately predict ERα status from a validation set of H&E‐stained WSIs of female BC but had poor discriminatory power when applied to male BC datasets prepared identically [75], these findings pose the intriguing possibility that ER function in BC may be sex‐specific. Indeed, sex specific features of ERα action have been identified in a genomic study. While most ERα binding sites were shared between male and female BCs in chromatin immunoprecipitation analysis, those associated with clinical outcome appeared sex specific [24]. ERα phosphorylation sites that can classify female outcome did not show predictive potential in male.

Limitations are acknowledged. The integrative bioinformatics focused solely on the genes common to all datasets. This was necessary to achieve data integration and minimise bias, but it leaves a potential gapviainformation loss through mutually exclusive genes. Like most studies on male BC, our results are limited by the number of patients available to analyse. These are small in comparison to similar work on female BC, but this is a common and well recognised obstacle when studying a rare pathology like male BC. Nevertheless, scope remains to validate these findings in larger datasets with comprehensive metadata. This was a retrospective study, which permitted accumulation of many more cases, which is an important consideration for a rare pathology. However, this meant that cases were identified and collected from a range of different sources, resulting in incomplete metadata for some cases as not all were systematically collected. However, scope remains to validate these findings in larger retrospective datasets such as the International Male Breast Cancer Program [7].

In conclusion, from the data presented herein, male BC is different biologically, in the same way that intrinsic subtypes of female BC are not identical. As BC heterogeneity is well recognised, we propose that male BC should not be separated from female BC but instead considered as a distinct and potentially unique subtype of BC. If verified mechanistically, the translational implications could potentially reshape how male BC is managed.

Author contributions statement

VSp and AHS conceived and designed the study. SC and RA‐E contributed to the study design. SC, RA‐E and VSp contributed to data analysis, interpretation and writing of the manuscript. SC, AD, MS, RA‐E and VSp contributed to immunohistochemical analysis and interpretation. SC, VSi, LO, CBM and PJvD contributed to data collection and interpretation. Collection, analyses and interpretation of mRNA gene expression datasets were performed by SC, MDM, VSi and CS. Co‐authors gave critical input, and all authors read and approved the final submitted version of the paper.

Acknowledgements

SC received an Elphinstone Scholarship awarded competitively from the University of Aberdeen and a Saltire Emerging Researcher European Exchange Award from the Scottish Funding Council. MS received an undergraduate elective bursary from the Pathological Society of Great Britain and Ireland. This work was additionally supported by project grants from NHS Grampian (19/027), Friends of ANCHOR (RS 20/2005), and the Cyril & Margaret Gates Charitable Trust (RG‐16653). Special thanks to Prof Felix Grassmann and Dr Fergal Waldron for their guidance and useful discussions about the bioinformatics work. We are grateful to staff at the NHS Grampian Biorepository, NHS Greater Glasgow and Clyde Biorepository, Northern Ireland Biobank, Wales Cancer Biobank, and the Breast Cancer Now Biobank for kindly providing tissue samples and data for the immunohistochemistry described. Thanks also to the Microscopy and Histology Core Facility in the Institute of Medical Sciences at the University of Aberdeen for their kind assistance with immunohistochemistry.

Data availability statement

The public gene expression datasets are available from NCBI‐GEOGSE31259and NCBI‐GEOGSE104730. Raw counts data in the TCGA‐BRCA dataset were downloaded individually from the Genomic Data Commons Data Portal and combined into a dataset. Other data that support the findings of this study are available from the corresponding author upon reasonable request.

Associated Data

Data Availability Statement

The public gene expression datasets are available from NCBI‐GEOGSE31259and NCBI‐GEOGSE104730. Raw counts data in the TCGA‐BRCA dataset were downloaded individually from the Genomic Data Commons Data Portal and combined into a dataset. Other data that support the findings of this study are available from the corresponding author upon reasonable request.

References

1. Rubin JB, Lagas JS, Broestl L,et al. Sex differences in cancer mechanisms. Biol Sex Differ 2020; 11: 17. doi.org/10.1186/s13293-020-00291-x
2. Chen Z, Xu L, Shi W,et al. Trends of female and male breast cancer incidence at the global, regional, and national levels, 1990–2017. Breast Cancer Res Treat 2020; 180: 481–490. doi.org/10.1007/s10549-020-05561-1
3. Humphries MP, Jordan VC, Speirs V. Obesity and male breast cancer: provocative parallels? BMC Med 2015; 13: 134. doi.org/10.1186/s12916-015-0380-x
4. Leone JP, Freedman RA, Leone J,et al. Survival in male breast cancer over the past 3 decades. J Natl Cancer Inst 2023; 115: 421–428. doi.org/10.1093/jnci/djac241
5. Murphy CE, Carder PJ, Lansdown MRJ,et al. Steroid hormone receptor expression in male breast cancer. Eur J Surg Oncol 2006; 32: 44–47. doi.org/10.1016/j.ejso.2005.09.013
6. Humphries MP, Sundara Rajan S, Honarpisheh H,et al. Characterisation of male breast cancer: a descriptive biomarker study from a large patient series. Sci Rep 2017; 7: 45293. doi.org/10.1038/srep45293
7. Cardoso F, Bartlett JMS, Slaets L,et al. Characterization of male breast cancer: results of the EORTC 10085/TBCRC/BIG/NABCG International Male Breast Cancer Program. Ann Oncol 2018; 29: 405–417. doi.org/10.1093/annonc/mdx651
8. Hassett MJ, Somerfield MR, Baker ER,et al. Management of male breast cancer: ASCO guideline. J Clin Oncol 2020; 38: 1849–1863. doi.org/10.1200/JCO.19.03120
9. Skliris GP, Rowan BG, al‐Dhaheri M,et al. Immunohistochemical validation of multiple phospho‐specific epitopes for estrogen receptor alpha (ERalpha) in tissue microarrays of ERalpha positive human breast carcinomas. Breast Cancer Res Treat 2009; 118: 443–453. doi.org/10.1007/s10549-008-0267-z
10. Kornegoor R, Verschuur‐Maes AHJ, Buerger H,et al. Molecular subtyping of male breast cancer by immunohistochemistry. Mod Pathol 2012; 25: 398–404. doi.org/10.1038/modpathol.2011.174
11. Perou CM, Sørlie T, Eisen MB,et al. Molecular portraits of human breast tumours. Nature 2000; 406: 747–752. doi.org/10.1038/35021093
12. Parker JS, Mullins M, Cheang MCU,et al. Supervised risk predictor of breast cancer based on intrinsic subtypes. J Clin Oncol 2009; 27: 1160–1167. doi.org/10.1200/JCO.2008.18.1370
13. Wallden B, Storhoff J, Nielsen T,et al. Development and verification of the PAM50‐based Prosigna breast cancer gene signature assay. BMC Med Genomics 2015; 8: 54. doi.org/10.1186/s12920-015-0129-6
14. Curtis C, Shah SP, Chin S‐F,et al. The genomic and transcriptomic architecture of 2,000 breast tumours reveals novel subgroups. Nature 2012; 486: 346–352. doi.org/10.1038/nature10983
15. Johansson I, Nilsson C, Berglund P,et al. High‐resolution genomic profiling of male breast cancer reveals differences hidden behind the similarities with female breast cancer. Breast Cancer Res Treat 2011; 129: 747–760. doi.org/10.1007/s10549-010-1262-8
16. Johansson I, Nilsson C, Berglund P,et al. Gene expression profiling of primary male breast cancers reveals two unique subgroups and identifies N‐acetyltransferase‐1 (NAT1) as a novel prognostic biomarker. Breast Cancer Res 2012; 14: R31. doi.org/10.1186/bcr3116
17. Johansson I, Ringnér M, Hedenfalk I. The landscape of candidate driver genes differs between male and female breast cancer. PLoS One 2013; 8: e78299. doi.org/10.1371/journal.pone.0078299
18. Zelli V, Silvestri V, Valentini V,et al. Transcriptome of male breast cancer matched with germline profiling reveals novel molecular subtypes with possible clinical relevance. Cancer 2021; 13: 4515–4529. doi.org/10.3390/cancers13184515
19. Lacle MM, Kornegoor R, Moelans CB,et al. Analysis of copy number changes on chromosome 16q in male breast cancer by multiplex ligation‐dependent probe amplification. Mod Pathol 2013; 26: 1461–1467. doi.org/10.1038/modpathol.2013.94
20. Kornegoor R, Moelans CB, Verschuur‐Maes AHJ,et al. Promoter hypermethylation in male breast cancer: analysis by multiplex ligation‐dependent probe amplification. Breast Cancer Res 2012; 14: R101. doi.org/10.1186/bcr3220
21. Cancer Genome Atlas Research Network , Weinstein JN, Collisson EA,et al. The cancer genome atlas pan‐cancer analysis project. Nat Genet 2013; 45: 1113–1120. doi.org/10.1038/ng.2764
22. Edgar R, Domrachev M, Lash AE. Gene expression omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Res 2002; 30: 207–210. doi.org/10.1093/nar/30.1.207
23. Parkinson H, Kapushesky M, Shojatalab M,et al. ArrayExpress – a public database of microarray experiments and gene expression profiles. Nucleic Acids Res 2007; 35: D747–D750. doi.org/10.1093/nar/gkl995
24. Severson TM, Kim Y, Joosten SEP,et al. Characterizing steroid hormone receptor chromatin binding landscapes in male and female breast cancer. Nat Commun 2018; 9: 482. doi.org/10.1038/s41467-018-02856-2
25. Tommasi S, Mangia A, Iannelli G,et al. Gene copy number variation in male breast cancer by aCGH. Cell Oncol (Dordr) 2011; 34: 467–473. doi.org/10.1007/s13402-011-0041-9
26. Callari M, Cappelletti V, de Cecco L,et al. Gene expression analysis reveals a different transcriptomic landscape in female and male breast cancer. Breast Cancer Res Treat 2011; 127: 601–610. doi.org/10.1007/s10549-010-1015-8
27. Biesma HD, Schouten PC, Lacle MM,et al. Copy number profiling by array comparative genomic hybridization identifies frequently occurring BRCA2‐like male breast cancer. Genes Chromosomes Cancer 2015; 54: 734–744. doi.org/10.1002/gcc.22284
28. Humphries MP, Sundara Rajan S, Droop A,et al. A case‐matched gender comparison transcriptomic screen identifies eIF4E and eIF5 as potential prognostic markers in male breast cancer. Clin Cancer Res 2017; 23: 2575–2583. doi.org/10.1158/1078-0432.CCR-16-1952
29. Durinck S, Spellman PT, Birney E,et al. Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt. Nat Protoc 2009; 4: 1184–1191. doi.org/10.1038/nprot.2009.97
30. Durinck S, Moreau Y, Kasprzyk A,et al. BioMart and Bioconductor: a powerful link between biological databases and microarray data analysis. Bioinformatics 2005; 21: 3439–3440. doi.org/10.1093/bioinformatics/bti525
31. Dobin A, Davis CA, Schlesinger F,et al. STAR: ultrafast universal RNA‐seq aligner. Bioinformatics 2013; 29: 15–21. doi.org/10.1093/bioinformatics/bts635
32. Anders S, Pyl PT, Huber W. HTSeq – a Python framework to work with high‐throughput sequencing data. Bioinformatics 2015; 31: 166–169. doi.org/10.1093/bioinformatics/btu638
33. Gendoo DM, Ratanasirigulchai N, Schröder MS,et al. Genefu: an R/Bioconductor package for computation of gene expression‐based signatures in breast cancer. Bioinformatics 2016; 32: 1097–1099. doi.org/10.1093/bioinformatics/btv693
34. Zhang Y, Parmigiani G, Johnson WE. ComBat‐seq: batch effect adjustment for RNA‐seq count data. NAR Genom Bioinform 2020; 2: lqaa078. doi.org/10.1093/nargab/lqaa078
35. Leek JT, Storey JD. Capturing heterogeneity in gene expression studies by surrogate variable analysis. PLoS Genet 2007; 3: 1724–1735. doi.org/10.1371/journal.pgen.0030161
36. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA‐seq data with DESeq2. Genome Biol 2014; 15: 550. doi.org/10.1186/s13059-014-0550-8
37. Lin SM, du P, Huber W,et al. Model‐based variance‐stabilizing transformation for Illumina microarray data. Nucleic Acids Res 2008; 36: e11. doi.org/10.1093/nar/gkm1075
38. Ward JH Jr. Hierarchical grouping to optimize an objective function. J Am Stat Assoc 1963; 58: 236–244.
39. Chen H, Tian T, Luo H,et al. Identification of differentially expressed genes at the single‐cell level and prognosis prediction through bulk RNA sequencing data in breast cancer. Front Genet 2022; 13: 979829. doi.org/10.3389/fgene.2022.979829
40. Alzubi MA, Turner TH, Olex AL,et al. Separation of breast cancer and organ microenvironment transcriptomes in metastases. Breast Cancer Res 2019; 21: 36. doi.org/10.1186/s13058-019-1123-2
41. Suzuki R, Shimodaira H. Pvclust: an R package for assessing the uncertainty in hierarchical clustering. Bioinformatics 2006; 22: 1540–1542. doi.org/10.1093/bioinformatics/btl117
42. Gao LL, Bien J, Witten D. Selective inference for hierarchical clustering. J Am Stat Assoc 2024; 119: 332–342. doi.org/10.1080/01621459.2022.2116331
43. Chen EY, Tan CM, Kou Y,et al. Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinform 2013; 14: 128. doi.org/10.1186/1471-2105-14-128
44. Kuleshov MV, Jones MR, Rouillard AD,et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res 2016; 44: W90–W97. doi.org/10.1093/nar/gkw377
45. Gene Ontology Consortium , Aleksander SA, Balhoff J,et al. The gene ontology knowledgebase in 2023. Genetics 2023; 224: 1–14. doi.org/10.1093/genetics/iyad031
46. Gillespie M, Jassal B, Stephan R,et al. The reactome pathway knowledgebase 2022. Nucleic Acids Res 2022; 50: D687–D692. doi.org/10.1093/nar/gkab1028
47. Huang R, Grishagin I, Wang Y,et al. The NCATS BioPlanet – an integrated platform for exploring the universe of cellular signaling pathways for toxicology, systems biology, and chemical genomics. Front Pharmacol 2019; 10: 445. doi.org/10.3389/fphar.2019.00445
48. Ogata H, Goto S, Sato K,et al. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res 1999; 27: 29–34. doi.org/10.1093/nar/27.1.29
49. Pico AR, Kelder T, van Iersel MP,et al. WikiPathways: pathway editing for the people. PLoS Biol 2008; 6: e184. doi.org/10.1371/journal.pbio.0060184
50. Liberzon A, Birger C, Thorvaldsdóttir H,et al. The molecular signatures database (MSigDB) hallmark gene set collection. Cell Syst 2015; 1: 417–425. doi.org/10.1016/j.cels.2015.12.004
51. Bankhead P, Loughrey MB, Fernández JA,et al. QuPath: open source software for digital pathology image analysis. Sci Rep 2017; 7: 16878. doi.org/10.1038/s41598-017-17204-5
52. Sørlie T, Perou CM, Tibshirani R,et al. Gene expression patterns of breast carcinomas distinguish tumor subclasses with clinical implications. Proc Natl Acad Sci U S A 2001; 98: 10869–10874. doi.org/10.1073/pnas.191367098
53. Kornegoor R, van Diest PJ, Buerger H,et al. Tracing differences between male and female breast cancer: both diseases own a different biology. Histopathology 2015; 67: 888–897. doi.org/10.1111/his.12727
54. Liu J, Sun Y, Qi P,et al. Gene expression profiling for the diagnosis of male breast cancer. BMC Cancer 2024; 24: 1584. doi.org/10.1186/s12885-024-13358-4
55. André S, Pinto AE, Silva GL,et al. Male breast cancer‐immunohistochemical patterns and clinical relevance of FASN, ATF3, and collagen IV. Breast Cancer 2021; 15: 11782234211002496. doi.org/10.1177/11782234211002496
56. Wörmann SM, Song L, Ai J,et al. Loss of P53 function activates JAK2‐STAT3 signaling to promote pancreatic tumor growth, stroma modification, and gemcitabine resistance in mice and is associated with patient survival. Gastroenterology 2016; 151: 180–193.e12. doi.org/10.1053/j.gastro.2016.03.010
57. Kenny TC, Schmidt H, Adelson K,et al. Patient‐derived interstitial fluids and predisposition to aggressive sporadic breast cancer through collagen remodeling and inactivation of p53. Clin Cancer Res 2017; 23: 5446–5459. doi.org/10.1158/1078-0432.CCR-17-0342
58. Shields MA, Ebine K, Sahai V,et al. Snail cooperates with KrasG12D to promote pancreatic fibrosis. Mol Cancer Res 2013; 11: 1078–1087. doi.org/10.1158/1541-7786.MCR-12-0637
59. Xu S, Xu H, Wang W,et al. The role of collagen in cancer: from bench to bedside. J Transl Med 2019; 17: 309. doi.org/10.1186/s12967-019-2058-1
60. Laklai H, Miroshnikova YA, Pickup MW,et al. Genotype tunes pancreatic ductal adenocarcinoma tissue tension to induce matricellular fibrosis and tumor progression. Nat Med 2016; 22: 497–505. doi.org/10.1038/nm.4082
61. Di Benedetto A, Mottolese M, Sperati F,et al. Association between AXL, hippo transducers, and survival outcomes in male breast cancer. J Cell Physiol 2017; 232: 2246–2252. doi.org/10.1002/jcp.25745
62. Di Benedetto A, Mottolese M, Sperati F,et al. The hippo transducers TAZ/YAP and their target CTGF in male breast cancer. Oncotarget 2016; 7: 43188–43198. doi.org/10.18632/oncotarget.9668
63. Christensen LG, Lautrup MD, Lyng MB,et al. Subtyping of male breast cancer by PAM50 and immunohistochemistry: a pilot study of a consecutive Danish cohort. Apmis 2020; 128: 523–530. doi.org/10.1111/apm.13068
64. Sánchez‐Muñoz A, Vicioso L, Santonja A,et al. Male breast cancer: correlation between immunohistochemical subtyping and PAM50 intrinsic subtypes, and the subsequent clinical outcomes. Mod Pathol 2018; 31: 299–306. doi.org/10.1038/modpathol.2017.129
65. Martín M, González‐Rivera M, Morales S,et al. Prospective study of the impact of the Prosigna assay on adjuvant clinical decision‐making in unselected patients with estrogen receptor positive, human epidermal growth factor receptor negative, node negative early‐stage breast cancer. Curr Med Res Opin 2015; 31: 1129–1137. doi.org/10.1185/03007995.2015.1037730
66. Bartlett JMS, Bayani J, Kornaga EN,et al. Computational approaches to support comparative analysis of multiparametric tests: modelling versus training. PLoS One 2020; 15: e0238593. doi.org/10.1371/journal.pone.0238593
67. Jagannathan V, Robinson‐Rechavi M. Meta‐analysis of estrogen response in MCF‐7 distinguishes early target genes involved in signaling and cell proliferation from later target genes involved in cell cycle and DNA repair. BMC Syst Biol 2011; 5: 138. doi.org/10.1186/1752-0509-5-138
68. Matthews J, Gustafsson JA. Estrogen receptor and aryl hydrocarbon receptor signaling pathways. Nucl Recept Signal 2006; 4: e016. doi.org/10.1621/nrs.04016
69. Thomas RS, Sarwar N, Phoenix F,et al. Phosphorylation at serines 104 and 106 by Erk1/2 MAPK is important for estrogen receptor‐alpha activity. J Mol Endocrinol 2008; 40: 173–184. doi.org/10.1677/JME-07-0165
70. Murphy L, Cherlet T, Adeyinka A,et al. Phospho‐serine‐118 estrogen receptor‐alpha detection in human breast tumors in vivo. Clin Cancer Res 2004; 10: 1354–1359. doi.org/10.1158/1078-0432.ccr-03-0112
71. Murphy LC, Niu Y, Snell L,et al. Phospho‐serine‐118 estrogen receptor‐alpha expression is associated with better disease outcome in women treated with tamoxifen. Clin Cancer Res 2004; 10: 5902–5906. doi.org/10.1158/1078-0432.CCR-04-0191
72. Yamashita H, Nishio M, Kobayashi S,et al. Phosphorylation of estrogen receptor alpha serine 167 is predictive of response to endocrine therapy and increases postrelapse survival in metastatic breast cancer. Breast Cancer Res 2005; 7: R753–R764. doi.org/10.1186/bcr1285
73. Jiang J, Sarwar N, Peston D,et al. Phosphorylation of estrogen receptor‐alpha at Ser167 is indicative of longer disease‐free and overall survival in breast cancer patients. Clin Cancer Res 2007; 13: 5769–5776. doi.org/10.1158/1078-0432.CCR-07-0822
74. Skliris GP, Nugent ZJ, Rowan BG,et al. A phosphorylation code for oestrogen receptor‐alpha predicts clinical outcome to endocrine therapy in breast cancer. Endocr Relat Cancer 2010; 17: 589–597. doi.org/10.1677/ERC-10-0030
75. Chatterji S, Niehues JM, van Treeck M,et al. Prediction models for hormone receptor status in female breast cancer do not extend to males: further evidence of sex‐based disparity in breast cancer. NPJ Breast Cancer 2023; 9: 91. doi.org/10.1038/s41523-023-00599-y

---

Republished from the open web under CC-BY. Authors: Chatterji S, Diack A, Szostok M, Morgan MD, Silvestri V, Ottini L, Moelans CB, van Diest PJ, Selli C, Sims AH, Abu-Eid R, Speirs V. Read the original.