An integrative analysis reveals cancer risk associated with artificial sweeteners

Abstract

Background: Artificial sweeteners (AS) have been widely utilized in the food, beverage, and pharmaceutical industries for decades. While numerous publications have suggested a potential link between AS and diseases, particularly cancer, controversy still surrounds this issue. This study aims to investigate the association between AS consumption and cancer risk.

Methods: Targets associated with commonly used AS were screened and validated using databases such as CTD, STITCH, Super-PRED, Swiss Target Prediction, SEA, PharmMapper, and GalaxySagittarius. Cancer-related targets were sourced from GeneCards, OMIM, and TTD databases. AS-cancer targets were identified through the intersection of these datasets. A network visualization ('AS-targets-cancer') was constructed using Cytoscape 3.9.0. Protein-protein interaction analysis was conducted using the STRING database to identify significant AS-cancer targets. GO and KEGG enrichment analyses were performed using the DAVID database. Core targets were identified from significant targets and genes involved in the 'Pathways in cancer' (map05200). Molecular docking and dynamics simulations were employed to verify interactions between AS and target proteins. Pan-cancer and univariate Cox regression analyses of core targets across 33 cancer types were conducted using GEPIA 2 and SangerBox, respectively. Gene chip datasets (GSE53757 for KIRC, GSE21354 for LGG, GSE42568 for BRCA, and GSE46602 for PRAD) were retrieved from the GEO database, while transcriptome and overall survival data were obtained from TCGA. Data normalization and identification of differentially expressed genes (DEGs) were performed on these datasets using R (version 4.3.2). Gene Set Enrichment Analysis (GSEA) was employed to identify critical pathways in the gene expression profiles between normal and cancer groups. A cancer risk prognostic model was constructed for key targets to further elucidate their significance in cancer initiation and progression. Finally, the HPA database was utilized to investigate variations in the expression of key AS-cancer target proteins across KIRC, LGG, BRCA, PRAD, and normal tissues.

Results: Seven commonly used AS (Aspartame, Acesulfame, Sucralose, NHDC, Cyclamate, Neotame, and Saccharin) were selected for study. A total of 368 AS-cancer intersection targets were identified, with 48 notable AS-cancer targets, including TP53, EGFR, SRC, PIK3R1, and EP300, retrieved. GO biological process analysis indicated that these targets are involved in the regulation of apoptosis, gene expression, and cell proliferation. Thirty-five core targets were identified from the intersection of the 48 significant AS-cancer targets and genes in the 'Pathways in cancer' (map05200). KEGG enrichment analysis of these core targets revealed associations with several cancer types and the PI3K-Akt signaling pathway. Molecular docking and dynamics simulations confirmed interactions between AS and these core targets. HSP90AA1 was found to be highly expressed across the 33 cancer types, while EGF showed the opposite trend. Univariate Cox regression analysis demonstrated strong associations of core targets with KIRC, LGG, BRCA, and PRAD. DEGs of AS-cancer core targets across these four cancers were analyzed. GSEA revealed upregulated and downregulated pathways enriched in KIRC, LGG, BRCA, and PRAD. Cancer risk prognostic models were constructed to elucidate the significant roles of key targets in cancer initiation and progression. Finally, the HPA database confirmed the crucial function of these targets in KIRC, LGG, BRCA, and PRAD.

Conclusion: This study integrated data mining, machine learning, network toxicology, molecular docking, molecular dynamics simulations, and clinical sample analysis to demonstrate that AS increases the risk of kidney cancer, low-grade glioma, breast cancer, and prostate cancer through multiple targets and signaling pathways. This paper provides a valuable reference for the safety assessment and cancer risk evaluation of food additives. It urges food safety regulatory agencies to strengthen oversight and encourages the public to reduce consumption of foods and beverages containing artificial sweeteners and other additives.

Keywords: Artificial sweetener; Cancer; Core targets; Machine learning; Molecular dynamics simulation; Network toxicology; Pathways in cancer.

Fig. 1

Flowchart of this study

Fig. 2

Artificial sweetener-cancer targets obtain. A Daily consumption of foods containing artificial sweeteners (AS). B Molecular structures of seven commonly used AS. C AS targets acquisition database. D AS targets Venn diagram. E Cancer targets Venn diagram. F AS-cancer targets Venn diagram

Fig. 3

AS-cancer core targets acquisition. A AS-target-cancer network. The pink octagon symbolizes AS, the green triangle denotes seven commonly used types of artificial sweeteners, the orange diamond represents cancer, the purple circle signifies 35 core targets, and the blue lines primarily indicate the relationships between the seven AS and the core targets. B PPI network diagram and visualization of AS-cancer core targets. In the PPI visualization diagram, core targets are sorted based on their degree value. A darker color and larger circle area indicate a higher Degree value, signifying stronger interactions among targets

Fig. 4

GO and KEGG enrichment analysis of AS-cancer core targets. A Triadic histogram depicting the GO biological processes (BP), molecular functions (MF), and cellular components (CC) of AS-cancer core targets. In this visualization, green denotes the BP terms, orange represents the CC terms, and purple signifies the MF terms. Each column reflects the GeneRatio associated with its respective terms. B KEGG enrichment analysis of top 10 pathways for AS-cancer core targets. Different colors on the right half of the ring represent various pathways. On the left half, different colors indicate PValue for the 35 core targets. The lines within the circle illustrate relationships between pathways and targets

Fig. 5

Molecular docking between AS and AS-cancer core targets. A The molecular docking heat map. The results display seven different circles in various colors, each representing a different AS molecule. The size of each circle correlates with the value of the binding energy. B–I 2D and 3D interacted structures of AS and AS-cancer core targets. CCND1-Aspartame, MAPK1-Acesulfame, CASP3-Sucralose, MAPK1-NHDC, CDK4-Cyclamate, KRAS-Neotame, KRAS-Saccharin, and BCL2L1-NHDC

Fig. 6

Molecular dynamics simulation. A and D Root mean square deviation (RMSD) of MAPK1-NHDC and CDK4-cyclamate. B and E Root mean square fluctuation (RMSF). C and F The radius of gyration (RG) and its values along the three axes (Rgx, Rgy, Rgz)

Fig. 7

AS-cancer core targets pan-cancer analysis expression difference heat map. The expression patterns of 35 AS-cancer core targets across 33 types of cancer are visualized. Red indicates higher expression levels, while blue indicates lower expression levels

Fig. 8

Univariate Cox regression analysis of AS-cancer core targets. A Visualization heat maps display the differential expression of 35 AS-cancer core targets across 33 types of cancer. In these maps, ‘ns’ denotes ‘0’, ‘*’ denotes ‘1’, ‘**’ denotes ‘2’, ‘***’ denotes ‘3’, and ‘****’ denotes ‘4’. B Forest plots of univariate Cox regression for CDK4 across 33 cancers show that an HR = 1 indicates no effect, HR < 1 indicates a favorable effect, and HR > 1 indicates an adverse effect. CDK4 demonstrates significant differences in LGG, LIHC, MESO, KIRP, KICH, SKCM, and LAML. C Forest maps of univariate Cox regression for CASP3 across 33 types of cancer reveal significant differences in LGG, ACC, SKCM, GBM, LIHC, and UVM

Fig. 9

Differentially expressed genes (DEGs) screening, target-pathway diagram and GSEA pathway enrichment analysis in clear cell renal cell carcinoma (GSE53757) and brain low-grade glioma (GSE21354). A Volcano plots depicting DEGs in a renal cancer dataset (GSE53757) reveal significant up-regulation of CCND1, MAPK1, MYC, CDKN1A, and PDGFRB genes in KIRC patients, while the EGF gene shows significant down-regulation in these patients. B Target-pathway diagrams in AS and KIRC. In these diagrams, red indicates genes significantly upregulated, blue denotes significantly down-regulated genes, and green represents genes with insignificant changes. Shapes correspond to different targets: triangles denote sweetener targets, squares denote disease targets, and circles denote targets common to both sweeteners and diseases. C The top 5 GSEA enrichment pathways up-regulated in KIRC. D The top 5 GSEA enrichment pathways down-regulated in KIRC. E Volcano plot of DEGs in a low-grade glioma dataset (GSE21354) shows significant up-regulation of RELA, BCL2L1, MYC, TP53, IGF1A, CCND1, and CDK4 in patients with LGG. F Target-pathway diagrams in AS and LGG. G The top 5 GSEA enrichment pathways up-regulated in LGG. H The top 5 GSEA enrichment pathways down-regulated in LGG

Fig. 10

Renal cancer prognosis targets screening, and a cancer risk prognostic model constructing. A Univariate Cox regression analysis initially identified 20 prognostic genes (P < 0.05). B A coefficient profile plot was generated against the log(λ) values in the LASSO regression analysis model. The optimal parameter (λ) was marked with a dashed line on the left. C LASSO regression analysis of the 20 prognostic genes. D Multivariate Cox analysis identified 5 prognostic genes. E Kaplan–Meier survival analysis in high and low risk groups. F ROC curves of overall survival at 1, 3 and 5 years of the prognostic model. G The relationship between the survival status of renal cancer patients and the risk score of the prognostic model, along with changes in gene expression levels as the risk score increases. The dashed line delineates the boundary between the low-risk and high-risk groups

Fig. 11

Brain low-grade glioma prognosis targets screening, and a glioma cancer risk prognostic model constructing. A Univariate Cox regression analysis initially identified 19 prognostic genes (P < 0.05). B A coefficient profile plot was generated against the log(λ) values in the LASSO regression analysis model. The optimal parameter (λ) was marked with a dashed line on the left. C LASSO regression analysis of the 19 prognostic genes. D Multivariate Cox analysis identified 8 prognostic genes. E Kaplan–Meier survival analysis in high and low risk groups. F ROC curves of overall survival at 1, 3 and 5 years of the prognostic model. G The relationship between the survival status of brain low-grade glioma patients and the risk score of the prognostic model, along with changes in gene expression levels as the risk score increases. The dashed line delineates the boundary between the low-risk and high-risk groups

Fig. 12

Validation of the expression of core target proteins in KIRC and LGG, as well as in their respective normal tissues. A–C Immunohistochemical staining of CDKN1A, ERBB2, and BCL2 proteins in KIRC and normal renal tissues. D–F Immunohistochemical staining of CDK4, JUN, and BCL2L1 in LGG and normal brain tissue