Prognostic significance of different molecular typing methods and immune status based on RNA sequencing in HR-positive and HER2-negative early-stage breast cancer

This study was conducted to evaluate the prognostic significance of different molecular typing methods and immune status based on RNA sequencing (RNA-seq) in hormone receptor (HR)-positive and human epidermal growth factor receptor 2 (HER2)-negative (HR + /HER2-) early-stage breast cancer and develop a modified immunohistochemistry (IHC)-based surrogate for intrinsic subtype analysis. The gene expression profiles of samples from 87 HR + /HER2- early-stage breast cancer patients were evaluated using the RNA-seq of Oncotype Dx recurrence score (RS), PAM50 risk of recurrence (ROR), and immune score. Intrinsic tumor subtypes were determined using both PAM50- and IHC-based detection of estrogen receptor, progesterone receptor, Ki-67, epidermal growth factor receptor, and cytokeratins 14 and 5/6. Prognostic variables were analyzed through Cox regression analysis of disease-free survival (DFS) and distant metastasis-free survival (DMFS). Survival analysis showed that ROR better predicted recurrence and distant metastasis compared to RS (for DFS: ROR, P = 0.000; RS, P = 0.027; for DMFS, ROR, P = 0.047; RS, P = 0.621). Patients with HR + /HER2- early-stage breast cancer was classified into the luminal A, luminal B, HER2-enriched, and basal-like subtypes by PAM50. Basal-like subgroups showed the shortest DFS and DMFS. A modified IHC-based surrogate for intrinsic subtype analysis improved the concordance with PAM50 from 66.7% to 73.6%, particularly for basal-like subtype identification. High level of TILs and high expression of immune genes predicted poor prognosis. Multi-factor Cox analysis showed that IHC-based basal-like markers were the only independent factors affecting DMFS. Prognosis is better evaluated by PAM50 ROR in early-stage HR + /HER2- breast cancer and significantly differs among intrinsic subtypes. The modified IHC-based subtype can improve the basal-like subtype identification of PAM50. High immunity status and IHC-based basal-like markers are negative prognostic factors.


Background
Breast cancer is a heterogeneous disease from both the clinical and pathological perspectives [1,2]. Although early detection is critical, some cases diagnosed in early stages may persist and recur locally, cause distant metastasis, or even result in mortality. In the last decade, determination of the tumor gene expression profiles has considerably improved the understanding of the biological heterogeneity of breast cancers [2,3]. Analysis of gene expression panels, such as Oncotype DX and PAM50, can provide prognostic information beyond traditional assessments of the risk of breast cancers [4,5]. Alternatively, the Oncotype DX panel (21 genes) calculates the recurrence score (RS), and the PAM50 (50 genes) provides the risk of recurrence (ROR) score, which also identifies intrinsic subtypes. Because of the limited use of molecular detection in Chinese patients, the abilities of Oncotype and PAM50 to accurately predict prognosis in patients with hormone receptor (HR)-positive human epidermal growth factor receptor 2 (HER2)-negative early breast cancer have not been compared.
Intrinsic molecular profiling provides additional clinically relevant insights compared to current pathology-based classifications. However, the technological complexity and high operating costs have limited its use in the clinic, leading to the development of immunohistochemistry (IHC)-based surrogates for identifying intrinsic subtypes of breast cancer (luminal A, luminal B, HER2-enriched, basal-like, and normal-like [6,7]. IHC analysis of the expression of ER, progesterone receptor (PR), HER2, and the proliferation marker Ki-67 have been used to replace genotyping in clinical applications [8]. This pathology-based classification is well-established and used clinically to determine treatment modalities (endocrine therapy, chemotherapy, and/or targeted therapy) and patient inclusion in clinical trials. The discordance between IHC-and PAM50-based intrinsic subtypes was previously reported; however, the discordance in the survival data of patients with early-stage breast cancer has not been established [9,10].
The immune status, particularly that of tumor-infiltrating lymphocytes (TILs), was recently identified as a useful marker for predicting prognosis and responses to neoadjuvant chemotherapy in patients with ER-positive breast cancer [11,12]. The immune-specific gene expression patterns of TILs can be evaluated by quantifying differences in the abundance of multiple immune infiltrates in various solid tumors [13,14]. However, the prognostic role of the immune status in HR-positive/HER2-negative (HR + /HER2-) early breast cancers and its relationship with different intrinsic subtypes are unclear.
This study was conducted to compare the prognosis prediction efficiency of two multi-gene expression panels and explore the prognostic effect of the immune status on HR + /HER2-early-stage breast cancer. Additionally, we improved the IHC-based surrogate method for predicting prognosis by incorporating the immune status into the model.

Patients
We recruited 87 patients with HR + /HER2-stage 1 breast cancer who underwent curative surgery at Peking Union Medical College Hospital (Beijing, China) between July 2012 and August 2017. The median age of the patients was 48.5 years (range, 30-78 years). The median follow-up duration was 60 months. Clinical information regarding age, TNM staging, disease recurrence, and death of all patients was collected for further analysis. The TNM stage was described according to the AJCC 8 th edition, and only patients with HR + /HER2-stage 1 breast cancer was included in the study. This study is approved by the ethics committee on Human Research of Peking Union Medical College Hospital (S-K1445).

Histological analysis
Sections of formalin-fixed paraffin-embedded breast tumor tissues were re-examined by evaluating the results of hematoxylin and eosin (HE) staining to confirm the diagnosis. The expression of ER, PR, HER2, Ki-67, epidermal growth factor receptor (EGFR), cytokeratin (CK) 14, and CK5/6 was analyzed via IHC using a Ventana BenchMark XT automated slide stainer (Ventana Medical Systems, Inc., Tucson, AZ, USA) according to the manufacturer's instructions. For each antibody, positive and negative control slides were included in each staining run. IHC slides were independently assessed by two qualified pathologists. The IHC experimental conditions for all antibodies and criteria for interpreting the results have been described previously [15]. Briefly, nuclear staining for ER and PR in > 1% of the tumor cells was considered as a positive result. The HER2 status was determined by either negative (0 or 1) or equivocal (2 +) HER2 staining in IHC, and the absence of HER2 amplification was confirmed by fluorescence in situ hybridization. The expression of CK5/6, EGFR, P53, and Ki-67 was detected at the time of diagnosis. Basal markers (CK5/6, EGFR, and CK14) were considered as positive when > 1% of the tumor cells displayed plasma or membrane staining of the basal markers.
Stromal TILs were also evaluated through HE staining by two experienced pathologists following the guidelines for TIL assessment recommended by the International TILs Working Group [16]. The percentage of stromal TILs were divided into low and high levels based on the receiver operating characteristic (ROC) curve, using a cutoff value of 13.5%.

Gene expression analysis
The tumor surface area was determined by HE staining before RNA extraction. RNA was extracted from five 10-μm formalin-fixed paraffin-embedded slides using an RNAstorm extraction kit (CD201, CELLDATA, Fremont, CA, USA) following the manufacturer's instructions, and macro-dissection was performed to protect the normal breast tissue from contamination. A minimum of 100 ng of purified total RNA was used to measure the expression of 50 tumor-related genes, 17 immune-related genes, and 5 housekeeping genes by RNA-seq on the iSeq platform (Illumina, San Diego, CA, USA), and a FASTQ file was generated for each sample. To guarantee the integrity of any downstream analysis, the sequencing data were further analyzed only when samples with less than 30% of missing genes among all genes and total reads larger than 10,000 were observed. In each FASTQ file, raw counts were generated using the ShortRead package in R; after a single read, one sequence was mapped to the targeted regions of the human genome. Next, the raw counts of all samples were normalized based on the sizes of the transcripts and library. A gene expression matrix was generated by calculating the counts per million per sample using the edgeR Bioconductor software package. Finally, the data were log2-transformed and evaluated using the K-nearest neighbor method with the median centered and column standardized gene expression data.

Calculation of Relapse (RS and ROR), intrinsic subtype analysis, and immunity evaluation
The RS was calculated using the Oncotype Dx panel, and the ROR score was calculated using the PAM50 panel [4]. The proliferation-weighted ROR score was calculated as described previously [7]. Patients were assigned to one of three risk groups (low, medium, or high) according to the RS and ROR-P scores, respectively as described before [4,7].
The four intrinsic subtypes of breast cancer, namely luminal A, luminal B, HER2-enriched, and basal-like, were identified using the PAM50 predictor based on the expression profiles [7].
Moreover, an immune score was determined based on the expression levels of 17 immune genes as previously described [17]. Patients were grouped into one of two immune groups, "iweak" or "istrong, " based on the immune score values. A suitable cutoff value, 45.5, was determined using the ROC curve (Supplementary Fig  S1). An immune score ≥ 45.5 was considered as "istrong, " whereas an immune score < 45.5 was considered as "iweak".
Modified IHC-based analysis of the intrinsic subtype was performed to identify basal-like subtypes. Basal-like cases confirmed by PAM50 were investigated to determine the IHC features of ER, PR, Ki-67, and basal markers (EGFR, CK14, and CK5/6). The criteria for identifying a basal-like subtype in this assay were as follows: any one of the basal markers was positive, Ki-67 was "high" (≥ 40%), and ER was "low" (≤ 10%).

Statistical analysis
SPSS version 17.0 (SPSS. Inc., Chicago, IL, USA) was used for statistical analysis. Qualitative variables were compared using the chi-square test or Fisher's exact test. The false discovery rate was applied to multiple testing corrections.
Disease-free survival (DFS) was defined as the time from surgery to the date of the first local or distant relapse. Distant metastasis-free survival (DMFS) was defined from the date of curative surgery to the date of distant metastasis or the last follow-up. Relapsed disease and metastasis were verified by diagnostic imaging and pathology during follow-up examinations. Three patients died during the follow-up period, of whom two died from the disease. Because the number of deaths was very small, they could not be evaluated by regression analysis; therefore, we only analyzed the relationship of factors with DFS and DMFS. Survival analyses were performed using Kaplan-Meier curves. Statistical significance was set at a two-tailed p-value < 0.05.

Clinical and histological characteristics
Of the resected tumors, 85 (97.7%) were invasive breast carcinoma of no special type, two were invasive lobular carcinoma. Majority of tumors (76 of 87, 87.4%) were low or media histologic grades (grade 1 or 2). Local recurrence was observed in 26 (29.9%) patients, whereas distant metastasis was observed in 15 (17.2%) patients. The median DFS time and DMFS time were 52 months and 55 months, respectively. Three patients (3.4%) died during follow-up, two of these deaths were attributable to breast cancer (Table 1).

ROR by PAM50 and its comparison to RS by Oncotype Dx
ROR and RS were categorized as low, medium, and high; the distributions of both methods are listed in Table 2. Kaplan-Meier survival analysis showed that the ROR measured by the two methods could predict actual recurrence ( Fig. 1a; P = 0.027, Oncotype Dx; Fig. 1b; P = 0.0000, PAM50). The low, medium, and high ROR detected by PAM50 was consistent with the actual prognosis in the survival curve. The survival curve for low and medium risks overlapped with each other detected by Oncotype Dx (Fig. 1a). The ROR of PAM50 more accurately predicted DFS compared to Oncotype Dx (Fig. 1b). Neither method could effectively predict the DMFS of the low-and medium-risk groups, as the curves of the two groups overlapped in PAM50 ( Fig. 1c and d). However, the ROR significantly differentiated the high-risk group from others when predicting DMFS (P = 0.621, Oncotype Dx; P = 0.047, PAM50) ( Fig. 1c and d). The concordance of recurrence risks provided by the two methods was poor (κ = 0.194).

Distribution of intrinsic subtypes and its prognostic significance
The following distributions of the intrinsic subtypes were detected by PAM50 analysis of the 87 analyzed breast tumors: 71% luminal A, 20% luminal B, 7% basal-like, and 2% HER2-enriched (Fig. 2a). Most HR + /HER2-breast cancers were of the luminal subtype. Six cases showed the intrinsic basal-like subtype and only 2 cases exhibited the HER2 enrichment subtype (too few for survival analysis). K-M analysis showed that the luminal subtypes had a significant longer DFS time than non-luminal type ( Fig. 2b; P = 0.003) but the DMFS time did not show statistical significance ( Fig. 2c; P = 0.233). According to the ROR of PAM50, the low-risk class consisted entirely of the luminal A subtypes, whereas the high-risk class comprised four basal-like and two luminal B subtypes (Table 2). Consistent with the ROR, the luminal A subtype displayed the longest DFS, and the basal-like  [19].
The concordance of subtype identification by PAM50 and the IHC-based method was poor (kappa = 0.246).   Fig. 4a and b, respectively). Moreover, the modified IHC-based intrinsic subtypes recognized all cases of the molecular basal-like subtype with shorter DFS and DMFS times than those of non-basal-like subtype cases (for DFS, P = 0.146; for DMFS, P = 0.021; Fig. 4c and d). Further analysis revealed that all HR + /HER2-cases expressing basal-like markers had a worse prognosis, regardless of the Ki-67 index (for DFS, P = 0.29; for DMFS, P = 0.014; Fig. 4e and f ).

Impact of immune status on prognosis
The immune status was evaluated by examining the TILs and calculating the immune scores. For the relatively small number of cases in our cohort, we calculated the best cutoff value of TILs using the ROC curve as 13.5%. According to this cutoff value, there were 75 cases of low-level TILs and 12 cases of high-level TILs. High-level TILs cases showed relatively shorter DFS and DMFS, and the difference for DMFS was significant (P = 0.058 for DFS, P = 0.018 for DMFS, Fig. 5a and b). Additionally, an immune score cutoff value of 45.5 was obtained using the ROC curve to differentiate the cases into i-strong and i-weak groups. The i-strong group also showed a worse prognosis (P = 0.105 for DFS, P = 0.041 for DMFS, Fig. 5c and d).
Furthermore, the prognostic value of immune status differed according to the intrinsic molecular subtype in cases expressing basal-like markers. The three molecular basal-like cases with metastasis were istrong and had high levels of TILs. Both the non-molecular basallike cases with metastasis were weak and had low levels of TILs.

Univariate and multivariate analysis of factors affecting DFS
Univariate analysis of 14 factors (age > 50 years, tumor size > 1 cm, histology high grade, Ki-67 index > 30%, IHC basal-like marker positive, modified IHC subtype, PAM50 molecular subtypes, ROR, RS, istrong, high level TILs infiltration, P53 positivity, chemotherapy, and radiotherapy) showed that four factors significantly affected DFS: tumor size, PAM50 subtypes, RS, and ROR (Table 5, Fig. 6a). Additionally, two factors significantly affected DMFS: IHC basal-like marker and high levels of TILs (Table 6, Fig. 6b). Multivariate analysis of the relevant factors revealed that a tumor size and ROR were independent factors affecting DFS (P = 0.028 and 0.030, respectively, Table 5, Fig. 7a). IHC basal-like markers positivity were the only independent factors affecting DFMS (P = 0.029, Table 6, Fig. 7b). Immune score was another factor that had great impact on DMFS in both univariate and multivariate analyses (P = 0.051 and 0.065 respectively, Table 6, Fig. 7). Other factors showed little effect on DFS and DMFS in both univariate and multivariate analyses.  Table 3 The clinicopathological characteristics of the molecular subtypes

Discussion
Increasing evidence has shown that multi-gene expression panels used for breast cancer subtyping, such as Oncotype DX, MammaPrint, EndoPredict, Prosigna, and BCI, can provide prognostic information [18]. The clinical utility of these screening tools to guide treatment has been evaluated in several international studies [20,21]. We compared the two most common multi-gene detection methods for ROR prediction in early-stage HR + / HER2-breast cancer. Our results showed that both the ROR of PAM50 and RS of Oncotype Dx analysis could be used to predict the prognosis of patients with HR + / HER2-breast cancer. However, the concordance of recurrence risks calculated for Oncotype Dx and PAM50 was poor. The ROR of PAM50 was found to more effective predict the risk of DFS and DMFS compared to using RS. Kaplan-Meier curves showed that RS could differentiate high RS levels from low and medium RS levels to predict  [22]. ROR scores were reported to be more accurate for predicting the risk of distant metastasis compared to using RS among endocrine-treated postmenopausal women with ER-positive breast cancer [23]. Similarly, a study of 774 postmenopausal patients with ER-positive/HER2-negative breast cancer from the ATAC trial compared the prognostic values of six multigene signatures and demonstrated that ROR has a greater prognostic value for late distant metastasis compared to RS in patients with lymph node-negative breast cancer  [24]. Compared to RS, ROR segregated more patients as medium-risk and fewer patients as low-risk. The accuracy of this categorization was confirmed by the marginally separated survival curve between the low-and medium-risk groups achieved using ROR rather than using RS. Of the 6 patients in the ROR-based high-risk class, 4 had the basal-like tumor subtype and 2 had the luminal B tumor subtype. The two cases of luminal B subtype had low Ki-67 levels (10% and 15%, respectively) and high ER expression, corresponding to the mediumrisk group by the 21-gene method and low-risk group by the IHC surrogate method; therefore, neither group was administered chemotherapy. However, both patients experienced tumor relapse within 12-24 months. These patients may have benefitted from chemotherapy if they had been assigned to the high-risk class based on the ROR. In addition to ROR, PAM50-based intrinsic molecular subtyping is commonly used to guide the clinical management of breast cancer treatments. Among HR + /HER2-early-stage breast cancers, most patients with luminal A breast cancer do not benefit from chemotherapy and only require endocrine therapy. However, luminal B breast cancer may require both chemotherapy and endocrine therapy. Furthermore, basal-like breast cancer may require stronger chemotherapy, and HER2-enriched cancers may benefit from targeted anti-HER2 therapy in combination with chemotherapy or endocrine therapy [18,19,25]. In our study sample, 9.2% (8/87) of patients with HR + /HER2-early-stage breast cancer in our cohort were non-luminal, which is within the previously reported range (8%-15%) [3]. Basal-like subtypes had the shortest DFS and DMFS times.
As for the poor prognosis of non-luminal breast cancers, it is important to differentiate between the nonluminal subtypes in HR + /HER2-patients. Although several international guidelines recommend using these     all breast cancers [10]. Using the St. Gallen IHC-based method for intrinsic subtyping (2015) in our study, only two subtypes, luminal A and luminal B, were detected, whereas no basal-like or HER2-enriched subtypes were identified. Additionally, poor concordance (66.7%) of intrinsic subtype identification between PAM50 and St. Gallen IHC-based surrogate subtyping was found in a previous study [3]. Unlike HR + /HER2-luminal breast cancer, the adjuvant therapy strategy in basal-like breast cancer is often very similar to that used in triple-negative breast cancers [3]. Identification of basal-like subtypes in HR + / HER2-breast cancers is important when a gene expression profiling assay cannot be performed. We found that all cases of basal-like intrinsic subtype defined by PAM50 expressed at least one of the basal-like markers and had a high Ki-67 index (≥ 40%) and low ER expression (≤ 10%). When the three basal markers (EGFR, cytokeratin CK14, and CK5/6) were added in conjunction with low ER expression and a high Ki-67 index in the IHC-based subtype assay, we identified all six PAM50 basal-like cases, and the concordance rate between the modified IHCbased intrinsic subtyping and those detected by PAM50 was increased to 73.6%.
Accumulating evidence has shown that TILs and immune-related gene signatures can be prognostic or predictive factors in breast cancer [27]. Although the ROR of PAM50 and its intrinsic molecular subtyping can provide valuable prognostic information, it does not include immune-related genes. The prognostic effect of the immune status is well-established in triple negative and HER2-enriched breast cancers [27,28]. Moreover, TILs and immune scores are prognostic in ER-negative and highly proliferative breast cancers [17,29]. However, little evidence exists regarding the role of immunity in HR + /HER2-breast cancers. We evaluated prognosis based on TILs and immune scores for all 87 cases and found that both TILs and the immune score were factors predicting worse prognosis (particularly a shorter DMFS) for HR + /HER2-early-stage breast cancers. Our results are similar to those of other studies that showed that high TILs were an adverse prognostic factor for survival in ERpositive/HER2-negative cases [30]. Blok et al. reported that patients with excess CD8-positive TILs did not benefit from exemestane treatment in early-stage breast cancers [31]. Most patients with HR + /HER2-early-stage breast cancer was administered adjuvant endocrine therapy, which may explain why those with high levels of TILs had shorter DFS and DMFS in our study. For breast cancer prognosis, the immune gene signatures expressed by TILs must be incorporated into current multi-gene panels to improve the prognostic value and guide treatment. Although many immune-related genes have emerged as prognostic or predictive biomarkers, a study of the expression level of 130 immune-related genes revealed a high level of heterogeneity in luminal breast tumors [32]. We tested a panel of 17 immune-related genes which had a strong predictive effect in triple-negative breast cancer in patients with HR-positive early-stage breast cancer [17]. Our results showed that the immune-related gene panel had prognostic value in HR + /HER2-breast cancer, particularly for DMFS. However, because of the small size of our cohort, further validation in larger cohorts is needed. We found that the effect of the immune status on the prognosis of HR + /HER2-cases varied according to the molecular subtype.
In conclusion, we found that prognosis significantly differed among intrinsic subtypes and was better evaluated by the PAM50 ROR than Oncotype DX RS in early HR + /HER2-breast cancer. Additionally, we found that a strong immune status negatively affected the prognosis of HR + /HER2-breast cancer and that immunerelated gene signatures should be incorporated into current multi-gene tests to enhance the prognostic value of current molecular subtyping methods. When multi-gene based intrinsic subtype analysis is not available, a modified IHC-based subtyping assay for identification of basal-like subtype is recommended in HR + / HER2-patients. The impact of IHC basal-like marker expression on prognosis of HR + /HER2-breast cancer also should be pay attention to. Our study is limited mainly by the relatively small cohort sample size, which should be reinforced in a larger cohort in further studies.