Augustin Ernoult
email: augustin.ernoult@inria.fr

Juliette Chabassier

BORE RECONSTRUCTION OF WOODWIND INSTRUMENTS USING THE FULL WAVEFORM INVERSION

Several techniques can be used to reconstruct the internal geometry of a wind instrument from acoustics measurements. In this study, the passive linear acoustic response of the instrument is simulated and an optimization process is used to fit the simulation to the measurements. This technique can be seen as a first step toward the design of wind instruments, where the targeted acoustics properties come no more longer from measurements but are imposed by the designer. The difficulties of this approach are to find the best acoustic observation allowing the reconstruction (impedance, reflection function, etc.) but also to have an efficient optimization process. The "full waveform inversion" (FWI) is a technique coming from the seismology community. It uses the knowledge of the equation modeling the wave propagation into the instrument (here the telegraphist equation) to have an explicit expression of the gradient of the function which is minimized. This gradient is evaluated with a low computational cost. The FWI methodology, along with 1D spectral finite element discretization in space, applied to the woodwind instruments (with tone holes, losses and radiation) is presented in this communication. The results obtained for the bore reconstruction with different acoustics observations are then compared and discussed.

INTRODUCTION

The aim of the bore reconstruction is to deduce the internal geometry of a wind instrument (bore profile) from acoustic measurement (typically input impedance measurement), giving a non-destructive tool to estimate the geometry of existing instrument.

Several previous studies conducted bore reconstruction of brass instruments, in which the resonator part is composed of one unique main tube, but to our knowledge, no reconstruction of woodwind instrument (with side holes) has been conducted. One approach is the pulse reflectometry technique [START_REF] Sharp | Acoustic pulse reflectometry for the measurement of musical wind instruments[END_REF][START_REF] Amir | A Discrete Model for Tubular Acoustic Systems with Varying Cross Section -The Direct and Inverse Problems. Part 1: Theory[END_REF], a signal approach in which the main tube is supposed to be composed of small conical parts, the length of which are related to highest frequency of the signal. However this technique can not be applied directly to woodwind instruments in which the resonator can be seen as a "pipes network" due to the side holes. A second possibility is to use optimization algorithms to inverse the impedance computation: find the geometry which minimizes the residual between measured and simulated impedance. This approach has been used by W.Kausel [START_REF] Kausel | Optimization of Brasswind Instruments and its Application in Bore Reconstruction[END_REF] on brass instruments by using a zero order optimization algorithm (Rosenbrock). This approach can be seen as a first step toward the design of wind instruments, where the targeted acoustics properties come no more longer from measurements but are imposed by the designer. Even if this principle can be applied to woodwind instruments, due to the complexity of the resonator geometry, it is necessary to use a more efficient optimization algorithm like quasi-Newton algorithm [START_REF] Nocedal | Numerical optimization[END_REF]Chap.3] for which the gradient of the objective function must be known.

The "full waveform inversion" (FWI) is a technique coming from the seismology community [START_REF] Virieux | An overview of fullwaveform inversion in exploration geophysics[END_REF]. It uses the knowledge of the equation modeling the wave propagation to have an explicit expression of the gradient of the acoustic fields inside the propagation medium. It gives the possibility to evaluate, with a low computational cost, the gradient of the objective function of the optimization problem as long as its relation to the acoustic field is differentiable.

The general principle of the FWI and its adaptation to woodwind instrument reconstruction are first presented. The reconstruction of an illustrative "instrument" is then discussed: a cylindrical pipe with four side holes.

All the results presented in this paper are obtained using the OpenWInD python GPLv3 toolbox, the latest version of which can be downloaded at openwind. gitlabpages.inria.fr.

INVERSE PROBLEM

Wave propagation model

To solve the wave equation in the resonator of the instrument and compute its input impedance, it is described as a pipes network. The main bore profile can be described by one or several pipes with any shape (with axial symmetry). At the side holes location, a T-joint junction is assumed between the two pipes of the main bore and the cylindrical pipe corresponding to the chimney hole. The design parameters necessary to describe this geometry are grouped in the vector m. The physical coefficients of the equations may therefore depend on m.

In each pipe n the scaled telegraphist equation is as-sumed:

     L n (m)Z v (x, m, ω)u n + dp n dx = 0 L n (m)Y t (x, m, ω)p n + du n dx = 0 ∀x ∈ [0, 1].
(1) where L n is the pipe length, p n is the acoustic pressure and u n the acoustic flow. The coefficients Z v and Y t account for thermo-viscous effects and depend on the geometry:

       Z v (x, m, ω) = jωρ S(x) [1 -J (k v (m, ω)R(x))] -1 Y t (x, m, ω) = jωS(x) ρc 2 [1 + (γ -1)J (k t (m, ω)R(x))] k v (m, ω) = jω ρ µ , k t (m, ω) = jωρ C p κ
where R is the section radius, S = πR 2 is the section area, Table 1 describes the air constants and the function J of a complex variable is defined as

J (z) = 2 z J 1 (z) J 0 (z) , ∀z ∈ C (3)
where J 0 and J 1 are the Bessel functions of the first kind [START_REF] Tijdeman | On the propagation of sound waves in cylindrical tubes[END_REF].

Air density (kg.m -3 )

ρ = 1.2929(T 0 /(T + T 0 )) Sound celerity (m.s -1 ) c = 331.45 (T + T 0 )/T 0 Viscosity (kg.m -1 .s -1 ) µ = 1.708e -5 (1 + 2.9e -3 T ) Therm. cond. (Cal/(m.s. • C)) κ = 5.77e -3 (1 + 3.3e -3 T ) Specific heat (Cal/(kg. • C)) Cp = 240 Ratio of spec. heat γ = 1.402
Table 1. Physical constant values [START_REF] Chaigne | Acoustics of Musical Instruments. Modern Acoustics and Signal Processing[END_REF]. T is the temperature in Celsius and T 0 = 273.15K. In this study T = 20 • C.

The boundary conditions of each pipe depend on its situation in the global instrument. At the entry (first point of the first pipe), a unitary input flow is imposed : u 1 (x = 0, ω) = 1. At the junctions between two main pipes (indexed by a and b) and one side hole (indexed by c), the entry flow is conserved and the entry pressures of each pipe follow a coupled equation [START_REF] Chaigne | Acoustics of Musical Instruments. Modern Acoustics and Signal Processing[END_REF]Chap.7]:

v a + v b + v c = 0, p a -p c p b -p c = m 11 m 12 m 21 m 22 jω v a v b
where the efficient masses can be obtained from measurements [START_REF] Chaigne | Acoustics of Musical Instruments. Modern Acoustics and Signal Processing[END_REF].

At each pipe end, a radiation impedance is imposed [START_REF] Rabiner | Digital Processing of speech signals[END_REF]:

(α + jωβ) p n (1, ω) = jω Z c v n (1, ω) (4)
where α and β are chosen to describe various radiation models, and can be put to (0, 0) to model a closed hole.

For each pipe the equations are discretized by the Finite Element Method (FEM) leading to a linear matrix equation [START_REF] Tournemenne | A comparison of a one-dimensional finite element method and the transfer matrix method for the computation of wind music instrument impedance[END_REF]. The junction and radiation equations are appended in order to close the discrete system. This procedure results in a global matricial formulation for each considered frequency f = ω/2π:

A tot (m, ω)U tot = E tot . (5)
where all the geometric coefficients contribute to the matrix A tot , the vector U tot contains the acoustic flow and pressure at each degree of freedom (dof) for all pipes and the vector E tot the source term (associated with the input boundary condition).

The system is then solved by inverting the matrix A tot . The input impedance is then the pressure at the dof corresponding to the entry of the instrument:

Z = RU tot ( 6
)
where R is the projection matrix to the corresponding dof.

Choice of observable

The aim of the reconstruction is to determine the value of the parameters in the design vector m by minimizing the residuals between an experimental observation and a simulated one. In our case, the measured data being the impedance, the observation must be a function of the impedance O(Z). To allow the use of quasi-Newton algorithm, this function must be C 1 . The identity, the modulus (as long as it is far from 0) and the phase are good candidates. The reflection function

R = Z -1 Z + 1 (7)
and its modulus are also interesting. The phase of the reflection function is wrapped and can not be used directly [START_REF] Ernoult | Woodwind instrument design optimization based on impedance characteristics with geometric constraints[END_REF]. The cost function J minimized by the optimization algorithm is:

J(m) = ω 1 2 |O(Z simu (ω)) -O(Z meas (ω))| 2 (8)
with Z simu and Z meas , the simulated and measured impedances respectively.

To perform correctly the reconstruction, the chosen observation must ideally be such that the cost function only has one local minimum in the design space at the point corresponding to the actual geometry. The evolution of the cost functions relying of previously evoked observations along one axis of the design space (the position of the first hole) are represented on Figure 1, for a specific scenario. In this case, the reflection function appears to be the observation with the least local minima (and the impedance phase to a lesser extent). It is the observation chosen for this study.

Gradient computation

To perform the optimization with a quasi-Newton algorithm it is necessary to estimate the gradient of the cost function with respect to each design variable m i . It is for this step that the approach used for the Full Waveform Inversion in geophysics is particularly interesting. The Frechet derivative is the mathematical framework used here. In some cases the adjoint state method can also be used efficiently [START_REF] Plessix | A review of the adjoint-state method for computing the gradient of a functional with geophysical applications[END_REF].

The gradient of the cost function can easily be deduced from the differentiation of the acoustic fields with respect to the design variables ∂U tot /∂m i , by formally using composed derivations with equations ( 6), ( 7) and [START_REF] Rabiner | Digital Processing of speech signals[END_REF].

The differentiation of the direct problem of Eq.( 5) yields that the differentiation of the acoustic fields can be obtained by solving a similar problem with another source term:

A tot ∂U tot ∂m i = - ∂A tot ∂m i U tot . (9)
This problem induces a low computational cost because the matrix A tot is already to be inverted for the direct problem.

The most technical part is maybe the computation of the derivative of A tot which necessitates to differentiate the coefficients Z v and Y t .

RECONSTRUCTION OF AN ILLUSTRATIVE CASE

This reconstruction approach is here applied to a simple "woodwind-like" tube, the geometry of which is well known. The simple "woodwind-like" tube is a cylinder (treated here as a specific cone) drilled by four holes. The dimensions are indicated in Figure 2. In this study, only the main bore radius at the entrance is supposed known. The aim of the reconstruction is to determine the values of the 14 other parameters: the length and the right-end radius of the main pipe and the locations, the radii and the chimney heights of the four holes (14 unknowns).

Description of the instrument

A preliminary study reveals that the most useful "fingering" for the reconstruction are the ones for which at the most one hole is open. The 5 corresponding input impedances are measured.

Results

The initial geometry is a 5 cm conical pipe with one hole at each centimeter. A strategy is used to have a rough estimation of the geometry before refining the values of the design variables. For this first step, the low frequencies are the most important. Only 10 frequency values evenly spread between 100 Hz and 200 Hz. First of all the "all closed" impedance is used to adjust only the main bore profile (length and right-end radius). The impedance for which only the fourth hole is open is then used to adjust this hole location, and so on to adjust the four holes locations. These first optimizations are very fast (around 20 sec.), but due to the bad dimensions of the holes (radius and chimney heights) these first estimations of the locations are quite bad.

To refine the values, all the fingerings are taken into account and all the design variables are adjusted together. This time more frequencies are considered (100 values evenly spread between 100 Hz and 3000 Hz). This step is longer (around 5 min). The final results and the comparison with the measured geometric values are indicated in the tables 2, 3, 4, 5. Finally, the errors on the reconstruction can come from uncertainties on the input impedance measurement. The relative uncertainties on the input impedance of a cylinder are around 8% with the home-made sensor used for this study and built following the principles described by Gibiat [START_REF] Gibiat | Acoustical impedance measurements by the two-microphone-three-calibration (TMTC) method[END_REF] with five microphones. A preliminary study of the error propagation in the studied geometry gives around 1mm of uncertainty on the locations, 0.4mm on the diameters and 0.5mm on the chimney heights. It appears therefore that the major identified defect of the reconstruction process presented here is the precision of the impedance measurements.

Variable

CONCLUSION

This article presents the adaptation of the Full Waveform Inversion to the bore reconstruction of woodwind instrument (pipe with side holes). The process is applied to the reconstruction of a cylinder pipe with four side holes. The results indicate that the process is able to reconstruct such an instrument but that it is very sensitive to the precision of the measured impedance used for the reconstruction. Performing the inversion directly from the measured pressure at the five microphones instead of using post processed quantities (impedance, reflection function) could reduce these uncertainties by providing redundant data to the algorithm.

Figure 1 .

 1 Figure 1. Variations of the residuals computed with different observations on the fingering where only the first hole is open, with respect to the relative location of the first hole. The reflection function has the smoother evolution.

Figure 2 .

 2 Figure 2. Sketch of the cylindrical tube with 4 side holes studied. It is described by 15 geometric parameters (3 for the main bore profile and 3 for each of the four holes).

Table 2 .

 2 Main bore dimensions (in mm): values obtained by geometric measurements and reconstruction.

		Measured Reconstr. Error
	Total length	287.5 ± 0.5	287.57	0.07
	Right-end diam.	4 ± 0.01	3.92	0.08
	Variable	Measured Reconstr. Error
	Hole 1 loc. 100 ± 0.5	100.50	0.5
	Hole 2 loc. 130 ± 0.5	130.46	0.46
	Hole 3 loc. 180 ± 0.5	180.12	0.12
	Hole 4 loc. 240 ± 0.5	240.18	0.18

Table 3 .

 3 Holes locations (in mm): values obtained by geometric measurements and reconstruction.It is important to notice that the measured locations have uncertainties around 0.5mm and the other dimensions around 0.01mm. The values obtained by the reconstruction process presented here are quite good for most of the design variables. At low frequency, the diameter and the chimney height of a hole have similar effect. The optimization algorithm can have difficulties to adjust together

	Variable	Measured Reconstr. Error
	Hole 1 diam.	3 ± 0.01	2.95	0.05
	Hole 2 diam. 3.5 ± 0.01	3.54	0.04
	Hole 3 diam. 3.5 ± 0.01	3.26	0.24
	Hole 4 diam. 2.5 ± 0.01	2.19	0.31

Table 4 .

 4 Holes diameters (in mm): values obtained by geometric measurements and reconstruction.

	Variable	Measured Reconstr. Error
	Chimney 1 1.7 ± 0.01	1.74	0.04
	Chimney 2 1.3 ± 0.01	1.46	0.16
	Chimney 3 1.5 ± 0.01	1.52	0.02
	Chimney 4 1.4 ± 0.01	1.14	0.26

Table 5 .

 5 Holes chimney heights (in mm): values obtained by geometric measurements and reconstruction. these two parameters. It can explain the "wide" deviation between measured and reconstructed values for the dimension of the hole 4. Using higher frequencies could improve the reconstruction. The radiation model could also be a source of uncertainties in the reconstruction, since a model error should impact the way physical coefficients influence the simulated impedances, especially close to the holes open end.


Daouda Ndiaye

Dr Maryame Sy

Agnès Pallier

Dr Sandra Même

Isidro De Silva

Sara Lacerda

Dr Aline M Nonat

Loic Charbonniere

Éva Tóth

Loïc J Charbonnière

A Ndiaye

Dr S Pallier

Dr É Lacerda

Tóth

Da Silva

Unprecedented kinetic inertness for a Mn2+-bispidine chelate: a novel structural entry for Mn2+-based imaging agents

Keywords: manganese, contrast agent, MRI, 52 Mn, kinetic inertness

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

Introduction

The high diagnostic value of contrast agents makes them indispensable for many Magnetic Resonance Imaging protocols. For 35 years now, Gd-complexes have been used in millions of human examinations and considered among the safest diagnostic drugs. However, the recent emergence of nephrogenic systemic fibrosis and its causal link to Gd-exposure, as well as the evidence on brain and bone accumulation of Gd have alerted the medical community. [1] In response to these safety concerns, some of the linear Gd 3+ chelates have been recently withdrawn and the use of others has been restricted. In this context, the replacement of Gd 3+ with more biocompatible, safer paramagnetic metal ions has become a major objective. Being an essential metal ion, as well as a good relaxation agent due to its five unpaired electrons (in the high spin state), slow electron spin relaxation and fast water exchange, Mn 2+ is an obvious alternative. [2] Indeed, recent years have witnessed increasing interest in Mn 2+ complexes as potential MRI agents. [3] In addition, manganese has a positron-emitting radioisotope, 52 Mn, with interesting decay properties for Positron Emission Tomography (PET) (t1/2 = 5.6 d, max. β + -energy: 575 keV).

Given its low β + decay intensity (29.6%), PET resolution as good as 1.2 mm could be achieved which is particularly important in small animal imaging. [4] The long t1/2 makes 52 Mn particularly adapted to image slow biological processes. Altogether, Mn 2+ is the unique metal ion to offer detection capabilities in both MRI and PET. [5] Thermodynamically stable and kinetically inert complexation of Mn 2+ is indispensable to avoid in vivo release of the free metal ion, which could potentially cause toxicity during MRI where large contrast agent quantities are needed, or lead to off-target signals in PET examinations. In addition, for good MRI efficiency, the complex should contain at least one inner sphere water molecule. However, the lower charge (with respect to Gd 3+ ) and the lack of ligand-field stabilization energy for the high-spin d 5 electron configuration are not favorable to achieve high thermodynamic stability. The highly labile nature of Mn 2+ sets an even more difficult challenge to meet. Generally, ligand rigidity and preorganization promote kinetic inertness, and indeed, the most promising Mn 2+ complexes are based on rigidified linear or macrocyclic chelators (Scheme 1). For instance, the trans-1,2-cyclohexylene backbone and the pyridine function ensure rigidity and are essential structural elements to provide slow dissociation of the MnPyC3A chelate (Scheme 1), a potential Mn 2+ -based agent which has successfully passed to in vivo experiments. [3a, 6] Among macrocyclic chelates, the non-hydrated MnDOTA was reported to have high inertness, [7] while the monohydrated [Mn(cisDO2A)] was found already ~20-fold more labile (dissociation half-lives, t1/2, estimated to 1061 h and 48 h, respectively, at pH 7.4). [8] Recently, derivatives of the pyridine-containing pyclen macrocycle were found to form stable and very inert Mn 2+ chelates. [3c] With t1/2 = 8×10 3 h estimated for pH 7.4, the Mn(PC2A-EA) complex bearing an ethyleneamine pendant arm and proposed as a pH-sensitive agent, has the highest resistance to dissociation ever reported for a monohydrated Mn 2+ complex. Nevertheless, this dissociation half-life remains a million times below the value for GdDOTA (t1/2 =2.7 × 10 9 h, pH 7.4 and 25 °C). In the aim of substantially improving the kinetic inertness of Mn 2+ chelates while ensuring good relaxation properties and water solubility, we are exploring bispidine (3,7diazabicyclo[3.3.1]nonane) derivatives. Bispidine-type ligands represent a versatile platform in coordination chemistry. [9] Pendant arms with varying coordinating functions can be attached in different positions to the bicyclic amine core in order to tune the coordination properties for metal ions of various sizes, coordination numbers or coordination geometries. Depending on the substituents, the ligand can adopt different conformations and configurations. [10] Among those, the chair-chair conformer is best suited for stable metal coordination and thanks to its exceptionally preorganized structure, bispidine derivatives tend to form very inert complexes. Bispidine complexes have been investigated with a large variety of metal ions and for different applications. Concerning Mn 2+ , the very scarce reports involve pyridine derivative bispidine ligands and are restricted to studies in the solid state or in non-aqueous solvents. [11] We hypothesized that the 2,4-pyridyl-disubstituted bispidol derivative L1 bearing one methylene carboxylic acid (Scheme 1) could be a good chelator for Mn 2+ . First, the very rigid bispidine scaffold can endow MnL1 with high kinetic inertness. The ligand is expected to coordinate in a five-dentate manner, involving two pyridine and two bispidine nitrogens as well as the methylene carboxylate. This leaves one coordination site for a hydration water, important for MRI application. The non-coordinating pending carboxylates might additionally create a second sphere contribution and increase proton relaxivity. Finally, L1 can be readily functionalized for biological targeting via the lysine terminal amine. This scaffold has been previously investigated for the complexation of 64/67 Cu 2+ [12] and 68 Ga 3+ . [13] A three dimensional representation of the putative structure of the MnL1 complex is shown in Scheme 2. We evidence here that despite a moderate thermodynamic stability, MnL1 has a resistance to acid catalyzed dissociation and transmetalation which is unprecedented for a Mn 2+ chelate. In addition, the relaxation efficiency of the complex was assessed by a combined 17 O NMR and NMRD study. The complex was injected into mice and its biodistribution and pharmacokinetics were monitored by MRI. Finally, the 52 Mn radiolabelling of L1 was explored and the stability of the radiocomplex was investigated.

Results and Discussion

L1 was synthesized in 4 steps as previously reported. [12]

Thermodynamic stability and kinetic inertness of MnL1

Ligand protonation constants were obtained by pH potentiometry and are reported in Table 1. By analogy to other bispidine ligands, [14] we attributed the first protonation constant to one of the tertiary amine groups, while the other amine of the bispidine skeleton remains likely unprotonated in the chair-chair conformation. [14b] Further protonation steps supposedly occur on the lysine amine (logKH2 = 10.31, typical of the lysine  amine protonation constant) [15] and then on the three carboxylates, while the pyridine nitrogens are not protonated above pH 2. Thermodynamic stability constants of MnL1 and ZnL1 were assessed by pH potentiometry (Table 1). In the case of Mn 2+ , slow complex formation prevented direct titration; therefore, batch samples were prepared (Fig. S1). MnL1 is three orders of magnitude less stable than ZnL1, in accordance with the Irving-Williams series. Figure 1 depicts the species distribution for MnL1 as a function of pH, supported by the pH-dependent relaxivities. Protonated complexes form with both Zn 2+ and Mn 2+ ; the first protonation constants are practically identical to logKH2 of the ligand and correspond to the protonation of the noncoordinating lysine amine. The similarity between the second protonation constant of the ligand and logKMnHL or logKZnHL further supports the attribution of the ligand protonation steps. At physiological pH, this species protonated on the lysine  amine is the only complex present. The stability of MnL1 is relatively modest, as shown by the pMn value of 6.65 calculated for 10 -5 M Mn 2+ and L1 concentration and pH 7.4 (Table 1). Nevertheless, once the chelate is formed, it has an exceptional resistance to dissociation. Kinetic inertness of metal complexes is typically assessed by transmetalation experiments with endogenously available metal ions, such as Zn 2+ , which is present at ~10 M concentration in the blood and represents the most abundant potential competitor of Mn 2+ . Different transmetalation protocols have been reported for Mn 2+ complexes, some allowing a simple comparison of the kinetic inertness under given conditions, others giving also access to the underlying dissociation mechanisms. We followed transmetalation in the presence of 10 or 50 equivalents of Zn 2+ , at pH 6 (50 mM MES buffer), 37°C, in a similar way to that proposed by Caravan et al. [3a] At either of the Zn 2+ concentrations, MnL1 did not show any relaxivity change for the 140 days of observation (Fig. 2). Under similar conditions (37°C, pH 6.0, 25 Zn 2+ equivalents), dissociation half-lives, t1/2, of 0.285 h and 54.4 h have been reported respectively for MnPyC3A [3a] and MnPC2A-EA [3c] , two reference compounds considered as potential MRI probes with good kinetic inertness. Analogous experiments have been carried out at more acidic pH values as well (pH 3.1-6.0, 50 equivalents Zn 2+ , 37 °C; Fig. 2). Under these pseudo-first order conditions, the dissociation rate constant, kobs, has been determined, and the estimated transmetalation halflives vary between 22.7 h at pH 3.13 and 130 days at pH 5.07. The inertness of MnL1 was also assessed in highly acidic solutions where the dissociation is much faster (0.01-1.0 M HCl, 25°C; Fig. 3). The dissociation rate constants, kobs, show a second order dependence on the proton concentration and could be fitted to Eq. 1.

𝑘 𝑜𝑏𝑠 = 𝑘 0 + 𝑘 1 [𝐻 + ] + 𝑘 2 [𝐻 + ] 2 (Eq. 1)
The fit resulted in k1 = (1.6±0.1)×10 -3 s -1 M -1 and k2 = (5.0±0.1)×10 -4 s -1 M -2 , while k0 was fixed to 0, otherwise the fit yielded a negative value with a large error (-3.5±2)×10 -5 s -1 . We should note that the kobs values determined above in the Zn 2+ transmetalation study in the pH 3.1-5.1 interval also show a second order dependence on [H + ] (Fig. S2), with k1 = (2.2±0.2)×10 -3 s -1 M -1 , in the same order of magnitude as the one above. This rate constant, characterizing the proton assisted dissociation pathway, is two orders of magnitude lower than that reported for MnPC2A-EA which is the most inert hydrated Mn 2+ chelate reported so far. [3c] At this point, it is difficult to extrapolate the dissociation half-life of the complex for physiological conditions, and in particular, more data are needed to describe the role of metal-assisted dissociation processes in the overall dissociation. These studies are in progress, but they are long due to the extremely low dissociation rates. Nevertheless, we can conclude that the rigid and preorganized bispidine skeleton provides very high kinetic inertness, unprecedented for a Mn 2+ chelate.

Relaxation properties of MnL1

Proton relaxivities have been recorded for MnL1 as a function of the magnetic field, at 25 and 37°C, in water (pH 7) and in human serum (Figs. 4a andS3). The NMRD data were complemented by variable temperature 17 O relaxation rate measurements (Fig. 4b) which give access to the hydration number and the water exchange parameters. The relaxivity is remarkably high for a Mn 2+ complex (r1 = 4.28 mM -1 s -1 and 3.37 mM -1 s -1 at 25°C and 37°C, respectively, 20 MHz, pH 7 in water). In comparison, relaxivities of 3.3 and 2.5 mM -1 s -1 (20 MHz, 25 and 37 °C) and 3.52 mM -1 s -1 (20 MHz, 25°C) were reported for the monohydrated complexes MnPCy3A and MnPC2A-EA. One hydration water molecule is confirmed for MnL1 from the temperature-dependent 17 O transverse relaxation rates using the method of Gale et al. (q = 0.85; Fig. S4). [17] The temperature-dependent transverse 17 O relaxation rates and the NMRD profiles have been analysed in the frame of the Solomon-Bloembergen-Morgan theory to determine the parameters describing water exchange and rotational dynamics. The fit yielded kex 298 = (5.1±0.7)×10 7 s -1 and H ǂ = 10.6 kJ/mol for the rate and the activation enthalpy of the water exchange. The water exchange rate is similar to that reported for MnPyC3A (kex 298 = 5.4×10 7 s -1 ) [3a] while it is about ten times lower than kex 298 on MnEDTA (kex 298 = 47×10 7 s -1 ). [18] The bispidine MnL1 chelate is structurally very different from other previously studied Mn 2+ complexes, and the reason for this relatively slow water exchange is difficult to identify. Nevertheless, this exchange rate remains high enough not to limit proton relaxivity. The rotational correlation time, R 298 = (100±2) ps, is higher than the values reported for small molecular weight Mn 2+ complexes (e.g. 56 ps for MnEDTA). [19] This might reflect that R is overestimated due to second sphere contribution to relaxivity from the noncoordinating carboxylates. Therefore, we also performed the analysis of the NMRD data by including an additional contribution from a second sphere relaxation mechanism to the overall relaxivity. To describe this second sphere mechanism, we assumed two 2 nd sphere water molecules (with an estimated Mn-H proton distance of 3.2 Å, 50 ps residence time and an activation enthalpy of 10 kJ/mol, based on previous simulations on Gd 3+ complexes; see Fig. S5 and ESI for all details). [20] This resulted in a rotational correlation time, R 298 = (72±4) ps which seems more reasonable on the basis of the molecular weight of the complex. This is obviously not a firm proof of a second sphere relaxivity contribution; nevertheless, it is evident that the relaxivities measured tend to be higher than the typical values reported for monohydrated Mn 2+ chelates. The relaxivities remain identical in the presence of up to 50 equivalents of citrate, carbonate and phosphate (Fig. S6). They increase to 6.46 and 5.12 mM -1 s -1 in human serum at 25 and 37 °C, respectively (20 MHz, Fig. 4a), a tendency similar to that observed for instance for MnPyC3A and attributed to low affinity binding to serum proteins. 52 Mn radiolabeling of L1

The need for long-lived   -emitting isotopes and the emergence of a combined PET-MRI technology have recently promoted an increasing interest in the production and complexation of 52 Mn. We have produced and purified 52 Mn at the cyclotron of the Orléans CNRS campus (see ESI for details). 52 Mn was recovered as 52 MnCl2 solution after purification through 2 anion exchange resin AG © 1-X8 columns. The radiolabelling of L1 was tested under different conditions of pH, temperature and 52 Mn:L1 ratio. The radiolabelling reactions were followed by radioTLC at 5, 15, 30 and 60 min. All reactions were performed in presence of ascorbic acid to prevent any oxidation reaction. Optimized reaction conditions were obtained for pH 7, 70 °C, 1h, achieving 99% labelling yield. 52 MnL1 is highly hydrophilic with a logP = -1.96±0.06, which falls within the range reported for other Mn complexes, e.g. MnDPDP or Mn(EDTA-BTA) with logP values of -3.07 and -1.84, respectively. [21] The in vitro stability of the radiocomplex was assessed in different media at pH 7.4 (water, PBS and saline (0.9% NaCl)) and in the presence of 0. HSA (Fig. 5), revealing its stability in all media up to 18h, and a slightly lower to moderate stability at 24h.

MRI biodistribution and pharmacokinetics of MnL1

To validate the potential in vivo use of MnL1, we have evaluated the biodistribution in wild type mice by MRI at 7 T. The signal intensity of the main organs (kidney, liver, muscle, lung) and the aorta was followed during 1 h post intravenous injection (0.06 mmol/kg dose; Fig. S7). Figure 6 depicts normalized signal intensity for the analyzed organs. The highest signal intensity was reached around 4 min after injection. MnL1 shows renal clearance and the uptake in the liver is limited. The estimated blood half-life is 21 min, typical of small molecules and similar to previously reported contrast agents. [3a, 22] This distribution profile was further confirmed by ex vivo ICP-OES measurements performed at the end of the MRI experiment. The data are compared to mice which did not receive any MnL1 injection (control) where the Mn-levels measured correspond to endogenous Mn (Fig. 7; Table S1). At 1.5 h post injection, higher Mn content is found in the kidneys and blood of injected mice. The blood uptake could be related to the low affinity binding of MnL1 to blood proteins, as also observed by relaxivity measurements. Compared to the Mn tissue content reported for MnPyC3A at 24 h post injection, our complex shows lower liver and kidney uptake. Data is presented in nmol/g tissue ± SD (n=3).

Conclusion

We demonstrate here that despite a relatively modest thermodynamic stability, the MnL1 complex has an exceptional resistance to dissociation. At 37°C, pH 6.0 and in the presence of up to 50 equivalents of Zn 2+ , MnL1 remained intact and did not show any relaxivity change for at least 140 days. MnL1 has one inner sphere water molecule. It has higher proton relaxivities than most monohydrated, small molecular weight Mn 2+ chelates, which might be related to a second sphere relaxivity contribution induced by the non-coordinating carboxylates functions.

In vivo MRI experiments performed in mice indicated quick renal clearance, which was also supported by ex vivo ICP determination of the Mn tissue content in different organs. L1 could be successfully labeled with 52 Mn and the radiocomplex has good stability in biological media. Based on these preliminary radiolabeling experiments, the conjugation of the bifunctional ligand L1 to biomolecules, such as antibodies can be envisaged. Indeed, for the imaging of such slowly circulating biomolecules, the long-lived 52 Mn PET isotope represents a clear advantage. In summary, the bispidine-based L1 chelator constitutes a very promising structural entry for the development of Mn-based imaging agents for MRI as well as for PET. Alternatively, this bispidine synthon can be also attractive for the design of combined PET-MRI probes. Such bimodal imaging applications based on Mn 2+ might be interesting when the in vivo fate of the imaging probe is to be followed for a longer period of time. Most importantly, the L1 ligand can provide excellent kinetic inertness which so far could not be achieved within the more traditional ligand families.

Experimental Section

Experimental details are described in detail in the Supporting Information.

Scheme 1 .

 1 Scheme 1. Ligands discussed.

Scheme 2 .

 2 Scheme 2. Three dimensional representation of the putative structure of MnL1.

4 Figure 1 .

 41 Figure 1. Species distribution curves calculated for MnL1 (1 mM) and pHdependent relaxivities (▲) measured at 25°C, 60 MHz.

Figure 2 .Figure 3 .

 23 Figure 2. Time-dependent variation of the longitudinal proton relaxation rates in a 1 mM MnL1 solution at 37°C and 60 MHz, 0.1 M KCl, in the presence of 50 equivalents of Zn 2+ . Empty symbols at pH 6.00 correspond to 10 equivalents of Zn 2+ . The curves represent the fit of the experimental data to yield the observed rate constants, kobs.

Figure 4 .

 4 Figure 4. a) Proton relaxivities of MnL1 measured in water at 25°C (♦) and 37 °C (▲) and in human serum at 37 °C (∆), and b) reduced transverse 17 O relaxation rates (9.4 T). The lines represent the simultaneous fit of the experimental points as explained in the text.

Figure 5 .

 5 Figure 5. In vitro stability of 52 MnL1 in different media at pH 7.4 (water, PBS and saline (0.9% NaCl)) and in the presence of 0.6 mM HSA, at different time points: 0, 1, 18 and 24 h (darker to lighter colour, respectively).

  [3a]

Figure 6 .

 6 Figure 6. a) Normalized signal intensity in the kidney, muscle, liver and aorta plotted as a function of time. Measurements were performed every 2 minutes for 4 mice. Standard deviations are not presented for better readability. b) Sagittal T1-weighted MR images prior to injection and at 4, 15 min and 1 h post intravenous injection of 0.06 mmol/kg of MnL1.

Figure 7 .

 7 Figure 7. Ex vivo ICP-OES quantification of the Mn tissue content in the main organs and blood of control and injected mice with MnL1 (1.5 h post injection).

Table 1 .

 1 Ligand protonation constants, stability constants of ML complexes and pMn values. Values in parenthesis correspond to one standard deviation.

		L1	PyC3A [3a]	trans-CDTA [16]	PC2A-EA [6]
	logKH1	11.44(1)	10.16	9.54	11.34
	logKH2	10.31(2)	6.39	5.97	8.93
	logKH3	4.71(5)	3.13	3.60	6.91
	logKH4	2.76(5)	-	2.52	1.97
	logKH5	2.22(4)	-	1.46	-

Acknowledgements

The authors thank the French National Research Agency (grant ANR-18-CE18-0008) for funding. We thank Dr. Frederic Szeremeta and Quentin Mura for their assistance with MRI experiments and Louis Frealle with the 52 Mn production.

Entry for the Table of Contents

The 2,4-pyridyl-disubstituted bispidol ligand L1 is a promising structural entry for the development of Mn-based imaging agents. Thanks to its preorganized, rigid structure, it provides exceptional inertness to the Mn 2+ complex. MnL1 is monohydrated, has good relaxation efficiency for MRI applications and undergoes rapid renal clearance when injected to mice. The 52 Mn-labeled radiocomplex could be obtained and was stable in biological media.


Gaétane Le Provost
email: gaetaneleprovost@free.fr

Isabelle Badenhausser

Yoann Le Bagousse-Pinguet

Yann Clough

Laura Henckel

Cyrille Violle

Vincent Bretagnolle

Marilyn Roncoroni

Peter Manning

Nicolas Gross

Gaëtane Le Provost

Gaëtane Le Provost

Land-use history impacts functional diversity across multiple trophic groups

Keywords:

most likely because they offer a stable and high quality habitat refuge for species with low mobility and specialized feeding niches. Our study shows that grassland-to-crop conversion has lasting impacts on the functional biodiversity of agricultural ecosystems. Accordingly, land-use legacy effects must be considered in conservation programs aiming to protect agricultural biodiversity. In particular, the retention of permanent grassland sanctuaries within intensive landscapes may offset ecological debts.

SIGNIFICANCE

The World's agricultural landscapes have seen dramatic declines in their biodiversity. Because biodiversity responses to land-use changes can be delayed, understanding how current biodiversity is impacted by historical land-use changes is crucial to protect and restore agricultural biodiversity.

By utilizing a long-term land-use record database, we show that historical land-use changes continue to shape current biodiversity. Our study expands our understanding of global biodiversity change by identifying land-use history effects to be as strong as the effects of current management and the amount of semi-natural land cover, often assumed to be the major drivers of biodiversity change in agricultural landscapes. Accordingly, land-use legacy effects must be considered in management programs that aim to protect and restore biodiversity.

INTRODUCTION

Habitat destruction caused by land-use change is a major driver of global biodiversity declines [START_REF] Newbold | Has land use pushed terrestrial biodiversity beyond the planetary boundary? A global assessment[END_REF].

While land-use changes often have time-delayed impacts on biodiversity [START_REF] Tilman | Habitat destruction and the extinction debt[END_REF][START_REF] Krauss | Habitat fragmentation causes immediate and time-delayed biodiversity loss at different trophic levels[END_REF], we know very little about the relationship between land-use history and present-day biodiversity. The presence of remnant and alternative habitats in the landscape can buffer biodiversity loss through time, if organisms are able to persist in, and disperse to these refuges (4). These meta-population processes can delay or prevent species extinction for years or decades [START_REF] Blitzer | Spillover of functionally important organisms between managed and natural habitats[END_REF], meaning we may have underestimated the impacts of past land-use change on biodiversity. Understanding the long-term response of biodiversity to land-use change is crucial if we are to limit species extinctions and biodiversity decline by implementing sound and timely conservation and restoration efforts.

Species can differ significantly in their responses to land-use changes [START_REF] Allan | Interannual variation in land-use intensity enhances grassland multidiversity[END_REF], making it hard to generalize biodiversity response as a whole, and to predict which species and taxa may suffer the most from habitat destruction in the long run. A solution to this problem is the use of species traits, as species sharing similar traits have been shown to respond consistently to land use change [START_REF] Gámez-Virués | Landscape simplification filters species traits and drives biotic homogenization[END_REF][START_REF] Flynn | Loss of functional diversity under land use intensification across multiple taxa[END_REF]. Furthermore, in agricultural landscapes, trait diversity plays an important role in determining the capacity of ecosystems to cope with future environmental changes [START_REF] Vogel | Chapter Three -Lost in trait space: species-poor communities are inflexible in properties that drive ecosystem functioning[END_REF] and secure the provision of key ecosystem services. However the long term response of functional diversity to land-use changes is unknown as previous studies investigating the effects of land-use history have focused on species richness [START_REF] Krauss | Habitat fragmentation causes immediate and time-delayed biodiversity loss at different trophic levels[END_REF][START_REF] Kuussaari | Extinction debt: a challenge for biodiversity conservation[END_REF], ignoring trait diversity. A general understanding of how the trait diversity of multiple trophic groups responds to land-use changes over time may therefore improve our capacity to manage agricultural biodiversity, and secure sustainable long-term agroecosystem functioning.

Here, we investigated how 20-years of past land-use changes affected the functional trait diversity of seven taxonomic groups in grasslands (vascular plants, wild bees, hoverflies, grasshoppers, carabid beetles, spiders and birds), belonging to five trophic groups (primary producers, pollinators, herbivores, invertebrate and vertebrate predators) (see for details SI Appendix, Table S2). We used a unique database of 220,000 land-use records in an agricultural region of 430 km² dominated by crop production in western France (11,000 fields recorded annually since 1994, ref. 10) to investigate the effect of multiple aspects of land-use history, operating at the field and landscape level, on functional diversity (SI Appendix, Fig. S1 and Table S1). Historically, the study region was a typical rural area composed of mixed crop-livestock systems. Fifty years ago, grassland was the dominant land use, covering 60% of the study area [START_REF] Bretagnolle | Towards sustainable and multifunctional agriculture in farmland landscapes: Lessons from the integrative approach of a French LTSER platform[END_REF]. These grasslands were either grazed or mown for livestock production. Since that time and up to the present day, shifts from livestock to annual crop production resulted in a strong decline in grassland cover and in 2014, grasslands covered only 12% of the area.

Within the study area, we sampled 75 grasslands with a wide range of land-use histories, while controlling for the effects of current landscape composition and the management of the sampled grassland field. Land-use history was assessed by the age of the sampled grassland field ("Field age") and by three independent landscape metrics (SI Appendix, Fig. S1 and Table S1): (i) the time elapsed since the first grassland-to-crop conversion for all fields within a landscape of 1km radius of the sampled grassland ("Time as cropland"). High value of time as cropland indicates that most grasslands in the landscape were converted into crops long ago; (ii) the permanency of the grassland cover in the landscape ("Grassland permanency"). A high value of grassland permanency corresponds to landscapes with old permanent grasslands; (iii) the turnover of rotation from crop to grassland in the landscape ("Crop-grassland turnover rate"). A high value of crop-grassland turnover rate indicates that sown grasslands are maintained for several years. These three metrics complement each other to inform on the loss, the stability and turnover of grasslands in the landscape.

We focused on a core set of independent organismal traits that mediate species responses to land-use changes (Fig. S2): body size, mobility (e.g. wing span) and resource-use traits (e.g. mandible strength or beak size). Such traits can directly determine a species' ability to respond rapidly to land-use intensification, e.g. frequent disturbances and resource homogenization in agricultural landscape may for species with smaller size, a generalist diet, shorter generation time and higher dispersal abilities [START_REF] Gámez-Virués | Landscape simplification filters species traits and drives biotic homogenization[END_REF][START_REF] Brown | Toward a metabolic theory of ecology[END_REF]. We summarized the functional diversity of the entire trophic chain using multitrophic and multi-trait diversity metrics [START_REF] Kuussaari | Extinction debt: a challenge for biodiversity conservation[END_REF]. We also investigated how land use affects different facets of functional trait diversity by compiling three multitrophic trait diversity indices for mobility, resource-acquisition and body size traits separately. We controlled for confounding effects of the local species pool (the multitrophic density and species richness, following ref. 9) in our analyses, to isolate net effects of land-use changes on functional diversity (see methods). We hypothesized that landscapes in which grassland was converted to crop long ago would support low grassland multitrophic trait diversity by selecting generalist species and reducing the diversity of body size and mobility traits. We further hypothesized that increasing grassland permanency and decreasing the crop-grassland turnover rate could mitigate biodiversity loss by providing stable and high quality habitat to low mobility organisms and feeding specialist species.

RESULTS

All models accounted for the influence of multitrophic species richness and density on trait diversity (Fig. 1). Multitrophic trait diversity increased linearly with multitrophic species richness and decreased with multitrophic density (Fig. 2), with the exception of resource-acquisition trait diversity for which multitrophic density had no effect. We also found multitrophic species richness to be more strongly affected by land-use history (43% of explained variance) than by current land use (18% of explained variance) (SI Appendix, Figs. S4 andS5). On the contrary, multitrophic density was mostly driven by current land use (62% vs. 38% of variance explained by land-use history).

The best-selected model explained 70% of the variance in multi-trait diversity (Fig. 1; see also SI Appendix, Table S5). Of this, current land use accounted for 16% of the explained variance (Fig. 1), where multi-trait diversity was higher in landscapes with high grassland cover (Fig. 2).

Land-use history explained 12% of multi-trait diversity, through the effects of grassland field age (5% of the explained variance) and grassland permanency in the landscape (accounted for 7% of the explained variance) (Fig. 2). These effects were positive: multi-trait diversity was higher in fields that had been grasslands for a long time and were surrounded by other permanent grasslands (Fig

2).

When considering each trait separately, we found that mobility, resource-acquisition and body size trait diversity responded differently to past and current land use. Land-use history was the main driver of the diversity of mobility traits (34% of explained variance) (Fig. 1). Mobility trait diversity was lower in fields embedded in landscapes converted to cropland long-ago (Fig. 2). In contrast, the diversity of mobility traits was higher in old grassland fields and in landscapes in which grassland cover was more permanent. Past and current land use accounted for an equal amount (28%) of explained variance for resource-acquisition trait diversity (Fig. 1). The diversity of resource-use strategies initially increased following grassland-to-cropland conversion in the landscape (during the first 10 years) but subsequently decreased (Fig. 2). Current land use could either mitigate this loss, as resource-acquisition trait diversity increased with present-day grassland cover, or accentuate it, as it decreased with grassland field productivity, a proxy of grassland management intensity, and the amount of mass flowering crops in the landscape. Our model explained less variation in body size diversity, although this did respond to current land use (16% of explained variance) (SI Appendix, Fig. S6).

When considering each trophic group separately, plant trait diversity was only impacted by current grassland field productivity (27% of explained variance). For animals, the importance of land-use history decreased with increasing trophic level, while the importance of current landscape-scale predictors increased (Figs. 3 and4, see also SI Appendix, Table S6). Herbivore (i.e. grasshoppers) and mutualist (i.e. wild bees and hoverflies) trait diversity was strongly affected by land-use history (respectively 39% and 27% of explained variance) and was higher in old grassland fields (Fig. 4). Land-use history accounted for 28% of the explained variance in invertebrate predator trait diversity, which increased in landscapes dominated by permanent grasslands.

Vertebrate predator trait diversity was mainly driven by current landscape composition, and increased in landscapes with a high cover of mass flowering crops (16% of explained variance) and forest (46% of explained variance).

DISCUSSION

We investigated how land-use history affects the functional diversity of multiple trophic groups in grasslands. We found that historical land use continues to play an important role in shaping the present-day functional diversity across multiple trophic groups, and that it is of equal importance to well established land-use drivers of biodiversity in agricultural landscapes, such as grassland management intensity [START_REF] Flynn | Loss of functional diversity under land use intensification across multiple taxa[END_REF] and current landscape composition [START_REF] Gámez-Virués | Landscape simplification filters species traits and drives biotic homogenization[END_REF]. Ignoring the role of land-use legacies on functional diversity may therefore hinder our ability to predict how biodiversity responds to land-use changes, and the functional consequences of biodiversity change on ecosystems.

Multi-trait diversity was most strongly influenced by the species pool, i.e. the density of individuals and the taxonomic diversity observed within the sampled grasslands. There was a negative relationship between multi-trait diversity and density. High population densities, together with low functional diversity, occur when a small number of functionally similar species benefit from land-use intensification in agricultural landscapes. These few species may benefit from the resources offered by particular crop monocultures and become hyper abundant, e.g. in pest outbreaks [START_REF] Hillebrand | Consequences of Dominance: A Review of Evenness Effects on Local and Regional Ecosystem Processes[END_REF]. However, by increasing functional diversity, we may limit dominance by a few species with potentially negative effects on crops, as suggested by the dilution effect hypothesis [START_REF] Bebber | The global spread of crop pests and pathogens[END_REF][START_REF] Schmidt | Biodiversity and the Dilution Effect in Disease Ecology[END_REF]. We also found that multi-trait diversity strongly increased with taxonomic diversity (Fig. 2), which highlights low functional redundancy among species across all the trophic groups studied. A low functional redundancy implies that species loss may be accompanied by strong declines in functional diversity and uncertain outcomes in agroecosystem functioning. While agricultural landscapes have the potential to host high trait diversity [START_REF] Provost | Trait-matching and mass effect determine the functional response of herbivore communities to land-use intensification[END_REF], the low functional redundancy observed in our study suggests that the functional biodiversity of agroecosystems is particularly vulnerable to land-use changes and species loss [START_REF] Laliberte | Land-use intensification reduces functional redundancy and response diversity in plant communities[END_REF].

Decomposing multi-trait diversity into three independent indices (Fig. S2) allowed us to investigate the mechanisms through which land-use history alters biodiversity in agricultural landscapes. The multitrophic diversity of mobility and resource-acquisition traits was lower in landscapes that were converted from grassland to cropland long ago (Fig. 2), highlighting that habitat destruction has not only immediate, but also long-lasting negative effects on the functional diversity of agroecosystems. The conversion of permanent grasslands into annual crops may have particularly strong effects on low mobility organisms and species with narrow feeding niches.

However, our results show that maintaining permanent grasslands in the surrounding landscape mitigates this loss (Fig. 2), probably by providing stable, heterogeneous and resource-rich habitats in which many species can persist in and migrate between. While increasing the quantity of seminatural habitats in the landscape is often viewed as an important method to conserve biodiversity in agricultural landscapes (e.g. 18), our results highlight that increasing the quality of these habitatsthrough the conservation of old permanent grasslands in the landscapemay be equally important. Finally, we observed that mobility trait diversity decreased linearly with time since landscape conversion to cropland, while resource-acquisition trait diversity showed a hump-shaped response: it increased during the first 10 years after grassland-to-crop conversion but then decreased (Fig. 2). Spatial and temporal disturbances can create a concentration of transient species in remnant habitats, thus locally increasing resource-acquisition trait diversity [START_REF] Tscharntke | Landscape moderation of biodiversity patterns and processes-eight hypotheses[END_REF]. However, resourceacquisition trait diversity will decrease after some time due to the scarcity of habitats and resources in the landscape [START_REF] Kuussaari | Extinction debt: a challenge for biodiversity conservation[END_REF]. Quantifying the duration of this effect may be essential to slow down species loss and delay extinction.

Resource-acquisition trait diversity was also strongly impacted by current land use. Locally, we found a strong negative effect of the grassland field productivity (Fig. 3). High grassland productivity is associated with high rates of fertilization, which decreases plant diversity, and in turn, the availability of feeding niches [START_REF] Harpole | Addition of multiple limiting resources reduces grassland diversity[END_REF]. This has been observed in both in experiments and realworld landscapes, and leads to a global homogenization of species resource-use strategies at higher trophic levels [START_REF] Gossner | Land-use intensification causes multitrophic homogenization of grassland communities[END_REF]. In addition, at the landscape scale, the proportion of mass flowering crops (ranging from 0% to 44% in the study area) negatively affected resource-acquisition trait diversity.

Mass flowering crops offer short pulses of uniform and homogeneous resources for herbivores and mutualists that may select for particular resource-acquisition traits and therefore reduce their diversity (e.g. [START_REF] Kormann | Local and landscape management drive trait-mediated biodiversity of nine taxa on small grassland fragments[END_REF]. As the diversity of resource-acquisition traits drives overall resource utilization [START_REF] Naeem | Declining biodiversity can alter the performance of ecosystems[END_REF], the community level homogenization of resource-acquisition traits following land-use intensification could threaten the delivery of essential agroecosystem services such as nutrient recycling [START_REF] Deraison | Functional trait diversity across trophic levels determines herbivore impact on plant community biomass[END_REF], pollination [START_REF] Holzschuh | Mass-flowering crops dilute pollinator abundance in agricultural landscapes across Europe[END_REF] and biological control [START_REF] Greenop | Functional diversity positively affects prey suppression by invertebrate predators: a meta-analysis[END_REF].

We did not find clear response of body size to land-use changes (Fig. 1 and SI Appendix, Fig. S6), in contrast to other studies focusing on single trophic groups (e.g. 37). Body size is an integrative trait related to many aspects of animal species physiology and ecology (e.g. metabolism, ref. [START_REF] Bebber | The global spread of crop pests and pathogens[END_REF]; mobility and dispersal, or else stoichiometry, ref. [START_REF] Mayfield | What does species richness tell us about functional trait diversity? Predictions and evidence for responses of species and functional trait diversity to land-use change[END_REF]) and may not have the same functional significance across multiple trophic groups [START_REF] Häussler | The biggest losers: Habitat isolation deconstructs complex food webs from top to bottom[END_REF]. As a result, our synthetic index of multitrophic body size may not show a consistent response to land-use changes. This may explain why our analyses predicted less variation in body size trait diversity, than they did for mobility or resourceacquisition traits.

We found that past and current land use had contrasting impacts on trait diversity of different trophic groups (Fig. 3). For instance, primary consumers (i.e. herbivores and mutualists) responded more strongly to local history (through the age of the sampled grasslands) while predators were more impacted by past and present-day landscape composition (Fig. 4 and SI Appendix, Table S6). Accordingly, delayed responses to land-use changes could create temporal mismatches between interacting trophic levels, if land-use changes disrupt interacting partners inconsistently. This may lead to the disruption of ecological interactions and further biodiversity loss [START_REF] Valiente-Banuet | Beyond species loss: the extinction of ecological interactions in a changing world[END_REF]. While plants are expected to experience time-delayed responses to habitat loss because of persistence in seed banks [START_REF] Purschke | Interactive effects of landscape history and current management on dispersal trait diversity in grassland plant communities[END_REF], plant trait diversity was not impacted by land-use history (Fig. 3). However, we found a strong effect of land-use history on plant species richness (SI Appendix, Fig. S7), suggesting that legacy effects of land use affect plant species independently of the traits investigated in this study. While the functional diversity of certain trophic groups showed divergent responses to land-use changes (Figs. 4), our study shows that there is a consistent pattern in the response of overall ecosystem trait diversity, thus suggesting that our integrated index of multi-trait diversity provides a simple quantitative measure of the whole ecosystem biodiversity response.

Our results revealed that present-day functional biodiversity was impacted by land-use changes that happened up to 20 years ago, thus demonstrating that land-use actions have longlasting impacts and that extinction debts (sensu 9) are commonplace in agricultural landscapes [START_REF] Krauss | Habitat fragmentation causes immediate and time-delayed biodiversity loss at different trophic levels[END_REF].

In an era of global biodiversity change, our results emphasize the need to consider land-use legacies in conservation programs that aim to protect the biodiversity of agroecosystems and its associated ecosystem functions and services. In particular, our results emphasize the need to preserve permanent grasslands in agricultural landscapes, as it provides a shelter to low mobility organisms and species with narrow feeding niches. While large-scale policy schemes encourage the retention of permanent grasslands, they have largely decreased in their cover over the last years [START_REF] Pe'er | EU agricultural reform fails on biodiversity[END_REF]. Our results call for immediate actions to conserve and restore permanent grasslands in order to preserve the functionality of agricultural landscapes and avoid future extinction debt.

MATERIALS AND METHODS

Study area. The study was conducted in 2014 in the Long Term Socio-Ecological Research (LTSER) "Zone atelier Plaine & Val de Sèvre" (ZAPVS) located in western France [START_REF] Bretagnolle | Description of long-term monitoring of farmland biodiversity in a LTSER[END_REF]. The LTSER covered approximately 430 km2 of an intensively managed agricultural plain. Historically, it was a typical rural area characterized by the presence of mixed crops-livestock systems (dairy goats and cows). Grassland was the dominant land use fifty years ago covering about 60% of the area [START_REF] Bretagnolle | Description of long-term monitoring of farmland biodiversity in a LTSER[END_REF]. Since that time, shifts from grazing livestock to feeding livestock, and shifts from livestock production to annual crop production have resulted in a strong decline in grassland cover.

In 2014, grasslands covered about 12% of the area and included artificial grasslands (i.e., alfalfa with 3% of the area), temporary (sown with pure grasses or in mixtures with legume species and ≤ 5 year old) and permanent grasslands (> 5 year old) managed by grazing, mowing or abandoned. The remaining areas were covered by crops (66% of landscape area). Soils are mostly composed of karst, with calcareous rocks providing shallow calcareous soils with low water retention and pH > 7. Since 1994, land cover of the study area has been monitored on a yearly basis at the field scale, by using about 30 land-use types (11,000 fields approximately, see ref. [START_REF] Diekötter | Oilseed rape crops distort plantpollinator interactions[END_REF] for methodological details on land cover monitoring), and has been stored in a Geographical Information System (GIS) database, running on QGIS v 2.14.

Grassland selection and land-use metrics.

We monitored 75 grassland fields within the study area. The grassland fields were selected among hay meadows of varying ages and vegetation types (pure legumes, pure or mixed grasses, legume and grass mixture, post-cultural vegetation). The average age of the grassland field was 8 year-old (sd = 6.55). We calculated for each of the 75 grassland fields the current landscape composition, and landscape metrics linked to land-use history, within a 1 km-radius landscape (SI Appendix, Fig. S1). This scale was chosen to approximate the dispersal distance of different taxa (e.g. 39). For landscape composition, we considered landscape elements known as favorable or resource rich habitats for the different groups of taxa, i.e. the proportion of the landscape covered by grasslands, forests and mass flowering crops (oilseed rape and sunflower). Landscape composition metrics varied between 0-35 % for grassland, 0-32 % for forests and 0-44 % for mass flowering crop covers. For land-use history metrics, we calculated for all field in the 1 km-radius landscape surrounding the focal grasslands: (i) the time elapsed since the first grassland-to-crop conversion (SI Appendix, Fig. S1); (ii) the time spent into current land use since the last tillage (hereafter called "Field age"). The age of crop-field was 1 year and the age of grassland field varied between 1 and 20 years (set to 20 years if it had not been ploughed since 1994); and (iii) the rate of turnover from crop to grassland since the last conversion of the field into cropland. This metric is linked to the lifespan of grasslands sown in crop rotation and is higher when the average lifespan of sown grasslands in the landscape is longer. To account for the size of the field, all these metrics were weighted by the field area. We then averaged these three metrics at the landscape level (hereafter called "Time as cropland", "Grassland permanency" and "Crop-grassland turnover rate") (SI Appendix, Fig. S1).

Species richness and density. Seven taxonomic groups were sampled on the 75 grassland fields.

The taxonomic groups were chosen for their association with key ecological functions: primary producers (vascular plants), pollinators (wild bees and hoverflies), herbivores involved in carbon and nitrogen cycling (grasshoppers), invertebrate predators that are natural enemies of pests (carabid beetles and spiders) and vertebrate predators (birds) of important cultural value (see SI Appendix, Table S2). Sampling was performed throughout the growing season from April to August 2014 in each grassland field following standardized protocols [START_REF] Bretagnolle | Towards sustainable and multifunctional agriculture in farmland landscapes: Lessons from the integrative approach of a French LTSER platform[END_REF]. Details of sampling methods are provided in SI Appendix Methods S1. Relative abundance per plant species was calculated as the sum of the species cover in the 10 quadrats divided by the total cover of all species. For animals, species density was estimated as the number of individuals captured (invertebrates) or recorded (birds) divided by sampling intensity (number of traps or point counts in the grassland). The number of species in each grassland field was determined in the lab for all invertebrate species.

Trait measurements and trait diversity. Functional trait data were collected for all taxonomic groups. Plant trait data came from a local database from the LSTER site [START_REF] Deraison | Functional trait diversity across trophic levels determines herbivore impact on plant community biomass[END_REF]. Bird morphological traits were compiled from the literature [START_REF] Cramp | Handbook of the birds of Europe, the Middle East and North Africa[END_REF]. For all other taxonomic groups, all measurements were performed using a stereo microscope (Leica Microsystems M50) equipped with an integrated high definition microscope camera (Leica IC80 HD). Details of trait measurements are provided in SI Appendix, Methods S2. In total, we collected trait data for 178 species distributed across the different trophic groups (SI Appendix, Methods S2).

We calculated the community abundance-weighted variance [START_REF] Enquist | Chapter Nine-Scaling from Traits to Ecosystems: Developing a General Trait Driver Theory via Integrating Trait-Based and Metabolic Scaling Theories[END_REF] for each taxonomic group separately and each trait separately. For single traits, the variance is a measure of the trait dispersion within a given community weighted by the abundance of each individual species and is a measure of functional trait diversity of a given community [START_REF] Enquist | Chapter Nine-Scaling from Traits to Ecosystems: Developing a General Trait Driver Theory via Integrating Trait-Based and Metabolic Scaling Theories[END_REF]. We also computed a multi-trait index of functional diversity based on trait dispersions for each taxonomic group separatelythe Functional Dispersion (FD) [START_REF] Laliberté | A distance-based framework for measuring functional diversity from multiple traits[END_REF].

Multitrophic trait diversity measurements. We used methods developed to study ecosystem multifunctionality to calculate a multi-trait diversity index considering on all traits and taxonomic groups. This approach (averaging approach, see 54) consists of calculating the average standardized values of multiple functional diversities. Our multi-trait functional diversity index was thus calculated as the average standardized FD values across the taxonomic groups (similar to multidiversity indices following 8). Using this index allows us to identify the environmental conditions, which maximize functional trait diversity across multiple taxonomic groups. Note that the averaged multi-trait diversity index was highly correlated with the threshold-50 multi-trait diversity index (37) -calculated as the percentage of community variance values that exceeds 50% of their maximum observed community variance (r=0.87, see SI Appendix, Table S3). Similarly, we calculated an average index of multitrophic densityaverage standardized values of the total density for the different animal taxonomic groups (i.e. grasshoppers, wild bees, hoverflies, carabids, spiders and birds)and an average index of multitrophic species richness, following 8. In addition, to test if land-use changes have similar effects on different axes of trait variation, we also calculated averaged indices for mobility trait diversity, resource-acquisition trait diversity and body size trait diversity values as the average standardized community variance values of these three types of traits, across the taxonomic groups. To calculate these trait indices, we excluded plant traits as the three selected plant traits are associated with different functions and other axes of variation. Finally, to test if the response of each trophic group to land-use changes is similar, we computed multi-trait diversity, species richness and density indices for five trophic groups separately: primary producers (vascular plants), mutualists (wild bees and hoverflies), herbivores (grasshoppers), invertebrate predators (carabids and spiders) and vertebrate predators.

Assessing the effects of past and current land use. We evaluated the effects of land-use changes on grassland trait diversity by using linear models. We ran separate analyses on: (i) the multi-trait diversity; (ii) the mobility, resource-acquisition and body size trait diversity; (iii) on the multi-trait diversity of the different trophic groups separately. We included in our models the effects of the grassland field age and of the three landscape land-use history metrics as predictors of functional trait diversity. We also controlled for current landscape composition (% grassland, % forest and % mass flowering crop areas). We used plant productivity as a proxy measure of grassland management, as it is related to fertilizer inputs. This was assessed by harvesting plant biomass each month (between February and August 2014) above a cutting height of 5 cm from the soil surface, within five 35 × 35 cm quadrats. Productivity was then calculated as the weight of dried-plant material (oven-dried at 60°C for 72h) product per square meter per day between the initial biomass measurement and peak biomass (end of May). We considered quadratic terms for grassland field age, landscape land-use history metrics, current grassland field productivity and landscape composition to assess potential non-linear effects of these variables. Further co-varying factors accounted for the number of mowing events and soil depth were included. Our models also integrated the coordinates of the centroid of each sampled grassland (latitude and longitude) to correct for additional spatial effects not accounted for by the field and landscape predictors [START_REF] Bagousse-Pinguet | Phylogenetic, functional, and taxonomic richness have both positive and negative effects on ecosystem multifunctionality[END_REF].

Finally, as the local species pool (defined by both the density and the species richness) may affect functional trait diversity metrics [START_REF] Mayfield | What does species richness tell us about functional trait diversity? Predictions and evidence for responses of species and functional trait diversity to land-use change[END_REF], we included in the models the multitrophic density and multitrophic species richness indices to isolate the net effects of land-use changes on multitrophic trait diversity. We ran separate analyses on the multitrophic density and the species richness, and results are presented in SI Appendix, Fig. S4 and Fig. S5.

We then performed a series of models for the multi-trait diversity, mobility trait diversity, resource-acquisition trait diversity and body size trait diversity; and for the multi-trait diversity of each taxonomic group separately. We first used a backward stepwise regression procedure using the software JMP 11 (SAS Institute) to select, between all models, the best-fitting models with lower AICc (Δ AICc < 2). Second, using the best selected models, we performed a model-averaging procedure based on AICc selection (delta AICc < 2) to determine parameter coefficients for the best final set of predictors of our response variables. This procedure was performed using the function dredge in the R package Multi-Model Inference (MuMIn) [START_REF] Bartoń | MuMIn: Model selection and model averaging based on information criteria (AICc and alike)[END_REF]. Model residuals were inspected for constant variance and normality. We standardized all variables (z-scored: mean-centred and divided by the standard deviation) to interpret parameter estimates on a comparable scale [START_REF] Schielzeth | Simple means to improve the interpretability of regression coefficients[END_REF]. Correlation among the predictors used was low (SI Appendix, Table S1) and did not induce multicollinearity issues in our analyses (SI Appendix, Table S4). The inclusion of many predictors in statistical models increases the chance of type I error (false positives). To account for this we used a Benjamini and Hochberg procedure to control for false discovery rates and adjust P-values [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF]. To evaluate the relative importance of the predictors as drivers of trait diversity, we expressed the importance of predictors as the percentage of variance they explain, based on the comparison between the absolute values of their standardized regression coefficients and the sum of all standardized regression coefficients from the predictors. This method is similar to a variance partitioning analysis because we previously transformed all predictors to z-scores [START_REF] Bagousse-Pinguet | Phylogenetic, functional, and taxonomic richness have both positive and negative effects on ecosystem multifunctionality[END_REF][START_REF] Gross | Functional trait diversity maximizes ecosystem multifunctionality[END_REF][START_REF] Sirami | Increasing crop heterogeneity enhances multitrophic diversity across agricultural regions[END_REF]. The following identifiable variance fractions were then examined: (i) land-use history, (ii) current land use, and (iii) the species pool.

FIGURE LEGENDS

Fig. 1. Importance of the drivers of multi-trait diversity, mobility trait diversity, resourceacquisition trait diversity and body size trait diversity. Relative effects (% R²), resulting from a model averaging procedure, were calculated for each group of predictors (i.e. land-use history, current land use and the species pool). All predictors were scaled to interpret parameter estimates on a comparable scale. Note that for mobility, resource-acquisition and body size trait diversity, we focused on animal traits and excluded plant traits from the analyses. Results were consistent considering spatial scales ranging from 500-1500 m radii surrounding the sampled grasslands (SI Appendix, Fig. S3). Fig. 2. Effects of land-use history (field age, time as cropland, grassland permanency and cropgrassland turnover rate), current land use (field productivity, % grassland, % mass flowering, % forest) and of the species pool (multitrophic density and richness) on (A) multi-trait diversity, (B) mobility trait diversity and (C) resource-acquisition trait diversity. Lines show model fits and shaded areas correspond to the 95% confidence band. Model predictions were calculated using a model averaging procedure. Note that for mobility trait and resource-acquisition trait diversity, we focused on animal traits and excluded plant traits from the analyses. All predictors were scaled to interpret parameter estimates on a comparable scale. P-values of the best selected models for each model parameter are given. P-values were adjusted following (57) to control for false discovery rates (see SI Appendix, Table S7). Gray dots (N=75) correspond to observed data. Fig. 3. Importance of the drivers of multi-trait diversity within the different trophic groups (primary producers, herbivores, mutualists, invertebrate predators and vertebrate predators). Relative effects (% R²), resulting from the model averaging procedure, were calculated for each group of predictors (i.e. land-use history, current land use and the species pool). All predictors were scaled to interpret parameter estimates on a comparable scale. Fig. 4. Effects of land-use history (field age, time as cropland, grassland permanency and cropgrassland turnover rate) and current land use (field productivity, % grassland, % mass flowering, % forest) on the multi-trait diversity of each trophic group. Lines show model fits for primary producers (green), herbivores (yellow), mutualists (orange), invertebrate predators (brown) and vertebrate predators (blue). Model predictions were calculated using the model averaging procedure. All predictors were scaled to interpret parameter estimates on a comparable scale. Pvalues of the best selected models for each model parameter are given, *, P-value < 0.05; ** Pvalue < 0.01; *** P-value < 0.001. P-values were adjusted following (58) to control for false discovery rates (see SI Appendix, Table S7). Dots (N=75) correspond to observed data and colors indicate the trophic group.

We thank Franck Coudray, Hélène Deraison, Yoann Erreca, Thierry Fanjas-Mercère, Jean-Luc Gautier, Louis Gross, Nadine Guillon, Florian Mezerette, Sophie Pillaud, Alexis Saintilan, Edoardo Tedesco for field assistance. The study was supported by the French government Initiatives d'Excellence -Initiatives Science/Innovation/Territoires/Économie (IDEX-ISITE) initiative 16-IDEX-0001 (CAP 20-25), the FarmLand research programme, an ERA-Net BiodivERsA project funded by the French Agence Nationale de la Recherche (ANR-11-EBID-0004), and the ECODEAL research programme, funded through the 2013-2014 BiodivERsA/FACCE-JPI joint call for research proposals. G.L.P was supported by a region Poitou-Charentes -department Deux-Sèvres PhD grant. N.G. was supported by the AgreenSkills+ fellowship programme, which has received funding from the European Union's Seventh Framework Programme under Grant Agreement FP7-609398 (AgreenSkills+ contract). Y.L.B.P. was supported by the European Research Council (BIODESERT) and by a Marie Sklodowska-Curie Actions Individual Fellowship (MSCA-IF) within the European Program Horizon 2020 (DRYFUN Project 656035). We thank Eric Allan, Sandra Lavorel, Hervé Jactel, Colin Fontaine and Pascal Carrère for their useful comments on the previous versions of the manuscript.

DATA AVAILABILITY

The data that support the findings of this study will be openly available at https://datadryad.org.

ACKOWLEDGMENTS