B Cells as Prognostic Biomarker After Surgery for Colorectal Liver Metastases

Background: The aim of this study was to identify more accurate variables to improve prognostication of individual patients with colorectal liver metastases (CRLM). Clinicopathological characteristics only partly explain the large range in survival rates.

Methods: MessengerRNA expression profiles of resected CRLM of two patient groups were analysed by mRNA sequencing: poor survivors (death from recurrent disease <30 months after surgery) and good survivors (no recurrent disease >60 months after surgery). Tumour and adjacent liver parenchyma samples were analysed.

Results: MessengerRNA expression profiling of the tumour samples identified 77 genes that were differentially expressed between the two survival groups at a False Discovery Rate (FDR) <0.1. In the adjacent liver parenchyma samples only one gene, MTRNR2L1, showed significantly higher expression in the good survivors. Pathway analysis showed higher expression of immune-related and stroma-related genes in tumour samples from good survivors. Expression data was then validated by immunohistochemistry in two cohorts comprising a total of 125 patients. Immunohistochemical markers that showed to be associated with good survival in the total cohort were: high K/L+ infiltration in tumour stroma [p = 0.029; OR 2.500 (95% CI 1.100–5.682)] and high CD79A+ infiltration in tumour stroma [p = 0.036; OR 2.428 (95%CI 1.062–5.552)].

Conclusions: A high stromal infiltration of CD79A+ B cells and K/L+ plasma cells might be favourable prognostic biomarkers after surgery for CRLM.

Introduction

Colorectal cancer is the cancer with the third highest incidence in Europe, and it often disseminates to the liver (1). Curative treatment by liver surgery is possible in about 10–20% of patients with colorectal liver metastases (CRLM) (2, 3) and survival rates after surgery have improved in the last decades, with reported 5-year survival rates ranging from 14 to 60% (4–6). Clinicopathological characteristics can partly explain this large range, but they are inconsistent in accurately determining the prognosis for the individual patient (7).

Recent advances in technology have led to the notion that molecular characteristics might outperform clinicopathological scoring systems in predicting patient treatment and survival (8). Four individual centres have published prognostic gene signatures for colorectal liver metastases based on expression microarrays, but no single gene was shared between the four signatures (9–12).

As early as 1889, Paget suggested that it is not only the metastatic tumour (the “seed”) but also the hosting organ or tissue (the “soil”) that might be important for metastatic tumour growth (13). Progression of remaining, clinically undetectable metastatic tumour cells in the remnant liver after liver surgery is dependent on both micrometastases and on a favourable micro-environment (14–16). For example, it has been shown that microRNA expression in adjacent liver parenchyma can be associated with recurrent disease and patient survival (17).

In this study, we applied genome-wide mRNA expression profiling by RNAseq in two patient groups selected on survival after liver surgery for CRLM: poor survivors, those who died from recurrences within 30 months after surgery, and good survivors, those alive without recurrences 60 months or more after surgery. Samples from both tumour tissue and adjacent liver parenchyma were included. Immunohistochemistry was performed in two cohorts to validate the RNAseq results on protein level. The goal of this study was to identify molecular markers for favourable patient survival.

Materials and Methods

Patients and Samples

In this study we included patients from two tertiary referral centres for liver surgery, the University Medical Centre Groningen (UMCG, cohort 1) and the Erasmus Medical Centre in Rotterdam (cohort 2). Patients were selected from prospectively maintained databases. Inclusion criteria were (1) R0 partial liver resection for CRLM, (2) no neoadjuvant chemotherapy before liver surgery and no adjuvant chemotherapy after liver surgery, (3) a Fong clinical risk score (5) of 3 or lower, (4) no detectable extrahepatic disease at time of surgery, (5) no other known malignant disease, and (6) availability of fresh frozen (−80°C) resected CRLM material. Of note, (neo)adjuvant chemotherapy is not a standard treatment in the Netherlands. None of the patients received perioperative chemotherapy; neither neoadjuvant nor adjuvant after liver surgery. Adjuvant chemotherapy after primary tumor resection was allowed provided that it was not within 6 months before liver surgery. Follow-up consisted of cross-sectional imaging and measurement of carcinoembryonic antigen (CEA) serum levels every 3–4 months during the first 2 years after liver resection and at 6-month intervals afterwards up to 5 years. Samples of two groups of patients with different survival rates were selected: poor survivors, who died of recurrent disease within 30 months after partial liver resection, and good survivors, who showed no evidence of recurrent disease at 60 months after liver resection. The samples of tumour tissue and adjacent liver parenchyma were reviewed by an experienced hepatopathologist to judge the quality of the tissue. Genome-wide mRNA expression profiling by RNAseq followed by immunohistochemistry was performed on a cohort of patients treated in the University Medical Centre Groningen (cohort 1). Corroboration of immunohistochemical results was performed using tissue samples from patients with similar inclusion criteria who underwent liver surgery in another tertiary referral centre for liver surgery in the Netherlands (cohort 2, Erasmus Medical Centre in Rotterdam). Baseline clinicopathological characteristics of both cohorts are summarised in Table 1.

TABLE 1

Table 1. Cohorts 1 and 2.

RNA/DNA Isolation

Macrodissected frozen samples of tumour tissue and adjacent liver parenchyma were included. Both genomic DNA and total RNA were isolated from 10 μm-tissue sections (RNA/DNA purification kit, Norgen Biotek Corporation, Thorold, Ontario, Canada). DNA/RNA isolation was performed according to the manufacturer's protocol. Quality check and RNA quantification of samples was carried out by capillary electrophoresis using the LabChip GX (Perkin Elmer, Waltham, Massachusetts, USA).

Microsatellite Instability and Mutational Hotspots

Microsatellite instability (MSI) was tested for all cases by amplifying 20 ng genomic DNA using primers for five polymorphic mononucleotide loci (NR21, NR24, BAT25, BAT26, MONO27). The resulting PCR products were analysed on the ABI 3730xl DNA Analyzer (Thermo Fisher, Waltham, Massachusetts, USA). The MSI status was assessed as MSI-instable when two out of five markers showed instability (18). For assessment of somatic driver mutations in KRAS (codon 12 and 13) and BRAF (V600E), genomic DNA was amplified by PCR and the resulting amplicons were analysed by Sanger sequencing. Primers are listed in Table S9.

mRNA Sequencing and Gene Expression Quantification

Sequence libraries were generated using the Quantseq 3' mRNA sample preparation kit (Lexogen, Vienna, Austria) starting from 500 ng total RNA of each sample. Library preparation was performed according to the manufacturer's protocol. Quality check and quantification of libraries was carried out by capillary electrophoresis using the LabChip GX (Perkin Elmer). Barcoded libraries were pooled in equimolar ratios aiming at a final concentration of 2–10 pM. Sequencing was performed on an Illumina HiSeq2500, applying a 50 bp single-read protocol and aiming at 5–10 million reads per sample.

The first 12 nucleotides (nt) were trimmed of the 50nt-reads to remove sequencing artefacts. Hisat version 0.1.5-beta (19) was used to align reads to human genome reference build 37 (20). SAMtools version 1.2 (21) was used to sort the aligned reads, and gene-level quantification was then performed by HTSeq version 0.6.1p1 (22) using “–mode=union,” enabling strandedness. Because of internal poly-A priming bias, only reads mapping within 500 bp upstream from the transcript termination sites—corrected for splicing—were counted, based on Ensembl version 75 gene annotation (23). Reads in genes with multiple transcript termination sites were summed up, resulting in a read count per gene. Genes with an average read count <20 in both survival groups were excluded from further analysis.

Data and Pathway Analysis mRNA Sequencing

Data analysis was performed in R (24). To discover possible bias, a principal component analysis was carried out on all 81 samples after read count normalization by the VSD-function of the DESeq2 package (25). Samples with <300,000 reads in the target regions were excluded from further analysis (Table S1). One sample was removed because of an undesirable patient characteristic (an additional oncological disease). Differential expression analysis using DESeq2 was performed separately for the tumour samples and adjacent liver parenchyma samples comparing the poor vs. good survivors (25). P-values were corrected for multiple testing by the Benjamini Hochberg method (False Discovery Rate, FDR). Heatmap visualizations were made after read count normalization by the regularized log transformation of the DESeq2 package, followed by mean centring per gene (25). Pathway analysis was performed by DAVID EASE software (26, 27), with p-values corrected for multiple testing by the Bonferroni method (Family Wise Error Rate, FWER).

Immunohistochemistry

We included formalin-fixed paraffin-embedded (FFPE) samples containing both tumour and adjacent liver tissue. Besides samples from patients treated at the University Medical Centre Groningen, Groningen, the Netherlands (cohort 1), we also included samples from patients with similar inclusion criteria who were treated at the Erasmus Medical Centre, Rotterdam, the Netherlands (cohort 2). Three μm sections from FFPE tissue blocks were deparaffinised in xylene, rehydrated in graded alcohol and hematoxylin and eosin (H&E) stained. H&E stained slides were examined to confirm the inclusion of the tumour-liver transition area and to score the histopathological growth pattern (28). The Ventana automated staining system (Roche, Basel, Switzerland) was used to stain the tissue sections for CD45 (RP2/18), CD4 (SP35), CD8 (SP57), CD79A (SP18), and Kappa/Lambda (double staining, polyclonal). Staining with the primary antibodies FOXP3 (236A/E7, ABCAM, 1/100 dilution) and SLAMF7 (HPA055945, Atlas Antibodies, 1/200 dilution) was done manually. Antigen retrieval in FOXP3-stained slides was performed in a pressure cooker using a Tris/EDTA buffer (PH 9.0). Antigen retrieval in SLAMF7-stained slides was performed in the microwave using a citrate buffer (PH 6.0). Primary antibodies were diluted with 1%BSA/PBS and incubated at room temperature for 60 min. Secondary and tertiary antibodies were diluted in 1%BSA/PBS and 1% AB serum and incubated at room temperature for 30 min. Appropriate positive and negative controls were used. Diaminobenzidine (DAB) was used as the chromogen, followed by a counterstaining by hematoxylin. The staining of the immunophenotypical markers was graded by microscopic inspection in a semi-quantitative scoring system for three different areas: the invasive margin, the tumour stroma and intra-tumoural region. The invasive margin was defined as the tumour-liver transition area, the tumour stroma was defined as the area of stroma surrounding the tumour cells, and intra-tumoural staining was defined as immunopositive intra-epithelial lymphocytes in tumour cell areas. The grading of immunoreactivity was supervised by two experienced researchers and consensus was achieved in all cases. All markers were scored in grades 1–3.

Statistical Analysis

Summary statistics were obtained using established methods and presented as percentages, median (interquartile range, IQR) or mean (standard deviation, SD). Correlation coefficients of ranked data were calculated using Spearman's rho. Variables associated with survival were first tested by univariable binary logistic regression analyses comparing poor vs. good survivors. Variables with a p-value <0.1 in univariable analysis were entered into the multivariable model. The used multivariable regression model was a binary logistic regression comparing poor vs. good survivors, entering all covariates simultaneously and including a constant in the model. Odds ratio (OR) and 95% confidence intervals (CI) were estimated, and a p-value <0.05 was considered significant. Statistical analyses were carried out with IBM SPSS Statistics V22 (IBM, Armonk, New York, USA).

Results

Quality Control

A summary of the mRNA sequencing quality control data is shown in Table S1. The mean (±SD) number of reads on target was 1.57 10e6 (±0.75 10e6). The variance in the dataset was explored by principal component analysis. Principal component 1 explained 78% of the total variance and showed a high concordance with tissue type (tumour vs. adjacent liver parenchyma, Figure S1). The expression of the liver-specific gene Albumin also resembled the differences in tissue type. Two adjacent liver samples that clustered with the tumour samples had very low levels of Albumin mRNA (0.10 and 0.37% of all reads in Albumin gene, respectively) compared to the mean Albumin mRNA level in the adjacent liver samples (6.21% ± 3.10) (Figure S1). Similarly, three tumour samples that clustered close to adjacent liver samples showed high levels of Albumin mRNA (10.69, 6.87, and 4.25%) compared to the mean Albumin expression in the tumour samples (0.65% ± 1.87) (Figure S1). This suggested a sampling error, and we therefore excluded these five samples from further analyses. Additionally, five samples with fewer than 300,000 reads in the target region were excluded. In total 70 samples were included for further analyses: 39 tumour samples and 31 adjacent liver samples. Of note, for 16 patients both a tumour sample and an adjacent liver sample were available. The additional 38 samples (70 minus 2^*16) are from different patients.

Clinicopathological Characteristics

The clinicopathological characteristics of the 39 patients of whom we included tumour samples are shown in Table S2, stratified by survival. High preoperative CEA was the only statistically different variable between both groups, with a high CEA in the poor survival group (p = 0.049). The clinical risk score, which is a combined score of five clinicopathological factors (5), was also significantly higher in the poor survivors (p = 0.031; Table S2).

The clinicopathological characteristics of the 31 patients with sampled adjacent liver parenchyma are shown in Table S3. There was no clinicopathological characteristic that was statistically different between the good and poor survivors.

Tumour-Specific Expression Profiling

After filtering for low read counts, 8,931 genes were included for further analysis by DESeq2 to identify differentially expressed genes (DEGs) between the poor vs. good survivors. Of these, 333 DEGs were differentially expressed between the two groups with a p <0.01, and for 77 of these DEGs the FDR was <0.1 (Table S4). Figure 1 depicts an unsupervised clustering-based heatmap of the 77 genes with a FDR <0.1, in which a clear separation of the survival groups is observed. To find out if specific biological pathways are over-represented in the 333 genes with p <0.01, we carried out a pathway analysis with DAVID EASE. Two biological entities appeared to be different between the survival groups (Table S5): the extracellular matrix and the immune system. The most significant GO-term-related pathway in each of these two entities was “extracellular matrix” (FWER = 1.17e-10) and “response to external stimulus” (FWER = 0.039), respectively. In general, a higher expression in both the extracellular matrix pathways and the immune-system-related pathways was observed in the good survivors.

FIGURE 1

Figure 1. Heatmap of tumour samples. Unsupervised clustering of the 77 genes with the lowest FDR values in DESeq2 analysis. Samples are shown on the x-axis and the 77 genes on the y-axis. A quantile colour scale is used with 10 different colours ranging from black (low expression) to yellow (high expression).

Correlation With Estimate

Pathway analysis suggested that the expression of extracellular-matrix-related genes and immune-related genes is different between the survival groups (Table S5). We further analysed this using the ESTIMATE (“Estimation of STromal and Immune cells in MAlignant Tumours using Expression data”) software (29). The goal of this software package is to estimate the amount of stromal cells and immune cells by calculating a stromal score and an immune score based on 282 genes. The stromal score was higher in the good survivors (42 ± 603 vs. −469 ± 601; p = 0.013), as was the immune score (725 ± 660 vs. 317 ± 495; p = 0.043). Overall, the differences we observed in the mRNA expression data of the poor vs. good survivors are likely to be caused by the presence of immune cells and stromal cells in the tumour samples.

Adjacent Liver Parenchyma Expression Profiling

Raw read counts of 8,772 genes were used as input for differential mRNA expression analysis by DESeq2. We compared the poor vs. good survivors to find DEGs (Table S6), and identified 109 genes with a p <0.01. Of these, only one gene was differentially expressed between the survival groups with a FDR <0.1: MTRNR2L1 was more highly expressed in the good survivors (Table S6). The relatively high FDRs in this analysis indicate that the results from the adjacent liver parenchyma analysis are not as reliable as the results from the tumour samples. A pathway analysis including the 109 DEGs was performed with DAVID EASE. The transfer RNA pathway was significantly enriched and had, in general, a higher expression in the poor survivors (Table S7). Unsupervised clustering did not show uniform DEGs across the adjacent liver samples, as can be seen in the heatmap (Figure S2). The separation of the survival groups is not as clear as that seen in the tumour samples (Figure 1). Therefore, in the validation, we focused on the mRNA expression data of the tumour samples.

Immunohistochemistry

Tissue slides of two cohorts were analysed, in which cohort 1 consisted of 47 patients (Groningen) and cohort 2 of 78 patients (Rotterdam). The samples in cohort 1 are mostly from the same patients included in the mRNA expression analysis. Table 1 shows that patients in cohort 1 had larger tumours, higher CEA levels, more frequent major liver surgery and a higher clinical risk score; Of note, all these variables are interrelated. Additionally, patients in cohort 2 had an older age.

To study the expression of stroma-related genes, a small pilot experiment was performed in which immunohistochemistry on 11 tumour tissue slides was performed using four stroma-related markers (FBLN1, MMP2, PRRX1, ABLIM1) that were differentially expressed in the mRNA expression analysis. The results did not validate the sequencing results (data not shown). We then hypothesized that the stroma-related differences might resemble the histopathological growth pattern (28). The two most prevalent growth patterns are a desmoplastic growth pattern, which is characterized by a rim of stromal cells between the tumour cells and the adjacent liver cells, and a replacement growth pattern, which lacks a desmoplastic rim. A 100% desmoplastic growth pattern was associated with good survival in both cohort 1 and 2 (p = 0.089 and p = 0.008, respectively; Table 2).

TABLE 2

Table 2. Clinicopathological and biological markers vs. survival in cohorts 1, 2 and combined.

The mRNA expression data prompted us to stain tumour sections to validate immune-related expression, by which we chose seven immunophenotypical markers [CD45, CD4, CD8, FOXP3, CD79A, Kappa/Lambda (K/L) and SLAMF7] to study both T cell and B cell expression. Besides T cell markers we also chose to stain B cell markers as multiple immunoglobulin genes were differentially expressed between the survival groups (Table S5). CD45+ and CD4+ cells generally had a higher abundancy compared to the other markers. Figure 2 shows detailed pictures of the staining. Table 2 shows the clinicopathological characteristics and the scoring of the seven immunophenotypical markers in the poor vs. good survivors of both cohorts separately and combined. The clinical risk score was a good predictor of survival in cohort 1 (p =0.006; Table 2), while in cohort 2 and the combined cohort the CRS was a less good predictor of survival [p = 0.931 and p = 0.055, respectively (Table 2)]. A high CD79A expression in the tumour stroma tended to be associated with good survival in both cohort 1 and 2 (p = 0.069 and p = 0.054, respectively). Similarly, K/L expression in the tumour stroma also tended to be associated with good survival in cohorts 1 and 2 (p = 0.081 and p = 0.070, respectively). A 100% desmoplastic growth pattern tended to be associated with good survival in cohort 1 (p = 0.089) and was significantly associated with good survival in cohort 2 (p = 0.008; Table 2).

FIGURE 2

Figure 2. Detailed pictures of immunohistochemistry. (A) x10 magnification image of high Kappa/Lambda + staining in the tumour stroma. Kappa light chains are stained brown, Lambda light chains are stained pink. (B) x4 magnification image of high CD79A+ staining (brown) in the invasive margin. (C) x10 magnification image of a replacement growth pattern by H&E staining. (D) x10 magnification image of a desmoplastic growth pattern by Masson's trichrome staining. Tumour and liver cells are stained purple, connective tissue (desmoplastic rim) is stained blue. Tu, tumour; Li, liver; St, tumour stroma.

Multivariable Analysis

A multivariable analysis was performed to test whether immunohistochemical markers have significant impact in predicting patient survival. The three most significant biomarkers were analysed in multivariable analysis (Table 3). Immunohistochemical markers that might be associated with good survival in cohorts 1 and cohort 2 were (respectively): high K/L+ infiltration in tumour stroma [p = 0.124; OR 2.682 (95%CI 0.762–15.248] and p = 0.067; OR 2.589 [95%CI 0.934–7.171)], high CD79A+ infiltration in tumour stroma [p = 0.110; OR 2.740 (95%CI 0.795–9.443)] and p = 0.095; OR 2.443 [95%CI 0.856–6.978)], and a 100% desmoplastic growth pattern [p = 0.157; OR 2.709 (95% CI 0.682–10.768) and p = 0.040; OR 5.858 (95%CI 1.083–31.693)]. These three biomarkers are significantly associated with patient survival (all p < 0.05) if we combine data from cohort 1 and 2 (Table 3). If we correct only for the clinicopathological factor CRS (high vs. low) in the total cohort, high K/L+ infiltration in tumour stroma [p = 0.026; OR 2.359 (95% CI 1.109–5.019) as high CD79A+ infiltration in tumour stroma (p = 0.031; OR 2.338 (95% CI 0.1.081–5.056)] as a 100% desmoplastic growth pattern [p = 0.007; OR 3.952 (95% CI 1.454–10.742)] were all significantly associated with favourable survival whereas the CRS was not (p = 0.168; p = 0.103; p = 0.163, respectively).

TABLE 3

Table 3. Multivariable analysis to predict survival.

Statistical correction for clinicopathological factors did not uniformly improve the prognostic value of biomarkers in multivariable analysis (30). This is partly explained by an association between the clinical marker tumour size and protein expression (Table S8). Of note, in patients with more than one metastasis (n = 36; 29% of all cases), the tumour with the largest diameter was used for immunohistochemical staining in 28 out of 36 patients. A desmoplastic growth pattern was associated with smaller tumour size [2.8 cm (1.5–4.9) vs. non-desmoplastic growth pattern 3.8 cm (2.6–5.1); p = 0.029]. In addition, high stromal CD79A+ infiltration was also associated with smaller tumour size [3.4 cm (2.0–4.5) vs. low CD79A+ infiltration 4.1 cm (3.0–7.0); p = 0.011]. Especially B cell-related markers scored in the invasive margin were associated with size of the tumour (Table S8).

Discussion

In this study we used mRNA sequencing to identify prognostic molecular markers in patients after liver surgery for colorectal liver metastases. We selected two patient groups: poor survivors (death due to recurrence within 30 months) and good survivors (disease-free survival > 60 months after liver surgery). We show that patients with good survival had a higher expression of immune-related and stroma-related genes. Additional analysis by ESTIMATE software indicated higher immune cell and stromal cell infiltration in the tumour samples of good survivors (29). Immunohistochemistry showed in two cohorts that high immune infiltration of CD79A+ B cells and K/L+ plasma cells in tumour stroma tended to be associated with good survival. Comparable results from these biomarkers were observed in multivariable analysis correcting for clinicopathological factors. In addition, the desmoplastic growth pattern, in which a stroma-rich pseudo-capsule around the tumour is present, might also be associated with good patient survival. The prognostic value of a desmoplastic growth pattern has been reported previously (31–36). In contrast, the associations between good patient survival and high stromal infiltration of CD79A+ cells and K/L+ cells are, to the best of our knowledge, novel findings.

The top differentially expressed genes in the tumour samples in this study did not overlap with those of the four previously published expression signatures (9–12). However, no individual gene was shared between these earlier studies either (9–12). This may partly be explained by the different experimental methods used by the studies. We used mRNA sequencing, while the earlier studies used different expression microarrays (9–12, 37). The inclusion criteria for patients were also very different between studies. Our study and that of van der Stok et al. (11) did not administer neoadjuvant or adjuvant chemotherapy to included patients, while the three other studies did (9, 10, 12). Neoadjuvant chemotherapy influences tumour biology and might thereby alter gene expression levels (38). Three of the studies aimed at identifying genes that are associated with recurrence rather than with overall survival (10–12), and the cut-off values to stratify rapid recurrence vs. late or no recurrence differed per study: Snoeren et al. (10) included patients with DFS ≤ 1 year vs. DFS >1 year, Snoeren et al. (13) included patients with DFS <6 months vs. DFS > 2 years, and van der Stok et al. included patients with DFS ≤ 1 vs. DFS >3 years. In our opinion we used the most straight-forward inclusion by selecting chemotherapy-naïve patients, the strongest endpoint (disease-related overall survival; death of recurrent disease within 30 months vs. alive and free of disease 5 years after liver surgery) and the most accurate gene expression technology. These inter-study differences might explain the lack of shared genes.

We also studied adjacent liver parenchyma and found a significantly higher expression of MTRNR2L1 in the good survivors. This gene is an isoform of MT-RNR (humanin), which is reported to have neuroprotective and anti-apoptotic functions and is mainly studied in age-related illness like Alzheimer's disease (39, 40).

In our study, the B-cell-related immunohistochemical markers CD79A and K/L suggest an association with patient survival. There are no previous studies that analysed CD79A and K/L in CRLM, but the prognostic value of B-cell marker CD20 was analysed by four other studies (41–44). Two studies reported an association of high CD20+ B cells in the invasive margin with a favourable overall survival (41, 44). In addition, one of these studies also showed an association between high CD20+ B cells in the tumour stroma and a favourable survival (41). Two other studies observed no associations with patient survival (42, 43). The results in literature are diverse when focussing on the role of the immune system on prognosis after surgery for CRLM (45). Although the role of tumour-infiltrating B cells and immunoglobulins is not as well-described as the role of tumour-infiltrating T cells, an association between B cell infiltration and favourable patient survival has been reported for other solid tumours (46–48). Anti-tumour immunity can be achieved through immunoglobulin production by plasma cells and, moreover, through increased cytotoxic T cell activation via antigen presenting B cells (46–48). In contrast, others have questioned the actual effect of the B cell response (47, 48). It is also proposed that the B cell infiltration is the result of IFNγ production and serves as a surrogate marker of the T-cell-mediated anti-tumour response (47).

Conclusions

In conclusion, our study of the liver metastasis of patients with colorectal cancer showed that high stromal infiltration of CD79A+ B cells and high stromal infiltration of K/L+ plasma cells might be favourable prognostic biomarkers after surgery for CRLM. Future evaluation on external cohorts is needed to prove whether these biomarkers truly are associated with patient survival and of practical clinical value.

Data Availability Statement

Datasets are in a publicly accessible repository: The generated and analysed RNA sequencing dataset for this study can be found in the European Nucleotide Archive (study Accession No: PRJEB35154).

Ethics Statement

This study was approved by the institutional review board of the University Medical Center Groningen, Groningen, the Netherlands (research registry no. 201501288). Written informed consent was not required in accordance with national legislation and the institutional requirements.

Author Contributions

JH, AG, DG, CV, RS, KJ, and KK: conceptualization and methodology. MT and KK: software. JH, LV, AG, DH, and PN: validation. JH and KK: formal analysis. JH, KK, and KJ: investigation and data curation. RS, KJ, and CV: resources. JH: writing—original draft preparation. LV, AG, DH, PN, AG, CV, RS, KJ, and KK: writing—review and editing. JH, KK, LV, and AG: visualization. LV, AG, CV, RS, KJ, and KK: supervision. KJ: project administration and funding acquisition. The samples of tumour tissue and adjacent liver parenchyma were reviewed on tissue quality by AG, experienced hepatopathologist and co-author in this study.

Funding

This study was supported by two grants from the JK de Cock foundation to JH.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Kate McIntyre for assistance in the editorial process.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2020.00249/full#supplementary-material

References

1. Ferlay J, Steliarova-Foucher E, Lortet-Tieulent J, Rosso S, Coebergh JWW, Comber H, et al. Cancer incidence and mortality patterns in Europe: estimates for 40 countries in (2012). Eur J Cancer. (2013) 49:1374–403. doi: 10.1016/j.ejca.2012.12.027

PubMed Abstract | CrossRef Full Text | Google Scholar

2. van der Pool AEM, Damhuis RA, IJzermans JNM, de Wilt JHW, Eggermont AMM, Kranse R, et al. Trends in incidence, treatment and survival of patients with stage IV colorectal cancer: a population-based series. Color Dis. (2012) 14:56–61. doi: 10.1111/j.1463-1318.2010.02539.x

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Manfredi S, Lepage C, Hatem C, Coatmeur O, Faivre J, Bouvier A-M. Epidemiology and management of liver metastases from colorectal cancer. Ann Surg. (2006) 244:254–9. doi: 10.1097/01.sla.0000217629.94941.cf

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Simmonds PC, Primrose JN, Colquitt JL, Garden OJ, Poston GJ, Rees M. Surgical resection of hepatic metastases from colorectal cancer: a systematic review of published studies. Br J Cancer. (2006) 94:982–99. doi: 10.1038/sj.bjc.6603033

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Fong Y, Fortner J, Sun RL, Brennan MF, Blumgart LH. Clinical score for predicting recurrence after hepatic resection for metastatic colorectal cancer: analysis of 1001 consecutive cases. Ann Surg. (1999) 230:309–18; discussion 318–21. doi: 10.1097/00000658-199909000-00004

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Hof J, Wertenbroek MWJLAE, Peeters PMJG, Widder J, Sieders E, de Jong KP. Outcomes after resection and/or radiofrequency ablation for recurrence after treatment of colorectal liver metastases. Br J Surg. (2016) 103:1055–62. doi: 10.1002/bjs.10162

CrossRef Full Text | Google Scholar

7. Moris D, Pawlik TM. Personalized treatment in patients with colorectal liver metastases. J Surg Res. (2017) 216:26–9. doi: 10.1016/j.jss.2017.04.013

PubMed Abstract | CrossRef Full Text | Google Scholar

8. El-Deiry WS, Goldberg RM, Lenz H, Shields AF, Gibney GT, Tan AR, et al. The current state of molecular testing in the treatment of patients with solid tumors, (2019). CA Cancer J Clin. (2019) 69:305–43. doi: 10.3322/caac.21560

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Ito H, Mo Q, Qin L-X, Viale A, Maithel SK, Maker AV, et al. Gene expression profiles accurately predict outcome following liver resection in patients with metastatic colorectal cancer. PLoS ONE. (2013) 8:e81680. doi: 10.1371/journal.pone.0081680

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Snoeren N, van Hooff SR, Adam R, van Hillegersberg R, Voest EE, Guettier C, et al. Exploring gene expression signatures for predicting disease free survival after resection of colorectal cancer liver metastases. PLoS ONE. (2012) 7:e49442. doi: 10.1371/journal.pone.0049442

PubMed Abstract | CrossRef Full Text | Google Scholar

11. van der Stok EP, Smid M, Sieuwerts AM, Vermeulen PB, Sleijfer S, Ayez N, et al. mRNA expression profiles of colorectal liver metastases as a novel biomarker for early recurrence after partial hepatectomy. Mol Oncol. (2016) 10:1542–50. doi: 10.1016/j.molonc.2016.09.002

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Snoeren N, Emmink BL, Koerkamp MJG, van Hooff SR, Goos JACM, van Houdt WJ, et al. Maspin is a marker for early recurrence in primary stage III and IV colorectal cancer. Br J Cancer. (2013) 109:1636–47. doi: 10.1038/bjc.2013.489

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Paget S. The distribution of secondary growths in cancer of the breast. 1889. Cancer Metastasis Rev. (1989) 8:98–101.

PubMed Abstract | Google Scholar

14. Kaplan RN, Riba RD, Zacharoulis S, Bramley AH, Vincent L, Costa C, et al. VEGFR1-positive haematopoietic bone marrow progenitors initiate the pre-metastatic niche. Nature. (2005) 438:820–7. doi: 10.1038/nature04186

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Marusawa H, Jenkins BJ. Inflammation and gastrointestinal cancer: an overview. Cancer Lett. (2014) 345:153–6. doi: 10.1016/j.canlet.2013.08.025

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Milette S, Sicklick JK, Lowy AM, Brodt P. Molecular pathways: targeting the microenvironment of liver metastases. Clin Cancer Res. (2017) 23:6390–9. doi: 10.1158/1078-0432.CCR-15-1636

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Kahlert C, Klupp F, Brand K, Lasitschka F, Diederichs S, Kirchberg J, et al. Invasion front-specific expression and prognostic significance of microRNA in colorectal liver metastases. Cancer Sci. (2011) 102:1799–807. doi: 10.1111/j.1349-7006.2011.02023.x

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Kawakami H, Zaanan A, Sinicrope FA. Microsatellite instability testing and its role in the management of colorectal cancer. Curr Treat Options Oncol. (2015) 16:30. doi: 10.1007/s11864-015-0348-2

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. (2015) 12:357–60. doi: 10.1038/nmeth.3317

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Genomes Project Consortium, Auton A, Brooks LD, Durbin RM, Garrison EP, Kang HM, et al. A global reference for human genetic variation. Nature. (2015) 526:68–74. doi: 10.1038/nature15393

CrossRef Full Text | Google Scholar

21. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. (2009) 25:2078–9. doi: 10.1093/bioinformatics/btp352

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Anders S, Pyl PT, Huber W. HTSeq–a Python framework to work with high-throughput sequencing data. Bioinformatics. (2015) 31:166–9. doi: 10.1093/bioinformatics/btu638

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Flicek P, Amode MR, Barrell D, Beal K, Billis K, Brent S, et al. Ensembl 2014. Nucleic Acids Res. (2014) 42:D749–55. doi: 10.1093/nar/gkt1196

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Team RC. R: A Language and Environment for Statistical Computing. R Found Stat Comput. Vienna, Austria (2016).

PubMed Abstract | Google Scholar

25. 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: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Huang DW, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. (2009) 37:1–13. doi: 10.1093/nar/gkn923

CrossRef Full Text | Google Scholar

27. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. (2009) 4:44–57. doi: 10.1038/nprot.2008.211

CrossRef Full Text | Google Scholar

28. van Dam P-J, van der Stok EP, Teuwen L-A, Van den Eynden GG, Illemann M, Frentzas S, et al. International consensus guidelines for scoring the histopathological growth patterns of liver metastasis. Br J Cancer. (2017) 117:1427–41. doi: 10.1038/bjc.2017.334

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. (2013) 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Peduzzi P, Concato J, Kemper E, Holford TR, Feinstem AR. A simulation study of the number of events per variable in logistic regression analysis. J Clin Epidemiol. (1996) 49:1373–9. doi: 10.1016/S0895-4356(96)00236-3

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Nielsen K, Rolff HC, Eefsen RL, Vainer B. The morphological growth patterns of colorectal liver metastases are prognostic for overall survival. Mod Pathol. (2014) 27:1641–8. doi: 10.1038/modpathol.2014.4

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Pinheiro RS, Herman P, Lupinacci RM, Lai Q, Mello ES, Coelho FF, et al. Tumor growth pattern as predictor of colorectal liver metastasis recurrence. Am J Surg. (2014) 207:493–8. doi: 10.1016/j.amjsurg.2013.05.015

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Brunner SM, Kesselring R, Rubner C, Martin M, Jeiter T, Boerner T, et al. Prognosis according to histochemical analysis of liver metastases removed at liver resection. Br J Surg. (2014) 101:1681–91. doi: 10.1002/bjs.9627

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Eefsen RL, Vermeulen PB, Christensen IJ, Laerum OD, Mogensen MB, Rolff HC, et al. Growth pattern of colorectal liver metastasis as a marker of recurrence risk. Clin Exp Metastasis. (2015) 32:369–81. doi: 10.1007/s10585-015-9715-4

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Nyström H, Tavelin B, Björklund M, Naredi P, Sund M. Improved tumour marker sensitivity in detecting colorectal liver metastases by combined type IV collagen and CEA measurement. Tumour Biol. (2015) 36:9839–47. doi: 10.1007/s13277-015-3729-z

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Van den Eynden GG, Bird NC, Majeed AW, Van Laere S, Dirix LY, Vermeulen PB. The histological growth pattern of colorectal cancer liver metastases has prognostic value. Clin Exp Metastasis. (2012) 29:541–9. doi: 10.1007/s10585-012-9469-1

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Balachandran VP, Arora A, Gönen M, Ito H, Turcotte S, Shia J, et al. A validated prognostic multigene expression assay for overall survival in resected colorectal cancer liver metastases. Clin Cancer Res. (2016) 22:2575–82. doi: 10.1158/1078-0432.CCR-15-1071

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Misale S, Di Nicolantonio F, Sartore-Bianchi A, Siena S, Bardelli A. Resistance to anti-EGFR therapy in colorectal cancer: from heterogeneity to convergent evolution. Cancer Discov. (2014) 4:1269–80. doi: 10.1158/2159-8290.CD-14-0462

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Gong Z, Tas E, Muzumdar R. Humanin and age-related diseases: a new link? Front Endocrinol. (2014) 5:210. doi: 10.3389/fendo.2014.00210

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Bodzioch M, Lapicka-Bodzioch K, Zapala B, Kamysz W, Kiec-Wilk B, Dembinska-Kiec A. Evidence for potential functionality of nuclearly-encoded humanin isoforms. Genomics. (2009) 94:247–56. doi: 10.1016/j.ygeno.2009.05.006

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Mlecnik B, Van den Eynde M, Bindea G, Church SE, Vasaturo A, Fredriksen T, et al. Comprehensive intrametastatic immune quantification and major impact of immunoscore on survival. J Natl Cancer Inst. (2018) 110:97–108. doi: 10.1093/jnci/djx123

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Berthel A, Zoernig I, Valous NA, Kahlert C, Klupp F, Ulrich A, et al. Detailed resolution analysis reveals spatial T cell heterogeneity in the invasive margin of colorectal cancer liver metastases associated with improved survival. Oncoimmunology. (2017) 6:e1286436. doi: 10.1080/2162402X.2017.1286436

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Tanis E, Julié C, Emile J-F, Mauer M, Nordlinger B, Aust D, et al. Prognostic impact of immune response in resectable colorectal liver metastases treated by surgery alone or surgery with perioperative FOLFOX in the randomised EORTC study 40983. Eur J Cancer. (2015) 51:2708–17. doi: 10.1016/j.ejca.2015.08.014

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Meshcheryakova A, Tamandl D, Bajna E, Stift J, Mittlboeck M, Svoboda M, et al. B cells and ectopic follicular structures: novel players in anti-tumor programming with prognostic power for patients with metastatic colorectal cancer. PLoS ONE. (2014) 9:e99008. doi: 10.1371/journal.pone.0099008

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Hof J, Kok K, Sijmons RH, de Jong KP. Systematic review of the prognostic role of the immune system after surgery of colorectal liver metastases. Front Oncol. (2019) doi: 10.3389/fonc.2019.00148

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Flynn NJ, Somasundaram R, Arnold KM, Sims-Mourtada J. The multifaceted roles of B cells in solid tumors: emerging treatment opportunities. Target Oncol. (2017) 12:139–152. doi: 10.1007/s11523-017-0481-x

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Guy TV, Terry AM, Bolton HA, Hancock DG, Shklovskaya E, Fazekas de St Groth B. Pro- and anti-tumour effects of B cells and antibodies in cancer: a comparison of clinical studies and preclinical models. Cancer Immunol Immunother. (2016) 65:885–96. doi: 10.1007/s00262-016-1848-z

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Yuen GJ, Demissie E, Pillai S. B lymphocytes and cancer: a love-hate relationship. Trends Cancer. (2016) 2:747–57. doi: 10.1016/j.trecan.2016.10.010

PubMed Abstract | CrossRef Full Text | Google Scholar