Pascal Petit
email: pascal.petit@univ-grenoble-alpes.fr

Dominique J Bicout

Health risk assessment with multiple reference indices 1 2

Keywords: Health risk assessment, probabilistic risk assessment, risk aggregation, cancer risk

Background: Conducting a risk assessment can be challenging, especially when dealing with several reference indices, which could lead to conflicting conclusions between studies. The common approach is to use a reference index from a single source based on the risk assessor's preference.

Objectives: To propose an approach for constructing a multi-reference index-based aggregated risk estimate using mathematical objectivity to reflect all of the available information.

Methods:

The aggregated risk estimate based on multiple reference indices (AREMRI) results from the weighted linear combination of risk distributions that were obtained with each reference index available. The weights were calculated using the degree of agreement among the reference indexbased risk distributions. The approach is illustrated through three practical cases of benzene inhalation cancer risk assessment using inhalation unit risks (IURs) from six different regulatory agencies.

Results:

The degrees of agreement between the reference index-based risk distribution, obtained with the six IURs, ranged from 0.7 to 92%. The highest weights were attributed to reference index-based risk distributions that had the highest degree of agreement with the maximum number of other reference index-based risk distributions. Regardless of the practical case considered, the AREMRI risk distribution resulted in the third highest risk compared to the six single risk distributions.

Conclusion:

Our approach can be useful in the presence of several reference indices by providing useful insights, consistency and direct comparisons between studies to support better-informed risk assessment and management decisions. This approach can shed some light on some of the uncertainties associated with the toxicological reference values in a risk assessment associated with the toxicological reference values. If the uncertainty is large, more detailed evaluation of the toxicological reference values would be needed.

Introduction

Risk assessments are paramount in ensuring the protection of human health against environmental pollutants [START_REF] Rotter | Overview on legislation and scientific approaches for risk assessment of combined exposure to multiple chemicals: the potential EuroMix contribution[END_REF]. Performing a risk assessment requires the use of toxicity reference indices. Reference indices are derived using rigorous methodological weights of evidence approaches by regulatory agencies based on available data generated in experimental toxicology and exposure assessment. For a given absorption pathway (oral, respiratory, cutaneous), dates when the index values were derived, the method of derivation, and knowledge or assumptions about substances' mechanisms of action are often different from one organization/agency/study to another (Petit et al. 2020). These differences arise mainly from the use of different data sources/key studies (e.g., the results from animal studies or epidemiological studies), extrapolation methods, correction factors, and commonly accepted assumptions and uncertainties [START_REF] Anses | Valeurs toxicologiques de référence (VTR). Guide d'élaboration de VTR. Rapport d'expertise collective[END_REF][START_REF] Beck | Approaches for describing and communicating overall uncertainty in toxicity characterizations: U.S. Environmental Protection Agency's Integrated Risk Information System (IRIS) as a case study[END_REF][START_REF] Melnick | Conflicting views on chemical carcinogenesis arising from the design and evaluation of rodent carcinogenicity studies[END_REF][START_REF] Epa | Guidance for applying quantitative data to develop data-derived extrapolation factors for interspecies and intraspecies extrapolation[END_REF]. Even when similar studies are considered for the same critical health effects and absorption pathway, organizations/agencies/studies can reach different conclusions due to different policies and/or expert judgments [START_REF] Beck | Approaches for describing and communicating overall uncertainty in toxicity characterizations: U.S. Environmental Protection Agency's Integrated Risk Information System (IRIS) as a case study[END_REF]. Depending on the reference index chosen, the results of the risk assessment can thus differ and lead to different conclusions and decisions, which may have broad and profound implications regarding the prevention of adverse effects and the protection of human health and the environment.

The double question that arises when carrying out a risk assessment is knowing which reference index to use and how to take into account the information provided by the other indices. Indeed, there are no criteria or benchmark/gold standard approaches guiding the choice of an index, especially since, by definition, all reference indices available from health/regulatory agencies are legitimate. In practice, we tend to use either an index from our own home country agency (when available) or an index from an agency to which we subjectively attach a certain importance. Thus, in this paper, we explored and suggest an alternative approach, which would be to construct and use an aggregated risk estimate based on multiple reference indices (AREMRI) created from mathematical objectivity. Such an AREMRI would allow for direct comparisons among studies, promote transparency and ensure a consistency of analysis during risk assessment. The aim here is not to judge how researchers and/or regulatory agencies derive reference indices or to criticize existing reference indices but rather to use them as is, starting from the fact that each of the indices carries a part of the truth.

The remainder of this paper is organized as follows. The general approach for constructing an AREMRI risk distribution (uncertain risks) is presented in Section 2. Section 3 is devoted to an illustrative application of the developed approach to inhalation cancer risks. Finally, the last section, Section 4, discusses the proposed approach and considers the practical implications and implementation of the developed analysis in an interactive web application.

Methods

Aggregated risk distribution construction based on multiple reference indices (AREMRI)

Consider a discrete ensemble of indicators of uncertain risk , each based on a reference index from a different source (e.g., organization/agency/study), and each being described by a distribution function normalized to one, , with = 1, 2, … , . It is assumed that all reference indices from sources are legitimate and relevant (developed for the same purpose) and that there is no objective and/or impartial reason to choose one over another for the risk assessment. As a first step, the ensemble of indicators can be characterized by the matrix of degrees of agreement, , given by (Eq.

1):

= ⋯ ⋯ ⋮ ⋮ ⋮ ⋮ ⋯ = 1 ⋯ 1 ⋯ ⋮ ⋮ ⋮ ⋮ ⋯ 1 (1)
where = is the degree of agreement between two reference index-based distributions (indicators) and such that 0 ≤ = ≤ 1 and = 1. = 0 refers to no agreement, while = 1 indicates total agreement between the two reference index-based risk distributions. A reference index-based distribution is defined as a risk distribution resulting from the use of a given/source-specific reference index (e.g., oral slope factor from US EPA). Mathematically, = overlap" , # is obtained by calculating the overlap between two (normalized to one) distribution functions and ; a general method for calculating the overlap (or overlapping area) between two functions is described elsewhere (Petit et al. 2020).

The content of is very informative to compare reference index-based risk distribution with each other and also to guide the development of a pragmatic choice of one reference index over the others, in possible combination with other techniques such as stochastic dominance or expected utility [START_REF] Verteramo Chiu | Ranking disease control strategies with stochastic outcomes[END_REF]. Whichever approach is used, the result would be to choose the indicator (i.e., reference index) from a single source (e.g., regulatory agency/study). Consequently, we are rather in favor of constructing a composite indicator resulting from the aggregation of all indicators to reflect all of the available information. Our aim is to develop such an aggregation based on the content of and, in this case, there will be no aggregation when = $ (identity matrix), i.e., when none of the reference index-based distributions would agree with the others.

To this end, we propose to construct an aggregated risk distribution based on multiple reference indices (AREMRI), % , as a weighted linear combination of reference index-based risk distributions , as follows (Eq. 2):

% = ∑ ' ( (2)
where ' (0 ≤ ' ≤ 1 such that ∑ ' ( = 1) is a contributing weight of the risk distribution generated with the reference index " " and is a risk distribution based on a single given reference index. Next, to determine the ' coefficients and obtain an objective aggregated distribution, we require that % has a maximum agreement with all risk distributions generated with all reference indices available , i.e., to maximize total agreement ) between % and all .

The total agreement ) is related to the agreement matrix as (Eq. 3):

) = * 1 1 ⋮ 1 + , 1 ⋯ 1 ⋯ ⋮ ⋮ ⋮ ⋮ ⋯ 1 -....../......0 1 ' ' ⋮ ' (3)
The agreement matrix 1 is symmetric with the trace, Tr 1 = ∑ 3 4 4( = , where 3 4 are the eigenvalues associated with (normalized to one) eigenvectors 5 4 of matrix 1 such that 0 ≤ 3 < 3 < ⋯ < 3 ≤ . When there is no agreement among all reference index-based risk distributions , 1

reduces to the identity matrix with = 0 for 7 , and the total agreement is always ) = 1 regardless of the weights ' . In this case, no aggregation based on the agreement matrix 1 is possible as there is no objective way of choosing one reference index over another. In the opposite situation when all = 1, i.e., perfect agreement among all reference index-based risk distributions, the total agreement is maximum, ) = . We will show below that, in this case, ' = 1 ⁄ (arithmetic mean).

Finally, when 0 < < 1, the maximum of the total agreement, ), is obtained by setting the ' coefficients equal to eigenvector 5 corresponding to the largest eigenvalue 3 as shown in (Eq. 4):

' ' ⋮ ' = "∑ 9 :,; : ;<= # > , > , ⋮ > , ⟹ ) = 3 "∑ 9 :,; : ;<= # * 1 1 ⋮ 1 + , > , > , ⋮ > , = 3 (4)
This determines the weights of % in Eq. 2. It is paramount to ensure that all > , have the same sign. In case of perfect agreement among all reference index-based risk distributions (indicators), i.e., all = 1, the highest eigenvalue of the agreement matrix 1 is 3 = , associated with the eigenvector

5 , = 1, 1, ⋯ ,1 , ; this leads to ' = ' = ⋯ = ' = 1 ⁄ .
When dealing with only two reference indices (indicators), i.e., = 2 and = 7 0, the highest eigenvalue of 1 is 3 = 1 + , associated with the eigenvector 5 , = 1, 1 , . In this case, the weights are ' = ' = 1 2 ⁄ , independent of , provided that 7 0.

To gain more understanding of these weights, consider an example of = 3 for which weights can be calculated explicitly as:

= B 1 0 1 C 0 C 1 D ⟹ B ' ' ' C D = E% =F GH% =F F G% FI F G% FI J * H + C C + (5)
where = < 1 and C = C < 1 such that C = C = 0, and the highest eigenvalue of 1 is

3 C = 1 + H + C
. Clearly, Eq. 5 shows how the weights are "proportional or a function" of the respective agreements. For instance, setting C = 0 leads to ' C = 0, i.e., the contribution of a reference index-based risk distribution (indicator) is zero when the reference index-based risk distribution does not agree with any other. In addition, the weight of a reference index-based risk distribution (indicator) increases with its total agreement with the other reference index-based risk distributions (indicators). Indeed, we have ' C < ' < ' when C < , and ' < ' C < ' for C > . These observations remain valid for > 3.

Application to inhalation cancer risks

To illustrate the usefulness of using an aggregated risk distribution based on multiple reference indices, we consider three practical cases/scenarios of risk assessment of benzene exposure where each exposure scenario is characterized by a distribution of benzene concentrations assumed constant over time and following a lognormal distribution with the geometric mean (GM) and geometric standard deviation (GSD) parameters given in Table 1 and an exposure duration of 4 hours a day, 120 days a year for 50 years. Regardless of the reference index considered, the concentration distributions used for a given scenario were exactly the same (same geometric mean and geometric standard deviation).

The concentration distribution parameters differed from one practical case to another with a GM set to 253, 253 and 58.4 µg/m 3 and a GSD set to 2.80, 1.03 and 2.80 for cases 1, 2 and 3, respectively. For all scenarios, both the GM and GSD were set at pragmatic random values generated with R software 4.0.5® (R Core Team, Vienna, Austria) for Windows 10©.

We used the InCaRisk app (Petit et al. 2020)

Results

For a given practical case, each IUR-based risk distribution resulted in different cancer estimates even though the exposure concentration distributions were the same. These differences came solely from the IUR values that differed from one agency to another (Table 1). Indeed, each IUR-based risk distribution had the same shape (height and width) because of similar concentration parameters (Table 1) but was shifted from one another (different mean) by a factor that depended on the ratio LM %NO PQ ; LM %NO PQ R S (Fig. 2C, 3C and 4C).

As seen in the matrix of agreement between sources (agencies) (Fig. 2A, 3A and 4A), the degrees of agreement were heterogeneous between IUR-based risk distributions, regardless of the practical case considered. The results were similar for cases 1 and 3. Degrees of agreement were lower for the second scenario than for the first and third scenarios because the GSDs were the lowest (narrower distributions) for case 2. For all scenarios, the highest degrees of agreement, 92% for cases 1 and 3 and 81% for case 2, were found between IUR-based risk distributions from RIVM and WHO, while the lowest degrees of agreement were found between IUR-based risk distributions from Health Canada and OEHHA (30% for cases 1 and 3, 0.7% for case 2). Regardless of the practical case considered, IUR-based risk distributions from OEHHA and ANSES had a high degree of agreement with each other (84% for cases 1 & 3, 80% for case 2), a medium degree of agreement with all other IUR-based risk distributions for cases 1 & 3 (30 to 56%), and a low degree of agreement with all other IUR-based risk distributions for case 2 (0.7 to 15%). IUR-based risk distributions from the WHO, US EPA, Health Canada and RIVM had a medium to high degree of agreement with each other (66 to 92% for cases 1 & 3 and 30 to 81% for case 2).

The weights were similar for cases 1 and 3, ranging from 0.136 to 0.189, and they ranged from 0.053 to 0.247 for case 2 (Fig. 2B, 3B and 4B). The highest weights were attributed to IUR-based risk distributions that had the highest degree of agreement with the maximum number of other IUR-based risk distributions. Conversely, IUR-based risk distributions that shared a high degree of agreement with the smallest number of other IUR-based risk distributions had the smallest weights.

For the first and third scenarios, the aggregated multi IUR-based risk distribution (AREMRI) followed a lognormal distribution of parameters GM = 1.65x10 -4 and GSD = 2.8 for case 1 and GM = 3.82x10 -5

and GSD = 2.84 for case 3. The AREMRI resulted in the third highest risk, with a leukemia risk of 6.56x10 -4 [IC 95%: 2.31x10 -5 -1.32x10 -3 ] for case 1 and 1.53x10 -4 [IC 95%: 5.30x10 -6 -3.07x10 -4 ] for case 3. For both scenarios, these leukemia risks were 2.5 and 2.2 times lower than the cancer risk estimated with the reference indices (IURs) of OEHHA and ANSES, respectively (Fig. 2C and4C).

Conversely, the aggregated multi-IUR-based risks were 1.5, 1.9, 2.3 and 3.5 times higher than the leukemia risk estimated with the reference indices from the US EPA, WHO, RIVM and Health Canada, respectively. The concentration GSDs used for cases 1 & 3 were the same (GSD = 2.80), which explains why the weights and degrees of agreement were similar for both scenarios. The only difference between case 1 and case 3 came from the estimated leukemia risk, which was 4.3 times lower for case 3 than for case 1. This difference corresponded to the ratio of their concentration GM (253 µg/m 3 for case 1 vs. 58.4 µg/m 3 for case 3). In other words, for a given concentration GSD, changing the concentration GM will shift the risk estimate toward higher values when the GM is increased and toward lower values when the GM is decreased, with little to no change at all of the degrees of agreement or the weighted values.

For the second scenario, the aggregated multi-IUR-based risk distribution (AREMRI) followed a lognormal distribution of parameters GM = 1.08x10 -4 and GSD = 1.52. The AREMRI resulted in the third highest risk, with a leukemia risk of 1.83x10 -4 [IC 95%: 4.77x10 -5 -2.55x10 -4 ], which was 3.5 and 3.1 times lower than the cancer risk estimated with the reference indices (IURs) of OEHHA and ANSES, respectively (Fig. 3C). Conversely, the aggregated multi-IUR-based risk was 1.1, 1.4, 1.6 and 2.5 times higher than the leukemia risk estimated with the reference indices from the US EPA, WHO, RIVM and Health Canada, respectively. For this scenario, the GM concentration was the same as that in the first scenario (GM = 253 µg/m 3 ), but the GSD concentration was lower (1.03 for case 2 vs. 2.80 for case 1). A lower GSD concentration for case 2 resulted in lower degrees of agreement (Fig. 3A) and better risk estimation (narrower 95% confidence interval) than for case 1. Indeed, there was a 5.3 factor between the lower and upper bounds of the 95% confidence interval of the risk distribution for 237 case 2 compared to a factor of 57 for case 1. 238

Discussion

Performing a risk assessment requires the use of toxicity reference indices. When dealing with several reference indices, risk assessment can be challenging, as there is no objective way to choose one reference index over the others. By definition, all indices available from sanitary/regulatory agencies are legitimate, but in practice, we tend to use an index from a source or agency to which we subjectively attach a certain importance. In this paper, we propose a novel approach consisting of constructing and using an aggregated multi-reference index-based risk estimate constructed only on mathematical objectivity, assuming that the reference indices considered are all relevant and suitable for the risk assessment to be conducted. An R code allowing for the estimation of risks (AREMRI) by inputting multiple reference indices is provided in the Supplemental Material.

The construction of the aggregated risk distribution was based on the weighted linear aggregation (arithmetic mean) of the probability risk distribution generated using a given reference index and using the degree of agreement among the reference index-based risk distributions for the assignment of weight. Other approaches could be used, such as stochastic dominance, which is a form of stochastic ordering that would favor one reference index-based risk distribution over the others based on a specific criterion. For example, using stochastic dominance, it would be possible to favor only the risk distribution resulting from the source or agency with the highest reference index, that is to say the most protective, or with the lowest reference index (weight of 1 while a weight of 0 for all the other reference index-based risk distribution) (Verteramo [START_REF] Verteramo Chiu | Ranking disease control strategies with stochastic outcomes[END_REF]. That kind of approach is, however, different from our goal, which was not to judge or provide arguments for choosing a reference index over the others but rather to objectively construct an aggregated risk indicator that takes into account several reference indices from different sources (e.g., sanitary agencies), as all reference indices are legitimate and carry part of the truth.

The proposed approach was illustrated using three practical cases where the exposure concentration distributions were assumed to follow lognormal distributions. Our approach is not limited to lognormal distributions and can work with any type of distribution. Two additional examples are available in the Supplemental Material (Fig. S1 andS2), with exposure concentrations following normal and Poisson distributions, respectively.

The approach consisting of constructing and using an aggregated multi-reference index-based risk estimate was illustrated for one kind of stochastic reference value (inhalation cancer unit risk). The presented approach could be used with other types of toxicological reference values, such as oral slope factors, drinking water unit risks or threshold/deterministic values. For some toxicological reference values, the distribution of the underlying toxicity data is available, and the reference value itself can sometimes correspond to a certain point on that distribution (e.g., for BMDL values). In that particular case, these reference values could potentially be set as distributions instead of single values using the distribution of the underlying toxicity data when available. While the presented approach can still be used with reference value distributions, it should be noted that the choice of the reference value distribution (e.g., gamma, lognormal, Poisson, Weibull) and the associated uncertainties (confidence intervals) will both have an important impact on the resulting risk estimate. In particular, in the presence of large confidence intervals (high data dispersion), the degrees of agreement (overlaps between distributions of risks) are so high that all of the reference indices contribute almost equally to the aggregate distribution and the notion of multi-reference indices becomes less relevant. This issue is particularly true with deterministic values (threshold values), for which the choice of the derivation method and, in particular, the choice and use of large safety factors may dominate the differences between the toxicological reference values. One possible solution to address this issue would be to replace traditional noncancer reference doses (threshold values) with probabilistic estimates [START_REF] Chiu | Beyond the RfD: Broad Application of a Probabilistic Approach to Improve Chemical Dose-Response Assessments for Noncancer Effects[END_REF].

Regardless of the practical case considered, IUR-based risk distributions from OEHHA and ANSES had a high degree of agreement with each other but a medium to low degree of agreement with other IUR-based risk distributions. The six benzene reference values used in this illustration were derived for the same purpose (lifetime benzene exposure), all based on cohort studies of benzene-exposed workers, but using different key studies and extrapolation/derivation methods among agencies, which may explain the observed differences (ANSES 2019;Health Canada 2010;[START_REF] Oehha | Consolidated table of OEHHA/ARB approved risk assessment health values[END_REF]RIVM 2001;[START_REF] Epa | Table[END_REF][START_REF] Who | Air Quality Guidelines for Europe[END_REF]. OEHHA and ANSES used the same extrapolation method, a relatively risk procedure, but chose two different key studies, [START_REF] Rinsky | Leukemia in benzene workers[END_REF] for OEHHA and [START_REF] Richardson | Temporal variation in the association between benzene and leukemia mortality[END_REF] for ANSES. RIVM chose the same key study as OEHHA but used a nonthreshold extrapolation method. The WHO used an average relative risk model with [START_REF] Crump | Risk of benzene-induced leukemia: a sensitivity analysis of the pliofilm cohort with additional follow-up and new exposure estimates[END_REF] as a key study. The US EPA applied a low-dose linearity utilizing the maximum likelihood estimates method based on six key studies [START_REF] Crump | Quantitative estimates of risk of leukemia from occupational exposure to benzene[END_REF][START_REF] Crump | Risk of benzene-induced leukemia: a sensitivity analysis of the pliofilm cohort with additional follow-up and new exposure estimates[END_REF][START_REF] Paustenbach | Benzene toxicity and risk assessment, 1972-1992: implications for future regulation[END_REF][START_REF] Rinsky | Leukemia in benzene workers[END_REF][START_REF] Rinsky | Benzene and leukemia -An epidemiologic risk assessment[END_REF][START_REF] Epa | Carcinogenic effects of benzene: an update[END_REF]. Health Canada applied a linear quadratic model of the exposure response relationship based on four key studies [START_REF] Bond | An update of mortality among chemical workers exposed to benzene[END_REF][START_REF] Rinsky | Benzene and leukemia -An epidemiologic risk assessment[END_REF]Wong 1987a,b).

Differences in policies and/or expert judgments could also have played a role in the observed differences [START_REF] Beck | Approaches for describing and communicating overall uncertainty in toxicity characterizations: U.S. Environmental Protection Agency's Integrated Risk Information System (IRIS) as a case study[END_REF].

The aggregated multi-reference index-based distribution (the degrees of agreement and weights) depends on the number of reference indices included in the aggregation as well as the concentration distribution parameters. The narrower the concentration distributions (lower GSD), the lower the degrees of agreement and the better the risk distribution estimation (narrower confidence interval). If the GSD concentration tends to 1, as in the second scenario, the AREMRI tends to the mean risk value and also tends to approach and equal a single reference index-based risk distribution (e.g., the risk distribution using the reference value from the US EPA for case 2).

Regarding the reference indices, we found, for our practical cases, six agencies that set/proposed an IUR for benzene and the risk of leukemia. There are most likely other agencies and studies around the world which we are not aware of that may also have proposed an IUR for benzene and the risk of leukemia for lifetime exposure. Integrating these potential other IURs into the aggregated multireference index-based risk distribution construction would change the leukemia risk obtained with the aggregated risk indicator. However, the approach presented in this paper does not remove the requirement of a prior screening of reference values to be included in the integrated approach. This screening should always be made to determine how reference indices were derived, which data they are based on, and their intended use (e.g., assessment of occupational exposure, exposure via foodstuffs and unintentional exposure). Indeed, it is paramount that the risk assessor has some knowledge regarding the reference values used in an assessment and understands the way the values were derived, their associated uncertainties and their relevance and validity for the risk assessment being carried out. If the approach presented in this paper is used to interpret the results of a risk assessment, then only reference values derived for the same purpose(s) should be included in the analysis.

Conclusions

The three practical cases highlighted the usefulness of comparing risk estimates generated using multiple reference indices from several sources (e.g., regulatory agencies) when assessing the human health risk associated with exposure to environmental pollutants.

Therefore, the proposed approach hinges on its ability to provide useful information when conducting risk assessment in the presence of several reference indices from various sources. These results are essential for risk managers, practitioners, and decision-makers by providing guidance and by supporting direct comparisons between studies. Indeed, our approach can provide a more complete, consistent and transparent risk assessment of pollutants that pose a threat to human health and support better-informed risk management decisions by giving some insight into the uncertainties associated with the choice of the reference index in the risk assessment. Our approach could also serve to integrate risk assessments, such as by obtaining an integrated view of risk assessments when trying to prioritize among different cases for remedial action.

To facilitate the use of our approach, the matrix of degrees of agreement, side-by-side comparison of reference index-based risk distributions, the choice of the reference values to consider and the construction of the aggregated multi-reference index-based risk distribution will be implemented in the InCaRisk app for inhalation cancer risk assessment. In addition to making information more readily understood and retained in a quicker time, the visual aids and representations provided by InCaRisk could be useful for improving risk communication and promoting transparency among studies when conducting risk assessments [START_REF] Beck | Approaches for describing and communicating overall uncertainty in toxicity characterizations: U.S. Environmental Protection Agency's Integrated Risk Information System (IRIS) as a case study[END_REF].

Regardless of the type of exposure (e.g., acute, subchronic, or chronic), absorption pathways, exposure setting (environmental or occupational), stages of the life cycle of a chemical (e.g., manufacture, use, disposal, consumer products or waste) and stages of the life of an individual (e.g., adulthood), the approach presented in this paper can be used. However, this approach can only be used with reference indices intended for the same purpose (e.g., the same absorption pathway, same type of exposure and same exposure setting), and some prior screening of reference indices has to be done to ensure that the aggregated approach does not look at a mixture of reference values that should not be integrated.

Finally, the proposed approach has applicability in other fields and studies with similar problematics, as long as there is some constraint on the choice of reference values to be included. Table 1: Summary of the lognormal exposure concentration distribution parameters used for three practical cases

Tables

Figure captions

Fig. 1

 1 presents a schematic summary of the different steps of the construction of the aggregated risk distribution based on multiple reference indices.

Fig. 1 :

 1 Fig. 1: Scheme of the general approach to construct an aggregated risk distribution based on multiple

Fig. 2 :

 2 Fig. 2: Side-by-side comparison of single IUR-based risk distributions for the construction of the aggregated multi-IUR-based risk distribution -case 1

Fig. 3 :

 3 Fig. 3: Side-by-side comparison of single IUR-based risk distributions for the construction of the aggregated multi-IUR-based risk distribution -case 2

Fig. 4 :

 4 Fig. 4: Side-by-side comparison of single IUR-based risk distributions for the construction of the aggregated multi-IUR-based risk distribution -case 3

Fig. 1 :

 1 Fig. 1: Scheme of the general approach to construct an aggregated risk distribution based on multiple

Fig. 2 :

 2 Fig. 2: Side-by-side comparison of single IUR-based risk distributions for the construction of the

Fig. 3 :

 3 Fig. 3: Side-by-side comparison of single IUR-based risk distributions for the construction of the

  with 1,000 simulations runs to generate all of the reference index-based risk distributions . Briefly, the InCaRisk app is a free user-friendly and interactive web application that allows for cancer risk estimation following inhalation exposure

	Environmental Protection Agency) and WHO (World Health Organization). The IURs for benzene set
	by these agencies regarding the risk of leukemia for lifetime exposure were 2.6x10 -5 (ANSES 2019),
	3.3x10 -6 (Health Canada 2010), 2.9x10 -5 (OEHHA 2019), 5x10 -6 (RIVM 2001), 7.8x10 -6 (US EPA
	2018) and 6x10 -6 m 3 /µg (WHO 2000).

(https://exporisk-timc.imag.fr/InCaRisk/). A total of 305 substances/chemicals with reference values from up to eight agencies are available. InCaRisk requires no programming knowledge to use. This app also incorporates a variety of features and options to make it easy to use, with the possibility of configuring the exposure settings (exposure scenario), the concentration distribution settings (choice of the type of distribution, distribution parameters, uploading an exposure concentration file) and to export the results as a report. This app offers the chance to have an immediate glimpse into the results and to see how the results change according to different setting configurations with interactive easy-toread graphs, which are downloadable. The reference index-based risk distributions were generated for six sanitary agencies for which the tumor site (leukemia) and tumor type (hematologic/immune) used in the inhalation unit risk (IUR) establishment were the same: ANSES (French Agency for Food, Environmental and Occupational Health & Safety), Health Canada, OEHHA (California Office of Environmental Health Hazard Assessment), RIVM (National Dutch Institute for Public Health and the Environment), US EPA (US Table 1 here (see Table 1 after the Reference section).

Acknowledgements: None.

Funding

This research did not receive a specific grant from agencies in the public, commercial, or not-for-profit sectors.

Author contributions

Pascal

Declaration of interests:

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.

Protection Agency, WHO: World Health Organization. Each agency name refers to the risk distribution generated using its benzene IUR and the exposure scenario considered. The term "Interagency" refers here to the aggregated risk estimate based on multiple reference indices (AREMRI).


Marion Cheron
email: cheron.marion@gmail.com

David Costantini

Frédéric Angelier

Cécile Ribout

François Brischoux

David Costantini

Aminomethylphosphonic acid (AMPA) alters oxidative status during embryonic development in an amphibian species

Keywords: aminomethylphosphonic acid, Bufo spinosus, oxidative stress, telomere

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Among the numerous sources of environmental contamination (e.g., fossil fuel combustion, industrial activities, [START_REF] Kendall | Wildlife Toxicology : Emerging Contaminant and Biodiversity Issues[END_REF], modern agriculture is well-known for the massive release of toxic substances in the environments [START_REF] Jayaraj | Review Article. Organochlorine pesticides, their toxic effects on living organisms and their fate in the environment[END_REF][START_REF] Solomon | Ecological Risk Assessment for Aquatic Organisms from Over-Water Uses of Glyphosate[END_REF]. Indeed, pesticides are specifically designed to alter the functioning of diverse organisms (weeds, fungi, and insects) which are detrimental to crop productivity [START_REF] Sharma | Worldwide pesticide usage and its impacts on ecosystem[END_REF]. Although traditional toxicology studies have convincingly demonstrated the potential toxicity of such compounds on non-target organisms [START_REF] Crump | Exposure to the Herbicide Acetochlor Alters Thyroid Hormone-Dependent Gene Expression and Metamorphosis in Xenopus Laevis[END_REF][START_REF] De Brito Rodrigues | Impact of the glyphosate-based commercial herbicide, its components and its metabolite AMPA on non-target aquatic organisms[END_REF][START_REF] Hasenbein | Environmentally relevant concentrations of herbicides impact non-target species at multiple sublethal endpoints[END_REF], they are not sufficient to assess their impact on wild animals for several reasons. First, they often rely on relatively high, environmentally irrelevant concentrations of agrochemicals in order to assess mortality thresholds (i.e., LD50, [START_REF] Relyea | What Doesn't Kill You Makes You Sluggish: How Sublethal Pesticides Alter Predator-Prey Interactions[END_REF]. As a consequence, reliance on LD50 limits our ability to identify more subtle, sublethal effects occurring at lower concentrations [START_REF] Beiras | Sublethal Toxicity at the Level of Organism[END_REF]. However, such sublethal effects have been shown to occur at extremely low concentrations and to be responsible for the disruption of several organismal traits that can affect population persistence [START_REF] Cauble | Sublethal Effects of the Herbicide Glyphosate on Amphibian Metamorphosis and Development[END_REF][START_REF] Cheron | Aminomethylphosphonic acid alters amphibian embryonic development at environmental concentrations[END_REF][START_REF] Lenkowski | Low concentrations of atrazine, glyphosate, 2, 4-dichlorophenoxyacetic acid, and triadimefon exposures have diverse effects on Xenopus laevis organ morphogenesis[END_REF].

Importantly, field monitoring often reveals that environmental contamination occurs at concentrations that are well below LD50 [START_REF] Awkerman | Simulated developmental and reproductive impacts on amphibian populations and implications for assessing long-term effects[END_REF][START_REF] Egea-Serrano | Understanding of the impact of chemicals on amphibians: a meta-analytic review: Impact of Pollution on Amphibians[END_REF]. Second, most toxicology studies use laboratory models which can differ from their wild counterparts in their response to environmental contamination [START_REF] Hahn | Casting a wide net: use of diverse model organisms to advance toxicology[END_REF][START_REF] Patisaul | Animal models of endocrine disruption[END_REF]. Indeed, laboratory conditions simplify the complexity that organisms face in natural environments and long-term selection of laboratory lineages are likely to influence the responses of these species to environmental constraints [START_REF] Hahn | Casting a wide net: use of diverse model organisms to advance toxicology[END_REF][START_REF] Segner | What constitutes a model organism in ecotoxicology?[END_REF]. Overall, one of the prominent current challenges is to assess the detrimental sublethal effects of contaminants at environmentally relevant concentrations on non-model organism's representative of exposure in the wild. This is particularly relevant for pesticides, because their breakdown after spread in agricultural fields can expose organisms to concentrations much lower than those originally used or tested in laboratory experiments [START_REF] Bonansea | The Fate of Glyphosate and AMPA in a Freshwater Endorheic Basin: An Ecotoxicological Risk Assessment[END_REF][START_REF] Wang | A review of pesticide fate and transport simulation at watershed level using SWAT: Current status and research concerns[END_REF].

Glyphosate remains one of the most widely used non-selective herbicide [START_REF] Clapp | Explaining Growing Glyphosate Use: The Political Economy of Herbicide-Dependent Agriculture[END_REF][START_REF] Helander | Glyphosate in northern ecosystems[END_REF]. Accordingly, the effects of glyphosate on wildlife, including investigations of commercial formulations and/or of different surfactants have been relatively thoroughly studied [START_REF] Berger | How Does Changing Pesticide Usage Over Time Affect Migrating Amphibians: A Case Study on the Use of Glyphosate-Based Herbicides in German Agriculture Over 20 Years[END_REF][START_REF] Bókony | Chronic exposure to a glyphosate-based herbicide makes toad larvae more toxic[END_REF][START_REF] Cauble | Sublethal Effects of the Herbicide Glyphosate on Amphibian Metamorphosis and Development[END_REF][START_REF] De Brito Rodrigues | Impact of the glyphosate-based commercial herbicide, its components and its metabolite AMPA on non-target aquatic organisms[END_REF][START_REF] Giaquinto | Effects of Glyphosate-Based Herbicide Sub-Lethal Concentrations on Fish Feeding Behavior[END_REF][START_REF] Gill | Glyphosate toxicity for animals[END_REF][START_REF] Guilherme | European eel (Anguilla anguilla) genotoxic and pro-oxidant responses following short-term exposure to Roundup(R)--a glyphosate-based herbicide[END_REF][START_REF] Herek | Can environmental concentrations of glyphosate affect survival and cause malformation in amphibians? Effects from a glyphosate-based herbicide on Physalaemus cuvieri and P. gracilis (Anura: Leptodactylidae)[END_REF][START_REF] Lenkowski | Low concentrations of atrazine, glyphosate, 2, 4-dichlorophenoxyacetic acid, and triadimefon exposures have diverse effects on Xenopus laevis organ morphogenesis[END_REF][START_REF] Matozzo | The Effects of Glyphosate and Its Commercial Formulations to Marine Invertebrates: A Review[END_REF][START_REF] Muñoz | Glyphosate and the key characteristics of an endocrine disruptor: A review[END_REF][START_REF] Solomon | Ecological Risk Assessment for Aquatic Organisms from Over-Water Uses of Glyphosate[END_REF]. However, its primary metabolite, aminomethylphosphonic acid (AMPA) is arguably the most common compound found in surface and groundwater worldwide [START_REF] Bonansea | The Fate of Glyphosate and AMPA in a Freshwater Endorheic Basin: An Ecotoxicological Risk Assessment[END_REF][START_REF] Medalie | Influence of land use and region on glyphosate and aminomethylphosphonic acid in streams in the USA[END_REF][START_REF] Okada | Glyphosate and aminomethylphosphonic acid (AMPA) are commonly found in urban streams and wetlands of Melbourne, Australia[END_REF][START_REF] Silva | Distribution of glyphosate and aminomethylphosphonic acid (AMPA) in agricultural topsoils of the European Union[END_REF].

The greater occurrence and abundance of AMPA in natural environments is likely linked to the extensive use of glyphosate and also because AMPA is an organic phosphonate derived from detergents [START_REF] Goto | AminoMethylPhosphonic acid (AMPA) in natural waters: Its sources, behavior and environmental fate[END_REF]. In North American wetland, AMPA occurred at 97% ranging from 0.25 to 10 µg.l -1 [START_REF] Okada | Glyphosate and aminomethylphosphonic acid (AMPA) are commonly found in urban streams and wetlands of Melbourne, Australia[END_REF]. Environmental concentrations found in aquatic environments in France range from 0.1 μg l -1 to 6.6 μg l -1 (data from Water Agencies, "Agence de l'eau Loire Bretagne" and "Agence de l'eau Adour-Garonne"). Despite its greater environmental occurrence and abundance than glyphosate, studies on the effects of AMPA at environmental concentrations on wildlife are comparatively fewer [START_REF] Cheron | Aminomethylphosphonic acid alters amphibian embryonic development at environmental concentrations[END_REF][START_REF] De Brito Rodrigues | Impact of the glyphosate-based commercial herbicide, its components and its metabolite AMPA on non-target aquatic organisms[END_REF][START_REF] Domínguez | Toxicity of AMPA to the earthworm Eisenia andrei Bouché, 1972 in tropical artificial soil[END_REF][START_REF] Guilherme | DNA and chromosomal damage induced in fish (Anguilla anguilla L.) by aminomethylphosphonic acid (AMPA)-the major environmental breakdown product of glyphosate[END_REF][START_REF] Martinez | Effects of glyphosate and aminomethylphosphonic acid on an isogeneic model of the human blood-brain barrier[END_REF][START_REF] Matozzo | Ecotoxicological hazard of a mixture of glyphosate and aminomethylphosphonic acid to the mussel Mytilus galloprovincialis (Lamarck 1819)[END_REF][START_REF] Matozzo | Effects of aminomethylphosphonic acid, the main breakdown product of glyphosate, on cellular and biochemical parameters of the mussel Mytilus galloprovincialis[END_REF]. Importantly, these studies have revealed sublethal effects on wild species, such as genotoxic effect [START_REF] Battaglin | Glyphosate and Its Degradation Product AMPA Occur Frequently and Widely in U.S. Soils, Surface Water, Groundwater, and Precipitation[END_REF][START_REF] Medalie | Influence of land use and region on glyphosate and aminomethylphosphonic acid in streams in the USA[END_REF][START_REF] Okada | Glyphosate and aminomethylphosphonic acid (AMPA) are commonly found in urban streams and wetlands of Melbourne, Australia[END_REF][START_REF] Silva | Distribution of glyphosate and aminomethylphosphonic acid (AMPA) in agricultural topsoils of the European Union[END_REF] and cellular and biochemical alterations [START_REF] Matozzo | Ecotoxicological hazard of a mixture of glyphosate and aminomethylphosphonic acid to the mussel Mytilus galloprovincialis (Lamarck 1819)[END_REF][START_REF] Matozzo | Effects of aminomethylphosphonic acid, the main breakdown product of glyphosate, on cellular and biochemical parameters of the mussel Mytilus galloprovincialis[END_REF]. AMPA has been shown to alter embryonic development in amphibians, increasing embryonic mortality and altering hatchling phenotype through processes occurring at low concentrations and suggesting nonmonotonicity [START_REF] Cheron | Aminomethylphosphonic acid alters amphibian embryonic development at environmental concentrations[END_REF].

Our goal was to identify potential physiological mechanisms that may be responsible for the alteration of embryonic development due to environmentally relevant concentrations of AMPA [START_REF] Cheron | Aminomethylphosphonic acid alters amphibian embryonic development at environmental concentrations[END_REF]. In this respect, we focused our investigations on oxidative stress, as such cellular biochemical process is known to reflect how environmental stressors affect organismal physiology [START_REF] Costantini | Oxidative Stress and Hormesis in Evolutionary Ecology and Physiology[END_REF]. The metabolic activity of cells is responsible for a continuous generation of reactive chemical species (e.g. reactive oxygen species ROS), which cause oxidation of numerous organic substrates (e.g. lipids, DNA).

Organisms are equipped with a complex array of molecular defence mechanisms, generally referred to as antioxidants. The cellular oxidative status is, therefore, determined by a complex network of molecules, including oxidative damages to organic substates, levels of both enzymatic and non-enzymatic antioxidants, or ROS production [START_REF] Costantini | Understanding diversity in oxidative status and oxidative stress: the opportunities and challenges ahead[END_REF]. Any changes in these molecular components of cellular oxidative status provide evidence for oxidative unbalance. There is growing consensus that oxidative status alteration is one relevant molecular mechanism that links environmental conditions (including contaminant exposure) to fitness, through effects on growth, fertility, reproduction and senescence [START_REF] Alonso-Alvarez | Increased susceptibility to oxidative stress as a proximate cost of reproduction: Oxidative stress as a cost of reproduction[END_REF][START_REF] Beaulieu | Biomarkers of oxidative status: missing tools in conservation physiology[END_REF][START_REF] Costantini | Oxidative stress in relation to reproduction, contaminants, gender and age in a long-lived seabird[END_REF][START_REF] Stier | Constraint and cost of oxidative stress on reproduction: correlative evidence in laboratory mice and review of the literature[END_REF]. In addition to markers of oxidative status, we included in our investigations telomere length that may reveal complementary effects of AMPA on DNA damage during amphibian embryonic development [START_REF] Burraco | Different effects of accelerated development and enhanced growth on oxidative stress and telomere shortening in amphibian larvae[END_REF][START_REF] Chatelain | The association between stressors and telomeres in non-human vertebrates: a meta-analysis[END_REF][START_REF] Cheron | Clutch quality is related to embryonic development duration, hatchling body size and telomere length in the spined toad ( Bufo spinosus )[END_REF]. Indeed, telomeres are repetitive nucleotide sequences located in each chromosome to ensure stability and viability, and telomere length is increasingly suggested to be an integrative biomarker of organismal, cellular, and molecular damages. Although telomere shortening occurs normally in cell division, this process can be magnified under stressful conditions especially when they are associated with oxidative stress [START_REF] Reichert | Does oxidative stress shorten telomeres in vivo ? A review[END_REF].

Accordingly, stress-induced telomere attrition may be linked to an intense cellular oxidative stress [START_REF] Chatelain | The association between stressors and telomeres in non-human vertebrates: a meta-analysis[END_REF] because guanine (the main nucleotide of telomeres) is particularly sensitive to oxidation [START_REF] Bjelland | Mutagenicity, toxicity and repair of DNA base damage induced by oxidation[END_REF] . Finally, telomere length might be seen as a proxy of individual quality although such relationship is not well establish in ectotherms species (Burraco et al., 2020a;[START_REF] Olsson | Ectothermic telomeres: it's time they came in from the cold[END_REF] Specifically, we hypothesized that environmentally relevant concentrations of AMPA, which have been shown to induce altered embryonic mortality, development duration and hatchling phenotype [START_REF] Cheron | Aminomethylphosphonic acid alters amphibian embryonic development at environmental concentrations[END_REF], will induce significant changes in oxidative status. We expected that exposition to AMPA concentrations will induce increase in antioxidant enzymes such as SOD, GPx and CAT and a depletion of thiols levels.

Moreover, we expected non monotonic relationships between oxidative markers and AMPA concentrations given that some studies showed nonlinear or hormetic effects of AMPA or glyphosate on organisms [START_REF] Brito | Hormetic effects of glyphosate on plants: Glyphosate hormesis[END_REF][START_REF] Cheron | Aminomethylphosphonic acid alters amphibian embryonic development at environmental concentrations[END_REF][START_REF] Domínguez | Toxicity of AMPA to the earthworm Eisenia andrei Bouché, 1972 in tropical artificial soil[END_REF]. If exposure to AMPA leads to increased oxidative stress, we also predicted that exposure to AMPA should be associated with shorter telomeres.

Materials and methods

Study species

Spined toad (Bufo spinosus) is a widespread amphibian species that persists in agricultural landscapes [START_REF] Gagné | Effect of landscape context on anuran communities in breeding ponds in the National Capital Region, Canada[END_REF][START_REF] Guillot | Landscape influences the morphology of male common toads ( Bufo bufo )[END_REF]. Breeding pond (46°090N, 0°240W, France) was monitored every day during the reproductive period. All eggs strings (N=10) were collected during a single day in January 2018 immediately after laying in order to avoid development in the natural environment. A subset of each string containing 120 eggs was randomly selected and kept for our experiment (see below), while the remaining eggs (i.e., 3000-5000) were released at their site of origin.

Treatment concentrations and chemical solutions

Environmental concentrations found in aquatic environments in France range from 0.1 μg l -1 to 6.6 μg l -1 (data from Water Agencies, "Agence de l'eau Loire Bretagne" and "Agence de l'eau Adour-Garonne") and are similar to concentrations found in aquatic environments in Northern America and Europe [START_REF] Bonansea | The Fate of Glyphosate and AMPA in a Freshwater Endorheic Basin: An Ecotoxicological Risk Assessment[END_REF][START_REF] Coupe | Fate and transport of glyphosate and aminomethylphosphonic acid in surface waters of agricultural basins[END_REF][START_REF] Goto | AminoMethylPhosphonic acid (AMPA) in natural waters: Its sources, behavior and environmental fate[END_REF]. Because our goal was to mimic the range of concentrations found in the wild, we aimed at three nominal concentrations spanning from low to relatively high concentrations (0.1 µg.l -1 ; 0.30µg.l -1 and 3.00 µg.l -1 ,Table 1). Stock solutions (0.1 g/L) of AMPA were obtained by dissolving commercial crystalline powder (aminomethylphosphonic acid, 99% purity, ACROS ORGANICS TM ) in dechlorinated tap water. These stock solutions were further diluted with dechlorinated tap water to reach the selected concentration for each treatment group (Table 1). We produced additional water samples for each treatment in order to perform analytical verifications of the concentrations. These analytical verifications of the actual concentrations were performed by an independent accredited analytical laboratory (QUALYSE, Champdeniers-Saint-Denis, France). Determination of AMPA in water samples was assessed using liquid chromatography-mass spectrometry (LC/MS/MS) with 9fluorenylmethyl chloroformate (FMOC-Cl) used as derivatization agent. The limit of quantification was 0.1 µg.l -1 for the determination of AMPA in water samples. Overall, these verifications showed that actual concentrations were similar to concentrations found in aquatic environments (see above) and that the differences between "low" (0.07±0.01 µg.l -1 ), "medium"(0.32±0.052 µg.l -1 ) and "high"(3.57±0.153 µg.l -1 ) treatments were significant (Table 1).

For clarity, we will refer to the treatments as Low, Medium and High hereafter.

Experimental design

Each subset of each clutch (containing 120 eggs) was further separated into segments containing 30 eggs (i.e., 120/4 treatments) that were randomly placed in different glassmolded aquariums containing 2 L of water and one of our experimental treatments (i.e., Control, Low, Medium or High) and monitored until hatching (total N = 4 treatments * 10 clutches representing 60 tanks containing 30 eggs each).

Eggs were subjected to different concentrations of AMPA (Low, Medium and High) and a control during the whole embryonic development until hatching. To maintain relatively constant exposure levels throughout the experiment but to avoid excessive mechanical disturbances to developing eggs, water was changed once a week according to the half-life of AMPA in water [START_REF] Battaglin | Glyphosate and Its Degradation Product AMPA Occur Frequently and Widely in U.S. Soils, Surface Water, Groundwater, and Precipitation[END_REF]. Our procedures did not alter egg jelly matrix surrounding the developing embryos throughout the experiment. We determined the stage of development according to the developmental table of [START_REF] Gosner | A Simplified Table for Staging Anuran Embryos and Larvae with Notes on Identification[END_REF]. Eggs were exposed from laying to hatching and hatching occurred at Gosner stage 22 (hereafter GS 22) after 16.10±0.16 days (see Table 2). All experiments took place in a thermal-controlled room with air and water temperature set at 17°C and with a natural day/night photoperiod (9h/15h).

We monitored egg segments and counted all the individuals that hatched (GS22) which we considered a metric of hatching success. Development duration was assessed as the time elapsed between the date of egg collection (egg laying) and the date of hatching. We took pictures of all live hatchlings using a camera (Panasonic Lumix DC-TZ55) in order to assess the length of individuals (N=80). Morphological measurements were performed with the software ImageJ [START_REF] Schneider | NIH Image to ImageJ: 25 years of Image Analysis[END_REF].

Upon hatching, for each clutch and within each treatment, we randomly collected 2 hatchlings that were pooled to assess oxidative status (because a single individual did not yield enough tissue for these assays) and 2 hatchlings that were used for telomere analyses.

Hence, we measured clutch levels of oxidative status parameters.

Oxidative status markers

We randomly selected 80 hatchlings (20/experimental group), including two hatchlings that hatched from each single clutch. We pooled together these two hatchlings because a smaller portion did not provide enough tissue for the analyses. Thus, we had 10 pools per experimental group (i.e. 40 samples in total). We homogenized hatchlings in Dulbecco's Phosphate Buffered Saline (Sigma-Aldrich, France) supplemented with 1 mM of phenylmethylsulfonyl fluoride (Sigma-Aldrich, France) as an inhibitor of proteases using a TissueLyser II (Qiagen) at 30 Hz for 4 minutes. Afterwards, we centrifuged tubes for 10 minutes at 4°C to obtain clean supernatants to be used for the assays. We measured the concentration of thiols using the -SHp test (Diacron International, Grosseto, Italy); the activity of the antioxidant enzyme superoxide dismutase (SOD) using the Ransod assay (RANDOX Laboratories, France); the activity of the antioxidant enzyme glutathione peroxidase (GPx) using the Ransel assay (RANDOX Laboratories, France); the activity of CAT using the OxiSelect CAT Activity Assay (Euromedex, France). We ran all analyses in duplicate and an average value of each duplicate was used for the statistical analyses (coefficients of variation < 10%). We normalized values of markers to the amount of proteins as quantified using the Bradford protein assay with albumin as reference standard (Sigma-Aldrich, France). All assays were run according to manufacturer's instructions and have been validated in amphibian species [START_REF] Liparoto | Biogeographic history moulds population differentiation in ageing of oxidative status in an amphibian[END_REF].

In addition to markers alone we chose to show balanced or unbalanced ratio of markers activities following de Haan et al. studies (1992[START_REF] De Haan | Elevation in the Ratio of Cu/Zn-Superoxide Dismutase to Glutathione Peroxidase Activity Induces Features of Cellular Senescence and This Effect Is Mediated by Hydrogen Peroxide[END_REF][START_REF] De Haan | Elevation in the Ratio of Cu/Zn-Superoxide Dismutase to Glutathione Peroxidase Activity Induces Features of Cellular Senescence and This Effect Is Mediated by Hydrogen Peroxide[END_REF], 1992). An unbalanced ratio of superoxide dismutase and glutathione peroxidase (SOD/GPx) can lead to a proliferation of H2O2 leading to accelerated cellular senescence and altered morphology [START_REF] De Haan | Elevation in the Ratio of Cu/Zn-Superoxide Dismutase to Glutathione Peroxidase Activity Induces Features of Cellular Senescence and This Effect Is Mediated by Hydrogen Peroxide[END_REF]. Moreover, since both GPx and catalase convert H2O2 into H2O, an unbalanced ratio of superoxide dismutase on GPx and catalase activities (SOD/(GPx+CAT) can lead to a proliferation of deleterious pro-oxidants [START_REF] De Haan | Cu/Zn superoxide dismutase mRNA and enzyme activity, and susceptibility to lipid peroxidation, increases with aging in murine brains[END_REF].

Relative telomere length

Similarly to sampling for oxidative stress markers, we randomly selected two hatchlings among the embryos that hatched of each clutch for telomere analysis (N=20 hatchlings from 10 clutches). Although, smaller portions of individuals (e.g., head, tail, muscle) can be used for telomere assays [START_REF] Burraco | Different effects of accelerated development and enhanced growth on oxidative stress and telomere shortening in amphibian larvae[END_REF][START_REF] Sánchez-Montes | Telomere attrition with age in a wild amphibian population[END_REF], hatchling spined toads were too small to provide enough tissue for these analyses and we used whole individuals [START_REF] Cheron | Clutch quality is related to embryonic development duration, hatchling body size and telomere length in the spined toad ( Bufo spinosus )[END_REF]). We used quantitative PCR (qPCR; BioRad CFX 96, Bio-Rad, USA) to determine telomere length using an established protocol [START_REF] Mclennan | Distinct telomere differences within a reproductively bimodal common lizard population[END_REF], which was adapted for spined toad [START_REF] Cheron | Clutch quality is related to embryonic development duration, hatchling body size and telomere length in the spined toad ( Bufo spinosus )[END_REF]. Whole hatchlings were digested with proteinase K and DNA was extracted using the Nucleospin Tissue Kit (Macherey-Nagel), following the manufacturer's instructions. DNA concentration and purity were assessed with a Nanodrop ND1000 spectrophotometer (Thermo Scientific). Universal telomere primers were used, and the control single-copy gene Recombination Activating Gene 1 (RAG-1) was selected and amplified using specific primers designed for the spined toad: RAG1-F 5'-GGGTCCTCTGATAGCCGAAA-3' and RAG1-R 5'-CATCATAACCTGTACCCCGGA-3'. This single-copy gene has previously been successfully used in multiple species (birds [START_REF] Sebastiano | Exposure to PFAS is Associated with Telomere Length Dynamics and Demographic Responses of an Arctic Top Predator[END_REF]; fish [START_REF] Dupoué | Chronic water restriction triggers sex-specific oxidative stress and telomere shortening in lizards[END_REF][START_REF] Mclennan | Distinct telomere differences within a reproductively bimodal common lizard population[END_REF][START_REF] Petitjean | Intraspecific variability of responses to combined metal contamination and immune challenge among wild fish populations[END_REF]) including amphibians [START_REF] Canestrelli | Biogeography of telomere dynamics in a vertebrate[END_REF]. All qPCRs were performed on three plates for each gene (RAG-1 and telomere) using 7.5 ng of DNA per reaction. The telomere and single-copy gene primers were respectively used at a concentration of 800 nM and 300 nM. To generate a six-point standard curve (from 50.0 ng to 1.5 ng) for controlling the amplifying efficiency of the reactions, serial dilutions of DNA from a pooled sampled of 10 tadpoles were included on the plate (in triplicates). For both telomere and RAG-1 amplification, the melt curves displayed single sharp curves, validating amplification specificity of these primers. A reference toad sample was run in triplicate in all plates to account for inter-plate variation. All samples were run in duplicates and were randomly distributed across the PCR plates. Amplification efficiencies reached Mean ± SE: RAG-1, 87.48 ± 5.92; TEL, 92.60 ± 6.00 and all R² were always very high (>0.99). The relative telomere length (expressed as T/S ratio) was calculated as the telomere copy number (T) relative to single-copy gene (S; RAG1), according to Cawthon (2002). Intra-plate variation for telomere and RAG1 were 2.07% and 2.53% respectively. Intra-plate variation for the T/S ratio was 3.46%. Inter-plate variations for telomere and RAG-1 were 2.07% and 1.89%, respectively. Inter-plate variation for the T/S ratio was 5.62%.

Statistical analysis

All statistical analyses were conducted with R statistical software v.4.0.0 (R Core Team, 2020) and RStudio v 1.2.5042 (RStudio, Inc.).

All data were tested for homogeneity of variance and normality with the Barlett's test and the Shapiro-Wilks test, respectively. We also checked the normality of the residuals using diagnostic plots. We checked which distribution fitted the best the model using Cullen & Frey plot (bootstrap at 500, "fitdistplus" package (Delignette-Muller and Dutang, 2015)). We computed z-scores in all predictor variables (scaled to mean of 0 and standard deviation of 1) to avoid multicollinearity between treatment when we tested interaction or relation with developmental-related traits but also to avoid large differences in scale in the model. We checked for multicollinearity using the variance inflation factor (VIF, non-multicollinearity when below 2. "car" package). Clutch identity was used as random effect in all our models.

Each clutch was assigned with the value of average telomere length ratio of the two tadpoles sampled. Because we did not measure oxidative status markers and telomere length on the same individual instead, we pooled individuals of the same family or used average values for telomere length, correlations between the two parameters were at clutch levels.

We used η 2 to measure the magnitude of the main effect (package "effect size").

η 2 : <0.01 considered negligible, between 0.01 and 0.06 considered small, between 0.06 and 0.14 considered medium and >0.14 considered large effect.

First to assess whether our AMPA treatment had an effect on oxidative status markers we fitted six linear mixed-effect models (LMERs, package lme4 , [START_REF] Bates | Fitting Linear Mixed-Effects Models Using lme4[END_REF] with "Treatment" as fixed effect and "SOD", "GPx", "CAT", thiols", "SOD/GPX" "SOD/(GPx+CAT)" as dependent variables. Four individuals displayed disproportionally elevated CAT (tested with "rosnertest" and "EnvStats" packages) and were excluded from the analyses. We calculated 95% confidence intervals, p value and coefficient for each levels of fixed effect factor from 1000 bootstraps iterations ("parameters", packages, Appendix A, table S1).). When p<0.05, the magnitude of the main effect was measured using η 2 (package "effectsize").

Second, to test whether developmental-related traits differed between AMPA treatment groups, we fitted three linear mixed-effect models ("lmer" function, "lme4" package [START_REF] Bates | Fitting Linear Mixed-Effects Models Using lme4[END_REF]) with "Treatment" as fixed effect and "embryonic development duration", "total length" and "telomere length" as dependent variables. Because hatching success followed binomial distribution, we fitted one generalized linear mixed models ("glmer", "lme4" package [START_REF] Bates | Fitting Linear Mixed-Effects Models Using lme4[END_REF] with "hatching success" as dependent variable.

Finally, to test if there was a relation between developmental-related traits and oxidative status and if this relation differed between treatment, we fitted LMER models with development-related trait ("embryonic development duration", "total length" and "telomere length") as dependent variables and z_scored oxidative status markers ("SOD", "GPx", "CAT", thiols", "SOD/GPX" "SOD/(GPx+CAT)") as fixed effect in interaction with "Treatment". We ran GLMER for binomial data ("hatching success"). To compare the trends and regression slopes we used "emtrends" function in "emmeans" package. Due to high correlation between oxidative status markers (SOD, GPx, CAT) and ratios, we performed statistical models separately (Appendix A, table S2).

Results

All oxidative status and developmental-related traits at each AMPA concentrations are summarized in Table 2.

Effect of AMPA treatment on oxidative status markers

CAT, SOD and GPx were not influenced by AMPA treatment (respectively: F3,25=2.05, p=0.125; F3,27=0.92 p=0.439; F3,27=0.24, p=0.869, Table 3). When including four outlier values for CAT, AMPA did not influence enzymatic activities (F3,27=0.70, p=0.560). AMPA treatment influenced concentration of thiols (F3,26=3.66, p=0.021, η 2 =0.21, Table 3). Hatchlings from low AMPA group displayed higher thiols than all other treatment groups (all p<0.028, Tukey HSD, Fig. 1).

AMPA treatment significantly influenced SOD/(GPx+CAT) ratio (F3,25=3.93, p=0.017, η 2 =0.27, Table 3). When including four outlier values for CAT, we found a marginal effect of AMPA in SOD/(GPx+CAT) ratio (F3,27=2.92, p=0.052). Ratios at Low and Medium concentrations were significantly lower compared to High concentration (respectively: -0.16[-0.30, -0.01], p=0.029; -0.21[-0.35, -0.06], p=0.006, Fig. 1). In Medium treatment, SOD/(GPx+CAT) ratio decreased compared to Control treatment (-0.13[-0.24, -0.02], p=0.026).

Finally, treatment did not affect SOD/GPx ratio (F3,36=2.02, p=0.127, Table 3) nor telomere length (F3,36=0.51, p=0.673)

Effect of AMPA treatment on developmental-related traits

Hatching success varied among treatments (χ 2 3=25.168, p<0.001); Low and Medium treatment yielded lower hatching success than Control or High treatment of AMPA (all p>0.003).

AMPA treatment did not influence embryonic development duration (F3,67=1.15, p=0.33).

However, total length of hatchlings differed between treatments (F3,67=5.57, p=0.002).

Hatchlings were larger at Low treatment compared to individuals from Control and High treatments (respectively, p=0.004, p=0.002).

Relationships among oxidative status markers and developmental-related traits

We did not find any relationships among any oxidative status marker and telomere length (all interaction p>0.26). Total length of hatchlings was not related to any oxidative status markers (p>0.25). Relation between total length of hatchlings and oxidative status markers in interaction with treatment was never significant (all p>0.201). We did not find any relation between oxidative status markers and embryonic development duration across or among groups (all p>0.186) Finally, although we found a significant interaction between SOD/GPx and SOD/(CAT+GPx) ratios, and hatching success, after close examination of the plots and removals of outliers, such interactions were no longer supported.

Discussion

Our study shows that the primary metabolite of glyphosate -AMPA -alters oxidative status, but not telomere length, during embryonic development in an amphibian species that occurs in agricultural landscapes. Interestingly, these physiological effects occurred at environmentally relevant-concentrations representative of actual concentrations of this contaminant in the wild. Furthermore, our results suggest non-monotonicity for some markers of oxidative status such as thiols and SOD/(CAT+GPx) ratio as shown for developmental-related traits in a related study [START_REF] Cheron | Aminomethylphosphonic acid alters amphibian embryonic development at environmental concentrations[END_REF]. SOD/(CAT+GPx) ratio is likely driven by CAT levels considering the trend for CAT to increase at Low and Medium level while GPx remained steady.

Exposure to low concentrations of AMPA influenced the levels of thiols but did not influence CAT levels, although both markers showed similar trends. For these markers, levels were higher at the lowest or the intermediate AMPA concentrations tested which suggest low concentration responses and non-monotonicity [START_REF] Hill | Nonmonotonic Dose-Response Curves Occur in Dose Ranges That Are Relevant to Regulatory Decision-Making[END_REF][START_REF] Shi | Progress in the studies on hormesis of low-dose pollutants[END_REF][START_REF] Vandenberg | Hormones and Endocrine-Disrupting Chemicals: Low-Dose Effects and Nonmonotonic Dose Responses[END_REF]. Interestingly, although we did not detect any influence of AMPA concentrations on SOD and GPx, we found that relevant ratios of these antioxidant enzymes (i.e., SOD/(GPx+CAT)) were influenced by AMPA treatments. Moreover, since both GPx and CAT convert H2O2 into H2O, an unbalanced ratio of superoxide dismutase on GPx and CAT activities (SOD/(GPx+CAT) can lead to a proliferation of deleterious pro-oxidants [START_REF] De Haan | Cu/Zn superoxide dismutase mRNA and enzyme activity, and susceptibility to lipid peroxidation, increases with aging in murine brains[END_REF]. SOD/(GPx+CAT) ratios showed low concentration responses and nonmonotonicity, with the lowest values obtained for the lower and intermediate AMPA concentrations. These results suggest a possible low-concentration stimulatory effect of AMPA with upregulation of both CAT and thiols. One important feature of the regulation of oxidative status involves, first the dismutation by superoxide dismutase (SOD) of . O2 -into H2O2 and then the conversion of H2O2 to H2O catalysed by CAT and/or GPx. However, when H2O2 is produced at high concentrations, CAT is more efficient in the conversion of H2O2 in H2O than GPx, which might explain why we observed upregulation of CAT but not of GPx [START_REF] Halliwell | Free Radicals in Biology and Medicine[END_REF]. It is important to stress that antioxidants play an important role in protecting early development and, accordingly, their levels are known to increase during embryonic development in amphibians (Dandapat et al., 2003;Anguiano et al., 2001). Moreover, during early development, CAT levels are high while GPx levels are low which show the important role of CAT in removing toxic hydrogen peroxide [START_REF] Aceto | Developmental Aspects of Detoxifying Enzymes in Fish ( Salmo Iridaeus )[END_REF], and GPx activities tend to gradually increase during growth [START_REF] Pašková | Teratogenicity and Embryotoxicity in Aquatic Organisms After Pesticide Exposure and the Role of Oxidative Stress[END_REF].

To our knowledge, only few studies investigated relevant concentration of AMPA on biochemical parameters [START_REF] Guilherme | DNA and chromosomal damage induced in fish (Anguilla anguilla L.) by aminomethylphosphonic acid (AMPA)-the major environmental breakdown product of glyphosate[END_REF][START_REF] Matozzo | Ecotoxicological hazard of a mixture of glyphosate and aminomethylphosphonic acid to the mussel Mytilus galloprovincialis (Lamarck 1819)[END_REF][START_REF] Matozzo | Effects of aminomethylphosphonic acid, the main breakdown product of glyphosate, on cellular and biochemical parameters of the mussel Mytilus galloprovincialis[END_REF] and most studies have investigated its parent compound, glyphosate ( [START_REF] De Brito Rodrigues | Impact of the glyphosate-based commercial herbicide, its components and its metabolite AMPA on non-target aquatic organisms[END_REF]. [START_REF] Guilherme | DNA and chromosomal damage induced in fish (Anguilla anguilla L.) by aminomethylphosphonic acid (AMPA)-the major environmental breakdown product of glyphosate[END_REF] found DNA and chromosomal damage when Anguilla Anguilla were exposed for 1 or 3 days at sublethal concentration (11.8 and 23.6 µg.l -1 ). In marine invertebrates, AMPA exposure altered oxidative status and enzymatic activities in Mytilus galloprovincialis when exposed for 7, 14 or 21 days to 1, 10 and 100 µg.l -1 [START_REF] Matozzo | Ecotoxicological hazard of a mixture of glyphosate and aminomethylphosphonic acid to the mussel Mytilus galloprovincialis (Lamarck 1819)[END_REF][START_REF] Matozzo | Effects of aminomethylphosphonic acid, the main breakdown product of glyphosate, on cellular and biochemical parameters of the mussel Mytilus galloprovincialis[END_REF]. In tadpole of another amphibian species (Pelobates cultripes), exposure to sublethal concentration of glyphosate induced altered antioxidant activity by decreasing glutathione reductase activity but did not affect SOD nor CAT levels [START_REF] Burraco | Physiological Stress Responses in Amphibian Larvae to Multiple Stressors Reveal Marked Anthropogenic Effects even below Lethal Levels[END_REF].

In Oncorhynchus mykiss , glyphosate-based formulations altered energy metabolism (Le Du-Carrée et al., 2021) when exposed to 1 µg.l -1 and in Anguilla Anguilla increased CAT levels in gills when exposed for 1 or 3 days at 58 and 116 µg.l -1 [START_REF] Guilherme | DNA damage in fish (Anguilla anguilla) exposed to a glyphosate-based herbicide -Elucidation of organ-specificity and the role of oxidative stress[END_REF]. Focusing on glyphosate's primary metabolite, our study reinforces the potential of this widely used herbicide and of its metabolite AMPA (the main contaminants detected in surface waters worldwide, [START_REF] Battaglin | Glyphosate and Its Degradation Product AMPA Occur Frequently and Widely in U.S. Soils, Surface Water, Groundwater, and Precipitation[END_REF][START_REF] Medalie | Influence of land use and region on glyphosate and aminomethylphosphonic acid in streams in the USA[END_REF][START_REF] Okada | Glyphosate and aminomethylphosphonic acid (AMPA) are commonly found in urban streams and wetlands of Melbourne, Australia[END_REF][START_REF] Silva | Distribution of glyphosate and aminomethylphosphonic acid (AMPA) in agricultural topsoils of the European Union[END_REF] to alter oxidative status in wildlife.

Interestingly, given the lower hatching success in Low and Medium groups, selective mortality of embryos susceptible to oxidative stress may have occurred in these two groups prior to hatching. According to this hypothesis, surviving embryos may be the individuals that are more resistant to both AMPA and AMPA-mediated oxidative stress. Such hypothesis dovetails with the lack of relationship between developmental traits and oxidative status markers. Indeed, embryonic development can be constrained by oxidative stress thought oxidation of the thiol group of proteins [START_REF] Goto | AminoMethylPhosphonic acid (AMPA) in natural waters: Its sources, behavior and environmental fate[END_REF], and accelerated development might enhance metabolism and induced cellular oxidative stress (Burraco et al., 2020b;[START_REF] Smith | Meta-analysis indicates that oxidative stress is both a constraint on and a cost of growth[END_REF]. This process might also explain why we did not detect any correlation between telomere length and oxidative status parameters at clutch level. Telomere dynamics are influenced by environmental conditions and telomere length generally shortens under environmental stress [START_REF] Chatelain | The association between stressors and telomeres in non-human vertebrates: a meta-analysis[END_REF]. We expected that a reduction in telomere length because prior work showed that exposure to environmentally relevant concentrations of AMPA may damage DNA ( [START_REF] De Brito Rodrigues | Impact of the glyphosate-based commercial herbicide, its components and its metabolite AMPA on non-target aquatic organisms[END_REF][START_REF] Guilherme | DNA and chromosomal damage induced in fish (Anguilla anguilla L.) by aminomethylphosphonic acid (AMPA)-the major environmental breakdown product of glyphosate[END_REF]. In our context, telomere length may be influenced by increased production of pro-oxidants or decreased of antioxidant defences [START_REF] Chatelain | The association between stressors and telomeres in non-human vertebrates: a meta-analysis[END_REF][START_REF] Reichert | Does oxidative stress shorten telomeres in vivo ? A review[END_REF]. However, if surviving embryos were more able to cope with oxidative stress, damages to their DNA may be reduced [START_REF] Coluzzi | Oxidative Stress Induces Persistent Telomeric DNA Damage Responsible for Nuclear Morphology Change in Mammalian Cells[END_REF]. Since the proliferation of H2O2 induces cellular senescence [START_REF] De Haan | Elevation in the Ratio of Cu/Zn-Superoxide Dismutase to Glutathione Peroxidase Activity Induces Features of Cellular Senescence and This Effect Is Mediated by Hydrogen Peroxide[END_REF] or telomere attrition in human primary fibroblasts [START_REF] Coluzzi | Oxidative Stress Induces Persistent Telomeric DNA Damage Responsible for Nuclear Morphology Change in Mammalian Cells[END_REF], selective mortality might select coping individuals with preserved telomere length. Early development constitutes a phase where oxidative stress levels are high because of the presumed high metabolic needs. However, our results on AMPA exposure suggest an alteration of the trade-off between oxidative stress and developmental traits, thus the potential for AMPA to affect expression of life-history traits.

Conclusions

Our results provide rare evidence for physiological effects of environmentally-relevant AMPA concentrations on amphibian embryos [START_REF] Cheron | Aminomethylphosphonic acid alters amphibian embryonic development at environmental concentrations[END_REF]. Interestingly, the integration of markers of oxidative status in our investigations allows to suggest potential physiological mechanisms responsible for the effects on developmental related traits detected in a previous study [START_REF] Cheron | Aminomethylphosphonic acid alters amphibian embryonic development at environmental concentrations[END_REF]. Our results also strengthen the idea that AMPA may also induce non-monotonic effects. These effects might indicate stimulation of the antioxidant machinery to buffer any pro-oxidant effect of AMPA, but also selected embryos that were more resistant to oxidative stress, two hypotheses which are not mutually exclusive. Further experiments using non-lethal concentration of AMPA are required in order to disentangle whether AMPA acts directly on the antioxidant machinery of tadpoles or if surviving individuals have higher ability to cope with oxidative stress. If our hypothesis on selective mortality holds true, it is possible that the long-term use of glyphosate may have acted as a potent selective pressure, thereby explaining why our study species can persist in agricultural landscapes [START_REF] Guillot | Landscape influences the morphology of male common toads ( Bufo bufo )[END_REF], but see Renoirt et al. 2021).

In light of our findings, we urge for regulatory toxicology to investigate effects of parent compounds but, more importantly, also of their metabolites at environmentally relevant concentrations on non-model species in order to comprehensively assess impacts of anthropogenic contaminants on the environment. Table 1. Concentrations (µg.l -1 ) measured in the experimental tanks (Mean± SD)

Tables

Acknowledgements

We thank Michael Bank and the anonymous reviewers for their helpful comments that improved previous versions of this manuscript. We thank the CNRS, the Agence de l'Eau Loire-Bretagne, the Agence de l'Eau Adour-Garonne, the Région Nouvelle-Aquitaine (Multistress 2017-1R20110, Aquastress 2018-1R20214, Amphitox 2019-1R20216), the ANSES (BiodiTox project # 2019/1/031) and the CPER Econat for funding our work.

Funding

Funding was provided by the CNRS, the Agence de l'Eau Loire-Bretagne, the Agence de l'Eau Adour-Garonne, the Région Nouvelle-Aquitaine (Multistress 2017-1R20110, Aquastress 2018-1R20214, Amphitox 2019-1R20216), the ANSES (BiodiTox project # 2019/1/031) and the CPER Econat.

Abbreviations

AMPA, aminomethylphosphonic acid; SOD, superoxide dismutase; GPx, glutathione peroxidase; CAT, catalase; SODGPx, ratio between SOD and GPx; SODCAT, ratio between SOD and CAT + GPx. GS: Gosner stage

Ethics approval

All applicable institutional and/or national guidelines for the care and use of animals were followed. This work was approved by the French authorities (COMETHEA ethic committee and Ministère de L'Enseignement Superieur, de la Recherche et de L'innovation) under permit #APAFIS#13477-2018032614077834 v7

Competing interests

We declare we have no competing interests.

Caption to figures

Figure 1: A/ Superoxide dismutase activity (units SOD/mg proteins), B/ Glutathione peroxidase activity (units GPx/mg proteins) C/ Thiols (mmol -1 of -SH groups/mg proteins,) D/ CAT activity (units CAT/mg proteins) E/ SOD/(GPx+CAT) ratio F/ SOD/GPx ratio relative to AMPA treatments in Bufo spinosus tadpoles. Data represent mean ± SE. Letters indicate significant differences.