<figure><img src="image_1.png" /></figure>

HAL
open science

# Modeling HepaRG metabolome responses to pyrrolizidine alkaloid exposure for insight into points of departure and modes of action

Estelle Dubreil, Keyvin Darney, Marie Laure Delignette-Muller, Audrey

Barranger, Sylvie Huet, Kevin Hogeveen, Thibaut Léger, Valérie Fessard,
Ludovic Le Hégarat

► To cite this version:

Estelle Dubreil, Keyvin Darney, Marie Laure Delignette-Muller, Audrey Barranger, Sylvie Huet, et al.. Modeling HepaRG metabolome responses to pyrrolizidine alkaloid exposure for insight into points of departure and modes of action. Journal of Hazardous Materials, 2024, 474, pp.134721. 10.1016/j.jhazmat.2024.134721. hal-04611387

HAL Id: hal-04611387

https://cnrs.hal.science/hal-04611387

Submitted on 17 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Modeling HepaRG metabolome responses to pyrrolizidine alkaloid exposure for insight into points of departure and modes of action

Estelle Dubreil, Keyvin Darney, Marie Laure Delignette-Muller, Audrey

Barranger, Sylvie Huet, Kevin Hogeveen, Thibaut Léger, Valérie Fessard,
Ludovic Le Hégarat

► To cite this version:

Estelle Dubreil, Keyvin Darney, Marie Laure Delignette-Muller, Audrey Barranger, Sylvie Huet, et al.. Modeling HepaRG metabolome responses to pyrrolizidine alkaloid exposure for insight into points of departure and modes of action. Journal of Hazardous Materials, 2024, 474, pp.134721. 10.1016/j.jhazmat.2024.134721. hal-04611387

HAL Id: hal-04611387

https://cnrs.hal.science/hal-04611387v1

Submitted on 17 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Modeling HepaRG metabolome responses to pyrrolizidine alkaloid exposure for insight into points of departure and modes of action

Estelle Dubreil<sup>a, *</sup>, Keyvin Darney<sup>b</sup>, Marie-Laure Delignette-Muller<sup>c</sup>, Audrey Barranger<sup>a</sup>, Sylvie Huet<sup>a</sup>, Kevin Hogeveen<sup>a</sup>, Thibaut Léger<sup>a</sup>, Valérie Fessard<sup>a</sup>, Ludovic Le Hégarat<sup>a</sup>

## HIGHLIGHTS

* PA effects on HepaRG cells were studied by concentration-response modelling and BMC.

* The endometabolome was found to be more sensitive than the exometabolome.

* LAS had the highest hepatotoxic potency, followed by HEL and RET.

* The bell shape is the most frequent trend in response to LAS.

* Lipids, bile acids and amino acids were especially affected by LAS.

## GRAPHICAL ABSTRACT

<figure><img src="image_4.png" /></figure>

## ARTICLE INFO

**Keywords:**

Alkaloid

Food contaminant

Dose-response model

Benchmark concentration

Metabolomics

## ABSTRACT

The new challenges in toxicology demand novel and innovative in vitro approaches for deriving points of departure (PODs) and determining the mode of action (MOA) of chemicals. Therefore, the aim of this original study was to couple in vitro studies with untargeted metabolomics to model the concentration-response of extra- and intracellular metabolome data on human HepaRG cells treated for 48 h with three pyrrolizidine alkaloids (PAs): heliotrine, retrorsine and lasiocarpine. Modeling revealed that the three PAs induced various monotonic and, importantly, biphasic curves of metabolite content. Based on unannotated metabolites, the endometabolome was more sensitive than the exometabolome in terms of metabolomic effects, and benchmark concentrations (BMCs) confirmed that lasiocarpine was the most hepatotoxic PA. Regarding its MOA, impairment of lipid metabolism was highlighted at a very low BMC (first quartile, 0.003 µM). Moreover, results confirmed that lasiocarpine targets bile acids, as well as amino acid and steroid metabolisms. Analysis of the endometabolome, based on coupling concentration-response and PODs, gave encouraging results for ranking toxins according to their hepatotoxic effects. Therefore, this novel approach is a promising tool for next-generation risk assessment, readily applicable to a broad range of compounds and toxic endpoints.

**Abbreviations:** NGRA, next-generation risk assessment; RP, reference point; NOAEL, no-observed-adverse-effect level; BMC, benchmark concentration; PA, pyrrolizidine alkaloid; LC-HRMS, liquid chromatography-high-resolution mass spectrometry; DMSO, dimethylsulfoxide; UHPLC, ultra-high-performance liquid chromatography; HEL, heliotrine; RET, retrorsine; LAS, lasiocarpine; EFSA, European Food Safety Agency; AOP, adverse outcome pathway; MOA, mode of action; POD, point of departure.

<footnote><sup>a</sup> ANSES, French Agency for Food, Environmental and Occupational Health & Safety, Fougères Laboratory, Toxicology of Contaminants Unit, 10 B rue Claude Bourgelat, 35306 Fougères, France</footnote>

<footnote><sup>b</sup> ANSES, French Agency for Food, Environmental and Occupational Health & Safety, Risk Assessment Department, 14 Rue Pierre et Marie Curie, 94701 Maisons-Alfort, France</footnote>

<footnote><sup>c</sup> University of Lyon 1, CNRS, VetAgro Sup, UMR 5558, Laboratoire de Biométrie et Biologie Evolutive, 69622 Villeurbanne, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: estelle.dubreil@anses.fr (E. Dubreil).</footnote>
<figure><img src="image_2.png" /></figure>

# Modeling HepaRG metabolome responses to pyrrolizidine alkaloid exposure for insight into points of departure and modes of action

Estelle Dubreil<sup>a, *</sup>, Keyvin Darney<sup>b</sup>, Marie-Laure Delignette-Muller<sup>c</sup>, Audrey Barranger<sup>a</sup>, Sylvie Huet<sup>a</sup>, Kevin Hogeveen<sup>a</sup>, Thibaut Léger<sup>a</sup>, Valérie Fessard<sup>a</sup>, Ludovic Le Hégarat<sup>a</sup>

## HIGHLIGHTS

* PA effects on HepaRG cells were studied by concentration-response modelling and BMC.

* The endometabolome was found to be more sensitive than the exometabolome.

* LAS had the highest hepatotoxic potency, followed by HEL and RET.

* The bell shape is the most frequent trend in response to LAS.

* Lipids, bile acids and amino acids were especially affected by LAS.

## GRAPHICAL ABSTRACT

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

**Keywords:**

Alkaloid

Food contaminant

Dose-response model

Benchmark concentration

Metabolomics

## ABSTRACT

The new challenges in toxicology demand novel and innovative in vitro approaches for deriving points of departure (PODs) and determining the mode of action (MOA) of chemicals. Therefore, the aim of this original study was to couple in vitro studies with untargeted metabolomics to model the concentration-response of extra- and intracellular metabolome data on human HepaRG cells treated for 48 h with three pyrrolizidine alkaloids (PAs): heliotrine, retrorsine and lasiocarpine. Modeling revealed that the three PAs induced various monotonic and, importantly, biphasic curves of metabolite content. Based on unannotated metabolites, the endometabolome was more sensitive than the exometabolome in terms of metabolomic effects, and benchmark concentrations (BMCs) confirmed that lasiocarpine was the most hepatotoxic PA. Regarding its MOA, impairment of lipid metabolism was highlighted at a very low BMC (first quartile, 0.003 µM). Moreover, results confirmed that lasiocarpine targets bile acids, as well as amino acid and steroid metabolisms. Analysis of the endometabolome, based on coupling concentration-response and PODs, gave encouraging results for ranking toxins according to their hepatotoxic effects. Therefore, this novel approach is a promising tool for next-generation risk assessment, readily applicable to a broad range of compounds and toxic endpoints.

**Abbreviations:** NGRA, next-generation risk assessment; RP, reference point; NOAEL, no-observed-adverse-effect level; BMC, benchmark concentration; PA, pyrrolizidine alkaloid; LC-HRMS, liquid chromatography-high-resolution mass spectrometry; DMSO, dimethylsulfoxide; UHPLC, ultra-high-performance liquid chromatography; HEL, heliotrine; RET, retrorsine; LAS, lasiocarpine; EFSA, European Food Safety Agency; AOP, adverse outcome pathway; MOA, mode of action; POD, point of departure.

* Corresponding author.

E-mail address: estelle.dubreil@anses.fr (E. Dubreil).

<footnote><sup>a</sup> ANSES, French Agency for Food, Environmental and Occupational Health & Safety, Fougères Laboratory, Toxicology of Contaminants Unit, 10 B rue Claude Bourgelat, 35306 Fougères, France</footnote>

<footnote><sup>b</sup> ANSES, French Agency for Food, Environmental and Occupational Health & Safety, Risk Assessment Department, 14 Rue Pierre et Marie Curie, 94701 Maisons-Alfort, France</footnote>

<footnote><sup>c</sup> University of Lyon 1, CNRS, VetAgro Sup, UMR 5558, Laboratoire de Biométrie et Biologie Evolutive, 69622 Villeurbanne, France</footnote>
# 1. Introduction

Pyrrolizidine alkaloids (PAs) are secondary phytochemicals, non-essential metabolites synthesized by plants for defense or reproduction [1]. PAs can cause severe food intoxication in humans and animals upon consumption of various natural foods and drinks (e.g. tea, honey, or edible plants) [2-7]. Their ubiquity in plants along with their potential toxicity make them a major food safety concern [8-10]. More than 660 different chemical structures have been identified in approximately 6000 plant species [10,11]. Teas and herbal infusions as well as other food items like honey can contain many PAs. A Europe wide survey of 1105 samples of animal- and plant-derived food collected in 2014 and 2015 showed that 91 % and 60 % of (herbal) teas and food supplements, respectively, contained at least one individual PA, with a mean concentration of 460 µg kg⁻¹ dry weight for teas [12]. Another study on 121 honey samples, mostly from Calabria (Italy), reported that 38 samples (31 %) contained PAs [13]. According to the European Commission's Rapid Alert System for Food and Feed (RASFF) portal (https://webgate.ec.europa.eu/rasff-window/screen/search), 25 notifications warned about PAs in 2022, mostly in cumin and oregano. Based on current data on exposure to PAs particularly in teas, herbal infusions and honey, the EFSA CONTAM Panel established a Reference Point of 237 µg/kg body weight per day to assess the carcinogenic risks of PAs; this RP is likely to be updated in the near future [9].

PAs share a common structure of a tertiary amine esterified with one or two necic acids, making it a necine base. PAs are classified depending on their type of necine base, which can be a retronecine, heliotridine, otonecine or platynecine type. Apart from platynecine types, all PAs possess a double bond at position 1,2 of the pyrrolizidine ring system within the necine base (1,2-unsatured PAs), which confer their hepatotoxicity [14]. Depending on their degree of esterification, PAs can be classified as monoesters (e.g. heliotrine, HEL), open diesters (e.g. lasiocarpine, LAS) or cyclic diesters (e.g. retrorsine, RET) [15]. The toxicity of PAs arises upon their activation, catalyzed by cytochrome P450 (CYP450) enzymes, resulting in dehydro-PAs and forming reactive pyrrolic esters that interact with nucleophilic sites in proteins and DNA [16]. In vivo studies have shown that several PAs can cause acute and chronic liver damage and have proven potency as hepatocarcinogens in rodents [15,17,18]. Some PAs (monocrotaline, riddelliine and lasiocarpine) have been classified as potential carcinogens to humans (group 2B) by the International Agency for Research on Cancer (IARC) [19].

Although the hepatotoxicity of PAs has been well demonstrated, risk assessment can be refined by identifying biomarkers or by describing the molecular mechanisms underlying PA-induced toxicity [15]. To do so, PA-induced metabolic changes in the liver can be examined with untargeted metabolomics through ultra-high-performance liquid chromatography (UPHLC) coupled with high-resolution mass spectrometry (HRMS), as already used for other hepatotoxic compounds [20-23]. A similar approach has been applied to determine PA hepatotoxicity in vivo. The study concluded that hepatotoxicity of the cyclic diester sen-ecionine is linked to disrupted homeostasis of bile acids in rat serum [24, 25]. Another study, on the PA-containing plant *Gynura segetum*, suggested that 10 metabolites involved in amino acids, lipids and energy metabolism pathways may be related to potential liver injury [26]. Recently, four different PAs, including heliotrine, were shown to impair bile acid homeostasis and secretion in HepaRG cells by downregulating hepatobiliary transporters — enzymes mediating bile acid synthesis and conjugation — and several transcription regulators. The results suggested that retronecine-type PAs are more potent apoptotic compounds than heliotridine- and otonecine-type PAs [27].

Omics data, such as provided by metabolomics, have been used to investigate the adverse outcome pathways (AOPs) and modes of action (MOAs) of hepatotoxicants [28,29]. Nowadays, omics is being explored as a promising approach for next-generation risk assessment (NGRA) with suggested thresholds such as reference points (RPs). Also known as the “point of departure” (POD), the RP is defined as the dose level

corresponding to the onset of the compound's effects [30,31]. The Scientific Committee of the European Food Safety Authority (EFSA) has advised using benchmark doses (BMDs), which are more relevant than the no-observed-adverse-effect level (NOAEL), to derive RPs [32,33]. A BMD is a concentration level corresponding to a predefined difference in response from the control. It is estimated from a fitted concentration-response curve, as the concentration at which the predefined benchmark response (BMR) is reached. However, to date, there are no recommendations on a concentration-response framework for omics data nor case studies demonstrating the application of metabolomics to derive RPs using the BMD approach [31,34]. The RP method based on omics data comes under the concept of new approach methods (NAMs) aiming to facilitate risk assessment without generating new in vivo data (US-EPA, https://www.epa.gov/chemical-research/epa-new-approach-methods-work-plan-reducing-use-vertebrate-animals-chemical). For instance, RPs derived from BMDs (in vivo) or BMCs (in vitro) established from transcriptomics data can predict RPs from conventional “endpoint” studies, thereby illustrating the benefit of transcriptomics for risk assessment [30,35]. Few studies have demonstrated that metabolomics can also predict RPs, although metabolomics-derived RPs have been shown to be more sensitive, quantitative indicators of human liver injury than transcriptomic-derived RPs [36]. More recently, one study showed that high-throughput metabolomics could be used to assess the BMCs of four chemicals on HepaRG cells [37]. It has also been suggested that metabolomic fingerprints can be aggregated to obtain a quantitative assessment of responses from which to derive RPs to substantiate the ecotoxicological risk assessment of chemicals [38,39]. In this approach, the continuous responses of metabolites are modeled over a gradient of concentrations using open-source software such as BMDExpress-2 (https://github.com/auerbachs/BMDExpress-2/releases), PROAST (https://www.rivm.nl/en/proast), EFSA software (https://efsa.openanalytics.eu/app/bmd), or the DRomics package developed by French teams supported by an MSC H2020 grant (https://github.com/aursiber/DRomics) [40].

Qualitative detection of early changes has demonstrated the relevance of using untargeted metabolomics to construct comprehensive MOAs of (hepato)toxic contaminants [41,42]. Here, we propose an innovative step forward by combining metabolomics with concentration-response modeling. These original results can benefit NGRA by identifying early response biomarkers or even molecular pathways at low BMCs, but also by providing insight into the biological mechanisms based on response trends. Moreover, a recent breakthrough in the quantitative assessment of AOPs has made it possible to incorporate complex biological mechanisms. Here, we compared and characterized the hepatotoxic responses of three PAs based on the exometabolome (or extracellular metabolome) and endometabolome (or intracellular metabolome) of HepaRG cells at a first-tier level. Using the DRomics turnkey tool, we modeled the responses of significant unannotated metabolites following exposure to increasing concentrations of three PAs to determine BMCs and describe their trends [40,43]. At a second-tier level, the identified hepatotoxicity biomarkers along with their associated pathways and BMCs were combined to suggest a MOA for LAS. The conclusions of this study help clarify dose-response modeling in metabolomics to advocate its use in risk assessment. This approach is promising for determining the hepatotoxicity of similar compounds belonging to the same family using in vitro systems and for improving the predictivity of potential hazards.

# 2. Materials and methods

## 2.1. Chemicals and reagents

High-purity PAs (> 97 % for LAS, and > 98 % for RET and HEL) were purchased from Phytoplan Diehm & Neuberger GmbH (Heidelberg, Germany). Stock solutions of PAs were prepared in dimethylsulfoxide (DMSO) and stored at -20 °C. The internal standards (ISs) diclofenac-d₄ and flunixine-d₃ were purchased from Witega Laboratorien Berlin-
# 1. Introduction

Pyrrolizidine alkaloids (PAs) are secondary phytochemicals, non-essential metabolites synthesized by plants for defense or reproduction [1]. PAs can cause severe food intoxication in humans and animals upon consumption of various natural foods and drinks (e.g. tea, honey, or edible plants) [2-7]. Their ubiquity in plants along with their potential toxicity make them a major food safety concern [8-10]. More than 660 different chemical structures have been identified in approximately 6000 plant species [10,11]. Teas and herbal infusions as well as other food items like honey can contain many PAs. A Europe wide survey of 1105 samples of animal- and plant-derived food collected in 2014 and 2015 showed that 91 % and 60 % of (herbal) teas and food supplements, respectively, contained at least one individual PA, with a mean concentration of 460 µg kg⁻¹ dry weight for teas [12]. Another study on 121 honey samples, mostly from Calabria (Italy), reported that 38 samples (31 %) contained PAs [13]. According to the European Commission's Rapid Alert System for Food and Feed (RASFF) portal (https://webgate.ec.europa.eu/rasff-window/screen/search), 25 notifications warned about PAs in 2022, mostly in cumin and oregano. Based on current data on exposure to PAs particularly in teas, herbal infusions and honey, the EFSA CONTAM Panel established a Reference Point of 237 µg/kg body weight per day to assess the carcinogenic risks of PAs; this RP is likely to be updated in the near future [9].

PAs share a common structure of a tertiary amine esterified with one or two necic acids, making it a necine base. PAs are classified depending on their type of necine base, which can be a retronecine, heliotridine, otonecine or platynecine type. Apart from platynecine types, all PAs possess a double bond at position 1,2 of the pyrrolizidine ring system within the necine base (1,2-unsatured PAs), which confer their hepatotoxicity [14]. Depending on their degree of esterification, PAs can be classified as monoesters (e.g. heliotrine, HEL), open diesters (e.g. lasiocarpine, LAS) or cyclic diesters (e.g. retrorsine, RET) [15]. The toxicity of PAs arises upon their activation, catalyzed by cytochrome P450 (CYP450) enzymes, resulting in dehydro-PAs and forming reactive pyrrolic esters that interact with nucleophilic sites in proteins and DNA [16]. In vivo studies have shown that several PAs can cause acute and chronic liver damage and have proven potency as hepatocarcinogens in rodents [15,17,18]. Some PAs (monocrotaline, riddelliine and lasiocarpine) have been classified as potential carcinogens to humans (group 2B) by the International Agency for Research on Cancer (IARC) [19].

Although the hepatotoxicity of PAs has been well demonstrated, risk assessment can be refined by identifying biomarkers or by describing the molecular mechanisms underlying PA-induced toxicity [15]. To do so, PA-induced metabolic changes in the liver can be examined with untargeted metabolomics through ultra-high-performance liquid chromatography (UPHLC) coupled with high-resolution mass spectrometry (HRMS), as already used for other hepatotoxic compounds [20-23]. A similar approach has been applied to determine PA hepatotoxicity in vivo. The study concluded that hepatotoxicity of the cyclic diester sen-ecionine is linked to disrupted homeostasis of bile acids in rat serum [24, 25]. Another study, on the PA-containing plant *Gynura segetum*, suggested that 10 metabolites involved in amino acids, lipids and energy metabolism pathways may be related to potential liver injury [26]. Recently, four different PAs, including heliotrine, were shown to impair bile acid homeostasis and secretion in HepaRG cells by downregulating hepatobiliary transporters — enzymes mediating bile acid synthesis and conjugation — and several transcription regulators. The results suggested that retronecine-type PAs are more potent apoptotic compounds than heliotridine- and otonecine-type PAs [27].

Omics data, such as provided by metabolomics, have been used to investigate the adverse outcome pathways (AOPs) and modes of action (MOAs) of hepatotoxicants [28,29]. Nowadays, omics is being explored as a promising approach for next-generation risk assessment (NGRA) with suggested thresholds such as reference points (RPs). Also known as the “point of departure” (POD), the RP is defined as the dose level

corresponding to the onset of the compound's effects [30,31]. The Scientific Committee of the European Food Safety Authority (EFSA) has advised using benchmark doses (BMDs), which are more relevant than the no-observed-adverse-effect level (NOAEL), to derive RPs [32,33]. A BMD is a concentration level corresponding to a predefined difference in response from the control. It is estimated from a fitted concentration-response curve, as the concentration at which the predefined benchmark response (BMR) is reached. However, to date, there are no recommendations on a concentration-response framework for omics data nor case studies demonstrating the application of metabolomics to derive RPs using the BMD approach [31,34]. The RP method based on omics data comes under the concept of new approach methods (NAMs) aiming to facilitate risk assessment without generating new in vivo data (US-EPA, https://www.epa.gov/chemical-research/epa-new-approach-methods-work-plan-reducing-use-vertebrate-animals-chemical). For instance, RPs derived from BMDs (in vivo) or BMCs (in vitro) established from transcriptomics data can predict RPs from conventional “endpoint” studies, thereby illustrating the benefit of transcriptomics for risk assessment [30,35]. Few studies have demonstrated that metabolomics can also predict RPs, although metabolomics-derived RPs have been shown to be more sensitive, quantitative indicators of human liver injury than transcriptomic-derived RPs [36]. More recently, one study showed that high-throughput metabolomics could be used to assess the BMCs of four chemicals on HepaRG cells [37]. It has also been suggested that metabolomic fingerprints can be aggregated to obtain a quantitative assessment of responses from which to derive RPs to substantiate the ecotoxicological risk assessment of chemicals [38,39]. In this approach, the continuous responses of metabolites are modeled over a gradient of concentrations using open-source software such as BMDExpress-2 (https://github.com/auerbachs/BMDExpress-2/releases), PROAST (https://www.rivm.nl/en/proast), EFSA software (https://efsa.openanalytics.eu/app/bmd), or the DRomics package developed by French teams supported by an MSC H2020 grant (https://github.com/aursiber/DRomics) [40].

Qualitative detection of early changes has demonstrated the relevance of using untargeted metabolomics to construct comprehensive MOAs of (hepato)toxic contaminants [41,42]. Here, we propose an innovative step forward by combining metabolomics with concentration-response modeling. These original results can benefit NGRA by identifying early response biomarkers or even molecular pathways at low BMCs, but also by providing insight into the biological mechanisms based on response trends. Moreover, a recent breakthrough in the quantitative assessment of AOPs has made it possible to incorporate complex biological mechanisms. Here, we compared and characterized the hepatotoxic responses of three PAs based on the exometabolome (or extracellular metabolome) and endometabolome (or intracellular metabolome) of HepaRG cells at a first-tier level. Using the DRomics turnkey tool, we modeled the responses of significant unannotated metabolites following exposure to increasing concentrations of three PAs to determine BMCs and describe their trends [40,43]. At a second-tier level, the identified hepatotoxicity biomarkers along with their associated pathways and BMCs were combined to suggest a MOA for LAS. The conclusions of this study help clarify dose-response modeling in metabolomics to advocate its use in risk assessment. This approach is promising for determining the hepatotoxicity of similar compounds belonging to the same family using in vitro systems and for improving the predictivity of potential hazards.

# 2. Materials and methods

## 2.1. Chemicals and reagents

High-purity PAs (> 97 % for LAS, and > 98 % for RET and HEL) were purchased from Phytoplan Diehm & Neuberger GmbH (Heidelberg, Germany). Stock solutions of PAs were prepared in dimethylsulfoxide (DMSO) and stored at -20 °C. The internal standards (ISs) diclofenac-d₄ and flunixine-d₃ were purchased from Witega Laboratorien Berlin-
Adlershof GmbH (Berlin, Germany). The stock solutions of ISs were prepared in methanol and stored at – 20 °C. DMSO was purchased from Sigma (St. Quentin-Fallavier, France). Liquid chromatography–mass spectrometry (LC-MS) grade acetonitrile, methanol, water, and formic acid were purchased from Thermo Fisher Scientific (Waltham, MA, USA).

## 2.2. HepaRG cell culture

Differentiated HepaRG cells, basal hepatic medium and supplements were acquired from Biopredic International (Saint-Grégoire, France). Cells were cultured in Williams E medium (Eurobio, Les Ulis, France) supplemented with 10 % fetal calf serum (FCS) (Perbio, Brebières, France), 100 units/mL penicillin (Invitrogen Corporation, Illkirch, France), 100 µg/mL streptomycin (Invitrogen Corporation), 5 µg/mL insulin (Sigma Aldrich, Lyon, France), 2 mM L-glutamine (Thermo Fisher Scientific, Waltham, MA, USA), and 25 µg/mL hydrocortisone succinate (Pharmacia & Upjohn, Guyancourt, France). Cells were seeded in a 75 cm² flask at a density of $2 \times 10^6$ then incubated at 37 °C with 5 % CO₂. For differentiation purposes, after 14 days of culture in a 75 cm² flask, the cells were seeded at 26,000 cells/cm² in either 96-well plates (cytotoxicity test) or 24-well plates (metabolomics), then incubated for two weeks before 1.7 % DMSO was added to the culture medium, and then incubated for two more weeks. In all cases, the medium was changed three times a week.

## 2.3. Cytotoxicity

### 2.3.1. MTT assay

To identify appropriate concentrations for metabolomics experiments, the subtoxic concentrations of PAs were determined using the MTT assay. Working solutions were freshly prepared by serial dilution in FCS-free medium, with a final concentration of 0.5 % DMSO. HepaRG cells were treated with a concentration range of each PA for 48 h, i.e. 0.39 to 200 µM for HEL and RET, 0.1 to 50 µM for LAS. At the end of the treatment, MTT (500 µg/mL in medium) was added to each well and incubated for 1 h at 37 °C. After solubilization with DMSO, and shaking for 10 min, absorbance was read at 570 nm using a FLUOstar OPTIMA

microplate reader (BMG Labtek, Champigny-sur-Marne, France). Viability was calculated as the percentage of mean absorbance relative to the solvent control condition (0.5 % DMSO). Technical triplicates per concentration and three independent experiments were performed.

### 2.3.2. Statistical analysis

For statistics, a one-way analysis of variance (ANOVA) and Dunnet's multiple comparison test were used; P values < 0.05 were considered statistically significant (Prism 9.0, GraphPad, La Jolla, CA, USA).

## 2.4. Experimental design for metabolomics

The experimental design is depicted in Fig. 1. For the analysis of the endo- and exometabolome, HepaRG cells (24-well plates in triplicate) were exposed for 48 h to 10 concentrations (N = 3 each) (0.006–0.016–0.041–0.102–0.256–0.64–1.6–4–10–25 µM) for LAS and RET, and 10 concentrations (N = 3 each) (0.026–0.065–0.16–0.41–0.102–2.56–6.4–16–40–100 µM) for HEL. The treatment medium contained 0.5 % DMSO without FCS. Twelve wells on each plate were used as a solvent control (0.5 % DMSO). For the exometabolome study, three control media without DMSO were also added.

## 2.5. Metabolome analysis

### 2.5.1. Exometabolome (EXO)

After treatment, the plates were immediately placed on ice. For each well, the medium was transferred into 2-mL centrifuge tubes and stored at – 80 °C until further extraction. After thawing, 300 µL of media was transferred into 5-mL polypropylene tubes. Then, 5 µL of IS solution (see 2.1) was added. Next, 1000 µL of a cold solution of methanol and acetonitrile (50:50, v:v) was added to each tube and vortexed for 20 s to extract the metabolites. After centrifugation (14,000 g, 4 °C, 5 min), the supernatants were transferred into new 5-mL polypropylene tubes and evaporated under a nitrogen flow at 30 °C. The extracts were dissolved in 80 µL of water and acetonitrile (95:5, v:v) and transferred into a vial for UHPLC-HRMS analysis. The quality control (QC) pool for EXO analysis was generated by collecting 5 µL aliquots of each sample.

<figure><img src="image_5.png" /><figcaption>Fig. 1. Experimental design.</figcaption></figure>
Adlershof GmbH (Berlin, Germany). The stock solutions of ISs were prepared in methanol and stored at – 20 °C. DMSO was purchased from Sigma (St. Quentin-Fallavier, France). Liquid chromatography–mass spectrometry (LC-MS) grade acetonitrile, methanol, water, and formic acid were purchased from Thermo Fisher Scientific (Waltham, MA, USA).

## 2.2. HepaRG cell culture

Differentiated HepaRG cells, basal hepatic medium and supplements were acquired from Biopredic International (Saint-Grégoire, France). Cells were cultured in Williams E medium (Eurobio, Les Ulis, France) supplemented with 10 % fetal calf serum (FCS) (Perbio, Brebières, France), 100 units/mL penicillin (Invitrogen Corporation, Illkirch, France), 100 µg/mL streptomycin (Invitrogen Corporation), 5 µg/mL insulin (Sigma Aldrich, Lyon, France), 2 mM L-glutamine (Thermo Fisher Scientific, Waltham, MA, USA), and 25 µg/mL hydrocortisone succinate (Pharmacia & Upjohn, Guyancourt, France). Cells were seeded in a 75 cm² flask at a density of $2 \times 10^6$ then incubated at 37 °C with 5 % CO₂. For differentiation purposes, after 14 days of culture in a 75 cm² flask, the cells were seeded at 26,000 cells/cm² in either 96-well plates (cytotoxicity test) or 24-well plates (metabolomics), then incubated for two weeks before 1.7 % DMSO was added to the culture medium, and then incubated for two more weeks. In all cases, the medium was changed three times a week.

## 2.3. Cytotoxicity

### 2.3.1. MTT assay

To identify appropriate concentrations for metabolomics experiments, the subtoxic concentrations of PAs were determined using the MTT assay. Working solutions were freshly prepared by serial dilution in FCS-free medium, with a final concentration of 0.5 % DMSO. HepaRG cells were treated with a concentration range of each PA for 48 h, i.e. 0.39 to 200 µM for HEL and RET, 0.1 to 50 µM for LAS. At the end of the treatment, MTT (500 µg/mL in medium) was added to each well and incubated for 1 h at 37 °C. After solubilization with DMSO, and shaking for 10 min, absorbance was read at 570 nm using a FLUOstar OPTIMA

microplate reader (BMG Labtek, Champigny-sur-Marne, France). Viability was calculated as the percentage of mean absorbance relative to the solvent control condition (0.5 % DMSO). Technical triplicates per concentration and three independent experiments were performed.

### 2.3.2. Statistical analysis

For statistics, a one-way analysis of variance (ANOVA) and Dunnet's multiple comparison test were used; P values < 0.05 were considered statistically significant (Prism 9.0, GraphPad, La Jolla, CA, USA).

## 2.4. Experimental design for metabolomics

The experimental design is depicted in Fig. 1. For the analysis of the endo- and exometabolome, HepaRG cells (24-well plates in triplicate) were exposed for 48 h to 10 concentrations (N = 3 each) (0.006–0.016–0.041–0.102–0.256–0.64–1.6–4–10–25 µM) for LAS and RET, and 10 concentrations (N = 3 each) (0.026–0.065–0.16–0.41–0.102–2.56–6.4–16–40–100 µM) for HEL. The treatment medium contained 0.5 % DMSO without FCS. Twelve wells on each plate were used as a solvent control (0.5 % DMSO). For the exometabolome study, three control media without DMSO were also added.

## 2.5. Metabolome analysis

### 2.5.1. Exometabolome (EXO)

After treatment, the plates were immediately placed on ice. For each well, the medium was transferred into 2-mL centrifuge tubes and stored at – 80 °C until further extraction. After thawing, 300 µL of media was transferred into 5-mL polypropylene tubes. Then, 5 µL of IS solution (see 2.1) was added. Next, 1000 µL of a cold solution of methanol and acetonitrile (50:50, v:v) was added to each tube and vortexed for 20 s to extract the metabolites. After centrifugation (14,000 g, 4 °C, 5 min), the supernatants were transferred into new 5-mL polypropylene tubes and evaporated under a nitrogen flow at 30 °C. The extracts were dissolved in 80 µL of water and acetonitrile (95:5, v:v) and transferred into a vial for UHPLC-HRMS analysis. The quality control (QC) pool for EXO analysis was generated by collecting 5 µL aliquots of each sample.

<figure><img src="image_4.png" /><figcaption>Fig. 1. Experimental design.</figcaption></figure>
## 2.5.2. Endometabolome (ENDO)

Once the cell culture medium was collected (for EXO analysis, see previous section), cells were washed twice with 300 µL of a cold 0.9 % NaCl buffer. Then 300 µL of cold methanol, 100 µL of cold water and 5 µL of IS solution was added to each well. The plates were ultra-sonicated (40 kHz) for 30 min in a water bath to detach the cells. Any adhering cells left over were scraped off. The extracted metabolites were transferred into 2-mL centrifuge tubes and centrifuged at 14,000 g at 4 °C for 5 min. The supernatants were transferred into 5-mL poly-propylene tubes prior to evaporation under a nitrogen flow at 30 °C. The extracts were dissolved in 40 µL of water and acetonitrile (95:5, v:v) and transferred into a vial for UHPLC-HRMS analysis. The QC pool for ENDO analysis was generated by collecting 5 µL aliquots of each sample.

## 2.5.3. UHPLC-HRMS analysis

EXO and ENDO triplicate samples, the 12 solvent controls, the three medium controls and the QC sample were all analyzed using a Q-Exactive Plus mass spectrometer coupled with a Vanquish Flex binary UHPLC system (both from Thermo Fisher Scientific, Waltham, MA, USA). The analysis was conducted according to Léger et al. [44]. Briefly, UHPLC was performed on an Acquity UPLC HSS T3 (1.8 µm, 1 mm × 150 mm, 100 Å) at 80 µL min⁻¹ flow. Mobile phases were composed of solvent A (water with 0.1 % formic acid) and solvent B (acetonitrile with 0.1 % formic acid). Mobile phase B was ramped linearly from 5 % to 10 % over 2.5 min, to 45 % over 7.5 min, to 65 % over 10 min, and to 100 % over 3 min. Finally, an isocratic flow at 100 % solvent B was applied for 3 min, and a decrease to 5 % in 1 min to return to initial conditions; this was followed by an equilibration step for 2 min, giving a total run time of 28 min. Metabolites were analyzed in the Orbitrap in full ion scan mode at a resolution of 70,000 (at *m/z* 200) with a mass range of *m/z* 70–1050. The samples were acquired in positive and negative electro-spray ionization modes (ESI+, ESI-). Fragments were obtained by higher energy collisional dissociation (HCD) activation with a collisional energy of 27 %, and an isolation width of 1.4 Da. MS/MS data were acquired in the Orbitrap cell in a data-dependent mode, in which the 20 most intense precursor ions were fragmented, with a dynamic exclusion of 30 s and a resolution of 17,500. Automatic gain control (AGC) targets were set to 3 × 10⁶ and 2 × 10⁵ for the MS and MS/MS steps, respectively. Maximum ion accumulation times were set to 50 and 45 ms for MS and MS/MS acquisitions, respectively. QC samples (for EXO and ENDO) were included every nine injections to check acquisition reproducibility. The ISs were visually inspected in QC samples regarding signal intensity, peak shape, and retention time.

Compound Discoverer version 3.2 (Thermo Fisher Scientific, Waltham, MA, USA) was used for data preprocessing. Results were filtered to minimize interference and background levels, specifically filtering for retention time (RT) greater than 2 min, a mass shift less than 5 ppm, a minimum peak intensity of 1000 in positive mode and 500 in negative mode, with pooled quality control samples (QC correction) to compensate for time-dependent batch effects. The relative standard deviation of detected peak areas in the QC samples before and after correction did not exceed the 30 % threshold. The minimum and maximum number of elements were as follows: 0 < N < 10; C < 90; 0 < 18; H > 190; P < 3; S < 5; Br < 3; Cl < 4. The maximum width of peaks selected was 0.4 min. The minimum number of scans per peak was set to 5. The fill gaps algorithm was applied for missing peaks (parameters of 5 ppm mass tolerance and of 1.5 S/N threshold). The measured areas were normalized using a constant median normalization method. The automatic subtraction of blanks was included for EXO datasets.

All area under the curve (AUC) values were transferred to Excel sheets for the 12 datasets: i.e. EXO analyzed in ESI+ and ESI- modes, ENDO analyzed in ESI+ and ESI- modes, and for each of the three PAs. The AUC values were log10-transformed.

## 2.6. Statistical analyses of EXO and ENDO

### 2.6.1. Data verification

The 12 EXO and ENDO datasets (in positive and negative ion modes) for the three PAs were analyzed using the DRomics tool (Larras et al., 2018) designed for concentration-response datasets from omics experiments. Data were analyzed using the DRomics shiny application (Shiny version, https://lbbe.univ-lyon1.fr/dromics.html) and the DRomics R package version 2.5-0 [45]. The data, i.e. AUC values log₁₀-transformed by data preprocessing, were imported into the DRomics shiny application (step 1). A principal component analysis (PCA) of each dataset was used to check the data after normalization and transformation, and to remove any outliers.

### 2.6.2. Selection of significantly responding metabolites

In step 2, the responding metabolites (named items in DRomics) (i.e. those displaying a significant change in signal with increasing PA concentrations) were selected using a quadratic trend test with a false discovery rate of 0.05 controlled by a Benjamini-Hochberg correction on the p-values, as suggested in Larras et al. (2018) [40].

### 2.6.3. Fitting the concentration-response models

In step 3, concentration-response was modeled for each of the responding metabolites. Some metabolites were also eliminated at this step when no model could be fitted (automatically carried out in DRomics). The monotonic and non-monotonic response trends were identified at this step based on the best-fit curves among the five models implemented in DRomics, using an information criterion (second-order Akaike criterion) to prevent overfitting. The best-fit model, its parameter values, the standard residual error (SDres) were reported as well as the trend of the curve: increasing and decreasing for monotonic curves, and U-shaped and bell-shaped for non-monotonic (i.e. biphasic) curves. A U shape corresponds to a decrease of the response followed by an increase, and a bell shape corresponds to an increase followed by a decrease of the response.

### 2.6.4. Calculation of benchmark concentrations (BMCs)

In step 4, the best-fit model for each metabolite was used to assess BMCs according to EFSA guidelines [32,33]. The BMC does not refer to a specific level of effect, but indicates the concentration leading to a level of change compared with the control response that takes into account the data variability around the modeled curve. It is calculated as the concentration corresponding to a benchmark response (BMR<sub>zSD</sub>) defined as follows: BMR<sub>zSD</sub>=y₀ + /- z * SD, where y₀ is the mean control response, SD is the residual standard deviation of the considered concentration-response fitted model, and z is a factor set to 1 in this work (default value proposed in DRomics and recommended by EFSA [32]). The individual BMC values were plotted with their confidence intervals using a bootstrap procedure. The empirical cumulative distribution functions (ECDF) of the BMCs were then plotted (i.e. ECDF of unannotated metabolites).

## 2.7. Annotation of endometabolites affected by lasiocarpine treatment (LAS ENDO)

The DRomics tool selected endometabolites affected by LAS treatment, which were then annotated. The identification was based on the more common adducts ([M+H]⁺, [M+Na]⁺, [M+K]⁺, [M+NH₄]⁺, [M]⁺) in ESI+ mode and ([M-H]⁻, [M+HCOOH-H]⁻, [M+CH₃COOH-H]⁻, [M-H-H₂O]⁻) in ESI- mode and took into consideration the *m/z*-value, the potential formula generated by the Compound Discoverer algorithm and checked in Chemspider, the isotope pattern, the measured retention time and the fragmentation spectra when acquired with DDA conditions. Endometabolites were first selected when identified against an in-house library (HepaRG DB) and an open-access database MTOX700 + for predicting toxicity [46] based on mass accuracy and retention times
## 2.5.2. Endometabolome (ENDO)

Once the cell culture medium was collected (for EXO analysis, see previous section), cells were washed twice with 300 µL of a cold 0.9 % NaCl buffer. Then 300 µL of cold methanol, 100 µL of cold water and 5 µL of IS solution was added to each well. The plates were ultra-sonicated (40 kHz) for 30 min in a water bath to detach the cells. Any adhering cells left over were scraped off. The extracted metabolites were transferred into 2-mL centrifuge tubes and centrifuged at 14,000 g at 4 °C for 5 min. The supernatants were transferred into 5-mL poly-propylene tubes prior to evaporation under a nitrogen flow at 30 °C. The extracts were dissolved in 40 µL of water and acetonitrile (95:5, v:v) and transferred into a vial for UHPLC-HRMS analysis. The QC pool for ENDO analysis was generated by collecting 5 µL aliquots of each sample.

## 2.5.3. UHPLC-HRMS analysis

EXO and ENDO triplicate samples, the 12 solvent controls, the three medium controls and the QC sample were all analyzed using a Q-Exactive Plus mass spectrometer coupled with a Vanquish Flex binary UHPLC system (both from Thermo Fisher Scientific, Waltham, MA, USA). The analysis was conducted according to Léger et al. [44]. Briefly, UHPLC was performed on an Acquity UPLC HSS T3 (1.8 µm, 1 mm × 150 mm, 100 Å) at 80 µL min⁻¹ flow. Mobile phases were composed of solvent A (water with 0.1 % formic acid) and solvent B (acetonitrile with 0.1 % formic acid). Mobile phase B was ramped linearly from 5 % to 10 % over 2.5 min, to 45 % over 7.5 min, to 65 % over 10 min, and to 100 % over 3 min. Finally, an isocratic flow at 100 % solvent B was applied for 3 min, and a decrease to 5 % in 1 min to return to initial conditions; this was followed by an equilibration step for 2 min, giving a total run time of 28 min. Metabolites were analyzed in the Orbitrap in full ion scan mode at a resolution of 70,000 (at *m/z* 200) with a mass range of *m/z* 70–1050. The samples were acquired in positive and negative electro-spray ionization modes (ESI+, ESI-). Fragments were obtained by higher energy collisional dissociation (HCD) activation with a collisional energy of 27 %, and an isolation width of 1.4 Da. MS/MS data were acquired in the Orbitrap cell in a data-dependent mode, in which the 20 most intense precursor ions were fragmented, with a dynamic exclusion of 30 s and a resolution of 17,500. Automatic gain control (AGC) targets were set to 3 × 10⁶ and 2 × 10⁵ for the MS and MS/MS steps, respectively. Maximum ion accumulation times were set to 50 and 45 ms for MS and MS/MS acquisitions, respectively. QC samples (for EXO and ENDO) were included every nine injections to check acquisition reproducibility. The ISs were visually inspected in QC samples regarding signal intensity, peak shape, and retention time.

Compound Discoverer version 3.2 (Thermo Fisher Scientific, Waltham, MA, USA) was used for data preprocessing. Results were filtered to minimize interference and background levels, specifically filtering for retention time (RT) greater than 2 min, a mass shift less than 5 ppm, a minimum peak intensity of 1000 in positive mode and 500 in negative mode, with pooled quality control samples (QC correction) to compensate for time-dependent batch effects. The relative standard deviation of detected peak areas in the QC samples before and after correction did not exceed the 30 % threshold. The minimum and maximum number of elements were as follows: 0 < N < 10; C < 90; 0 < 18; H > 190; P < 3; S < 5; Br < 3; Cl < 4. The maximum width of peaks selected was 0.4 min. The minimum number of scans per peak was set to 5. The fill gaps algorithm was applied for missing peaks (parameters of 5 ppm mass tolerance and of 1.5 S/N threshold). The measured areas were normalized using a constant median normalization method. The automatic subtraction of blanks was included for EXO datasets.

All area under the curve (AUC) values were transferred to Excel sheets for the 12 datasets: i.e. EXO analyzed in ESI+ and ESI- modes, ENDO analyzed in ESI+ and ESI- modes, and for each of the three PAs. The AUC values were log10-transformed.

## 2.6. Statistical analyses of EXO and ENDO

### 2.6.1. Data verification

The 12 EXO and ENDO datasets (in positive and negative ion modes) for the three PAs were analyzed using the DRomics tool (Larras et al., 2018) designed for concentration-response datasets from omics experiments. Data were analyzed using the DRomics shiny application (Shiny version, https://lbbe.univ-lyon1.fr/dromics.html) and the DRomics R package version 2.5-0 [45]. The data, i.e. AUC values log₁₀-transformed by data preprocessing, were imported into the DRomics shiny application (step 1). A principal component analysis (PCA) of each dataset was used to check the data after normalization and transformation, and to remove any outliers.

### 2.6.2. Selection of significantly responding metabolites

In step 2, the responding metabolites (named items in DRomics) (i.e. those displaying a significant change in signal with increasing PA concentrations) were selected using a quadratic trend test with a false discovery rate of 0.05 controlled by a Benjamini-Hochberg correction on the p-values, as suggested in Larras et al. (2018) [40].

### 2.6.3. Fitting the concentration-response models

In step 3, concentration-response was modeled for each of the responding metabolites. Some metabolites were also eliminated at this step when no model could be fitted (automatically carried out in DRomics). The monotonic and non-monotonic response trends were identified at this step based on the best-fit curves among the five models implemented in DRomics, using an information criterion (second-order Akaike criterion) to prevent overfitting. The best-fit model, its parameter values, the standard residual error (SDres) were reported as well as the trend of the curve: increasing and decreasing for monotonic curves, and U-shaped and bell-shaped for non-monotonic (i.e. biphasic) curves. A U shape corresponds to a decrease of the response followed by an increase, and a bell shape corresponds to an increase followed by a decrease of the response.

### 2.6.4. Calculation of benchmark concentrations (BMCs)

In step 4, the best-fit model for each metabolite was used to assess BMCs according to EFSA guidelines [32,33]. The BMC does not refer to a specific level of effect, but indicates the concentration leading to a level of change compared with the control response that takes into account the data variability around the modeled curve. It is calculated as the concentration corresponding to a benchmark response (BMR<sub>zSD</sub>) defined as follows: BMR<sub>zSD</sub>=y₀ + /- z * SD, where y₀ is the mean control response, SD is the residual standard deviation of the considered concentration-response fitted model, and z is a factor set to 1 in this work (default value proposed in DRomics and recommended by EFSA [32]). The individual BMC values were plotted with their confidence intervals using a bootstrap procedure. The empirical cumulative distribution functions (ECDF) of the BMCs were then plotted (i.e. ECDF of unannotated metabolites).

## 2.7. Annotation of endometabolites affected by lasiocarpine treatment (LAS ENDO)

The DRomics tool selected endometabolites affected by LAS treatment, which were then annotated. The identification was based on the more common adducts ([M+H]⁺, [M+Na]⁺, [M+K]⁺, [M+NH₄]⁺, [M]⁺) in ESI+ mode and ([M-H]⁻, [M+HCOOH-H]⁻, [M+CH₃COOH-H]⁻, [M-H-H₂O]⁻) in ESI- mode and took into consideration the *m/z*-value, the potential formula generated by the Compound Discoverer algorithm and checked in Chemspider, the isotope pattern, the measured retention time and the fragmentation spectra when acquired with DDA conditions. Endometabolites were first selected when identified against an in-house library (HepaRG DB) and an open-access database MTOX700 + for predicting toxicity [46] based on mass accuracy and retention times
when available, then checked in human open-access libraries such as the Human Metabolome Database (HMDB) (www.hmdb.ca) or KEGG Compound Database (https://www.kegg.jp/kegg/compound/). The MS/MS fragmentation spectra using HMDB and mzCloud (www.mzcloud.org) were studied manually to improve annotation confidence. Only metabolites that could be annotated with a confidence level 3 (L3) or higher, according to the annotation confidence system of Schymanski et al., were considered consistent for the study of affected pathways and assigned BMCs [47,48]. L3 refers to tentative candidates (i.e. insufficient information) for one exact structure, e.g. annotation to lipid class level), and level 2 (L2) refers to a probable structure based on an unambiguous library spectrum match or diagnostic MS/MS fragments, without any

available standard or literature for confirmation.

Metabolites identified within pathways in both positive and negative modes were grouped and processed in the DRomics interpreter (Shiny version, https://lbbe-shiny.univ-lyon1.fr/DRomics/inst/DRomicsInterpreter-shiny/). Two main HMDB and KEGG annotation levels were considered: pathway class (e.g. lipid metabolism) and pathway subclass (e.g. fatty acid metabolism, energy metabolism). Only pathways involving at least three responding metabolites were studied for the BMC assessment and characterization of both the trend plot and sensitivity plot by pathway.

<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Number of metabolites (endometabolome (ENDO), exometabolome (EXO), in ES1 + and ESI- modes) in percentages associated with each response trend for HepaRG cells treated for 48 h with pyrrolizidine alkaloids (PAs). Trend: Increase Decrease Bell U.</figcaption></figure>
when available, then checked in human open-access libraries such as the Human Metabolome Database (HMDB) (www.hmdb.ca) or KEGG Compound Database (https://www.kegg.jp/kegg/compound/). The MS/MS fragmentation spectra using HMDB and mzCloud (www.mzcloud.org) were studied manually to improve annotation confidence. Only metabolites that could be annotated with a confidence level 3 (L3) or higher, according to the annotation confidence system of Schymanski et al., were considered consistent for the study of affected pathways and assigned BMCs [47,48]. L3 refers to tentative candidates (i.e. insufficient information) for one exact structure, e.g. annotation to lipid class level), and level 2 (L2) refers to a probable structure based on an unambiguous library spectrum match or diagnostic MS/MS fragments, without any

available standard or literature for confirmation.

Metabolites identified within pathways in both positive and negative modes were grouped and processed in the DRomics interpreter (Shiny version, https://lbbe-shiny.univ-lyon1.fr/DRomics/inst/DRomicsInterpreter-shiny/). Two main HMDB and KEGG annotation levels were considered: pathway class (e.g. lipid metabolism) and pathway subclass (e.g. fatty acid metabolism, energy metabolism). Only pathways involving at least three responding metabolites were studied for the BMC assessment and characterization of both the trend plot and sensitivity plot by pathway.

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Number of metabolites (endometabolome (ENDO), exometabolome (EXO), in ES1 + and ESI- modes) in percentages associated with each response trend for HepaRG cells treated for 48 h with pyrrolizidine alkaloids (PAs). Trend: Increase Decrease Bell U.</figcaption></figure>
# 3. Results

## 3.1. Cytotoxic effects of PAs in HepaRG cells

To determine the appropriate PA concentrations for metabolomic experiments, the viability of HepaRG cells after a 48 h treatment was determined using the MTT assay (Supplementary Fig. 1).

LAS was the most cytotoxic compound on HepaRG cells, inducing a decrease in viability of 46 ± 12 % at the highest concentration tested, 50 µM. A 41 ± 3 % decrease in viability was also observed at 200 µM RET, whereas HEL had only a slight effect on the viability of HepaRG cells.

Based on these results, we selected the highest subtoxic concentrations, not inducing more than 30 % toxicity, for the metabolomic experiments: 25 µM for LAS and RET, and 100 µM for HEL.

## 3.2. ENDO and EXO concentration-response characterization

### 3.2.1. Responses of unannotated metabolites

Monotonic (i.e. increase or decrease) and non-monotonic (i.e. U-shaped or bell-shaped pattern) response trends for increasing concentration gradients were examined for each PA, for ENDO and EXO conditions and in ESI+ and ESI- modes (Fig. 2). The responsive metabolites were hereafter examined without annotation (named unannotated metabolites). First, the total number of responsive ENDO metabolites was higher than the EXO metabolites for LAS (2398 and 1127 metabolites, respectively) and HEL (1537 and 252 metabolites, respectively), and almost equivalent for RET (around 450 metabolites). The metabolites following a non-monotonic curve response (yellow U-shaped curve, and green bell-shaped curve, Fig. 2) ranged from 25 % to 75 % of responsive metabolites. For LAS, a substantial number of the responsive metabolites showed a bell-shaped trend (51 % for ENDO and 39 % for EXO) and approximately a quarter showed an increasing trend. For RET, 34 % showed a linear decreasing trend for ENDO, and a U-shaped trend was predominant for EXO. For HEL, 47 % of the ENDO metabolites showed U-shaped curves; fewer metabolites (252) were obtained for EXO, mostly modeled with a linear decreasing curve.

### 3.2.2. BMCs of the endometabolome (ENDO) from unannotated metabolites

For ENDO, for each PA, an ECDF plot with a color gradient coding for the response was constructed from the BMCs of all responsive unannotated metabolites (Fig. 3A). The BMCs of metabolites identified in HepaRG cells following LAS exposure ranged from almost 0 (0.00006) to 20.5 µM for ESI+ and ESI- modes combined, with the first quartile of metabolites having a BMC (BMC Q1) below 0.01 µM for both modes. For RET, a comparable range of BMCs from 0.0007 to 18.5 µM was found for both modes and BMC Q1 values below 0.5 µM and 0.6 µM for ESI+ and ESI- modes, respectively. For HEL, the BMCs ranged from approximately 0.0003 to 80 µM for both modes, and BMC Q1 values below 0.06 µM and 0.09 µM for ESI+ and ESI- modes, respectively.

### 3.2.3. BMCs of the exometabolome (EXO) from unannotated metabolites

For EXO, a cumulative distribution plot for each PA was constructed from the BMCs of all unannotated metabolites (Fig. 3B). The BMCs of metabolites highlighted after LAS exposure ranged from almost 0 (0.00006) to approximately 20 µM for the ESI+ and ESI- modes combined, with BMC Q1 values below approximately 1.2 µM for both modes. The BMC distribution for RET exposure showed a comparable range, from almost 0 (0.0001) to 19.5 µM for ESI+ and ESI- modes with BMC Q1 values below 0.14 µM and 0.007 µM for the ESI+ and ESI- modes, respectively. Higher BMCs were determined for HEL, with a distribution from almost 0 (0.0003) to approximately 73 µM for ESI+ and ESI- modes, respectively, with BMC Q1 values below 0.2 µM and 3.1 µM for the ESI+ and ESI- modes, respectively.

## 3.3. Metabolic pathway responses of HepaRG cells to lasiocarpine (LAS)

### 3.3.1. Metabolites and associated pathways as markers of hepatotoxicity

Compared with RET and HEL, LAS had a higher number of affected unannotated metabolites for ENDO as well as a lower BMC Q1, indicating a higher sensitivity of HepaRG cells in response to LAS exposure. For these reasons, the data of unannotated endometabolites following LAS treatment (LAS ENDO) were selected for tentative annotation and BMC pathway assignment.

Significant modeled and annotated metabolites associated with biological pathways disrupted in response to 48 h treatment with a LAS concentration range are reported in Table 1. Additional information on annotated metabolites (m/z, ESI mode, adduct, delta mass, MW, proposed formula, database attribution, main fragments, mzCloud scoring) are provided in Supplementary Table 1. Some fitted concentration-response curves are presented in Fig. 4. Based on HMDB classification, 87 metabolites were annotated. Lipids and organic acids were the two main affected pathway superclasses. In the lipids superclass, 59 metabolites were annotated and assigned to four pathway classes with a large section of fatty acyls (40 metabolites, mainly fatty acids). Steroids were the second class of affected metabolites, with 16 metabolites annotated. Among them, five bile acids were disrupted, and 11 more metabolites distributed in five pathway subclasses of steroids such as androstane and pregnane. Finally, three metabolites were annotated as glycerophospholipids and flavonoids. In the organic acids superclass, 23 metabolites were annotated and were assigned to three pathways. The class that was most affected was carboxylic acids and derivatives, with 19 annotated metabolites mainly belonging to amino acids. The last three metabolites belonged to three other classes. Finally, four other superclasses grouping five compounds, including two nucleotides, were affected.

The responses of the main subclasses, i.e. fatty acids, fatty amides, bile acids, steroid derivatives and amino acids, were further analyzed regarding response trends and BMCs.

### 3.3.2. In-depth pathway determination using metabolic responses

The concentration-response relationship of the pathways was based on trends in the responses of metabolites involved in each pathway subclass and the corresponding BMC Q1.

The responses obtained for the main pathway subclasses, i.e. fatty acids, fatty amides, bile acids, steroid derivatives and amino acids, are presented in the trend plot (Fig. 5A). Overall, biphasic responses were common, with two of the pathway subclasses showing U-shaped curves and all pathway subclasses showing bell-shaped curves. Therefore, because the responses were mostly biphasic, it is crucial to have dose-response tools that can analyze models that are more complex than simple linear responses, especially when a wide exposure (concentration) range is applied (here, for instance, the lowest and highest concentration for each PA differed by a factor of 4000). The bell shape was the most frequently found model and both fatty acids and bile acids had this type of curve almost exclusively. About one third of amino acids and fatty acids showed a monotonic increasing curve. For steroids, monotonic linear decreasing curves, and biphasic U-shaped curves were observed following LAS treatment.

All biological pathway subclasses appeared to be highly sensitive to LAS exposure, as reflected by their BMC Q1 values (Fig. 5B). The BMC Q1 value represents the concentration based on ECDF where 25 % of metabolites in the pathway are affected. The 25 % percentile (Q1) is a threshold previously suggested for BMDs assessed from transcriptomic or metabolomic data [39,49]. The ECDF specific to each pathway subclass are presented in Supplementary Fig. 2. Lipids were the most sensitive metabolites, showing impairment below a threshold of 0.05 µM (BMC Q1). In contrast, amino acids responded above this value. The most sensitive metabolites were mainly associated with the pathways of fatty amides and fatty acids. However, very few metabolites were classified as fatty amides. Bile acids were more sensitive than the other lipids
# 3. Results

## 3.1. Cytotoxic effects of PAs in HepaRG cells

To determine the appropriate PA concentrations for metabolomic experiments, the viability of HepaRG cells after a 48 h treatment was determined using the MTT assay (Supplementary Fig. 1).

LAS was the most cytotoxic compound on HepaRG cells, inducing a decrease in viability of 46 ± 12 % at the highest concentration tested, 50 µM. A 41 ± 3 % decrease in viability was also observed at 200 µM RET, whereas HEL had only a slight effect on the viability of HepaRG cells.

Based on these results, we selected the highest subtoxic concentrations, not inducing more than 30 % toxicity, for the metabolomic experiments: 25 µM for LAS and RET, and 100 µM for HEL.

## 3.2. ENDO and EXO concentration-response characterization

### 3.2.1. Responses of unannotated metabolites

Monotonic (i.e. increase or decrease) and non-monotonic (i.e. U-shaped or bell-shaped pattern) response trends for increasing concentration gradients were examined for each PA, for ENDO and EXO conditions and in ESI+ and ESI- modes (Fig. 2). The responsive metabolites were hereafter examined without annotation (named unannotated metabolites). First, the total number of responsive ENDO metabolites was higher than the EXO metabolites for LAS (2398 and 1127 metabolites, respectively) and HEL (1537 and 252 metabolites, respectively), and almost equivalent for RET (around 450 metabolites). The metabolites following a non-monotonic curve response (yellow U-shaped curve, and green bell-shaped curve, Fig. 2) ranged from 25 % to 75 % of responsive metabolites. For LAS, a substantial number of the responsive metabolites showed a bell-shaped trend (51 % for ENDO and 39 % for EXO) and approximately a quarter showed an increasing trend. For RET, 34 % showed a linear decreasing trend for ENDO, and a U-shaped trend was predominant for EXO. For HEL, 47 % of the ENDO metabolites showed U-shaped curves; fewer metabolites (252) were obtained for EXO, mostly modeled with a linear decreasing curve.

### 3.2.2. BMCs of the endometabolome (ENDO) from unannotated metabolites

For ENDO, for each PA, an ECDF plot with a color gradient coding for the response was constructed from the BMCs of all responsive unannotated metabolites (Fig. 3A). The BMCs of metabolites identified in HepaRG cells following LAS exposure ranged from almost 0 (0.00006) to 20.5 µM for ESI+ and ESI- modes combined, with the first quartile of metabolites having a BMC (BMC Q1) below 0.01 µM for both modes. For RET, a comparable range of BMCs from 0.0007 to 18.5 µM was found for both modes and BMC Q1 values below 0.5 µM and 0.6 µM for ESI+ and ESI- modes, respectively. For HEL, the BMCs ranged from approximately 0.0003 to 80 µM for both modes, and BMC Q1 values below 0.06 µM and 0.09 µM for ESI+ and ESI- modes, respectively.

### 3.2.3. BMCs of the exometabolome (EXO) from unannotated metabolites

For EXO, a cumulative distribution plot for each PA was constructed from the BMCs of all unannotated metabolites (Fig. 3B). The BMCs of metabolites highlighted after LAS exposure ranged from almost 0 (0.00006) to approximately 20 µM for the ESI+ and ESI- modes combined, with BMC Q1 values below approximately 1.2 µM for both modes. The BMC distribution for RET exposure showed a comparable range, from almost 0 (0.0001) to 19.5 µM for ESI+ and ESI- modes with BMC Q1 values below 0.14 µM and 0.007 µM for the ESI+ and ESI- modes, respectively. Higher BMCs were determined for HEL, with a distribution from almost 0 (0.0003) to approximately 73 µM for ESI+ and ESI- modes, respectively, with BMC Q1 values below 0.2 µM and 3.1 µM for the ESI+ and ESI- modes, respectively.

## 3.3. Metabolic pathway responses of HepaRG cells to lasiocarpine (LAS)

### 3.3.1. Metabolites and associated pathways as markers of hepatotoxicity

Compared with RET and HEL, LAS had a higher number of affected unannotated metabolites for ENDO as well as a lower BMC Q1, indicating a higher sensitivity of HepaRG cells in response to LAS exposure. For these reasons, the data of unannotated endometabolites following LAS treatment (LAS ENDO) were selected for tentative annotation and BMC pathway assignment.

Significant modeled and annotated metabolites associated with biological pathways disrupted in response to 48 h treatment with a LAS concentration range are reported in Table 1. Additional information on annotated metabolites (m/z, ESI mode, adduct, delta mass, MW, proposed formula, database attribution, main fragments, mzCloud scoring) are provided in Supplementary Table 1. Some fitted concentration-response curves are presented in Fig. 4. Based on HMDB classification, 87 metabolites were annotated. Lipids and organic acids were the two main affected pathway superclasses. In the lipids superclass, 59 metabolites were annotated and assigned to four pathway classes with a large section of fatty acyls (40 metabolites, mainly fatty acids). Steroids were the second class of affected metabolites, with 16 metabolites annotated. Among them, five bile acids were disrupted, and 11 more metabolites distributed in five pathway subclasses of steroids such as androstane and pregnane. Finally, three metabolites were annotated as glycerophospholipids and flavonoids. In the organic acids superclass, 23 metabolites were annotated and were assigned to three pathways. The class that was most affected was carboxylic acids and derivatives, with 19 annotated metabolites mainly belonging to amino acids. The last three metabolites belonged to three other classes. Finally, four other superclasses grouping five compounds, including two nucleotides, were affected.

The responses of the main subclasses, i.e. fatty acids, fatty amides, bile acids, steroid derivatives and amino acids, were further analyzed regarding response trends and BMCs.

### 3.3.2. In-depth pathway determination using metabolic responses

The concentration-response relationship of the pathways was based on trends in the responses of metabolites involved in each pathway subclass and the corresponding BMC Q1.

The responses obtained for the main pathway subclasses, i.e. fatty acids, fatty amides, bile acids, steroid derivatives and amino acids, are presented in the trend plot (Fig. 5A). Overall, biphasic responses were common, with two of the pathway subclasses showing U-shaped curves and all pathway subclasses showing bell-shaped curves. Therefore, because the responses were mostly biphasic, it is crucial to have dose-response tools that can analyze models that are more complex than simple linear responses, especially when a wide exposure (concentration) range is applied (here, for instance, the lowest and highest concentration for each PA differed by a factor of 4000). The bell shape was the most frequently found model and both fatty acids and bile acids had this type of curve almost exclusively. About one third of amino acids and fatty acids showed a monotonic increasing curve. For steroids, monotonic linear decreasing curves, and biphasic U-shaped curves were observed following LAS treatment.

All biological pathway subclasses appeared to be highly sensitive to LAS exposure, as reflected by their BMC Q1 values (Fig. 5B). The BMC Q1 value represents the concentration based on ECDF where 25 % of metabolites in the pathway are affected. The 25 % percentile (Q1) is a threshold previously suggested for BMDs assessed from transcriptomic or metabolomic data [39,49]. The ECDF specific to each pathway subclass are presented in Supplementary Fig. 2. Lipids were the most sensitive metabolites, showing impairment below a threshold of 0.05 µM (BMC Q1). In contrast, amino acids responded above this value. The most sensitive metabolites were mainly associated with the pathways of fatty amides and fatty acids. However, very few metabolites were classified as fatty amides. Bile acids were more sensitive than the other lipids
A

<figure><img src="image_7.png" /><figcaption>**Fig. 3A.** Cumulative distributions of the benchmark concentration (BMC) values for the endometabolites provided and exported from the DRomics tool, with color gradient coding for the response along the concentration. The x-axis represents the BMC values (in µM). The y-axis represents the empirical cumulative distribution function (ECDF) (cumulative proportions of metabolites). Very low estimated BMC values, lower than a threshold defined as the smallest non-null tested concentration divided by 100, were set to this threshold. The signal is scaled between -1 (blue) and +1 (red) to focus on the shape of each response rather than on its amplitude.</figcaption></figure>
A

<figure><img src="image_6.png" /><figcaption>**Fig. 3A.** Cumulative distributions of the benchmark concentration (BMC) values for the endometabolites provided and exported from the DRomics tool, with color gradient coding for the response along the concentration. The x-axis represents the BMC values (in µM). The y-axis represents the empirical cumulative distribution function (ECDF) (cumulative proportions of metabolites). Very low estimated BMC values, lower than a threshold defined as the smallest non-null tested concentration divided by 100, were set to this threshold. The signal is scaled between -1 (blue) and +1 (red) to focus on the shape of each response rather than on its amplitude.</figcaption></figure>
B

<figure><img src="image_8.png" /><figcaption>**Fig. 3B.** Cumulative distributions of the BMC values for the exometabolites provided and exported from the DRomics tool. The x-axis represents the BMC values (in µM). The y-axis represents the empirical cumulative distribution function (ECDF) (as a proportion of metabolites). Very low estimated BMC values, lower than a threshold defined as the smallest non-null tested concentration divided by 100, were set to this threshold. The signal is scaled between -1 (blue) and +1 (red).</figcaption></figure>
B

<figure><img src="image_7.png" /><figcaption>**Fig. 3B.** Cumulative distributions of the BMC values for the exometabolites provided and exported from the DRomics tool. The x-axis represents the BMC values (in µM). The y-axis represents the empirical cumulative distribution function (ECDF) (as a proportion of metabolites). Very low estimated BMC values, lower than a threshold defined as the smallest non-null tested concentration divided by 100, were set to this threshold. The signal is scaled between -1 (blue) and +1 (red).</figcaption></figure>
<caption>**Table 1**</caption>

Main pathway superclasses, classes and subclasses (HMDB classification) corresponding to proposed metabolites (lasiocarpine (LAS) endometabolome (ENDO), ESI+ and ESI- modes), associated BMC₁SD values, models, trends, and adjusted probability values (p-val).

<table><thead><tr><th>Pathway super-class</th><th>Pathway class</th><th>Pathway sub-class</th><th>Putative name</th><th>BMC<sub>1 SD</sub> (µM)</th><th>Model</th><th>Trend</th><th>adjpvalue</th></tr></thead><tbody><tr><td rowspan="15">Lipids and lipid-like molecules</td><td rowspan="15">Fatty Acyls</td><td rowspan="15">Fatty acids and conjugates</td><td>2-Amino-tetradecanoic acid</td><td>0.0001</td><td>log-Gauss-probit</td><td>bell</td><td>2.6E-05</td></tr><tr><td>5-Oxo-6E,8Z-tetradecadienoic acid</td><td>0.0006</td><td>log-Gauss-probit</td><td>bell</td><td>3.5E-05</td></tr><tr><td>Elaidic acid</td><td>0.013</td><td>log-Gauss-probit</td><td>bell</td><td>2.7E-04</td></tr><tr><td>Tetranor-12-HETE</td><td>10.20</td><td>linear</td><td>inc</td><td>9.3E-04</td></tr><tr><td>4-Hydroxy-tridecanoic acid</td><td>0.001</td><td>log-Gauss-probit</td><td>U</td><td>1.3E-03</td></tr><tr><td>8-Amino-7-oxononanoic acid</td><td>0.031</td><td>log-Gauss-probit</td><td>bell</td><td>1.4E-03</td></tr><tr><td>Heptadecanoic acid</td><td>0.019</td><td>log-Gauss-probit</td><td>bell</td><td>1.1E-02</td></tr><tr><td>Palmitelaidic acid</td><td>14.280</td><td>linear</td><td>inc</td><td>2.6E-02</td></tr><tr><td>Linoleic acid</td><td>0.720</td><td>log-Gauss-probit</td><td>bell</td><td>2.8E-02</td></tr><tr><td>Palmitic acid</td><td>17.28</td><td>linear</td><td>inc</td><td>3.5E-02</td></tr><tr><td>2-Hydroxystearate</td><td>19.71</td><td>linear</td><td>inc</td><td>4.6E-02</td></tr><tr><td>3-Carboxy-4-methyl-5-propyl-2-furanpropionate</td><td>20.53</td><td>linear</td><td>dec</td><td>3.8E-02</td></tr><tr><td>Sebacic acid</td><td>0.002</td><td>log-Gauss-probit</td><td>bell</td><td>9.8E-05</td></tr><tr><td>6,9-Dioxo-decanoic acid</td><td>2.102</td><td>Gauss-probit</td><td>bell</td><td>1.0E-04</td></tr><tr><td>Asperitaconic acid A</td><td>0.012</td><td>log-Gauss-probit</td><td>bell</td><td>1.7E-03</td></tr><tr><td>8E-Heptadecenedioic acid</td><td>10.36</td><td>linear</td><td>inc</td><td>1.9E-03</td></tr><tr><td>N-Palmitoyl alanine</td><td>0.006</td><td>log-Gauss-probit</td><td>bell</td><td>2.3E-03</td></tr><tr><td>N-Palmitoyl glutamine</td><td>0.007</td><td>log-Gauss-probit</td><td>bell</td><td>3.6E-03</td></tr><tr><td>Isopalmitic acid</td><td>0.015</td><td>log-Gauss-probit</td><td>bell</td><td>3.6E-03</td></tr><tr><td>Pantheric acid C</td><td>0.001</td><td>log-Gauss-probit</td><td>bell</td><td>6.5E-03</td></tr><tr><td>Suberic acid</td><td>0.003</td><td>log-Gauss-probit</td><td>bell</td><td>1.0E-02</td></tr><tr><td rowspan="15">Fatty amides</td><td rowspan="15">Fatty acids and conjugates</td><td>2E,4E,6E-Nonatrienal</td><td>19.32</td><td>linear</td><td>dec</td><td>2.0E-02</td></tr><tr><td>N-Oleoyl glutamine</td><td>0.007</td><td>log-Gauss-probit</td><td>bell</td><td>2.0E-02</td></tr><tr><td>6,9,12,15,18,21-Tetracosahexaynoic acid</td><td>0.053</td><td>log-Gauss-probit</td><td>bell</td><td>2.1E-02</td></tr><tr><td>Collimonin C</td><td>0.029</td><td>log-Gauss-probit</td><td>U</td><td>2.3E-02</td></tr><tr><td>Aleprestic acid</td><td>0.002</td><td>log-Gauss-probit</td><td>bell</td><td>2.4E-02</td></tr><tr><td>3Z-Dodecenedioic acid</td><td>0.002</td><td>log-Gauss-probit</td><td>bell</td><td>2.4E-02</td></tr><tr><td>Mycinonic acid III</td><td>0.003</td><td>log-Gauss-probit</td><td>U</td><td>2.7E-02</td></tr><tr><td>Ambrettolic acid</td><td>15.39</td><td>linear</td><td>inc</td><td>2.9E-02</td></tr><tr><td>Undecylenic acid</td><td>0.001</td><td>log-Gauss-probit</td><td>bell</td><td>3.5E-02</td></tr><tr><td>16-Hydroxypalmitate</td><td>19.32</td><td>linear</td><td>inc</td><td>3.7E-02</td></tr><tr><td>Oleamide</td><td>0.001</td><td>log-Gauss-probit</td><td>bell</td><td>3.5E-08</td></tr><tr><td>N-Decanoyl-homoserine lactone</td><td>0.000</td><td>log-Gauss-probit</td><td>bell</td><td>5.8E-04</td></tr><tr><td>Decanamide</td><td>0.001</td><td>log-Gauss-probit</td><td>bell</td><td>1.0E-03</td></tr><tr><td>Palmitamide</td><td>0.011</td><td>log-Gauss-probit</td><td>bell</td><td>1.5E-02</td></tr><tr><td>Lineolic acids and derivatives</td><td>9-HODE</td><td>0.009</td><td>log-Gauss-probit</td><td>bell</td><td>2.4E-03</td></tr><tr><td>Fatty acyl glycosides</td><td>Methyl 9-(alpha-D-galactosyloxy)nonanoate</td><td>0.002</td><td>log-Gauss-probit</td><td>bell</td><td>2.4E-08</td></tr><tr><td>Lineolic acids and derivatives</td><td>Linoleic acid</td><td>0.168</td><td>log-Gauss-probit</td><td>bell</td><td>2.1E-02</td></tr><tr><td>Fatty alcohols</td><td>1,2-Eicosanediol</td><td>0.006</td><td>log-Gauss-probit</td><td>bell</td><td>3.9E-04</td></tr><tr><td>Steroids and steroid derivatives</td><td>Bile acids, alcohols and derivatives</td><td>Glycodeoxycholic acid</td><td>0.010</td><td>log-Gauss-probit</td><td>bell</td><td>9.2E-16</td></tr></tbody></table>

(continued on next page)
<caption>**Table 1**</caption>

Main pathway superclasses, classes and subclasses (HMDB classification) corresponding to proposed metabolites (lasiocarpine (LAS) endometabolome (ENDO), ESI+ and ESI- modes), associated BMC<sub>1 SD</sub> values, models, trends, and adjusted probability values (p-val).

<table><thead><tr><th>Pathway super-class</th><th>Pathway class</th><th>Pathway sub-class</th><th>Putative name</th><th>BMC<sub>1 SD</sub> (µM)</th><th>Model</th><th>Trend</th><th>adjpvalue</th></tr></thead><tbody><tr><td rowspan="12">Lipids and lipid-like molecules</td><td rowspan="12">Fatty Acyls</td><td rowspan="12">Fatty acids and conjugates</td><td>2-Amino-tetradecanoic acid</td><td>0.0001</td><td>log-Gauss-probit</td><td>bell</td><td>2.6E-05</td></tr><tr><td>5-Oxo-6E,8Z-tetradecadienoic acid</td><td>0.0006</td><td>log-Gauss-probit</td><td>bell</td><td>3.5E-05</td></tr><tr><td>Elaidic acid</td><td>0.013</td><td>log-Gauss-probit</td><td>bell</td><td>2.7E-04</td></tr><tr><td>Tetranor-12-HETE</td><td>10.20</td><td>linear</td><td>inc</td><td>9.3E-04</td></tr><tr><td>4-Hydroxy-tridecanoic acid</td><td>0.001</td><td>log-Gauss-probit</td><td>U</td><td>1.3E-03</td></tr><tr><td>8-Amino-7-oxononanoic acid</td><td>0.031</td><td>log-Gauss-probit</td><td>bell</td><td>1.4E-03</td></tr><tr><td>Heptadecanoic acid</td><td>0.019</td><td>log-Gauss-probit</td><td>bell</td><td>1.1E-02</td></tr><tr><td>Palmitelaidic acid</td><td>14.280</td><td>linear</td><td>inc</td><td>2.6E-02</td></tr><tr><td>Linoleic acid</td><td>0.720</td><td>log-Gauss-probit</td><td>bell</td><td>2.8E-02</td></tr><tr><td>Palmitic acid</td><td>17.28</td><td>linear</td><td>inc</td><td>3.5E-02</td></tr><tr><td>2-Hydroxystearate</td><td>19.71</td><td>linear</td><td>inc</td><td>4.6E-02</td></tr><tr><td>3-Carboxy-4-methyl-5-propyl-2-furanpropionate</td><td>20.53</td><td>linear</td><td>dec</td><td>3.8E-02</td></tr><tr><td>Sebacic acid</td><td>0.002</td><td>log-Gauss-probit</td><td>bell</td><td>9.8E-05</td></tr><tr><td>6,9-Dioxo-decanoic acid</td><td>2.102</td><td>Gauss-probit</td><td>bell</td><td>1.0E-04</td></tr><tr><td>Asperitaconic acid A</td><td>0.012</td><td>log-Gauss-probit</td><td>bell</td><td>1.7E-03</td></tr><tr><td>8E-Heptadecenedioic acid</td><td>10.36</td><td>linear</td><td>inc</td><td>1.9E-03</td></tr><tr><td>N-Palmitoyl alanine</td><td>0.006</td><td>log-Gauss-probit</td><td>bell</td><td>2.3E-03</td></tr><tr><td>N-Palmitoyl glutamine</td><td>0.007</td><td>log-Gauss-probit</td><td>bell</td><td>3.6E-03</td></tr><tr><td>Isopalmitic acid</td><td>0.015</td><td>log-Gauss-probit</td><td>bell</td><td>3.6E-03</td></tr><tr><td>Pantheric acid C</td><td>0.001</td><td>log-Gauss-probit</td><td>bell</td><td>6.5E-03</td></tr><tr><td>Suberic acid</td><td>0.003</td><td>log-Gauss-probit</td><td>bell</td><td>1.0E-02</td></tr><tr><td rowspan="12">Fatty amides</td><td rowspan="12">Fatty acids and conjugates</td><td>2E,4E,6E-Nonatrienal</td><td>19.32</td><td>linear</td><td>dec</td><td>2.0E-02</td></tr><tr><td>N-Oleoyl glutamine</td><td>0.007</td><td>log-Gauss-probit</td><td>bell</td><td>2.0E-02</td></tr><tr><td>6,9,12,15,18,21-Tetracosahexaynoic acid</td><td>0.053</td><td>log-Gauss-probit</td><td>bell</td><td>2.1E-02</td></tr><tr><td>Collimonin C</td><td>0.029</td><td>log-Gauss-probit</td><td>U</td><td>2.3E-02</td></tr><tr><td>Aleprestic acid</td><td>0.002</td><td>log-Gauss-probit</td><td>bell</td><td>2.4E-02</td></tr><tr><td>3Z-Dodecenedioic acid</td><td>0.002</td><td>log-Gauss-probit</td><td>bell</td><td>2.4E-02</td></tr><tr><td>Mycinonic acid III</td><td>0.003</td><td>log-Gauss-probit</td><td>U</td><td>2.7E-02</td></tr><tr><td>Ambrettolic acid</td><td>15.39</td><td>linear</td><td>inc</td><td>2.9E-02</td></tr><tr><td>Undecylenic acid</td><td>0.001</td><td>log-Gauss-probit</td><td>bell</td><td>3.5E-02</td></tr><tr><td>16-Hydroxypalmitate</td><td>19.32</td><td>linear</td><td>inc</td><td>3.7E-02</td></tr><tr><td>Oleamide</td><td>0.001</td><td>log-Gauss-probit</td><td>bell</td><td>3.5E-08</td></tr><tr><td>N-Decanoyl-homoserine lactone</td><td>0.000</td><td>log-Gauss-probit</td><td>bell</td><td>5.8E-04</td></tr><tr><td>Decanamide</td><td>0.001</td><td>log-Gauss-probit</td><td>bell</td><td>1.0E-03</td></tr><tr><td>Palmitamide</td><td>0.011</td><td>log-Gauss-probit</td><td>bell</td><td>1.5E-02</td></tr><tr><td>Lineolic acids and derivatives</td><td>9-HODE</td><td>0.009</td><td>log-Gauss-probit</td><td>bell</td><td>2.4E-03</td></tr><tr><td>Fatty acyl glycosides</td><td>Methyl 9-(alpha-D-galactosyloxy)nonanoate</td><td>0.002</td><td>log-Gauss-probit</td><td>bell</td><td>2.4E-08</td></tr><tr><td>Lineolic acids and derivatives</td><td>Linoleic acid</td><td>0.168</td><td>log-Gauss-probit</td><td>bell</td><td>2.1E-02</td></tr><tr><td>Fatty alcohols</td><td>1,2-Eicosanediol</td><td>0.006</td><td>log-Gauss-probit</td><td>bell</td><td>3.9E-04</td></tr><tr><td>Steroids and steroid derivatives</td><td>Bile acids, alcohols and derivatives</td><td>Glycodeoxycholic acid</td><td>0.010</td><td>log-Gauss-probit</td><td>bell</td><td>9.2E-16</td></tr></tbody></table>

(continued on next page)
<table><caption>Table 1 (continued)</caption><thead><tr><th>Pathway super-class</th><th>Pathway class</th><th>Pathway sub-class</th><th>Putative name</th><th>BMC<sub>1</sub><sup>SD</sup><br>(µM)</th><th>Model</th><th>Trend</th><th>adjpvalue</th></tr></thead><tbody><tr><td rowspan="12"></td><td rowspan="12"></td><td rowspan="12"></td><td>7-ketodeoxycholic acid</td><td>0.010</td><td>log-Gauss-probit</td><td>bell</td><td>2.4E-02</td></tr><tr><td>7-Alpha-Hydroxy-3-oxo-4-cholestenoate</td><td>3.822</td><td>Gauss-probit</td><td>bell</td><td>3.3E-07</td></tr><tr><td>Glycochenodeoxycholic acid</td><td>8.104</td><td>linear</td><td>dec</td><td>2.0E-05</td></tr><tr><td>Taurocholic acid</td><td>2.344</td><td>log-Gauss-probit</td><td>bell</td><td>1.2E-04</td></tr><tr><td>Sterols</td><td></td><td>Sarmentogenin</td><td>1.338</td><td>Gauss-probit</td><td>bell</td><td>1.0E-07</td></tr><tr><td rowspan="4">Hydroxysteroids</td><td>Pinnasterol</td><td></td><td>Gauss-probit</td><td>dec</td><td>3.0E-16</td></tr><tr><td>Corticosterone</td><td></td><td>Gauss-probit</td><td>U</td><td>1.8E-05</td></tr><tr><td>Cortisone</td><td></td><td>Gauss-probit</td><td>dec</td><td>1.1E-05</td></tr><tr><td>11-Dehydrocorticosterone</td><td>0.009</td><td>log-Gauss-probit</td><td>U</td><td>3.1E-04</td></tr><tr><td></td><td>18-Oxocortisol</td><td></td><td>log-Gauss-probit</td><td>dec</td><td>3.2E-13</td></tr><tr><td>Pregnan e steroids</td><td></td><td>Megestrol</td><td>3.529</td><td>Gauss-probit</td><td>U</td><td>2.7E-05</td></tr><tr><td>Androstane steroids</td><td></td><td>11-Hydroxyandrosterone</td><td>0.002</td><td>log-Gauss-probit</td><td>bell</td><td>1.2E-06</td></tr><tr><td rowspan="3"></td><td rowspan="3"></td><td>11-Ketotestosterone</td><td>18.32</td><td>linear</td><td>dec</td><td>4.1E-02</td></tr><tr><td>Androstenedione</td><td>3.415</td><td>Gauss-probit</td><td>U</td><td>3.1E-04</td></tr><tr><td>2-Hydroxyestrone</td><td>0.001</td><td>log-Gauss-probit</td><td>U</td><td>6.2E-08</td></tr><tr><td>Flavonoids</td><td></td><td>4'-Hydroxy-4-(4-hydroxystyryl)-7-methoxyflavan</td><td>0.039</td><td>log-Gauss-probit</td><td>bell</td><td>5.4E-03</td></tr><tr><td rowspan="2">Glycerophospholipids</td><td rowspan="2"></td><td>1-(2-Methoxy-octadecanyl)-sn-glycero-3-phosphoserine</td><td>0.005</td><td>log-Gauss-probit</td><td>bell</td><td>1.2E-03</td></tr><tr><td>PS(18:2(9Z,12Z)/0:0)</td><td>0.007</td><td>log-Gauss-probit</td><td>bell</td><td>1.0E-03</td></tr><tr><td rowspan="18">Organic acids and derivatives</td><td rowspan="18">Carboxylic acids and derivatives</td><td rowspan="18">Amino acids, peptides and analogues</td><td>Leucylalanine</td><td>3.765</td><td>Gauss-probit</td><td>inc</td><td>4.8E-11</td></tr><tr><td>Leucylserine</td><td>0.055</td><td>log-Gauss-probit</td><td>bell</td><td>1.9E-04</td></tr><tr><td>Valylvaline</td><td>0.072</td><td>log-Gauss-probit</td><td>bell</td><td>1.7E-03</td></tr><tr><td>Prolylvaline</td><td>0.000</td><td>log-Gauss-probit</td><td>bell</td><td>2.7E-03</td></tr><tr><td>Alanyltyrosine</td><td>0.118</td><td>log-Gauss-probit</td><td>bell</td><td>5.3E-03</td></tr><tr><td>Phosphocreatine</td><td>9.598</td><td>linear</td><td>dec</td><td>8.3E-03</td></tr><tr><td>Methionine</td><td>4.032</td><td>Gauss-probit</td><td>inc</td><td>8.6E-03</td></tr><tr><td>Phenylalanylhistidine</td><td>0.573</td><td>log-Gauss-probit</td><td>bell</td><td>1.2E-02</td></tr><tr><td>Phenylalanylvaline</td><td>0.104</td><td>log-Gauss-probit</td><td>bell</td><td>2.3E-02</td></tr><tr><td>N-Acetylglutamic acid</td><td>18.87</td><td>linear</td><td>inc</td><td>4.9E-02</td></tr><tr><td>Nicotinuric acid</td><td>3.344</td><td>Gauss-probit</td><td>dec</td><td>2.3E-12</td></tr><tr><td>Arginylphenylalanine</td><td>0.001</td><td>log-Gauss-probit</td><td>bell</td><td>1.3E-05</td></tr><tr><td>Phenylalanylalanine</td><td>0.049</td><td>log-Gauss-probit</td><td>bell</td><td>6.0E-04</td></tr><tr><td>Gamma-Glutamylcysteine</td><td>0.067</td><td>exponential</td><td>inc</td><td>9.5E-04</td></tr><tr><td>Leucylleucine</td><td>0.044</td><td>log-Gauss-probit</td><td>bell</td><td>2.4E-03</td></tr><tr><td>N-Acetylglutamine</td><td>11.12</td><td>linear</td><td>dec</td><td>7.4E-03</td></tr><tr><td>N-Acetylleucine</td><td>0.007</td><td>log-Gauss-probit</td><td>bell</td><td>6.3E-06</td></tr><tr><td>Capryloylglycine</td><td>10.49</td><td>linear</td><td>dec</td><td>7.7E-04</td></tr><tr><td>Alanylglutamate</td><td>1.983</td><td>log-Gauss-probit</td><td>dec</td><td>2.3E-07</td></tr><tr><td></td><td></td><td></td><td>Cis or trans-Aconitic acid</td><td>16.82</td><td>linear</td><td>inc</td><td>3.5E-02</td></tr><tr><td></td><td>Carboximidic acids and derivatives</td><td>Carboximidic acids</td><td>N1-Acetylspermidine or N8-</td><td>3.041</td><td>Gauss-probit</td><td>U</td><td>3.2E-06</td></tr><tr><td></td><td>Organooxygen compounds</td><td>Carbohydrates and conjugates</td><td>D-Ribonic acid</td><td>0.227</td><td>log-Gauss-probit</td><td>U</td><td>8.3E-03</td></tr><tr><td></td><td>Hydroxy acids and derivatives</td><td>Medium-chain hydroxy acids and derivatives</td><td>2-Hydroxydecanedioic acid</td><td>12.88</td><td>linear</td><td>dec</td><td>4.3E-03</td></tr><tr><td></td><td></td><td>Beta hydroxy acids and derivatives</td><td>Malic acid</td><td>3.8</td><td>Gauss-probit</td><td>inc</td><td>4.8E-03</td></tr><tr><td>Organic oxygen compounds</td><td>Organooxygen compounds</td><td>Carbonyl compounds</td><td>3-Ketosphingosine</td><td>0.008</td><td>log-Gauss-probit</td><td>bell</td><td>2.7E-02</td></tr><tr><td>Organic nitrogen compounds</td><td>Organonitrogen compounds</td><td>Amines</td><td>Hexadecasphingosine</td><td>0.0004</td><td>log-Gauss-probit</td><td>bell</td><td>7.5E-04</td></tr></tbody></table>

(continued on next page)
<table><caption>Table 1 (continued)</caption><thead><tr><th>Pathway super-class</th><th>Pathway class</th><th>Pathway sub-class</th><th>Putative name</th><th>BMC<sub>1 SD</sub><br>(µM)</th><th>Model</th><th>Trend</th><th>adjpvalue</th></tr></thead><tbody><tr><td rowspan="12"></td><td rowspan="12"></td><td rowspan="12"></td><td>7-ketodeoxycholic acid</td><td>0.010</td><td>log-Gauss-probit</td><td>bell</td><td>2.4E-02</td></tr><tr><td>7-Alpha-Hydroxy-3-oxo-4-cholestenoate</td><td>3.822</td><td>Gauss-probit</td><td>bell</td><td>3.3E-07</td></tr><tr><td>Glycochenodeoxycholic acid</td><td>8.104</td><td>linear</td><td>dec</td><td>2.0E-05</td></tr><tr><td>Taurocholic acid</td><td>2.344</td><td>log-Gauss-probit</td><td>bell</td><td>1.2E-04</td></tr><tr><td>Sterols</td><td></td><td>Sarmentogenin</td><td>1.338</td><td>Gauss-probit</td><td>bell</td><td>1.0E-07</td></tr><tr><td rowspan="4">Hydroxysteroids</td><td>Pinnasterol</td><td></td><td>Gauss-probit</td><td>dec</td><td>3.0E-16</td></tr><tr><td>Corticosterone</td><td></td><td>Gauss-probit</td><td>U</td><td>1.8E-05</td></tr><tr><td>Cortisone</td><td></td><td>Gauss-probit</td><td>dec</td><td>1.1E-05</td></tr><tr><td>11-Dehydrocorticosterone</td><td>0.009</td><td>log-Gauss-probit</td><td>U</td><td>3.1E-04</td></tr><tr><td></td><td>18-Oxocortisol</td><td></td><td>log-Gauss-probit</td><td>dec</td><td>3.2E-13</td></tr><tr><td>Pregnan e steroids</td><td></td><td>Megestrol</td><td>3.529</td><td>Gauss-probit</td><td>U</td><td>2.7E-05</td></tr><tr><td>Androstane steroids</td><td></td><td>11-Hydroxyandrosterone</td><td>0.002</td><td>log-Gauss-probit</td><td>bell</td><td>1.2E-06</td></tr><tr><td rowspan="3">Flavonoids</td><td></td><td>11-Ketotestosterone</td><td>18.32</td><td>linear</td><td>dec</td><td>4.1E-02</td></tr><tr><td>Androstenedione</td><td></td><td>Gauss-probit</td><td>U</td><td>3.1E-04</td></tr><tr><td>2-Hydroxyestrone</td><td>0.001</td><td>log-Gauss-probit</td><td>U</td><td>6.2E-08</td></tr><tr><td rowspan="2">Glycerophospholipids</td><td>Flavonoids</td><td></td><td>4'-Hydroxy-4-(4-hydroxystyryl)-7-methoxyflavan</td><td>0.039</td><td>log-Gauss-probit</td><td>bell</td><td>5.4E-03</td></tr><tr><td>Glycerophosphoserines</td><td></td><td>1-(2-Methoxy-octadecanyl)-sn-glycero-3-phosphoserine</td><td>0.005</td><td>log-Gauss-probit</td><td>bell</td><td>1.2E-03</td></tr><tr><td rowspan="22">Organic acids and derivatives</td><td rowspan="22"></td><td>Carboxylic acids and derivatives</td><td></td><td>PS(18:2(9Z,12Z)/0:0)</td><td>0.007</td><td>log-Gauss-probit</td><td>bell</td><td>1.0E-03</td></tr><tr><td>Organooxygen compounds</td><td></td><td>Leucylalanine</td><td>3.765</td><td>Gauss-probit</td><td>inc</td><td>4.8E-11</td></tr><tr><td rowspan="12">Carboximicids and derivatives</td><td rowspan="12">Organooxygen compounds</td><td>Carboximicids</td><td></td><td>Leucylserine</td><td>0.055</td><td>log-Gauss-probit</td><td>bell</td><td>1.9E-04</td></tr><tr><td>Carbohydrates and conjugates</td><td></td><td>Valylvaline</td><td>0.072</td><td>log-Gauss-probit</td><td>bell</td><td>1.7E-03</td></tr><tr><td>Medium-chain hydroxy acids and derivatives</td><td></td><td>Prolylvaline</td><td>0.000</td><td>log-Gauss-probit</td><td>bell</td><td>2.7E-03</td></tr><tr><td>Beta hydroxy acids and derivatives</td><td></td><td>Alanyltyrosine</td><td>0.118</td><td>log-Gauss-probit</td><td>bell</td><td>5.3E-03</td></tr><tr><td>Carbonyl compounds</td><td></td><td>Phosphocreatine</td><td>9.598</td><td>linear</td><td>dec</td><td>8.3E-03</td></tr><tr><td rowspan="2">Organonitrogen compounds</td><td rowspan="2">Organonitrogen compounds</td><td>Methionine</td><td>4.032</td><td>Gauss-probit</td><td>inc</td><td>8.6E-03</td></tr><tr><td>Phenylalanylhistidine</td><td>0.573</td><td>log-Gauss-probit</td><td>bell</td><td>1.2E-02</td></tr><tr><td rowspan="12">Organic nitrogen compounds</td><td rowspan="12">Organic nitrogen compounds</td><td>Phenylalanylvaline</td><td>0.104</td><td>log-Gauss-probit</td><td>bell</td><td>2.3E-02</td></tr><tr><td>N-Acetylglutamic acid</td><td>18.87</td><td>linear</td><td>inc</td><td>4.9E-02</td></tr><tr><td>Nicotinuric acid</td><td>3.344</td><td>Gauss-probit</td><td>dec</td><td>2.3E-12</td></tr><tr><td>Arginylphenylalanine</td><td>0.001</td><td>log-Gauss-probit</td><td>bell</td><td>1.3E-05</td></tr><tr><td>Phenylalanylalanine</td><td>0.049</td><td>log-Gauss-probit</td><td>bell</td><td>6.0E-04</td></tr><tr><td>Gamma-Glutamylcysteine</td><td>0.067</td><td>exponential</td><td>inc</td><td>9.5E-04</td></tr><tr><td>Leucylleucine</td><td>0.044</td><td>log-Gauss-probit</td><td>bell</td><td>2.4E-03</td></tr><tr><td>N-Acetylglutamine</td><td>11.12</td><td>linear</td><td>dec</td><td>7.4E-03</td></tr><tr><td>N-Acetylleucine</td><td>0.007</td><td>log-Gauss-probit</td><td>bell</td><td>6.3E-06</td></tr><tr><td>Capryloylglycine</td><td>10.49</td><td>linear</td><td>dec</td><td>7.7E-04</td></tr><tr><td>Amines</td><td></td><td>Alanylglutamate</td><td>1.983</td><td>log-Gauss-probit</td><td>dec</td><td>2.3E-07</td></tr><tr><td>Tricarboxylic acids and derivatives</td><td></td><td>Cis or trans-Aconitic acid</td><td>16.82</td><td>linear</td><td>inc</td><td>3.5E-02</td></tr><tr><td>Carboximicids and derivatives</td><td>Carboximicids</td><td></td><td>N1-Acetylspermidine or N8-</td><td>3.041</td><td>Gauss-probit</td><td>U</td><td>3.2E-06</td></tr><tr><td>Organooxygen compounds</td><td>Carbohydrates and conjugates</td><td></td><td>D-Ribonic acid</td><td>0.227</td><td>log-Gauss-probit</td><td>U</td><td>8.3E-03</td></tr><tr><td>Hydroxy acids and derivatives</td><td>Medium-chain hydroxy acids and derivatives</td><td></td><td>2-Hydroxydecanedioic acid</td><td>12.88</td><td>linear</td><td>dec</td><td>4.3E-03</td></tr><tr><td></td><td>Beta hydroxy acids and derivatives</td><td></td><td>Malic acid</td><td>3.8</td><td>Gauss-probit</td><td>inc</td><td>4.8E-03</td></tr><tr><td>Organic oxygen compounds</td><td>Organooxygen compounds</td><td>Carbonyl compounds</td><td>3-Ketosphingosine</td><td>0.008</td><td>log-Gauss-probit</td><td>bell</td><td>2.7E-02</td></tr><tr><td></td><td>Organonitrogen compounds</td><td>Amines</td><td>Hexadecasphingosine</td><td>0.0004</td><td>log-Gauss-probit</td><td>bell</td><td>7.5E-04</td></tr></tbody></table>

(continued on next page)
<table><caption>**Table 1 (continued)**</caption><thead><tr><th>Pathway super-class</th><th>Pathway class</th><th>Pathway sub-class</th><th>Putative name</th><th>BMC<sub>1</sub> SD (µM)</th><th>Model</th><th>Trend</th><th>adjpvalue</th></tr></thead><tbody><tr><td>Organoheterocyclic compounds</td><td>Biotin and derivatives</td><td>NA</td><td>Biotin</td><td>16.66</td><td>linear</td><td>inc</td><td>2.0E-02</td></tr><tr><td rowspan="2">Nucleosides, nucleotides and analogues</td><td rowspan="2">Pyrimidine nucleosides</td><td>Pyrimidine 2'-deoxyribonucleosides</td><td>Deoxycytidine</td><td>0.2777</td><td>log-Gauss-probit</td><td>bell</td><td>3.3E-03</td></tr><tr><td>Pyrimidine ribonucleotides</td><td>UTP</td><td>0.0044</td><td>log-Gauss-probit</td><td>U</td><td>2.8E-06</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Fig. 4.** Fitted response-concentration curves for four annotated metabolites.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 5.** Trend plot (A) and sensitivity plot (B) based on the first quartile (Q1) of BMD values associated with each pathway subclass for the lasiocarpine endometabolome.</figcaption></figure>

at 0.1 µM. The underlying steroid pathway was affected at a threshold (BMC Q1) of 0.4 µM.

### 3.4. BMC comparison for lasiocarpine (LAS)

The BMCs for LAS obtained from metabolomics in this study and those retrieved in the literature from genotoxicity potencies measured in the HepaRG cell line are summarized in Table 2. The BMCs from unannotated metabolites and their confidence interval are provided in Supplementary Fig. 3. BMC Q1 values for unannotated and annotated metabolites in our metabolomic experiment were comparable for ENDO. For annotated endometabolites, the lipid pathway was the most sensitive, with a BMC median at 0.01 µM. The BMC of genotoxic potencies using the micronucleus assay was around 1 µM, which does not differ

much from the BMC median of the organic acid pathway calculated in this study. Another in vitro genotoxicity assay using γH2AX induction generated a high BMC value at 5.3 µM, closer to our BMC median value of 3.22 µM for unannotated extracellular metabolites (EXO) than for our BMC median value of 0.29 µM for unannotated intracellular metabolites (ENDO).

## 4. Discussion

### 4.1. Comparison of the sensitivity of the endometabolome (ENDO) and exometabolome (EXO) responses

In this study, we assessed the metabolomic responses of HepaRG cells following 48 h of exposure to increasing concentrations of LAS, RET, or
<table><caption>**Table 1 (continued)**</caption><thead><tr><th>Pathway super-class</th><th>Pathway class</th><th>Pathway sub-class</th><th>Putative name</th><th>BMC<sub>1</sub> SD (µM)</th><th>Model</th><th>Trend</th><th>adjpvalue</th></tr></thead><tbody><tr><td>Organoheterocyclic compounds</td><td>Biotin and derivatives</td><td>NA</td><td>Biotin</td><td>16.66</td><td>linear</td><td>inc</td><td>2.0E-02</td></tr><tr><td rowspan="2">Nucleosides, nucleotides and analogues</td><td rowspan="2">Pyrimidine nucleosides</td><td>Pyrimidine 2'-deoxyribonucleosides</td><td>Deoxycytidine</td><td>0.2777</td><td>log-Gauss-probit</td><td>bell</td><td>3.3E-03</td></tr><tr><td>Pyrimidine ribonucleotides</td><td>UTP</td><td>0.0044</td><td>log-Gauss-probit</td><td>U</td><td>2.8E-06</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Fig. 4.** Fitted response-concentration curves for four annotated metabolites.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 5.** Trend plot (A) and sensitivity plot (B) based on the first quartile (Q1) of BMD values associated with each pathway subclass for the lasiocarpine endometabolome.</figcaption></figure>

at 0.1 µM. The underlying steroid pathway was affected at a threshold (BMC Q1) of 0.4 µM.

### 3.4. BMC comparison for lasiocarpine (LAS)

The BMCs for LAS obtained from metabolomics in this study and those retrieved in the literature from genotoxicity potencies measured in the HepaRG cell line are summarized in Table 2. The BMCs from unannotated metabolites and their confidence interval are provided in Supplementary Fig. 3. BMC Q1 values for unannotated and annotated metabolites in our metabolomic experiment were comparable for ENDO. For annotated endometabolites, the lipid pathway was the most sensitive, with a BMC median at 0.01 µM. The BMC of genotoxic potencies using the micronucleus assay was around 1 µM, which does not differ

much from the BMC median of the organic acid pathway calculated in this study. Another in vitro genotoxicity assay using γH2AX induction generated a high BMC value at 5.3 µM, closer to our BMC median value of 3.22 µM for unannotated extracellular metabolites (EXO) than for our BMC median value of 0.29 µM for unannotated intracellular metabolites (ENDO).

## 4. Discussion

### 4.1. Comparison of the sensitivity of the endometabolome (ENDO) and exometabolome (EXO) responses

In this study, we assessed the metabolomic responses of HepaRG cells following 48 h of exposure to increasing concentrations of LAS, RET, or
<table><caption>**Table 2**
Benchmark concentrations (BMCs) from different experiments with HepaRG cells exposed to lasiocarpine (LAS).</caption><thead><tr><th>BMC<br>μM</th><th>Unannotated metabolites</th><th>Unannotated metabolites</th><th>Annotated metabolites</th><th>Lipid pathway</th><th>Organic acid pathway</th><th>Literature</th><th>Literature</th></tr><tr><th></th><th>EXO</th><th>ENDO</th><th>ENDO BMC<sub>1</sub> SD</th><th>ENDO</th><th>ENDO</th><th>BMC<sub>50</sub></th><th>BMD<sub>100</sub></th></tr><tr><th></th><th>BMC<sub>1</sub> SD</th><th></th><th>a</th><th></th><th>ENDO BMC<sub>1</sub> SD</th><th>γH2AX</th><th>(BMDL-BMDU)</th></tr></thead><tbody><tr><th>BMC</th><td>/</td><td>/</td><td>/</td><td>/</td><td>/</td><td>5.3</td><td>0.8-1.1</td></tr><tr><th>BMC 25th percentile</th><td>1.23</td><td>0.01</td><td>0.007</td><td>0.003</td><td>0.06</td><td>/</td><td>/</td></tr><tr><th>BMC median</th><td>3.22</td><td>0.29</td><td>0.06</td><td>0.01</td><td>0.57</td><td>/</td><td>/</td></tr></tbody></table>

a This study, results presented in Section 3.2.2 for ENDO and 3.2.3 for EXO

b This study, results from annotated metabolites in Supplementary Table 1

c This study, results from annotated metabolites presented in Table 1, Section 3.3.1

d Louisse et al., 2019 [57] (BMC<sub>50</sub> for Hill model: BMC50; 1.5-fold γH2AX induction compared with the background)

e Allemang et al., 2018 [65] Critical Effect Size at 1, representing a 100 % increase over the background level (doubling of the micronucleus count of the corresponding control). The calculated lower and upper confidence limits (CEDL, CEDU) of the critical effect dose (CED) resulting from the exponential and Hill model 3 were compiled

HEL, and we investigated them at two levels: ENDO (i.e. intracellular metabolites) and EXO (i.e. extracellular metabolites).

To compare and analyze the response trends in ENDO versus EXO, we used unannotated metabolites that gave significant models with DRomics. Our results on unannotated metabolites demonstrated that HepaRG cells are sensitive to all three PAs, as shown by the number of metabolites affected in both ENDO and EXO analyses. Owing to its differentiation and biotransformation capacities and its genetic stability, the human HepaRG cell model has numerous advantages for studying hepatotoxic effects of compounds [47,50,51]. Using unannotated metabolites to retrieve information on cellular responses can be a useful method when the metabolome is poorly annotated. Furthermore, it can be used at a first-tier level to sort compounds in the same family that cause hepatotoxic effects. To investigate temporal trends between chemicals, Malinowska et al. (year) showed that different rank-based approaches on unannotated and annotated features are more or less equivalent when assessing the BMC, varying less than tenfold across the approaches [37]. We applied BMCs based on the ranked first quartile (BMC Q1) value to compare the sensitivities of the metabolome in response to the three PAs. For LAS, the difference between the BMC Q1 values assessed from unannotated and annotated metabolites was less than twofold.

Cellular responses to xenobiotics such as PAs can be revealed by analyzing the intracellular and the extracellular metabolites to identify the perturbations induced in cell metabolism [42]. The intracellular metabolome (ENDO) extracted from cell or tissue content [22,23] is much more frequently studied than the extracellular metabolome (EXO) extracted from cell culture media [36]. However, EXO can provide information on the secretion of key molecules [42], and it is relatively easy to prepare samples from cell culture media. EXO is also relevant for investigating the xenobiotic metabolism and the metabolites produced from the parent compound, such as PAs, particularly when toxifying reactions occur [15,52].

The discrepancy between ENDO and EXO metabolomes regarding unannotated metabolites suggest that ENDO is of primary interest for studying concentration-responses related to hepatotoxicity. In this study, more intracellular metabolites (ENDO) than extracellular metabolites (EXO) were perturbed for at least two PAs: LAS and HEL. For LAS, these metabolites showed higher sensitivity in ENDO than EXO experiments regarding the Q1, and the response curves of the EXO and ENDO metabolites primarily followed a biphasic bell shape, with fewer metabolites following a linear response model. Bell-shaped curves can indeed describe the behavior of metabolites, which do not necessarily have a linear response to increasing concentrations. For this reason, biphasic shapes are included in several tools for the dose-response analysis of high-throughput omics data [45]. They can describe the behavior of metabolites in response to low and high concentrations. For example, adaptive behavior or induced pathways may be elicited at low

concentrations, but inhibited or not elicited at higher concentrations [53]. Overall, biphasic models accounted for at least half of the models, except for HEL EXO for which fewer metabolites were impaired.

## 4.2. Potency of lasiocarpine (LAS) on HepaRG cells compared with HEL and RET

The BMC assessment based on cumulative distribution served as an indicator to classify the hepatotoxic effect of three PAs. As described above, the metabolome of HepaRG cells was most affected by LAS treatment, which led to a higher number of modulated metabolites and cumulative distributions of BMCs generated than the other two PAs. Notably, most ENDO metabolites had a BMC Q1 below 0.01 µM. Calculated thresholds such as BMC/BMD can serve as reference values to classify various compounds depending on their potential to induce liver injury. Consequently, compounds with a higher potential for hepatotoxicity may require in-depth investigation. To identify those of toxicological concern, one study applied concentration-response modeling to transcriptomics data generated after exposure of both differentiated and non-differentiated human HepaRG cells to 24 reference compounds [54]. High-throughput systems are particularly useful for studying concentration-response models of multiple compounds in omics. For example, Reardon et al. applied high-throughput transcriptional profiling to assess the effects of 23 PFASs on human 3D spheroids [55]. In this study, we used a low-throughput approach to ensure a sufficient amount of analyzable material as a proof-of-concept for using metabolomics-derived BMCs. At a second-tier level, this strategy could be further extended to yet more compounds, depending on the sensitivity of the HRMS instrument. The use of such BMCs can then foster read-across features for some families of compounds, based on structure-dependent toxicity. From our results on the ENDO metabolome, we inferred the order of toxicity LAS>HEL>RET based on BMC Q1, i.e. 0.01 µM (LAS ENDO, ESI+ mode) < 0.06 µM (HEL ENDO, ESI+ mode) < 0.5 µM (RET ENDO, ESI+ mode). Therefore, it can be hypothesized that structures with an open-chain diester like LAS cause more effects than other types of structures. Moreover, open-chain monoesters like HEL appear to induce stronger molecular toxic effects on HepaRG cells than cyclic-chain structures like RET, a cyclic diester. Some studies have explored PA-mediated structure toxicity in relation to endpoints. Impairment of bile acid homeostasis and cell death were the two main effects observed in HepaRG cells treated with retronecine-type PAs, regardless of their diester structure (open chain or cyclic chain) [27,56]. The toxicity of the diester-type PAs has been confirmed in other studies [14,57,58]. Out of 22 structurally different PAs, diesters — whether they had open-chain or cyclic-chain structures — had the greatest potential for deregulating the expression of genes associated with cholestasis in HepaRG cells and for decreasing bile acid concentrations [59]. However, in our metabolomic experiment, the monoester
<table><caption>**Table 2**
Benchmark concentrations (BMCs) from different experiments with HepaRG cells exposed to lasiocarpine (LAS).</caption><thead><tr><th>BMC<br>μM</th><th>Unannotated metabolites</th><th>Unannotated metabolites</th><th>Annotated metabolites</th><th>Lipid pathway</th><th>Organic acid pathway</th><th>Literature</th><th>Literature</th></tr><tr><td></td><td>EXO</td><td>ENDO</td><td>ENDO BMC<sub>1</sub> SD</td><td>ENDO</td><td>ENDO</td><td>BMC<sub>50</sub></td><td>BMD<sub>100</sub></td></tr><tr><td></td><td>BMC<sub>1</sub> SD</td><td>BMC<sub>1</sub> SD</td><td></td><td></td><td></td><td>γH2AX</td><td>(BMDL-BMDU)</td></tr><tr><td></td><td><sup>a</sup></td><td><sup>a</sup></td><td><sup>b</sup></td><td><sup>c</sup></td><td><sup>c</sup></td><td><sup>d</sup></td><td>Micronucleus</td></tr></thead><tbody><tr><td>BMC</td><td>/</td><td>/</td><td>/</td><td>/</td><td>/</td><td>5.3</td><td>0.8-1.1</td></tr><tr><td>BMC 25th percentile</td><td>1.23</td><td>0.01</td><td>0.007</td><td>0.003</td><td>0.06</td><td>/</td><td>/</td></tr><tr><td>BMC median</td><td>3.22</td><td>0.29</td><td>0.06</td><td>0.01</td><td>0.57</td><td>/</td><td>/</td></tr></tbody></table>

<sup>a</sup> This study, results presented in Section 3.2.2 for ENDO and 3.2.3 for EXO

<sup>b</sup> This study, results from annotated metabolites in Supplementary Table 1

<sup>c</sup> This study, results from annotated metabolites presented in Table 1, Section 3.3.1

<sup>d</sup> Louisse et al., 2019 [57] (BMC<sub>50</sub> for Hill model: BMC50; 1.5-fold γH2AX induction compared with the background)

<sup>e</sup> Allemang et al., 2018 [65] Critical Effect Size at 1, representing a 100 % increase over the background level (doubling of the micronucleus count of the corresponding control). The calculated lower and upper confidence limits (CEDL, CEDU) of the critical effect dose (CED) resulting from the exponential and Hill model 3 were compiled

HEL, and we investigated them at two levels: ENDO (i.e. intracellular metabolites) and EXO (i.e. extracellular metabolites).

To compare and analyze the response trends in ENDO versus EXO, we used unannotated metabolites that gave significant models with DRomics. Our results on unannotated metabolites demonstrated that HepaRG cells are sensitive to all three PAs, as shown by the number of metabolites affected in both ENDO and EXO analyses. Owing to its differentiation and biotransformation capacities and its genetic stability, the human HepaRG cell model has numerous advantages for studying hepatotoxic effects of compounds [47,50,51]. Using unannotated metabolites to retrieve information on cellular responses can be a useful method when the metabolome is poorly annotated. Furthermore, it can be used at a first-tier level to sort compounds in the same family that cause hepatotoxic effects. To investigate temporal trends between chemicals, Malinowska et al. (year) showed that different rank-based approaches on unannotated and annotated features are more or less equivalent when assessing the BMC, varying less than tenfold across the approaches [37]. We applied BMCs based on the ranked first quartile (BMC Q1) value to compare the sensitivities of the metabolome in response to the three PAs. For LAS, the difference between the BMC Q1 values assessed from unannotated and annotated metabolites was less than twofold.

Cellular responses to xenobiotics such as PAs can be revealed by analyzing the intracellular and the extracellular metabolites to identify the perturbations induced in cell metabolism [42]. The intracellular metabolome (ENDO) extracted from cell or tissue content [22,23] is much more frequently studied than the extracellular metabolome (EXO) extracted from cell culture media [36]. However, EXO can provide information on the secretion of key molecules [42], and it is relatively easy to prepare samples from cell culture media. EXO is also relevant for investigating the xenobiotic metabolism and the metabolites produced from the parent compound, such as PAs, particularly when toxifying reactions occur [15,52].

The discrepancy between ENDO and EXO metabolomes regarding unannotated metabolites suggest that ENDO is of primary interest for studying concentration-responses related to hepatotoxicity. In this study, more intracellular metabolites (ENDO) than extracellular metabolites (EXO) were perturbed for at least two PAs: LAS and HEL. For LAS, these metabolites showed higher sensitivity in ENDO than EXO experiments regarding the Q1, and the response curves of the EXO and ENDO metabolites primarily followed a biphasic bell shape, with fewer metabolites following a linear response model. Bell-shaped curves can indeed describe the behavior of metabolites, which do not necessarily have a linear response to increasing concentrations. For this reason, biphasic shapes are included in several tools for the dose-response analysis of high-throughput omics data [45]. They can describe the behavior of metabolites in response to low and high concentrations. For example, adaptive behavior or induced pathways may be elicited at low

concentrations, but inhibited or not elicited at higher concentrations [53]. Overall, biphasic models accounted for at least half of the models, except for HEL EXO for which fewer metabolites were impaired.

## 4.2. Potency of lasiocarpine (LAS) on HepaRG cells compared with HEL and RET

The BMC assessment based on cumulative distribution served as an indicator to classify the hepatotoxic effect of three PAs. As described above, the metabolome of HepaRG cells was most affected by LAS treatment, which led to a higher number of modulated metabolites and cumulative distributions of BMCs generated than the other two PAs. Notably, most ENDO metabolites had a BMC Q1 below 0.01 µM. Calculated thresholds such as BMC/BMD can serve as reference values to classify various compounds depending on their potential to induce liver injury. Consequently, compounds with a higher potential for hepatotoxicity may require in-depth investigation. To identify those of toxicological concern, one study applied concentration-response modeling to transcriptomics data generated after exposure of both differentiated and non-differentiated human HepaRG cells to 24 reference compounds [54]. High-throughput systems are particularly useful for studying concentration-response models of multiple compounds in omics. For example, Reardon et al. applied high-throughput transcriptional profiling to assess the effects of 23 PFASs on human 3D spheroids [55]. In this study, we used a low-throughput approach to ensure a sufficient amount of analyzable material as a proof-of-concept for using metabolomics-derived BMCs. At a second-tier level, this strategy could be further extended to yet more compounds, depending on the sensitivity of the HRMS instrument. The use of such BMCs can then foster read-across features for some families of compounds, based on structure-dependent toxicity. From our results on the ENDO metabolome, we inferred the order of toxicity LAS>HEL>RET based on BMC Q1, i.e. 0.01 µM (LAS ENDO, ESI+ mode) < 0.06 µM (HEL ENDO, ESI+ mode) < 0.5 µM (RET ENDO, ESI+ mode). Therefore, it can be hypothesized that structures with an open-chain diester like LAS cause more effects than other types of structures. Moreover, open-chain monoesters like HEL appear to induce stronger molecular toxic effects on HepaRG cells than cyclic-chain structures like RET, a cyclic diester. Some studies have explored PA-mediated structure toxicity in relation to endpoints. Impairment of bile acid homeostasis and cell death were the two main effects observed in HepaRG cells treated with retronecine-type PAs, regardless of their diester structure (open chain or cyclic chain) [27,56]. The toxicity of the diester-type PAs has been confirmed in other studies [14,57,58]. Out of 22 structurally different PAs, diesters — whether they had open-chain or cyclic-chain structures — had the greatest potential for deregulating the expression of genes associated with cholestasis in HepaRG cells and for decreasing bile acid concentrations [59]. However, in our metabolomic experiment, the monoester
HEL appeared to be more toxic than the cyclic diester RET, because more ENDO metabolites were affected by HEL than by RET, and the BMC Q1 for HEL was lower. Our results can be attributed to the OH position: LAS (diester) and HEL (monoester) are a heliotridine (7 S) necine-base type, whereas RET (diester) is an aetronecine (7 R) necine-base type.

Given that LAS was the PA that affected the HepaRG metabolome the most and gave a low BMC as well as a higher number of affected ENDO metabolites, we investigated LAS further.

### 4.3. Metabolic pathways impaired by LAS exposure

Our study explored for the first time LAS molecular mechanisms of hepatotoxicity at the intracellular level and associated BMC. The most sensitive pathways impaired by LAS were related to lipid (fatty amides and fatty acids) metabolism (BMC Q1 below 0.01 µM), and then to bile acid metabolism (BMC Q1 at 0.01 µM). Next, many impaired metabolites related to amino acid metabolism were affected (BMC Q1 at 0.05 µM). Finally, the metabolism of other steroid derivatives was impaired at higher concentrations (BMC Q1 at 0.4 µM).

In exposure of the liver to PAs, three main molecular MOAs underlying hepatotoxicity have been described: (i) DNA damage response, (ii) induction of apoptosis by cellular oxidative stress through the intrinsic mitochondrial pathway or excessive production of reactive oxygen species (ROS), and (iii) effects on bile acid homeostasis (Gluck et al., 2020; Wang et al., 2021). Among the pathways that were impaired by LAS, the metabolism of bile acid has been extensively studied in the literature and shown to be greatly affected by different PAs [24,25,27, 59]. In particular, Waizenegger et al. revealed that the hepatotoxicity of four PAs is associated with a significant decrease in bile acids in both the intracellular and extracellular contents of HepaRG cells [27]. In our study, the responses of fatty acids and bile acids followed a bell-shaped trend, indicating an impact on bile acid homeostasis induced by both low and high concentrations of LAS. Indeed, in response to low concentrations of hepatotoxicants, cells may compensate for the mechanistic effects by increasing their lipid content, but this compensation may not occur at higher concentrations of the toxicants because lipid synthesis is also affected [20]. The amino acid metabolism was impaired at higher BMCs than the lipid metabolism, but still at very low concentrations. Like lipids and bile acids, most of the amino acids responded with a bell-shaped trend throughout the exposure range. In support of this observation, primary bile acids are very rapidly conjugated with glycine and taurine after their synthesis [59]. The primary purpose of bile acids is to facilitate digestion of fat through its surfactant properties, with emulsification of fats into micelles [60]. Amino acids, lipids and the bile acid glycocholic acid have also been described as being impaired in a recent in vivo study in the plasma and urine of rats exposed to a decoction of *Gynura segetum*, an herbal medicine containing PAs [26]. The authors identified biomarkers of hepatotoxicity, including some amino acids and lipids, using a metabolic profiling analysis [26]. In HepaRG cells, we also observed impairment of metabolites conjugated with glutamic acid, mainly with a linear increasing trend. This finding can be associated with disruptions in the tricarboxylic acid (TCA) cycle [47], potentially causing oxidative stress, and is further? corroborated by the impaired level of cis- or trans-aconitic acid. Among other impaired amino acids, nicotinuric acid showed a linear decreasing trend, partly explaining the increased production of amino acids involved in glutamate metabolism. In addition, impairment of the mitochondrial function can be revealed by disorders observed in lipids, organic acids and nucleotides, such as modifications in the metabolism of sphingolipids describing a bell-shaped trend, and an effect on UTP production at a very low BMCs. Finally, the steroid class including steroids other than bile acids, such as hydroxysteroids, pregnane, androstane and estrane steroids, were impaired at higher concentrations than lipids, and followed a U-shaped trend. To our knowledge, this is the first evidence that PA may affect steroid metabolism. Their concentration-response trend was contrary to that of bile acids and lipids. Therefore, LAS may act as a

pregnane X receptor (PXR) agonist or antagonist, because this transcription factor is known to act as a physiological sensor in the regulation of bile acid and lipid metabolisms and of steroid and endocrine homeostasis [61]. A previous study showed that LAS displayed PXR-activating properties in PXR-transfected HepG2 cells, and this activation may participate in liver steatosis ([62]; https://aopwiki.org/aops/60).

### 4.4. BMC assessment through metabolomics

The assessment of PODs such as BMCs, derived from concentration-response relationship of in vitro omics data is fully compliant with NAMs. However, there is a need for further knowledge and consensus among the scientific community before it can be useful and fully accepted for use in risk assessment. Genomic and transcriptomic data have been explored for this purpose, and have led to guidelines and recommendations [33,49,63], but application to metabolomics requires further investigation. In this study, we explored the concentration-response trends and derived BMDs for intra- and extracellular metabolites of HepaRG cells produced after treatment with three PAs. The time of exposure (48 h) was deemed appropriate for obtaining a sufficient number of affected metabolites. A previous comparative study indicated that the time points of 24 h and 48 h produced consistent BMC values for three out of the four tested chemicals [37]. We were also able to compare trends in the modeled responses of unannotated metabolites as well as the BMCs from both the exo- and endometabolome. We concluded that, for unannotated metabolites, EXO provided a less sensitive BMC Q1 than ENDO. This result may depend on the class of compounds and their toxicokinetic profile. Moreover, the comparison of different BMCs for LAS shows that lipid metabolism is by far the most sensitive class of metabolites, even compared with BMCs derived for genotoxic potency (Table 2). Nevertheless, BMCs derived from specific pathway or subpathway classes must be supported by a sufficient level of confidence in metabolite annotation, particularly when very few metabolites are annotated in a given (sub)pathway. The metabolites that could be annotated mostly showed L3 confidence levels, with the rest being L2 (mzCloud matching). Confidence levels are also crucial when defining the metabolites that can be used as input for regulatory toxicology [64]. Finally, as for transcriptomics, a framework needs to be established to define, in view of the thousands of metabolites involved, which BMCs must be examined (1st percentile, 25th percentile (Q1), median) so as to be sufficiently protective and comparable to traditional apical PODs [34,49,55].

## 5. Conclusion

Overall, our findings indicate that exposure to PAs via food and beverages may have an impact on health, with hepatotoxic effects observed in HepaRG cells at very low concentrations of PAs in our investigation of impaired metabolism pathways. This study applied concentration-response modeling and BMC derivation to understand and dissociate the hepatotoxic effects of three different PAs. To our knowledge, this is the first study showing the benefit of such an approach on several contaminants from the same family, but with different chemical structures. Using unannotated metabolites revealed that the endometabolome was more sensitive than the exometabolome regarding effects on the metabolome, and that LAS was confirmed as the most hepatotoxic PA. Interestingly, our trend analysis highlighted the relevance of bell-shaped responses as indicators of early adaptive changes. At a second-tier level, the MOA of LAS was enriched by highlighting specific perturbed biochemical pathways such as lipid metabolism, which was the most sensitive impaired pathway, followed by amino acid metabolism. In particular, bile acid deregulation was confirmed, which is consistent with several previous studies; however, our is the first report of impairment of steroid metabolism, necessitating further studies. The response of HepaRG cells to PA exposure led to
HEL appeared to be more toxic than the cyclic diester RET, because more ENDO metabolites were affected by HEL than by RET, and the BMC Q1 for HEL was lower. Our results can be attributed to the OH position: LAS (diester) and HEL (monoester) are a heliotridine (7 S) necine-base type, whereas RET (diester) is an aetronecine (7 R) necine-base type.

Given that LAS was the PA that affected the HepaRG metabolome the most and gave a low BMC as well as a higher number of affected ENDO metabolites, we investigated LAS further.

### 4.3. Metabolic pathways impaired by LAS exposure

Our study explored for the first time LAS molecular mechanisms of hepatotoxicity at the intracellular level and associated BMC. The most sensitive pathways impaired by LAS were related to lipid (fatty amides and fatty acids) metabolism (BMC Q1 below 0.01 µM), and then to bile acid metabolism (BMC Q1 at 0.01 µM). Next, many impaired metabolites related to amino acid metabolism were affected (BMC Q1 at 0.05 µM). Finally, the metabolism of other steroid derivatives was impaired at higher concentrations (BMC Q1 at 0.4 µM).

In exposure of the liver to PAs, three main molecular MOAs underlying hepatotoxicity have been described: (i) DNA damage response, (ii) induction of apoptosis by cellular oxidative stress through the intrinsic mitochondrial pathway or excessive production of reactive oxygen species (ROS), and (iii) effects on bile acid homeostasis (Gluck et al., 2020; Wang et al., 2021). Among the pathways that were impaired by LAS, the metabolism of bile acid has been extensively studied in the literature and shown to be greatly affected by different PAs [24,25,27, 59]. In particular, Waizenegger et al. revealed that the hepatotoxicity of four PAs is associated with a significant decrease in bile acids in both the intracellular and extracellular contents of HepaRG cells [27]. In our study, the responses of fatty acids and bile acids followed a bell-shaped trend, indicating an impact on bile acid homeostasis induced by both low and high concentrations of LAS. Indeed, in response to low concentrations of hepatotoxicants, cells may compensate for the mechanistic effects by increasing their lipid content, but this compensation may not occur at higher concentrations of the toxicants because lipid synthesis is also affected [20]. The amino acid metabolism was impaired at higher BMCs than the lipid metabolism, but still at very low concentrations. Like lipids and bile acids, most of the amino acids responded with a bell-shaped trend throughout the exposure range. In support of this observation, primary bile acids are very rapidly conjugated with glycine and taurine after their synthesis [59]. The primary purpose of bile acids is to facilitate digestion of fat through its surfactant properties, with emulsification of fats into micelles [60]. Amino acids, lipids and the bile acid glycocholic acid have also been described as being impaired in a recent in vivo study in the plasma and urine of rats exposed to a decoction of *Gynura segetum*, an herbal medicine containing PAs [26]. The authors identified biomarkers of hepatotoxicity, including some amino acids and lipids, using a metabolic profiling analysis [26]. In HepaRG cells, we also observed impairment of metabolites conjugated with glutamic acid, mainly with a linear increasing trend. This finding can be associated with disruptions in the tricarboxylic acid (TCA) cycle [47], potentially causing oxidative stress, and is further? corroborated by the impaired level of cis- or trans-aconitic acid. Among other impaired amino acids, nicotinuric acid showed a linear decreasing trend, partly explaining the increased production of amino acids involved in glutamate metabolism. In addition, impairment of the mitochondrial function can be revealed by disorders observed in lipids, organic acids and nucleotides, such as modifications in the metabolism of sphingolipids describing a bell-shaped trend, and an effect on UTP production at a very low BMCs. Finally, the steroid class including steroids other than bile acids, such as hydroxysteroids, pregnane, androstane and estrane steroids, were impaired at higher concentrations than lipids, and followed a U-shaped trend. To our knowledge, this is the first evidence that PA may affect steroid metabolism. Their concentration-response trend was contrary to that of bile acids and lipids. Therefore, LAS may act as a

pregnane X receptor (PXR) agonist or antagonist, because this transcription factor is known to act as a physiological sensor in the regulation of bile acid and lipid metabolisms and of steroid and endocrine homeostasis [61]. A previous study showed that LAS displayed PXR-activating properties in PXR-transfected HepG2 cells, and this activation may participate in liver steatosis ([62]; https://aopwiki.org/aops/60).

### 4.4. BMC assessment through metabolomics

The assessment of PODs such as BMCs, derived from concentration-response relationship of in vitro omics data is fully compliant with NAMs. However, there is a need for further knowledge and consensus among the scientific community before it can be useful and fully accepted for use in risk assessment. Genomic and transcriptomic data have been explored for this purpose, and have led to guidelines and recommendations [33,49,63], but application to metabolomics requires further investigation. In this study, we explored the concentration-response trends and derived BMDs for intra- and extracellular metabolites of HepaRG cells produced after treatment with three PAs. The time of exposure (48 h) was deemed appropriate for obtaining a sufficient number of affected metabolites. A previous comparative study indicated that the time points of 24 h and 48 h produced consistent BMC values for three out of the four tested chemicals [37]. We were also able to compare trends in the modeled responses of unannotated metabolites as well as the BMCs from both the exo- and endometabolome. We concluded that, for unannotated metabolites, EXO provided a less sensitive BMC Q1 than ENDO. This result may depend on the class of compounds and their toxicokinetic profile. Moreover, the comparison of different BMCs for LAS shows that lipid metabolism is by far the most sensitive class of metabolites, even compared with BMCs derived for genotoxic potency (Table 2). Nevertheless, BMCs derived from specific pathway or subpathway classes must be supported by a sufficient level of confidence in metabolite annotation, particularly when very few metabolites are annotated in a given (sub)pathway. The metabolites that could be annotated mostly showed L3 confidence levels, with the rest being L2 (mzCloud matching). Confidence levels are also crucial when defining the metabolites that can be used as input for regulatory toxicology [64]. Finally, as for transcriptomics, a framework needs to be established to define, in view of the thousands of metabolites involved, which BMCs must be examined (1st percentile, 25th percentile (Q1), median) so as to be sufficiently protective and comparable to traditional apical PODs [34,49,55].

## 5. Conclusion

Overall, our findings indicate that exposure to PAs via food and beverages may have an impact on health, with hepatotoxic effects observed in HepaRG cells at very low concentrations of PAs in our investigation of impaired metabolism pathways. This study applied concentration-response modeling and BMC derivation to understand and dissociate the hepatotoxic effects of three different PAs. To our knowledge, this is the first study showing the benefit of such an approach on several contaminants from the same family, but with different chemical structures. Using unannotated metabolites revealed that the endometabolome was more sensitive than the exometabolome regarding effects on the metabolome, and that LAS was confirmed as the most hepatotoxic PA. Interestingly, our trend analysis highlighted the relevance of bell-shaped responses as indicators of early adaptive changes. At a second-tier level, the MOA of LAS was enriched by highlighting specific perturbed biochemical pathways such as lipid metabolism, which was the most sensitive impaired pathway, followed by amino acid metabolism. In particular, bile acid deregulation was confirmed, which is consistent with several previous studies; however, our is the first report of impairment of steroid metabolism, necessitating further studies. The response of HepaRG cells to PA exposure led to
substantial changes in metabolites involved in toxicity, as well as those with a role in adaptation mechanisms. In particular, fatty acids and bile acids are implicated in the most sensitive pathways that were impaired by PAs, suggesting that they may be biomarkers of hepatotoxic response.

Improving the predictivity of this original method now requires testing new advanced liver models such as 3D spheroids of HepaRG cells or primary human hepatocytes and co-cultures with non-parenchymal cells to translate in vitro PODs into equivalent exposures for humans. In conclusion, as already shown in previous studies, our work confirms that the combination of metabolomics with a concentration-response tool to derive PODs is a promising approach for NGRA, and highlights the urgent need for recommendations on deriving in vitro PODs from metabolomic experiments.

### CRediT authorship contribution statement

**Estelle Dubreil:** Writing – review & editing, Writing – original draft, Supervision, Resources, Project administration, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Keyvin Darney:** Writing – review & editing, Validation, Formal analysis. **Marie-Laure Delignette-Muller:** Writing – review & editing, Software, Formal analysis. **Barranger Audrey:** Writing – review & editing, Investigation, Formal analysis. **Sylvie Huet:** Resources, Investigation, **Kevin Hogeveen:** Writing – review & editing, Validation, Formal analysis. **Thibaut Léger:** Writing – review & editing, Methodology, Investigation, Conceptualization. **Valérie Fessard:** Writing – review & editing, Supervision, Resources, Project administration. **Ludovic Le Hégarat:** Writing – review & editing, Supervision, Resources, Project administration, Formal analysis, Conceptualization.

### 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.

### Appendix A. Supporting information

Supplementary data associated with this article can be found in the online version at doi:10.1016/j.jhazmat.2024.134721.

### References

[1] Akinboye, A.J., Kim, K., Choi, S., Yang, I., Lee, J.G., 2023. Alkaloids in food: a review of toxicity, analytical methods, occurrence and risk assessments. Food Sci Biotechnol.

[2] Arseculeratne, S.N., Gunatilaka, A.A.L., Panabokke, R.G., 1981. Studies on medicinal plants of Sri Lanka: occurrence of pyrrolizidine alkaloids and hepatotoxic properties in some traditional medicinal herbs. J Ethnopharmacol 4, 159-177.

[3] Deinzer, M.L., Thomson, P.A., Burgett, D.M., Isaacson, D.L., 1977. Pyrrolizidine alkaloids: their occurrence in honey from tansy ragwort (Senecio jacobaea L.). Science 195, 497-499.

[4] Mattocks, A.R., 1971. The occurrence and analysis of pyrrolizidine alkaloid n-oxides. Xenobiotica 1, 451-453.

[5] Parton, K., Bruere, A.N., 2002. Plant poisoning of livestock in New Zealand. NZ Vet J 50, 22-27.

[6] Stewart, M.J., Steenkamp, V., 2001. Pyrrolizidine poiso²ning: a neglected area in human toxicology. Ther Drug Monit 23, 698-708.

[7] Tandon, B.N., Tandon, H.D., Tandon, R.K., Narndranathan, M., Joshi, Y.K., 1976. An epidemic of veno-occlusive disease of liver in central India. Lancet 308, 271-272.

[8] Chain, E.Po.CitF., 2011. Scientific opinion on pyrrolizidine alkaloids in food and feed. EFSA J 9, 2406.

[9] Chain, E Panel oCitF., Knutsen, H.K., Alexander, J., Barregård, L., Bignami, M., Brüscheiler, B., et al., 2017. Risks for human health related to the presence of pyrrolizidine alkaloids in honey, tea, herbal infusions and food supplements. EFSA J 15, e04908.

[10] Schrenk, D., Gao, L., Lin, G., Mahony, C., Mulder, P.P.J., Peijnenburg, A., et al., 2020. Pyrrolizidine alkaloids in food and phytomedicine: Occurrence, exposure, toxicity, mechanisms, and risk assessment - A review. Food Chem Toxicol 136.

[11] Buchmueller, J., Kaltner, F., Gottschalk, C., Maares, M., Braeuning, A., Hessel-Pras, S., 2022. Structure-dependent toxicokinetics of selected pyrrolizidine alkaloids in vitro. Int J Mol Sci 23.

[12] Mulder, P.P.J., López, P., Castelari, M., Bodi, D., Ronczka, S., Preiss-Weigert, A., et al., 2018. Occurrence of pyrrolizidine alkaloids in animal- and plant-derived food: results of a survey across Europe. Food Addit Contam - Part A Chem, Anal, Control, Expo Risk Assess 35, 118-133.

[13] Roncada, P., Isani, G., Peloso, M., Dalmonte, T., Bonan, S., Caprai, E., 2023. Pyrrolizidine alkaloids from monofloral and multifloral Italian honey. Int J Environ Res Public Health 20.

[14] Haas, M., Wirachowski, K., Thibol, L., Küpper, J.H., Schrenk, D., Fahrer, J., 2023. Potency ranking of pyrrolizidine alkaloids in metabolically competent human liver cancer cells and primary human hepatocytes using a genotoxicity test battery. Arch Toxicol 97, 1413-1428.

[15] He, Y., Zhu, L., Ma, J., Lin, G., 2021. Metabolism-mediated cytotoxicity and genotoxicity of pyrrolizidine alkaloids. Arch Toxicol 95, 1917-1942.

[16] Ma, J., Li, M., Li, N., Chan, W.Y., Lin, G., 2021. Pyrrolizidine alkaloid-induced hepatotoxicity associated with the formation of reactive metabolite-derived pyrrole-protein adducts. Toxins 13.

[17] Chen, X., Ma, J., He, Y., Xue, J., Song, Z., Xu, Q., et al., 2021. Characterization of liver injury induced by a pyrrolizidine alkaloid in rats. Phytomedicine 89.

[18] Ebmeyer, J., Rasinger, J.D., Hengstler, J.G., Schaudien, D., Creutzenberg, O., Lampen, A., et al., 2020. Hepatotoxic pyrrolizidine alkaloids induce DNA damage response in rat liver in a 28-day feeding study. Arch Toxicol 94, 1739-1751.

[19] Rivera-Pérez, A., Romero-González, R., Garrido Frenich, A., 2023. Determination and occurrence of alkenylbenzenes, pyrrolizidine and tropane alkaloids in spices, herbs, teas, and other plant-derived food products using chromatographic methods: review from 2010-2020. Food Rev Int 39, 1110-1136.

[20] Cuykx, M., Beirnaert, C., Rodrigues, R.M., Laukens, K., Vanhaecke, T., Covaci, A., 2019. Untargeted liquid chromatography-mass spectrometry metabolomics to assess drug-induced cholestatic features in HepaRG® cells. Toxicol Appl Pharmacol 379.

[21] Iturrope, E., Da Silva, K.M., Robeyns, R., Van De Lavoir, M., Boeckmans, J., Vanhaecke, T., et al., 2022. Metabolic signature of ethanol-induced hepatotoxicity in HepaRG cells by liquid chromatography-mass spectrometry-based untargeted metabolomics. J Proteome Res 21, 1153-1166.

[22] Martínez-Sena, T., Moro, E., Moreno-Torres, M., Quintás, G., Hengstler, J., Castell, J.V., 2023. Metabolomics-based strategy to assess drug hepatotoxicity and uncover the mechanisms of hepatotoxicity involved. Arch Toxicol.

[23] Rodrigues, R.M., Kollipara, L., Chaudhari, U., Sachinidis, A., Zahedi, R.P., Sickmann, A., et al., 2018. Omics-based responses induced by bosentan in human hepatoma HepaRG cell cultures. Arch Toxicol 92, 1939-1952.

[24] Xiong, A., Yang, F., Fang, L., Yang, L., He, Y., Wan, Y.J.Y., et al., 2014. Erratum: metabolomic and genomic evidence for compromised bile acid homeostasis by senecionine, a hepatotoxic pyrrolizidine alkaloid. Chem Res Toxicol 27 (5) (775-786) DOI: 10.1021/tx400451q), Chemical Research in Toxicology, 27 (2014) 1458.

[25] Xiong, A., Yang, F., Fang, L., Yang, L., He, Y., Wan, Y.Y.J., et al., 2014. Metabolomic and genomic evidence for compromised bile acid homeostasis by senecionine, a hepatotoxic pyrrolizidine alkaloid. Chem Res Toxicol 27, 775-786.

[26] Qiu, S., Zhang, H., Fei, Q., Zhu, F., Wang, J., Jia, X., et al., 2018. Urine and plasma metabolomics study on potential hepatotoxic biomarkers identification in rats induced by Gynura segetum. J Ethnopharmacol 216, 37-46.

[27] Waizenegger, J., Glück, J., Henricsson, M., Luckert, C., Braeuning, A., Hessel-Pras, S., 2021. Pyrrolizidine alkaloids disturb bile acid homeostasis in the human hepatoma cell line heparg. Foods 10.

[28] Olesti, E., González-Ruiz, V., Wilks, M.F., Boccard, J., Rudaz, S., 2021. Approaches in metabolomics for regulatory toxicology applications. Analyst 146, 1820-1834.

[29] Vinken, M., 2019. Omics-based input and output in the development and use of adverse outcome pathways. Curr Opin Toxicol 18, 8-12.

[30] Black, M.B., Stern, A., Efremenko, A., Mallick, P., Moreau, M., Hartman, J.K., et al., 2022. Biological system considerations for application of toxicogenomics in next-generation risk assessment and predictive toxicology. Toxicol Vitr 80.

[31] Viant, M.R., Ebbels, T.M.D., Beger, R.D., Ekman, D.R., Epps, D.J.T., Kamp, H., et al., 2019. Use cases, best practice and reporting standards for metabolomics in regulatory toxicology. Nat Commun 10.

[32] Committee, E.S., Hardy, A., Benford, D., Halldorsson, T., Jeger, M.J., Knutsen, K. H., et al., 2017. Update: use of the benchmark dose approach in risk assessment. EFSA J 15, e04658.

[33] Committee, E.S., More, S.J., Bampidis, V., Benford, D., Bragard, C., Halldorsson, T. L., et al., 2022. Guidance on the use of the benchmark dose approach in risk assessment. EFSA J 20, e07584.

[34] Harrill, J.A., Viant, M.R., Yauk, C.L., Sachana, M., Gant, T.W., Auerbach, S.S., et al., 2021. Progress towards an OECD reporting framework for transcriptomics and metabolomics in regulatory toxicology. Regul Toxicol Pharmacol 125.

[35] Gwinn, W.M., Auerbach, S.S., Parham, F., Stout, M.D., Waidyanatha, S., Mutlu, E., et al., 2020. Evaluation of 5-day in vivo rat liver and kidney with high-throughput transcriptomics for estimating benchmark doses of apical outcomes. Toxicol Sci 176, 343-354.

[36] Crizer, D.M., Ramaiahgari, S.C., Ferguson, S.S., Rice, J.R., Dunlap, P.E., Sipes, N.S., et al., 2021. Benchmark concentrations for untargeted metabolomics versus transcriptomics for liver injury compounds in in vitro liver models. Toxicol Sci 181, 175-186.
substantial changes in metabolites involved in toxicity, as well as those with a role in adaptation mechanisms. In particular, fatty acids and bile acids are implicated in the most sensitive pathways that were impaired by PAs, suggesting that they may be biomarkers of hepatotoxic response.

Improving the predictivity of this original method now requires testing new advanced liver models such as 3D spheroids of HepaRG cells or primary human hepatocytes and co-cultures with non-parenchymal cells to translate in vitro PODs into equivalent exposures for humans. In conclusion, as already shown in previous studies, our work confirms that the combination of metabolomics with a concentration-response tool to derive PODs is a promising approach for NGRA, and highlights the urgent need for recommendations on deriving in vitro PODs from metabolomic experiments.

### CRediT authorship contribution statement

**Estelle Dubreil:** Writing – review & editing, Writing – original draft, Supervision, Resources, Project administration, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Keyvin Darney:** Writing – review & editing, Validation, Formal analysis. **Marie-Laure Delignette-Muller:** Writing – review & editing, Software, Formal analysis. **Barranger Audrey:** Writing – review & editing, Investigation, Formal analysis. **Sylvie Huet:** Resources, Investigation, **Kevin Hogeveen:** Writing – review & editing, Validation, Formal analysis. **Thibaut Léger:** Writing – review & editing, Methodology, Investigation, Conceptualization. **Valérie Fessard:** Writing – review & editing, Supervision, Resources, Project administration. **Ludovic Le Hégarat:** Writing – review & editing, Supervision, Resources, Project administration, Formal analysis, Conceptualization.

### 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.

### Appendix A. Supporting information

Supplementary data associated with this article can be found in the online version at doi:10.1016/j.jhazmat.2024.134721.

### References

[1] Akinboye, A.J., Kim, K., Choi, S., Yang, I., Lee, J.G., 2023. Alkaloids in food: a review of toxicity, analytical methods, occurrence and risk assessments. Food Sci Biotechnol.

[2] Arseculeratne, S.N., Gunatilaka, A.A.L., Panabokke, R.G., 1981. Studies on medicinal plants of Sri Lanka: occurrence of pyrrolizidine alkaloids and hepatotoxic properties in some traditional medicinal herbs. J Ethnopharmacol 4, 159-177.

[3] Deinzer, M.L., Thomson, P.A., Burgett, D.M., Isaacson, D.L., 1977. Pyrrolizidine alkaloids: their occurrence in honey from tansy ragwort (Senecio jacobaea L.). Science 195, 497-499.

[4] Mattocks, A.R., 1971. The occurrence and analysis of pyrrolizidine alkaloid n-oxides. Xenobiotica 1, 451-453.

[5] Parton, K., Bruere, A.N., 2002. Plant poisoning of livestock in New Zealand. NZ Vet J 50, 22-27.

[6] Stewart, M.J., Steenkamp, V., 2001. Pyrrolizidine poiso²ning: a neglected area in human toxicology. Ther Drug Monit 23, 698-708.

[7] Tandon, B.N., Tandon, H.D., Tandon, R.K., Narndranathan, M., Joshi, Y.K., 1976. An epidemic of veno-occlusive disease of liver in central India. Lancet 308, 271-272.

[8] Chain, E.Po.CitF., 2011. Scientific opinion on pyrrolizidine alkaloids in food and feed. EFSA J 9, 2406.

[9] Chain, E Panel oCitF., Knutsen, H.K., Alexander, J., Barregård, L., Bignami, M., Brüscheiler, B., et al., 2017. Risks for human health related to the presence of pyrrolizidine alkaloids in honey, tea, herbal infusions and food supplements. EFSA J 15, e04908.

[10] Schrenk, D., Gao, L., Lin, G., Mahony, C., Mulder, P.P.J., Peijnenburg, A., et al., 2020. Pyrrolizidine alkaloids in food and phytomedicine: Occurrence, exposure, toxicity, mechanisms, and risk assessment - A review. Food Chem Toxicol 136.

[11] Buchmueller, J., Kaltner, F., Gottschalk, C., Maares, M., Braeuning, A., Hessel-Pras, S., 2022. Structure-dependent toxicokinetics of selected pyrrolizidine alkaloids in vitro. Int J Mol Sci 23.

[12] Mulder, P.P.J., López, P., Castelari, M., Bodi, D., Ronczka, S., Preiss-Weigert, A., et al., 2018. Occurrence of pyrrolizidine alkaloids in animal- and plant-derived food: results of a survey across Europe. Food Addit Contam - Part A Chem, Anal, Control, Expo Risk Assess 35, 118-133.

[13] Roncada, P., Isani, G., Peloso, M., Dalmonte, T., Bonan, S., Caprai, E., 2023. Pyrrolizidine alkaloids from monofloral and multifloral Italian honey. Int J Environ Res Public Health 20.

[14] Haas, M., Wirachowski, K., Thibol, L., Küpper, J.H., Schrenk, D., Fahrer, J., 2023. Potency ranking of pyrrolizidine alkaloids in metabolically competent human liver cancer cells and primary human hepatocytes using a genotoxicity test battery. Arch Toxicol 97, 1413-1428.

[15] He, Y., Zhu, L., Ma, J., Lin, G., 2021. Metabolism-mediated cytotoxicity and genotoxicity of pyrrolizidine alkaloids. Arch Toxicol 95, 1917-1942.

[16] Ma, J., Li, M., Li, N., Chan, W.Y., Lin, G., 2021. Pyrrolizidine alkaloid-induced hepatotoxicity associated with the formation of reactive metabolite-derived pyrrole-protein adducts. Toxins 13.

[17] Chen, X., Ma, J., He, Y., Xue, J., Song, Z., Xu, Q., et al., 2021. Characterization of liver injury induced by a pyrrolizidine alkaloid in rats. Phytomedicine 89.

[18] Ebmeyer, J., Rasinger, J.D., Hengstler, J.G., Schaudien, D., Creutzenberg, O., Lampen, A., et al., 2020. Hepatotoxic pyrrolizidine alkaloids induce DNA damage response in rat liver in a 28-day feeding study. Arch Toxicol 94, 1739-1751.

[19] Rivera-Pérez, A., Romero-González, R., Garrido Frenich, A., 2023. Determination and occurrence of alkenylbenzenes, pyrrolizidine and tropane alkaloids in spices, herbs, teas, and other plant-derived food products using chromatographic methods: review from 2010-2020. Food Rev Int 39, 1110-1136.

[20] Cuykx, M., Beirnaert, C., Rodrigues, R.M., Laukens, K., Vanhaecke, T., Covaci, A., 2019. Untargeted liquid chromatography-mass spectrometry metabolomics to assess drug-induced cholestatic features in HepaRG® cells. Toxicol Appl Pharmacol 379.

[21] Iturrope, E., Da Silva, K.M., Robeyns, R., Van De Lavoir, M., Boeckmans, J., Vanhaecke, T., et al., 2022. Metabolic signature of ethanol-induced hepatotoxicity in HepaRG cells by liquid chromatography-mass spectrometry-based untargeted metabolomics. J Proteome Res 21, 1153-1166.

[22] Martínez-Sena, T., Moro, E., Moreno-Torres, M., Quintás, G., Hengstler, J., Castell, J.V., 2023. Metabolomics-based strategy to assess drug hepatotoxicity and uncover the mechanisms of hepatotoxicity involved. Arch Toxicol.

[23] Rodrigues, R.M., Kollipara, L., Chaudhari, U., Sachinidis, A., Zahedi, R.P., Sickmann, A., et al., 2018. Omics-based responses induced by bosentan in human hepatoma HepaRG cell cultures. Arch Toxicol 92, 1939-1952.

[24] Xiong, A., Yang, F., Fang, L., Yang, L., He, Y., Wan, Y.J.Y., et al., 2014. Erratum: metabolomic and genomic evidence for compromised bile acid homeostasis by senecionine, a hepatotoxic pyrrolizidine alkaloid. Chem Res Toxicol 27 (5) (775-786) DOI: 10.1021/tx400451q), Chemical Research in Toxicology, 27 (2014) 1458.

[25] Xiong, A., Yang, F., Fang, L., Yang, L., He, Y., Wan, Y.Y.J., et al., 2014. Metabolomic and genomic evidence for compromised bile acid homeostasis by senecionine, a hepatotoxic pyrrolizidine alkaloid. Chem Res Toxicol 27, 775-786.

[26] Qiu, S., Zhang, H., Fei, Q., Zhu, F., Wang, J., Jia, X., et al., 2018. Urine and plasma metabolomics study on potential hepatotoxic biomarkers identification in rats induced by Gynura segetum. J Ethnopharmacol 216, 37-46.

[27] Waizenegger, J., Glück, J., Henricsson, M., Luckert, C., Braeuning, A., Hessel-Pras, S., 2021. Pyrrolizidine alkaloids disturb bile acid homeostasis in the human hepatoma cell line heparg. Foods 10.

[28] Olesti, E., González-Ruiz, V., Wilks, M.F., Boccard, J., Rudaz, S., 2021. Approaches in metabolomics for regulatory toxicology applications. Analyst 146, 1820-1834.

[29] Vinken, M., 2019. Omics-based input and output in the development and use of adverse outcome pathways. Curr Opin Toxicol 18, 8-12.

[30] Black, M.B., Stern, A., Efremenko, A., Mallick, P., Moreau, M., Hartman, J.K., et al., 2022. Biological system considerations for application of toxicogenomics in next-generation risk assessment and predictive toxicology. Toxicol Vitr 80.

[31] Viant, M.R., Ebbels, T.M.D., Beger, R.D., Ekman, D.R., Epps, D.J.T., Kamp, H., et al., 2019. Use cases, best practice and reporting standards for metabolomics in regulatory toxicology. Nat Commun 10.

[32] Committee, E.S., Hardy, A., Benford, D., Halldorsson, T., Jeger, M.J., Knutsen, K. H., et al., 2017. Update: use of the benchmark dose approach in risk assessment. EFSA J 15, e04658.

[33] Committee, E.S., More, S.J., Bampidis, V., Benford, D., Bragard, C., Halldorsson, T. L., et al., 2022. Guidance on the use of the benchmark dose approach in risk assessment. EFSA J 20, e07584.

[34] Harrill, J.A., Viant, M.R., Yauk, C.L., Sachana, M., Gant, T.W., Auerbach, S.S., et al., 2021. Progress towards an OECD reporting framework for transcriptomics and metabolomics in regulatory toxicology. Regul Toxicol Pharmacol 125.

[35] Gwinn, W.M., Auerbach, S.S., Parham, F., Stout, M.D., Waidyanatha, S., Mutlu, E., et al., 2020. Evaluation of 5-day in vivo rat liver and kidney with high-throughput transcriptomics for estimating benchmark doses of apical outcomes. Toxicol Sci 176, 343-354.

[36] Crizer, D.M., Ramaiahgari, S.C., Ferguson, S.S., Rice, J.R., Dunlap, P.E., Sipes, N.S., et al., 2021. Benchmark concentrations for untargeted metabolomics versus transcriptomics for liver injury compounds in in vitro liver models. Toxicol Sci 181, 175-186.
[37] Malinowska, J.M., Palosaari, T., Sund, J., Carpi, D., Weber, R.J.M., Lloyd, G.R., et al., 2023. Derivation of metabolic point of departure using high-throughput in vitro metabolomics: investigating the importance of sampling time points on benchmark concentration values in the HepaRG cell line. Arch Toxicol 97, 721-735.

[38] Creusot, N., Chaumet, B., Eon, M., Mazzella, N., Moreira, A., Morin, S., 2022. Metabolomics insight into the influence of environmental factors in responses of freshwater biofilms to the model herbicide diuron. Environ Sci Pollut Res 29, 29332-29347.

[39] Larras, F., Billoir, E., Scholz, S., Tarkka, M., Wubet, T., Delignette-Muller, M.L., et al., 2020. A multi-omics concentration-response framework uncovers novel understanding of triclosan effects in the chlorophyte Scenedesmus vacuolatus. J Hazard Mater 397.

[40] Larras, F., Billoir, E., Baillard, V., Siberchicot, A., Scholz, S., Wubet, T., et al., 2018. DRomics: a Turnkey TOOL to Support the Use of the Dose-response Framework for Omics Data in Ecological Risk Assessment. Environ Sci Technol 52, 14461-14468.

[41] Araújo, A.M., Carvalho, F., De Pinho, P.G., Carvalho, M., 2021. Toxicometabolomics: small molecules to answer big toxicological questions. Metabolites 11.

[42] Balerczyk, A., Damblon, C., Elena-Herrmann, B., Panthu, B., Ratureau, G.J.P., 2020. Metabolomic approaches to study chemical exposure-related metabolism alterations in mammalian cell cultures. Int J Mol Sci 21, 1-22.

[43] D.-M. Marie Laure, S. Aurélie, L. Floriane, B. Elise, DRomics, a workflow to model and make sense of dose-response (multi-)omics data in (eco)toxicology, bioRxiv, (2023) 2023.02.2002.527852.

[44] Léger, T., Balaguer, P., Le Hégarat, L., Fessard, V., 2023. Fate and PPARγ and STATs-driven effects of the mitochondrial complex I inhibitor tebufenpyrad in liver cells revealed with multi-omics. J Hazard Mater 442.

[45] Delignette-Muller, M.L., Siberchicot, A., Larras, F., Billoir, E., 2023. DRomics, a workflow to exploit dose-response omics data in ecotoxicology. Peer Community J 3.

[46] Sostare, E., Lawson, T.N., Saunders, L.R., Colbourne, J.K., Weber, R.J.M., Sobanski, T., et al., 2022. Knowledge-driven approaches to create the MTox700+ Metabolite panel for predicting toxicity. Toxicol Sci 186, 208-220.

[47] Cuykx, M., Rodrigues, R.M., Laukens, K., Vanhaecke, T., Covaci, A., 2018. In vitro assessment of hepatotoxicity by metabolomics: a review. Arch Toxicol 92, 3007-3029.

[48] Schymanski, E.L., Jeon, J., Gulde, R., Fenner, K., Ruff, M., Singer, H.P., et al., 2014. Identifying small molecules via high-resolution mass spectrometry: Communicating confidence. Environ Sci Technol 48, 2097-2098.

[49] Farmahin, R., Williams, A., Kuo, B., Chepelev, N.L., Thomas, R.S., Barton-Maclaren, T.S., et al., 2017. Recommended approaches in the application of toxicogenomics to derive points of departure for chemical risk assessment. Arch Toxicol 91, 2045-2065.

[50] Anthérieu, S., Chesné, C., Li, R., Guguen-Guillouzo, C., Guillouzo, A., 2012. Optimization of the HepaRG cell model for drug metabolism and toxicity studies. Toxicol Vitr 26, 1278-1285.

[51] Jiang, J., Pieterman, C.D., Ertaylan, G., Peeters, R.L.M., de Kok, T.M.C.M., 2019. The application of omics-based human liver platforms for investigating the mechanism of drug-induced hepatotoxicity in vitro. Arch Toxicol 93, 3067-3098.

[52] Ebmeyer, J., Braeuning, A., Glatt, H., These, A., Hessel-Pras, S., Lampen, A., 2019. Human CYP3A4-mediated toxification of the pyrrolizidine alkaloid lasiocarpine. Food Chem Toxicol 130, 79-88.

[53] Colas, S., Le Faucheur, S., 2024. How do biomarkers dance? Specific moves of defense and damage biomarkers for biological interpretation of dose-response model trends. J Hazard Mater 465.

[54] Ramaiahgari, S.C., Auerbach, S.S., Saddler, T.O., Rice, J.R., Dunlap, P.E., Sipes, N. S., et al., 2019. The power of resolution: Contextualized understanding of biological responses to liver injury chemicals using high-throughput transcriptomics and benchmark concentration modeling. Toxicol Sci 169, 553-566.

[55] Reardon, A.J.F., Rowan-Carroll, A., Ferguson, S.S., Leingartner, K., Gagne, R., Kuo, B., et al., 2021. Potency ranking of per- and polyfluoroalkyl substances using high-throughput transcriptomic analysis of human liver spheroids. Toxicol Sci 184, 154-169.

[56] Waizenegger, J., Braeuning, A., Templin, M., Lampen, A., Hessel-Pras, S., 2018. Structure-dependent induction of apoptosis by hepatotoxic pyrrolizidine alkaloids in the human hepatoma cell line HepaRG: Single versus repeated exposure. Food Chem Toxicol 114, 215-226.

[57] Louisse, J., Rijkers, D., Stoopen, G., Holleboom, W.J., Delagrange, M., Molthof, E., et al., 2019. Determination of genotoxic potencies of pyrrolizidine alkaloids in HepaRG cells using the γH2AX assay. Food Chem Toxicol 131, 110532.

[58] Rutz, L., Gao, L., Küpper, J.H., Schrenk, D., 2020. Structure-dependent genotoxic potencies of selected pyrrolizidine alkaloids in metabolically competent HepG2 cells. Arch Toxicol 94, 4159-4172.

[59] Glück, J., Waizenegger, J., Braeuning, A., Hessel-Pras, S., 2021. Pyrrolizidine alkaloids induce cell death in human hepatic cells in a structure-dependent manner. Int J Mol Sci 22, 1-12.

[60] Di Ciaula, A., Garruti, G., Baccetto, R.L., Molina-Molina, E., Bonfrate, L., Wang, D. Q.H., et al., 2017. Bile acid physiology. Ann Hepatol 16, s4-s14.

[61] Liang, Y., Gong, Y., Jiang, Q., Yu, Y., Zhang, J., 2023. Environmental endocrine disruptors and pregnane X receptor action: a review. Food Chem Toxicol 179, 113976.

[62] Luckert, C., Braeuning, A., Lampen, A., Hessel-Pras, S., 2018. PXR: structure-specific activation by hepatotoxic pyrrolizidine alkaloids. Chem-Biol Interact 288, 38-48.

[63] P. National Toxicology, NTP Research Reports, in: NTP Research Report on National Toxicology Program Approach to Genomic Dose-Response Modeling: Research Report 5, National Toxicology Program, Durham (NC), 2018.

[64] Malinowska, J.M., Viant, M.R., 2019. Confidence in metabolite identification dictates the applicability of metabolomics to regulatory toxicology. Curr Opin Toxicol 16, 32-38.

[65] Allemang, A., Mahony, C., Lester, C., Pfuhler, S., 2018. Relative potency of fifteen pyrrolizidine alkaloids to induce DNA damage as measured by micronucleus induction in HepaRG human liver cells. Food Chem Toxicol 121, 72-81.
[37] Malinowska, J.M., Palosaari, T., Sund, J., Carpi, D., Weber, R.J.M., Lloyd, G.R., et al., 2023. Derivation of metabolic point of departure using high-throughput in vitro metabolomics: investigating the importance of sampling time points on benchmark concentration values in the HepaRG cell line. Arch Toxicol 97, 721-735.

[38] Creusot, N., Chaumet, B., Eon, M., Mazzella, N., Moreira, A., Morin, S., 2022. Metabolomics insight into the influence of environmental factors in responses of freshwater biofilms to the model herbicide diuron. Environ Sci Pollut Res 29, 29332-29347.

[39] Larras, F., Billoir, E., Scholz, S., Tarkka, M., Wubet, T., Delignette-Muller, M.L., et al., 2020. A multi-omics concentration-response framework uncovers novel understanding of triclosan effects in the chlorophyte Scenedesmus vacuolatus. J Hazard Mater 397.

[40] Larras, F., Billoir, E., Baillard, V., Siberchicot, A., Scholz, S., Wubet, T., et al., 2018. DRomics: a Turnkey TOOL to Support the Use of the Dose-response Framework for Omics Data in Ecological Risk Assessment. Environ Sci Technol 52, 14461-14468.

[41] Araújo, A.M., Carvalho, F., De Pinho, P.G., Carvalho, M., 2021. Toxicometabolomics: small molecules to answer big toxicological questions. Metabolites 11.

[42] Balerczyk, A., Damblon, C., Elena-Herrmann, B., Panthu, B., Ratureau, G.J.P., 2020. Metabolomic approaches to study chemical exposure-related metabolism alterations in mammalian cell cultures. Int J Mol Sci 21, 1-22.

[43] D.-M. Marie Laure, S. Aurélie, L. Floriane, B. Elise, DRomics, a workflow to model and make sense of dose-response (multi-)omics data in (eco)toxicology, bioRxiv, (2023) 2023.02.2002.527852.

[44] Léger, T., Balaguer, P., Le Hégarat, L., Fessard, V., 2023. Fate and PPARγ and STATs-driven effects of the mitochondrial complex I inhibitor tebufenpyrad in liver cells revealed with multi-omics. J Hazard Mater 442.

[45] Delignette-Muller, M.L., Siberchicot, A., Larras, F., Billoir, E., 2023. DRomics, a workflow to exploit dose-response omics data in ecotoxicology. Peer Community J 3.

[46] Sostare, E., Lawson, T.N., Saunders, L.R., Colbourne, J.K., Weber, R.J.M., Sobanski, T., et al., 2022. Knowledge-driven approaches to create the MTox700+ Metabolite panel for predicting toxicity. Toxicol Sci 186, 208-220.

[47] Cuykx, M., Rodrigues, R.M., Laukens, K., Vanhaecke, T., Covaci, A., 2018. In vitro assessment of hepatotoxicity by metabolomics: a review. Arch Toxicol 92, 3007-3029.

[48] Schymanski, E.L., Jeon, J., Gulde, R., Fenner, K., Ruff, M., Singer, H.P., et al., 2014. Identifying small molecules via high-resolution mass spectrometry: Communicating confidence. Environ Sci Technol 48, 2097-2098.

[49] Farmahin, R., Williams, A., Kuo, B., Chepelev, N.L., Thomas, R.S., Barton-Maclaren, T.S., et al., 2017. Recommended approaches in the application of toxicogenomics to derive points of departure for chemical risk assessment. Arch Toxicol 91, 2045-2065.

[50] Anthérieu, S., Chesné, C., Li, R., Guguen-Guillouzo, C., Guillouzo, A., 2012. Optimization of the HepaRG cell model for drug metabolism and toxicity studies. Toxicol Vitr 26, 1278-1285.

[51] Jiang, J., Pieterman, C.D., Ertaylan, G., Peeters, R.L.M., de Kok, T.M.C.M., 2019. The application of omics-based human liver platforms for investigating the mechanism of drug-induced hepatotoxicity in vitro. Arch Toxicol 93, 3067-3098.

[52] Ebmeyer, J., Braeuning, A., Glatt, H., These, A., Hessel-Pras, S., Lampen, A., 2019. Human CYP3A4-mediated toxification of the pyrrolizidine alkaloid lasiocarpine. Food Chem Toxicol 130, 79-88.

[53] Colas, S., Le Faucheur, S., 2024. How do biomarkers dance? Specific moves of defense and damage biomarkers for biological interpretation of dose-response model trends. J Hazard Mater 465.

[54] Ramaiahgari, S.C., Auerbach, S.S., Saddler, T.O., Rice, J.R., Dunlap, P.E., Sipes, N. S., et al., 2019. The power of resolution: Contextualized understanding of biological responses to liver injury chemicals using high-throughput transcriptomics and benchmark concentration modeling. Toxicol Sci 169, 553-566.

[55] Reardon, A.J.F., Rowan-Carroll, A., Ferguson, S.S., Leingartner, K., Gagne, R., Kuo, B., et al., 2021. Potency ranking of per- and polyfluoroalkyl substances using high-throughput transcriptomic analysis of human liver spheroids. Toxicol Sci 184, 154-169.

[56] Waizenegger, J., Braeuning, A., Templin, M., Lampen, A., Hessel-Pras, S., 2018. Structure-dependent induction of apoptosis by hepatotoxic pyrrolizidine alkaloids in the human hepatoma cell line HepaRG: Single versus repeated exposure. Food Chem Toxicol 114, 215-226.

[57] Louisse, J., Rijkers, D., Stoopen, G., Holleboom, W.J., Delagrange, M., Molthof, E., et al., 2019. Determination of genotoxic potencies of pyrrolizidine alkaloids in HepaRG cells using the γH2AX assay. Food Chem Toxicol 131, 110532.

[58] Rutz, L., Gao, L., Küpper, J.H., Schrenk, D., 2020. Structure-dependent genotoxic potencies of selected pyrrolizidine alkaloids in metabolically competent HepG2 cells. Arch Toxicol 94, 4159-4172.

[59] Glück, J., Waizenegger, J., Braeuning, A., Hessel-Pras, S., 2021. Pyrrolizidine alkaloids induce cell death in human hepatic cells in a structure-dependent manner. Int J Mol Sci 22, 1-12.

[60] Di Ciaula, A., Garruti, G., Baccetto, R.L., Molina-Molina, E., Bonfrate, L., Wang, D. Q.H., et al., 2017. Bile acid physiology. Ann Hepatol 16, s4-s14.

[61] Liang, Y., Gong, Y., Jiang, Q., Yu, Y., Zhang, J., 2023. Environmental endocrine disruptors and pregnane X receptor action: a review. Food Chem Toxicol 179, 113976.

[62] Luckert, C., Braeuning, A., Lampen, A., Hessel-Pras, S., 2018. PXR: structure-specific activation by hepatotoxic pyrrolizidine alkaloids. Chem-Biol Interact 288, 38-48.

[63] P. National Toxicology, NTP Research Reports, in: NTP Research Report on National Toxicology Program Approach to Genomic Dose-Response Modeling: Research Report 5, National Toxicology Program, Durham (NC), 2018.

[64] Malinowska, J.M., Viant, M.R., 2019. Confidence in metabolite identification dictates the applicability of metabolomics to regulatory toxicology. Curr Opin Toxicol 16, 32-38.

[65] Allemang, A., Mahony, C., Lester, C., Pfuhler, S., 2018. Relative potency of fifteen pyrrolizidine alkaloids to induce DNA damage as measured by micronucleus induction in HepaRG human liver cells. Food Chem Toxicol 121, 72-81.