A total of 224 fluoropyrimidine-treated patients (n = 120 in the “toxicity group”; n = 104 in the “no-toxicity group”) were sequenced by NGS. After quality control, 11 of the 224 sequenced samples were excluded from analysis because they did not pass the quality check of sequencing data by fastQC. All of these samples were characterized by a low number of produced paired-end reads (number of reads < 1,000,000) and a sequence length > 220 bp versus the expected 150–200 bp, indicating a failure in the sequencing step. An additional 35 samples were excluded from analysis based on the result of unsupervised cluster analysis, which classified the sample as “good” or “poor” based on coverage. Thus, the final “toxicity group” included 82 patients, whereas the “no-toxicity group” included 96 patients. Nine genes (i.e., UPRT, TYMP, RXRA, TK1, MIR27B, MIR27A, MIR23A, KDM6A, KDM6B) were excluded from the analysis because of unfavorable coverage (i.e., not meeting the criteria of median depth coverage of 50X for the genes defined by the panel in at least 95% of the samples).

Overall, the coverage analysis showed a median number of 542,162 mapped reads (range: 61,074–1563800) with a median percentage of reads at 1 × and 10x (depth of coverage) of 86% and 2%, respectively. Looking at coverage at the genomic level, the median depth at the gene level was homogeneous between the “toxicity” and “no-toxicity” groups, allowing comparison of the two cohorts.

The main clinical–demographic features of the two study populations are shown in Table 1. The two groups were balanced overall in terms of clinical–demographic characteristics, with the exception of treatment setting, as the adjuvant setting was more frequent among the controls while first-line or further lines of treatment were more frequent among the cases. Of the 82 patients in the “toxicity group,” 50 (61.0%) developed grade 3 toxicity and 32 (39.0%) developed grade 4 toxicity as the maximum grade of hematologic or non-hematologic toxicity experienced during the entire course of treatment. Forty-seven of 82 patients (57.3%) developed grade ≥ 3 hematologic toxicity, with neutropenia being the most common adverse event (35/47; 74.5%). Fifty-eight of 82 patients (70.7%) developed grade ≥ 3 non-hematologic toxicity, with diarrhea being the most common adverse event (27/58; 46.6%). Thirty-five of 82 patients (42.7%) experienced severe hematologic or non-hematologic toxicity within the first cycle of treatment (acute toxicity, cycle ≤ 1), and forty-seven patients (57.3%) after the first cycle.

A total of 7,420 and 7,896 germline variants were called against the reference genome in the “toxicity” and “no-toxicity” groups, respectively. The mean coverage (read depth) of the identified genetic variants was 97 (range: 50–530) and 75 (range: 50–214) for the “toxicity” and “no-toxicity” groups, respectively.

A total of 471 unique genetic variants were identified (426 single-nucleotide polymorphisms [SNP], 25 deletions [DEL], 17 insertions [INS], 3 double-nucleotide polymorphisms [DNP]) in the “toxicity group” (Fig. 2A). Of these, 275/471 (58.4%) were common (MAF ≥ 1%) and 196/471 (41.6%) were rare/very rare (MAF < 1%) or novel variants. Thirty-six of the 275 (13.1%) common variants and 45 of the 196 (23.0%) rare/very rare or novel variants were predicted as deleterious. Most of the identified variants were in the 3’ region (213/471, 45.2%). Of the remaining variants, 105/471 (22.3%) were synonymous, 91/471 (19.3%) were missense and 43/471 (9.1%) were in the 5’ region.

Pie chart visualizing the type of genetic germline variants identified in groups. A “toxicity” and B “no-toxicity” groups. “Other” variants include: De novo Start OutOfFrame, Frame Shift Del, Frame Shift Ins, In Frame Ins, and Start Codon SNP

In the “no-toxicity” group, 466 unique genetic variants were identified (414 SNP, 31 DEL, 17 INS, 3 DNP, 1 oligonucleotide polymorphism [ONP]). Of these 274/466 (58.8%) were common (MAF ≥ 1%) and 192/466 (41.2%) were rare/very rare (MAF < 1%) or novel variants. Thirty-five of the 274 (12.8%) common variants and 44 of the 192 (22.9%) rare/very rare or novel variants were classified as deleterious. Most of the identified variants were in the 3’ region (220/466, 47.2%). Of the remaining variants, 104/466 (22.3%) were synonymous variants, 85/466 (18.2%) were missense variants and 41/466 (8.8%) were in the 5’ region (Fig. 2B).

The frequencies of common and rare (including novel) variants by gene and subpathway in the “toxicity” and “no-toxicity” groups are shown in Figs. 3 and 4, respectively. As can be seen from the figures, the frequency distribution of the common variants in the two groups is similar at both at gene and subpathway level, while the frequency distribution of the rare variants in the two groups is remarkably different.

Frequency of common (MAF ≥ 1%) and rare (MAF < 1%, including novel) variants according to gene in patients with and without toxicity

Frequency of common (MAF ≥ 1%) and rare (MAF < 1%, including novel) variants according to subpathway in patients with and without toxicity

The role of rare/novel and common variants and pathways in the risk of developing severe fluoropyrimidine-related toxicity was assessed by SKATjoint analysis. The list of genes and related subpathways is reported in Additional file 1: Table S1. The SKAT statistics were calculated considering all variants, variants in the UTRs, and missense variants, grouping them also for the functional effect. The results of the SKATs analysis are shown in Table 2.

The DPYS (all variants, P = 0.024) and the PPARD (all variants deleterious, P = 0.039; 3’UTR/5’UTR variants deleterious, P = 0.022) genes were identified to significantly impact the risk of developing fluoropyrimidine-related toxicity by applying the Bonferroni correction. At subpathway level, “nuclear receptor” (all variants, p = 0.033) and “membrane transporter” (all variants, P = 0.037) resulted the classes of genes significantly associated with the risk of severe toxicity.

Based on the results of the SKAT analysis (Table 2), the role of the burden of rare and novel variants in the DPYS and PPARD genes on the risk of severe fluoropyrimidine-related toxicity was further investigated. The complete list of rare (MAF < 1%) and novel variants identified in the DPYS and PPARD genes in the “toxicity” and “no-toxicity” groups is shown in the Additional file 1: Tables S2 and S3, respectively.

As for the DPYS gene, it was significant in the “all variants” group: seven and four rare/very rare variants were detected in the case and control groups, respectively (Additional file 1: Table S2). The mean number of DPYS variants per 100 patients was significantly higher in the “toxicity” group that in the “no-toxicity” group (P = 0.047). The number of patients with at least one rare DPYS variant differed between groups (12.2% and 4.2% in the “toxicity” and “no-toxicity” groups, respectively), albeit with borderline significance (P = 0.055). Carrying at least one rare DPYS variant was associated with an approximately fourfold higher risk of severe cumulative toxicity (P = 0.030) after accounting for potential confounders in the multivariable model (Table 3). The same effect was detected for acute severe toxicity, although not significantly (P = 0.082). It was confirmed that the presence of at least one rare DPYD missense variant in the present study population increased the risk of cumulative severe toxicity, as previously published [2], although not significantly due to the small sample size (OR = 6.03; 95% CI 0.65–56.27; P = 0.115). After excluding five cases and one control with rare DPYD missense mutations from the analysis, the rare DPYS variants maintained their negative impact on the risk of severe cumulative toxicity (OR = 4.06; 95% CI 1.12–14.81; P = 0.115).

The PPARD gene was significantly associated with toxicity in both “all variants deleterious” and “3’UTR/5’UTR variants deleterious” groups; the two groups overlapped because all deleterious variants identified were in the UTR region (Additional file 1: Table S3). The mean number of PPARD variants per 100 patients was not significantly higher in the “toxicity” group than in the “no-toxicity” group (P = 0.472). The number of patients with at least one rare PPARD variant was not significantly different between groups (P = 0.595). Carrying at least one rare PPARD variant was not associated with the risk of experiencing severe cumulative (P = 0.759) or acute (P = 0.777) severe toxicity in multivariate analysis (Table 3).

The role of DPYS and PPARD polymorphisms (MAF > 1%) as predictors of severe toxicity was further investigated. The toxicity risk for each common DPYS (all variants) and PPARD (deleterious 3’UTR/5’UTR variants) variant identified by NGS analysis was compared between the “toxicity” and “no-toxicity” groups and is summarized in Table 4.

Two common polymorphisms were detected in the DPYS gene, the synonymous rs2298840 (c.216C > T; p.Phe72Phe) and the 3’UTR rs143004875 (g.chr8:105391734_105391735insT). Remarkably, these variants were identified only in the “toxicity” group, whereas they were not detected in the “no-toxicity” group. Although statistical analysis is hampered by the small number of cases, the rs143004875-T variant allele was associated with an increased risk of severe toxicity for both cumulative (P = 0.002) and acute toxicity (P = 0.005). A similar trend was observed for rs2298840, with all patients carrying the polymorphic A-allele developing severe toxicity.

As for the PPARD gene, all three identified deleterious variants are located in the UTR region: rs2016520 (5’UTR, g.chr6:35378778C > T), rs9658170 (3’UTR, g.chr6:35394504G > A), and rs9658167 (3’UTR, g.chr6:35394080G > A). Among these variants, rs2016520 was associated with the risk of developing severe cumulative (P = 0.001) and acute toxicity (P = 0.0001). Particularly, patients with the rs2016520-TT genotype had about threefold and sixfold higher risk of severe cumulative and acute toxicity, respectively.

To investigate in detail the functional role of the DPYS-rs2298840, DPYS-rs143004875, and PPARD-rs2016520 polymorphisms, bioinformatic analysis was performed using HaploReg v4.1 [24], RegulomeDB v2.0.3 [25], and Ensembl’s VEP-Ensembl GRCh37release release 110—July 2023 [22]. The methods and detailed results are summarized in Additional file 1: Table S4. In brief, the synonymous DPYS-rs2298840 variant could have a moderate impact on gene functionality and/or expression, as it broadly alters regulatory chromatin status and consensus sequence for transcription factors and exhibits eQTL hits. An impact on the splicing mechanism could also not be excluded. All these effects were summarized by a RegulomeDB rank score of 3a (i.e., transcription factors binding + any motif + DNase peak) and a probability score of 0.61235. The CADD score is 14.91. A similar functional prediction was obtained for the 3’UTR DPYS-rs143004875 polymorphism, although the effects appear to be smaller (RegulomeDB rank score = 6 [i.e., motif hit], probability score = 0.22365; CADD score not available). The 5’UTR PPARD-rs2016520 variant, that is in linkage (r2 ≥ 0.8) with 17 additional polymorphisms, could have a moderate functional consequence since it affects chromatin architecture and DNA accessibility for gene transcription, it is located in a transcriptional binding element with a resulting impact on the regulation of protein expression (eQTL hits) and regulates the splicing mechanism. All of these effects were summarized by a RegulomeDB rank score of 5 (i.e., transcription factors binding or DNase peak) and a probability score of 0.13454. The CADD score is 16.81.