The study was conducted on 125 patients aged between 17 and 50 from the Republic of Kazakhstan with early-onset of CRC. All the patients signed the informed consent and filled in the detailed questionnaire which included questions related to the FHC. The consent form and the questionnaire were approved by the IRB/IEC of Kazakh Institute of Oncology and Radiology (KazIOR) and Asfendiyarov Kazakh National Medical University (Almaty, Kazakhstan). In cases observed before 2017, cancer staging was done according to the 7th edition of the AJCC Cancer Staging System; in cases observed after 2017—according to the 8th edition of the AJCC Cancer Staging System (The AJCC Cancer Staging Manual, Eighth Edition, https://cancerstaging.org/references-tools/deskreferences/Pages/8EUpdates.aspx#).

Genomic DNA (gDNA) was isolated from EDTA treated peripheral blood samples using the GeneJET Genomic DNA Purification Kit (Thermo Scientific, USA). Qualitative and quantitative characteristics of the DNA samples were evaluated by the spectrophotometer (Eppendorf Biophotometer plus, Germany) and fluorometer (Quantus™ Promega, USA). At least 50 ng of double-stranded nuclear DNA was obtained per sample.

Next-generation sequencing (NGS) was performed using TruSight Rapid Capture Kit (Illumina, USA) in combination with the TruSight Cancer Sequencing Panel (Illumina, USA) according to the manufacturer's recommendations.

DNA libraries with 4 mM molarities were subjected to clustering using a standard flow cell and were sequenced on the MiSeq platform (Illumina) using the MiSeq Reagent Kit v3 (600 cycles).

Bioinformatics data analysis was carried out using two methodological approaches with different algorithms of analysis. In the first approach, the NGS data were analyzed using MiSeq Reporter v.3.0 software (Illumina). This software allowed mapping and aligning the sequenced sequences and comparing them with the reference sequence of the human genome (GRCH37/hg19) using the Burrows-Wheeler Aligner (BWA) algorithm. Search and detection of variants for specific regions of the genome were carried out using the Genome Analysis Toolkit (GATK, Broad Institute, Cambridge, USA).

In the second approach, various bioinformatics methods and software packages were used to optimize the workflow. The mapping and alignment of a sequence read against the human reference genome of the GRCH37/hg19 version were conducted using the Bowtie2 algorithm with a very sensitive local parameter. FastQC and MultiQC were used to evaluate the quality of reads of sequences. Picard Tools and SAM Tools/BCF tools were used to convert SAM files into BAM files, sort, and index the mapped reads, remove intermediate files, merge the BAM files, and identify duplicates. Java Runtime Environment and R Bioconductor were used to execute software scripts. GATK was used to re-align the mapped reads around the areas with insertions/deletions. GATK “Haplotype caller” strategy was used to filter and detect genomic variants. Finally, all VCF (variant call format) files containing only alternative variants obtained by GATK were imported into the Variant Studio 3.0 software (Illumina). In the study, only the variants identified by both approaches were considered.

Variant Studio was used to annotate and interpret genetic variants. A Q-score of 30, corresponding to an error rate of 0.1%, was selected as the acceptance threshold value. The following filtration parameters were set to exclude poor quality nucleotide variants: read depth >50×, alternative read depth >20×, quality value >100. Variants that matched the filtering parameters were selected for the study.

Genetic variants were annotated in accordance with the nomenclature of the Human Genome Variation Society (HGVS) (17); the interpretation was done using the Single Nucleotide Polymorphism Database (dbSNP, http://www.ncbi.nlm.nih.gov/projects/SNP/), ClinVar database, InSIGHT/LOVD database, and COSMIC. Integrative Genomics Viewer was applied to visualize the variants (18, 19).

Identified variants were categorized in accordance with the guidelines of the American College of Medical Genetics and Genomics (ACMG) (20) as pathogenic, likely pathogenic, uncertain significance (VUS), benign, and likely benign. Pathogenic mutations included those resulting in a premature stop codon (frameshift mutations and nonsense), variants with uncorrected splicing and variants affecting protein function according to the results of functional studies. Novel mutations such as frameshift mutations and nonsense mutations leading to protein truncation, as well as mutations without well-established functional studies, were classified as likely pathogenic. Synonymous and intronic variants not affecting splicing were regarded to be benign/probably benign, as well as the variants with a minor allele frequency (MAF) found in 5% or more of the general population according to the 1000 Genomes database (1000G) (21), the Exome Sequencing Project database (ESP6500, https://esp.gs.washington.edu) or the Exome Aggregation Consortium database (ExAC) (20). The remaining variants with no functional data and not fulfilling the criteria of classification as pathogenic/likely pathogenic or benign/likely benign, or with conflicting evidence of their benign and pathogenic nature, were defined as VUS.

In silico, bioinformatics tools such as Sorting Intolerant From Tolerant (SIFT) (22) and Polymorphism Phenotyping-2 (PolyPhen-2, http://genetics.bwh.harvard.edu/pph2/) (23) were used to predict potential pathogenic effects of missense variants on protein structure and function. The variants that matched the following criteria: a SIFT score of <0.05 and a PolyPhen-2 score between 0.95 and 1.0, were considered to be strongly suspected of being deleterious (24).

Demographics, clinical-pathological characteristics, and personal and family history variables were analyzed. Descriptive statistics were summarized as frequency distributions for categorical variables and as means for continuous variables. Wilcoxon rank sum tests were used to compare continuous variables, and the chi-square test was performed to compare categorical variables. The results were considered statistically significant at a p < 0.05. All statistical tests were two-sided at the 5% level of significance. Data were processed and analyzed by using R software, version 3.6.0 (W.N. Venables, D.M. Smith, and the R Core Team).

Sequencing data have been deposited at the NCBI SRA archive with BioProject record PRJNA548114 and SRA accession SUB5747000.

Study cohort included 62 (49.6%) men and 63 (50.4%) women with early onset CRC. Th ethnic composition of the patient group included Kazakhs−61 (48.8%), Russians−36 (28.8%), Uyghurs−7 (5.6%), Azerbajanians−3 (2.4%), Koreans−3 (2.4%), Germans−3 (2.4 %), Ukrainians−2 (1.6%), Tatars−2 (1.6%), Belarusians, Dungans, Turks−1 patient each (0.8% each), patients with mixed ethnic background−5 (4%).

The median age of patients at CRC diagnosis was 40.7 (39.6 in men, 41.7 in women). 25 patients (20.1% of the early CRC onset cohort, of them, 15 men and 10 women) were diagnosed with CRC before the age of 35 (median age 29.08 ± 4.75 for both sexes; 28.4 ± 45.2 in men and 30.1 ± 5.15 in women).

In terms of clinical pathology, patients with early CRC onset differed by TNM staging, tumor localization, and histological types. The pre-dominant histological type was moderately differentiated adenocarcinoma detected in 81 (64.8%) cases. Half of the patients (57.2%) were diagnosed with stage III cancer. The rectum was a more common tumor site in both men and women compared to other sites. No significant statistical difference was found in all clinical characteristics between males and females (Table 1).

Of all the subjects, 19 (15.2%) patients with a family history of cancer (FHC), 6 had the CRC syndrome (Table 2). The characterization of all the patients with early onset CRC is shown in Figure 1.

A total of 125 unrelated patients with early onset CRC were analyzed by NGS using the TruSight Cancer Sequencing Panel which targeted a set of 94 genes known to play a role in cancer pre-disposition. Bioinformatics analysis of NGS data has revealed 11,152 variants (Supplementary Table 1) in 85 genes, of them, 3,790 missense, 6,254 synonymous variants, 44 3′UTR variants, 10 frameshift variants, five stop-gain variants, four in-frame deletions, two splice donors, one splice acceptor variant, and 1,042 intron or non-coding variants. APC, BRCA2/1, ALK, BRIP1, EGFR, FANCA, FANCD2, FANCI, HNF1A, MEN1, NSD1, PMS2, RECQL4, RET, SLX4, WRN, and XPC genes mutated most often.

The approach to the interpretation of NGS data described in the corresponding databases and guidelines mentioned in the Materials and Methods section has revealed 24 pathogenic/likely pathogenic mutations. 472 variants were classified as VUS; 10,656 variants—as benign/likely benign.

Twenty-four pathogenic and probable pathogenic mutations were detected in 20 patients (16% of the entire group of patients) with cancer of the cecum, ascending colon, hepatic flexure of the colon, splenic flexure of the colon, sigmoid colon, rectosigmoid colon, and rectum. Of them, five patients had a family history of FAP, CRC, gastric cancer, or ovarian cancer. Six patients (four cases of synchronous cancers and two cases of metachronous cancer) had primary multiple tumors (PMTs). Out of those pathogenic and likely pathogenic variants, five were detected in the APC gene, three in FANCI, three in BRCA2, two in BRCA1, two in MLH1, and one mutation each, in MSH6, MUTYH, BLM, NBN, ATM, BMPR1A, CHEK2, AIP, and DICER1 genes. As shown in Table 3, the analysis of mutation type has revealed 10 frameshift mutations, five missense mutations, five stop-gain mutations, one in-frame deletion, and three mutations involving uncorrected splicing. All detected pathogenic and likely pathogenic mutations were in the heterozygous state. In total, 23 mutations were unique, eight of them represented novel variants. Those new mutations have not been previously mentioned in the LOVD and ClinVar databases and have not been described in publications. All mutations were rare, and for 17 of them, frequency data was not available in the 1000G, ESP6500 and ExAC databases (Table 3). 62.5% of the mutations were available in the dbSNP database, 29.1% in the COSMIC catalog, and 38.4% in the ClinVar and/or LOVD databases. Multiple alterations of genes were noticed in three patients. Three frameshift variants including two novel variants (c.907delA in APC gene and c.6743dupA in ATM) and one previously reported mutation (c.657_661delACAAA in NBN) were detected in one patient (male, Kazakh) at the age of 40 with rectum cancer. Two more patients (one male at the age of 26 with FAP and sigmoid colon cancer and one male at the age of 46 with ascending colon cancer and rectum cancer) had one novel frameshift variant each (FANCI c.3340delA and BMPR1A c.152delC, respectively) in combination with other previously reported frameshift variants (APC c.3613delA and CHEK2 c.599T > C, respectively, Table 3).

Pathogenic mutations were most common in the APC gene that encoded tumor suppressor protein acting as an antagonist of the Wnt signaling pathway. APC gene mutations were strongly associated with CRC carcinogenesis (12, 13). Most of them were found in CRC patients with polyposis, FAP, and PJS that could be considered as disease-causing mutations. Pathogenic mutations were also common in Fanconi Anemia associated gene (FANCI, 12.5%) known to play a significant role in repairing double-stranded DNA breaks by homologous recombination and restoring inter-chain DNA crosslinks (26). Fanconi Anemia genes have been reported to play an important role in CRC inherited pre-disposition (27). The mutations of one of the mismatch repair (MMR) genes, MLH1 associated with Lynch syndrome (28, 29) were also quite frequent (8.3%).

Mutations were also found in non-CRC susceptibility genes such as BRCA1 (8.3%) and BRCA2 (12.5%). Those tumor suppressor genes controlling genome integrity were strongly associated with breast and ovarian cancer (30, 31). Pearlman et al. have shown a high prevalence of mutations in the non-CRC genes (BRCA1/2−8.3%) among patients with early onset of CRC (32). However, some studies have demonstrated a modest link between CRC risk and mutations in BRCA1/2 (33, 34). Notably, 1 female (Patient ID: CRC629) with a BRCA1 mutation (c.5329dupC (also known as c.5382insC)/p.Gln1777ProfsTer74, 0.8%) with primary multiple cancer (CRC, ovarian, and breast cancer) had a family history of ovarian cancer. In the studies of CRC patients from different populations, BRCA1 c.5382insC mutation was found in unselected Jewish Ashkenazi patients (0.44%) (25, 35), in unselected Polish patients (0.29%) and in patients with a reported family history of CRC (0.93%) (34).

Of all the 472 identified VUS, 131 were potentially deleterious according to at least one in silico prediction tool (PolyPhen-2 and SIFT) (Supplementary Table 2). Forty-five of all VUS were selected according to the classification shown in Figure 2 and they were categorized as being strongly deleterious (Supplementary Table 3). Of these 45 VUS (42 missense variants in 35 patients and three in-frame deletions in three patients) were identified. Three patients (Patient ID: CRC530, CRC597, and CRC625) had in-frame deletions in the APC, FANCE, and TSC2 gene, respectively. The largest number of VUSs were found in MSH2 involved in mismatch repair and BRCA2 which is a tumor suppressor gene involved in the maintenance of genome stability, specifically the homologous recombination pathway for double-strand DNA repair (Figure 3).

Eight novel alterations were identified. Two of the novel alterations, c.7429G>A/p.Gly2477Arg in the ATM gene and c.1306C > A/p.Leu436Met in the FANCD2 gene were identified in patients with FHC (Patient ID: CRC335 and CRC587). Other six novel variants were detected in patients with no FHC, c.1135G > A/p.Ala379Thr in NSD1 gene, c.2777A > G/p.Glu926Gly in RB1 gene, c.2693G > A/p.Arg898Lys in the BLM gene, c.1031A > C/p.Gln344Pro in MSH2 gene, c.1493T > G p.Phe498Cys in DICER1 gene and c.278G > A/p.Arg93His in PMS1 gene in a patient (Patient ID: CRC380, CRC442, CRC599, CRC600, and CRC635).

The three novel variants in patients (Patient ID: CRC599, CRC600, CRC635) were found in high and moderate penetrance genes (BLM, MSH2, and PMS1) which are known to be associated with CRC so these mutations can be regarded as one of the risk factors for developing CRC.

The patient cohort under study was conditionally divided into three subgroups by their clinical and molecular characteristics: patients with FHC, patients with PMTs, and patients with sporadic cases (SC) of the disease had the average age 38.9 ± 9.1, 44.6 ± 3.3, and 40.7 ± 7.3, respectively. Detailed information about the most significant mutations with deleterious effects found in patients with FHC (n = 19), PMTs (n = 10), and SCs (n = 98) are represented in Supplementary Tables 4–6, respectively. For patients with FHC, PMTs and SC, 32, 9, and 93 variants were classified as VUS with potentially damaging effects according to the criteria of PolyPhen and SIFT databases and six, five, and eleven variants identified were considered as pathogenic in accordance with the ACMG guidelines, respectively. Early-onset CRC patients with pathogenic variants and VUS were evaluated according to the level of cancer pre-disposition genes (high and moderate penetrance genes listed in the National Comprehensive Cancer Network (NCCN) (36, 37), which are shown in Figures 4, 5, respectively. Pathogenic mutations in high penetrance CRC genes were higher in patients with FHC (21.1 vs. 3.1%; P = 0.0002) and in patients with PMTs (20.0 vs. 3.1%; P = 0.0004) compared with patients without/unknown FHC (3.1%) (Figure 4).

Segregation analysis was done only for two patients with FAP and CRC with the information of relatives (Patient ID: CRC613: sister and daughter, Patient ID: CRC2 brother) but not for the whole cohort because of the unavailability of data (Figure 6). These patients had three frameshifts and two stop codon mutations in APC, MLH1, and FANCI genes. One of those mutations in the APC gene c.3613delA (p.Ser1205AlafsTer60) was previously described in the Polish population among family members with FAP (38). Other mutations [one in FANCI c.3340delA (p.Thr1114ProfsTer50), two in APC c.4128T > G (p.Tyr1376Ter) and c.2835delG (p.Thr946HisfsTer9), and one in MLH1gene c.436C > T (p.Gln146Ter)] have not been previously described in the literature. However, APC c.4128T > G somatic mutation has been reported in the COSMIC database.

DNA sequencing has revealed identical genetic defects in the proband's brother aged 28: germline hereditary mutations [one missense (c.2146G > A) in the MLH1 gene and one pathogenic frameshift mutation (c.3613delA) in the APC gene] and daughter of a FAP (CRC613) proband with primary multiple CRC tumors had identical APC 4128T > G/p.Tyr1376Ter pathogenic stop-gain and FANCM c.4881T > G missense mutations, while the same proband's sister had another APC c.2835delG/p.Thr946HisfsTer9 pathogenic frameshift mutation (novel) and a MLH1 c.436C > T/p.Gln146Ter pathogenic stop-gain mutation (described in ClinVar and LOVD) (Table 3).

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.