Nadia Bouain

Gabriel Krouk

Benoît Lacombe

Hatem Rouached

published or not. The documents may come   L'archive ouverte pluridisciplinaire


Sara Maqrot

Simon De Givry

Marc Tchamitchian

Gauthier Quesnel
email: gauthier.quesnel@inrae.fr

Learning parameters of the Wedelin heuristic with application to crew and bus driver scheduling $

Keywords: integer programming, heuristics, Lagrangian relaxation, sensitivity analysis, parameter tuning, crew and bus driver scheduling, set partitioning problem 2010 MSC: 90C10, 90C59

Heuristics are important techniques designed to find quickly good feasible solutions for hard integer programs. Most heuristics depend on a solution of the relaxed linear program. Another approach based on Lagrangian relaxation offers a number of advantages over linear programming, namely it is extremely fast for solving large problems. Wedelin heuristic is such a Lagrangian based heuristic, initially developed to solve airline crew scheduling problems. The performance of this method depends crucially on the choice of its numerous parameters. To adjust them and learn which ones have important influence on whether a solution is found and its quality, we propose to conduct a sensitivity analysis followed by an automatic tuning of the most influential parameters.

We have implemented a C++ open-source solver called baryonyx which is a parallel version of a (generalized) Wedelin heuristic. We used the Morris method to find useful continuous parameters. Once found, we fixed other parameters and let a genetic optimization algorithm using derivatives adjust the useful ones in order to get the best solutions for a given time limit and training instance set. Our experimental results done mostly on crew and bus driver scheduling benchmarks tackled as set partitioning problems show the significant improvements obtained by tuning the parameters and the good performances of our approach compared to state-of-the-art exact and inexact integer programming solvers.

Introduction

Heuristics are important techniques designed to find quickly good feasible solutions for hard integer programs. Most heuristics depend on a solution of the relaxed linear program. Another approach based on Lagrangian relaxation offers a number of important advantages over linear programming [START_REF] Marshall | The Lagrangian relaxation method for solving integer programming problems[END_REF], namely it is extremely fast for solving large problems.

Wedelin heuristic [START_REF] Wedelin | An algorithm for large scale 0/1 integer programming with application to airline crew scheduling[END_REF] is such a Lagrangian based heuristic. Wedelin [START_REF] Wedelin | The design of a 0/1 integer optimizer and its application in the Carmen system[END_REF] described the basic principle of the algorithm and its application in the Carmen system for scheduling crew rotations for airlines. Alefragis et al. [START_REF] Alefragis | Parallel integer optimization for crew scheduling[END_REF] presented a scalable parallelization of the original algorithm used in the Carmen system. Grohe and Wedelin [START_REF] Grohe | Cost Propagation-Numerical Propagation for Optimization Problems[END_REF] and [START_REF] Wedelin | Revisiting the in-the-middle algorithm and heuristic for integer programming and the max-sum problem[END_REF] introduced a similar algorithm for the max-sum problem. Ernst et al. [START_REF] Ernst | A new Lagrangian heuristic for the task allocation problem[END_REF] described a variation of the Wedelin algorithm and applied it to the staff planning problem. Bastert et al. [START_REF] Bastert | A generalized Wedelin heuristic for integer programming[END_REF] presented many extensions and generalizations of Wedelin algorithm with various improvements. They evaluated the performance of their variant on a set of instances from different sources, where the results were favorable compared to FICO Xpress, a commercial optimization software. Starting from [START_REF] Bastert | A generalized Wedelin heuristic for integer programming[END_REF], we propose an extension performing multiple runs in parallel.

A major difficulty in the use of optimization methods is the parameter setting. It is important for each problem to find a set of parameter values that leads to optimal performances. Choosing the best values manually requires a lot of experimentation. There are several methods for automatic configuration of parameter values [START_REF] Ágoston | Parameter control in evolutionary algorithms[END_REF][START_REF] Hamadi | Autonomous search[END_REF]. Eiben et al. [START_REF] Ágoston | Parameter control in evolutionary algorithms[END_REF] classify methods into two categories according to the manner of use, before running the optimization algorithm (parameters tuning) and during its execution (parameter control).

The control methods include self-adjust the parameter values for the resolution and dynamically control these values to improve the solution search. These techniques are widely applied to self-adjust the parameters of evolutionary algorithms [START_REF] Hamadi | Autonomous search[END_REF], such as crossover rate [START_REF] Yang | Adaptive crossover in genetic algorithms using statistics mechanism[END_REF], mutation rate [START_REF] Smith | Self adaptation of mutation rates in a steady state genetic algorithm[END_REF] and population size [START_REF] Endre Eiben | Evolutionary algorithms with on-the-fly population size adjustment[END_REF] for genetic algorithms. These methods can be subdivided into two branches: deterministic and (self-)adaptive [START_REF] Hinterding | Adaptation in evolutionary computation: A survey[END_REF]. Deterministic methods are based on deterministic rules, which do not change during the execution of the algorithm. Their goal is to calculate approximate values of parameters and adjust them according to the problem [START_REF] Barbulescu | Dynamic representations and escaping local optima: Improving genetic algorithms and local search[END_REF]. Adaptive methods use information on the current state of the search to change parameter values. Adaptation is effected by changing the objective function, by increasing or decreasing the penalty coefficients for violations of constraints, from one generation to another. This prevents poor settings to conduct future generations [START_REF] Harik | A parameter-less genetic algorithm[END_REF]. These control methods are developed in an automated framework for setting the parameters of a specific problem (one instance). If the goal is to solve different instances, these methods can be costly in terms of computing time, given the number of parameters and instances.

The principle of parameter tuning methods is simpler. Parameters do not change values during search. They are adjusted before the execution of the algorithm, and remain unchanged after. The settings obtained by learning on a subset of instances is used to solve all instances of a problem. Different parameter tuning strategies exist in the literature and depend on the type of parameters: discrete, continuous or categorical. Two famous examples are CALIBRA [START_REF] Adenso | Fine-tuning of algorithms using fractional experimental designs and local search[END_REF] and ParamILS [START_REF] Hutter | Automatic algorithm configuration based on local search[END_REF] parameter tuning methods for the discrete case. They explore the parameter configuration search space using (iterated) local search. CALIBRA limits the number of parameter configuration evaluations using partial statistical designs. ParamILS saves computation time by doing partial evaluations of the training set.

A sensitivity analysis can facilitate parameter tuning by focusing on important parameters for a problem. The goal is to reduce the time and effort in resolving sensitive parameters. In the literature, the work on the use of both techniques (sensitivity analysis to reduce the parameter search space and automatically adjusting their values) is not much discussed.

Kim et al. [START_REF] Kim | Sensitivity analysis and parameter tuning scheme for global sea-ice modeling[END_REF] perform a sensitivity analysis on the dynamic parameters of sea ice model, using automatic differentiation. Information on the gradient provided by the latter are used in a parameter optimization algorithm based on quasi-Newton method [START_REF] Zhu | Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization[END_REF]. Teodoro et al. [START_REF] Teodoro | Algorithm sensitivity analysis and parameter tuning for tissue image segmentation pipelines[END_REF] combine sensitivity analysis and automatic parameter tuning for an image segmentation problem. This approach permits to identify the least influential parameters and reduce the parameter configuration search space (100 points instead of trillion points).

In these previous works, parameter tuning is performed for a single instance. In this paper we propose a protocol to generate a universal set of parameter values from a subset of instances of a given problem to be solved by our multi-run Wedelin algorithm. The basic idea is to i) select a subset of instances, ii) investigate the sensitivity of parameters for fixing the values of the non-sensitive parameters to their default values, iii) automatically adjust the values of the significant parameters by black-box optimization, and iv) apply the learned parameter setting on all instances of the problem.

In Section 2, we describe the Wedelin heuristic and our multi-run extension for parallelism. In Section 3, we present our protocol for parameter tuning. In Section 4, we give experimental results and conclude.

A Parallel Version of Wedelin Heuristic

We are interested in minimizing 0/1 linear programs of the following form:

min cx s.t. Ax = b x ∈ {0, 1} n
where c ∈ R n is a vector of n costs, b ∈ N m a column vector of constant terms, and A ∈ {0, 1} m×n a coefficient matrix for the constraints. In the following, constraint index i refers to the ith row in equation Ax = b, and index j to the jth column in A associated to variable x j . Let J = {1, . . . , n} be the set of column indices, and J(i) = {j : a ij = 1} ⊆ J the subset of column indices occurring in constraint i. Similarly, let I = {1, . . . , m} be the set of row indices, and I(j) = {i : a ij = 1} ⊆ I the subset of constraint indices involving x j .

In this paper we adopt an approximate method, called Wedelin algorithm or Inthe-middle algorithm. This method is designed for solving linear problems on Boolean variables [START_REF] Wedelin | An algorithm for large scale 0/1 integer programming with application to airline crew scheduling[END_REF] and uses the Lagrangian relaxation.

The Lagrangian relaxation

The Lagrangian relaxation consists in moving into the objective function all the constraints. These constraints are built into the objective function as linear combinations where the coefficients are Lagrangian multipliers. They penalize the objective function if one of the integrated constraint is violated. This relaxation has the advantage over linear relaxation to directly provide integer solutions and to tackle large size problems.

The Wedelin Heuristic

In this paper, we are particularly interested in a specific heuristic based on the Lagrangian relaxation known as Wedelin heuristic. This method solves linear programming problems with Boolean variables with a specific mathematical form such as the set partitioning problem. This heuristic tries to directly solve the linear problem:

min x cx -π(Ax -b) (1)
s.c.

x ∈ {0, 1} n

Where π = {π 1 , π 2 , ..., π m } ∈ R m represent the Lagrangian multipliers or dual variables. The resolution is trivial when π i elements in π are fixed to π:

πb+ min x (c -πA)x (2) s.c. x ∈ {0, 1} n
The solution xj is constructed with the reduced cost sign with (c j -πa j ) and a j the associated column in the coefficient matrix A.

xj =    1 if (c j -πa j ) < 0, 0/1 if (c j -πa j ) = 0, 0 if (c j -πa j ) > 0. (3)
If the reduced cost c = (c -πA) is non-zero, the solution is unique. However, if one or a few elements of c are zero, it will be difficult to find a feasible solution for the not relaxed problem. The goal of this algorithm is to find a π where all elements are non-zero and where the single solution of the relaxed problem is realizable for the primal problem.

We note that the change of the value of a single component πi modifies the values (and signs) of the reduced costs c, and therefore affects the x solution.

The main idea of the Wedelin algorithm (see Algorithm 1) is to consider iteratively a single constraint i, updating the element πi such that the x vector satisfies the i constraint of the primal problem. There is always a selection range of πi . The algorithm chooses the value of πi in the middle of this interval. This corresponds to the dual search for descent by coordinates.

Algorithm 1 represents the core of the solver. It takes as input a definition of the problem: the matrix A and the vectors b and c. It expects an output vector of Boolean, x, a solution of the problem. The algorithm starts at line 1 with the construction of an initial x vector which permits to build the list of violated constraints R. It uses the Bastert et al. [START_REF] Bastert | A generalized Wedelin heuristic for integer programming[END_REF] proposal to penalize variables with positive costs. The main loop begins at line 2 and run until the loop limit p and/or the κ max is reached. For each violated constraints (line 3), it (i) decreases the preferences p, (ii) produces the reduced cost vector r and sort it according to the reduced costs, (iii) updates the Lagrangian multipliers at line 4, (iv) affects the x vector and the preference matrix p. Finally, at line 6, we update the list of violating constraints and exit with the solution found or adjust the κ adjustment and the iteration l.

Input : A ∈ {0, 1} m×n , b ∈ N m , c ∈ R n Output:
x ∈ {0, 1} n with Ax = b and cx small or message no solution for j ∈ {1, . . . , n} do // Build an initial variable assignment

if c j ≤ 0 then xj ← 1 else xj ← 0 end Let π ∈ R m , π ← 0, p ∈ R m×n , p ← 0, l ← 1, κ ← κ min , R ← {i : n j=1 a ij xj = b i } // List of violated constraints while l ≤ ρ and κ ≤ κ max do for i ∈ R do // Update every violated constraint for j ∈ J(i) = {j : a ij = 0} do pij ← θ × pij // History exp. decay for j ∈ J(i) do // Build reduced costs r j ← c j -i∈I(j) a ij πi -i∈I(j) pij end r ϕ[1] < r ϕ[2] < ... < r ϕ[|J(i)|]
// Sort variables by increasing red. costs

πi ← πi + 1 2 (r ϕ[bi+1] + r ϕ[bi] ) // Update Lagrangian multipliers if l ≤ ω then ∆ ← 0 // Do not perturb reduced costs during warmup else ∆ ← κ 1-κ (r ϕ[bi+1] -r ϕ[bi] ) + δ for j ∈ {1, . . . , b i } do // Update variables and preferences positively piϕ[j] ← piϕ[j] + ∆ ; xϕ[j] ← 1 end for j ∈ {b i + 1, . . . , |J(i)|} do // or negatively piϕ[j] ← piϕ[j] -∆ ; xϕ[j] ← 0 end end R ← {i : n j=1 a ij xj = b i } // Update the list of violated constraints if R = ∅ then return x // If empty, exit with solution found κ ← κ + κ step ( |R| m ) α // Adaptive adjustment of step size l ← l + 1 // Next iteration
end return no solution found Algorithm 1: Wedelin algorithm with local preferences and adaptive step size.

The Parallel Solver Baryonyx

In this paper, we introduce a new integer and Boolean linear programming solver called baryonyx. This solver is largely based on the algorithms provided in Bastert et al. [START_REF] Bastert | A generalized Wedelin heuristic for integer programming[END_REF]. We have introduced several modifications to the proposed algorithms to (i) allow the reuse of previous solutions found, (ii) diversify the search, and also (iii) exploit todays symmetric multiprocessor CPUs.

Baryonyx accepts two modes: solver and optimizer. In the solver mode, it runs once trying to satisfy all the constraints (the exact implementation of Algorithm 1). In the optimizer mode, it runs in parallel according to the number of processors, and tries to satisfy all the constraints and to optimize the solution at each run, reporting the best solution found for all runs when it reaches its time limit.

The Wedelin heuristic depends only on the δ parameter (line 5 in Algorithm 1) to diversify the search. To improve diversification, we develop several random processes (see the technical documentation of baryonyx). The most important process of diversification is the initialization mechanism of x. Indeed, x is used to determine violated constraints before any computation. Default, a deterministic initialization is proposed by Bastert et al. [START_REF] Bastert | A generalized Wedelin heuristic for integer programming[END_REF] where xj equals 1 if c j ≤ 0 otherwise 0. We propose to extend this part by combining a random process and several different algorithms. We use the Bernoulli's law and its parameter p ∈ [0, 1] to provide random Boolean.

random Each xi are initialized with the Bernoulli law.

bastert For each xi , either the variable is initialized by the cost variable c i (See section 1 in Algorithm 1) or by a random Boolean.

best For each xi , either the variable is initialized by the best solution found previously xi or by random Boolean.

Since Wedelin's algorithm runs several times in baryonyx, we changed the initialization of x for each iteration. The default baryonyx choice is the best-cycle policy. It starts with the bastert policy. If is fails to find a solution, it restarts with the random policy, else it uses the best solution found as initial solution and it tries to improve this solution with the best policy three times. If it fails it switches to the first step of the algorithm (bastert policy) otherwise it keeps the best policy. Figure 1 shows the diagram of the best-cycle policy.

To further increase the diversification, we exploit the symmetric multiprocessor CPUs. Each core shares matrices A, vectors b and c. The vector x is local for each process. Each core have its own pseudo-random number generator initialized with different random seed.

Learning Continuous Parameter Values

In this section, we present our protocol for tuning continuous parameters of Wedelin heuristic automatically. Algorithm 1 has eleven numerical parameters as described in Table 1. Parameters τ, Γ, ρ, ω are integers whereas the others are continuous. We assume a fixed time limit τ and a fixed number of cores Γ, and treat the other integer parameters as continuous values ] Bernoulli distribution with success probability p using a rounding function. So, we have nine continuous parameters to tune. The Wedelin heuristic is viewed as a complex function with K = 9 input parameters (x 1 , . . . , x K ) or factors and a single output value y(x 1 , . . . , x K ) which corresponds to the quality of the parameter configuration over a set of training instances. For each benchmark category, we select an evenly distributed subset (20%) of instances to be part of the training set.

Quality of a parameter configuration

The quality of a parameter configuration (x 1 , . . . , x K ) is computed as follow. First, we execute our parallel baryonyx solver on every training instance using a modified Wedelin-Good configuration1 . Each execution uses Γ cores2 in parallel during 2τ time limit. We collect the best l init

y(x 1 , . . . , x K ) = 1 - 1 N N e=1 l e -l init e 10u init e -l init e
If no solution is found for a particular instance and parameter configuration then we assume l e = 10u init e . By multiplying u init e by 10, we assume not finding a solution is at least ten times worse than finding a good solution close to the best initial solution. When y > 1, we have found a better configuration than the modified WedelinGood configuration.

Selection of important parameters by sensitivity analysis

We used the Morris method [START_REF] Morris | Factorial sampling plans for preliminary computational experiments[END_REF][START_REF] Saltelli | Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models[END_REF] to determine which inputs have important effects on the output. The goal is to discover which parameters are important and cannot be ignored. These parameters will be fine tuned later on (see Section 3.3). For that a factorial sampling plan is built from individually randomized one-factor-at-a-time designs. The Morris method randomly selects L initial configurations within a regular K-dimensional d-level grid. Each parameter is discretized into d levels including its bounds (given in Table 1). Starting from each initial configuration, a configuration trace is constructed by changing the value of one parameter at a time until all the parameters have been modified. The step value is usually chosen as ∆ = (d -1)/2. See an example in Figure 2(left). Each trace corresponds to K +1 configuration evaluations, resulting in a total of L×(K +1) evaluations.The Morris method performs a limited amount of configuration evaluations which is a linear function of the number of parameters K.

The Morris method provides qualitative sensitivity measures allowing to rank the input factors in order of importance, but not to quantify by how much one given factor is more important than another [START_REF] Saltelli | Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models[END_REF]. In Figure 2

δy l k = y(x 1 , . . . , x k-1 , x k + ∆, x k+1 , . . . , x K ) -y(x 1 , . . . , x k , . . . , x K ) ∆
When µ * k is large but not σ k , parameter x k has an important overall influence on the output. If both measures are large, it corresponds to a non-linear effect on the output or an input involved in interaction with other factors [START_REF] Saltelli | Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models[END_REF]. This is the case for parameter κ step in Figure 2(right). We use µ * k to rank the input factors. In our experiments, we have K = 9, L = 50, d = 10. We selected the four most important factors, κ min , κ step , δ, θ, which was the same set of parameters found by the Morris method in all our benchmark categories.

Optimization of Selected Parameter Values

genoud (GENetic Optimization Using Derivatives) [START_REF] Mebane | R version of genetic optimization using derivatives[END_REF] is a black-box optimization method for solving nonlinear, nonsmooth, and even discontinuous functions. It combines a genetic algorithm with a quasi-Newton method. The quasi-Newton method is the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method [START_REF] Gill | Practical Optimization[END_REF] using a built-in numerical derivative. It helps to quickly find a local optimum when the current parameter configuration is "in a smooth neighborhood of the local optimum point" [START_REF] Mebane | R version of genetic optimization using derivatives[END_REF]. It may also prevent to find the global optimum if used too early or too aggressively. In practice, a number of initial warmup generations are explored before BFGS is applied on the best individual in the current population.

genoud starts with a random population including the modified WedelinGood configuration. Then, next generations are built from 8 genetic algorithm operators (used in equal proportion) dedicated to continuous parameters: Cloning, Uniform Mutation, Boundary Mutation, Non-Uniform Mutation, Polytope Crossover, Simple Crossover, Whole Non-Uniform Mutation, and Heuristic Crossover [START_REF] Mebane | R version of genetic optimization using derivatives[END_REF]. For instance, the Polytope Crossover computes a new parameter configuration that is a convex combination of as many individuals as there are parameters to tune.

In our experiments, a (hard) maximum number of 10 generations is performed with a population size of 100. BFGS is applied on the best individual at each generation after the eighth generation. We observed between 803 to 1157 evaluations of the output value y. It can be more than 10 × 100 due to BFGS evaluations and less because genoud will not evaluate the same configuration twice.

Experimental Results

Baryonyx is a free software (MIT license) in C++17 for solving Boolean and integer linear programming problems. It is provided as a command line program, as a shared library, and has an encapsulation to the R software. Following results were achieved using the version 0.4 of baryonyx built with gcc-7.2.

The baryonyx wrapper for R is called Rbaryonyx. It relies on the rcpp package to facilitate exchanges between the two libraries. Rbaryonyx can read lp files, solves the problem and returns a list of solution(s) and solving information data to easily link with other packages such as R sensitivity [START_REF] Iooss | A review on global sensitivity analysis methods[END_REF] for the Morris method and rgenoud [START_REF] Mebane | R version of genetic optimization using derivatives[END_REF].

All computations were performed on a cluster of 32-physical-core Intel Xeon CPU E5-2683 v4 at 2.1 GHz and 4 GB of RAM per core. In order to speed up the parameter tuning process, we parallelize the evaluation of the training instances, by taking care that the actual number of executions of baryonyx multiplied by Γ is less than the available number of cores. Depending on the size of the training set, each parameter configuration evaluation took between 1 to 1.5 minute. The Morris method took between 9 to 13 hours per benchmark family. The genoud method took from 14 to 28 hours.

In the following, we compare baryonyx against IBM ILOG cplex version 12.8 and LocalSolver3 version 8.0. cplex is a state-of-the-art exact MIP solver. LocalSolver is a mathematical programming local search solver using a simulated annealing based on ejection chain moves specialized for maintaining the feasibility of Boolean constraints and an efficient incremental evaluation using a directed acyclic graph [START_REF] Benoist | Localsolver 1. x: a black-box local-search solver for 0-1 programming. 4OR[END_REF]. Every solver ran in parallel mode using 30 cores. The solving time limit is 60 seconds (except for VCS where it is 1, 800 sec.). cplex and LocalSolver use their default parameter configurations 4 . baryonyx ran with three different static parameter configurations (Supplementary Table .7) plus the one found by rgenoud (Supp. Table .8, where the corresponding generation of the best configuration found is also mentioned).

When available, we also report results from other publications [START_REF] Borndörfer | Aspects of set packing, partitioning, and covering[END_REF][START_REF] Suniel D Curtis | Constructing driver schedules using iterative repair[END_REF][START_REF] Bastert | A generalized Wedelin heuristic for integer programming[END_REF][START_REF] Umetani | Exploiting variable associations to configure efficient local search algorithms in large-scale binary integer programs[END_REF], but we must treat these results with care as they do not correspond to the same computer, nor time limit, and were obtained on a sequential machine.

All the instances have been preprocessed by cplex. A direct translation of 0/1 linear programming lp file format into LocalSolver lsp format has been done.

SPP benchmark

This test set consists in 55 instances of airline crew scheduling problems expressed as set partitioning problems. These problems are obtained from the OR-library [5]5 . SPP instances were provided by four different airline corporations where a subset of these problems was initially solved in [START_REF] Hoffman | Solving airline crew scheduling problems by branch-andcut[END_REF], and further experimented by other exact [START_REF] Borndörfer | Aspects of set packing, partitioning, and covering[END_REF] and local search methods [START_REF] Chu | A genetic algorithm for the set partitioning problem[END_REF][START_REF] Umetani | Exploiting variable associations to configure efficient local search algorithms in large-scale binary integer programs[END_REF].

These instances are easily solved by cplex within the 60-second time limit. The same optimum values were also reported in [START_REF] Borndörfer | Aspects of set packing, partitioning, and covering[END_REF] and [START_REF] Umetani | Exploiting variable associations to configure efficient local search algorithms in large-scale binary integer programs[END_REF] (when available). LocalSolver could not find a solution on four instances, whereas all baryonyx configurations always found a solution. Using default or optimized by genoud configurations were the best options, respectively at 0.06% and 0.01% to the optimum.

Telebus benchmark

Telebus6 [START_REF] Borndörfer | Aspects of set packing, partitioning, and covering[END_REF] is a scheduling problem to program tour vehicles for disabled persons in Berlin. The goal is to provide a one-off service with minimum costs, while respecting a set of constraints, such as vehicle capacity and mandatory breaks. The problem was modeled in two steps. The first step, "clustering", identifies all possible bus circuits that can carry several people at a time. The goal is to select a set of controls with a minimal vehicle travel distance. In the second step, "chaining", the selected commands are chained to generate bus circuits that respect all the constraints. The objective is to minimize the total distance traveled by the vehicles. These two steps represent 28 instances of set partitioning problems divided into two equal-size subfamilies (14 v/clustering, 14 t/chaining) corresponding to different periods in the year 1996. t/chaining instances are more difficult than v/clustering instances and their optimum is mostly unknown (except for t0415, t0420, t0421).

cplex solved all v/clustering and t0415 instances within the time limit, doing slightly better than the original branch-and-cut method of [START_REF] Borndörfer | Aspects of set packing, partitioning, and covering[END_REF] and the 4-flip neighborhood local search algorithm of [START_REF] Umetani | Exploiting variable associations to configure efficient local search algorithms in large-scale binary integer programs[END_REF]. baryonyx and LocalSolver were a few percent below.

baryonyx rgn performed extremely well on t/chaining instances, obtaining the best results among all tested methods except on three instances where cplex and 4-flip local search found better solutions. In particular, for two open instances (t1717/t1722), it improved MIPLIB 2017 best reported solutions (from 184271/114245 to 165881/167523).

Looking at our training set, we observed it does not include any v-instances because the initial modified WedelinGood evaluation (see Section 3.1) failed to produce any solutions.

0.61%

1.88%

0.95%

2.72% (28)

0.44% (28) 2.95% (28) 35.82% (18) 2.90% (28) 3.01% (27)
We also compared our solver using WedelinFast and WedelinGood static configurations with the original version developed by Bastert [START_REF] Bastert | A generalized Wedelin heuristic for integer programming[END_REF]. Our methods were able to find a solution for all the instances whereas the original approach could not for at least six instances. It demonstrates the robustness of doing multiple runs in parallel compared to a single run (using a longer 600-second time limit) as done in [START_REF] Bastert | A generalized Wedelin heuristic for integer programming[END_REF].

CSPLib022 benchmark

CSPLib0227 is a library of 12 bus driver scheduling problems reformulated as set partitioning problems. The problems come from different bus companies: Reading (r1 to r5a), Centre West Ealing area (c1, c1a, c2), the former London Transport (t1 and t2). These problems are relatively small and easy. cplex took less than 7.7 seconds to solve its most difficult instance c1a. baryonyx rgn found all the optima in less than 2.2 seconds (solving time) for the largest instance r3. LocalSolver did not find the optimum for this instance, neither found a solution for two other instances. The iterative repair local search method GENET [START_REF] Suniel D Curtis | Constructing driver schedules using iterative repair[END_REF] got the worst results.

VCS benchmark

VCS instances are randomly generated bus and driver scheduling problems. VCS1200 is a medium-size instance ( 1, 200 constraints, 130, 000 variables), and VCS1600 is a large problem ( 1, 600 constraints, 500, 000 variables) [START_REF] Zaghrouti | Integral simplex using decomposition for the set partitioning problem[END_REF].

With a 1, 800-second time limit, cplex did not find any integral solution for the largest instance whereas baryonyx with its default setting or the one found by genoud produced

0.00%

0.00%

0.00%

0.00%

1.54%

relatively good solutions (5.3% to the optimum for VCS1600). A first solution at 5.7% to the optimum for VCS1600 was found by default baryonyx in 280 seconds (solving time).

Table 5: Computational results (relative distance to best-known solutions for solved instances and in parentheses number of solved instances) on VCS instances [START_REF] Zaghrouti | Integral simplex using decomposition for the set partitioning problem[END_REF].

Instances baryonyx cplex default fast good genoud VCS (2) 9.25% (2) (0) (0) 7.32% (2)
12.68%

In order to find the optimum values, we also ran cplex without any time limit. It took 1 hour for VCS1200 and 4.3 hours for VCS1600 to be solved to optimality by cplex version 12.7.1 on a 4-core Intel CPU i7-4600U at 2.1 GHz.

Nqueens benchmark

Nqueens8 instances represent the weighted n-queen problem where the goal is to place n queens on a chessboard n × n so that none of them attack each other (exactly one queen in each row and each column, at most one queen in each ascending diagonal and each descending diagonal). There are n 2 0/1 variables x i,j . The objective function is to minimize i,j c i,j x i,j with c i,j a random value uniformly sampled from [1, n]. The largest instance among 8 has 1 million variables.

This problem combines set partitioning linear constraints (for rows and columns) and set packing linear constraints (for diagonals). It shows the ability of baryonyx to tackle efficiently a larger class of integer linear programs, with better results than cplex and LocalSolver thanks to parameter tuning by genoud.

Conclusion

baryonyx is a parallel multi-start meta-heuristic which offers good results on large crew and bus driver scheduling problems expressed as set partitioning problems in a relatively short amount of time. It can be used just after preprocessing to provide a good initial upper bound for a complete branch-and-bound solver. Depending on the usage context and time available, when off-line tuning of the parameters is allowed, we show a methodological process to select the important factors and optimize them. Even using a small training set, we could significantly improve the results. This methodology is readily available as supplementary R scripts next to baryonyx source code. Concerning the comparison with the other solvers, it is important to note that we did not try to tune their parameters and better results could be obtained as reported in [START_REF] Hutter | Automated configuration of mixed integer programming solvers[END_REF]. Concerning the set partitioning problems (all our benchmarks except nqueens), more preprocessing rules could be applied.We made some preliminary experiments with set covering problems but the results were not as good as with set partitioning or set packing problems. It remains as future work to evaluate our solver on a larger set of benchmarks.

We plan to exploit the Lagrangian multipliers to give dual bound information which could make our solver exact in some cases, a feature already available in LocalSolver using linear relaxation. Dealing with a quadratic objective function as done in [START_REF] Grohe | Cost Propagation-Numerical Propagation for Optimization Problems[END_REF] is also an important topic we would like to work on and apply to planning problems in agronomy [START_REF] Maqrot | A Mixed Integer Programming Reformulation of the Mixed Fruit-Vegetable Crop Allocation Problem[END_REF]. baryonyx rgn κ min κ step δ θ SPP (VCS) [START_REF] Ágoston | Parameter control in evolutionary algorithms[END_REF] 0.00 1.88 × 10 -4 3.19 × 10 -3 2.85 × 10 -1 telebus [START_REF] Adenso | Fine-tuning of algorithms using fractional experimental designs and local search[END_REF] 1.05 × 10 -1 4.10 × 10 -4 1.13 × 10 -2 3.54 × 10 -1 CSPLib022 [0] 0.00 2.00 × 10 -4 1.00 × 10 -2 6.00 × 10 -1 nqueens [START_REF] Adenso | Fine-tuning of algorithms using fractional experimental designs and local search[END_REF] 3.37 × 10 -2 7.45 × 10 -3 9.04 × 10 -3 7.89 × 10 -1 Umetani [START_REF] Umetani | Exploiting variable associations to configure efficient local search algorithms in large-scale binary integer programs[END_REF] baryonyx Table .12: Computational results (objective value) on VCS instances (using SPP parameter configuration for baryonyx rgn ). "-": no solution found.

Figure 1 :

 1 Figure1: The best-cycle algorithm is designed to both improve diversification and to look for a better solution when a solution is found. Each box represents a run of the complete Wedelin algorithm. The solver stops when time limit is reached.

e

  and worst u init e solution objective values returned by Algorithm 1 for each instance e. If no solution is found we remove this instance from the training set. During the training phase, the normalized quality of a parameter configuration is obtained from the mean over N valid training instances of the relative distance gap to the best initial solutions:

Figure 2 :

 2 Figure 2: (Left) Example of L = 6 traces in the Morris design for K = 2 parameters discretized into d = 8 levels. The step value here is ∆ = 1. (Right) Sensitivity measures of the elementary effects for the CSPLib022 benchmark.

Table . 10 :

 10 Computational results (objective value) on telebus instances (in bold: training set and best solutions, "-": no solution found)

20 Table . 13 :

 2013 Computational results (objective value) on nqueens instances (in bold: training set and best solutions, "-": no solution found).

Table 1 :

 1 baryonyx numerical parameters with their typical domain ranges and descriptions.

	Parameter Range	Description
	τ	[1, +∞[	CPU time limit (in seconds)
	Γ	[1, +∞[	Number of parallel runs
	ρ	[10 2 , 10 5 ]	Number of iterations inside a run
	ω	[0, 100]	Number of warmup iterations before using κ
	κ min	[0, 0.5]	Minimum value for κ approximation
	κ step	[10 -4 , 10 -2 ]	Step value for κ approximation
	κ max	[0.6, 1]	Maximum value for κ approximation
	α	[0, 2]	Adaptive adjustment of κ based on the number
			of violated constraints [4]
	δ	[10 -3 , 10 -1 ]	Random perturbation on reduced costs
	θ	[0, 1]	Strength of history of local preferences p [4]
	p	[10 -4 , 1 -10 -4

Table 2 :

 2 Computational results (relative distance to best-known solutions for solved instances and in parentheses number of solved instances) on SPP instances[START_REF] Beasley | Or-library: distributing test problems by electronic mail[END_REF].

							Local-
	Instances		baryonyx		cplex	Solver	[7]	[31]
		default	fast	good	genoud
	aa	0.23 %	0.80 %	1.19%	0.08%	0.00%	23.02%	0.00%	1.64%
	(6)	(6)	(6)	(6)	(6)	(6)	(3)	(6)	(6)
	us	0.09 %	1.48%	2.60%	0.02%	0.00%	0.00%	0.00%	0.04%
	(4)	(4)	(4)	(4)	(4)	(4)	(3)	(4)	(4)
	nw	0.00%	2.64%	3.96%	0.00%	0.00%	8.53%	0.00%
	(43)	(43)	(43)	(43)	(43)	(43)	(43)	(43)	N/A
	kl	0.69%	0.00%	0.00%	0.00%	0.00%	0.00%	0.00%
	(2)	(2)	(2)	(2)	(2)	(2)	(2)	(2)	N/A
	Mean	0.06%	2.26%	3.41%	0.01%	0.00%	8.54%	0.00%
	(55)	(55)	(55)	(55)	(55)	(55)	(51)	(55)	N/A

Table 3 :

 3 Computational results (relative distance to best-known solutions for solved instances and in parentheses number of solved instances) on telebus instances[START_REF] Borndörfer | Aspects of set packing, partitioning, and covering[END_REF].

							Local-
	Instances		baryonyx		cplex	Solver	[7]	[31]
		default	fast	good	genoud
	v0415-0421	0.14%	0.12%	0.11%	0.30%	0.00%	0.04%	0.00%	0.00%
	(7)	(7)	(7)	(7)	(7)	(7)	(7)	(7)	(7)
	v1616-1622	0.43%	1.08%	1.45%	1.43%	0.00%	6.21%	0.01%	0.09%
	(7)	(7)	(7)	(7)	(7)	(7)	(7)	(7)	(7)
	t0415-0421	0.02 %	0.01 %	0.03%	0.03%
	(7)	(7)	(7)	(7)

Table 4 :

 4 Computational results (relative distance to best-known solutions for solved instances and in parentheses number of solved instances) on CSPLib022 instances[START_REF] Suniel | CSPLib problem 022: Bus driver scheduling[END_REF].

	Instances		baryonyx	cplex	LocalSolver	[9]
		default	fast	good	genoud
	CSPLib022	6.87%
	(12)

Table 6 :

 6 Computational results (relative distance to best-known solutions for solved instances and in parentheses number of solved instances) on nqueens instances.

	Instances		baryonyx		cplex	LocalSolver
		default	fast	good genoud
	nqueens	4.33%	8.11%	5.44%	0.61%	453.26%	36.34%
	(8)	(6)	(7)	(8)	(8)	(8)	(8)

Table . 7

 . : Static parameter configurations for baryonyx.

	Parameters	Configurations WedelinFast WedelinGood baryonyx default
	τ	60 seconds (1800 for VCS)
	Γ	30 parallel runs (only 3 for training)
	ρ		10 5 iterations
	ω	20 warmup iterations
	κ min		0
	κ step	10 -3	2 × 10 -4		10 -3
	κ max		0.6
	α δ	0 0.01	(1 -θ)	1 min n j=1 |cj |,cj =0 max n j=1 |cj |
	θ	0.4	0.6		0.5
	p	0			0.5
	Violated constraint order	none		random
	Initialization policy	bastert	best-cycle

Table . 8

 . : Learnt parameter configurations for baryonyx rgn . Remaining parameters use baryonyx default configuration. After every benchmark name, we put between square brackets the generation number where this configuration has been found by rgenoud.

Table . 9

 . : Computational results (objective value) on SPP instances (in bold: training set and best solutions, "-": no solution found)

	Beasley[5]
	LocalSolver
	cplex
	baryonyx
	Rows
	Columns
	Instances

We set α = 1, p = 0.5 with the best-cycle strategy, and a randomized order of violated constraints at each iteration.

In our experiments, we used a smaller number of cores Γ =

during the training phase than during the test phase (Γ = 30) in order to save computation time.

https://www.localsolver.com

For cplex, we set EP AGAP = EP GAP = EP IN T = 0 to avoid premature stop.

http://people.brunel.ac.uk/ ~mastjjb/jeb/orlib/sppinfo.html

http://www.zib.de/opt-long_projects/TrafficLogistics/Telebus/index.en.html

http://www.csplib.org/Problems/prob022

https://forgemia.inra.fr/thomas.schiex/cost-function-library/tree/master/random/ wqueens

Acknowledgments. We are grateful to the GenoToul (Toulouse) Bioinformatic platform for providing us computational support for this work.


Lauris Evariste
email: lauris.evariste@gmail.com

Emmanuel Flahaut

Clément Baratange

Maialen Barret

Florence Mouchet

Eric Pinelli

Anne Marie Galibert

Brigitte Soula

Laury Gauthief

Ecotoxicological assessment of commercial boron nitride nanotubes toward Xenopus laevis tadpoles and host-associated gut microbiota

Keywords: Boron n itride nanotubes, Xenopus, genotoxicity, gut microbiota, 16S

Despite the growing interest for boron nitride nanotubes (BNNT) due to their unique properties, data on the evaluation of the environmental risk potential of this emerging engineered nanoma terial are currently lacking. Therefore, the ecotoxicity of a commercial form of BNNT (containing tubes, hexagonal-boron nitride, and boron) was assessed in vivo toward larvae of the amphibian Xenopus /oevis. Following the exposure, multiple endpoints were measured in the tadpoles as well as in bacterial communities associated to the host gut. Exposure to BNNT led to boron accumulation in host tissues and was not associated to genotoxic effects. However, the growth of the tadpoles increased due to BNNT exposure. This parameter was associated to remodeling of gut microbiome, benefiting to taxa from the phylum Bacteroidetes. Changes in relative abun dance of this phylum were positively correlated to larval growth. The obtained results support the finding that BNNT are biocompatible as indicated by the absence of toxic effect from the tested nanomaterials. ln addition, byproducts, especially free boron present in the tested prod uct, were overall beneficial for the metabolism of the tadpoles.

1 ntroduction

Boron nitride nanotubes (BNNTI are nanoparticles structurally analogous to carbon nanotubes (CNTs) in which boron and nitrogen replace carbon atoms [START_REF] Chopra | soron Nitride Nanotubes: Sdence[END_REF]. lncreasing attention is being paid to these BNNT due to their unique properties such as high thermo-mechanical stability [START_REF] Suryavanshi | Elastic Modulus and Resonance Behavior of Boron Nitride Nanotubes[END_REF][START_REF] Chen | Boron Nitride Nanotubes: Pronounced Resistance to Oxidation[END_REF]), electrical insulation [START_REF] Radosavljevié | uElectrical Properties and Transport in Boron Nitride Nanotubes[END_REF]) and thermal conductivity (Chang et al. 2006). As high quality materials are increasingly available due to the improvement of synthesis techniques, BNNT are triggering great interest for the development of applications in a wide range of area such as com posite improvement, medicine, hydrogen storage as well as aquatic pollution remediation [START_REF] Kim | Boron Nitride Nanotubes: Synthesis and Applications[END_REF][START_REF] Kalay | synthesis of Boron Nitride Nanotubes and Their applications[END_REF][START_REF] Zhi | Boron Nitride Nanotubes: functionalization and Composites[END_REF][START_REF] Merlo | Boron Nitride Nanomaterials: Biocompatibility and Bio-Applications[END_REF][START_REF] Laie | Boron Nitride for Hydrogen Storage[END_REF][START_REF] Yu | Boron Nitride-Based Materials for the Removal of Pollutants from Aqueous Solutions: A Review[END_REF]. Such keen interest for this material predicts a large-scale production in the near future. As any engineered nanomaterial, BNNT could be released into the environment during the whole material life cycle, from the production, the use and the waste disposai or recycling to reach significant levels in aquatic ecosystems [START_REF] Yu | Boron Nitride-Based Materials for the Removal of Pollutants from Aqueous Solutions: A Review[END_REF][START_REF] Mottier | Environmental Impact of Engineered Carbon Nanoparticles: From Releases to Effects on the Aquatic Biota[END_REF]. For these reasons, this implies to evaluate its biocompatibility and its potential hazard for the environment before considering the mass use of this nanomaterial.

Studies related to BNNT toxicity remain scarce and inconsistent. lndeed, although many of the published works indicate biocompatibility of raw and functionalized BNNT in vitro as well as in vivo [START_REF] Chen | soron Nitride Nanotubes are Noncytotoxic and Can be Functionalized for Interaction with Proteins and Cells•[END_REF][START_REF] Emanet | uchapter 3 -Biocompatibility Evaluation of Boron Nitride Nanotubes[END_REF]Ciofani et al. 2010;[START_REF] Ciofani | cytocompatibility, Interactions, and Uptake of Polyethyleneimine-Coated Boron Nitride Nanotubes by Living Cells: Confirmation of Their Potential for Biomedical Applications[END_REF][START_REF] Ciofani | cytocompatibility Evaluation of Gum Arabic-Coated Ultra-Pure Boron Nitride Nanotubes on Human Cells[END_REF][START_REF] Fernandez-Yague | Effects of Polydopamine Functional ization on Boron Nitride Nanotube Dispersion and Cytocompatibility[END_REF][START_REF] Rocca | uPectin-Coated Boron Nitride Nanotubes: ln Vitro Cyto-/ lmmune-Compatibility on RAW 264.7 Macrophages[END_REF][START_REF] Salvetti | uln Vivo Biocompatibility of Boron Nitride Nanotubes: Effects on Stem Cell Biology and Tissue Regeneration in Planarians[END_REF], some studies highlighted cell type dependent cytotoxic and genotoxic effects [START_REF] Augustine | Assessing Size-Dependent Cytotoxicity of Boron Nitride Nanotubes Using a Novel Cardiomyocyte AFM Assay[END_REF][START_REF] Horvath | In Vitro Investigation of the Cellular Toxicity of Boron Nitride Nanotubes[END_REF][START_REF] Çal | In Vitro Investigation of the Effects of Boron Nitride Nanotubes and Curcumin on DNA Damage[END_REF]. Ecotoxicological data are lacking while it is necessary for the sustainable develop ment of BNNT-based products. lndeed, as for other nanoparticles, a large scale usage of BNNT would lead to releases into the environment at different stages of their life-cycle [START_REF] Mottier | Environmental Impact of Engineered Carbon Nanoparticles: From Releases to Effects on the Aquatic Biota[END_REF]Sun ét al. 2016;Kéllér and Lazaréva 2014;[START_REF] Bundschuh | Nanoparticles in the Environment: Where Do We Come From, Where Do We Go to?[END_REF]. Thus, there is a need to fill knowledge gaps concerning the ecotoxicity of BNNT, particu larly toward aquatic ecosystems as they constitute a receptacle for most contaminants. For this purpose, the amphibian Xenopus laevis was chosen as bio logical model to assess the ecotoxic potential of BNNT as larval stages were shown to be sensitive to contaminants, including nanoparticles [START_REF] Evariste | Gut Microbiota of Aquatic Organisms: A Key Endpoint for Ecotoxicological Studies[END_REF]Mouchet et al. 2010;[START_REF] Bourdiol | siocompatible Polymer-Assisted Dispersion of Multi Walled Carbon Nanotubes in Water, Application to the Investigation of Their Ecotoxicity Using Xenopus laevis Amphibian Larvae[END_REF][START_REF] Muzi | uExamining the Impact of Multi-Layer Graphene Using Cellular and Amphibian Models[END_REF][START_REF] Saria | ushort Term Exposure to Multi-Walled Carbon Nanotubes lnduce Oxidative Stress and DNA Damage in Xenopus /aevis Tadpoles[END_REF][START_REF] Colombo | reratogenic Hazard of BP8-Coated Silver Nanoparticles to Xenopus /aevis[END_REF][START_REF] Marfn-Barba | unravelling the Mechanisms That Determine the Uptake and Metabolism of Magnetic Single and Multicore Nanoparticles in a Xenopus laevis Model[END_REF]. Parameters related to the larval growth were shown to constitute endpoints of interests as they reflect the overall health status of the organisms. ln a con text of exposure to nanoparticles of carbon allo tropes such as CNTs, toxicity leading to larval growth inhibition was described by the specific sur face area of the tested nanoparticles [START_REF] Mottier | surface Area of Carbon Nanoparticles: A Dose Metric for a More Realistic Ecotoxicological Assessment[END_REF][START_REF] Lagier | surface Area of Carbon Based Nanoparticles Prevai ls on Dispersion for Growth Inhibition in Amphibians[END_REF]. ln addition, genotoxic effects could lead to non-negligible consequences as damaged DNA may cause cellular dysfunctions [START_REF] Jackson | The DNA-Damage Response in Human Biology and Disease[END_REF] leading to the death of organisms and further impact the ability of popula tions to maintain [START_REF] Sukumaran | Effects of Genotoxicity and lts Consequences at the Population Level in Sexual and Asexual Artemia Assessed by Analysis of Inter-Simple Sequence Repeats (ISSR)[END_REF]. For this reasons, biomarkers related to genotoxicity are considered as a pertinent ecological endpoint at the population, community and ecosystem level [START_REF] Anderson | Genetie and Molecular Ecotoxicology: A Research Framework[END_REF]. Thus, alterations of genetic materials in amphibians would contribute to the decline of populations while this class contains many endangered species [START_REF] Oertli | uconservation and Monitoring of Pond Biodiversity: introduction[END_REF].

lt is now widely accepted that the growth of organisms is dependent from its capacity to acquire energy from trophic resources while this energy metabolism is can be strongly influenced by the metabolic capacities of the gut microbiota [START_REF] Mithieux | Gut Microbiota and Host Metabolism: What Relationship[END_REF][START_REF] Musso | #Interactions between Gut Microbiota and Host Metabolism Predisposing to Obesity and Diabetes[END_REF][START_REF] Cani | Metabolism in 2013: The Gut Microbiota Manages Host Metabolism: Nature Reviews[END_REF]. ln addition, the gut microbiome is emerging as a central target in environmental toxi cology studies due to its contribution in the regula tion of multiple physiological processes of the host but also in the metabolization of environmental pollutants (Claus, Guillou, and Ellero-Simatos 2016;[START_REF] Evariste | Gut Microbiota of Aquatic Organisms: A Key Endpoint for Ecotoxicological Studies[END_REF][START_REF] Adamovsky | The Gut Microbiome and Aquatic Toxicology: An Emerging Concept for Environmental Health[END_REF]. Alterations of the gut microbiome composition may dysregulate the normal physiological functioning of the host and lead to various diseases [START_REF] Durack | rhe Gut Microbiome: Relationships with Disease and Opportunities for Therapy[END_REF]. Furthermore, in the amphibians, changes in the gut microbiota composition during ear1y stages of life were shown to influence later life resistance to infections [START_REF] Knutie | Early-Life Disruption of Amphibian Microbiota Decreases Later-Life Resistance to Parasites[END_REF]. Environmental pollution constitutes a major factor influencing the composition of the gut flora and previous studies reported that a wide range of con taminants were able to induce gut dysbiosis in association to deleterious issues on host [START_REF] Evariste | Gut Microbiota of Aquatic Organisms: A Key Endpoint for Ecotoxicological Studies[END_REF][START_REF] Jin | Effects of Environ mental Pollutants on Gut Microbiota[END_REF]. Thus, the present work aims to investigate the consequences of an in vivo exposure of X. laevis tadpoles to BNNT com bining monitoring of host physiology and gut bac terial community composition. The former endpoint constitutes an integrator of the overall host-metab olism while the latter one would allow to better understand the causes of potential metabolic changes induced by BNNT exposure. Overall, the results obtained in this study suggest that BNNT are rather biocompatible. lndeed, exposure to the nanoparticles did not induced irreversible genotoxic effects and stimulated the growth of the tadpoles. This increased growth was associated to remodeling of the host-associated gut microbiota benefiting to the phylum Bacteroidetes. Measurement of boron released in the exposure media and accumulated in the organisms suggest that the free boron associ ated to the nanoparticles are responsible for the beneficial effects observed.

Materials and methods

1. Synthesis and characterization of boron nitride nanotubes

Boron nitride nanotubes were obtained from BNNT LLC (USA). They correspond to the 'BNNT Pl-Beta Products' range of materials as classified on the website of the company (http://www.bnnt.com/ products). This material is described as containing hexagonal boron nitride (h-BN) and boron as well as being catalyst free (elemental B and N > 99.90/2) with ca. 25% of Boron in addition to BNNT. According to the manufacturer, the specific surface area is >200 m 2 /g. Sample was used as-received, after simple grinding (agate mortar and pestle) in order to make it suitable for the experiments and characterization. The sample was characterized by X-ray diffraction (XRD), Brunauer-Emmet-Teller (BET), TGA, TEM, Raman and chemical elemental analysis.

X-ray diffraction analysis was performed using a Bruker D4 ENDEAVOR X-ray Diffractometer (Cu.Kcx. = 1.5406 Â). Specific surface area was meas ured with a Micrometrics Flow Sorb Il 2300 by N 2 adsorption using the BET theory. Dry BNNT powder sample was first degassed at 100 °C under N 2 atmosphere for 2 h before being cooled down to 77 K (liquid nitrogen temperature) for adsorption of a monolayer of N 2 • After heating the sample back to room temperature, the desorption peak was recorded. A calibration was done by injecting a known quantity of N 2 (in the same conditions as for the sample measurement). The uncertainty of the measurement was 3%. Thermogravimetric analysis was performed in air atmosphere at a heating rate of 1 °C/min using a SETARAM TAG 16. Transmission Electron Microscopy observations were performed using a TEM JEOL 1400 with an acceleration voltage of 120 kV after deposition of a few drops of a dilute suspension in ethanol on Copper TEM grids (Lacey carbon). Outer diameter of the BNNTs was meas ured from HRTEM images using lmageJ software. For chemical elemental analysis of Boron, the sam ple was digested in an open system with a mixture of sulfuric and nitric acids du ring 8 h. Elemental analysis was performed on an ICP-AES ICAP 6500 from Thermofisher Scientific, Bremen. Both the digestion and analysis were performed by CREALINS, 6Napse group (uncertainty of the meas urement: 3%). Nitrogen and Oxygen elemental anal yses were performed by ISA (CNRS, Lyon, France) by organic micro-analysis after total combustion at 1050 °C under helium/oxygen flux for nitrogen; total pyrolysis at 1080 °C under nitrogen flux for Oxygen. Raman signature was analyzed to get information on the structural quality of the nanotubes (Labram HR800 Horiba Yvon Jobin, À= 532 nm). X-Ray photoelectron spectroscopy (XPS) was used to determine the quantitative atomic composition of the BNNT (XPS Kalpha ThermoScientific). IR analysis was performed using a Perkin Elmer Spectrum One FT-IR spectrometer and inclusion of the sample in a dry KBr pellet.

The stability of the BNNT suspension was eval uated using a Turbiscan LAB (Formulation) equip ment, at room temperature. Suspensions were prepared at 10 mg L 1 in exposure media before performing analysis. The suspension was placed in a glass vial and both transmission and reflection were monitored vs time for 24 h (1 scan per minute).

Xenopus laevis rearing, breeding and exposure conditions

Xenopus laevis originated from our certified breed ing facilities (under the approval number A31113002) and the experimental procedure was approved by an ethic committee (CEEA-073). Rearing and breeding conditions of Xenopus were performed as previously described [START_REF] Mouchet | ucharacterisation and in Vivo Ecotoxicity Evaluation of Double-Wall Carbon Nanotubes in Larvae of the Amphibian Xenopus /aevis[END_REF]Mouchet et al. 2010). Pregnant mare's gonadotropin was injected in sexually mature indi viduals to induce spawning. Fecundated eggs obtained were bred in active charcoal filtered tap water at 22 ± 2 °C and fed with ground aquarium fish food (TetraPhyll®) until they reach stage 50 according to Nieuwkoop & Faber development table [START_REF] Nieuwkoop | uNormal Table of Xenopus /aevis (Daudin). a Systematical and Chronological Survey of the Develop ment from the Fertilized Egg till the End of Metamorphosis[END_REF]. Following the international standard ISO 21427-1 procedure (ISO/ FDIS 21427-1. [START_REF] Harrison | Quantification of hexagonal boron nitride impurities in boron nitride nanotubes via FTIR spectroscopy[END_REF]2006), 20 larvae per experimental condition were exposed for 12 days under semi static conditions with daily exposure media renewal. Larvae were fed daily ad libitum with ground aquar ium fish food (TetraPhylr®, Tetra, Melle, Germany). Food was added twice a day, directly after media renewal and at the end of the day. Unconsumed food was removed daily during the media renewal process. Negative control (NC) condition was com posed of reconstituted water devoid of contaminant (294 mg L 1 CaCl 2 .2H 2 O; 123.25 mg L 1 MgSO 4 .7H 2 O; 64.75 mg L 1 NaHCO 3 ; 5.75 mg L 1 KCI) [START_REF] Evariste | Gut Microbiota of Aquatic Organisms: A Key Endpoint for Ecotoxicological Studies[END_REF][START_REF] Mouchet | ucharacterisation and in Vivo Ecotoxicity Evaluation of Double-Wall Carbon Nanotubes in Larvae of the Amphibian Xenopus /aevis[END_REF]) while posi tive control (PC) for genotoxicity assessment con tained cyclophosphamide monohydrate ([6055-19-2], Sigma, France) at 40 mg L 1

• BNNT were dis persed by sonication before contamination of the exposure medium. Tested concentrations of BNNT were 0, 0.1, 1, and 10 mg L 1 •

Analys is of boron clearance in the water and accumulation in tadpoles

Boron concentration was measured in reconstituted 1 water from the control and the 1 0 mg L of BNNT conditions after 24 h of incubation in the presence and absence of tadpoles in the media. Water sam ples were centrifuged twice at 1500 x g during 20 min and filtered at 0.22 µm to remove nanotubes from the samples before boron analysis. After 12 days of exposure, boron concentration was measured in larvae from the control and the 1 0 mg L 1 of BNNT gr oups. Guts of tadpoles were removed before lyophilization of organisms, in order to avoid measurement of BNNT-related boron from the gut content. Boron concentration in the water and in the larvae was determined by ICP-AES at CREALINS laboratory (Lyon, France).

Micronucleus test and ce// cycle analysis

The micronucleus assay was performed in accord ance with the ISO 21427-1 guidelines (ISO/FDIS 21427- 1.1, 2006). The micronucleus formation con stitutes a good indicator of irreversible genotoxic effects, integrating aneugenic and clastogenic effects occurring after exposure to a genotoxic compound. For this purpose, after 12 days of expos ure, Xenopus larvae were anesthetized by bathing in MS222 solution at 100 mg L 1 before collecting blood samples from cardiac puncture. Blood smears were prepared, fixed in methanol for 10 minutes before performing hematoxylin and eosin staining.

Micronucleated erythrocytes were accounted over a total of 1000 cells (MNE %0) using optical microscopy (oil immersion lens, xl 500). For cell cycle analysis, blood sub-samples were fixed using cold ethanol (70% v/v). Prior to the flow cytometric analysis, cells were rinsed using PBS and labeled with FxCycle ™ PI/RNase Staining Solution (Life Technologies SAS) following manufacturer's recom mendations. Propidium iodide fluorescence was measured using Cytoflex (Beckman Coulter, USA) equipped with a 488-nm excitation laser. For each sample, 10,000 events were acquired in a region corresponding to cells of interest after doublet dis crimination using FSC-H versus FSC-A.

La,va/ growth measurement

The total length of larvae was measured at the beginning (d0) and at the end of the 12 days of exposure (dl 2) using lmageJ 1.49 software. Normalized growth rate was determined as previ ously described [START_REF] Mottier | surface Area of Carbon Nanoparticles: A Dose Metric for a More Realistic Ecotoxicological Assessment[END_REF][START_REF] Lagier | surface Area of Carbon Based Nanoparticles Prevai ls on Dispersion for Growth Inhibition in Amphibians[END_REF]) using the following formula:

. .

(Ldl 2-MLd0
)

Normal,zed SIZe (%) =

MLd0

x 100

X (M::12)

Ld12 corresponds to the length of one larvae at 12 days, MLd0 is the mean length at day 0 of larvae from the exposure condition and MLCdl 2 is the mean length of larvae from the negative control at day 12.

Analysis of sequences from gut microbiota survey

Genomic DNA was extracted from the whole intes tine of larvae using DNeasy PowerSoil kit (QIAGEN) according to manufacturer recommendations with the following adjustments: samples were incubated 1 0 min at 65 °C after adding the solution Cl and elution buffer C6 was incubated during 10 min before performing the last centrifugation. The V4-V5 region of the 16S rRNA gene was amplified using 515 F (5'-GTGYCAGCMGCCGCGGTA-3')/928R (5'-CCCCGYCAATTCMTTTRAGT-3') primer pair [START_REF] Wang | conservative Fragments in Bacterial 16S rRNA Genes and Primer Design for 16S Ribosomal DNA Amplicons in Metagenomic Studies[END_REF] and the following PCR protocol: 94 °C for 120 s, 30 cycles of 94 °C for 60 s, 65 °C for 40s, 72 °C 30s and 72 °C for l0min. Amplicon sequencing was performed using an Illumina MiSeq (2 x 250 pb) by the Get PlaGe platform (Genotoul, Toulouse, France).

Demultiplexed data were processed using FROGS (Find Rapidly OTU with Galaxy Solution) pipeline on Galaxy [START_REF] Escudié | FROGS: Find, Rapi dly, OTUs with Galaxy Solution[END_REF]. Contigs with a length between 380 and 500 pb were kept, clustered with Swarm [START_REF] Mahé | swarm: robust and Fast Clustering Method for Amplicon-Based Studies[END_REF]) with an aggregation dis tance of 3 and chimeras were removed. Filters were applied to remov e singletons and keep for analysis OTUs with a minimum abundance of 0.005% of the sequences [START_REF] Bokulich | ouality-Filtering Vastly lmproves Diversity Estimates from Illumina Amplicon Sequencing[END_REF]). The taxonomy affiliation was performed using Blastn against the Silva 132 database (pintail 80).

Statistical analysis

Results from micronucleus assay were analyzed using McGill non-parametric test [START_REF] Mcgill | Variations of Box Plots[END_REF] on median values of each group of larvae. This test consists in comparing medians of samples of size n (where n 2: 7) and in determining their 95% confidence intervals (95% Cl). 95% Cl are expressed by M ± 1.57 x IQR/ ✓n, where M is the median and IQR is the inter-quartile range [START_REF] Mcgill | Variations of Box Plots[END_REF]. The difference between the medians of the test groups and the median of the NC group is significant with 95% certainty if there is no overlap. For cell-cycle data, normality was assessed with Kolmogorov-Smirnov test and homogeneity of variances with Levene's test. One way analysis of variance (ANOVA) followed by Tukey test were used to compare cell-cycle phase distribution among conditions.

For data manipulations of gut microbiota survey, OTUs counts, alpha diversity indexes and Weighted Unifrac Distances calculations as well as multidi mentional scaling (MDS) plot were carried out using 'Phyloseq' R package [START_REF] Mcmurdie | Phyloseq: An R Package for Reproducible Interactive Analysis and Graphies of Microbiome Census Data[END_REF]. Graph visualization of OTUs relative abundances was performed using 'ggplot2' package [START_REF] Wickham | ggplot2: Elegant Graphies for Data Analysis[END_REF] while differential abundance of bacterial gen era between exposed conditions compared to the control group was performed using 'Deseq2' R package [START_REF] Love | Moderated Estimation of Fold Change and Dispersion for RNA-seq data with DESeq2[END_REF]. For multi variate analysis of variance between groups, PERMANOVA was performed using Adonis function from the 'vegan' R package [START_REF] Oksanen | Package 'vegan[END_REF]. 10-15nm) (Figure l(c,d)) as well as some more disor ganized material. Part of the nanotubes are damaged or even broken (Figure 1 (c,d), black arrows), which may be related to the grinding treatment that we have applied. On Figure 1 (c,d), it is possible to see that while most nanopartides look amorphous, a few are crystalized. Ali observed nanoparticles were coated with similar polyhedral graphitic-like h-BN shells. Outer diameter of the BNNT ranged from 2 to 14 nm with a mean outer diameter of 6 ± 2.6 nm (Figure 2).

Results and discussion

3.

We obtained through measurement of elemental an alysis a total boron content between 38.6 and 42.5 wt. % (ca. 50.5 at. % using the highest weight value), 46.3 wt. % of Nitrogen (ca. 42.5 at. %) and 8.7 wt. % of Oxygen (ca. 7 at. %). If we make a first hypothesis that all the Nitrogen in the sample is in BN, then the excess of Boron is of ca. 8 at. %.

Assuming that all the oxygen is involved in B 2 0 3 (the most stable boron oxide), then the excess of Boron (corresponding to elemental boron) would be ca. 3.3 wt.%. This is rather far from the informa tion from the provider, although this explains that variations are likely between different batches. Details results from the XRD, XPS, TGA, Raman and IR analysis are detailed in the Supplementary information .

The variation of the transmission data (compared to t0) vs time at the top and at the bottom of the vial is presented in Figure 3(A,B) respectively. At the top of the vial, transmission progressively increased as the BNNTs were settling down. On the opposite, the trans mission progressively decreased at the bottom of the vial as the sediment was forming. ln both cases, it is obvious that most of the clarification/sedimentation took place in 8-10 h, at a rather constant speed (before slowing down) (Supplementary Figure SS).

Boron concentration in the exposure media and larvae

Analysis of boron concentration in the water indicated that reconstituted water from the control condition was devoid of boron. ln the 10 mg L 1 of BNNT con dition, 0.32 mg L 1 of boron were measured after 24 h of incubation in the media in absence of Xenopus tadpoles while after 24 h in presence of larvae, the concentration of boron decreased below 0.1 mg L 1 (quantification limit). During the exposure, larvae pro duced BNNT-containing feces according to the color differences with the control group. At the end of the 12 days of exposure, 2 ± 4 µg g 1 of boron were meas ured in the whole tissues of the larvae from the con trai group. Boron concentration in larvae exposed to 10 mg L 1 of BNNT reached 183 ± 83 µg g 1 of tissue.

BNNT were shown to be structurally very stable [START_REF] Chen | Quantitative Characterization of Structural and Mechanical Properties of Boron Nitride Nanotubes in High Temperature Environments[END_REF][START_REF] Song | #Excellent Oxidation Endurance of Boron Nitride Nanotube Field Electron Emitters[END_REF]) and there is few chance for the digestive tract of tadpoles to consti tute drastic enough conditions allowing BNNT deg radation as strong acidic conditions are needed to dissolve this material [START_REF] Kleinerman | Dissolution and Characterization of Boron Nitride Nanotubes in Superacid: Langmuir[END_REF]. ln addition, carbon-based struc tural homologues of BNNT were shown to accumulate in the lumen of Xenopus tadpoles during the exposure while a transfer across the intestinal wall was not evidenced [START_REF] Mouchet | ulnternational Amphibian Micronucleus Standarclized Procedure (ISO 21427-1) for in Vivo Evaluation of double-walled carbon nanotubes toxicity and genotoxicity in water[END_REF]. Thus, crossing of intestinal barrier is not likely to occur in the case of BNNT exposure and it is highly probable that the boron concentration measured in the water as well as in the tissues of the larvae originated from free boron as byproduc t found in the commercial form of the tested BNNT. However to confirm this hypothesis, the use of imaging technics would be needed to determine if BNNT are present or not in the bio logical tissues. lndeed, amphibians were previously shown to accumulate free boron in polluted sites and are able to tolerate boron concentration in the water up to 10 mg L 1 for extended period without adverse effects [START_REF] Eisler | Boron hazards to fish, wildlife, and inverte brates: a synoptic review[END_REF][START_REF] Emiroglu | soron Concentration in Water, Sediment and Different Organisms around Large Borate Deposits of Turkey[END_REF].

Measurement performed under the experimental conditions indicated a full clearance of the water column by the larvae over 24 hours, leading to rejection of BNNT in the animal feces.

Micronucleus assay and ce// cycle analysis

Xenopus larvae exposed to cyclophosphamide at 40 mg L 1 (positive control = PC) exhibited signifi cantly higher MNE compared to the control group (negative control = NC), validating results obtained from the micronucleus assay. Among experimental conditions containing BNNT, although a trend to an increasing number of micronucleated cells was observed after exposure to 10 mg L 1 of the nano particle, no statistically significant increase of MNE %0 was accounted (Figure 4(A)). Analysis of erythro cyte cell cycle high ligh ted a significant decrease in S-phase cells in conditions containing BNNT com pared to the negative control (ANOVA, S-phase: p < 0.001 ), while the percentage of G0/G 1 ce lis increased and G2/M decreased only after exposure to 1 mg L 1 of BNNT (ANOVA, G0/G 1: p < 0.001; G2/ M: p < 0.001) (Figure 4(B)).

The results obtained indicated that the commer cial form of BNNT tested was not leading to irre versible genotoxic effects at the tested concentrations as no significant induction of micro nucleated erythrocytes was observed. Thus, even if we cannot exclude that reversible DNA damages could occur through single strand break DNA, the obtained results indicated that BNNT exposure do not lead to permanent DNA alteration. This obser vation is contrary to in vitro studies indicating dis turbances of cell proliferation and genotoxicity of BNNT [START_REF] Fernandez-Yague | Effects of Polydopamine Functional ization on Boron Nitride Nanotube Dispersion and Cytocompatibility[END_REF][START_REF] Horvath | In Vitro Investigation of the Cellular Toxicity of Boron Nitride Nanotubes[END_REF][START_REF] Çal | In Vitro Investigation of the Effects of Boron Nitride Nanotubes and Curcumin on DNA Damage[END_REF][START_REF] Emanet | interaction of Carbohydrate Modified Boron Nitride Nanotubes with Living Cells[END_REF], but consistent with the few available works per formed in vivo. lndeed, no genotoxicity was observed in the planaria [START_REF] Salvetti | uln Vivo Biocompatibility of Boron Nitride Nanotubes: Effects on Stem Cell Biology and Tissue Regeneration in Planarians[END_REF] as well as in drosophila [START_REF] Demir | Antigenotoxic Potential of Boron Nitride Nanotubes[END_REF] after BNNT exposure. ln the latter study, BNNT showed antioxidant and antigenotoxic properties against a known genotoxic compound. ln addition, similar effects were observed in Xenopus laevis exposed to other non-oxidized 1 D carbon counterparts of BNNT [START_REF] Mouchet | ucharacterisation and in Vivo Ecotoxicity Evaluation of Double-Wall Carbon Nanotubes in Larvae of the Amphibian Xenopus /aevis[END_REF]Mouchet et al. 2010). lt was previously indicated that baron exert protective effects against genotoxic compounds in other bio logical models [START_REF] Turkez | Effects of Borie Acid and Borax on Titanium Dioxide Genotoxicity[END_REF][START_REF] Ince | Protective Effects of Boron on Cyclophosphamide lnduced Lipid Peroxidation and Genotoxicity in Rats[END_REF][START_REF] Ql | #Evaluation of Genotoxic and Antigenotoxic Effects of Boron by the Somatic Mutation and Recombination Test (SMARn on Drosophila[END_REF][START_REF] Tepedelen | Borie Acid Reduces the Formation of DNA Double Strand Breaks and Accelerates Wound Healing Process[END_REF][START_REF] Alak | Borax Supplementation Alleviates Hematotoxicity and DNA Damage in Rainbow Trout (Oncorhynchus mykiss) Exposed to Copper[END_REF]. For this reason, it is not pos sible yet to determine in these conditions whether free baron associated to the commercial mixture is protective against the genotoxic potential of the different types of BN nanoparticles or if BN are devoid of genotoxic potential.

Larval growth rate

Exposure to BNNT led to a significant increase of Xenopus tadpole growth in a dose-dependent man ner (ANOVA, p < 0.001) (Figure S(A)). Growth rate increased from 20.54 ± 9.53% at 0.1 mg L 1 to 44 ± 13.97% at the highest tested concentration compared to negative contrai. However, mean weight/length ratio remained unchanged in larvae exposed to the tested conditions (ANOVA,

p = 0.143) (Figure S(B)).
Exposure to BNNT led to a dose-dependent increase of larval growth that is counterintuitive regarding to the literature indicating that no matter the number of dimension of the carbon-based nanomaterial (CBN) considered, growth inhibition of Xenopus tadpoles was described by the surface area of exposure [START_REF] Mottier | surface Area of Carbon Nanoparticles: A Dose Metric for a More Realistic Ecotoxicological Assessment[END_REF][START_REF] Lagier | surface Area of Carbon Based Nanoparticles Prevai ls on Dispersion for Growth Inhibition in Amphibians[END_REF]. As the surface area of the tested BNNT is compar able to those of the CBNs used in these studies and the observed effects are oppos ite, we can empha size that this growth inhibition mo del is specific to CBNs and is thus not applicable for BNNT. The main hypothesis for mechanisms underlying growth inhibition following CBNs exposure was associated to nutrient depletion and reduction of nutrient intake due to the presence of agglomerated CBNs in the gut [START_REF] Lagier | surface Area of Carbon Based Nanoparticles Prevai ls on Dispersion for Growth Inhibition in Amphibians[END_REF]. The accumulation of BNNT in the tadpole gut was not as clear as previ ously observed with carbon-based nanomaterials (CBNs) that is consistent with the absence of growth alteration measured. Thus, despite the cap acity of BNNT to adsorb nutrients [START_REF] Farmanzadeh | interaction of Vitamins A, 81, C, 83 and D with Zigzag and Armchair Boron Nitride Nanotubes: A DFT Study[END_REF], the hypothesis of nutrient intake limitation is not likely to occur in the case of BNNT exposure.

Previous studies demonstrated that a deficiency (< 0.003 mg L 1 ) as well as too high concentrations of boron (<50 mg L 1) were detrimental for the development of amphibian embryo and impaired the reproduction of adults [START_REF] Fort | Adverse Reproductive and Developmental Effects in Xenopus from lnsufficient Boron[END_REF][START_REF] Fort | chronic Feeding of a Low Boron Diet Adversely Affects Reproduction and Development in Xenopus /aevis[END_REF][START_REF] Laposata | Effects of Boron and Nitrate on Hatching Success of Amphibian Eggs[END_REF][START_REF] Fort | impact of Boron Deficiency on Xenopus /aevis[END_REF]. According to these studies, the boron con centrations measured in the exposure media of this work (0.32 mg L 1

) do not represent a critical con centration impairing tadpole's physiology that is consistent. Despite the lack of literature regarding the enhanced growth of amphibian larvae exposed to boron, these effects were observed in fishes with no explanation of the mechanism involved in such growth increase [START_REF] Eckhert | soron Stimulates Embryonic Trout Growth[END_REF][START_REF] Rowe | The Response of Trout and Zebrafish Embryos to Low and High Boron Concentrations is U-Shaped[END_REF]Oz, lnanan, and Dikel 2018). However, it was sug gested that boron was playing a role in the thyroid axis of X. laevis, probably in the synthesis of T3 [START_REF] Fort | impact of Boron Deficiency on Xenopus /aevis[END_REF]). Thus, among the possible path ways involved in this observed growth stimulation, it may in part occur from a stimulation of the meta morphosis by the boron.

Gut microbiota survey

Among the different tested conditions, richness was not affected (Observed species: ANOVA p =0.505; Chao 1: ANOVA p =0.826). However, evenness calculated with Shannon index was shown to be significantly decreased after 12 days of exposure to l0 mgL 1 of BNNT (ANOVA p<0.001) (Figure 6(A)). Bacterial communities were shown to be signifi cantly affected by the BNNT concentration as révéaléd by MDS pérforméd with Wéightéd-Unifrac distances and PERMANOVA analysis (F = 14.146; r2 = 0.726, p = 0.001) (Figure 6(B)). Pairwise compari sons indicated that gut bacterial communities were similar between the control group and the lowest tested concentration 0.1 mg L 1, while the two other tested concentrations 1 and 1 0 mg L 1 were different from each other and from any other conditions.

Three major phyla composed the gut microbiota of X. laevis tadpoles: Bacteroidetes, Proteobacteria and Firmicutes with relative abundance in the con trol group representing 62 ± 8.5%, 24.7 ± 8.3% and 12.2 ± 5.5% of the whole microbial communities respectively (Figure 7). Exposure to BNNT led to a significant decrease of Proteobacteria relative abun dance at 1 and 10 mg L 1 of BNNT (ANOVA p < 0.001 ). On the opposite, a significant increase of Bacteroidetes is noticed at these concentrations (ANOVA p <0.001), leading to a phylum relative abundance of 74 ± 1.17% at the highest tested concentration.

At the family level of the phylum Bacteroidetes, Bacteroidaceae were shown to significantly increase (ANOVA p < 0.001) at 1 and 10 mg L 1 reaching 97.36 ± 1.16% of the overall Bacteroidetes in the higher BNNT concentration compared to the 82.24 ± So/2 observed in the control group (Supplementary Figure S6A ). On the contrary, other families such as Spirosomaceae, Flavobacteriaceae, Crodnitomicaceae or Weeksellaceae were shown to significantly decrease (ANOVA p <0.001; p <0.001; p = 0.0356; p = 0.0023 respectively). Similarly, inside the phylum Proteobacteria, taxa from the family Magnetospirillaceae were shown to significantly increase after exposure to 10 mg L 1 of BNNT (ANOVA p =0.00116) (Supplementary Figure S6B), while Moraxellaceae and Burkholderiaceae were shown to significantly decrease at this concentration (ANOVA p < 0.001 and p = 0.00811 respectively).

At the genus scale, 12, 27, and 39 genera were differentially observed between the control group and the BNNT concentration of 0.1, 1, and 10 mg L 1 respectively (excluding multi-affiliated and unknown genera). Thus, even if the Firmicutes phylum was not globally affected, some genera were shown to be differentially observed between the control and the 1 and 1 O mg L 1 of BNNT conditions (Supplementary Figure S7). Genera which relative abundance was shown to be unaffected, decreased or increased were clas sified as resistant, sensitive or opportunistic respectively. Thus, each category represented about 58.2% of the overall taxa for the resistant, 21.5 and 20.3% for the sensitive and opportunis tic respectively. Even if it is difficult to differentiate the part of the effects associated to indirect or direct effects of BNNT and/or to byproducts exposure, several possi bilities can be considered concerning the effects observed on gut bacterial communities after the exposure of tadpoles to this commercial BNNT.

Considering the potential indirect effects, as pre viously mentioned, boron could play a role in the :.. -:.. :.. thyroid axis that is involved in the development and maturation of the gut [START_REF] Chalmers | isotope Effect on the Thermal Conductivity of Boron Nitride Nanotubes[END_REF][START_REF] Lshizuya-Oka | Amphibian Organ Remodeling during Metamorphosis: insight into Thyroid Hormone-lnduced Apoptosis[END_REF]. lt was also demonstrated that gut bacterial communities were changing dur ing the metamorphosis of amphibians, notably through an increase of the phyla Bacteroidetes, a decrease of Proteobacteria and decrease of diversity [START_REF] Kohl | Restructuring of the Amphibian Gut Microbiota through Metamorphosis: The Amphibian Gut Microbiota[END_REF][START_REF] Chai | changes in Intestinal Microbiota of Bufo Gargarizans and lts Association with Body Weight During Metamorphosis[END_REF]) which is consist ent with our data. Thus, we cannot exclude that the changes occurring in the gut microbial communities could be in part associated to the stimulation of thyroid axis leading to the acceleration of gut maturation. These changes in bacterial consortium could also be associated to direct effects of BNNT. Surprisingly, antibacterial properties of BN nanoparticles have not been thoroughly investigated yet [START_REF] Merlo | Boron Nitride Nanomaterials: Biocompatibility and Bio-Applications[END_REF] and results available from the literature regarding the antibacterial activity are contradictory. lndeed, despite its potential to interact with bio logical membranes and to form heteroaggregates with bacteria [START_REF] Thomas | insertion Mechanism and Stability of Boron Nitride Nanotubes in Lipid Bilayers[END_REF][START_REF] Wang | Heteroaggregation of Nanoparticles with Biocolloids and Geocolloids[END_REF], it was indicated that raw BNNT at concentrations reaching 1 mg L 1 were not inducing E.coli and S.aureus gr owth inhibition [START_REF] Nithya | uAqueous Dispersion of Polymer Coated Boron Nitride Nanotubes and Their Antibacterial and Cytotoxicity Studies[END_REF]. ln addition, expos ure of P. aeruginosa strain to boron nitride nanosheets at a concentration of 10 mg L induced damages to bacterial membranes in E.coli [START_REF] Zhang | Nanotoxicity of Boron Nitride Nanosheet to Bacterial Membranes: Langmuir: The[END_REF]. These studies were performed through exposure of single bacterial species in sus pension in a liquid medium that is very different from the exposure conditions found in the gut. However, few data emphasizing the effects of BN nanoparticles on complex bacterial communities are available. Strong antibiofilm activities of 2 D BN were observed in established biofilm composed of single bacterial species (Klvanç et al. 2018).

� -=. � 3 !2 � � 3 3 3 3 � � � � 1 mgJL 0 -- --=. � !2 � � 3 3 3 3 3 � � � � � 10 mg/L -� � � � � !2 � :è 3 � � � � � � r r r r � � � r Phylum _J Bacteroldetes D Firmicutes ■ Fusobacteria ■ Proteobacteria
As boron compounds were shown to exert weak toxicity against multiple bacterial strains (Minimum inhibitory concentration over 1 g L 1) that possess a good tolerance to this element [START_REF] Ahmed | Mechanism of Boron Tolerance in Soil Bacteria[END_REF][START_REF] Yilmaz | Minimum lnhibitory and Minimum Bactericidal Concentrations of Boron Compounds against Several Bacterial Strains, Turk[END_REF][START_REF] Sayin | uAntibacterial and Antibiofilm Effects of Boron on Different Bacteria[END_REF], few direct effects from free boron could be expected on gut microbial com munities. Boron was shown to constitute an essen tial trace element for bacteria, contributing to physiological and metabolic activities [START_REF] Kabu | Biological Effects of Boron: Reviews of Environmental Contamination and Toxîcology[END_REF]Uluisik, Karakaya, and Koc 2018). Furthermore, boron was shown to be implied in bacterial communication through the activation of an extracellular signaling molecule (autoinducer Al-2), involved in quorum sensing [START_REF] Chen | structural Identification of a Bacterial Quorum-Sensing Signal Containing Boron[END_REF]. This process allows the regulation of gene expres sion and diverse physiological activities in response to fluctuations in bacterial density [START_REF] Miller | Quorum Sensing in Bacteria[END_REF][START_REF] Federle | Interspecies Communication in Bacteria[END_REF][START_REF] Waters | Quorum Sensing: cell to-Cell Communication in Bacteria[END_REF][START_REF] Papenfort | #Quorum Sensing Signal-Response Systems in Gram-Negative Bacteria[END_REF]. Exposure to nanoparticles were previously shown to influence quorum sensing process [START_REF] Singh | uBactericidal, Quorum Quenching and anti-Biofilm Nanofactories: A New Niche for Nanotechnologists[END_REF][START_REF] Xiao | Impairment of Biofilm Formation by TiO2 Photocatalysis through Quorum Quenching[END_REF][START_REF] Hayat | Quorum Quenching: Role of Nanoparticles as Signal Jammers in Gram-Negative Bacteria[END_REF]). This could also be possible for the BN nanoparticle exposure. Thus, disturbances of quorum sensing from the gut bac teria could lead to changes in microbial commun ities [START_REF] Thompson | Manipulation of the Quorum Sensing Signal Al-2 Affects the Antibiotic-Treated Gut Microbiota[END_REF].

Statistical correlations

Correlation analysis indicated strong positive corre lations among tadpole biometric parameters. Significant negative correlation between the phy lum Proteobacteria and Bacteroidetes are noticed (Pearson, r = -0.91 p < 0.001). Firmicutes/ Bacteroidetes ratio was shown not to be correlated with any of the larvae biometric parameters or growth rate. However, phylum Bacteroidetes was shown to be significantly and positively correlated with larval length (r = 0.62, p = 0.0034), growth rate (r = 0.51, p = 0.0211 ). On the contrary, phylum Proteobacteria was shown to be significantly nega tively correlated with larval length (r = -0.69, p=0.0008) and growth (r = -0.55, p=0.012).

Changes in the Firmicutes/Bacteroidetes ratio of gut microbiota was shown to be associated to metabolic disorders in several species [START_REF] Ley | obesity Alters Gut Microbial Ecology[END_REF][START_REF] Ley | Ecological and Evolutionary Forces Shaping Microbial Diversity in the Hu man Intestine[END_REF][START_REF] Haiser | 0eveloping a Metagenomic View of Xenobiotic Metabolism[END_REF][START_REF] Li | Gut Microbiota Contributes to the Growth of Fast-Growing Transgenic Common Carp (Cyprinus carpio L.)[END_REF]. ln this study, the F/ B ratio was not significantly affected (ANOVA, p =0.2) and was not correlated to growth parame ters following exposure to BN. Thus, we can suggest that the growth stimulation observed is not associ ated to the potential induction of a metabolic disorder in host. On the contrary, phylum Bacteroidetes was positively correlated with growth parameters of the tadpoles. Nevertheless, the most known biological function of members of the phy lum Bacteroidetes from the gut is the degradation of biopolymers such as polysaccharides to produce carbohydrates [START_REF] Thomas | Environmental and Gut Bacteroidetes: The Food Connection[END_REF][START_REF] Johnson | Microbiome and Metabolic Disease: revisiting the Bacterial Phylum Bacteroidetes[END_REF]. Thus, this increase of Bacteroidetes relative abundance can be associated to improved produc tion of carbohydrates and production of energy, leading to an improved fitness for the host. The use of other omic techniques such as transcriptomic would be needed in order to deterrnine the links between the changes in gut flora and modifications of functional capacities.

Conclusion

The aim of this study is to fill lacking data concern ing the ecotoxicity of a commercial form of boron nitride nanotubes toward the aquatic compartment. We used an original approach based on the meas urement of toxicological endpoints including geno toxicity and growth parameters of an amphibian species as well as bacterial communities associated to the host gut. The obtained results indicate an overall biocompatibility of the tested BN mixture toward X. laevis tadpoles. Significant induction of larval growth was shown to be correlated with changes in the gut microbial communities of the host. These changes in host physiology are most probably due to indirect effects of byproducts, especially free boron that could stimulate the mat uration of the gut, benefiting to bacteria favoring the host metabolism. Thus, BNNT alone represent a minor threat for amphibians in aquatic environ ments. However, due to its sorption capacities, interactive effects with other common contaminants such as PAHs, heavy metals or pesticides are remaining to be assessed to fully characterize its ecotoxic potential.

Figure 1 .

 1 Figure 1. High-resolution TEM images of the BNNT sample. (a) overview at low magnification illustrating the fibrous structure. (b-d) High-resolution images showing on (b) the presence of multlwalled nanotubes (from 2 to 10 concentric walls), (c) damaged walls of some nanotubes (black arrow) and two isolated amorphous particles surrounded by some graphitic-like shells (white arrows), (d) some filling inside some nanotubes (dotted white arrow), a broken and open nanotube (black arrow) as well as an amorphous particle which seems to be at the end of a nanotube. The inset (same scale) of (d) also shows two isolated particles surrounded by some graphitlc-like shells, one being amorphous while the other one (top) is crystallized.

7 . 2 /

 72 Bnnt characteristicsBET measurement indicated a specific surface area of 163 m g, a little bit lower than the > 200 m 2 /g claimed by the provider. TEM images (Figure 1 (a)) shows that the sample is composed of nanotubes (2-10 walls, Figure l(b)) and small nanoparticles encapsulated in graphitic-like h-BN shells (size ca.

Figure 2 .

 2 Figure 2. BNNT outer diameter distributions evaluated from the HRTEM images (n = 91 ).

Figure 3 .

 3 Figure 3. Variation of the transmission vs time for suspensions of BNNT at 10 mg L -1 in exposure media measured at the top of the via! (A) and at the bottom of the via! (B).

Figure 4 .

 4 Figure 4. Micronucleus induction (A) and cell cycle (B) measured in erythrocytes of Xenopus tadpoles exposed for 12 days to increasing concentrations of BNNT. MNE: micronucleated erythrocytes; NC: negative control; PC: positive control; *significant differ ence compared to the NC (McGill test). For cell cycle analysis, N = 12, A NOVA p < 0.001 followed by Tukey test. Letters indicate significant differences between concentrations tested for each phase of the cell cycle.

Figure S .

 S Figure S. Normalized growth of X. laevis larvae measured after 12 days of exposure to increasing concentrations of BNNT. ANOVA p < 0.001 followed by Tukey test. Letters indicate significant differences between concentrations.

Figure 6 .

 6 Figure6. Effects of 12 days of exposure to BNNT on the structure of gut bacteri al communities of X. laevis tadpoles. Observed richness, estirnated richness (Chao1 index) and evenness (Shannon index) are compared be tween the exposure conditions (0, 0.1, 1 and 10 mg L-1 ) (A). MDS plot of bacterial communities based on Unweighted unifr ac distances (B).

Figure 7 .

 7 Figure 7. Phylum composing the gut microbiota after 12 days of exposure to BNNT concentrations ranging from O to 10 mg L _,.

Acknowledgments

The authors thank BNNT, LLC for providing us with the material used in this study. J. Esvan is acknowledged for his help with the acquisition and analysis of XPS data.

Disclosure statement

No potential confl ict of interest was reported by the author(s).

Ecotoxicological assessment of commercial Boron Nitride Nanotubes towards Xenopus laevis tadpoles and host-associated gut microbiota

Supplementary information

Lauris Evariste 1 , Emmanuel Flahaut 2 , Clément Baratange 1 , Maialen Barret 1 , Florence Mouchet 1 , Eric Pinelli 1 , Anne Marie Galibert 2 , Brigitte Soula 2 , Laury Gauthier 1 1 EcoLab, Université de Toulouse, CNRS, INPT, UPS, Toulouse, France

Detailed BNNT characteristics

XRD pattern showed 2 main peaks at ca. 2 = 25.9° and 42.4°. The first one may correspond to the (002) line of h-BN (although the interlayer spacing is larger than expected at 3.44 Å instead of 3.36 Å), while the second one may be a superimposition of the (100) and (101) lines. XPS analysis (Figure S1 and Table S1) revealed 40.2 at.% of Boron involved in B-N bond, and 39.9 at.% of Nitrogen involved in the same bond, thus BN represents at least 80 at. % of the total. It is however impossible from the XPS analysis to make the difference between, BN involved in nanotubes and BN particles. The only exogenous element evidenced by XPS analysis was Carbon, which is likely to come from contamination during sample handling and preparation. 4.4 at.% of Boron would be involved in O-B-N bonds and 1.6 at. % would be involved in direct B-O bonds. This confirms that the sample contains a large majority of BN but that some boron oxide is also likely to be present, while pure elemental boron seems excluded (no B-B bond identified). TGA analysis (Figure S2) indicated a good thermal stability of the BNNTs in air until 425°C, temperature at which some oxidation started to occur in 2 steps. The first step was observed between 425°C and 825°C, while the second one, much more pronounced, began at 825°C and stopped only close to 990°C. The first oxidation step was earlier attributed to the oxidation of elemental boron [START_REF] Harrison | Quantification of hexagonal boron nitride impurities in boron nitride nanotubes via FTIR spectroscopy[END_REF] but should have a different origin here as elemental boron was ruled out from the XPS analysis. The comparison with purified BNNTs (prepared by air oxidation treatment followed by washing in boiling water, for characterization purpose only) revealed that the first oxidation step was absent in the purified material, indicating that it may be related either to BN not included in BNNTs, or to boron compounds (oxide, boric acid) present in the raw material. Raman analysis at 532 nm (Figure S3) revealed mainly a single peak at ca. 1370 cm -1 attributed to the E2g mode related to in-plane vibration of the B-N bond. However, it seems that the same peak may be observed both for BNNTs and h-BN [START_REF] Harrison | Quantification of hexagonal boron nitride impurities in boron nitride nanotubes via FTIR spectroscopy[END_REF]. XRD analysis revealed mostly a main broad peak just below 26° corresponding to the (002) line in BNNTs [START_REF] Harrison | Quantification of hexagonal boron nitride impurities in boron nitride nanotubes via FTIR spectroscopy[END_REF]2], while the characteristic peak of h-BN for the same line was obviously absent from the BNNT sample. This analysis also ruled out the presence of detectable amounts of crystallized boric acid or boron oxide. The other peaks cannot be used to make a clear difference between h-BN and BNNTs.

Figure S3: Raman analysis of BNNTs (532 nm).

Finally, from IR analysis (Figure S4), it appears that this is not easy to clearly differentiate between h-BN and BNNTs (although H. Harrison et al have recently proposed a method based on FTIR spectroscopy to do so [START_REF] Harrison | Quantification of hexagonal boron nitride impurities in boron nitride nanotubes via FTIR spectroscopy[END_REF]. For BNNTs, a longitudinal mode is expected at 1369 cm-1, while a tangential circumferential one is expected at 1545 cm -1 [START_REF] Lee | Effective growth of boron nitride nanotubes by thermal chemical vapor deposition[END_REF]. These 2 peaks are observed but are very broad and strongly overlap in our sample. In h-BN, a transverse optical mode and a longitudinal mode both resonate near 1350 cm -1 , which are not clearly resolved in our sample. However, the peak at 1540 cm -1 was also experimentally reported for h-BN [START_REF] Geick | Normal Modes in Hexagonal Boron Nitride[END_REF]. We observe a narrow peak at 811 cm -1 , attributed to an out-of-plane vibration (Radial) vibration of B-N bonds. Peaks in the 970 -1170 cm -1 range are often attributed to oxidized forms such as boron oxide, boric acid or metaboric acid [5] but are observed at different places in our sample (925, 1032, 1102 cm-1) and exhibit a very low intensity.