Prognosis and Genomic Landscape of Liver Metastasis in Patients With Breast Cancer

Objective The prognosis of breast cancer liver metastasis (BCLM) is poor, and its molecular mechanism is unclear. We aimed to determine the factors that affect the prognosis of patients with BCLM and investigate the genomic landscape of liver metastasis (LM). Methods We described the prognosis of patients with BCLM and focused on prognosis prediction for these patients based on clinicopathological factors. Nomogram models were constructed for progression-free survival (PFS) and overall survival (OS) by using a cohort of 231 patients with BCLM who underwent treatment at Shandong Cancer Hospital and Institute (SCHI). We explored the molecular mechanism of LM and constructed driver genes, mutation signatures by using a targeted sequencing dataset of 217 samples of LM and 479 unpaired samples of primary breast cancer (pBC) from Memorial Sloan Kettering Cancer Center (MSKCC). Results The median follow-up time for 231 patients with BCLM in the SCHI cohort was 46 months. The cumulative incidence of LM at 1, 2, and 5 years was 17.5%, 45.0%, and 86.8%, respectively. The median PFS and OS were 7 months (95% CI, 6–8) and 22 months (95% CI, 19–25), respectively. The independent factors that increased the progression risk of patients with LM were Karnofsky performance status (KPS) ≤ 80, TNBC subtype, grade III, increasing trend of CA153, and disease-free interval (DFS) ≤ 1 year. Simultaneously, the independent factors that increased the mortality risk of patients with LM were Ki-67 ≥ 30%, grade III, increasing trend of CA153, pain with initial LM, diabetes, and DFI ≤ 1 year. In the MSKCC dataset, the LM driver genes were ESR1, AKT1, ERBB2, and FGFR4, and LM matched three prominent mutation signatures: APOBEC cytidine deaminase, ultraviolet exposure, and defective DNA mismatch repair. Conclusion This study systematically describes the survival prognosis and characteristics of LM from the clinicopathological factors to the genetic level. These results not only enable clinicians to assess the risk of disease progression in patients with BCLM to optimize treatment options, but also help us better understand the underlying mechanisms of tumor metastasis and evolution and provide new therapeutic targets with potential benefits for drug-resistant patients.


INTRODUCTION
Breast cancer (BC) is a malignant tumor with the highest incidence in women, and the trend of rejuvenation is significant (1). Metastatic breast cancer (mBC) is constantly diagnosed, and the position of metastatic organs strongly correlates with survival time, despite progress has been made in the treatment and prognosis of early BC (2)(3)(4). The prognosis of liver recurrence is the second-worst outcome after brain metastasis (2,4). About half of the patients with mBC eventually develop liver metastasis (LM) (5), and this probability is increasing every year (6).
The American Joint Committee on Cancer Staging is most commonly used to predict the prognosis of cancer patients. It includes TNM staging, pathological grade, and tumor expression status of biological indicators such as estrogen receptor (ER), progesterone receptor (PR), and human epidermal growth factor receptor 2 (HER-2). Additionally, some investigations have pointed out that biological and pathological parameters could be used for predicting the recurrence or prognosis of patients with BC (7,8). However, an authoritative prediction model for the prognosis of patients with breast cancer liver metastasis (BCLM) has not yet been developed.
In the early years, the mechanisms of BCLM were not explored extensively, and mainly the inherent structure and microenvironment of the liver were focused on. For example, the high expression of Claudin-2 can enhance the adhesion to extracellular matrix proteins, thereby increasing the potential of BC cells to transfer to the liver (9). Another example is that the interaction between chemokine receptors on tumor cell membranes and chemokines in the microenvironment, such as CCL2 (10,11) and CXCR4-CXCL12/SDF-1 (12)(13)(14), plays important roles in LM. Moreover, with the rapid development of sequencing technology in the recent years, a large number of early BC gene maps have been reported, along with several reports on the overall mutation characteristics of mBC (15)(16)(17)(18). However, the genetic landscape related to LM alone has not yet been characterized.
In this study, we aimed to focus on prognosis prediction and construct nomogram models for progression-free survival (PFS) and overall survival (OS) by using clinicopathological characteristics, especially the innovative application of the dynamic changes of hematological indicators, in a cohort of patients with BCLM who underwent treatment at Shandong Cancer Hospital and Institute (SCHI). We investigated the molecular mechanism of LM and constructed driver genes, mutation signatures (MSs) by using a targeted sequencing dataset from the Memorial Sloan Kettering Cancer Center (MSKCC), which was verified with the dataset of the Predictive Oncology team of the University of Aix-Marseille (POTUAM) (Figure 1).

SCHI Cohort
Patients who met the following criteria were eligible for this study: (1) diagnosed with BCLM and treated continuously with SCHI ( Figure 1) from December 2008 to December 2018, (2) complete clinical pathology records, (3) age ≥ 18 and ≤ 75, and (4) Karnofsky performance status (KPS) ≥ 60. Patients who met any of the following criteria were excluded from the study: (1) bilateral BC, (2) primary and/or metastatic liver cancer, (3) other invasive malignant diseases within five years, and (4) medical records deemed unqualified according to the investigator's opinion. All the medical records of patients with BCLM diagnosed by pathology or imaging were collected retrospectively. This study was approved by the Institutional Review Board, and the personal information of all patients and attending doctors was also deleted from the data set. Clinical features, pathological features, imaging examinations, treatment methods, and survival information were collected by two independent researchers according to a standardized process, and disagreements were resolved through discussions with a third expert.
The date of LM diagnosis was based on the date of the imaging examination or the date of biopsy pathology report. Disease-free interval (DFI) refers to the length of time from BC diagnosis to relapse or metastasis. PFS refers to the length of time from the diagnosis of BCLM to disease progression. OS refers to the length of time from BCLM diagnosis to death from any cause or the last follow-up. All surviving patients at the time of analysis were censored at the date of their last follow-up.

Significant Mutant Genes and Mutational Signatures
Significantly mutant genes were identified with Maftools (20) across the entire cohort of pBCs and LMs. Oncodrive based on the algorithm oncodriveCLUST (21) was used to identify cancer genes (driver) from a given mutant allele fractions. We then used the NMF (22) algorithm to decompose MSs based on the set of known signatures from the COSMIC database (23), and calculated the cosine similarity to identify the best match.

Statistical Analysis
Chi-square test or Fisher's exact test was used to evaluate differences in count data. Comparisons of tumor mutational burden were performed using Mann-Whitney U test. We used the Kaplan-Meier method for survival analysis and evaluated the difference between the Kaplan-Meier curves by applying the logrank test. p < 0.05 was considered statistically significant.
The Cox proportional hazards model was used to analyze the univariate and multivariate factors associated with survival. Finally, we used the RMS r package to draw the nomograms. All statistical analyses were performed using SPSS software for mac v26 and R v4.0.0.

Clinicopathological Characteristics of the SCHI Cohort
Between December 2008 and December 2018, 410 patients with BCLM were admitted to our hospital. A total of 231 patients were eligible. Of the 231 patients, 213 had first recurrent LM and 18 had subsequent recurrent LM. Among the patients with first recurrent LM, 69 had only LM, and 144 had metastasis to other organs. (Figure 1). The median follow-up time in this study was 46 months (range, 12-118 months). The clinicopathological characteristics of the patients at baseline are shown in Table 1.
Notably, we observed highly accumulated ESR1 mutations in LM, and the ESR1 mutation rate of LM (20%) was significantly higher than that of pBC (3%) (Fisher's exact test, p = 6.20 × 10 −12 ) ( Figure 5D). TP53, as the gene with the highest mutation rate, mutated exclusively with ESR1 in LM (pair-wise Fisher's exact test, p = 6.11 × 10 −5 ) ( Figures 5B, 6A). ERBB2, which also had a high mutation frequency in LM, was another gene that was mutually exclusive with the ESR1 mutation ( Figures 5B, 4A). Moreover, we found that almost all GATA3 mutations were accompanied by ESR1 mutations (Figure 6A), and the incidence of GATA3 mutations in LM was slightly higher than that in pBC. ESR1 mutations in LM patients were not only more frequent but also more concentrated in mutation sites, in contrast to the low frequency and scattered mutation patterns in patients with pBC. Thirty-five LM samples in total carried mutations were resistant to ESR1 aromatase inhibitor (AI) (D538G:19, Y537S/N:16), compared with only three pBC samples (D538G:1, Y537S/N:2) (Fisher's exact test, p = 8.37 × 10 −16 ) ( Figure 6B). To discover the biological pathways that play a key role in LM, we enriched the mutation matrix with known oncogenic signaling pathways in TCGA cohorts (24). Four oncogenes (FGFR4, ERBB2, ROS1, and IGF1R) and one suppressor gene (NF1) in the RTK-RAS pathway, which ranked first, were mutated more frequently in LM (Figures 6C, D).

Three MSs in Patients With BCLM
For the purpose of clarifying the etiological mechanism of the different mutation rates between LM and pBC in the MSKCC dataset, and explaining the potential mechanism of BC metastasis to the liver, non-negative matrix factorization (NMF) created by Nik-Zainal S et al. (22) was used to extract MSs from 96 subtypes of three-base context of mutations. A total of three prominent signatures were matched ( Figure 7A): MS1, best matched to COSMIC signature 13 (cosine-similarity: 0.822), is characterized by C > T mutations at the TpC dinucleotide and an APOBEC Cytidine Deaminase (C > G) phenotype; MS2, best matched to COSMIC signature 7 (cosine-similarity: 0.605), is characterized mainly by C > T and T > C mutations with other types of base substitutions contributing less and intricate patterns formed due to exposure to ultraviolet; and MS3, best matched to COSMIC signature 6 (cosine-similarity: 0.724), is mainly characterized by C > T mutations and a map caused by defective DNA mismatch repair ( Figure 7B).

DISCUSSION
We described the prognosis of patients with BCLM, established a prognostic prediction model based on clinical pathological factors, and characterized the genomic landscape of patients with LM with external dataset in this study. We innovatively incorporated the dynamic changes of traditional tumor markers, such as CEA, CA153, and blood indicators, such as granulocytelymphocyte ratio, in the early stage of treatment into the prognostic analysis of LM. It is encouraging to note that the dynamic changes of CA153 greatly improved the prediction performance of the model. The median follow-up time for the 231 patients with BCLM in the SCHI cohort was 46 months. In this cohort, the cumulative incidence of LM within 5 years was as high as 86.8%. The median PFS and OS were 7 months (95% CI, 6-8) and 22 months (95% CI, [19][20][21][22][23][24][25], respectively (Figures 2A, B). PFS at 1 year was 25.1%, and the OS at 1 and 2 years was 77.1% and 45.8%, respectively. We found that the OS of patients with first recurrent only LM (n = 69) was significantly longer than that of patients with first recurrent LM with other organs (n = 144) and subsequent recurrent LM (n = 18) (p = 0.0036). This may be because patients with LM as the only first site of metastasis had better KPS (chi-square test, p = 0.005). Conversely, patients who are accompanied by metastasis to other organs tend to have heavier tumor burden and higher tumor heterogeneity, and are more susceptible to drug resistance. In addition, our results support the view that the natural process of LM is also strongly influenced by the biology of BC subtypes, which is consistent with previous studies (6,25,26). Patients with TNBC subtype were more prone to LM in the early postoperative period ( Figure 2G) and had shorter PFS ( Figure 2E) and OS ( Figure  2F). Considering that the molecular subtypes of 64 LMs obtained a 32.8% (21/64) mutation rate compared with the paired pBCs (Figure 4), biopsies of metastasis are useful for the reassessment of the metastatic sites to define a more effective treatment strategy for patients with BCLM (27). Thus, the ER, PR, and HER-2 statuses need to be reassessed by biopsy when LM occurs. The indicators of KPS, Ki-67, and grade were recorded in previous reports of mBC. KPS ≤ 70 (28), Ki-67 ≥ 20% (29,30), and grade III (31) were independent risk factors for the prognosis of mBC. These results are generally consistent with those of our study. Moreover, in a study by Nishimura et al. (29), DFI is inversely correlated with the Ki-67 values. We also observed that the shorter the DFI, the worse the prognosis of patients with LM. Finally, our data showed, for the first time, that patients with diabetes and initial pain in the liver have a shorter OS.
In addition to the above factors that have an impact on the prognosis of LM, we have included the dynamic changes of CA153, CEA, and granulocyte-lymphocyte ratio after two cycles of LM treatment as prognostic factors into consideration for the first time. The increasing trend of CA153 showed a strong correlation with shorter PFS (HR: 2.79, 95% CI: 1. 86-4.16; p < 0.001) and OS (HR: 2.32, 95% CI: 1. 60-3.37; p < 0.001). More importantly, the dynamic changes of CA153 could reflect the curative effect in a timelier and more accurate manner, and help decision makers adjust the treatment plan according to the curative effect. To summarize, we established a prognostic model for BCLM based on these clinical and pathological factors (Figure 3). Compared with the prognostic nomogram model without CA153_trend (Supplementary Figure 2), the prognostic nomogram model including the dynamic changes of CA153 had greatly improved predictive ability (C-index: 0.743  In the SCHI cohort, the critical cause of the 32.8% molecular subtype mutation rate of LMs relative to pBCs is changes at the genome level. The genomic landscape of early BC has been reported many times. In addition, there is evidence that genomic changes are obtained in the process of cancer metastasis and progression and the genomic pattern of early cancers cannot represent lethal cancers (32)(33)(34)(35)(36). Recently, there have been subsequent reports on the genomic landscape of mBC (16)(17)(18). However, the molecular mechanism for LM remains an area that has not yet been fully developed (37). In view of the above, we decoded the genomic landscape of LM by using the MSKCC dataset (19). Among the 14 SMGs in the LM samples, we identified four driver genes: ESR1 (20%), AKT1 (8%), ERBB2 (7%), and FGFR4 (4%). Compared with the 21 potential driver genes (TP53, ESR1, CDH1, MAP3K1, GATA3, CBFB, ARID1A,  (16), most of the four driver genes (3/4, ESR1, AKT1, and ERBB2) of LM were more consistent, verifying that ESR1, ERBB2, and AKT1 can drive LM. The other 18 inconsistent driver genes of mBC also indicate that the driver genes of LM are not exactly the same as those of other metastases, that is to say these inconsistent driver genes of mBC may play a role in metastasis in other organs, but not in liver. To further verify our results, we quoted the POTUAM dataset. The mutation frequencies of the four driver genes in the MSKCC dataset were highly consistent with those in the POTUAM dataset (Supplementary Figure 1). Furthermore, in the POTUAM dataset, except for FGFR4, which had a low mutation rate, the mutation frequencies of ESR1, AKT1, and ERBB2 in the LM were higher than those in other metastases. Moreover, LM had a significantly different mutation spectrum from other metastases (Supplementary Figure 3). Our study characterized the driver genes of LM isolated from other metastases, making the molecular targets of LM clearer and more focused. Interestingly, ESR1 and TP53 were the most mutually exclusive mutant gene pair in LMs (pairwise Fisher's exact test, p = 6.11 × 10 −5 ). As previously reported (15,38), ESR1 mutations were enriched in patients with LMs (Fisher's exact test, P < 0.001). In contrast to the low frequency and scattered mutation patterns in the pBC cohort, the ESR1 mutations were mainly concentrated in D538G and Y537S/N. This is consistent with a previous study published in the journal JAMA Oncology (39). Another notable gene is NF1. In a recent study, Bertucci et al. (15) performed whole-exome sequencing of 617 tumor samples from patients with metastatic BC and revealed that the mutation frequency of NF1 is negatively correlated with the prognosis of patients with HR +/HER-2-mBC. In the current study, we found that as a tumor suppressor gene, NF1 had the most frequent mutation in the RTK-RAS pathway in patients with LM. To better understand which mutation processes facilitated the progression of LM, we further evaluated the distribution of MSs. LMs matched three prominent signatures: S13 (APOBEC cytidine deaminase), S7 (ultraviolet exposure), and S6 (defective DNA mismatch repair). According to previous reports, APOBEC activation can mediate secondary resistance to endocrine therapy (40). Evaluating BCLM using a next-generation sequencing approach helps us better understand the underlying mechanisms of tumor metastasis and evolution and provide new therapeutic targets with potential benefits for drugresistant patients. Moreover, owing to the improvement of the current circulating tumor DNA detection technology, if the mutant LM driver genes can be detected in the circulating blood, it may be a powerful tool for early warning of LM. Finally, by comparing the TMB between LMs and pBCs (2.77 vs. 1.98, U test, p < 0.001), we also confirmed that metastatic foci are more genetically complex than primary foci, as in previous studies (32,33,35).
Our study has several limitations. First, this research was retrospective, so there may be some unavoidable biases. Second, the construction of the prognosis prediction model was based on a SCHI single-center cohort; thus a prospective multi-center verification is required before subsequent promotion. Third, we used an external sequencing dataset to characterize the genomic landscape of BCLM, and further exploration requires the support of more cell and animal experiments.
In summary, BCLM is a complex process that involves many factors. This study systematically describes the survival prognosis and the characteristics of LM from clinicopathological factors to the genetic level. The independent factors that increased the progression risk of patients with LM were KPS ≤ 80, TNBC subtype, grade III, increasing trend of CA153, and DFI ≤ 1 year. Simultaneously, the independent factors that increased the mortality risk of patients with LM were Ki-67 ≥ 30%, grade III, increasing trend of CA153, pain with initial LM, diabetes, and DFI ≤ 1 year. Mutations in the driver genes ESR1, AKT1, ERBB2, and FGFR4, and the MS APOBEC cytidine deaminase, ultraviolet exposure, and defective DNA mismatch repair may provide convenience for LM. We believe that our study makes a significant contribution to the literature because it can help clinicians evaluate the risk of disease progression in patients with BCLM to determine the optimal treatment strategy. Our results also provide a better understanding of the mechanisms underlying BCLM progression and evolution and provide new therapeutic targets that can potentially benefit drug-resistant patients or be eligible for clinical trials.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

ETHICS STATEMENT
The study was approved by the institutional review board of Shandong Cancer Hospital and Institute, and the personal information of all patients and attending doctors was deleted from the data set. The ethics committee waived the requirement of written informed consent for participation.

AUTHOR CONTRIBUTIONS
CT collected the patient's clinicopathological data and downloaded external data sets. SL followed up the patient's survival information. CT and SL jointly analyzed the data and prepared the manuscript. XS and YW supervised the research and revised the manuscript. All authors contributed to the article and approved the submitted version. study design, data collection and analysis, decision to publish, or preparation of the manuscript.

SUPPLEMENTARY MATERIAL
The