Exploration of predictive biomarkers for postoperative recurrence of stage II/III colorectal cancer using genomic sequencing

Abstract Postoperative recurrence of colorectal cancer (CRC) eventually leads to therapeutic failure; therefore, treatment strategies based on accurate prediction of recurrence are urgently required. To identify biomarkers that can predict treatment outcomes, we compared the mutational profiles of surgically resected specimens from patients with recurrent cancer with those from patients with non‐recurrent cancer. Target sequencing, whole‐exome sequencing (WES), or whole‐genome sequencing (WGS) was performed on 89 and 58 tumors from recurrent and non‐recurrent cases, respectively. WGS revealed the driver mutations that were not detected with target sequencing or WES, including the structural variations affecting ZFP36L2. Loss of function of ZFP36L2 was frequently observed in primary tumors from recurrent cases. Furthermore, the recurrence‐free survival of patients with loss of function of ZFP36L2 was significantly shorter relative to patients with no loss of ZFP36L2 function. In summary, the study demonstrated that detailed genomic analysis could help improve precision medicine for CRC.


| INTRODUCTION
Colorectal cancer (CRC) is the second leading cause of cancer-related deaths in Japan (Vital Statistics Japan, Ministry of Health, Labor and Welfare, https://www. mhlw.go.jp/engli sh/datab ase/db-hw/vs01.html). While radical resection is the first line of treatment, postoperative adjuvant chemotherapy is usually administered to patients with stage III CRC. Although the postoperative recurrence rate of stage II CRC is approximately 15%, administration of chemotherapy in such cases is generally left to the discretion of the clinician, since its efficacy for stage II CRC is still unknown. [1][2][3] Extensive research has elucidated the landscape of genomic aberrations in CRC, [4][5][6][7] enabling the implementation of precision medicine through therapeutic stratification; however, the application of prognostic prediction based on mutational profiles remains limited. For instance, low expression of CDX2 is a predictive marker for chemosensitivity, 8 mutations in PTPRT serve as a potential predictive marker for metastatic CRC, 7 mutations in BRAF V600E is a prognostic factor for stage II microsatellite stable (MSS) CRC, 9 and mutations in KRAS or BRAF are associated with worse outcomes for patients with stage III MSS CRC. 10,11 However, the prediction of postoperative recurrence is unsatisfactory. Many challenges remain in improving treatment strategies based on such prediction.
In this study, we aimed to identify additional biomarkers by comparing the mutational profiles of surgically resected specimens from patients with recurrent cancer relative to those from patients with non-recurrent cancer. Our findings could help improve the prediction of treatment outcomes in CRC and increase the efficacy of therapeutic strategies.

| Clinical specimens
Tissues containing tumor and corresponding normal mucosa were collected from surgically resected specimens obtained from patients treated at The University of Tokyo Hospital between 2006 and 2019. Specimens were either snap-frozen in liquid nitrogen immediately after resection and stored at −80°C or immersed in RNAlater in TissueProtect Tubes (Qiagen) overnight at 4°C, and stored at −20°C until use. Genomic DNA was extracted from tissue samples using the DNeasy Blood & Tissue Kit or the QIAamp DNA Mini Kit (Qiagen). Total RNA was extracted with RNA-Bee (Tel-Test) and treated with DNase I (Qiagen) using the RNeasy Mini Kit (Qiagen).

| Target sequencing
Two hundred and fifty ng genomic DNA from each sample was fragmented to 200-bp lengths using the LE220 Focused ultrasonicator (Covaris). The End Prep Enzyme Mix (New England Biolabs) was used for fragment-end repair and the NEBNext Adaptor for Illumina (New England Biolabs) was used for adaptor ligation. Adaptor-ligated DNA fragments were amplified for nine cycles. Seven hundred and fifty ng of each PCR product was subjected to target capture with custom-made probes, as described previously. 12 The captured library was amplified over 11 cycles, and massive parallel sequencing of the isolated fragments was performed with the HiSeq2500 (Illumina,) using the paired-end option.

| Whole-exome sequencing
Two hundred and fifty ng of genomic DNA from each sample was fragmented to 200-bp lengths using the Focused ultrasonicator (Covaris). Adaptor-ligated DNA fragments, obtained as described in the previous sub-section, were amplified over nine cycles, and 750 ng of each PCR product was used to capture the exonic region using the SureSelect All Exon Kit 6 (Agilent Technologies) and the Sciclone G3 automation system (PerkinElmer). The captured library was amplified over 11 cycles, and massive parallel sequencing of isolated fragments was performed with the HiSeq2500 (Illumina) or the NovaSeq6000 (Illumina).

| Whole-genome sequencing
Five hundred ng of genomic DNA from each sample was used for library preparation using the automated library preparation run on Bravo (Agilent Technologies) with the TruSeq DNA PCR-Free Library Prep Kit (Illumina). IDT for Illumina-TruSeq UD Indexes (Illumina) was used as an unrelated adapter for index reads. Massive parallel sequencing of the isolated fragments was performed with the NovaSeq6000 (Illumina).

| Sequencing analysis
Reads from target sequencing, paired-end WES, and WGS were independently aligned to the human reference genome (hg38) using the Burrows-Wheeler Aligner (http://bio-bwa. sourc eforge.net/) and Bowtie 2 (http://bowti e-bio.sourc eforge.net/bowti e2/index.shtml). Somatic (synonymous and non-synonymous) mutations were called using an in-house caller and two publicly available mutation callers: Genome Analysis Toolkit (https://gatk.broad insti tute.org/hc/en-us) MuTect 2 (https://gatk.broad insti tute.org/hc/en-us/artic les/36003 75938 51-Mutect2) and VarScan 2 (http://varsc an.sourc eforge.net/). Mutations were discarded if any of the following criteria were met: total read number < 20, variant allele frequency in tumor samples <0.05, mutant read number in the germline control samples >2, a mutation in only one genome strand, a variant present in normal human genomes of the 1000 Genomes Project dataset (https:// www.inter natio nalge nome.org/) or the in-house database. Gene mutations were annotated using SnpEff (http://snpeff.sourc eforge.net). The copy number status was analyzed using an in-house pipeline that determined the logR ratio (LRR) as follows: (I) homozygous (variant allele frequency [VAF] ≤ 0.05 or ≥ 0.95) or heterozygous (VAF 0.4-0.6) SNPs were selected from the genomes of the related normal samples in the 1000 Genomes Project database; (II) normal and tumor read depths for the selected SNPs were adjusted based on the G + C percentage of a 100-bp window flanking position; (III) LRR was calculated as log 2 t i n i , where n i and t i are normal and tumor-adjusted depths at position i, respectively; and (IV) each representative LRR was determined by using the median of a moving window (1 Mb) centered at position i. The LRR of the copy number for both alleles, i.e., the major and the minor alleles were determined for every region of the genome. The p values for copy number gain or loss of the respective genomic regions were determined from LRRs with a permutation test (100,000 iterations) in GISTIC2. (https://softw are.broad insti tute.org/cance r/cga/ gistic). Q values were calculated from the P values using the 'qvalue' package in R (http://github.com/jdsto rey/qvalue). The allele-specific copy number status was inferred using the LRR in FACETS. 13 Structural variations (SVs) were detected by analyzing WGS data with the Genomon 2.6.3 pipeline (https://github.com/Genom on-Project).

| Signature analysis
Mutational signatures were analyzed using the Wellcome Trust Sanger Institute Mutational Signature Framework (http://jp.mathw orks.com/matla bcent ral/filee xchan ge/38724 -wtsi-mutat ional -signa ture-frame work). The optimal number of signatures was determined based on the signature stabilities and average Frobenius reconstruction errors.

| Clinical information and survival analysis
Clinical information, including sex, age, preoperative diagnosis, site of the lesion, the experience of preoperative chemoradiotherapy, pathological findings, recurrence, overall survival, and recurrence-free survival, were obtained from all patients included in the study. Clinical data were anonymized at the Department of Healthcare Information Management, The University of Tokyo.
Patients were assessed after operation at every 3 months, for as long as possible. Patients underwent physical examination, review of symptoms suggestive of recurrence, semi-yearly computed tomography of chest, abdomen, and pelvis, measuring carcinoembryonic antigen and carbohydrate antigen 19-9 at each visit, and colonoscopy at years 1, 2, 3, 4 and 5. RFS is determined at the first event of recurrent disease (local, distant or regional). OS is determined at the time of decease. Patients without a recurrence or decease were censored at their last assessment for follow-up.

| Clonal analysis
The cellular prevalence of clones carrying individual nonsynonymous mutations was determined using PyClone. 16 The determined cellular prevalence was used as the input data for LiCHeE 17 to identify phylogenetic relationships between the clones. Clinicians estimated the mode of recurrence by considering the location of the lesion, the time elapsed between the first and second surgeries, and the surgical procedure. Concordance between the results of clonal analysis and inference based on clinical information was evaluated.

| Statistics
Numerical variables were summarized by the median and range. Comparisons of numerical variables between groups were performed using the two-tailed Wilcoxon rank-sum test. Comparisons of categorical variables between groups were performed using Fisher's exact test. The survival curve was established by the Kaplan-Meier method and compared using the log-rank test with Bonferroni adjustment. The Cox proportional-hazards model was also used to evaluate the effects of multiple variables. Statistical analysis was performed using R 3.5.1 (R Core Team, 2015) or SciPy 1.1.0 (Community Library Project, 2021).

| Overview of patients
From 2537 colorectal cancer (CRC) specimens surgically resected in the University of Tokyo hospital between January 5, 2006, and January 31, 2019, recurrent CRCs (rCRCs) and stage IV CRCs for which fresh frozen tumor tissue was available were selected for this study ( Figure 1); 89 tumors from 68 patients that showed recurrence within 5 years after radical surgery, were selected as rCRCs. As a control group, 58 tumors from 49 patients, who remained free from recurrence for at least 5 years, were selected as non-recurrent CRCs (nrCRCs). The nrCRCs were chosen in order from newest to oldest, so that they would have the same number as the rCRCs. Three samples (one tumor with high microsatellite instability, one associated with Lynch syndrome, and one associated with ulcerative colitis) were excluded from the analysis. In addition, we included 43 tumor tissues from 42 patients with stage IV CRC at the pre-operative diagnosis. In total, we analyzed 190 samples from 159 patients. Multiple samples were collected from 22 patients (52 samples).
Patients' characteristics are presented in Table 1; the mean age at surgery was 63 years, and 50.9% of the patients were male. The median follow-up period was 1811 days (67-4620 days), whereas the median period for recurrence after surgery was 906 days (44-2537 days). Compared to those in the cohort from The Cancer Genome Atlas (TCGA) program, the clinical features in our study did not show any difference.

| Genomic analysis of rCRC
Target sequencing was conducted using 43 tumors from 28 patients with rCRC and six tumors from three patients with nrCRC. We selected 26 tumors from 18 patients with rCRC and five tumors from three patients with nrCRC, based on the relatively higher tumor content estimated from target sequencing, for analysis via WGS. We selected 46 tumors from 40 patients with rCRC, 52 tumors from 46 patients with nrCRC, and 43 tumors from 42 patients with stage IV CRC for analysis via WES, and 78 tumors from 59 patients with rCRC, 53 tumors from 43 patients with nrCRC, and 39 tumors from 39 patients with stage IV CRC for analysis via transcriptome sequencing ( Figure 1).
The results of target sequencing, WES, and WGS are shown in Figure 2 and Figure S1. The distribution of CMS categories was similar across patients with nrCRC, rCRC, and stage IV CRC. Our findings did not agree very well with a previous report of prognosis based on CMS classification in that the outcomes of patients with CMS4 tumors in this study were not particularly worse than those reported in the original report. 14 The discrepancy might be due to ethnic differences or how the cohort in this study was selected from patients who were eligible for first-line curable resection.
The detected single-nucleotide variants were concordant with those reported in previous genomic studies of CRCs. 5,6 In our cohort, we did not detect any hypermutator (tumors having ≥10 mutations/Mb), microsatellite instability or mutation in POLD1 and POLE, or germline mutation in APC and MUTYH. We identified somatic mutations in APC (

| Copy number alterations in rCRC
On investigating the copy number alterations in nr-CRCs and rCRCs, we noticed that copy number gains in chr 1: 144,000,000-164,000,000 were more frequent in nrCRCs ( Figure 3A). Survival analysis indicated that CNAs in this genomic region may be prognostic factors ( Figure 3B). Although there were several genes in this region related to colorectal cancer including BCL9 20 and DDR2, 21,22 we could not identify upregulated genes as a result of copy number gain in this region. Next, we analyzed the molecular pathways related to DEGs between cases with and without copy number gain in this region ( Figures 3C,D) and found immune response to be negatively correlated with copy number gain in this region, and likely related to worse outcomes. By contrast, biological pathways associated with protein processing were upregulated in tumors with copy number gain in this region.

| Structural variations (SVs) in rCRC
Analysis of WGS data revealed a median of 76 (7-298) SVs, a fraction of which affected the major oncogenes ( Figure 4A). SVs in APC were identified in two cases (deletion and inversion in tumor number 1003 and deletion in tumor number 1061); SV in KRAS was found in one case (tandem duplication in tumor number 1018); SVs in FBXW7 were found in two cases (deletion in tumor number 1006 and tandem duplication in tumor number 1052) and SV in PIK3CA was found in one case (deletion in tumor number 1052). In tumor number 1003, a large deletion was found in one allele of APC, resulting in loss of heterozygosity, whereas the other allele had multiple SVs, including a deletion in the exonic region. In tumor number 1061, a large deletion was found in one allele of APC, resulting in loss of heterozygosity, whereas the other allele had a disruptive frameshift deletion that probably eliminated APC function. In tumor number 1006, a deletion was found in the exon region of FBXW7. In tumor number 1052, a deletion was found in PIK3CA that resulted in the truncation of its transcript. In tumor number 1019, a deletion was found in the second exon of ZFP36L2 that resulted in truncation, whereas in tumor number 1031, we detected a deletion encompassing the entire region of ZFP36L2 and another in the promoter. Candidates for fusion genes were not identified in this study.    Figure 4B). By contrast, TP53 was more frequently altered by point mutations accompanied by loss of heterozygosity (88/97, 90.7%). In one tumor analyzed via WGS, two distinct SVs were found to affect the two APC alleles. Since a substantial number of SVs affecting major oncogenic events were not identified by WES, analysis of SVs via WGS was considered to facilitate the detection of oncogenic abnormalities that were not identified by target sequencing or WES.

| Clinical significance of mutations in ZFP36L2
As described above (Figure 2), mutations and SVs in ZFP36L2 were recurrently observed. We confirmed the sequencing results with manual inspection ( Figure S2 and Table S1). Since a substantial fraction of mutations were frameshifts ( Figure 5A), the detected mutations and SVs were considered to impair the function of ZFP36L2 in most cases, suggesting ZFP36L2 as a tumor suppressor in CRC. Since it appeared that mutations in ZFP36L2 were enriched among rCRCs compared with among nr-CRCs, we performed a survival analysis to explore the possibility that mutations in ZFP36L2 were associated with tumor recurrence. Interestingly, RFS of patients with ZFP36L2 mutations was significantly shorter than that of patients without the mutations ( Figure 5B, Table S2). Non-synonymous mutations in ZFP36L2 were found in 30 of 594 cases from the TCGA cohort. Of the 223 cases with disease-free survival data available, 10 had nonsynonymous mutations in ZFP36L2. However, no significant difference was detected in the disease-free survival of patients with or without ZFP36L2 mutations in the TCGA cohort ( Figure S3).
DEGs detected on comparing tumors with or without ZFP36L2 mutations were subjected to pathway analysis. Interestingly, tumors with ZFP36L2 mutations had increased expression of genes associated with metabolic pathways such as cholesterol metabolism and mitochondrial function ( Figure 5C). This finding suggests that tumors with ZFP36L2 mutations might be characterized by a different metabolic state, making metabolic pathways a potential therapeutic target. However, further study is required. Expression of genes associated with immune response was suppressed in tumors with ZFP36L2 mutations ( Figure 5D), which might contribute to the poor outcome of patients with these tumors.
Although mutations in KRAS and BRAF are reportedly related to inferior treatment outcomes, 10,11 we did not observe significantly shorter survival in patients with KRAS or BRAF mutations, probably due to relatively small cohort size ( Figure 5E,F). While 71% (10 of 14) of mutations in ZFP36L2 were detected in patients with KRAS mutations (Figure 5G), inclusion of ZFP36L2 mutations in patient stratification might increase the sensitivity of recurrence prediction with minimal loss of specificity ( Table 2).

| Clonal analysis
When tumors recurred after curative tumor resection, clinicians inferred their origin based on clinical information or pathological examination. Postoperative therapy is determined based on the tumor origin, i.e. whether cancer has emerged from the original tumor ("recurrence") or de novo ("double cancer"). Thus, it would be important to know whether discrimination based on clinical and pathological information is correct.
In total, 31 tumors from 12 patients were considered to be recurrent or peritoneally disseminated and sharing the same origin as the primary tumors, whereas 23 tumors from 10 patients were considered to be double cancers. The assessments were performed independent of this study in the course of routine clinical practice, considering the tumor site, time intervals, and likelihood of tumor recurrence from a pathological viewpoint.
Using clonal analysis of mutational information, we inferred that five cases of synchronous and two cases of   metachronous cancer developed de novo, whereas four cases of metachronous cancer occurred via local metastasis or peritoneal dissemination ( Figure 6). In all pairs of synchronous tumors, we found no shared mutations, including those in APC, TP53, or copy number alterations, indicating that the tumors may have initiated independently. By contrast, most metachronous lesions probably originated from a common ancestor as the initial sample and shared mutations in APC and TP53. In two cases, the clinical assessment was discordant with the clonal analysis. In the first case, the primary cancer (tumor number 1011) was a rectal adenocarcinoma, and rectal resection was performed as part of radical surgery. In the postoperative periodic examination, a new tumor was found on the oral side of the anastomosis of the former surgery (tumor number 1042). The tumor was not on the anastomosis site, and whether it was a local recurrence or double cancer was difficult to

| DISCUSSION
Precise prediction of postoperative recurrence would help improve treatment strategies for stage II/III CRC. In this study, genomic analysis was conducted to identify predictive biomarkers for the detection of stage II/III CRC recurrence and to explore genomic abnormalities. For selected cases, we performed WGS to explore genomic abnormalities in more detail, and several new findings were obtained.
Our results showed that LOF mutations in ZFP36L2 may be predictive of recurrence. ZFP36L2 encodes a zincfinger RNA-binding protein, which is reportedly involved in regulating mRNA stabilization 23 and protein synthesis. 23,24 Previous studies have shown that ZFP36L2 suppresses the proliferation of colon cancer cell lines in vitro 25 and that LOF mutations in ZFP36L2 are enriched in metastatic CRCs. 18,19 It was also reported that LOF of ZFP36L2 promotes tumor progression by impairing apoptosis due to DNA damage. 26 Therefore, ZFP36L2 could be a tumor suppressor for CRC, and its deficiency may predict recurrence or worse outcomes. However, the findings were inconsistent with previously reported data that revealed the expression of ZFP36L2 to be increased in 10% of gastric cancer, and ZFP36L2 to promote tumor proliferation. 27 In addition, the analysis of publicly available TCGA data was not consistent with that of our cohort. Such discrepancies may be attributed to differences in ethnicity, sequencing procedures (including library preparation), or mutational analysis methods, which affect the detection sensitivity for ZFP36L2 mutations. Indeed, mutations in ZFP36L2 were not described in the original report of the TCGA cohort. Moreover, we found gene expression associated with the immune reaction to be altered in tumors with ZFP36L2, which might explain the mechanistic basis of recurrence. However, as the difference in analytical procedures may affect the detection of ZFP36L2 mutations, a unified procedure along with a larger cohort would be required for elucidating the functions of ZFP36L2.
ZFP36L2 mutation status was associated with differences in RFS, but OS did not differ by ZFP36L2 mutation status; ZFP36L2 mutation status may predict tumor recurrence, but not response to second-line therapy. Further  analysis is needed to clarify the prognostic impact of ZFP36L2 mutation. Our data showed that chromosomal SVs often involve major oncogenes related to CRC. For instance, we found both alleles of APC to be affected by SV in one case and a relatively large deletion leading to loss of heterozygosity in other cases. Such aberrations are likely to escape target sequencing or WES detection. Therefore, a comprehensive analysis via WGS and long-read sequencing could help reveal genomic abnormalities and identify the underlying mechanism of action in CRC, such as mutations in ZFP36L2.
Clonal analysis was performed to examine the degree of concordance between clinical and genetic data. Our results showed that clinical assessment is precise, except for cases where the origin of tumors is difficult to infer based on clinical information alone. Clonal analysis based on sequencing data could contribute to improving clinical practice.
This study had the following limitations: (1) cases were selected based on the availability of tumor samples, and thus, comparisons between control groups were not possible; (2) our cohort was a mixture of patients with stage F I G U R E 6 Clonal analysis of metachronous and synchronous tumors. Metachronous tumors are shown on the left side, and synchronous tumors are shown on the right side. For individual cases, the anatomical locations are shown on the left, and the clonal evolutionary process is shown on the right. Cases in which tumors were derived from the same clone are marked by a blue box and those derived from independent clones by an orange box. Cases where the origin of tumors was difficult to infer clinically are marked by dottedline boxes. Numbers in circles indicate the number of mutations considered for clonal analysis. Tumor identification is indicated by fourdigit numbers. Anast, anastomosis; GL, germ line; Rec, recurrence II and stage III CRC, and the influence of confounding factors was unavoidable. However, we selected the nrCRC group such that the percentage of stage II tumors would be the same in the rCRC and nrCRC groups, which, we believe, justifies our analysis. Further investigation will be required to validate the significance of ZFP36L2 deficiency via prospective studies, especially for stage II CRC.
In conclusion, our study showed that detailed genomic analysis could improve the precise clinical characterization of CRC. In particular, ZFP36L2 was identified as a candidate biomarker for rCRC. We also believe that our data will provide a clue to introduce WGS into clinical genomics in the future.