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)(5)(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)(10)(11)(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)(15)(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.

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

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 MSIinstable 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 tumourliver 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 intraepithelial 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).

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 Frontiers in Oncology | www.frontiersin.org 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.

Correlation With Estimate
Pathway analysis suggested that the expression of extracellularmatrix-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 stromarelated 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).
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).

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 (   Three biomarkers were analysed in multivariable analysis. Each biomarker is adjusted for the most significant prognostic clinical factor(s) within the respective cohort according to the "one variable per 10 events" rule (30). In cohort 1, the variables were corrected for "major liver surgery (≥ 3 segments)". In cohort 2, correction was performed using the variables: age, sex and tumour size. In the combined cohort, correction was performed using the variables: age, sex, solitary tumour, tumour size, and cohort. An odds ratio >1 corresponds with good survival. CI, confidence interval; K/L, Kappa Lamda. in multivariable analysis (30). This is partly explained by an association between the clinical marker tumour size and protein expression (   (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)(32)(33)(34)(35)(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)(10)(11)(12). However, no individual gene was shared between these earlier studies either (9)(10)(11)(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)(10)(11)(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)(11)(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 straightforward 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)(42)(43)(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)(47)(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)(47)(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.

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