Machine learning prioritization of antibiotic residues in aquatic foods reveals exposure-driven genotoxic risk mediated by BCL2

TL;DR

This study uses machine learning to identify antibiotic residues in aquatic foods that pose genotoxic risks, showing how exposure levels drive these risks and how BCL2 can mitigate them.

Contribution

A novel ML framework that integrates surveillance data and mechanistic validation to prioritize genotoxic risks from antibiotic mixtures in food.

Findings

Prioritized mixtures of antibiotics reduced cell viability and induced DNA damage in hepatocytes.

BCL2 overexpression significantly reduced DNA damage, confirming its protective role.

ML revealed that exposure levels, not just toxicity, are key drivers of risk.

Abstract

Antibiotic residues in aquatic foods pose genotoxic risks. Traditional monitoring focuses on individual Maximum Residue Limit (MRL) compliance, often overlooking cumulative multi-residue risks. We developed an interpretable machine learning (ML) framework integrating surveillance data (3719 samples, 17 antibiotics) with mechanistic validation to prioritize risks in Guangzhou (2021−2023). Exposure metrics and in silico hazard predictions were analyzed via clustering and ranking. Enrofloxacin, sulfamethoxazole, and 3-amino-2-oxazolidinone were prioritized, with risk drivers deconstructed via Explainable AI (SHAP). In L-02 hepatocytes, prioritized mixtures reduced viability (70.21 ± 7.49%), increased apoptosis (7.12 ± 2.75%), and induced DNA damage (tail DNA% 5.25 ± 1.03%) (all p < 0.05). BCL2 overexpression significantly attenuated this damage (tail DNA% 2.90 ± 2.65%, p < 0.05), confirming its role as a key functional mediator. This surveillance-to-mechanism workflow provides a data-driven paradigm for identifying mechanistic biomarkers and prioritizing food safety interventions, surpassing traditional compliance-based monitoring.

•Integrated exposure-toxicity metrics via ML to prioritize residues in aquatic foods.•SHAP revealed exposure outweighs intrinsic toxicity as primary risk driver.•Validated antibiotic mixtures induce DNA damage via BCL2-mediated apoptosis.•ML uncovered “risk-shifting” from restricted ENR to emergent SMZ violations.•A replicable paradigm translating big data into actionable policies.

Integrated exposure-toxicity metrics via ML to prioritize residues in aquatic foods.

SHAP revealed exposure outweighs intrinsic toxicity as primary risk driver.

Validated antibiotic mixtures induce DNA damage via BCL2-mediated apoptosis.

ML uncovered “risk-shifting” from restricted ENR to emergent SMZ violations.

A replicable paradigm translating big data into actionable policies.

Figures

Overall risk profile and regulatory landscape of antibiotic residues in Guangzhou aquatic foods (2021–2023). (A) Risk landscape. A bubble chart illustrating the detection rate (x-axis), exceedance rate (y-axis), and mean concentration (bubble size) for each antibiotic. (B) Contribution to total residue load. A treemap visualizing the percentage contribution of individual antibiotics to the total quantified residue mass across all positive detections. (C) Detection vs. conditional exceedance rates. A dumbbell plot comparing the overall detection rate with the conditional exceedance rate (the pr

Spatiotemporal distribution and contamination hotspots of antibiotic residues. This figure illustrates the spatial distribution of antibiotic contamination across different aquatic food categories and species, alongside its temporal evolution. (A) Overall Contamination Rate by Category. A bar chart showing the overall detection rate of antibiotics for each aquatic food category. (B) Detection Rate Heatmap. A heatmap displaying the detection rates (%) of the top 10 antibiotics across key aquatic food categories. (C) ENR in Freshwater Fish Species. A lollipop chart showing the number of ENR dete

Integrated risk prioritization and driver analysis of antibiotic residues. (A) Exposure–hazard quadrant visualizing *k*-means clustering (*k* = 4) of 14 antibiotics based on standardized features (mean concentration, detection rate, ADMET-predicted toxicity, bioaccumulation); ENR (red) occupies the high-exposure/high-hazard quadrant. (B) TOPSIS risk scores ranking 12 antibiotics with complete datasets, identifying ENR, SMZ, and AOZ as top three. (C) Radar plot of standardized indicator profiles for the top five TOPSIS-ranked antibiotics. (D) Confusion matrix comparing domain knowledge (quadran

Toxicological risk assessment and in vitro toxicity verification of prioritized antibiotics. (A)Toxic equivalent quantity (*TEQ*) of 6 antibiotics interacting with receptors. (B) L-02 cell viability, (C) apoptosis levels, (D) tail DNA percentage (Tail DNA%), (E) tail length (μm), and (F) Olive tail moment in L-02 cells following 24-h treatment with antibiotic mixtures. (G) Representative flow cytometry plots illustrating apoptosis. (H) Representative fluorescence microscopy images of DNA damage from the comet assay. Values are expressed as mean ± SD (n = 3 independent biological plicates). Sta

Network toxicology and mechanistic validation reveal *BCL2* as a key mediator of antibiotic-induced genotoxicity. (A) Relationship diagram of the three antibiotics and their predicted protein targets. (B) Venn diagram showing target gene overlap among the three antibiotics. (C) Venn diagram identifying the 274 genes common to both antibiotic targets and the genotoxicity-related gene set. (D–F) Functional enrichment analysis of the 274 intersecting genes, showing top GO terms and KEGG pathways. (G) PPI network used to identify 10 core hub genes. (H) A representative molecular docking model of a

Graphical Abstract. ML-prioritized antibiotic residues (ENR, SMZ, AOZ) in Guangzhou aquatic foods (2021–2023) drive BCL2-mediated genotoxicity; surveillance-to-mechanism paradigm for targeted food safety regulation.Fig. 6

Full text

Introduction

1

The escalating global demand for high-quality protein has driven intensive aquaculture practices, significantly increasing reliance on antibiotics for disease control in high-density farming systems (Liu et al., 2017; Wang et al., 2022). Global antibiotic consumption in aquaculture exceeded 10,000 metric tons in 2017 and is projected to rise by 33% by 2030, with Asia, particularly China, accounting for the dominant share (Schar et al., 2020; Shao et al., 2021). While essential for production, this widespread use leads to unavoidable antibiotic residue accumulation in aquatic foods and the environment (Ben et al., 2019; Li et al., 2021). Numerous studies have documented concerning levels of antibiotics, such as chloramphenicol (CPL) and nitrofuran metabolites, in fish and shellfish, often exceeding regulatory limits (Luo et al., 2021).

These residues pose significant public health threats. Beyond the well-established concern of fostering antibiotic-resistant bacteria that can transfer to humans (Arsène et al., 2022), chronic dietary exposure to antibiotic residues has been linked to potential carcinogenicity, allergic reactions, developmental toxicity in children, and gut microbiota disruption (Cazer et al., 2020; Ravindra et al., 2023). Critically, the long-term health implications of low-dose, multi-residue exposure remain poorly characterized. While international standards exist (e.g., Codex Alimentarius), effectively managing the complex mixture of antibiotic residues in aquatic foods remains a formidable challenge (Okocha et al., 2018). Conventional monitoring, often focused on compliance with Maximum Residue Limits (MRLs) for individual compounds, struggles to objectively stratify risks arising from diverse exposure patterns and varying intrinsic toxicities within complex environmental matrices.

Recent advancements in data science and bioinformatics offer transformative potential. Machine learning (ML), particularly unsupervised clustering, enables unbiased, data-driven stratification of chemical risks (Li et al., 2025; Revelou et al., 2025; Singh et al., 2024). However, most current environmental ML applications are confined to the ‘black-box’ prediction of antibiotic occurrences or concentrations, often lacking a systematic bridge between large-scale monitoring data and mechanistic toxicology. Concurrently, Explainable Artificial Intelligence (XAI) techniques, such as SHapley Additive exPlanations (SHAP), facilitate a transition from mere prediction to robust ‘driver analysis.’ By quantifying the specific contributions of exposure metrics versus hazard features, XAI serves as a critical bridge that links macroscopic monitoring data with microscopic mechanistic toxicology (Marzidovšek et al., 2024; Vega García & Aznarte, 2020; Xiong et al., 2024). Furthermore, network toxicology and molecular docking allow for the systematic identification of molecular targets and pathways underlying toxicity (Chu & Zi, 2024; He et al., 2024). However, a critical gap exists in translating large-scale environmental surveillance data into actionable, risk-prioritized regulatory strategies using these integrated computational and mechanistic approaches.

Guangzhou, a major hub for aquaculture and consumption in China, exemplifies these challenges. Assessing the dynamic public health risk from antibiotic residues here requires moving beyond simple residue screening to address fundamental questions: (1) How do contamination patterns and key risk drivers (e.g., high exposure vs. high intrinsic toxicity) vary spatiotemporally and across species in response to practices and regulations? (2) What are the co-occurrence patterns of residues, and which combinations drive the highest cumulative risk? Simply targeting the highest-concentration compounds or predefined toxicant lists is insufficient, as it overlooks the complex interplay of exposure prevalence, bioaccumulation potential, and mixture effects that determine actual biological impact.

Therefore, this study pioneers an integrated, multi-tiered framework designed to bridge these gaps, systematically transforming surveillance data into biomarker-guided regulatory insights. Our approach progresses from macro-level spatiotemporal monitoring to micro-level mechanistic understanding: First, to comprehensively map the contamination landscape, we performed 23,483 analytical tests for 17 antibiotic residues across a large-scale collection of 3719 aquatic food samples obtained from Guangzhou over a three-year period (2021–2023). Second, we developed and applied an interpretable ML pipeline, integrating empirical exposure metrics (concentration, detection rate) with in silico hazard predictions (bioaccumulation, toxicity), to objectively prioritize antibiotic residues based on their holistic risk profile. Leveraging XAI (SHAP), we conducted a robust driver analysis to deconstruct the factors influencing this prioritization. Finally, for the ML-prioritized compounds, we employed network toxicology, molecular docking, and in vitro validation to elucidate their underlying genotoxic mechanisms, identifying a key biomarker (BCL2) for potential monitoring.

This study aims not only to provide a robust scientific basis for targeted, risk-based management of antibiotic residues in aquaculture but also to establish a replicable “ML-to-mechanism” paradigm for evidence-based chemical risk governance in the food chain. By demonstrating how ML-driven risk quadrants and biomarker insights can directly guide regulatory resource allocation (e.g., focusing inspections on high-priority species or compounds), we address a critical need in food safety governance: transforming big data into prioritized regulatory action.

Materials and methods

2

Chemicals and materials

2.1

The 17 antibiotics (Table S1) were obtained from MedChemExpress (MCE, China). Dimethyl sulfoxide (DMSO; CAS 67–68-5), used as solvent, was supplied by MP Biomedicals (Shanghai, China). Stocks were maintained at 4 °C. Cellular assay reagents included Cell Counting Kit-8 (APExBIO, Shanghai, China), GoTaq® qPCR Master Mix (Promega, Madison, WI, USA), and the Annexin V-PE/7-AAD Kit (KeyGEN BioTECH, Jiangsu, China).

Sample collection

2.2

A surveillance program (2021–2023) collected 3719 commercial aquatic samples from Guangzhou markets using a stratified random sampling strategy. The samples were immediately transported to the laboratory under refrigeration. Sample preparation was performed according to routine National standards. For analysis, samples were peeled, shelled, or deboned, and edible parts were homogenized. Since all samples were bought as deceased products for consumption, no ethical approval was required. Samples were collected monthly using a stratified random sampling strategy based on market size and geographic distribution across all 11 administrative districts.

The comprehensive sampling covered four major categories (freshwater foods, marine foods, shellfish, and other aquatic products) and encompassed 814 economically important species, providing a snapshot of the region's antibiotic residue landscape.

Chemical analysis

2.3

Fish muscle tissues were homogenized and accurately weighed (2.0 g per sample). Antibiotics were extracted using a mixture of acetonitrile and water (80:20, v/v) containing 0.1% formic acid, followed by vortexing and ultrasonic-assisted extraction for 20 min. The extracts were centrifuged at 10,000 ×g for 10 min at 4 °C. The supernatants were purified using Oasis HLB solid-phase extraction (SPE) cartridges preconditioned with methanol and water. After elution and drying under nitrogen, residues were reconstituted in 1 mL of 10% methanol in water and filtered through a 0.22 μm membrane filter. Multiple antibiotics were quantified via ultra-high-performance liquid chromatography coupled with tandem mass spectrometry (UHPLC-MS/MS) on a C18 reverse-phase column (2.1 × 100 mm, 1.7 μm, Milford, MA, USA) at 40 °C. The mobile phase consisted of (A) 0.1% formic acid in water and (B) 0.1% formic acid in methanol, delivered at a flow rate of 0.3 mL/min under gradient elution. The injection volume was 5 μL. MS/MS detection was conducted using an electrospray ionization source (ESI) in positive mode with multiple reaction monitoring (MRM). Instrumental parameters (ion spray voltage, gas flows, and source temperature) were optimized for each analyte. Quantification was performed using external calibration curves prepared in matrix-matched standards.

For quality assurance and control, procedural blanks, matrix spikes, and quality control (QC) samples were included in each batch. Isotope-labeled internal standards were added to all samples prior to extraction to correct for matrix effects and signal variability. Method recovery, linearity, limit of detection (LOD), and limit of quantification (LOQ) were validated according to standard protocols. Recovery rates for antibiotics ranged between 70% and 120%, with relative standard deviations (RSDs) below 15%.

Risk prioritization and machine learning modeling

2.4

Antibiotics were prioritized using a tiered evaluation of exposure and hazard. Risk profiles were initially constructed using four primary indicators: mean concentration and detection rate derived from surveillance, alongside bioaccumulation factor and toxicity score predicted by ADMETlab 3.0. These features were standardized via Z-score transformation, and K-means clustering (k = 4) was applied to partition residues into four risk quadrants.

Because quadrant assignment provided only a qualitative stratification, the Technique for Order of Preference by Similarity to Ideal Solution (TOPSIS) was subsequently applied to derive a more granular composite risk ranking. To achieve this high-resolution prioritization, the evaluation was expanded to eight standardized indicators: mean concentration, detection rate, exceedance rate, species richness, regulatory status, bioaccumulation factor, and specific toxicity endpoints (carcinogenicity and genotoxicity). TOPSIS was executed in R as follows: (i) a decision matrix was constructed using antibiotics as rows and the eight indicators as columns; (ii) equal weights were assigned to reflect balanced contributions; (iii) the positive (PIS) and negative (NIS) ideal solutions were defined; (iv) normalized Euclidean distances were calculated; and (v) relative closeness coefficients were computed, yielding scores from 0 to 1.

To ensure the objectivity of the prioritization and avoid synthetic bias, no data imputation was performed. A fundamental technical distinction was established between non-detections and missing database parameters: non-detectable concentrations were treated as numerical zeros, which are valid inputs for calculation, whereas compounds lacking structure-specific descriptors were treated as mathematical nulls. Consequently, a stepwise completeness-driven exclusion process was applied: entries lacking primary structural or toxicological descriptors were deterministically excluded from initial clustering, while those with gaps in secondary hazard metrics were further omitted from the final quantitative ranking due to the mathematical requirement for a strictly complete data chain.

An extreme gradient boosting (XGBoost) model was implemented using the caret and xgboost packages in R to approximate K-means labels. Standardized features served as predictors, and cluster assignments as categorical outcomes. The model was trained on the full dataset using 5-fold cross-validation repeated three times (‘repeatedcv’, number = 5, repeats = 3). Hyperparameters were optimized via grid search across the following search space: nrounds (50, 100, 150), max_depth (1,2,3), eta (0.1, 0.3, 0.4). To control model complexity and prevent over-tuning, other parameters were held constant: gamma (0), colsample_bytree (0.8), min_child_weight (1), and subsample (1). The optimal parameter combination (nrounds = 100, max_depth = 3, and eta = 0.3) was selected based on the highest average cross-validation accuracy. The trained XGBoost model was utilized as a descriptive framework to facilitate subsequent interpretability analysis.

Model interpretability analysis

2.5

The model's decision-making process was deconstructed using SHapley Additive exPlanations (SHAP), a game-theoretic approach implemented via the kernelshap package in R. For each antibiotic, SHAP values were calculated for every feature, with the training data serving as the background dataset to define the baseline model prediction. The resulting values were visualized using the shapviz package to generate a suite of interpretability plots: (i) a global feature importance plot based on mean absolute SHAP values; (ii) SHAP summary (beeswarm) plots showing the distribution and direction of feature impacts; (iii) feature dependence plots revealing non-linear relationships; and (iv) a waterfall plot for a representative sample, providing a detailed, additive explanation of a single prediction. This comprehensive SHAP analysis enabled a multi-layered interpretation of the antibiotic risk drivers at both global and local levels.

Toxicological risk assessment and in silico toxicity prediction

2.6

We began by calculating the Estimated Daily Intake (EDI) of antibiotic residues from consuming aquatic foods to evaluate the potential health risks. Acknowledging variations in dietary habits and physiological tolerance, the assessment was stratified for distinct population subgroups: children (<18 years) and adults (≥18 years), as well as urban and rural residents, based on parameters from the Exposure Factors Handbook of the Chinese Population. The EDI was calculated based on the following equation (Gao et al., 2025):

$$[eqn]$$

where Df represents daily fish consumption (37.86 g/day for children and 32.73 g/day for adults), C represents the median level of the specific antibiotic in aquatic foods (μg/kg), and BW refers to the average body weight (60.8 kg for rural populations and 63.4 kg for urban populations).

To translate intake into a biologically relevant risk metric, a Toxicity Equivalency Quantity (TEQ) was calculated. First, the chemical structures of antibiotics were initially obtained from the PubChem database. Subsequently, six toxicity endpoints (NR-ER, SR-ATAD5, NR-AR, NR-AR-LBD, NR-PPAR-γ, SR-ARE) were predicted using the ADMETlab 3.0 platform (https://admetmesh.scbdd.com/) to obtain their corresponding Toxicity Equivalency Factors (TEFs) (Table S2). The TEQ was then calculated as follows (Gao et al., 2025):

$$[eqn]$$

Based on the above results, three priority antibiotics were screened out, and then the organ toxicity was predicted by two independent platforms, ADMETlab 3.0 (https://admetlab3.scbdd.com/) and Deep-PK Predictions (https://biosig.lab.uq.edu.au/deeppk/).

Antibiotic grouping

2.7

To simulate human exposure to antibiotics via aquatic foods, three dose groups were designed. The steady-state plasma concentration (Css) was calculated as follows (Flynn, 2007; Kum et al., 2024):

$$[eqn]$$

In this equation, the parameter F represents the bioavailability, which was predicted using the SwissADME database. The maximum daily intake, D, was derived from the highest measured antibiotic residue (as per Eq. 1), while CL corresponds to the clearance rate predicted by ADMETlab. To establish biologically relevant in vitro exposure levels, the calculated Css was multiplied by a composite uncertainty factor (UF) of 1000, consistent with established risk assessment guidelines (Stedeford et al., 2007) and recent methodological precedents in in vitro-to-in vivo extrapolation (Liu et al., 2024; Romeo et al., 2022). This 1000-fold amplification was employed to bridge the gap between acute cellular response and chronic human exposure, accounting for inter-individual variability (10-fold), uncertainties in predictive models (10-fold), and chronic bioaccumulation effects (10-fold). The final experimental dose groups are detailed in Supplementary Table S3.

Cell culture and cell viability assessment

2.8

Human normal hepatocyte L-02 cells (also known as HL-7702) were obtained from the Guangdong Provincial Center for Disease Control and Prevention (Guangzhou, China). Cells were cultured in RPMI-1640 medium supplemented with 12% fetal bovine serum (FBS) and maintained at 37 °C in a humidified incubator with 5% CO₂. Antibiotic stock solutions were prepared in DMSO and further diluted with culture medium to the required concentrations, ensuring a final DMSO concentration of ≤0.1% (v/v).

Cell viability was assessed using the Cell Counting Kit-8 (CCK-8) assay. L-02 cells were seeded at 8000 cells per well in 96-well plates and treated for 12, 24, or 48 h. CCK-8 reagent was added to the culture medium (1,10 ratio), incubated at 37 °C for 1 h in the dark, and absorbance was measured at 450 nm using a BioTek microplate reader. Each condition was performed in six technical replicates and repeated independently three times.

Apoptosis assay

2.9

L-02 cells were seeded at 5 × 10^5^ cells per well in 6-well plates and incubated overnight. Cells were treated with antibiotic mixtures for 24 h, harvested using trypsin without ethylenediaminetetraacetic acid (EDTA), and collected in 1.5 mL centrifuge tubes. Cell pellets were rinsed twice with precooled phosphate-buffered saline (PBS) by centrifugation to remove supernatant.

For staining, cell pellets were resuspended in 100 μL binding buffer, followed by addition of 5 μL 7-aminoactinomycin D (7-AAD) and 1 μL Annexin V-phycoerythrin (PE). Suspensions were vortexed gently and incubated in the dark for 5–10 min, then diluted with 400 μL binding buffer. Samples were analyzed by flow cytometry. Assays were performed in three biological replicates.

Comet assay

2.10

L-02 cells were suspended in low-melting-point agarose and applied onto slides precoated with normal-melting-point agarose. After solidification at 4 °C, cells were lysed at 4 °C in the dark for 1 h, rinsed, and equilibrated in cold alkaline electrophoresis buffer for 20 min to allow for DNA unwinding. Electrophoresis was performed at 25 V for 30 min. Subsequently, the slides were neutralized with Tris-HCl and stained with NA-Red (1500 dilution).

DNA damage was quantified using fluorescence microscopy (150 cells randomly selected per group) via CASP software. Parameters assessed included tail DNA percentage (tail DNA%), tail length, and olive tail moment. Assays were performed in three biological replicates.

Construction of the antibiotic target library

2.11

To construct a comprehensive target library for the selected antibiotics, we systematically mined five well-established databases: ChEMBL (https://www.ebi.ac.uk/chembl/), STITCH (http://stitch.embl.de/), PharmMapper (https://www.lilab-ecust.cn/pharmmapper/), DrugBank (https://go.drugbank.com/), and SwissTargetPrediction (http://www.swisstargetprediction.ch/). We limited the searches to Homo sapiens to identify possible human protein targets. We then aggregated the target lists retrieved from all sources. This combined dataset underwent a rigorous curation process, which included the removal of duplicate entries and the standardization of all target gene names against the UniProt database to ensure nomenclature consistency. The final, curated set of unique proteins constituted the antibiotic target library used for subsequent analyses.

Construction of the genotoxicity-related gene set

2.12

To establish a comprehensive gene set associated with genotoxicity, genes related to “Genotoxicity” were obtained from GeneCards (https://www.genecards.org/), while additional relevant genes were retrieved from OMIM (https://omim.org/) using the terms “Genetic Diseases” and “Genetic Toxicity”. The retrieved gene lists were merged, and duplicate entries were removed to ensure a unique reference set. This final, curated collection constituted the genotoxicity-related gene library, serving as the foundation for identifying the overlapping targets between antibiotic exposure and genotoxic risk.

Enrichment analyses of gene ontology (GO) and Kyoto Encyclopedia of genes and genomes (KEGG) pathways

2.13

We conducted GO and KEGG pathway enrichment analyses on the genes shared between the antibiotics' targets and genotoxicity-related genes. The analyses were carried out in the R software utilizing the ClusterProfiler and Org.Hs.eg.db packages. For biological annotations, terms and pathways with p < 0.05 were selected.

Protein-protein interaction (PPI) network construction

2.14

We constructed a PPI network to map the connections among the overlapping targets identified between the antibiotic mixtures and genotoxicity. We generated this network using the STRING database (https://string-db.org/). After exporting the interaction data, we visualized and refined the PPI network using Cytoscape version 3.9.1, with its MCODE plugin employed to establish the network's structural layout.

Core target screening

2.15

Based on the intersection of genes ranked first in the MCODE clustering and those identified using the CytoNCA and cytoHubba plugins, we obtained the hub genes in the PPI network. The key genes were screened out based on the median of “Betweenness”, the median of “Closeness”, and twice the median of “Degree”.

Molecule-protein binding assay

2.16

We obtained the 3D coordinates of proteins from PDB database (http://www.rcsb.org/pdb). AutoDock Tools 1.5.6 software was used to preprocess both the ligand molecules and protein. The data were input into AutoDock Vina software to carry out molecular docking calculations. The binding energy of protein and ligand of each pose was calculated by the software automatically and the binding poses of ligand and protein was viewed and analyzed by PyMOL 3.1.3 software.

Cell transfection

2.17

The coding sequence of human BCL2 was cloned into the pEZ-Lv105 vector (EX-H3307-Lv105; GeneCopoeia) (Table S6). An empty vector (EX-NEG-Lv105) served as the negative control. Cells were seeded in six-well plates (70%–90% confluency), and transfected using Lipofectamine 2000 and Opti-MEM according to a 1:1 ratio with plasmid DNA. After 30 min at room temperature, the complexes were added to cells and incubated for 8 h. Puromycin selection was applied post-transfection to obtain stable cell lines. At 72 h after transfection, the overexpression was confirmed by quantitative reverse transcription-polymerase chain reaction (qRT-PCR).

qRT-PCR

2.18

Using TRIzol reagent, total RNA was extracted from L-02 cells. After lysis, chloroform was added to the mixture, and the aqueous phase containing RNA was obtained by centrifugation. RNA was subsequently precipitated with isopropanol, and the resulting pellet was collected and subjected to washes with 75% chilled ethanol and anhydrous ethanol. The purified RNA was then redissolved for subsequent steps.

cDNA was synthesized from the extracted RNA using the GoScript™ Reverse Transcription System for subsequent gene expression analysis. Real-time PCR was carried out on the Applied Biosystems™ 7500 Fast Dx instrument (Thermo Fisher Scientific) with GoTaq® qPCR Master Mix, and relative expression levels were quantified using the 2^−ΔΔCT^ method with GAPDH as the internal reference gene. All target gene primer sequences are detailed in Supplementary Table S4.

Statistical analysis

2.19

Descriptive statistical analyses of the surveillance data, including the characterization of spatiotemporal, co-occurrence, and regulatory compliance patterns, were conducted in R (version 4.4.3). All visualizations were created using RStudio and the ggplot2 ecosystem, including packages such as treemapify, ggalluvial, and UpSetR.

For the in vitro experimental data, statistical analyses were carried out using GraphPad Prism version 8.0. Results are expressed as mean values accompanied by standard deviation (SD). To evaluate differences between two groups, independent-samples t-tests were applied, while one-way ANOVA was used for comparisons involving more than two groups. Statistical significance was defined as a p-value <0.05. Each experiment was independently repeated three times.

Results

3

Spatiotemporal distribution, risk profile, and co-occurrence patterns of antibiotic residues

3.1

A total of 3719 aquatic food samples, collected from Guangzhou markets between 2021 and 2023 and representing 814 species across four categories (freshwater foods, marine foods, shellfish, and other aquatic species), underwent 23,483 tests for 17 antibiotics. Nine antibiotic residues were detected, with concentrations ranging from not detected to a maximum of 9426 μg/kg (Table 1).Table 1. Concentration (μg/kg) of antibiotics in different fish species collected from Guangzhou.Table 1. AntibioticsFreshwater foodsMarine foodsShellfishOther aquatic foodsConcentrationDetected rate (%)ConcentrationDetected rate (%)ConcentrationDetected rate (%)ConcentrationDetected rate (%)MinMaxMinMaxMinMaxMinMaxFF3.5108.80.79–––0.12824.150.1119.313.1ENR1.38942641.721.98107411.92.5217622.233.6248240.74TMP13280.2414390.69––––––SMZ2.28714.019.711510.824.824.820.244.8845813.33OTC4944941.3228.32522.38––––––CPL0.350.540.190.641900.212.882.880.1NDND0OTC/CTC/TC4944940.813.53661.09––––––SEMNDND01.22120.67NDND0NDND0AOZ1.3633.50.581.6343.80.64NDND012.212.20.47AMOZNDND0NDND0NDND0–––AHDNDND0NDND0NDND0NDND0OFLNDND0NDND0NDND0NDND0MNZNDND0NDND0––––––NORNDND0NDND0NDND0NDND0PEFNDND0NDND0NDND0NDND0CTC–––NDND0––––––TC–––NDND0––––––Sum (Total Tests)ND94266.04ND11511.75ND17621.14ND4826.14ND: Not Detected.

Detection rates differed by category (Table S8). Freshwater foods showed an overall detection rate of 6.04% (n = 462/7645), which exceeded rates in marine foods (1.75%, n = 161/9195; χ^2^ test, p < 0.001) and shellfish (1.14%, n = 65/5698; χ^2^ test, p < 0.001). Enrofloxacin (ENR) accounted for much of the variation, with a detection rate of 41.72% (n = 413/990) in freshwater foods (concentration range: 1.38–9426 μg/kg) compared to 11.9% (n = 123/1034) in marine foods (range: 1.98–1074 μg/kg; Fisher's exact test, p < 0.001). At the species level, ENR reached 9426 μg/kg in bass, while sulfamethoxazole (SMZ) attained 1151 μg/kg in leopard coral grouper.

Across all 746 positive detections, ENR comprised 602 instances (80.7%), corresponding to an overall detection rate of 19.3% (n = 602/3119 tests for ENR). This rate surpassed those for florfenicol (FF; 3.8%, n = 55/1452; χ^2^ test, p < 0.001) and SMZ (2.1%, n = 43/2006; χ^2^ test, p < 0.001). Of 84 non-compliant samples (0.36% of total; one-sample proportion test, p < 0.001), ENR contributed 42 (50.0%), yielding an exceedance rate of 1.35% (n = 42/3119 samples tested against ENR maximum residue limits; one-sample proportion test, p < 0.001) (Fig. 1A). Notably, ENR accounted for a dominant 84.5% of the total residue mass load across detections (weighted mean concentration: 82.55 μg/kg) (Fig. 1B).Fig. 1. Overall risk profile and regulatory landscape of antibiotic residues in Guangzhou aquatic foods (2021–2023). (A) Risk landscape. A bubble chart illustrating the detection rate (x-axis), exceedance rate (y-axis), and mean concentration (bubble size) for each antibiotic. (B) Contribution to total residue load. A treemap visualizing the percentage contribution of individual antibiotics to the total quantified residue mass across all positive detections. (C) Detection vs. conditional exceedance rates. A dumbbell plot comparing the overall detection rate with the conditional exceedance rate (the proportion of detected samples that are non-compliant) for the top 10 most detected antibiotics. (D) Detection rate by regulatory status. A jittered boxplot comparing the distribution of detection rates for ‘Banned’ versus ‘Restricted Use’ antibiotics (log-scaled y-axis). (E) Co-occurrence patterns. An UpSet plot displaying the frequency of single and intersecting sets of detected antibiotics.Fig. 1

Regulatory status affected detection and exceedance patterns (Fig. 1B–D). Banned substances, including nitrofuran metabolites (semicarbazide [SEM], 3-amino-2-oxazolidinone [AOZ]) and chloramphenicol (CPL), exhibited low detection rates (e.g., CPL: 0.16%, n = 6/3656) but 100% conditional exceedance upon detection (n = 6/6). Restricted-use antibiotics accounted for 96.5% of detections (n = 720/746) and residue mass (Supplementary Fig. S1). While specific restricted-use compounds showed high prevalence, the median detection rates for restricted-use antibiotics did not significantly differ from those for banned substances (Wilcoxon rank-sum test, p = 0.642) (Fig. 1D).

Co-occurrences appeared in 12.3% of positive samples (n = 76/616). Single-residue detections predominated (n = 540/616, 87.7%). The ENR-based combinations were the most frequent co-occurrence patterns, with one quadruple co-occurrence of SMZ-CPL-AOZ-ENR (n = 1) (Fig. 1E). Overall, surveillance detected nine antibiotic residues in 746 instances across 3719 samples, with ENR comprising 80.7% of positives and 84.5% of total quantified mass.

Spatiotemporal distribution reveals critical hotspots and dynamic shifts

3.2

Spatially, the contamination was highly heterogeneous across different food categories. Freshwater Fish emerged as the primary reservoir for antibiotic residues, exhibiting the highest overall detection rate of 6.9%, significantly higher than Other Aquatic (6.1%), Marine Shrimp (2.3%), Marine Fish (1.9%), Freshwater Shrimp (1.3%), and Shellfish (1.1%) (Fig. 2A). Marine Crab and Freshwater Crab showed 0% detection. This disparity underscores that freshwater environments are critical areas for contamination control.Fig. 2. Spatiotemporal distribution and contamination hotspots of antibiotic residues. This figure illustrates the spatial distribution of antibiotic contamination across different aquatic food categories and species, alongside its temporal evolution. (A) Overall Contamination Rate by Category. A bar chart showing the overall detection rate of antibiotics for each aquatic food category. (B) Detection Rate Heatmap. A heatmap displaying the detection rates (%) of the top 10 antibiotics across key aquatic food categories. (C) ENR in Freshwater Fish Species. A lollipop chart showing the number of ENR detections in the top 10 contaminated freshwater fish species, with point size indicating maximum concentration. (D) Overall Annual Trends. A line chart illustrating the annual overall antibiotic detection and exceedance rates from 2021 to 2023. (E) ENR Trend in Bass. A line chart showing the annual detection rate of ENR specifically in Bass from 2021 to 2023. (F) Shift in Exceedance Composition. An alluvial plot depicting the annual changes in the composition (number of exceedances) of different antibiotics that caused non-compliance.Fig. 2

A detailed heatmap (Fig. 2B) further pinpointed the epicenter of this risk at the intersection of specific pollutants and aquatic food categories. The analysis overwhelmingly highlighted ENR in Freshwater Fish as the most critical risk hotspot, with a staggering detection rate of 47.9%. Other notable hotspots included ENR in Other aquatic foods (40.7%), SMZ in Other aquatic foods (13.3%), and FF in Other aquatic foods (13.1%). This comprehensive matrix provides a granular view of where specific antibiotic contamination issues are most concentrated.

Drilling down to the species level within these hotspots (Fig. 2C), a detailed investigation identified specific aquatic species as principal culprits for ENR contamination in freshwater fish. Bass was the leading species, accounting for 87 instances of ENR detections. Grass carp followed closely as the second major contributor with 82 instances, after consolidating synonymous entries and various tissue samples (e.g., meat and belly). Other frequently contaminated species identified through this precise mapping included crucian carp (34 detections), loach (25 detections), and grouper fish (22 detections). The size of the bubbles in Fig. 2C indicates that some of these species, particularly bass, could carry extremely high concentrations (up to 9426 μg/kg), which poses a significant direct exposure risk and identifying them as high-priority targets for subsequent mechanistic validation.

Temporally, the contamination landscape proved to be dynamic (Fig. 2D and F). The overall annual trends in detection and exceedance rates showed fluctuations (Fig. 2D). While the detection rate was 3.5% in 2021, it slightly increased to 3.7% in 2022 before dropping to 1.7% in 2023. Similarly, the exceedance rate was 0.36% in 2021, peaked at 0.50% in 2022, then declined to 0.20% in 2023. More profoundly, the composition of non-compliant antibiotics underwent a significant shift (Fig. 2F). The proportion of exceedances attributed to ENR decreased steadily over the three years (from 54% in 2021 to 40% in 2023). Concurrently, the relative contribution of SMZ to the total number of exceedances surged in 2023, suggesting a potential “risk-shifting” phenomenon. This dynamic indicates that regulatory pressure on one dominant antibiotic may inadvertently drive increased reliance on others, posing an evolving challenge for food safety management.

The annual trend of ENR specifically in Bass (Fig. 2E) provides a focused view of the primary culprit's dynamics. The detection rate for ENR in Bass consistently remained high, showing an increase from 86% in 2021 to 87.5% in 2022, and reaching 100% in 2023 among the samples tested for Bass, highlighting a persistent and escalating issue in this specific high-risk species.

Integrated risk prioritization of antibiotic residues using interpretable machine learning

3.3

Seventeen antibiotics were initially profiled across the surveillance dataset. Following the screening criteria defined in Section 2.4, 14 antibiotics were visualized in the risk quadrants; AMOZ, MNZ, and the composite group OTC/CTC/TC were excluded due to the absence of primary structural identifiers and validated hazard descriptors. Furthermore, AHD and TC were omitted from the subsequent TOPSIS ranking due to gaps in secondary hazard metrics—specifically, missing bioaccumulation factors and specific carcinogenicity/genotoxicity predictions in the database. This systematic refinement yielded a final set of 12 antibiotics with complete, multidimensional data chains for quantitative prioritization.

K-means clustering (k = 4) was applied to four standardized features (mean concentration, detection rate, ADMET-predicted toxicity, predicted bioaccumulation) and cluster assignments were projected onto a two-dimensional exposure–hazard quadrant (Fig. 3A). Among the 14 antibiotics shown in Fig. 3A, ENR was the only compound positioned in the high-exposure / high-hazard quadrant (Quadrant I; 7.1% of visualized compounds). Eight visualized compounds (8/14, 57.1%) comprised Quadrant II (Low Exposure, High Hazard), and the remaining five compounds (35.7%) occupied the low-hazard quadrants (Quadrants III and IV) .Fig. 3. Integrated risk prioritization and driver analysis of antibiotic residues. (A) Exposure–hazard quadrant visualizing k-means clustering (k = 4) of 14 antibiotics based on standardized features (mean concentration, detection rate, ADMET-predicted toxicity, bioaccumulation); ENR (red) occupies the high-exposure/high-hazard quadrant. (B) TOPSIS risk scores ranking 12 antibiotics with complete datasets, identifying ENR, SMZ, and AOZ as top three. (C) Radar plot of standardized indicator profiles for the top five TOPSIS-ranked antibiotics. (D) Confusion matrix comparing domain knowledge (quadrant) vs. TOPSIS high-risk classification. (E–H) Explainable AI (XAI) results for risk driver analysis, including: (E) mean absolute SHAP values ranking feature importance. (F) SHAP beeswarm plots showing the distribution and direction of feature impacts; (G) SHAP waterfall plot for a representative high-priority sample (total SHAP = 0.432); and (H) SHAP interaction matrix showing correlations between risk drivers. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)Fig. 3

Because quadrant assignment alone did not uniquely identify all candidates, a TOPSIS ranking was computed from the same standardized inputs (Fig. 3B). TOPSIS ranked ENR, SMZ, and AOZ as the top three of 12 compounds. While ENR's ranking was primarily associated with its dominant exposure (84.5% of total residue mass), the high prioritization of SMZ and AOZ was driven by their significantly elevated hazard and toxicity scores. The substantial margin in TOPSIS scores between these top three compounds and the remainder supported their selection for downstream mechanistic analysis.

A radar plot summarized standardized indicator profiles for the visualized antibiotics (Fig. 3C). ENR exhibited the highest standardized values for exposure metrics (detection, exceedance, and concentration); AOZ and SMZ displayed elevated but sub-maximal standardized profiles across exposure and in silico hazard axes relative to the visualized panel, illustrating how the four input features combined to shape each compound's risk signature. To evaluate concordance between the two complementary prioritization approaches, a cross-tabulation confusion matrix compared domain knowledge-based quadrant classifications (from K-means clustering) against TOPSIS-derived risk categorizations (Fig. 3D). The matrix demonstrated clear concordance between the two approaches; the five TOPSIS high-risk compounds were strategically distributed across the prioritized risk zones (one in Quadrant I, two in Quadrant II, and two in Quadrant III). Similarly, the medium-high risk compounds were effectively mapped to Quadrant II (five compounds) and Quadrant III (two compounds). This distribution underscores the synergy between clustering-based stratification and distance-based composite ranking in identifying high-threat contaminants.

An XGBoost surrogate classifier (optimized via grid search as detailed in Section 2.4) was trained to approximate K-means labels, enabling a transparent transition from predictive modeling to robust driver analysis. SHAP analyses (kernelshap, shapviz) quantified feature importance, identifying mean concentration (0.171) as the primary risk driver, followed by toxicity score (0.108) and detection rate (0.103) (Fig. 3E). For the representative highest-risk sample (Fig. 3G; total SHAP sum = 0.432), mean concentration (+0.190) and toxicity score (+0.133) provided the most significant contributions to its risk classification. The SHAP interaction matrix (Fig. 3H) further elucidated the synergy between risk features, showing strong pairwise correlations between mean concentration and detection rate (r = 0.91) and detection rate and toxicity score (r = 0.64). This interpretable pipeline bridges the gap between macroscopic surveillance data and microscopic toxicological risk, providing a data-driven basis for subsequent mechanistic validation.

Based on the integration of quadrant visualization, TOPSIS ranking and XAI-driven driver analysis, ENR, SMZ and AOZ were selected as the high-priority antibiotics for subsequent mechanistic experiments. This selection effectively bridges the gap between large-scale environmental surveillance and the investigation of underlying toxicological responses.

Toxicological risk assessment and genotoxicity prediction of prioritized antibiotics

3.4

Building on the machine learning prioritization of enrofloxacin (ENR), sulfamethoxazole (SMZ), and 3-amino-2-oxazolidinone (AOZ) as high-risk residues from the surveillance dataset, a toxicological risk assessment was performed, with comparisons to other non-compliant antibiotics (semicarbazide [SEM], chloramphenicol [CPL], and oxytetracycline [OTC]) to contextualize their profiles.

Estimated daily intakes (EDI) were calculated for the six antibiotics using median surveillance concentrations and subgroup-specific consumption from the Exposure Factors Handbook of the Chinese Population (Supplementary Table S7 and Fig. S4). EDIs were higher in children (<18 years) and rural residents than adults (≥18 years) and urban residents. The median EDI of all antibiotics was 11.96 ng/kg BW/day for rural child and 9.91 ng/kg BW/day for urban adults.

Toxicity equivalency quantities (TEQ) were derived for six endpoints (NR-ER, SR-ATAD5, NR-AR, NR-AR-LBD, NR-PPAR-γ, SR-ARE) using toxicity equivalency factors (TEFs) predicted by ADMETlab 3.0 based on PubChem structures (Table S2). TEFs for prioritized antibiotics exceeded comparators in NR-ER (ENR 0.541, AOZ 0.672, SMZ 0.034 vs. SEM 0.273, CPL 0.125, OTC 0.00004; one-way ANOVA on log-transformed TEFs, p < 0.001) and SR-ARE (ENR 0.068, AOZ 0.049, SMZ 0.003 vs. SEM 0.004, CPL 0.435, OTC 0.006; p < 0.01). The mean TEFs across endpoints were 1.8-fold higher for ENR, SMZ, and AOZ than comparators (Fig. 4A).Fig. 4. Toxicological risk assessment and in vitro toxicity verification of prioritized antibiotics. (A)Toxic equivalent quantity (TEQ) of 6 antibiotics interacting with receptors. (B) L-02 cell viability, (C) apoptosis levels, (D) tail DNA percentage (Tail DNA%), (E) tail length (μm), and (F) Olive tail moment in L-02 cells following 24-h treatment with antibiotic mixtures. (G) Representative flow cytometry plots illustrating apoptosis. (H) Representative fluorescence microscopy images of DNA damage from the comet assay. Values are expressed as mean ± SD (n = 3 independent biological plicates). Statistical significance was determined by one-way ANOVA: *p < 0.05, **p < 0.01, and ***p < 0.001 versus the control group.Fig. 4

In silico predictions for ENR, SMZ, and AOZ across six endpoints used ADMETlab 3.0 and Deep-PK Predictions (Table 2). High-risk scores (≥0.85) aligned for genotoxicity (ADMETlab 3.0: ENR 1.00, SMZ 1.00, AOZ 0.897) and micronucleus formation (Deep-PK: ENR 0.99, SMZ 0.989, AOZ 0.945). Other endpoints scored lower (e.g., carcinogenicity means: ADMETlab 3.0 0.514 ± 0.235; Deep-PK 0.450 ± 0.366; paired t-test, p = 0.124; no significant platform difference). Hepatotoxicity was elevated (ADMETlab 3.0: ENR 0.98, SMZ 1.00, AOZ 0.633).Table 2. Predictive value of antibiotics toxicity endpoint.Table 2. DatabaseToxicityAssessmentENRSMZAOZADMETlab 3.0Carcinogenicity0.3380.4250.779Drug-Induced Liver Injury0.9810.633Respiratory0.9890.0440.532Skin Sensitization0.1930.0650.762Eye Irritation0.0150.4040.996Genotoxicity110.897DeeP-PK PredictionsCarcinogenicity0.30.1660.883Drug-Induced Liver Injury0.5760.7250.144Respiratory0.9590.3950.91Skin Sensitization0.350.5270.578Eye Irritation0.0040.4740.997Micronucleus0.990.9890.945

L-02 hepatocytes were treated with low-, medium-, and high-dose mixtures (Supplementary Table S3) for 12, 24, and 48 h. Viability (CCK-8 assay; n = 6 technical replicates/condition) decreased at 24 h: medium dose 77.53 ± 1.70% and high dose 70.21 ± 7.49%. At 48 h, high-dose viability was 74.82 ± 1.08% (Fig. 4B; Supplementary Fig. S2).

Apoptosis (Annexin V-PE/7-AAD flow cytometry; n = 3 biological replicates) at 24 h increased: high dose 7.12 ± 2.75% (p = 0.012 vs. control 3.44 ± 0.09%) (Fig. 4C, G).

Comet assay (n = 150 cells/group; n = 3 replicates) at 24 h high dose showed tail DNA% 5.25 ± 1.03% (p < 0.001 vs. control 1.31 ± 0.45%), tail length 7.50 ± 0.79 μm, and Olive tail moment 1.21 ± 0.28% (Fig. 4D–F, H).

Collectively, in silico predictions and in vitro validation provided strong evidence of genotoxicity for the ENR, SMZ, and AOZ mixtures. The 24 h high-dose exposure resulted in both reduced cell viability (to 70.21 ± 7.49%) and significant DNA damage (tail DNA% of 5.25 ± 1.03%), prompting further network toxicology to elucidate the underlying molecular mechanisms.

Network toxicology and mechanistic validation identify BCL2 as a key mediator of antibiotic-induced genotoxicity

3.5

Following the in vitro confirmation that the prioritized antibiotic mixture induces genotoxicity, we employed an integrated network toxicology and molecular validation approach to elucidate the underlying molecular mechanisms.

A comprehensive target library was constructed for the three antibiotics by integrating data from five established databases (Fig. 5A). An analysis of these predictions revealed that 37 targets were common to all three antibiotics (SMZ, ENR, and AOZ), suggesting potential shared mechanisms of toxicity (Fig. 5B). To specifically focus on genotoxicity, this combined target list was intersected with a curated set of 5050 genotoxicity-associated genes. This critical filtering step identified 274 overlapping genes, which are considered the putative mediators of antibiotic-induced genotoxicity (Fig. 5C).Fig. 5. Network toxicology and mechanistic validation reveal BCL2 as a key mediator of antibiotic-induced genotoxicity. (A) Relationship diagram of the three antibiotics and their predicted protein targets. (B) Venn diagram showing target gene overlap among the three antibiotics. (C) Venn diagram identifying the 274 genes common to both antibiotic targets and the genotoxicity-related gene set. (D–F) Functional enrichment analysis of the 274 intersecting genes, showing top GO terms and KEGG pathways. (G) PPI network used to identify 10 core hub genes. (H) A representative molecular docking model of an antibiotic binding to BCL2. (I) Treemap comparing the binding energies of the antibiotics with the ten core hub proteins. (J) Relative mRNA expression of the hub genes in L-02 cells after antibiotic exposure. (K) qRT-PCR confirmation of BCL2 overexpression in stably transfected L-02 cells. (L) Representative comet assay images showing DNA damage in control (Vector) and BCL2-overexpressing (OE-BCL2) cells after treatment. (M) Quantification of comet assay parameters (Tail DNA, Tail Length, Olive Tail Moment). Values are expressed as mean ± SD from three independent experiments. *p < 0.05, **p < 0.01.Fig. 5

To understand the biological functions of these 274 genes, GO and KEGG enrichment analyses were performed. The results revealed significant enrichment in GO terms related to the regulation of cell proliferation and response to steroid hormones (Fig. 5D, E), as well as KEGG pathways including the critical “PI3K-Akt signaling pathway” (Fig. 5F). To pinpoint the central regulators within this network, a PPI network was constructed from the 274 intersecting genes (Fig. 5G). Through centrality analysis of this network, ten core hub genes were identified: AKT1, BCL2, CASP3, CCND1, ESR1, HSP90AA1, MTOR, PARP1, SRC, and STAT3.

To investigate whether the antibiotics could directly interact with these pivotal proteins, molecular docking was performed for all three antibiotics against all ten hub targets. The docking poses for all interactions are provided in the supplementary materials (Fig. S3), with the binding mode for BCL2 presented as a representative example (Fig. 5H). The docking analyses predicted robust binding affinities for all complexes, with binding energies ranging from −11.54 to −16.22 kcal/mol, well below the threshold of −7.2 kcal/mol that indicates strong interaction (Supplementary Table S5). The relative binding energies for each ligand-receptor complex are compared in the tree diagram in Fig. 5I, which graphically illustrates the high interaction strength across all core targets.

To experimentally validate these in silico findings, the transcriptional expression profiles of the ten genes were examined in L-02 cells after exposure to the mixed antibiotics. The qRT-PCR results demonstrated that the mixture significantly downregulated the expression of BCL2, CASP3, and CCND1 (Fig. 5J). Notably, BCL2 was selected for further validation due to its significant downregulation and robust interaction in docking analyses. To definitively establish this functional role, an L-02 cell line stably overexpressing BCL2 (OE-BCL2) was constructed, with qRT-PCR confirming a relative expression increase of 80.89 ± 20.02 (Fig. 5K).

Subsequently, a comet assay was used to assess DNA damage. As expected, the antibiotic mixture induced severe DNA damage in vector+High cells, resulting in a Tail DNA % of 10.99 ± 3.28%. In striking contrast, BCL2 overexpression provided a significant protective effect, with Tail DNA % notably reduced to 2.90 ± 2.65% (Fig. 5L, M). Similar significant reductions were also observed for Tail Length (from 26.00 ± 5.77 μm to 5.50 ± 3.11 μm) and Olive Tail Moment (from 4.64 ± 1.49% to 0.67 ± 0.45%). These results provide strong functional evidence that BCL2 is not merely a correlated biomarker but a key mediator of antibiotic-induced genotoxicity.

Discussion

4

The escalating global demand for aquaculture foods, coupled with intensive farming practices, has led to an alarming dependence on antibiotics, posing a significant environmental health challenge worldwide (Griboff et al., 2020). While the pervasive occurrence of antibiotic residues in aquatic foods is well-documented (Adel et al., 2017; Chen et al., 2018; Kang et al., 2018; Wang et al., 2017), effective risk management remains hindered by the complexity of residue mixtures and the lack of robust, data-driven methodologies for prioritizing compounds based on integrated exposure and hazard profiles. This study directly addresses this critical gap by establishing an integrated “surveillance-to-mechanism” framework. Our large-scale spatiotemporal monitoring, involving 23,483 analyses on aquatic foods from Guangzhou (2021–2023), revealed a contamination landscape dominated by legally permitted but frequently non-compliant antibiotics, particularly ENR. This finding underscores a significant regulatory paradox: while outright bans effectively suppress the detection of prohibited substances like pefloxacin and norfloxacin, controlling the misuse of permitted antibiotics remains a formidable challenge for food safety governance. The exceptionally high concentrations detected for ENR (up to 9426 μg/kg) and SMZ (up to 1151 μg/kg), significantly exceeding their respective MRLs and often surpassing levels reported in previous studies (Chang et al., 2020; Fang et al., 2021; Jansomboon et al., 2016), highlight the inadequacy of compliance-focused monitoring alone. This emphasizes the pivotal need for risk-based prioritization strategies that account for the interplay of prevalence, concentration, bioaccumulation potential, and intrinsic toxicity to guide targeted interventions and resource allocation within complex environmental matrices.

Our detailed spatiotemporal analysis provided critical insights into the geographical and biological specificities driving contamination heterogeneity. Freshwater fish emerged unequivocally as the primary reservoir for antibiotic residues, exhibiting significantly higher detection rates than marine foods. This pronounced disparity likely stems from the intensive production modes, longer farming cycles, and limited self-purification capacity characteristic of freshwater aquaculture systems. Within this high-risk category, our analysis pinpointed ENR in freshwater fish, particularly in bass (Lateolabrax spp.), as the contamination epicenter, characterized by staggering detection rates (up to 100% in bass samples by 2023) and the highest single concentration measured (9426 μg/kg). The co-occurrence patterns, revealing significant associations like the frequent pairing of ENR and SMZ, further illuminate the complexity of real-world mixture exposures consumers face, necessitating consideration of cumulative or synergistic effects in risk assessment.

Temporally, while the overall non-compliance rate showed a decline, a concerning “risk-shifting” phenomenon was evident. The proportion of violations attributed to ENR decreased over the three years, coinciding with a dramatic surge in exceedances caused by sulfonamides, primarily SMZ, in 2023. This dynamic trend strongly suggests that strengthened regulatory pressure on one dominant antibiotic (ENR) inadvertently drove increased reliance on others (SMZ), posing an evolving and adaptive challenge for food safety management. This finding underscores the critical importance of dynamic, data-informed surveillance systems capable of detecting such shifts in near real-time to enable proactive regulatory adjustments.

Moving beyond descriptive statistics, our interpretable ML framework provided a crucial paradigm shift for risk assessment. By integrating empirical exposure metrics (concentration, detection rate) with in silico hazard predictions (toxicity, bioaccumulation), the framework objectively stratified antibiotic residues. Crucially, the application of SHAP deconstructed the model's decision-making, revealing a key insight: exposure metrics (concentration and prevalence) were the dominant drivers of the final risk ranking in this specific environmental context, outweighing intrinsic toxicity scores for many compounds. This implies that, for the studied Guangzhou aquatic ecosystem, reducing the overall environmental load and prevalence of antibiotics would be a more effective immediate risk mitigation strategy than focusing solely on compounds with the highest theoretical toxicity. The unequivocal identification of ENR, SMZ, and AOZ as the top-tier risks was therefore not merely a reflection of their inherent hazard but a data-driven verdict based on their integrated risk profile, quantitatively demonstrating their preeminence as threats due to widespread occurrence and high concentrations. This approach represents a significant advancement from compliance-focused towards risk-focused regulation. The ML-derived risk quadrants and SHAP insights offer a transparent, updatable decision-support tool for regulators. For instance, identifying ENR's risk as heavily driven by its prevalence in freshwater fish (like bass) directly informs targeted monitoring strategies, shifting focus upstream to aquaculture ponds rather than solely relying on market screening. Furthermore, the framework's ability to flag the emergent SMZ risk through feature dependence analysis exemplifies its potential for predictive regulation in response to dynamic contamination patterns like risk-shifting.

A pivotal strength of our integrated framework is bridging the gap between risk identification and mechanistic understanding. For the ML-prioritized compounds (ENR, SMZ, and AOZ), our in silico toxicological risk assessment using TEQ confirmed their higher potential for interacting with critical biological pathways (e.g., estrogen receptor NR-ER, antioxidant response SR-ARE) compared to other non-compliant residues (SEM, CPL, OTC). More importantly, a robust consensus between two independent prediction platforms (ADMETlab 3.0 and Deep-PK Predictions) pinpointed genotoxicity and micronucleus formation as the primary predicted toxicological endpoints for all three prioritized antibiotics. This compelling in silico hypothesis was experimentally validated in vitro. Exposure of human L-02 hepatocytes to a mixture reflecting the prioritized residues induced significant, dose-dependent DNA damage and apoptosis, alongside reduced cell viability. This finding is particularly concerning as it demonstrates genotoxic potential not just for individual known genotoxins (Akhtara & Bharali, 2025; Bariweni et al., 2022; Bhattacharya et al., 2020), but for a realistic mixture of high-priority residues identified from actual field surveillance, suggesting potential additive or synergistic effects relevant to human dietary exposure.

Leveraging network toxicology, we progressed from observing what to elucidating how. The identification of a hub-gene network centered around apoptosis regulation (CASP3, BCL2) and cell cycle control (CCND1) provided a coherent mechanistic blueprint. The significant downregulation of BCL2, CASP3, and CCND1 mRNA levels in antibiotic-exposed L-02 cells aligned with the observed cellular phenotype (DNA damage, apoptosis). Functionally, stable overexpression of BCL2 in L-02 cells significantly rescued them from antibiotic mixture-induced DNA damage, demonstrating that BCL2 is not merely a correlative biomarker but a key functional mediator of the genotoxic response. Mechanistically, the observed suppression of BCL2 likely disrupts the critical BAX/BCL2 rheostat balance, triggering mitochondrial outer membrane permeabilization and the subsequent activation of the downstream caspase cascade. This toxicological pathway mirrors the mechanisms of other established foodborne genotoxins, such as heterocyclic aromatic amines, which similarly induce DNA-damage-dependent apoptosis by altering the BAX/BCL2 ratio (Deng et al., 2022; El-Hefny et al., 2020; Pezdirc et al., 2013). Such parallels suggest that the BCL2-mediated apoptotic axis serves as a convergent and sensitive target for diverse chemical stressors in the aquatic food chain. This mechanistic insight highlights BCL2 as a potential biomarker for monitoring the genotoxic risk associated with antibiotic residues in aquatic foods.

While this study provides a comprehensive framework and significant insights, several limitations warrant acknowledgment to guide future research. Our validation was confined to in vitro models using a single human hepatocyte line, which, while relevant for initial toxicity screening and metabolism, necessitates in vivo studies to confirm genotoxicity and assess chronic low-dose effects in whole organisms. Additionally, expanding beyond genotoxicity to other mechanisms—such as endocrine disruption, immunotoxicity, or gut microbiota dysbiosis—would yield a more holistic risk assessment. At the conclusion of this discussion, Fig. 6 provides a graphical abstract encapsulating the ‘surveillance-to-mechanism’ framework, from ML-driven risk prioritization to BCL2-mediated genotoxicity validation, underscoring its implications for targeted food safety governance.Fig. 6. Graphical Abstract. ML-prioritized antibiotic residues (ENR, SMZ, AOZ) in Guangzhou aquatic foods (2021–2023) drive BCL2-mediated genotoxicity; surveillance-to-mechanism paradigm for targeted food safety regulation.Fig. 6

Conclusion

5

This study establishes an integrated “surveillance-to-mechanism” framework that transforms big data into actionable food safety policies. By harmonizing large-scale monitoring, interpretable machine learning, and mechanistic toxicology, we demonstrate that ENR, SMZ, and AOZ pose the highest threat in Guangzhou's aquatic foods due to their pervasive detection in key species (e.g., bass), not solely intrinsic toxicity. ML-enabled real-time tracking revealed “risk-shifting” to sulfonamides post-ENR restrictions, demanding adaptive governance, implying that dynamic risk management is critical. Furthermore, we validated that prioritized antibiotics induce DNA damage through the dysregulation of the BCL2-centered apoptotic pathway, suggesting the potential of BCL2 as a biomonitoring target for antibiotic genotoxicity. Our paradigm-linking ML risk quadrants to biomarker-validated mechanisms-enables precision regulation (e.g., hotspot monitoring, mixture assessment) for global food safety systems.

CRediT authorship contribution statement

Huangqu Zhu: Writing – review & editing, Writing – original draft, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. Kaili Zhou: Investigation, Formal analysis, Data curation. Yuanzhi Li: Investigation, Data curation. Qingqiong Zhou: Investigation, Data curation. Mingjun Peng: Investigation, Data curation. Xinlan Wu: Project administration, Methodology. Xinwu Mao: Project administration, Methodology. Qiaoyuan Yang: Writing – review & editing, Writing – original draft, Supervision, Resources, Project administration, Methodology, Funding acquisition.

Informed consent statement

Not applicable.

Institutional review board statement

Not applicable.

Funding

This work was supported by the 10.13039/501100003453Natural Science Foundation of Guangdong Province (Grant Number: 2022A1515010727 to Y.Q.) and Key Scientific Research Project of Guangdong Provincial Department of Education (Grant Number: 2022ZDZX2047 to Y.Q.).

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.