Similar to findings in environmental studies [68,69,70,71,72], results of PCA (Figs. S1–S3) illustrate whether statistically significant differences exist between non-contaminated and contaminated samples for both earthworms and lettuce. The most pronounced differences between contaminated and non-contaminated samples are observed in lettuce roots (Fig. S2), where the sample groups are distinctly separated. This is followed by lettuce leaves (Fig. S3), while in the case of earthworm samples, the groups are not clearly separated (Fig. S1); this does not imply the absence of differences; rather, it may be due to the limitations of PCA in capturing subtle variations in this dataset. Notably, most of the significant features underlying the differences between sample groups are linked to alterations in the metabolomes of earthworms and lettuce, involving a diverse range of compounds, including lipids, amino acids, saccharides, and others. Although these changes are pertinent to the overall understanding of the system, they fall outside the scope of this manuscript, which focuses on pharmaceutical metabolites. Furthermore, the assurance of high data quality and reliability in various measuring modes is reinforced by the close proximity and tight clustering of the QC pooled samples. Additionally, the graphical placement of QC samples within the figures, interspersed between different sample types, further contributes to the robustness of the data validation process.

In the process of feature picking, a volcano plot analysis was employed with parameters (p < 0.05 and positive fold change (FC > 3.0)). The use of positive fold change was deliberate, aligning with the focus on identifying only drug-related metabolites, that should not be in blank samples. As in previous studies [7, 19], the FC value was set at 3, mirroring the definition of the limit of detection (LoD), and low concentrations of substances. Subsequently, only features deemed statistically significant and showing a positive fold change were subjected to further evaluation (the quantity of significant features in Table S5). These significant features were cross-referenced with our in silico spectral library, which was designed for the matching of pharmaceutical metabolites. A total of 26 drug-related metabolites were successfully identified (Table 1, respectively Tables S6–S7 with SMILES codes).

The prediction process identified a total of 3762 metabolites using three different approaches: 2704 metabolites were identified through software prediction (SP), 194 metabolites were discovered via literature search (LS), and 864 metabolites were identified using common metabolic pathways (CMP). After conducting the predictions, it became apparent that almost 90% of the predicted metabolites exhibited m/z values lower than 800 (3387 out of 3762). Hence, this threshold was utilized as the upper limit for LC-qTOF analysis to optimize the scan time and enhance the sensitivity of the measurements.

The overlap between the three prediction approaches is evaluated in Fig. 2 using a Venn diagram. Software predictions and literature searches identify exact molecular structures, treating structural isomers as distinct metabolites. In contrast, predictions based on common metabolic pathways provide only tentative structures, without differentiating between structural isomers. To accurately compare the approaches, several structural isomers were often grouped together under a single common metabolic pathway (reaction), as illustrated in Fig. 2. For example, when all three approaches overlapped, the software prediction identified 68 structures, while the literature search yielded 45 structures. All of these metabolites were accounted for by 25 common reaction pathways, due to the generation of structural isomers. This grouping was necessary to ensure a representative and meaningful comparison. The prediction of exact structures through SP can simplify metabolite identification, as it allows for the generation of more precise in silico MS2 spectra, unlike CMP, which only provides tentative structures. Furthermore, the partial overlap between CMP and SP (or also LS) suggests that the available software for metabolite prediction is not limited to simple oxidation, reduction reactions, hydrolysis, or conjugation with amino acids, saccharides, etc. It also predicts the cleavage of parent compounds into smaller molecules, which are not described by the list of common metabolic pathways, although sometimes it can be analytically quite challenging to link these small degradation product to specific small-drug molecules, as these compounds may also naturally occur within the organism. Nevertheless, the software prediction did not encompass all the common metabolic pathways described in the available literature [4, 8, 9, 18, 20, 21, 36, 52,53,54,55,56,57,58,59], contrary to our expectations.

Venn diagram illustrating the overlap among three approaches for metabolite structure prediction, highlighting the number of metabolites predicted by each method

This weak overlap may be due to several factors. The software’s prediction algorithms might be optimized for specific types of reactions or tailored to particular organisms, leading to the omission of pathways described in the literature for different organisms. Additionally, the complexity and variability of biological systems could result in metabolites that are not predicted by the software, even though they are reported in the literature. The algorithms might also focus on common metabolic routes, potentially missing less frequent or more complex pathways. Lastly, there may be a lag in updating prediction models with the latest scientific discoveries.

Following successful metabolite identification, a comparison of the different prediction methods for Met-ID success was visualized in a Venn diagram (Fig. 3 and Table S7). Breaking down the successful identifications, software prediction identified 11 compounds, literature search identified 12 compounds, and common reaction pathways identified 16 compounds. Notably, common reaction pathways resulted in a slightly higher number of successfully identified metabolites. However, each method produced a unique set of compounds, and only 3 of 26 substances were predicted by all three approaches. Given that no single approach demonstrated superiority, it can be concluded that the strategy employing all three methods for prediction was successful. Meanwhile, studies [7, 11, 13, 18,19,20,21,22, 53, 56, 74, 75] rely on either literature search, common metabolic pathways, or a simple blank subtraction in mass spectrometry. The comprehensive use of diverse prediction methods enhances the likelihood of capturing a broader spectrum of metabolites, thus contributing to a more thorough and reliable Met-ID process. Moreover, the overlap between the three approaches for successful metabolite identification (Fig. 3) is quite low, much like the overlap observed in the number of metabolites predicted by each method (Fig. 2). If we examine whether the most intense metabolites were identified by the SP, LS, or CMP approaches, we find that in the case of earthworms (Fig. S43), both metabolites (Atenolol-LS1 and Sulfamethoxazole-LS1) were predicted by all three approaches. In contrast, for L. sativa major metabolites (Figs. S44–S52), the major metabolites (such as Ketoprofen-R87, Ketoprofen-R81, Ketoprofen-R63, Sulfamethoxazole-R63, Enrofloxacin-R63) are formed through glucosidation, a process described in common metabolic pathways, but are not predicted by software predictions. Conversely, metabolites identified solely by SP (such as Enrofloxacin-M242, Erythromycin-M46, Erythromycin-M361, Ketoprofen-M648, and Ketoprofen-M835) exhibited lower ion intensities.

Venn diagram comparing metabolite identification efficiency across different prediction methods, showing the number of metabolites identified by each method

In this analysis, various MS/MS modes were compared on the basis of the number of library matches, which were visualized using a Venn diagram (Fig. 4). Each metabolite:sample type combination was considered as one element, totaling 67 unique elements. ESI + demonstrated slightly more matches (74) than ESI − (55). We postulate that, despite the typical measurement of pharmaceuticals in ESI + , efficient ionization in ESI − was facilitated by conjugation with various endogenous molecules, such as glucose. Both ESI + and ESI − modes yielded both diverse and unique compounds in Met-ID, highlighting the efficiency of utilizing both modes and their complementarity. In ESI + , DDA + (16 unique metabolites) outperformed DIA + (3 unique metabolites), except in the case of earthworm samples, where both MS/MS modes successfully identified 2 pharmaceutical metabolites. This difference may be attributed to the use of a precursor ion list and better MS2 quality due to the low concentrations in the complex matrix. In ESI − , the difference between DIA − (7) and DDA − (3) was less significant, given the overall lower number of identified metabolites. Consequently, combining DDA + and DDA − resulted in 56 of 67 successful identifications while halving the analysis time. Previous studies on the identification of drug-related metabolites have commonly focused on ESI + [11, 13, 18,19,20, 22, 55, 56, 75], whereas in the case of ibuprofen [53] or triclocarban [74], ESI − is utilized. Rarely, both ionization modes are employed [7, 21]. MS/MS mode studies typically use either DDA [7, 11, 19] or DIA [18, 22, 74, 75] mode.

Venn diagram comparing metabolite identification efficiency across different MS/MS modes, showing the number of metabolites identified by each mode

The majority of environmental Met-ID studies [18, 19, 21, 22, 55, 56, 74] typically focus on the identification of metabolites derived from a single parent drug per experiment. To enhance the study throughput, we simultaneously investigated the metabolites of a mixture comprising six pharmaceuticals. In total, 26 drug-related metabolites were successfully identified (Table 1). The distribution among parent substances is as follows: atenolol (1), enrofloxacin (5), ketoprofen (12), sulfamethoxazole (3), tetracycline (2), and erythromycin (3). A pivotal study [76] characterized five identification confidence levels in HRMS for elucidating small molecules and their metabolites. These confidence levels are commonly employed in Met-ID studies [8, 11, 19, 22]. Some studies even synthesized metabolites after successful structure elucidation to attain level 1 confidence (confirmed structure by reference standard) [8, 19]. Within the scope of this study, only pharmaceutical metabolites with identification confidence level 1 (confirmed structure by reference standard), level 2 (probable structure by library spectrum match or experimental data), and level 3 (tentative structures supported with MS2 experimental data or, as in our case, by a library match with in silico MS/MS predictions) were reported. Metabolites with confidence levels 4 and 5 were excluded because of complex matrix, low concentrations, missing MS2 spectra, and overall low data reliability. The overall measured MS2 quality is influenced by the MS/MS measurement mode, with DDA spectral quality being higher than in the case of DIA due to differences in functionality. Additionally, sample matrix and metabolite concentration can negatively impact spectral quality [77]. In some instances, MS2 spectra for both parent drugs and metabolites were obtained in both ESI + and ESI − . All MS2 spectra of the identified metabolites in this study are illustrated in Figs. S4–S42, and their corresponding product ions are listed in Table S6.

For metabolites with known structures, the product ions were found to be in agreement with in silico predictions by CFM-ID 4.0 [60]. When the structure of a metabolite was already described, the observed product ions were consistent with available literature, such as acetyl-sulfamethoxazole (Fig. S6 with study [40]), conjugation of sulfamethoxazole with glucose (Fig. S8 with study [40]), conjugation of sulfamethoxazole with pterin (Fig. S10 with study [41]), atenolol acid (Fig. S12 with study [23]), conjugation of ketoprofen with glucose (Fig. S19 with study [36]), conjugation of ketoprofen with malonic acid and glucose (Fig. S21 with study [36]), and double conjugation of ketoprofen with glucose (Fig. S22 with study [36]).

In instances where only the tentative structure was known, derived from the prediction of the metabolite through common metabolic pathways, the product ions typically exhibited alignment with those of the parent substances. This alignment often included the formation of [M + H]+ (or [M-H]−) ions of the parent substances as intense product ions for the metabolite, especially in cases of phase II metabolism involving conjugation with saccharides or amino acids [21]. Phase III processes may encompass the formation of secondary conjugations, such as malonyl, which serves as a signal for transport from the cytoplasm to the vacuole [19]. Regarding metabolism phase I reactions (e.g., hydroxylation, methylation), the presence of precursor ion of the parent drug in MS2 spectra of metabolites is less common due to fragmentation rules, which are dependent on molecular structure. Despite this, the product ions of parent substances are commonly shared with drug-related metabolites.

This observation aligns with a study [18] that found that diclofenac metabolites tend to cleave off parts of the molecule attached during conjugation metabolism, even at low collision energies. This suggests that the ion is fragmented to yield characteristic fragments while retaining information about the intact precursor. Similar fragmentation patterns were also noted in other studies [11, 18, 20, 36, 53, 55, 75] for different compounds such as diclofenac, trimethoprim, ibuprofen, clarithromycin, naproxen, and mefenamic acid. From these findings in both ionization modes, specific fragmentation patterns emerge, such as Δm/z 162.0528 for conjugation with glucose [