Discovering potential biomarkers for Ochratoxin A production by Penicillium nordicum in dry-cured meat matrices through untargeted metabolomics

Food Control 161 (2024) 110343 Available online 30 January 2024 0956-7135/© 2024 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/). Discovering potential biomarkers for Ochratoxin A production by Penicillium nordicum in dry-cured meat matrices through untargeted metabolomics David Garrido-Rodríguez a , María J. Andrade a , Josu´ e Delgado a , * , Eva Cebri´ an a , Inmaculada Barranco-Chamorro b a Higiene y Seguridad Alimentaria, Instituto Universitario de Investigaci´ on de Carne y Productos C´ arnicos, Facultad de Veterinaria, Universidad de Extremadura, Avda. de Las Ciencias S/n, 10003, C´ aceres, Spain b Departamento de Estadística e Investigaci´ on Operativa, Facultad de Matem´ aticas, Universidad de Sevilla, 41012, Sevilla, Spain ARTICLE INFO Keywords: Ochratoxin A (OTA) Untargeted metabolomics mixOmics MetaboAnalyst Penicillium nordicum Ripened meat product sPLS-DA ABSTRACT Ochratoxin A (OTA) is a mycotoxin produced by Penicillium nordicum, being this species a notable producer in dry-cured meat products. This toxin is one of the main fungal contaminants found in a variety of foods, including ripened meat products. Due to its harmful effects on health, it is essential to develop strategies from different approaches to identify whether products are contaminated and try to predict mycotoxin production during maturation to establish necessary corrective measures. In this study, targeted and untargeted metabolomic analyses were conducted to identify differential metabolites produced by three P. nordicum strains. These metabolites could be related to OTA production in agar-ham and dry-fermented sausage agar media under different water activity conditions. A Q Exactive™ Plus Hybrid Quadrupole-Orbitrap™ Mass Spectrometer was used to identify potential biomarkers and predictor-variables associated with OTA, employing: six comparatives that covered various aspects before and after OTA emergence, two analysis platforms (RStudio and MetaboAnalyst) and two approaches (imputation and non-imputation). The analysis was performed following two strategies, a Univariate and a Multivariate Statistical Analysis. As a result, several significant metabolites were identified as significantly linked to OTA production, including Asenjonamide C, (S,S)-Anacine, SPF-32629 A and Pseudouridine, which could potentially serve as biomarkers for OTA production. The study also highlighted the importance of using bioinformatic platforms, such as RStudio and MetaboAnalyst, for comprehensive analysis when integrating multiple outputs. This study paves the way for OTA prevention and taking corrective actions within the Hazard Analysis and Critical Control Points (HACCP) framework, after pointing out metabolites associated with OTA production, which deserve to be further studied in hams undergone to industrial ripening. 1. Introduction Ochratoxin A (OTA) is a pentaketyde with a molecular weight of 438.7 g/mol mainly produced by molds of the genera Penicillium and Aspergillus, among which Penicillium nordicum can be highlighted on drycured meat products. This toxin, firstly described in 1965 (Van der Merwe et al., 1965), is one of the main fungal metabolite contaminant in a great diversity of foods, including dry-cured meat products, where the high salt and protein concentration create the perfect substrate for this microorganism proliferation (Wang et al., 2016). Furthermore, it has been observed that its growth ideal conditions are at temperatures higher than 10 ◦C and a water activithe meat products have notty (a w ) around 0.90 (Rodríguez et al., 2015), conditions that take place during maturation of dry-cured meat products. About its toxic potential, the International Agency of Research on Cancer (IARC) classified the OTA as a possible carcinogen (2 B) in 1993, being nephrotoxic, hepatotoxic, neurotoxic, genotoxic, etc. The harmful effects on health of this toxin can be explained by its toxicokinetic, summarized on the scientific opinion published by the European Food Safety Authority (EFSA) (Schrenk et al., 2020). Due to its harmful effects on human health, the first regulation of OTA was established in 1990. In 2006, the EFSA set the maximum * Corresponding author. E-mail address: [email protected] (J. Delgado). Contents lists available at ScienceDirect Food Control journal homepage: www.elsevier.com/locate/foodcont https://doi.org/10.1016/j.foodcont.2024.110343 Received 13 November 2023; Received in revised form 7 January 2024; Accepted 24 January 2024 Food Control 161 (2024) 110343 2 tolerable intake of OTA at 14 ng/kg body weight per day, although other institutions, such as the U.S. Food and Drug Administration (FDA), have an even lower limit of 5 ng/kg body weight per day. The most up-to-date legislation is the Commission Regulation (EU) 2023/915 of April 25, 2023 on maximum levels for certain contaminants in food and repealing Regulation (EC) No 1881/2006 (European Commission, 2023), although the meat products have not been included with regard to OTA yet. Despite all this, nowadays the only legislation that establishes permitted maximum levels of OTA in meat products is Circolare 9 giugno 1999 No. 10 from Italy, setting a limit of 1 μ g/kg in pork meat (Ministero della Sanit` a, 1999). Because of all the aforementioned, it is essential to develop strategies, from different approaches, that allow the identification of contaminated products and, what is more interesting for the meat industry, to foresee the production of the mycotoxin during maturation in order to establish the necessary corrective measures, when required. While there are several available methods to identify products contaminated with OTA or detect OTA-producing strains, as described by Eshelli et al. (2018), there is a lack of studies focusing on metabolomics to gain a better understanding of its metabolism and predict its production in meat products. In the field of food microbiology, metabolomic studies serve as a powerful tool to evaluate various aspects related to food safety. These studies contribute to identify biomarkers of contaminants, such as the mentioned by Gallart-Ayala et al. (2015). These approaches also aid in identifying compounds produced by microorganisms, as highlighted by He et al. (2021). Additionally, metabolomics can be utilized for food authentication, quality characterization and other related purposes (Li et al., 2021). Through metabolomic analysis, primary and secondary metabolites present in a biological sample can be identified, leading to the generation of a metabolic profile. This profile reveals the metabolites that are either depleted or overproduced in response to various environmental, genetic, pathological and developmental conditions. Consequently, specific metabolic profiles can be associated with the production of toxins or other substances of interest. In this study, an untargeted metabolomic analysis was carried out to identify the differential metabolites produced by three P. nordicum strains, which could be associated with their OTA production in agarham and dry-fermented sausage agar (agar-dfs) matrices under different a w conditions (Supplementary Figs. 1A & B). The targeted and untargeted metabolomic results were analyzed by a Univariate Statistical Analysis (USA) and a Multivariate Statistical Analysis (MSA), being both performed by two approaches: Imputed, where all the metabolites detected were analyzed; and Common, where only the metabolites shared by both compared conditions for each comparative were analyzed to unveil those metabolites linked to OTA, as biomarkers or predictor-variables. 2. Materials and methods 2.1. Penicillium nordicum strains and culturing conditions Three P. nordicum strains were used: P. nordicum BFE856 (Pn856), from the Federal Research Centre of Nutrition and Food (Germany); P. nordicum PN15 (Pn15), from Food Hygiene and Safety Group at the University of Extremadura (Spain); and P. nordicum CBS 323.92 (Pn92), from Centraalbureau voor Schimmelcultures (The Netherlands). The use of behavioral different strains allows to extract global information considering multitude of ochratoxigenic strains, which enhance robustness of the results. These strains were inoculated and grown under 4 different culturing conditions: 1) agar-ham 0.94 a w (17 % w/w glycerol (Fisher Bio- Reagents™), 3 % w/w freeze-dried ham (purchased from a local establishment) and 2 % w/w European Bacteriological Agar (Condalab S.A.)), incubated at 20 ◦C. 2) agar-ham 0.90 a w (25 % w/w glycerol, 3 % w/w freeze-dried ham, 2 % w/w European Bacteriological Agar) incubated at 20 ◦C. 3) agar-dfs 0.94 a w (17 % w/w glycerol, 3 % w/w freeze-dried “salchich´ on” (purchased from a local establishment), 2 % w/w European Bacteriological Agar) incubated at 15 ◦C. 4) agar-dfs 0.90 a w (25 % w/w glycerol, 3 % w/w freeze-dried “salchich´ on”, 2 % w/w European Bacteriological Agar) incubated at 15 ◦C (Supplementary Fig. 1A). All the components were diluted using distilled water (Elix Technology Inside Milli-Q® Integral 5 Water System). The aim was to simulate the ones that take place during the drying stage of dry-cured meat products. 2.2. Sample preparation and collection P. nordicum strains were inoculated on yeast extract sucrose agar medium (Condalab S.A.) and incubated for 10 days at 25 ◦C. Spores were extracted with phosphate saline buffer and glass rod to be counted in a Thoma chamber (BlauBrand®). Spores were then inoculated by depositing 3 points of 4 μ L/point of a 10 6 spores/mL solution on the culture media and incubated under the conditions mentioned in Section 2.1. in culture incubators (FOC 225 E, Refrigerated Thermostat, VELP® Scientifica). Incubation lasted about 15 days, collecting colony samples from the first day of observable growth (c.a. ≥3 mm diameter) for 12 days. Three biological replicates were collected for each sampling, being weighted in a precision balance Shimadzu AUX220 and stored at −20 ◦C. At the end of sample collection, each collected mycelium was treated with acetonitrile Optima® HPLC (High Performance Liquid Chromatography) quality (Fisher Chemical) 1:5 (w/v). Then, samples were incubated at 200 rpm on the Unimax 2010 Heidolph® for an hour in the dark. Supernatants were collected after 10 min of centrifugation at 10.000 g and kept at −20 ◦C. The day when metabolomic analysis was performed, samples stored at −20 ◦C were diluted at 50 % with Milli-Q water, to enhance the chromatographic resolution, and set in the autosampler to be injected into the HPLC system of the Q Exactive™ Plus Hybrid Quadrupole-Orbitrap™ Mass Spectrometer (Supplementary Figs. 1A and 2A). 2.3. Q Exactive™ Plus Hybrid Quadrupole-Orbitrap™ mass spectrometer & software The metabolomic analysis, prior to the bioinformatics treatment, was conducted using the Q Exactive™ Plus Hybrid Quadrupole-Orbitrap™ Mass Spectrometer (MS; hereinafter referred to as Orbitrap-MS) (Supplementary Fig. 2A) in conjunction with the Free Style® and Compound Discover® software, for targeted and untargeted metabolomic assays, respectively. Both software were developed by Thermo Fisher Scientific. 2.3.1. OTA-targeted metabolomic analysis A C18 reverse-phase lipophilic chromatographic column (Thermo Fisher Accucore Aq C18 (150 ×2.1 mm) column 2.6 μ m particle size) was used in conjunction with the Free Style® software. The analysis was performed following the previously described by Cebri´ an et al. (2020). A 10-point calibration curve spanning concentrations from 0 to 90 ng/mL was injected and analyzed. The parameters used for it and to perform the targeted analyses can be found in Supplementary Table 1. The Limit of Detection (LOD) and Limit of Quantification (LOQ) for OTA were 0.0012 ng/mL and 0.01 ng/mL, respectively. The OTA targeted metabolomic results were summarized in an .xlsx file as Area Under the Curve (AUC) for each sample and replicate. These results were used in the conception and design of the untargeted analysis and deposited in the following GitHub repository: https://github. com/11DavidGR11/Metabolomic_analysis_from_AUC_datasets.git, specifically in folder “1_Data\“. Once these results were analyzed (Section 3.1.), several samples were selected to be processed by the untargeted approach to try to detect those metabolites linked to OTA production. Different time points were chosen to compare OTA-production with non-OTA production over time. D. Garrido-Rodríguez et al. Food Control 161 (2024) 110343 3 2.3.2. Untargeted metabolomics analysis Based on the OTA production, four comparatives were initially developed: comparative A to D (Supplementary Fig. 1B; Table 1). After obtaining untargeted datasets two additional comparatives (E and G) were designed (Supplementary Fig. 1B; Table 1). Four different types of samples were also used (Supplementary Table 1). The analyses performed with the Orbitrap-MS followed the parameters reported in the Supplementary Table 1. The Compound Discover® software followed the workflow given in Supplementary Fig. 1C. The process involves generating a mass trace for each ionized and detected metabolite, aligning RT with quality control samples, annotating identified compounds through various databases (ChemSpider, mzCloud, Natural Products Atlas 2020_06, Endogenous metabolites database with 4400 compounds), and distinguishing background compounds. As a result of this process an .xlsx was created for each condition of the four original comparatives and deposited in folder “1_Data\Untargeted _metabolomic_data\Individual_conditions” of the GitHub repository mentioned in Section 2.3.1. The data sets for comparatives E and G were obtained by duplicating their respective source data. 2.4. Bioinformatic analyses of metabolomic data The bioinformatic analysis was performed using two analytical platforms: RStudio, where we used the “Metabolomic_Analysis.Rmd” script, and the MetaboAnalyst.ca software. For each platform and for every comparison, we adopted two distinct approaches. The first approach involved analyzing all the detected metabolites for each condition in each comparison, irrespective of whether they were exclusively found in one condition (Imputed approach). The second approach aimed to streamline the study by excluding metabolites not detected in both conditions of a given comparison (Common approach). We refrained from conducting a global comparative analysis using these approaches due to the limited number of metabolites common to all conditions. In conclusion, we carried out two different statistical analyses, including Univariate and Multivariate analyses (Supplementary Fig. 2C). 2.4.1. OTA-targeted metabolomic analysis After splitting the raw dataset, of AUC values (three biological replicates per measured day), into different sheets based on the strain and growth conditions (a total of 12), the calibration curve was used to convert the units to ng/mL (Fig. 1A). Since the lowest point in the standard curve had an AUC intensity of 8.04 ×10 3 , corresponding to 0.01 ng/mL of OTA, values of AUC lowest than 6 ×10 3 were fixed as 0 ng/mL. Moreover, AUC values higher than the final point of the calibration curve (3.41 ×10 7 ) were eliminated from the dataset. Later, results were plotted in histograms (Fig. 1B), and the comparative design was performed based on the presence or absence of OTA in the samples. For this process, the “Targeted_metabolomic_data.Rmd” script, available in folder “1_Scripts\” of the GitHub repository previously mentioned, was used. 2.4.2. Raw data pre-processing by using RStudio Simple comparatives raw data were pre-processed with the script “Prepare_Datasets_to_Metabolomic_Analysis. Rmd”, deposited in folder “1_Scripts\” of the GitHub repository mentioned, and file “Correspondences.xlsx” in folder “1_Data\” (Supplementary Fig. 1D). With these files, AUC values of duplicates metabolites were aggregated, leaving a single row for each metabolite; unidentified compounds were renamed with their assigned chemical formula or their m/z and retention time (RT) values; row group information was added; and both conditions of each comparative were merged into a single dataset according to the compound name. For the initial approach, non-common metabolites (NA values) were imputed to include them into the analysis and identify their relationship with OTA production. However, in the second approach, they were removed from the dataset to study the differences just in common metabolites profiles. Final datasets were saved as .csv files in folder “1_Data\Untargeted_metabolomic_data\*ready_to_use”. Comparative G pre-processing script is in folder “1_Scripts\Comparative_G_scripts” and its correspondence file can be found in “1_Data \Comparative_G”. In this case, raw data files were treated as previous comparison, been merged all of them in a single final dataset with groups noOTA/OTA. Metabolites that were not found in all samples of a condition/group were removed from the data set, as we were only interested in metabolites that could be detected in all samples. Final dataset was saved as .csv file in folder “1_Data\Comparative_G\Untargeted_metabolomic_data\Ready_to_use”. 2.4.3. Metabolomic analysis by using RStudio The pre-processed data were analyzed with the script “Metabolomic_Analysis.Rmd” deposited in folder “1_Scripts\” of the GitHub repository. To perform this analysis the workflow provided in Supplementary Fig. 1D was followed. First, AUC values were normalized with the “NormalyzerDE” package (Willforss et al., 2019). The Quantile method was applied to single comparatives (Quantile-normalized.txt files obtained), whereas Loess method seemed more appropriate for comparative F (CycLoess-normalized.txt). The normalized data were analyzed by applying Principal Component Analysis (PCA) and Hierarchical Clustering on Principal Component (HCPC) based on Ward’s criterion from the “FactoMineR” (Lˆ e et al., 2008), and “factoextra” Table 1 Conditions compared on each comparative and the aim of each one. Condition 1 Condition 2 Aim Comparative A a Day 2 of growth of Pn15 at dryfermented sausage agar (agar-dfs) 0.90 water activity (a w ) (no future OTA production) Day 2 of Pn15 at agar-dfs 0.94 a w (future OTA production) Identify differential metabolites prior to OTA production Comparative B a Day 8 of growth of Pn15 at agar-dfs 0.90 a w (no detected OTA) Day 8 of Pn15 at agar-dfs 0.94 a w (detected OTA) Identify differential metabolites associated with OTA production and their relationship with a w Comparative C a Day 2 of growth of Pn15 at agar-ham 0.90 and 0.94 a w (no detected OTA) Day 10 of Pn15 at agar-ham 0.90 and 0.94 a w (detected OTA) Identify differential metabolites associated with OTA production, regardless of the a w levels Comparative D a Day 11 of growth of the three strains (Pn856, Pn15 and Pn92) at agar-dfs 0.90 a w (no detected OTA) Day 11 of growth of the three strains at agarham 0.90 a w (detected OTA) Identify the matrix effect in the OTA production and the metabolites associated with it Comparative E b Day 2 of growth of Pn15 at agar-dfs 0.94 a w (no detected OTA) Day 8 of growth of Pn15 at agardfs 0.94 a w (detected OTA) Identify differential metabolites associated with OTA production along the growth time at 0.94 a w Comparative G b (Global OTA samples) All samples and replicates where OTA was not detected All samples and replicates where OTA was detected Examine the overall differences in metabolite profiles between these two groups, OTA and non-OTA detected, independent of growth conditions a Comparatives developed before the untargeted metabolomic analysis were performed with the Q Exactive™ Plus Hybrid Quadrupole-Orbitrap™ Mass Spectrometer, and the dataset were obtained. b Comparatives developed after the untargeted metabolomic analysis were performed with the Q Exactive™ Plus Hybrid Quadrupole-Orbitrap™ Mass Spectrometer, and the dataset were obtained. D. Garrido-Rodríguez et al. Food Control 161 (2024) 110343 4 package (Kassambara & Mundt, 2020) to check the right groups separation (Supplementary Fig. 3). The USA was performed with “Limma” package (Ritchie et al., 2015). Although this package was suited for linear models and testing differential expression in microarray datasets, its use in metabolomic data has been previously described by Yu (2023). This package uses a moderated t-test and applies empirical Bayes methods to fit linear models and evaluate differential expressions. To compare the expression level mean ( μ ) under two conditions, for every metabolite the hypothesis test is (Equation (1)): H0: μ 1 i= μ 2 ivs H1: μ 1 i∕= μ 2 i.(1) Since μ 2 i, μ 1 i>0,the test can be carried out in terms of Log 2 Fold Change (LFC) (Equation (2)): H0:log2 μ 2 i μ 1 i =0vs H1:log2 μ 2 i μ 1 i ∕= 0.(2) As results, p-values and LFC values were obtained, performing a subsequent adjustment of the p-values by different methods with the aim of controlling the overall type I error. The methods used were Bonferroni and Benjamini-Hochberg (“p.adjust” function from the stats package in R was used). To identify the metabolites with significant differences for each comparative analyzed, specific thresholds were set. Absolute value of LFC (|LFC|) greater than or equal to 2 and a significance level ( α ) of 0.01 were used. These criteria aim to increase the confidence of the results by identifying metabolites that exhibit significant changes in relative abundance of at least a factor of 4.0 in each comparative, which were represented by upset plot (“UpSetR”; Conway & Gehlenborg, 2020). These plots helped to visualize and explore the significant metabolites, providing support to choose the Benjamini-Hochberg as the adjusting p-values method. After that decision, Volcano plots (“ggplot2” from package collection “Tidyverse”; Wickham et al., 2019) were performed, along with tables where the 10 most significant compounds, for each condition in each comparative, are listed. On the other hand, a sparse Partial Least Squares Discriminant Analysis (sPLS-DA) was conducted as a MSA using the function “splsda” from “mixOmics” package (Rohart et al., 2017). sPLS-DA allows us to carry out variable selection in metabolites when predicting the response variable (presence or absence of OTA), coded as 1/0, respectively. Based on the challenges faced during the selection of the optimal number of variables and components, the analysis utilized 2 components and 25 variables as suggested by the PCA results. The results of the sPLS-DA analysis include loadings plots, individual samples plots, variable plots, and heatmaps. Also, the top 10 most important predictor-variables for each comparative were presented in tables. Finally, a sensitivity analysis of our results was performed by using the function “predict” from “mixOmics” package (Rohart et al., 2017). Similarly to Tarazona et al. (2023), a few samples from the original datasets were subjected to slight perturbations, and the prediction data Fig. 1. Ochratoxin A (OTA) targeted results. A) Calibration pattern curve performed with known concentrations of OTA samples (ng/mL). B) Error bar plots of the OTA concentration detected in each condition and Penicillium nordicum strain with the conditions used for each comparative. Each plot shows the amount of OTA produced (ng/g) per sampling day in each growth condition of the experiment. The P. nordicum strains used in the experiment were Pn15, Pn856 and Pn92, and the growth conditions were two culture medium (agarham (J), dry-fermented sausage-agar (agar-dfs; S)), and two water activity (a w ) 0.90 and 0.94). Magenta boxes: they are the conditions compared in comparative A, second sampling day of Pn15 at agar-dfs 0.90 a w , against second sampling day of Pn15 at agar-dfs 0.94 a w . Blue boxes: comparative B conditions, day 8 of sampling of Pn15 at agar-dfs 0.90 aw, against eighth sampling day of Pn15 at agar-dfs 0.94 a w . Red boxes: indicate the conditions used for comparative C design. Day 2 of sampling of Pn15 at agar-ham at both a w , versus day 10 of sampling of Pn15 at agar-ham at both aw. Orange boxes: correspond to the conditions compared in comparative D. Eleventh sampling day of the three strains at agar-dfs 0.90 a w , against eleventh sampling day of the three strains at agar-ham 0.90 a w . Filled purple boxes: they are the conditions compared in the comparative E. Second sampling day of Pn15 at agar-dfs 0.94 a w , against eighth sampling day of Pn15 at agar-dfs 0.94 a w . C) Error bar plots of OTA concentration detected in each condition and P. nordicum strain with the conditions not used in the comparatives. Each plot shows the amount of OTA produced per sampling day in each growth condition of the experiment, displaying sampling days on the X-axis and OTA concentration on the Y-axis. Created with BioRender.com. D. Garrido-Rodríguez et al. Food Control 161 (2024) 110343 5 were compared with the original data using the Mahalanobis distance (Reprint of: Mahalanobis (2018)). The results of these analyses are summarized in the Supplementary Table 8. 2.4.4. Metabolomic analysis by using MetaboAnalyst.ca As a comparative study, the preprocessed .csv files were analyzed using the Statistical Analysis [one factor] function of the web-based tool MetaboAnalyst.ca (Pang et al., 2021), following a workflow adapted from Herraiz GilBioInfoGRX (2022) and that can be referred in Supplementary Fig. 1E. The configuration used for data loading was peak intensities as data type and samples in columns (unpaired) as format. While no data filtering was performed, samples were normalized by sum and pareto scaling. Subsequently, the data were visualized using PCA samples plots, heatmaps displaying the top 25 metabolites, and a dendrogram representation based on the Euclidean distance measure and Ward clustering algorithm. USA consisted of a fold change analysis with a threshold value of 4.0 (|LFC| >2.0) and direction of comparison of group 2 (in most of them OTA was detected) versus group 1 (OTA was not detected). A two-sample t-tests was performed with an α =0.01 and the p-value adjustment method False Discovery Rate (FDR). With both results, a volcano plot was performed and the corresponding data file with significant metabolites was obtained. Again, the MSA method used in the analysis was sPLS-DA, using 2 components and 25 variables. To assess the performance of models, Leave One Out Cross Validation (LOOCV) was used. 3. Results and discussion 3.1. OTA production The OTA production by the three P. nordicum strains evaluated under each growth condition was highly heterogeneous (Fig. 1B and C). The use of multitude of ochratoxigenic strains, behavioral different among them, allows to extract global information, which improves the results’ robustness. This mycotoxin was detected in agar-ham at 0.94 a w for every P. nordicum strain, whereas no strain produced OTA in agar-dfs at 0.90 a w . For agar-dfs 0.94 a w , Pn15 was the only strain capable of producing OTA (Fig. 1B.3, 1C.2 and 1C.4). These results are not in agreement with those obtained by Rodríguez et al. (2015) and S´ anchez-Montero et al. (2019), who detected OTA production by P. nordicum in the range of temperature 10–30 ◦C and 0.87–0.99 a w . In both studies, higher OTA production was detected at the lowest temperatures and a w values within the referred ranges. These differences could be due to the fact that in our study, glycerol was used to modulate the a w , which could be used as a substrate to resist stress (Magan, 2007), while the two mentioned studies used common salt (NaCl) for this purpose. In this regard, NaCl has been postulated as a potential inducer of OTA production by P. nordicum through an osmotic compensation mechanism in environments with high salt concentrations (Delgado et al., 2018; Schmidt-Heydt et al., 2012), being a plausible explanation of these apparent discrepancies. When OTA production is compared among the three P. nordicum strains, Pn15 showed the highest values (Fig. 1B). It also produced OTA under the broadest range of conditions and exhibited the most anticipated distribution pattern of product accumulation over production time. In fact, OTA production exhibited a continuous increase over time, with little variation (Fig. 1B1, 1B3, and 1B4), in contrast to that observed for Pn92 (Fig. 1B5 and 1C3). For Pn92 the OTA production in agar-ham showed a highly heterogenous distribution which increases and decreases along the 12 sampling days, reaching the highest concentration on day 6 for agar-ham at 0.94 a w (Fig. 1C3), and on day 11 for agar-ham at 0.90 a w (Fig. 1B5). On the other hand, from the fifth day of sampling, OTA production by Pn856 in agar-ham 0.94 a w decreased until concentrations even lower than the first days (Fig. 1C1). Although at first it may seem puzzling, the literature led us to argue that it could be due to the degradation of the mycotoxin in the medium, a fact previously described by other authors (Bailly et al., 2005). Similar results were found in the studies with P. nordicum on dry-fermented sausages during their industrial ripening by Delgado et al. (2021). Unfortunately, due to the absence of untargeted metabolomic analysis conducted in the observed conditions of reduced OTA accumulation, the assessment of potential OTA degradation metabolites was not carried out, and as a result, this study could not confirm this phenomenon. Based on these findings, we employed several comparative and proposed some studies to identify the metabolites associated with OTA production under different growth conditions using untargeted metabolomic analysis. These comparatives are those mentioned in Section 2.3.2. And displayed in the Supplementary Fig. 1B. 3.2. Number of significant metabolites by the univariate statistical analysis Although the significant number of metabolites was too low in the Common approaches, the statistical thresholds (Sections 2.4.3. and 2.4.4.) were maintained to allow the comparison of results and the rigor of the analysis. As mentioned in Section 2.4.3., the hypothesis test was carried out by comparing group 2 (OTA detected group for all comparatives and group cA_S17 for comparative A) against group 1 (no-OTA detected group for all comparatives and group cA_S25 for comparative A) for each comparative. Using these thresholds, the number of significant metabolites was quite different for each platform, being usually detected a greater number of them in the RStudio analysis. E.g., in comparative A, 162 increased and 188 decreased metabolites were obtained using RStudio, whilst MetaboAnalyst resulted in 132 increased and 174 decreased metabolites. As expected, the number of significant metabolites detected by the Imputed approaches was notably higher than that obtained by the Common approaches. E.g., comparative C using RStudio, showed a total of 851 significant metabolites detected by the Imputed approach, in contrast of the 154 significant metabolites obtained by the Common one. The summary can be found in Table 2, while the metabolites identification is contained in Supplementary Tables 2, 3, 4, 5 and 6. As it can be seen in the upset plots of Fig. 2A, where significant metabolites obtained for each simple comparative were contrasted along the four different analyses followed (both platforms and approaches), most of metabolites classified as significant by the Common approaches are included in the significant metabolites obtained by the Imputed approaches. Comparative A is the only one in which all common significant metabolites have also been detected as significant in the Imputed approaches (Fig. 2A). Whereas, in the other comparisons, there is at least one metabolite that has only been detected as significant in Common approaches in one or both platforms used for the analysis. An example of this last situation is the (S,S)-Anacine from comparative B (Supplementary Table 3). These differences in metabolites that are determined to be significant in a Common approach but not in an Imputed one could be due to the pvalue adjustment method. With Benjamini-Hochberg method, when we have smaller number of metabolites (test), as in Common approaches (Table 2), the p-value adjustment may be too conservative and reduce the statistical power of the analysis, which may make difficult to detect truly significant relationships (Pounds & Cheng, 2005). However, the advantages of this method with large number of tests, as in Imputed approaches, make it the ideal method for this study. We are not only trying to identify significant metabolites linked to OTA, but also trying to make a comparison between the results obtained from different approaches. Fig. 2A displays some compounds classified as significant metabolites by all the analysis performed on both platforms (RStudio and MetaboAnalyst) and approaches (Imputed and Common). In comparative A there are 7 shared significant metabolites (Supplementary Table 2), comparative B shares 9 significant metabolites (Supplementary Table D. Garrido-Rodríguez et al. Food Control 161 (2024) 110343 6 Table 2 Number of significant metabolites detected in each condition and analytical approach, based on the Benjamini-Hochberg p-value adjustment method. The statistical threshold used were significance level ( α ) of 0.01 and absolute value of the Log 2 Fold Change (|LFC|) greater or equal to 2.0. As the hypothesis test was H0:log2 μ 2 i μ 1 i =0|H1:log2 μ 2 i μ 1 i ∕= 0, results are expressed according to the second group of each comparative. D. Garrido-Rodríguez et al. Food Control 161 (2024) 110343 7 3), comparative C shares 87 significant metabolites (Supplementary Table 4), comparative D shares 33 significant metabolites (Supplementary Table 5), and comparative E shares 19 significant metabolites (Supplementary Table 6). Since this study introduces a novel approach by employing untargeted metabolomics for analyzing P. nordicum-derived metabolites, the novelty of this study makes impossible the direct comparisons with previous studies exploring this ochratoxigenic mold in similar terms. Instead, comparisons were made with other metabolomic studies using different techniques. The results obtained using the Orbitrap-MS system and Compound Discoverer software yielded a much higher total number of metabolites detected compared to similar studies. In the present study between 534 and 1,560 different metabolites were detected in each comparative analyzed by the Imputed approach. For instance, a semi-quantitative UHPLC-MS/MS analysis in the study of mycotoxin-associated microbiome in nuts resulted in 191 metabolites (Varga et al., 2013). Additionally, a study on the metabolome associated with Aspergillus flavus using an HPLC-MS/QTOF system identified 14 potential biomarker metabolites for aflatoxin production (Saldan et al., 2018). While in our study, more than 800 significantly produced metabolites were detected in some comparatives (Table 2). 3.3. Multivariate Statistical Analysis results and main predictor-variables Hypothesis tests are useful to identify differentially behavioral metabolites and may have an explanatory or exploratory purpose. Following this analytical strategy, metabolites with statistically significant differences between the conditions compared in each comparative were considered as possible biomarkers (Supplementary Fig. 1B; Fig. 2B and C). However, in this context, we have many variables (metabolites) highly correlated, which is a serious problem for traditional statistical techniques. The high correlation between the explanatory variables demands the use of multivariate techniques to provide a more complete comprehension of biochemical processes. Since metabolites are continuous variables and the variable response is 0/1, sparse PLS-DA (sPLS-DA) is appropriate. Using this technique to predict the presence or absence of OTA allows us to carry out a selection of variables involved in the process (Figs. 3 and 4). From the initial number of compounds in each comparative (Fig. 3A), 25 were selected based on their high correlation with one sample group and very low correlation with the opposite group (Fig. 3B–F). This technique was used in previous studies (Ruiz-Perez et al., 2020; Todorov et al., 2020) in different fields but with similar aims, getting as conclusions that this analytical technique offers high sensitivity in categorical classifications. Based on that previously mentioned, it was expected to obtain some possible predictor-variables (metabolites) different from those detected as significant by the USA. The hypothesis tested between these groups is referenced in Supplementary Fig. 4. The 25 possible predictor-variables of the Imputed approach of each comparative, are all included in the statistically significant metabolites. Whereas, in the Common approaches there are some predictor-variables which had not been previously detected as statistically significant. In comparatives A and B these differences could be due to the fact that we had selected 25 variables to perform the sPLS-DA, whilst the number of statistically significant metabolites is smaller. However, in comparative D (Supplementary Fig. 4), Fig. 2. Plots of the Univariate Statistical Analysis results. A) Upset-plots comparing the significant metabolites obtained by each study. These studies were performed in two platforms (analysis on RStudio and with MetaboAnalyst.ca software), and two approaches (Imputed and Common). The metabolites that will be analyzed are those that have been detected in both condition but in different concentrations. B) Volcano-plots of all significant metabolites, including those that were only detected in one of the compared conditions. These plots are made by plotting the LFC value, of each detected metabolite, on X-axis; and the Log 10 of its p-values adjusted by Benjamini-Hochberg on Y-axis. The conditions compared to perform the analysis and graphs are shown in each individual main title plot. Those metabolites that appear in orange are increased in the first group indicated in the plot-main, while the blue ones are increased in the second group. C) Volcano-plots of common significant metabolites, deleted those that were not detected in the two compared conditions. For each individual plot-main the conditions compared to perform the analysis and graphs are shown. Metabolites shown in orange are increased in the first group indicated, while the blue ones are increased in the second group. Created with BioRender.com. D. Garrido-Rodríguez et al. Food Control 161 (2024) 110343 8 the Cytisine has been detected as a possible predictor-variable only in the Common approach. Moreover, 2 ′ -Deoxyinosine, Leucyl-Proline, Maniwamycin C, and Purine are shared in such comparative as significant metabolites among the possible predictor-variables of the Common approach and the Imputed approach of the USA (Supplementary Fig. 4). 3.4. OTA as a result of untargeted metabolomic studies The major aim of this study was to identify metabolites linked to OTA and its production. Therefore, as mentioned in Sections 2.3.1. and 2.3.2., OTA-targeted metabolomic analysis was followed by untargeted metabolomic analysis. Firstly, OTA was measured to determine if it was present or absent in the samples (Fig. 1), and then comparatives were performed and analyzed by untargeted metabolomics. Based on this, OTA should be detected as a significant ‘potential predictor-variable’ metabolite in all comparatives where it was targeted-detected in one condition and no targeted-detected in the other one. 3.4.1. OTA as a differential metabolite by the univariate statistical analysis When referring to untargeted metabolomics, as it can be seen in the volcano plots of Fig. 2B and C, OTA was detected as significant metabolite only in comparative D, where it was identified as significant metabolite by both approaches (Imputed and Common). This contradiction can be explained by the mass spectrometer full scan configuration, without SIM operational mode, that generated these results. In this mode, the less abundant metabolites can be masked by the most abundant ones owing to the limited dynamic range of the MS. This is one of the reasons for using a combination of targeted and untargeted metabolomic approaches in this study. Moreover, in comparative C and D, ochratoxin B (OTB), an intermediate metabolite of OTA biosynthesis (Wang et al., 2018), was identified as significant metabolite (Fig. 2B and C). Regarding the OTB, it was only detected as significant in comparative D by using RStudio (Supplementary Table 5). The disparities observed between the results obtained in each analysis platform (RStudio and MetaboAnalyst) may be attributed to the normalization step, conducted prior to the statistical analysis (quantile versus pareto scaling), or the statistical test (moderated t-test versus two-sample t-test). From the observations made, quantile normalization (performed in RStudio) seems to be the best method in this case, as its posterior USA successfully identified OTB as a statistically significant metabolite. This outcome was expected, given that OTB is a precursor in the biosynthesis of OTA (Wang et al., 2018). 3.4.2. OTA among the main biomarkers obtained by the Multivariate Statistical Analysis Tables 3 and 4 summarize the 10 principal metabolites classified as predictor-variables of each condition in each comparative. The first one includes metabolites obtained by the Imputed approach, whilst the second one includes those obtained by the Common approach. As can be seen, most of them are not shared along both analyses. Thus, to reinforce the predictive character of the possible biomarkers, the results will be focused on the most significant of the shared ones. In Fig. 4 the most important predictor-variables obtained by the RStudio analysis for each comparative and approach are represented. OTA has not been detected as a statistically significant compound, being only detected in comparative D by the Common approach (Fig. 4D–Table 4). These results suggest that the analysis conducted with the “Metabolomic_Analysis.Rmd” script in RStudio is superior to the one performed on MetaboAnalyst. Alternatively, parameters used in the MetaboAnalyst study should be reconsidered. 3.5. Most outstanding metabolites of the achieved results Metabolites detected as statistically significant by the Imputed approach can be directly considered differential metabolites in the production of OTA, with a likely biomarker character of its production or absence. It is crucial to point out that the metabolite identification, based on exact mass matches and pattern of MS2 (tandem mass spectrometry) fragmentation, when available, allows in a relatively accurate way the identification of every metabolite. Nonetheless, it depends on the library availability and a final confirmation by magnetic resonance would be the further step to definitively confirm every potential biomarker tentatively identified in this study. Although their RT, under the same chromatographic conditions and m/z with the MS2 pattern when available, could be used with the purpose of biomarker detection with no additional identifying requirements. Additionally, it is important to perform a properly background noise deletion; if not, differences in the used culture conditions could be mistaken with really differences. It is the case of the Glycerol, the metabolite with highest significance differences in comparative A (Fig. 2A). Though it can be considered a strong biomarker for the OTA absence, it is important to note that it is increased at day 2 of sampling for Pn15 at agar-dfs 0.90 a w (Supplementary Fig. 1B), which was designed with a 25 % of glycerol, against the 17 % used in the compared condition. These results can be also observed by the sPLS-DA test (Fig. 4A; Table 3). Fig. 3. Heatmaps showing the correlation between samples and the 25 principal metabolites classified as predictor-variables by the Multivariate Statistical Analysis with sparse Partial Least Square Discriminant Analysis (sPLS-DA). A) Example of a heatmap showing the correlation between all metabolites in comparative A). The metabolites with a higher or lower correlation with the samples are distributed on the sides of the graph; whilst those with no statistically significant differences are in the center of the graph. B) Heatmap for the analysis of comparative A. C) Heatmap for the analysis of comparative B. D) Heatmap for the analysis of comparative C. E) Heatmap for the analysis of comparative D. F) Heatmap for the analysis of comparative E. Use Supplementary Table 9 to decode the IDs of the metabolites displayed in each heatmap. D. Garrido-Rodríguez et al. Food Control 161 (2024) 110343 9 (caption on next page) D. Garrido-Rodríguez et al. Food Control 161 (2024) 110343 16 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. Data availability Data will be made available on request. Acknowledgments R. Geisen is acknowledged for the P. nordicum BFE856 donation. Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi. org/10.1016/j.foodcont.2024.110343. References Abdelkader, M., Philippon, T., Asenjo, J., Bull, A., Goodfellow, M., Ebel, R., Jaspars, M., & Rateb, M. (2018). Asenjonamides A–C, antibacterial metabolites isolated from Streptomyces asenjonii strain KNN 42.f from an extreme-hyper arid Atacama Desert soil. Journal of Antibiotics, 71. https://doi.org/10.1038/s41429-017-0012-0 Bailly, J. D., Tabuc, C., Qu´ erin, A., & Guerre, P. (2005). Production and stability of patulin, ochratoxin A, citrinin, and cyclopiazonic acid on dry cured ham. Journal of Food Protection, 68(7), 1516–1520. https://doi.org/10.4315/0362-028x-68.7.1516 Boyes-Korkis, J. M., Gurney, K. A., Penn, J., Mantle, P. G., Bilton, J. N., & Sheppard, R. N. (1993). Anacine, a new benzodiazepine metabolite of Penicillium aurantiogriseum produced with other alkaloids in submerged fermentation. Journal of Natural Products, 56(10), 1707–1717. https://doi.org/10.1021/np50100a008 Cebri´ an, E., Nú˜ nez, F., G´ alvez, F. J., Delgado, J., Bermúdez, E., & Rodríguez, M. (2020). Selection and evaluation of Staphylococcus xylosus as a biocontrol agent against toxigenic moulds in a dry-cured ham model system. Microorganisms, 8(6). https:// doi.org/10.3390/microorganisms8060793 Conway, J., & Gehlenborg, N. (2020). UpSetR: A more scalable alternative to venn and euler diagrams for visualizing intersecting sets. (1.4.0). Delgado, J., da Cruz Cabral, L., Rodríguez, M., & Rodríguez, A. (2018). Influence of ochratoxin A on adaptation of Penicillium nordicum on a NaCl-rich dry-cured hambased medium. International Journal of Food Microbiology, 272, 22–28. https://doi. org/10.1016/j.ijfoodmicro.2018.02.020 Delgado, J., Rond´ an, J. J., Nú˜ nez, F., & Rodríguez, A. (2021). Influence of an industrial dry-fermented sausage processing on ochratoxin A production by Penicillium nordicum. International Journal of Food Microbiology, 339, Article 109016. https:// doi.org/10.1016/j.ijfoodmicro.2020.109016 Eshelli, M., Qader, M., Jambi, E., Hursthouse, A., & Rateb, M. (2018). Current status and future opportunities of omics tools in mycotoxin research. Toxins, 10(11), 433. https://doi.org/10.3390/toxins10110433 European Commission. (2023). Commission Regulation (EU) 2023/915 of 25 April 2023 on maximum levels for certain contaminants in food. Official Journal of the European Union. Gallart-Ayala, H., Ch´ ereau, S., Dervilly-Pinel, G., & Bizec, B. Le (2015). Potential of mass spectrometry metabolomics for chemical food safety. Bioanalysis, 7(1), 133–146. https://doi.org/10.4155/bio.14.267 Gil, H., & BioInfoGRX, S. (2022). Mini-Tutorial: An´ alisis de Datos Metabol´ omicos con MetaboAnalyst 5.0. YouTube https://youtu.be/fgWbj0IokjI.(2023-06-01). He, T., Yuan, C., & Zhao, C. (2021). Long intragenic non-coding RNA p53-induced transcript (LINC-PINT) as a novel prognosis indicator and therapeutic target in cancer. Biomedicine & Pharmacotherapy, 143, Article 112127. https://doi.org/ 10.1016/j.biopha.2021.112127 Kassambara, A., & Mundt, F. (2020). Factoextra : Extract and visualize the results of multivariate data analyses (1.0.7). Larsen, T. O., Svendsen, A., & Smedsgaard, J. (2001). Biochemical characterization of ochratoxin A-producing strains of the genus Penicillium. Applied and Environmental Microbiology, 67(8), 3630–3635. https://doi.org/10.1128/AEM.67.8.3630- 3635.2001 Lˆ e, S., Josse, J., & Husson, F. (2008). FactoMineR: An R package for multivariate analysis. Journal of Statistical Software, 25(1 SE-Articles), 1–18. https://doi.org/ 10.18637/jss.v025.i01 Li, S., Tian, Y., Jiang, P., Lin, Y., Liu, X., & Yang, H. (2021). Recent advances in the application of metabolomics for food safety control and food quality analyses. Critical Reviews in Food Science and Nutrition, 61(9), 1448–1469. https://doi.org/ 10.1080/10408398.2020.1761287 Magan, N. (2007). Fungi in extreme environments. Mycotaxon, 4, 85–103. Mahalanobis, P. C. (2018). Reprint of: Mahalanobis, P.C. (1936) “on the generalised distance in statistics.”. Sankhya - Series A, 80(1), 1–7. https://doi.org/10.1007/ s13171-019-00164-5 Ministero della Sanit` a. (1999). Direttive in materia di controllo ufficiale sui prodotti alimentari: valori massimi ammissibili di micotossine nelle derrate alimentari di origine nazionale, comunitaria e paesi terzi. Gazz. Uff. della Repubb. Ital., 135, 52–57. Nú˜ nez, F., Westphal, C. D., Bermúdez, E., & Asensio, M. A. (2007). Production of secondary metabolites by some terverticillate penicillia on carbohydrate-rich and meat substrates. Journal of Food Protection, 70(12), 2829–2836. https://doi.org/ 10.4315/0362-028X-70.12.2829 Pang, Z., Chong, J., Zhou, G., de Lima Morais, D. A., Chang, L., Barrette, M., Gauthier, C., Jacques, P.-´ E., Li, S., & Xia, J. (2021). MetaboAnalyst 5.0: Narrowing the gap between raw spectra and functional insights. Nucleic Acids Research, 49(W1), W388–W396. https://doi.org/10.1093/nar/gkab382 Pounds, S., & Cheng, C. (2005). Sample size determination for the false discovery rate. Bioinformatics, 21(23), 4263–4271. https://doi.org/10.1093/bioinformatics/bti699 Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., & Smyth, G. K. (2015). Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research, 43(7), e47. https://doi.org/10.1093/nar/gkv007. e47. Rodríguez, A., Capela, D., Medina, ´ A., C´ ordoba, J. J., & Magan, N. (2015). Relationship between ecophysiological factors, growth and ochratoxin A contamination of drycured sausage based matrices. International Journal of Food Microbiology, 194, 71–77. https://doi.org/10.1016/j.ijfoodmicro.2014.11.014 Rohart, F., Gautier, B., Singh, A., & Lˆ e Cao, K.-A. (2017). mixOmics: An R package for ’omics feature selection and multiple data integration. PLoS Computational Biology, 13(11), Article e1005752. https://doi.org/10.1371/journal.pcbi.1005752 Ruiz-Perez, D., Guan, H., Madhivanan, P., Mathee, K., & Narasimhan, G. (2020). So you think you can PLS-DA? BMC Bioinformatics, 21(1), 2. https://doi.org/10.1186/ s12859-019-3310-7 Saldan, N. C., Almeida, R. T. R., Avíncola, A., Porto, C., Galuch, M. B., Magon, T. F. S., Pilau, E. J., Svidzinski, T. I. E., & Oliveira, C. C. (2018). Development of an analytical method for identification of Aspergillus flavus based on chemical markers using HPLC-MS. Food Chemistry, 241, 113–121. https://doi.org/10.1016/j. foodchem.2017.08.065 S´ anchez-Montero, L., C´ ordoba, J. J., Peromingo, B., ´ Alvarez, M., & Nú˜ nez, F. (2019). Effects of environmental conditions and substrate on growth and ochratoxin A production by Penicillium verrucosum and Penicillium nordicum: Relative risk assessment of OTA in dry-cured meat products. Food Research International, 121, 604–611. https://doi.org/10.1016/j.foodres.2018.12.025 Schmidt-Heydt, M., Graf, E., Stoll, D., & Geisen, R. (2012). The biosynthesis of ochratoxin A by Penicillium as one mechanism for adaptation to NaCl rich foods. Food Microbiology, 29(2), 233–241. https://doi.org/10.1016/j.fm.2011.08.003 CONTAM (EFSA), Schrenk, D., Bodin, L., Chipman, J. K., del Mazo, J., Grasl-Kraupp, B., Hogstrand, C., Hoogenboom, L.(R.), Leblanc, J.-C., Nebbia, C. S., Nielsen, E., Ntzani, E., Petersen, A., Sand, S., Schwerdtle, T., Vleminckx, C., Wallace, H., Alexander, J., Dall’Asta, C., … Bignami, M. (2020). Risk assessment of ochratoxin A in food. EFSA Journal, 18(5), Article e06113. https://doi.org/10.2903/j. efsa.2020.6113 Shah, Z. A., Khan, K., Rashid, H. U., Shah, T., Jaremko, M., & Iqbal, Z. (2022). Insights into metabolic and pharmacological profiling of Aspergillus ficuum through bioinformatics and experimental techniques. BMC Microbiology, 22(1), 295. https:// doi.org/10.1186/s12866-022-02693-w Shimatani, T., Hosotani, N., Ohnishi, M., Kumagai, K., & Saji, I. (2006). SPF-32629 A and B, novel human chymase inhibitors produced by Penicillium sp. Journal of Antibiotics, 59(1), 29–34. https://doi.org/10.1038/ja.2006.5 Tarazona, S., Ugidos, M., Salguero, P., & Arzalluz, ´ A. (2023). Multivariate methods. Chapter sparse PLS-DA (sPLS-DA). Prediction. Todorov, H., Searle-White, E., & Gerber, S. (2020). Applying univariate vs. multivariate statistics to investigate therapeutic efficacy in (pre)clinical trials: A Monte Carlo simulation study on the example of a controlled preclinical neurotrauma trial. PLoS One, 15(3), Article e0230798. https://doi.org/10.1371/journal.pone.0230798 Van der Merwe, K. J., Steyn, P. S., Fourie, L., Scott, D. B., & Theron, J. J. (1965). Ochratoxin A, a toxic metabolite produced by Aspergillus ochraceus Wilh. Nature, 205 (976), 1112–1113. https://doi.org/10.1038/2051112a0 Varga, E., Glauner, T., Berthiller, F., Krska, R., Schuhmacher, R., & Sulyok, M. (2013). Development and validation of a (semi-)quantitative UHPLC-MS/MS method for the determination of 191 mycotoxins and other fungal metabolites in almonds, hazelnuts, peanuts and pistachios. Analytical and Bioanalytical Chemistry, 405(15), 5087–5104. https://doi.org/10.1007/s00216-013-6831-3 Wang, Y., Wang, L., Liu, F., Wang, Q., Selvaraj, J. N., Xing, F., Zhao, Y., & Liu, Y. (2016). Ochratoxin A producing fungi, biosynthetic pathway and regulatory mechanisms. Toxins, 8(3), 83. https://doi.org/10.3390/toxins8030083 Wang, Y., Wang, L., Wu, F., Liu, F., Wang, Q., Zhang, X., Selvaraj, J. N., Zhao, Y., Xing, F., Yin, W.-B., & Liu, Y. (2018). A consensus ochratoxin A biosynthetic pathway: Insights from the genome sequence of Aspergillus ochraceus and a comparative genomic analysis. Applied and Environmental Microbiology, 84(19), e01009–e01018. https://doi.org/10.1128/AEM.01009-18 D. Garrido-Rodríguez et al. Food Control 161 (2024) 110343 17 Wickham, H, Averick, M, Bryan, J, Chang, W, McGowan, LD, François, R, … Yutani, H (2019). Welcome to the tidyverse. Journal of Open Source Software, 4(43), 1686. https://doi.org/10.21105/joss.01686 Willforss, J., Chawade, A., & Levander, F. (2019). NormalyzerDE: Online tool for improved normalization of omics expression data and high-sensitivity differential expression analysis. Journal of Proteome Research, 18(2), 732–740. https://doi.org/ 10.1021/acs.jproteome.8b00523 Xia, K., He, X., Dai, Q., Cheng, W.-H., Qi, X., Guo, M., Luo, Y., Huang, K., Zhao, C., & Xu, W. (2014). Discovery of systematic responses and potential biomarkers induced by ochratoxin A using metabolomics. Food Additives & Contaminants Part A, Chemistry, Analysis, Control, Exposure & Risk Assessment, 31(11), 1904–1913. https:// doi.org/10.1080/19440049.2014.957249 Yu, M. (2023). Meta-workflow. Chapter 10: Statistical analysis. D. Garrido-Rodríguez et al.