From January 2020 to June 2021, a total of 85 colorectal cancer patients with or without organ metastasis who underwent surgical resections in the Jiangsu Cancer Hospital were included in this study. Surgical samples of the primary lesion were collected for genomic sequencing tests. Matched peripheral blood samples were also profiled as germline controls. The clinicopathologic features of each patient were retrospectively reviewed. All patients provided informed consents for taking part in this study.

The tissue and blood samples were tested in a CLIA- and CAP-accredited central laboratory (Nanjing Geneseeq Technology Inc., Nanjing, China). In brief, genomic DNA (gDNA) was extracted from formalin-fixed paraffin-embedded tissue samples by using QIAamp DNA FFPE Tissue Kit (Qiagen) according to the manufacturer’s protocol. Peripheral blood was centrifuged at 1,900g for 10 min and gDNA from the white blood cells was also extracted as negative germline control. DNA quality was examined using the NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific Waltham, MA, USA), and DNA quantification was performed on the Qubit 3.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA) with the dsDNA HS Assay Kit (Life Technologies, Carlsbad, CA, USA). The purified DNA was fragmented, and DNA libraries were constructed with the KAPA Hyper Prep kit (KAPA Biosystems) according to the manufacturer’s instruction. A panel targeting 425 cancer related genes (GeneseeqPrime™) using customized xGen lockdown probes (Integrated DNA Technologies) was used for hybridization capture. The hybridization reaction was performed with Dynabeads M-279 (Life Technologies) and xGen lockdown hybridization and wash kit (Integrated DNA Technologies). The target-enriched libraries were sequenced on the Illumina Hiseq4000 NGS platform (Illumina, San Diego, CA, USA) using the 150 bp paired-end reading according to the manufacturer’s protocol. The average coverage depth for the tumor tissue and the white blood cells were 1,000× and 60× respectively.

The sequence data analysis was carried out as previously described (24).In brief, the adapter reads, low quality reads and N bases were trimmed from the raw FASTQ files using Trimmomatic (25). Paired-end reads were aligned to the reference human genome (hg19). PCR deduplication was performed with Picard (http://broadinstitute.github.io/picard/). The base quality score recalibration and indel realignment were performed using the Genome Analysis Tool Kit (25). Somatic single nucleotide variants (SNVs) and small insertion/deletions (indels) were called by VarScan2 (26). Sequencing results from the tissue samples were compared with the matched white blood cell to identify the somatic mutations. SNVs and INDELs were further filtered using the following criteria: i) minimum ≥5 variant supporting reads and variant allele fraction (VAF) of ≥1%, ii) filtered out if present in >1% population frequency in the 1000 Genomes Project or ExAC database, iii) filtered out through an internal database of recurrent sequencing errors (≥3 variant reads and ≤20% VAF in at least 30 out of 2,000 normal control samples) on the same sequencing platform. Copy number variation (CNV) was screened using CNV kit. A fold change of more than 1.50 was considered to be a copy number gain and less than 0.65 was defined as copy number loss. Tumor mutational burden (TMB) was defined as the total number of base substitutions and indels in the coding region of targeted genes, namely, synonymous variants to reduce sampling noise and excluding known driver mutations as they are over-represented in the Panel, as previously described (27, 28).

MSI analysis was performed as previously described (29). The stability of a total of 52 microsatellite sites covered by the 425-cancer gene panel with a minimum of 15 bp repeats, including the classic MSI sites BAT-25, BAT-26, NR-21, NR-24, and MONO-27, was estimated and compiled into an overall MSI score. A site is considered qualified if it is covered by at least 101× depth of coverage. A sample is identified as MSI if more than 40% of the qualified sites with at least 100× coverage displayed instability.

Comparisons of proportion between groups were carried out using the Fisher’s exact test. Correlation analyses among gene alterations in different group were calculated by Spearman’s rank test. A two-sided P-value of <0.05 was considered statistically significant. All statistical analyses were performed using R (v.3.4.1).

Baseline clinicopathological features of the 85 CRC patients are summarized in Table 1 . The median age of the study cohort was 55 years (range 29-78) and 50% were males. All patients were ECOG PS 0-1. Of these, 34 (40%) of patients had their primary lesions localized to the left colon, 23 (27.1%) right colon, and 12 (22.2%) rectum. The majority (77/85, 90.4%) of the cohort had adenocarcinoma, followed by the mixed histological type of six (7.1%) patients. Based on the status of metastasis, patients were subdivided into non-organ metastasis (Non-OM, n = 59) and organ-metastasis (OM, n = 26) groups. Notably, two (2.4%) cases were classified as signet-ring cell carcinoma. Both were stage IV CRC and belonged to the OM group, with metastasis to the very uncommon organ peritoneum plus bladder and peritoneum plus bone, respectively. Within the OM group, 20 (76.9%) were with liver metastasis, followed by 17 (65.4%) peritoneum metastasis, five (19.2%) lung metastasis and three other organ metastasis. Moreover, there are ten (38.5%), 13 (50%), and three (11.5%) cases possessing one, two, and three organ metastatic sites, respectively.

A total of 1,981 somatic mutations in 334 genes were detected. The top frequently altered genes in this cohort included TP53 (82.4%), APC (69.4%), KRAS (55.3%), SMAD4 (25.9%), PIK3CA (22.4%), and FBXW7 (21.2%) ( Figure 1 ). The majority of variant types were missense mutations (1,427/1,981, 72.0%), which was followed by nonsense mutations (236/1,981, 11.9%) and frame-shift mutations (205/1,981, 10.3%; Figure 2 ). The two most common mutated genes in CRC were analyzed. A total of 80 TP53 alterations in 70 patients were detected ( Supplementary Figure S1A ). The p.R175H alteration was most common and 90% (72/80) of the alterations in TP53 were oncogenic or likely oncogenic according to OncoKB annotation. Mutations in LRP1B were found in 12 (14.1%) patients, and were mutually exclusive with TP53 alterations (P <0.01; Figure 3 ). Meanwhile, we detected 105 APC alterations belonging to 59 patients ( Supplementary Figure S1B ) and 88.6% (93/105) of the mutations were oncogenic and likely oncogenic.

Next, we examined the frequencies of oncogenic driver mutations in our cohort. Oncogenic mutations in KRAS, NRAS, and BRAF were detected in 55.3% (47/85), 7.1% (6/85), and 11.8% (10/85) of the cohort. KRAS G12 (82%, 41/50) was the most commonly mutated codon, and included G12D (30%, 15/50), G12V (16%, 8/50), G12R (8%, 4/85), G12C (6%, 3/50), G12A (4%, 2/50), and G12S (4%, 2/50) variants. The second most frequently altered codon was G13, consisting of G13D (14%, 7/50) and G13R (2%, 1/50). Furthermore, we also detected another eight (16%) types of KRAS variants located in the 2nd, 3rd, and 4th exons ( Supplementary Figure S1C ).

We also detected relatively high frequencies of alterations in EGFR (4.7%, 4/85), ALK (8.2%, 7/85), MET (2.4%, 2/85), ROS1 (5.9%, 5/85), ERBB2 (7.1%, 6/85), and NTRK1/2/3 (10.1%, 10/85), which may potentially serve as targets for targeted therapies. None of these mutations were known oncogenic/targetable mutations, except for an irregular SPTBN1~ALK fusion gene that retained the complete kinase domain of ALK. The presence of ALK fusion suggested that CRC might benefit from ALK-targeted therapies.

In addition, we performed mutational analysis comparing patients with various clinical characteristics. As compared to T1-3, T4 tumors had significantly higher mutational frequency in three genes, BRAF (5/15 vs. 5/69, p = 0.01), IFNE (2/15 vs. 0/69, p = 0.03), and RNF43 (5/69 vs. 7/69, p = 0.03). By comparing the patients with and without lymph metastasis, we found that alteration frequencies of 79 genes were significantly different between two groups, of those six genes had mutational frequencies more than 10% in this cohort ( Supplementary Table 1 ). Specifically, mutations of five genes favored non-lymph metastasis, namely ATM (8/28 vs. 2/56, p <0.01), CHD4 (7/28 vs. 2/56, p <0.01), AMER1 (9/28 vs. 5/56, p = 0.01), KMT2A (7/28 vs. 3/56, p = 0.01), and NOTCH2 (7/28 vs. 3/56, p = 0.01). On the other hand, TP53 mutations were enriched in the lymph metastasis group (19/28 vs. 51/56, p = 0.01). Moreover, a comparison between patients older or younger than median age was performed and showed that mutations in KDR (7/44 vs. 0/41, p = 0.01) were more common in the younger patients while mutations in ARAF (0/44 vs. 5/41, p = 0.02), AXL (0/44 vs. 5/41, p = 0.02) and APC (26/44 vs. 33/41, p = 0.04) were more common in the older patients. These results could potentially help to reveal the pathogenesis of CRC patients with diverse clinical features.

To identify potential genetic risk factors of organ metastasis, we analyzed the genetic features between the Non-OM and OM subgroups. To examine if the results were reliable, a TCGA CRC dataset (23) was used as a validation cohort.

To assess the influences of single gene aberrations on organ metastasis, we compared the genetic mutations (with mutational frequencies >10%) between the two groups. We observed a higher frequency of FBXW7 alterations in Non-OM than that in the OM group (27.1% vs. 7.7%; p = 0.05; Figure 1 ), which was validated in the TCGA CRC dataset (Non-OM, 19.3%; OM, 0%; p <0.01; Supplementary Figure S2A ). To further examine if FBXW7 mutations act independently of other genes, we performed a mutual exclusivity or co-occurrence analysis. Notably, the FBXW7 mutations occurred independently from mutations in other genes ( Figure 3 ). Functional annotations of all the FBXW7 variants with OncoKB revealed that the oncogenic or likely oncogenic FBXW7 mutations were significantly enriched in the Non-OM group than that in the OM group (20.3% vs. 3.8%; p = 0.05; Supplementary Figures S2B, C ; Supplementary Table 2 ). We also evaluated the functional variants in various oncogenic signaling pathways that might be related to CRC organ metastasis. Interestingly, the NOTCH pathway locating downstream of FBXW7 was more frequently altered in the Non-OM than in the OM group in both our cohort (22% vs. 3.8%, p = 0.05) and the TCGA cohort (10.9% vs. 0, p = 0.05). Thus, we speculated that FBXW7 and its downstream NOTCH pathway might play an independent and suppressive role in CRC organ metastasis.

In addition, we studied the distribution of all the altered genes in both our and the TCGA cohorts and found that there are 54 genes (all with mutational frequencies less than 10%) only mutated in the Non-OM groups of both cohorts but not in the OM groups. A pathway analysis implicated that these mutations were significantly enriched in several important cancer-related pathways, including the p38 MAPK cascade, response to cAMP, and leukocyte activation, etc. ( Supplementary Figure S3 ). Thus, these pathways might play vital roles in the regulation/suppression of CRC organ metastases.

Apart from the feature of single genes and pathways, distinct mutational features were identified among different groups. We analyzed the mutational type proportion between the two groups. The results suggested that in both groups, the top common mutational type were missense, frame-shift and nonsense variants, with minor trends of differences between the Non-OM subgroup [68.4% (793/1,160), 13.7% (159/1,160), and 11.7% (136/1,160) respectively] and the OM subgroup [77.2% (634/821), 5.6% (46/821), and 12.2% (100/821) ( Supplementary Figure S4A ). Similar results were observed in the TCGA cohort ( Supplementary Figure S4B ). Moreover, the numbers of the somatic mutations of each patient between the groups were also investigated, though no significant difference was observed in either our cohort or the validation cohort ( Supplementary Figures S5A, B ).

Since there are evidences that microsatellite (MS) instable CRC and tumor mutation burden (TMB) high CRC may possess disparate biological characters (30, 31), we subsequently inspected the associations of TMB and MS status with respect to the organ metastasis status ( Table 1 ). In our cohort, 22 patients had high TMB (TMB ≥10 mut/Mb), of which seven patients were MSI-H ( Supplementary Figure S6 ). Given the moderate sample size, both the TMB-H (n = 22, 25.9%, p = 0.43) and MSI-H (n = 7, 8.2%, p = 0.43) patients showed no statistically significant associations with OM vs. Non-OM compared with the TMB-L and MSS patients, respectively ( Supplementary Figures S7A, B ). In FBXW7-mutated tumors, there was a trend towards higher TMB compared with those with wild-type FBXW7 (p = 0.12; Supplementary Figure S8 ). In line with the previous result, within the TMB-L population, FBXW7 mutations more commonly occurred in the Non-OM group than that of OM group (26.2% vs. 4.8%, p = 0.05). Interestingly, within the TMB-H population, it was CTNNB1, but not FBXW7, mutations that significantly differed between the non-OM and OM groups, with the frequency of 23.5% (4/17) in the Non-OM and 80% (4/5) in the OM (p = 0.04) groups. This finding suggested that TMB-H CRC and TMB-L CRC might be driven by distinct organ metastatic mechanisms. Due to the low number of the MSI-H cases, we could not perform the similar analysis.