The study included a cohort of female patients who were diagnosed with CRC and operated between January 1, 2008, and June 30, 2012. This investigation included 269 patients with available data on clinical information, tumor characteristics, hormonal status as well as ER, ER, KRAS, CysLT₁R, CysLT₂R, COX-2, 15-PGDH, β-catenin, Mucin-2 and PGD2 synthase expression in CRC tissue. The study population is briefly described in the Supplementary Materials. Details about the study design, patient follow-up and data collection are provided elsewhere (7).

Tumor samples were retrieved and incorporated into tissue microarray (TMA) blocks based on the protocol described earlier (7). The tissues were stained with specific antibodies for the expression of ERα ERβ and other proteins of interest (Supplementary Material). Two independent investigators (GT and RE), blinded to the patient and tumor characteristics, evaluated the staining immunoreactivity using the immunoreactive score (IRS) with a range 0–9, which was calculated as a multiplication of staining intensity (0 = negative, 1 = weak, 2 = moderate and 3 = strong) with percentage of positive stained cells (1 = <10%, 2 = 11–50% and 3 = >50%) (7). The staining intensity was determined based on the criteria of Konstantinopoulos et al. (4), which are described in the Supplementary Materials. For ERα and ERβ expression, only the nuclear staining intensity was taken into consideration, based on which they were also scored as categorical variables, respectively low/high and negative/positive expression (Figure 2A). Briefly, negative and weak ERβ staining were grouped as low expression and moderate and strong ERβ staining as high expression (7). Because ERα is very little expressed in the normal colonic mucosa (1, 2), we defined its expression as positive if more than 10% of the nuclei were stained, regardless the staining intensity. All the other tumor samples that had <10% of the nuclei stained, regardless the staining intensity, were considered to have negative ERα expression. Each tumor sample was in duplicate. Cores with loss of tissue or with only stromal tissue were excluded from the analysis.

Normalized RNA sequencing data in transcripts per million (TPM), reverse phase protein array (RPPA) data, and the associated clinical information of the colon adenocarcinoma (COAD) samples were downloaded from the TCGA dataset (https://portal.gdc.cancer.gov/; https://tcpaportal.org/tcpa/; ≤ June 20, 2020). Out of 361 patients, 12 patients missing pathological information, 16 patients with a follow-up period of ≤30 days, and 52 patients with metastasis (stage IV) were eliminated. Thus, 282 patients with clinical information were included in the study. Normalized gene expression and protein expression data from the TCGA-COAD dataset were log2-transformed for further analysis.

To identify independent prognostic parameters and to validate the independent prognostic value of ERα and ERβ, univariate and multivariate Cox regression analyses were performed in the TCGA-COAD dataset on the ERα and ERβ gene and protein signature and clinicopathological parameters. Parameters with P < 0.05 in the univariate analysis were further included in the multivariate Cox regression analysis. The TCGA samples were divided into high- and low-risk groups according to the optimal cutoffs determined by the Youden Index association criteria and analyzed using Circos visualization package (24).

The variables were compared between the group of interest using Pearson's χ² test or Fisher's exact test for categorical variables and the Mann-Whitney U test or t-test for continuous variables. Survival curves, generated via the Kaplan-Meier method, were compared between the groups using the log-rank test. Univariate and multivariate Cox proportional hazards regression models were applied, and hazard ratios (HRs) together with 95% confidence intervals (CIs) were calculated to determine the risk of death or cancer recurrence. Receiver operating characteristic (ROC) curves were used to calculate the area under the curve (AUC) to determine the predictive ability of the final model with combined ERβ + ERα expression compared to models with only one ER expression or the basic model. Binary logistic regression model was used to determine the odds ratios (ORs) of having a metastatic event for each unit increase in ERα and ERβ intensity. The estimates with their corresponding 95% CIs were used to build forest plots by the ggplot2 package in R. Statistical analyses were performed using SPSS version 23.0 (SPSS, IBM, Armonk, NY, USA) and GraphPad Prism version 8.0a (GraphPad Software, Inc., San Diego, CA, USA). A two-sided P < 0.05 was considered statistically significant.

We had 306 primary CRC samples available for the evaluation of ERα and ERβ expression. Fourteen patients, who were previously operated and treated for breast cancer, were excluded from the study due to the risk of ERα alterations from the anti-estrogen therapies (Figure 1). We successfully evaluated ERβ in 300 CRC patients and ERα in 270 CRC patients. Based on the staining intensity assessed with IHC, ERβ expression was categorized as low and high, while ERα expression was categorized as negative and positive (Figure 2A). We next compared the expression of these receptors between normal and matched cancer tissues and found that compared to ERα expression levels, ERβ expression levels were higher in both normal and cancer tissues (Figure 2B). However, compared to normal tissues, a downregulation of ERβ and an upregulation of ERα were observed in the matched CRC tissues (Figure 2B, see violin bar graph). Since we previously reported that high ERβ expression correlated with better prognosis in CRC (7), we investigated the distribution of ERα expression in patients with low and high ERβ expression. We grouped the patients into four categories based on ERα and ERβ expression (Figure 2C). We found that 79% of patients with high ERβ expression had also negative ERα expression compared with 63% in the low ERβ group (Figure 2D). Likewise, the percentage of patients with positive ERα expression was higher in the low ERβ expression group (37%) than in the high ERβ expression group (21%) (Figure 2D). For representative IHC images of matched pairs of patients for both ERα and ERβ expression, see Supplementary Figure 1A.

Next, we used ESR1 (ERα) and ESR2 (ERβ) mRNA levels from the TCGA-COAD database to investigate the differential expression of ERα and ERβ in CRC patients with TNM stage I disease and TNM stage IIIc+IV disease. Compared to those with stage I disease, a smaller percentage of patients with stage IIIc+IV disease had upregulated ESR2 mRNA levels (Figure 2E). Additionally, ESR2 levels were lower in patients with stage IIIc+IV disease than in those with stage I disease (Figure 2E). Furthermore, ESR1 mRNA levels were obviously higher in patients with stage IIIc+IV disease than in those with stage I disease (Figure 2E).

Because the role of ERα expression in CRC is very little studied and all our results are based on antibody staining, we tested the specificity of the antibodies we used, in order to validate the antibodies. First, we stained the normal breast tissue, which is known to abundantly express ERα (positive control), and normal kidney, prostate, and skin tissues, which are known to lack ERα expression (negative controls, Supplementary Figure 1B) (25–27). Next, the same tissues were also stained with another anti-ERα antibody, D12 (Supplementary Figure 1C), which is widely used for the detection of ERα expression (28–30). We randomly stained 59 patients from the Female cohort with the D12 antibody. As shown in Supplementary Figure 1D the distribution of the IRS for nuclear ERα expression for each patient (n = 59) was the same for both antibodies. Likewise, when the patients were grouped as positive and negative nuclear ERα expression, no significant difference was observed between the two antibodies (P = 0.11, Supplementary Figure 1E). Out of 59 patients randomly stained with D12 antibody, 13 patients (22%) were positive for ERα expression, while 19 patients (32%) were detected as positive using the cocktail antibody (Supplementary Figure 1E). This could be explained by the fact that the cocktail antibody 1D5 + 6F11 was created by mixing two monoclonal antibodies that detect two different epitopes (31, 32). Representative IHC images of matched-pair CRC tissues for both antibodies are shown in the Supplementary Figure 1F.

Out of 252 patients with successful staining for the KRAS mutation, only 31 (12.3%) had positive staining (Figure 2F). Patients with a KRAS mutation had a significantly higher intensity of ERα expression (P < 0.05) and a tendency to have lower ERβ expression (P = 0.06) than patients with wild-type (WT) KRAS (Figure 2F). Additionally, we observed that 19% of patients with positive ERα expression had KRAS mutations, while 9% of patients with negative ERα expression had KRAS mutations (Figure 2G). An opposite tendency was observed when looking at the distribution of KRAS mutations in patients with low and high ERβ expression. While 15% of patients with low ERβ expression had KRAS mutations, only 7% of patients with high ERβ expression had KRAS mutations (Figure 1G). However, no statistical significance was reached. To further validate these findings, we used mRNA data from the TCGA-COAD public database and found a strong and significant positive correlation between the mRNA levels of ESR1 (ERα) and KRAS mutations, while no correlation was found with ESR2 mRNA levels (ERβ) (Figure 2H).

Previously we reported that high nuclear ERβ expression is independently associated with better OS and DFS in female CRC patients (7). Herein, we report that CRC patients with negative nuclear ERα expression have 19% lower risk for 5-years overall mortality (HR = 0.81; 95% CI, 0.68-0.94; P = 0.042, Figure 3A). Likewise, in the TCGA-COAD cohort, low ERα protein expression (HR = 0.73; 95% CI, 0.62-0.92; P = 0.035, Figure 3B) and high ERβ protein expression (HR = 0.78; 95% CI, 0.68-0.89; P = 0.001, Figure 3C) are associated with better prognosis of CRC patients. Additionally, we investigated the predicting ability of ERα and ERβ expression in our female patient's cohort calculating the ROC curves. We found that ERα expression predicts the 5-years OS with higher specificity (AUC = 0.720, Sensitivity = 65.22 and Specificity = 79.37, Figure 3D), while ERβ expression with higher sensitivity (AUC = 0.674, Sensitivity = 71.05 and Specificity = 49.42, Figure 3E). When we combined the ERα and ERβ expression, the predicting ability for 5-years OS in CRC patients was significantly improved with higher sensitivity and higher specificity (AUC = 0.842, Sensitivity = 71.53 and Specificity = 82.90, Figure 3F). Next, we looked at the risk score profile with TNM-stage and 5-years OS event by combining the ERα and ERβ expression in four groups as described above (Figure 2C). As shown in Figure 3G, the subgroups with positive ERα expression had the highest risk score profile, while the patients with negative ERα expression had the lowest risk score profile, despite the ERβ expression levels.

Next, we investigated the combined role of ERα and ERβ expression in CRC OS and DFS (Figure 4). The Cox regression analysis showed that patients with combined high ERβ + negative ERα expression were independently associated with better OS and had a 77% reduction in overall mortality (Figures 4A,B, Supplementary Table 1), as well as better DFS with a 90% reduction in cancer recurrence (Figures 4C,D, Supplementary Table 1) after adjustment for age, TNM stage and tumor vascular invasion, compared to patients with combined low ERβ + positive ERα expression, which were taken as the reference group. This finding was consistent even for the subgroups of patients with stage I-III cancer (Figure 4E), patients with colon cancer (Supplementary Figures 2A,B) and patients who did not receive adjuvant treatment (Figure 4F and Supplementary Figure 2C). In the second group of patients with low ERβ expression, even though the expression of ERα remained negative, the risk was increased by 14% for overall mortality and 33% for cancer recurrence compared to patients with combined high ERβ negative ERα expression (Supplementary Table 1). In addition, in the third group of patients with positive ERα expression, even though the expression of ERβ was high, the increase in the risks of overall mortality and cancer recurrence was much lower than that in the first group with combined high ERβ + negative ERα expression (3 and 22% lower, respectively: Supplementary Table 1, multivariate analysis). It is difficult to draw any conclusions about the subgroup of patients with rectal cancer due to the very small number of patients in each category, especially the category with combined high ERβ + positive ERα expression that has only one patient, n = 1 (Supplementary Table 1, Supplementary Figures 2D,E). These results clearly show that CRC patients with combined high ERβ + negative ERα expression have the best prognosis and that the subgroup with combined low ERβ + positive ERα expression has the worst prognosis.

To further investigate the role of the combined ERα and ERβ expressions in predicting CRC prognosis, we evaluated the ROC curves for the basic model (adjusted for age, TNM stage and tumor vascular invasion), the model extended with only ERβ expression, the model extended with only ERα expression, and the model that included the combined ERβ + ERα expressions. As shown in Figures 4G,H, the AUC was significantly higher for the model with the combined ERβ + ERα expressions than for all the other models for both OS and DFS. However, the predictive ability of the combined ERβ + ERα extended model was higher for DFS (AUC = 0.812, Figure 4H') than for OS (AUC = 0.801, Figure 4G'). The same results were obtained using the TCGA-COAD external cohort, where the combined expression of ERs had the best predictive ability for DFS compared with the other models (Figures 4I,I'). These results clearly show that the combined expression of ERα and ERβ plays an important role in predicting the prognosis of CRC patients.

We aimed to evaluate the distribution of clinical parameters and tumor characteristics between patients with combined high ERβ + negative ERα expression, considered to be the best prognostic group, and those with combined low ERβ + positive ERα expression, considered to be the worst prognostic group. As shown in Table 1, patients with combined high ERβ + negative ERα expression had a significantly lower number of overall deaths and cancer recurrence events, smaller tumor extent, fewer tumor metastases in the regional lymph nodes and distant organs, predominantly stage I and II disease, and were less likely to receive adjuvant treatment after the operation. Additionally, tumors with combined high ERβ + negative ERα expression had a higher frequency of the mucinous type of COAD and a never smoking status (Table 1).

We explored the hormonal characteristics of CRC female patients in relation to the combined ERα and ERβ expression. We found that female patients with combined high ERβ + negative ERα expression had a lower number of pregnancies (mean ± standard error of the mean, 1.8 ± 0.13, P = 0.04; Figure 5A) and shorter breastfeeding times (calculated as the total breastfeeding months for all the children a woman had; 8.2 ± 0.95, P = 0.08; Figure 5B) than female patients with combined low ERβ + positive ERα expression (2.2 ± 0.14 and 10.8 ± 1.2, respectively). No significant differences were observed between the two groups regarding the age of menopause and age of menarche (Figures 5C,D). Next, we examined how the use of hormonal contraception (HC) differed between the two groups. We found that most of the female patients with combined high ERβ + negative ERα expression never used HC compared with women with combined low ERβ + positive ERα expression (63% vs. 37%, P = 0.02, Figure 5E). When we looked at the type of HC, we found that 61% of female patients with combined high ERβ + negative ERα expression had never used combined (estrogen and progesterone) HC and 48% of them had used combined HC. In the subgroup of women with combined low ERβ + positive ERα expression 39% had never used combined HC and 52% had used combined HC (P = 0.07, Figure 5F). However, no difference was observed between the two groups regarding the use of progesterone HC (Figure 5G). We also looked at the use of hormone replacement therapy (HRT) and found that most of the female patients with combined high ERβ + negative ERα expression had used HRT for more than 5 years, while very few female patients with combined low ERβ + positive ERα expression had used HRT for a long time (71 and 29%, respectively, P = 0.02, Figure 5H). All the female patients who had used combined (estrogen and progesterone) HRT had combined high ERβ + negative ERα expression (P < 0.0001; Figure 5I). No significant results were found regarding the use of estrogen HRT (Figure 5J).

To further explore the prognostic role of combined ERα and ERβ expression in CRC patients, we correlated the patient with combined high ERβ + negative ERα expression or combined low ERβ + positive ERα expression with proteins important in CRC development and progression (Figure 6A). We noticed that patients with combined low ERβ + positive ERα expression had lower IRSs for CysLT₁R (P < 0.01), COX-2 (P < 0.001) and nuclear β-catenin (P < 0.001), which are connected to enhanced cell proliferation and poor patient outcome (18, 19, 23), compared to patients with combined high ERβ + negative ERα expression (Figure 6A, Supplementary Figure 3 for IHC images). On the other hand, patients with combined high ERβ + negative ERα expression had higher IRSs for CysLT₂R (P < 0.001), membrane β-catenin (P < 0.001), 15-PGDH (P < 0.01) and PGD2 synthase (P < 0.001), which are associated with a better outcome in CRC (20, 23, 33, 34) (Figure 6A, Supplementary Figure 3). Since we observed a higher frequency of mucinous adenocarcinomas in the group of patients with combined high ERβ + negative ERα expression, we investigated the association with Mucin-2 expression known to be reduced in CRC tissues compared to the normal mucosa (35, 36). We found that patients with combined high ERβ + negative ERα expression had significantly higher IRSs for Mucin-2 expression levels (P < 0.05) than patients with combined low ERβ + positive ERα expression (Figure 6A, Supplementary Figure 3). In the TCGA-COAD cohort, the same correlations were observed between the combined protein expression of ERs and CysLT₁R, COX-2, CysLT₂R and PGD2 synthase, whereas no correlation was found for combined ERs expression with 15-PGDH and Mucin-2 expression levels (Figure 6B).

We investigated the risk of having a metastatic event for each unit increase in the ERβ and ERα staining intensity, evaluated by IHC. We found that for each unit increase in the ERβ intensity, the risk of having a metastatic event were significantly and independently decreased by 60% after adjustment for age, TNM stage and tumor vascular invasion (OR = 0.40; 95% CI: 0.19–0.82; P = 0.012; Figure 7A). In addition, for each unit increase in the ERα intensity, the risk of having a metastatic event increased almost 2.5-fold (OR = 2.47; 95% CI: 1.15–5.32; P = 0.021; Figures 7A,B). The ERα intensity was strongly associated with liver metastasis, where for each unit increase in the ERα intensity, the risk of liver metastasis independently increased almost 4-fold (OR = 3.72; 95% CI: 1.36–10.17; P = 0.01; Figures 7A,B). However, no role of ERβ was found in lung metastasis and the promoting effect of increased ERα staining intensity (OR = 3.48; 95% CI: 1.38–8.77; P = 0.008) disappeared after adjustment for other confounding factors (OR = 3.05; 95% CI: 0.99–9.42; P = 0.052; Figure 7A). Importantly, each unit increase in the ERβ intensity significantly and independently decreased the risk of local recurrence and abdominal metastasis by 79% (OR = 0.21; 95% CI: 0.06–0.67; P = 0.009; Figures 7A,B). These results were summarized graphically using the forest plots, where the increased risk is shown in red, and the decreased risk is shown in blue (Figure 7B).