Metabolomics approach reveals disruption of metabolic pathways in the marine bivalve Mytilus galloprovincialis exposed to WWTP effluent

Introduction

The development of high resolution mass spectrometry (HRMS) has enabled to screen, in an untargeted way, hundreds of compounds at trace levels in any environmental matrices (e.g. water, sediment, organism). The power and sensitivity of this technology make it relevant and suitable to answer ecotoxicological issues, and especially through environmental metabolomics [START_REF] Lankadurai | Environmental metabolomics: an emerging approach to study organism responses to environmental stressors[END_REF]. Indeed, HRMS-based chemical profiles collected on fluids or biological tissues of exposed organisms are likely to reflect exposure by the direct non-targeted detection of parent compounds and their metabolites (exposome). They may also reflect effects by detection and characterization of the modulation of endogenous metabolites (endometabolome) that may lead to physiological disturbances induced by the chemical exposure [START_REF] Courant | Basics of mass spectrometry based metabolomics[END_REF]. The aim of the present study was to perform an untargeted metabolomics approach to highlight endogenous metabolites disrupted by an urban wastewater treatment plant (WWTP) effluent exposure and then elucidate the biological effects of such exposure on Mediterranean mussels (Mytilus galloprovincialis). An untargeted analysis was also conducted on the WWTP effluent extract in order to characterise the exposure. Since pharmaceutical and personal care products consumption has been rising during the last decades, these compounds may constitute a large class of contaminants present in the urban WWTP effluent (470,000 inhabitants equivalent). Hence, the following WWTP effluent preparation was focused on the extraction of dissolved hydrophilic compounds.

Materials and methods

The WWTP effluent was prepared as follows. 15L of effluent was deposited on 30 cartridges (Oasis HLB® 500 mg, 6 cc) conditioned with 5 mL MeOH and 5 mL distilled water. Cartridges were then washed (5 mL Milli Q water), dried and eluted each with 2 x 5 mL of MeOH. The eluate was evaporated until 3 mL of extract was obtained. The exposure was performed in controlled laboratory conditions. After seven days of acclimatisation, mussels from a treated group (n=25; 5 tanks containing 5 mussels each per 2L seawater) were exposed to 100 µL of the WWTP effluent extract corresponding to an environmental dilution of the effluent of 5%. Mussels from the control group (n=25) were exposed to solvent (100µL MeOH into each tank of 2L of seawater). During the acclimatisation and exposure periods, filtered seawater was continuously aerated and renewed daily (static renewal). Mussels were fed once daily with the marine green alga Tetraselmis suecica at constant density (10,000 cells/mL). Seven days following the exposure, mussels were dissected and their sex determined. No mortality was observed during the exposure period. Metabolic fingerprints were generated from digestive gland (control group: n=9 males and n=10 females; exposed group: n=10 males and n=10 females) by LC-HRMS analysis (reverse phase PFPP column, Exactive Orbitrap) operating in positive and negative ionisation mode. Afterwards, univariate and multivariate tests were applied to highlight disrupted metabolites in male and female following exposure. A gender specific response was then assessed quantitatively and qualitatively by an ANOVA Multiblock OPLS (AMOPLS). This method takes advantage from the structured design of multifactorial experiments to decompose and quantify the sources of variability, and investigate main effects (exposure, gender) or interaction(s) (exposure×gender) using dedicated components [START_REF] Boccard | Exploring Omics data from designed experiments using analysis of variance multiblock Orthogonal Partial Least Squares[END_REF].

Results and discussion

3.1 Characterisation of the WWTP effluent extract Untargeted analysis of the WWTP effluent extract revealed numerous signals in both negative (ESI-; n=8194) and positive (ESI+; n=5773) electrospray ionisation modes. Among those signals, the Toxin and Toxin Targeted Database [START_REF] Wishart | Direct interactions between metabolites are indicated with full arrows[END_REF] allowed the chemical annotation of 2912 and 1622 signals in ESI-and ESI+ respectively, based on an accurate mass tolerance of 0.002 Da. Majority of annotated compounds belonged to the classes "Drug", "Pesticide" or "Household compound", confirming the main urban nature of the WWTP effluent. Raw formulae of chlorine or sulphur containing contaminants were confirmed on the basis of chlorine or sulphur isotope ratio in mass spectrometry (e.g. sulfamethoxazole, nordazepam, fenofibric acid, fipronil, etc). When available in the laboratory, annotations were confirmed by the injection of the corresponding analytical standards in the same analytical conditions (diclofenac, venlafaxine, ofloxacin, oxazepam, 4-nonylphenol, etc.).

Disruption of male mussel endometabolome

The study of the endometabolome showed numerous metabolites modulated on treated male mussels in response to WWTP effluent exposure. The highlighted metabolites belong mainly to amino acids metabolism (e.g. tyrosine, phenylalanine, proline, etc.), neurohormones (dopamine and a serotonin metabolite), purine and pyrimidine metabolism (e.g. adenosine, guanine, uracil etc.), citric acid cycle intermediates (e.g. malate, fumarate), and a component involved in oxidative stress defence (oxidised glutathione). Figure 2 shows an example of altered metabolic pathways, i.e. phenylalanine and tyrosine metabolism. Modulation of these metabolites could have effects on several biological processes, which may lead to a negative impact of organism fitness, such as energy metabolism, DNA, RNA and protein synthesis, osmoregulation, byssus formation and adhesion, and reproduction.

Do males and females respond identically?

Female mussels also showed many impacted metabolites. Some of them were common to male, while many others seemed to be specific of female response. Thanks to the AMOPLS analysis, we were able to decompose the variability: the exposure effect (independently of gender) represented 7% (p<0.01) of the total variability, and the exposure×gender interaction, i.e. the gender specific response to WWTP effluent exposure, represented 4% (p<0.05) of the total variability. Meaningful biochemical information was obtained from the component related to each effect. The former highlighted metabolites modulated in response to exposure whatever the gender. Those relevant metabolites may constitute a set of biomarkers after further validation processes. The latter highlighted key metabolites involved in gender specific response (for example dopamine modulation was only observed in male in response to the WWTP exposure). This may provide a better understanding of response variabilities to the exposure.

Conclusion

Our work demonstrates the relevance of an untargeted metabolomics approach to characterise both exposome and endometabolome. We highlighted key molecular events triggered by chemical exposure, which could lead to adverse outcomes on individual fitness. In addition, it underlines that gender has to be taken into account for further research on response of organisms to stress. AMOPLS is a promising approach to better characterise multiple factors of variability in response to a stress. Keywords: metabolomics, LC-HRMS, wastewater treatment plant effluent, metabolic pathways, Mytilus galloprovincialis


F Calvo

P Poulain

Transitions between secondary structures in isolated polyalanines

Keywords: 64, 60, Cn Order-disorder transformations; statistical mechanics of model systems -02, 70, Uu Applications of Monte Carlo methods -87, 14, Ee Proteins

Monte Carlo simulations of gas-phase polyalanine peptides have been carried out with the Amber ff96 force field. A low-temperature structural transition takes place between the α-helix stable conformation and β-sheet structures, followed by the unfolding phase change. The transition state ensembles connecting the helix and sheet conformations are investigated by sampling the energy landscape along specific geometric order parameters as putative reaction coordinates, namely the electric dipole µ, the end-to-end distance d, and the gyration radius Rg. By performing series of shooting trajectories, the committor probabilities and their distributions are obtained, revealing that only the electric dipole provides a satisfactory transition coordinate for the α ↔ β interconversion. The nucleus at the transition is found to have a high helical content.

Introduction

The structure of proteins is closely related to their biological function. In many cases, the possibility to undergo conformational transitions can be another important property for the function to be fully acquired, e.g. for transport [1] or motor [2,3] proteins, or in enzymes [4]. Misfolding of some proteins is due to an unwanted conformational change into a wrong, non-native structure often leading to fatal diseases. Prion proteins and the native Aβ peptide in Alzheimer's disease are two important examples of αhelix-rich proteins which change to β-sheet upon aggregation into the so-called amyloid fibrils [5]. The competition between α helices and β sheets is also observed in proteins with nonhierarchical folding, as in β-lactoglobulin where α-helix intermediates are seen as transient conformations during folding into the β-sheet native structure [6].

While experimental methods such as high-resolution X-ray or NMR can often characterize the end states of a conformational change in solution, the entire pathway is much harder to study. Fluorescence resonance energy transfer provides detailed information about specific intermolecular distances and their time evolution, however these partial data are far from full atomistic accuracy. In the gas phase, several methods have been developed during the last decade to get insight into the structural properties of proteins as an alternative to the solvent or matrix characterizations. Studying isolated biomolecules also has the advantage of putting aside the intermolecular interactions, thus simplifying greatly the analysis [START_REF] De Vries | Special issue on Molecular Physics of Building Blocks of Life under Isolated or Defined Conditions[END_REF]. Infrared spectroscopy [8][9][10], ion mobility [11,12], H/D exchange reactions [START_REF] Oh | Proc. Natl. Acad. Sci. USA[END_REF], peptide fragmentation [START_REF] Oh | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Horn | [END_REF][15] or electric dipole measurements [16] have all contributed significantly to our understanding of gas-phase peptides. In particular, the existence of the α-helix and β-sheet conformations has been attested in a number of cases [17]. Yet, determining transition pathways remains an issue in vacuo, and the results of computer simulations can be invaluable.

Conformational changes in silico can be addressed in a number of ways, depending on the precise aims. Direct molecular dynamics [18][19][20][START_REF] Ohkubo | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Yang | [END_REF][23] and Monte Carlo (MC) simulations [24][25][26][27] provide useful insight into the energetics and kinetics of the folding transition, especially in the case of polyalanines. However, in their basic form they are usually not appropriate for studying rare events involving large free energy barriers. From the structural point of view, one would like to characterize not only the two end states, but also the optimal set of reaction coordinates along which the transition is best visualized. The common approach here is intuition, by selecting a priori order parameters defined according to the molecular geometry, and by calculating the probability distribution of finding these variables and the corresponding free energy maps. Upon identification of barriers between the two states, additional trajectories can be run to calculate the rate constants via transition path sampling (TPS) [28,29] or related schemes [30][31][32]. More recently, several groups [33][34][START_REF] Best | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Peters | [END_REF] have proposed improvements of the pioneering TPS method in which the quality of the reaction coordinate is assessed by computing the committor probability p c , that is the chance that the system falls into one state earlier than in the other state. This quantity, though it may not appear immediately physical, allows a proper definition of a transition state as the set of conformations for which p c = 1/2.

In a previous work, we have theoretically shown that gas-phase polyalanines can exhibit a thermally stable βsheet intermediate between the low-temperature α-helical native state and the high-temperature unfolded state [37]. These results are consistent with experimental measurements [38] and also with previous theoretical suggestions from other groups [39][START_REF] Levy | Proc. Natl. Acad. Sci[END_REF][START_REF] Ding | [END_REF]. Here we focus on the α-helix ↔ β-sheet interconversion, which we attempt to characterize in terms of specific order parameters. More precisely, we have chosen the electric dipole µ, the end-to-end distance d between the N-terminal atom and the H atom from the hydroxyl group in the carbonyl terminal, and finally the gyration radius R g . Near the α-β equilibrium transition temperature, the free energies along these three parameters show stable basins corresponding to both conformations. However, as will be seen below, dividing surfaces are not straightforward to identify and, more importantly, they may not reflect the true transition states. By shooting many trajectories from several regions of the energy landscape, the committor probabilities can be calculated, thus helping in assessing the relevances of these order parameters as reaction coordinates for the α ↔ β interconversion.

The article is organized as follows. In the next section, we describe all the computational methods used to sample the potential energy surface of our model system, the isolated octa-alanine. The equilibrium results and the transition state ensembles associated with each of the three order parameters are critically discussed in Sec. 3, and some concluding remarks are finally given in Sec. 4.

Methods

We have chosen a short enough polyalanine that is found stable in both the α-helical and β-hairpin configurations, namely the octa-alanine system. This peptide is modelled using the Amber ff96 force field [42]. Because the original parameters of Amber ff96 were originally fitted to reproduce the properties of hydrated biomolecules, the partial charges had to be increased to account for polarization effects. When used in the gas phase, this parameter set overstabilizes β-sheet structures with respect to helices [44]. To circumvent this problem and shield the charges, a dielectric constant ε r = 2 was used in the simulations, thus satisfactorily reproducing the electric dipole measured in deflection experiments [43].

We have used two kinds of methods in this work. A first goal is to sample the energy landscape of our polypeptide at thermal equilibrium, in order to characterize the α ↔ β interconversion and to calculate Landau free-energy curves corresponding to specific order parameters. In a second part, the quality of these order parameters as suitable reaction coordinates will be assessed by shooting trajectories from putative transition states and estimating the committor averages and their probability distributions.

Sampling the energy landscape

The potential energy surface of the octa-alanine was explored using parallel tempering (PT) Monte Carlo simulations as a reference. 32 replicas were propagated simultaneously in the temperature range 50 K ≤ T ≤ 1000 K with a geometric progression. Ten series of 10 6 Monte Carlo cycles per replica were propagated consecutively, the last configurations of each simulation providing the starting conditions for the next one. Here one MC move consisted of randomly selecting one torsion angle of the backbone or side chains, and randomly rotating it by an angle δθ drawn from a range [-θ max , θ max ], with θ max adjusted at each temperature to get approximately 50% acceptance rate according to the Metropolis rule. Here one MC cycle is a series of 23 consecutive MC moves. After each cycle, one exchange between two random adjacent replicas was attempted with 10% probability.

While the average potential energy E or its fluctuations ∆E 2 suffice to monitor the unfolding transition, other quantities are required to locate the presence of α-helix or β-sheet secondary structures. We thus calculated the electric dipole vector µ and its modulus µ, the end-to-end distance d between the nitrogen atom in Nter position and the hydrogen atom from the hydroxyl group in C-ter position, as well as the gyration radius R g . In addition to the standard thermal averages µ , d and R g , we recorded the two-dimensional histograms h T (E, A) as a function of potential energy and geometric parameter A = µ, d or R g at the temperature T . These histograms were then processed into the joint densities of states Ω(E, A), and subsequently into the thermal averages E , A or their fluctuations as a continuous function of temperature, via a multiple histogram analysis.

Complementary Monte Carlo simulations were performed using the Wang-Landau (WL) procedure [45] for joint densities of states, recently improved using our annealing scheme [46]. Here the two-dimensional density of states Ω(E, A) is numerically calculated with increasing accuracy by penalizing each visited state more and more as it is being visited during the exploration. This algorithm eventually converges to the true density of states, which in turns provides all thermodynamical properties in the canonical ensemble after Laplace transformation.

The two-dimensional Wang-Landau method, previously used with the electric dipole as the extra variable A, was found to perform very well for small peptides [46] when compared with parallel tempering Monte Carlo. Here we have chosen to repeat these simulations in the same conditions, but with the end-to-end distance d and the gyration radius R g as the two extra variables to the potential energy E. The simulation details for these Wang-Landau simulations are the same as in Ref. [46], and will not be further detailed here. In particular, the inverval in A was discretized into 200 bins for the three variables and the energy interval consisted of 300 bins. The ranges for the three geometric variables were 0 ≤ µ ≤ 40 D, 0 ≤ d ≤ 35 Å, and 4 Å≤ R g ≤ 10 Å, respectively.

The microcanonical densities of states Ω(E, A) can also be used to calculate the canonical probability of find-ing the system with parameter A within some error dA at temperature T , again via a simple Laplace transformation:

p(A, T )dA ∝ dA Ω(E, A) exp(-E/k B T )dE, (1)
with k B the Boltzmann constant. From the probability p(A, T ), the Landau free energy F (A, T ) is defined at temperature T as

F (A, T ) = -k B T ln p(A, T ). (2)
Bistable systems are expected to show two minima in the Landau free energy as a function of a suitable function A. At the transition temperature, and for a first-order like phase transition, the two minima should have the same depth.

Transition state ensembles

The shape of the probability p(A, T ) or the free energy F (A, T ) provides a great deal of information about the stable regions of the energy landscape, after projection onto the physically meaningful order parameters A. However, even the presence of two comparable minima in F (A, T * ) near the equilibrium temperature T * does not qualify the parameter A as a good reaction coordinate.

In order to determine whether the parameter A is an adequate transition coordinate, one must calculate the committor probability p c (A) at the barrier A † , that is the probability that a trajectory initiated at A † will fall into one minimum earlier than in the other minimum. By definition, the transition region between the two minima is the isocommittor surface, or the set of points in phase space for which p c = 1/2. Therefore a necessary condition for A to be a good transition coordinate is that the average committor probability p c taken for points having A ≃ A † is close to one half. However, as discussed in particular by Best and Hummer [START_REF] Best | Proc. Natl. Acad. Sci. USA[END_REF] and by Geissler and coworkers [47], this necessary condition on the average p c over a set of configurations is not sufficient. It may recover situations where this set actually comprises of state points already belonging to the two end states, resulting in an apparent free energy barrier only after projecting onto a wrong coordinate and further statistical averaging. A more rigorous approach is to calculate the distribution P(p c ) of the probability p c after repeating the calculation of p c from many different starting points at A † . If P(p c ) has a peaked shape centered at 1/2, then the variable A can be considered as a suitable reaction coordinate for the transition [START_REF] Best | Proc. Natl. Acad. Sci. USA[END_REF]47]. Conversely, a bimodal distribution P(p c ) peaked at the extremities p c = 0 and p c = 1 cannot be considered as appropriate. This situation was illustrated by Du and coworkers [48] who found that neither the number of native contacts nor the looplength distributions could describe the folding transition of lattice models of proteins.

For the present system, we have first sampled sets of transition points by saving configurations corresponding to the desired value of A 0 within some small range ±δA/2, selecting A 0 near the maximum of the free energy barrier A † separating the α-helix and β-sheet minima. This sampling was carried out using the parallel tempering Monte Carlo simulations without any bias for the parameters d and R g . Sampling configurations according to their electric dipole turned out to be more difficult, due to the much higher energy barriers reflecting very low occurence probabilities. For this parameter, and following the suggestion of Best and Hummer [START_REF] Best | Proc. Natl. Acad. Sci. USA[END_REF], an harmonic umbrella potential W (µ) = k(µ -µ 0 ) 2 around the desired value µ 0 was added to the potential energy, with the spring constant k = 5 kcal/mol/D 2 .

Once a set {R}(A 0 ) of N configurations characterized by their value of A(R) ∈ [A 0 -δA/2, A 0 + δA/2] had been selected, the committor probabilities p c (R) were estimated by shooting independent Monte Carlo trajectories initiated from each of these configurations. These 'aimless' shooting trajectories started at the temperature T = 150 K, which is below the folding temperature but above the α ↔ β interconversion, and after 5 × 10 4 MC cycles the temperature was decreased to 80 K for another quenching round of 5 × 10 4 MC cycles. For each starting configuration R, the number M of trajectories converging to either of the α-helical or β-sheet states after the 10 5 MC cycles was required to be 50. If one simulation did not end into one of these states before the 10 5 MC cycles, then it was just ignored and a new independent simulation was restarted. The committor probability p c is here given by the number M β /M of trajectories ending in the β-sheet state, divided by the total number M = M α + M β of trajectories. The different values p c obtained for N = 250 representative configurations of the order parameters A close to A 0 subsequently provide the distribution P(p c ) used to assess the quality of A as a good transition coordinate.

Results and discussion

Equilibrium properties

The variations of the heat capacities of the octa-alanine obtained using the parallel tempering and Wang-Landau Monte Carlo simulations are represented in Fig. 1. The thermally averaged electric dipole µ , end-to-end distance d and radius of gyration R g found for the reference parallel tempering Monte Carlo simulations are also shown in this figure. All heat capacity curves exhibit two peaks, the large peak centered near 215 K being robust for all other WL simulations. At this temperature, the gyration radius and the end-to-end distance sharply increase, as a signature of the unfolding transition.

The low-temperature peak centered near T * = 80 K is characterized by a drop in both µ and d , but a small increase in R g . The most stable conformation of the peptide is an α-helix, which is associated with a high electric dipole µ ≃ 24.5 D resulting from the alignment of intramolecular hydrogen bonds. The helical state has a significant end-to-end distance d ≃ 10 Å, but a reduced R g ≃ 4.6 Å reflecting its compactness. β-hairpin conformations, on the other hand, have a lower dipole and end-to-end distance, but a slightly higher radius of gyration. The stable state of the peptide thus changes from α-helical to β-sheet at the temperature where the heat capacity reaches its first peak. This stable β-sheet state can be further evidenced by calculating other structural indicators, such as the helical fraction or the fluctuations in the overlap function with the native state [49]. The presence of an intermediate β-sheet phase between the α-helical and unfolded random coil states is not specific to the octa-alanine, as longer alanine-rich peptides also exhibit this peculiar behaviour [37]. The β-sheet state is stabilized due to its higher entropy related to softer vibrational modes, whereas α-helices are energetically favored but more tightly bound [37,38]. If the peptide is hydrated, the helical conformation is stabilized by its favorable electrostatic interactions with the polarisable solvent. Most computer simulations of polyalanines in explicit or implicit solvent usually reported a single helix-coil transition induced by temperature [18,[START_REF] Ohkubo | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Yang | [END_REF]24,25], matching the picture of the Zimm-Bragg model [50].

Having a closer examination at the low-temperature heat capacity peak, clear differences between the performances of the three Wang-Landau Monte Carlo simulations are found in comparison with the reference parallel tempering data. While using the (E, µ) order parameters correctly reproduces the location, height and width of this peak, the same cannot be said about the two other flavors (E, d) and (E, R g ). Repeating the Wang-Landau simulations from independent initial conditions, the results obtained with the latter two sets of variables are found to fluctuate significantly, whereas the density of states in E and µ is reproducible. These findings show that sampling with the d or R g variables in addition to E is less efficient than using the electric dipole. This, in turn, already suggests that neither d or R g are reliable order parameters to explore the conformation space corresponding to the α and β states and their transition. The relative efficacy of these order parameters will now be assessed by looking at free energy profiles and computing the committor probabilities.

Electric dipole as the transition coordinate

We first consider the Landau free energy as a function of the order parameter µ. The variations of F (µ) at the temperatures of 60, 80, and 100 K, near the α ↔ β equilibrium, are represented in Fig. 2. Here and in all other free energy calculations, the curves have been shifted such that their minimum is exactly at F = 0. At the three selected temperatures, the free energy curves exhibit two main minima near µ = 9 D and µ = 24.5 D, which correspond to the β-sheet and α-helical conformations, respectively. These minima have the same depth at T = T * , but at 60 K (100 K) the α state (β state) is more stable. The equal free energies of the two states at the transition temperature confirm the first-order character of the α ↔ β interconversion. The barrier between the two states is located in the range µ ∼ 18-20 D, and is rather high (about 9 kcal/mol). Several sets of conformations have been sampled according to their electric dipole. Each set is characterized by its center µ and its width δµ = 0.2 D used as a tolerance factor during the selection. For a given set of conformations, Monte Carlo trajectories were shot until either of the α or β states was reached. In practice, the α-helical state is defined by its electric dipole µ, which must be larger than 22 D. Conformations are considered to be in the β-sheet state if µ ≤ 15 D and d ≤ 7 Å.

Varying µ in the range 17 D ≤ µ ≤ 23 D, the committor probabilities p c are calculated by monitoring the number of trajectories ending up in the β-sheet state earlier than in the α helical state. The averaged committor probability p c is plotted in the upper panel of Fig. 2 as a function of µ. This average value shows a steady decrease from 1 (for µ ≤ 17 D) to 0 (for µ ≥ 22 D). According to our definitions, any trajectory starting at a conformation with µ ≥ 22 D is already in the α state. For dipoles near µ ≃ 20.0 ± 0.1 D, the committor probability is about 0.5, indicating that conformations with this value of µ belong to the isocommittor surface, that is the optimal transition coordinate.

More insight into the quality of the electric dipole as reaction coordinate is provided by looking at the distributions of committor probabilities, shown in Fig. 3 for three values of µ. These distributions are narrowly located near p c = 0 or 1 for µ = 21.0 ± 0.1 D and µ = 18.0 ± 0.1 D, respectively. For µ = 20.0 ± 0.1 D, the distribution is symmetrically peaked near p c = 0.5. This is the expected behavior for a good reaction coordinate at the transition state.

Radius of gyration as the reaction coordinate

The Landau free energy profiles F (R g ) as a function of the gyration radius are displayed in Fig. 4 for the same three temperatures of 60, 80, and 100 K. Only the interval 4.4 Å ≤ R g ≤ 6 Å is shown as it relevant for both the α and β states. The free energies have a minimum at each of these states. Contrary to the free energy curves as a function of the electric dipole, the two minima are not equally deep at the equilibrium temperature. However this is partly compensated by the much broader minimum of the β-sheet state (approximately 4.9 Å ≤ R g ≤ 5.1 Å versus 4.58 Å ≤ R g ≤ 4.62 Å for the helical state). The free energy barrier is close to R g = 4.76 Å but is only moderate, around 1 kcal/mol. This significant difference in the energy barrier when changing the order parameter indicates that important orthogonal variables are missing in the description of the kinetics when selecting conformations based on their gyration radius only [29]. Conformations have been sampled based on their value of R g , for gyration radii in the range 4.66-4.85 Å and δR g = 0.03 Å. The average committor probabilities, represented in the upper panel of Fig. 4 against R g , sharply increase as R g crosses the barrier. The distributions of committor probabilities corresponding to 4.75 Å ≤ R g ≤ 4.78 Å and 4.79 Å ≤ R g ≤ 4.82 Å, for which p c = 0.09 and 0.65, respectively, are shown in Fig. 5.

Both distributions are strongly bimodal with nonzero values only for p c = 0 and p c = 1. Hence the starting conformations of the shooting trajectories already lie in the α or β states, rathers than being actually located inbetween. The gyration radius appears as a very poor reaction coordinate for characterizing the α ↔ β transition.

End-to-end distance as the reaction coordinate

The Landau free energy profiles F (d) obtained using the end-to-end distance d have been represented in Fig. 6, again at the same three temperatures. Conversely to the previous choices of order parameters, more than two minima are found in the variations of the free energy curves at the α/β equilibrium. The stable α-helical and β-sheet states are the deepest minima at the transition temperature of 80 K, they are located near d = 10.2 Å and d = 5.3 Å, respectively. The presence of extra metastable states at d = 8.1 Å and d = 11.5 Å makes it difficult to identify a clear barrier. Furthermore it would be desirable to characterize the conformations corresponding to the metastable minima, especially given that they are separated by rather small barriers (∼ 0.5 kcal/mol), hence being rather probable. In this purpose, sets of conformations at selected values of d within δd = 0.1 Å have been gathered near the three barriers in F (d), namely d = 6.9 Å, d = 9.0 Å, and d = 10.9 Å, as well as near the metastable minima themselves. The average committor probabilities p c calculated after repeating the shooting trajectories are represented in Fig. 6 as a function of d. Quite strikingly, their variations are not monotonic with the order parameter: the chances of ending up in the β-hairpin conformations are high not only for d ≤ 9.0 Å, but also for d ≥ 10.9 Å. This is a straightforward consequence of how d is connected to the conformation. Low d distances usually belong to the β-hairpin state. Torsions along the backbone may break the alignment between the two strands, as seen on the typical conformation found at the barrier near d = 6.9 Å on Fig. 6. Very large values of d correspond to extended conformations, and are thus more likely to belong also to the basin of attraction of the β-sheet state. In the intermediate range 7.0 Å ≤ d ≤ 11.0 Å, conformations can be of the α-helical type with some degree of stretching, or of the β-sheet type with the two strands not perfectly parallel. Therefore the end-to-end distance is a good parameter only for characterizing β-hairpin conformations, but it is not able to discriminate efficiently α helices from extended conformations, which more likely fold into β sheets due to their higher entropy. The distributions P(p c ) of the committor probabilities for which the averages are closest to 50%, that is d = 9.0 Å and 10.9 Å, are displayed in Fig. 7. As was the case for the gyration radius, the two distributions are essentially peaked at p c = 0 and p c = 1, indicating that the conformations sampled according to their end-to-end distance often belong to either of the α or β states. It can be noticed, though, that the two peaks at p c = 0 and 1 are less sharp for d = 9.0 Å. Therefore this set of conformations captures the true transition coordinate less approximately than the gyration radius.

Discussion

The Landau free energies obtained from Monte Carlo simulations show important differences depending on the order parameter considered. While the α-helix and β-sheet stable states appear on these curves as two well-defined minima, the barriers are not always easy to infer. Conformations belonging to either of the α or β states do not differ significantly by their gyration radius, and the barrier between the free energy minima is particularly narrow. In the case of the end-to-end distance order parameter, additional minima corresponding to deformed sheet conformations and multiple barriers blur the picture of a simple two-state system. Only the electric dipole seems to resolve unambiguously helix and strand conformations.

The relative performances of the order parameters µ, R g , and d to describe the α ↔ β interconversion, as assessed by the distributions of committor probabilities at the isocommittor surface, clearly show that only the electric dipole can play the role of a proper reaction coordinate for this transition. The strongly bimodal distributions of the committor probabilities found for the gyration radius and, albeit to a lesser extent, for the end-to-end distance, are related to the proximity of the selected conformations with respect to the α and β states. In contrast, the structures sampled at the top of the free energy barrier in the electric dipole case are not as clearly associated to the end states. These observations can be further rationalized by considering the correlations between the two order parameters remaining when selecting one of them among the three quantities (µ, R g , d). In Fig. 8 these correlations are represented for three order parameters near the isocommittor condition p c = 1/2.

The correlation diagram obtained when selecting conformations based on their electric dipole exhibits one single peak near R g ≃ 4.7-4.8 Å and d = 8.9-9.2 Å. Conversely, conformations selected according to their radius of gyration or their end-to-end distance lead to two peaks which can be easily identified to either of the α or β states. Therefore, neither R g nor d alone provide good reaction coordinates for the α ↔ β transition. On the other hand, the electric dipole seems to be a satisfactory transition coordinate, and the Landau free energy F (µ) at T = 80 K can be used to determine the properties of the transition state.

The location of the barrier, near µ ≃ 20 D, is rather close to the helical state, rather than to the β-hairpin basin. This is reflected on the strong assymetry in the free energy curves F (µ) shown in Fig. 2, in which the α-helix basin is much narrower than the broad β basin. Hence the conformations at the transition state have a high helical fraction. The formation of this helix core from a β-strand conformation, or the breaking of a small part of the full helix, are the rate limiting steps of the transition. On both sides of the barrier, the dipole varies over only about 2 D, with the substantial energy cost of 5-6 kcal/mol to reach the barrier top. This high energy barrier makes the α ↔ β interconversion unlikely to be observed in molecular dynamics simulations, but could well take place in molecular beams with time-of-flights longer than a few microseconds.

Conclusion

As the elementary building blocks of proteins, helices and sheets are of primary importance in determining the folding or misfolding properties and associated functions of many biologically relevant systems. The archetypal polyalanines studied here naturally fold into α helices, but can also form β-hairpin conformations by entropic stabilization. In the present work, the competition between the α and β secondary structures was studied from the perspective of the transition pathways connecting them. Using dedicated Monte Carlo simulations, the α ↔ β interconversion was first characterized by the Landau free energies, which provide a picture of the energy landscape projected onto specific order parameters. The electric dipole µ, the gyration radius R g , and the end-to-end distance d were used to monitor the transition near the equilibrium temperature for octa-alanine as found by our simulations based on the Amber ff96 force field. The variations of the free energy profiles usually show two clear minima associated with each of the α and β states, from which a transition region is easily located at the barrier. However, the end-to-end distance is not such a good coordinate for structural characterization, because many conformations without any helical content share the same value of d as the helix state.

By sampling conformations at the putative transition barrier, trajectories were shot to estimate the committor probability p c to the β-hairpin basin. For all order parameters, suitable ranges were located leading to the isocommittor property p c = 1/2 which is required for the order parameter to be a suitable reaction coordinate. However, the distributions of committor probabilities at the isocommittor surface reveal contrasting behaviors. Most conformations sampled at the barrier in R g or d actually belong to the α or β minima, whereas conformations at the barrier in µ are closer to real transition states. Therefore, only the latter order parameter provides a good reaction coordinate to describe the α ↔ β interconversion. The present results also showed that the transition state has a rather large extent of helicity, which further requires a high potential energy change of more than 5 kcal/mol.

The three order parameters used here were selected based on physical intuition. One possible extension of the present work could be to refine the reaction coordinate following the methods of previous groups [33][34][START_REF] Best | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Peters | [END_REF], particularly the automated schemes proposed by Ma and Dinner [33] or Peters and Trout [START_REF] Peters | [END_REF]. Optimizing the reaction coordinate would also improve the efficiency of sampling methods such as the Wang-Landau algorithm for joint densities of states, because the bottleneck of sampling precisely lies in the crossing of the highest energy barriers. It would also be of interest to investigate larger polyalanines, to determine the extent of helicity in the transition states between α and β conformations. Beyond the competition between these secondary motifs, the folding transition itself could be studied with the same set of methods and dedicated order parameters.

Fig. 1 .

 1 Fig. 1. Upper panel: thermally averaged dipole moment µ , gyration radius Rg , and end-to-end distance d of the Ala8 peptide as a function of temperature, obtained from parallel tempering MC simulations. Lower panel: canonical heat capacity of the same system obtained from MC simulations, or from Wang-Landau simulations using joint densities of states on energy E and µ, Rg, or d as an extra order parameter.

Fig. 2 . 1 Fig. 3 .

 213 Fig. 2. Upper panel: average committor probability pc (µ) of ending first in the sheet conformation, starting from a sample of configurations with prescribed dipole moment µ. Lower panel: Landau free energy profiles along the µ order parameter, at temperatures close to the helix/sheet transition.

Fig. 4 .

 4 Fig. 4. Upper panel: average committor probability pc (Rg) of ending first in the sheet conformation, starting from a sample of configurations with prescribed gyration radius Rg. Lower panel: Landau free energy profiles along the Rg order parameter, at temperatures close to the helix/sheet transition.

Fig. 5 .

 5 Fig. 5. Distribution of committor probabilities pc for ending first in the sheet conformation, for trajectories initiated at conformations with prescribed values of the gyration radius Rg.

Fig. 6 .

 6 Fig. 6. Upper panel: average committor probability pc (d) of ending first in the sheet conformation, starting from a sample of configurations with prescribed end-to-end distance d. Lower panel: Landau free energy profiles along the d order parameter, at temperatures close to the helix/sheet transition.

Fig. 7 .

 7 Fig. 7. Distribution of committor probabilities pc for ending first in the sheet conformation, for trajectories initiated at conformations with prescribed values of the end-to-end distance d.

Fig. 8 .

 8 Fig. 8. Correlation between the two remaining order parameters, for conformations taken at the transition state defined by the reaction coordinate in (a) dipole moment; (b) gyration radius; (c) end-to-end distance.

Acknowledgements

Computer simulations were carried out at the IDRIS and CINES computer centers, which we gratefully acknowledge. We also thank GDR 2758 for financial support. We thank Ph. Dugourd, R. Antoine, and M. Broyer for fruitful discussions.


Thibaut Dumas
email: thibaut.dumas@umontpellier.fr

Bénilde Bonnefille

Elena Gomez

Julien Boccard

Hélène Fenet

Frédérique Courant

Metabolomics approach reveals disruption of metabolic pathways in the marine bivalve Mytilus galloprovincialis exposed to a WWTP effluent extract

Introduction

Conventional wastewater treatment plants (WWTP) discharge a highly diverse range of organic contaminants in aquatic environments [START_REF] Deblonde | Emerging pollutants in wastewater: A review of the literature[END_REF] [START_REF] Loos | EU-wide monitoring survey on emerging polar organic contaminants in wastewater treatment plant effluents[END_REF] . Because of the high demographic pressure, coastlines are affected by contaminant discharges through the development of WWTP sea outfalls, but also by contributions from catchment areas. The health of those important ecosystems could be threatened by contaminant release. In this context, new approaches can be helpful to assess the effects of multicontamination on marine organisms without a priori knowledge. Environmental metabolomics could be one of those approaches since it is able to provide meaningful information on the biochemical response of organisms to a stress [START_REF] Lankadurai | Environmental metabolomics: an emerging approach to study organism responses to environmental stressors[END_REF] . The aim of the present study was to develop a metabolomics approach to highlight key metabolites disrupted by WWTP effluent exposure and then elucidate the biological effects of such exposure on Mediterranean mussels (Mytilus galloprovincialis).

Materials and methods

An exposure has been performed in controlled laboratory conditions. Mussels (n=50, 5 mussels per becher of 2L) were acclimatised for seven days. Then, 25 mussels were exposed to a WWTP effluent extract corresponding to an environmental dilution of 5% while 25 mussels were exposed to solvent (100µL MeOH into 2L of seawater). During the acclimatisation and exposure periods, filtered seawater was continuously aerated and renewed daily (static renewal) and mussels were fed once daily with the marine green alga Tetraselmis suecica at constant density (10,000 cells/mL). Seven days following the exposure, mussels were dissected and their sex determined. No mortality was observed during the exposure period. Metabolic fingerprints were generated from digestive gland by liquid chromatography-high resolution mass spectrometry (control group: n=9 males, n=10 females and n=6 undifferentiated; exposed group: n=10 males, n=10 females and n=5 undifferentiated). Afterwards, univariate and multivariate tests were applied to highlight disrupted metabolites in males and females following exposure. A gender specific response was also assessed quantitatively and qualitatively by an ANOVA Multiblock OPLS (AMOPLS). This method takes advantage from the structured design of multifactorial experiments to decompose and quantify the sources of variability, and investigate main effects (exposure, gender) or interaction(s) (exposure×gender) using dedicated components [START_REF] Boccard | Exploring Omics data from designed experiments using analysis of variance multiblock Orthogonal Partial Least Squares[END_REF] . Moreover, to better understand these observed effects, the WWTP effluent extract was characterised based on a suspect screening approach of 80 contaminants in an attempt to provide a first understanding of the relationship between observed effects and contaminants.

Results and discussion

WWTP effluent effects on male mussels

Treated male mussels showed numerous metabolites altered in response to WWTP effluent exposure. The highlighted metabolites belong mainly to amino acids metabolism (e.g. tyrosine, phenylalanine, leucine, proline, etc.), neurohormones (dopamine and a serotonin metabolite), purine and pyrimidine metabolism (e.g. adenosine, adenine, guanine, uracil etc.), citric acid cycle intermediates (e.g. malate, fumarate), and a component involved in oxidative stress defence (oxidised glutathione). Figure 1 shows an example of altered metabolic pathways, i.e. phenylalanine and tyrosine metabolism. Modulation of these metabolites could have effects on several biological processes, which may lead to a negative impact of organism fitness, such as energy metabolism, DNA, RNA and protein synthesis, osmoregulation, byssus formation and reproduction. The caracterisation of the WWTP effluent extract has allowed the detection of 38 PhACs and their metabolites from various therapeutic classes, and 4 pesticides, occurring in the WWTP effluent extract. Based on previously published literature, some of them can be related to the observed effects, which supports our hypotheses about the impacted biological processes.

Do males and females respond identically?

Female mussels also showed many impacted metabolites. Some of them were common to male, while many others seemed to be specific of female response. Thanks to the AMOPLS analysis, we were able to decompose the variability: the exposure effect (independently of gender) represented 7% (p<0.01) of the total variability, and the exposure×gender interaction, i.e. the gender specific response to WWTP effluent exposure, represented 4% (p<0.05) of the total variability. Meaningful biochemical information was obtained from the components related to each effect. The former highlighted metabolites modulated in response to exposure whatever the gender. The latter highlighted key metabolites involved in gender specific response (for example dopamine modulation was only observed in male in response to the WWTP exposure). This may provide a better understanding of response variabilities to the exposure.

Conclusions

In addition to enriching knowledge about the impact of WWTP effluents in marine environment, our work demonstrates the relevance of a metabolomics approach to highlight key molecular events triggered by chemical exposure, which could lead to adverse outcomes at the individual level. In addition, it underlines that gender has to be taken into account for further research on response of organisms to stress. AMOPLS is a promising approach to better characterise multiple factors of variability in response to a stress.

Figure 1 :

 1 Figure 1: Partial map of phenylalanine and tyrosine metabolism in the response of male mussels to WWTP effluent exposure. Metabolites modulated by exposure are shown in bold. Red arrows correspond to down-modulations and green arrows to up-modulations (exposed compared to control). Direct connections between metabolites are indicated with full arrows.

Acknowledgement -The Author thank the Agence Nationale de Sécurité Sanitaire de l'alimentation, de l'environnement et du travail (AMeCE 2015/1/091) and the Agence Nationale de la Recherche (IMAP ANR-16-CE34-0006-01) for their financial support. The doctoral fellowship of Thibaut Dumas is financially supported by the Agence Nationale de la Recherche (IMAP ANR-16-CE34-0006-01).


Cédric Durand

Fabrice Casset

Bernard Legrand

Marc Faucher

Philippe Renaux

Denis Mercier

Denis Renaud

Didier Dutartre

Eric Ollier

P Ancey

L Buchaillot

Characterization of IN-IC integrable in-plane nanometer scale resonators fabricated by a silicon on nothing advanced CMOS technology

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

INTRODUCTION

Devices are fabricated on 8 inch tools, hence proving the feasibility of MEMS/NEMS with a Front-End process in an industrial clean-room facility. The fabrication is based on a Silicon On Nothing (SON) Technology [START_REF] Monfray | Technical Digest International Electron Devices Meeting[END_REF] under development at STMicroelectronics for future CMOS generations. This work paves the way for the very large scale integration (VLSI) of micro-electro-mechanical resonators for RF applications.

The aim in a near future is to be able to replace spaceconsuming functions such as time references made by off chip Quartz oscillators, by small sizes MEMS/NEMS resonators that could be totally integrated in-IC, resulting in cost, consumption and sizes reductions [START_REF] Nguyen | Proceedings of the 13 th IEEE International Conference On Solid-state Sensors, Actuators and Microsystems Transducers'05[END_REF].

The paper reports on the fabrication process, measurement set-up, electrical equivalent circuit of the device, and gives experimental results such as devices electrical responses, temperature and pressure behavior.

Analytic calculation, finite element and electrical simulations have been performed and fit well measurement results.

DEVICE FABRICATION

Devices were fabricated using the SON technology [START_REF] Monfray | Technical Digest International Electron Devices Meeting[END_REF], to achieve sub-100nm gaps and 400nm thick single-crystal silicon resonators using only pure front-end process [START_REF] Durand | Symposium on Design Test Integration Packaging of MEMS/MOEMS[END_REF] and materials, ensuring in-IC integration capability. The fabrication process flow presented in Figure 1 starts with the patterning of active areas through a thermal SiO 2 layer. A SiGe sacrificial layer is then grown by selective epitaxy, followed by a high boron doped (5.10 19 at/cm 3 ) non-selective single-crystal silicon epitaxy (Figure 1a). E-beam lithography is used to define gaps and the resonator structure leading to a 47nm gap resolution (Figure 1b). The 400nm thick silicon structural layer is then etched by an anisotropic plasma to define air gaps (Figure 1c). Then structures are released by isotropic plasma etching (Figure 1d). The released structures were protected from the metallization by a non-conformal SiO 2 deposition (Figure 1e). Pads are then formed by a NiSi salicidation and aluminum deposition (Figure 1f). A L=4µm w=0.2µm beam with a d=100nm gap is shown in Figure 2. Several dimensions are available to demonstrate multi frequency possibilities on a same chip in the 10-500 MHz range.

MEASUREMENT SET-UP

RF characterizations are performed using a HP8753D network analyzer as depicted in Figure 3. The transmitted signal through the resonator is measured between the vibrating beam and the input port of the device. A bias voltage V P is applied on the input port; the beam and the substrate being grounded to avoid any pull-in effect of the structure with the substrate.

ELECTRICAL EQUIVALENT CIRCUIT

We consider a device based on a clamped-clamped (CC) beam with the following dimensions: L=4.9µm, w=0.26µm and d=95nm. We propose in Figure 4 an equivalent circuit of the device around the first resonance mode of the beam. The circuit will be then used in the ADS software to simulate the frequency response of the device.

Figure 4: Electrical equivalent circuit of the device dynamic response around the resonance frequency obtained on ADS

The electrical model is composed of two 140 Ω access resistances in series, corresponding to the resistance of the electrical pads and the 10µm long doped silicon access lines. The two parallel branches [START_REF] Clark | Proceedings of the 11 th IEEE International Conference On Solid-state Sensors, Actuators and Microsystems Transducers'01[END_REF] correspond respectively to the path of the motional current arising from the mechanical vibration of the beam, and to the path of the parasitic current induced by the direct capacitive coupling between the input and output of the device.

The motional branch is classically described as a series R x L x C x circuit where R x corresponds to the motional resistance. The electrical R x L x C x parameters are linked to mechanical properties of the device by the following relations:

2 . . η Q m k R x = , 2 η m L x = , k C x 2 η = (1) (2) (3) with 2 0 0 d S V p ε η = (4)
m and k are respectively the effective mass and effective stiffness of the resonating beam for the vibration mode under consideration, Q the quality factor, η the electromechanical coupling, V p the bias voltage, ε 0 the permittivity, S the surface, d 0 the initial gap. Simulations were performed using the following values: R x =40.085kΩ,

L x =0.176H, C x =15.894aF.
The second branch is composed of a capacitance. This branch corresponds to the path of the direct coupling parasitic current that is superimposed to the motional current. It induces an anti-resonance peak in the electrical frequency response of the device. The Value for simulation was the following: C 0 =102fF. Figure 5 shows a comparison between the electrical ADS simulation based on the electrical circuit of Figure 4 and the frequency response of the device around the fundamental resonance frequency. The transmission peak is in good correlation with measurement as well as the antiresonance peak shape and frequency. Work is in progress to improve simulations so as to be able to better predict the electrical responses for future devices fabrication.

EXPERIMENTALS RESULTS

Devices dynamic electrical responses

Beam structures with dimensions L=4 to 10µm and w=0.2 to 0.5µm have shown resonance frequencies from 15 to 100MHz for the fundamental mode of vibration under vacuum. Figure 6 gives the transmission and phase response of a clamped-clamped beam (L=4.9µm, w=0.26µm, d=95nm). The fundamental resonance was measured at 95.10MHz for V p =15V. This is in good agreement with analytic calculations that give 94.71MHz and FEM simulations with 93.90MHz. A motional resistance of 40kΩ was extracted from measurements. As quality factor of 2630 under vacuum was calculated from the measured response, considering the bandwidth at the middle height of the peak. The quality factor was also estimated by analytical calculation as detailed in [START_REF] Hao | [END_REF] for the same structure. The analytical calculations give a quality factor of 3015, with a limitation due to support losses. This is in good correlation with measurement. Only a 15% discrepancy between calculation and measurements is observed, probably due to some losses not taken into account in calculations. The fundamental resonance was measured at 67.35MHz for V p =19V. This is in good agreement with analytic calculations that give 66.47MHz and FEM simulations with 67.86MHz. A quality factor Q of 1375 and a motional resistance of 119kΩ are obtained under vacuum. The resonance frequency shifts linearly with V p 2 with a tuning range better than 1%. This enables a frequency tuning for temperature, process or aging compensation for time reference perspectives. An effective gap value of 150nm was extracted from measurements of the resonance frequency shift as a function of V p (curves Figure 7). We also observed third order resonances on a clampedclamped beam (L=10µm, w=0.35µm, d=115nm). The fundamental and third order modes are respectively at 29.7MHz and 154.8MHz (Figure 8), in good agreement with FEM simulations giving 156.9MHz. Figure 9 shows the effects of temperature on both resonant frequency and quality factor for a clampedclamped beam resonator. The study was made on a 150°C temperature range (from -50°C to 100°C), corresponding to specifications for cellular phones.

Dimensions µm

f 0 meas MHz f 0 ana MHz Q meas R m kΩ V p V 4.

Resonators behavior with environment

The resonant frequency shifts down quasi-linearly with temperature, leading to an average of -32ppm/°C, without any compensation. The negative frequency shift when temperature is growing up is attributed to the negative stiffness induced by the beam elongation, implying a compressive stress. We calculated that a V p variation of 2.5V is enough to compensate the frequency shift in a 100°C temperature range. This electrical tuning could be used to compensate the frequency shift for time reference applications.

Figure 9 also shows that the quality factor is growing up with temperature, from 1300 at -40°C to 5500 at 100°C (2630 at 23°C). We think that this tendency could be attributed to the residual stress on the structure at ambient temperature. When temperature is growing up, the induced stress can compensate the residual stress and limit losses. We suppose that a maximum of quality factor can be observed at higher temperatures than 100°C when the induced stress due to temperature will compensate a maximum of residual stress. More measurements should be investigated to verify this hypothesis. Figure 10 shows quality factor dependence on the pressure extracted from measurements. It varies from 1330 under vacuum to 220 at 140 Torr. Even with small beam dimensions, air damping is still the major contributor of energy losses when working at higher pressure than 1 Torr.

CONCLUSION AND PERSPECTIVES

In-IC nanometer scale resonators have been fabricated using SON technology on 8 inch wafers, and successfully characterized up to 154 MHz, with a study of temperature and pressure effects on resonators.

Future work will focus on improving the f.Q product and the low voltage operation of the resonators, improving the output detection by the use of integrated MOSFET transistors [START_REF] Grogg | Proceedings of the 14 th IEEE International Conference On Solid-State Sensors, Actuators and Microsystems Transducers'07 and Eurosensors[END_REF], and working on in-IC integration of devices.

  1a) Active areas after selective SiGe epitaxy followed by a nonselective Si epitaxy 1b) E-beam step leading to a 47 gap resolution 1c) 400nm deep and 88nm wide air gap after anisotropic dry etching 1d) Released silicon structure after SiGe dry isotropic etching 1e) Non-conformal SiO 2 deposition as a protection for metallization on released structures 1f) Metallization pads: NiSi silicide and AlCu stack isolated from the substrate by a thick SiO 2

Figure 1 :

 1 Figure 1: Fabrication process flow and SEM illustrations

Figure 2 :

 2 Figure 2: Perspective SEM picture of a clamped-clamped beam (L=4µm, w=0.2µm and d=100nm)

Figure 3 :

 3 Figure 3: Schematic of the electrical characterization set-up (All tests were performed with devices under a vacuum chamber)

Figure 5 :

 5 Figure 5: Frequency response comparison of electrical ADS simulation and device measurement on a clamped-clamped beam resonator (L=4.9µm, w=0.26µm and d=95nm)

Figure 6 :

 6 Figure 6: Frequency and phase response measured on a clamped-clamped beam resonator (L=4.9µm, w=0.26µm and d=95nm)

Figure 7 :

 7 Figure 7: Frequency response measured on a clampedclamped beam resonator as a function of V p (L=5.0µm, w=0.19µm and d=110nm) Figure 7 gives four resonance curves versus V p for a clamped-clamped beam (L=5.0µm, w=0.19µm, d=105nm).The fundamental resonance was measured at 67.35MHz for V p =19V. This is in good agreement with analytic calculations that give 66.47MHz and FEM simulations with 67.86MHz. A quality factor Q of 1375 and a motional resistance of 119kΩ are obtained under vacuum. The resonance frequency shifts linearly with V p 2 with a tuning range better than 1%. This enables a frequency tuning for temperature, process or aging compensation for time reference perspectives. An effective gap value of 150nm was extracted from measurements of the resonance frequency shift as a function of V p (curves Figure7).

Figure 8 :

 8 Figure 8: Frequency response showing the third resonance mode of a clamped-clamped beam resonator (L=10µm, w=0.35µm and d=115nm)

Figure 9 :

 9 Figure 9: Measured effects of temperature on resonance frequency and Q factor of a clamped-clamped beam resonator (L=4.9µm, w=0.26µm and d=95nm)

Figure 10 :

 10 Figure 10: Evolution of Q as a function of pressure for a CC beam resonator (L=5.0µm, w=0.19µm and d=110nm)

Table 1 :

 1 RF characteristics of the presented devicesTable1resumes the RF characteristics of the devices.

	9 x 0.26 d 0 =95nm	95.10 94.71 2630	40	15
	5.0 x 0.19 d 0 =110nm	67.35 66.47 1375	119	19

ACKNOWLEDGMENT

The authors wish to acknowledge Daniel Grogg and Adrian Ionescu from EPFL for their collaboration on temperature measurements, made in their laboratory.


Johann Mourier
email: johann.mourier@gmail.com

Joachim Claudet

Serge Planes

Human-induced shifts in habitat use and behaviour of a marine predator: the effects of bait provisioning in the blacktip reef shark

Keywords: Provisioning modifies habitat use of sharks wildlife provisioning, shark, coral reef, habitat use, sexual segregation

published or not. The documents may come

Human-induced shifts in habitat use and behaviour of a marine predator: the effects of bait provisioning in the blacktip reef shark

Introduction

Wildlife watching is often presented as a win-win scenario for both people and nature as it can provide new sources of jobs and income to local economies and can enhance environmental awareness [START_REF] Orams | Feeding wildlife as a tourism attraction: a review of issues and impacts[END_REF][START_REF] Knight | Making Wildlife Viewable: Habituation and Attraction[END_REF]. However, the economic viability of commercial wildlife watching depends on the predictable and constant sightings of wild animals [START_REF] Whittaker | Capacity norms on bear viewing platforms[END_REF], which can be enhanced by attracting them with food. These practices have rapidly developed worldwide in different terrestrial and marine ecosystems but have also raised concerns regarding possible negative consequences to the targeted animals and their ecosystems. A number of studies have now demonstrated potential alterations of behavioural and physiological parameters in terrestrial vertebrates, cetaceans and fish, including changes in overall health and increased aggression toward humans (see [START_REF] Brena | Shark and ray provisioning: functional insights into behavioral, ecological and physiological responses across multiple scales[END_REF][START_REF] Gallagher | Biological effects, conservation potential, and research priorities of shark diving tourism[END_REF]Orams, 2002 for reviews). These can act as additional stressors on populations of megafauna that are frequently already endangered or threatened.

For sharks, one of the two taxonomic groups most concerned by marine wildlife watching together with cetaceans [START_REF] Christiansen | Food provisioning increases the risk of injury in a long-lived marine top predator[END_REF], the use of bait is common practice as these species are often sparse or elusive [START_REF] Brena | Shark and ray provisioning: functional insights into behavioral, ecological and physiological responses across multiple scales[END_REF][START_REF] Gallagher | Biological effects, conservation potential, and research priorities of shark diving tourism[END_REF]. An increasing number of studies over the last decade show that provisioning can have a number of impacts on individuals, including on their abundance and residency [START_REF] Clua | Behavioural response of sicklefin lemon sharks Negaprion acutidens to underwater feeding for ecotourism purposes[END_REF][START_REF] Brunnschweiler | Long-Term Changes in Species Composition and Relative Abundances of Sharks at a Provisioning Site[END_REF][START_REF] Kiszka | Using unmanned aerial vehicles (UAVs) to investigate shark and ray densities in a shallow coral lagoon[END_REF], movement patterns and activity space [START_REF] Fitzpatrick | Variation in depth of whitetip reef sharks: does provisioning ecotourism change their behaviour?[END_REF][START_REF] Huveneers | The effects of cage-diving activities on the fine-scale swimming behaviour and space use of white sharks[END_REF]or physiology [START_REF] Barnett | Ecotourism increases the field metabolic rate of whitetip reef sharks[END_REF][START_REF] Huveneers | Interacting with wildlife tourism increases activity of white sharks[END_REF]. In their review, [START_REF] Brena | Shark and ray provisioning: functional insights into behavioral, ecological and physiological responses across multiple scales[END_REF] further highlighted that these most commonly described alterations of individual behaviour can have cascading effects at the group and community scales [START_REF] Drew | Shark-based tourism presents opportunities for facultative dietary shift in coral reef fish[END_REF][START_REF] Meyer | The effects of wildlife tourism provisioning on non-target species[END_REF].

Alteration of the size of individual home range is one of the changes most commonly observed after provisioning sharks or rays in the same location for long periods of time (several years) [START_REF] Huveneers | The effects of cage-diving activities on the fine-scale swimming behaviour and space use of white sharks[END_REF]. Sharks and rays tend to increase their time residency and decrease their activity space following provisioning activities, at least for subsets of the populations [START_REF] Clua | Behavioural response of sicklefin lemon sharks Negaprion acutidens to underwater feeding for ecotourism purposes[END_REF][START_REF] Fitzpatrick | Variation in depth of whitetip reef sharks: does provisioning ecotourism change their behaviour?[END_REF][START_REF] Corcoran | Supplemental Feeding for Ecotourism Reverses Diel Activity and Alters Movement Patterns and Spatial Distribution of the Southern Stingray, Dasyatis americana[END_REF][START_REF] Kiszka | Using unmanned aerial vehicles (UAVs) to investigate shark and ray densities in a shallow coral lagoon[END_REF]. These modifications in space utilization question whether they promote substantive trade-offs in activity budgets and in turn alter energy budget [START_REF] Barnett | Ecotourism increases the field metabolic rate of whitetip reef sharks[END_REF], individual fitness and overall the structuring role these top predators play in their ecosystem. While most studies have shown that provisioning unlikely affects the long-term and overall population behaviour of large species of sharks [START_REF] Hammerschlag | Don't bite the hand that feeds: assessing ecological impacts of provisioning ecotourism on an apex marine predator[END_REF][START_REF] Meyer | The impact of wildlife tourism on the foraging ecology and nutritional condition of an apex predator[END_REF], less evidence is available on smaller species (< 3 m total length) that are potentially more vulnerable to such activities due to restricted movements [START_REF] Maljković | Effects of tourism-related provisioning on the trophic signatures and movement patterns of an apex predator, the Caribbean reef shark[END_REF][START_REF] Kiszka | Using unmanned aerial vehicles (UAVs) to investigate shark and ray densities in a shallow coral lagoon[END_REF].

Here, using a Control-Impact design, we investigated potential changes induced by provisioning in behaviour at the individual and population levels of a common reef shark species. In particular, using acoustic telemetry, we assessed how shark provisioning activities restricted to specific locations for the last two decades can affect habitat use of both females and males for a period of more than two years.

Materials and methods

Study site and species

The study was conducted at Moorea Island (17°30'S; 149°51'W), French Polynesia. The tourism industry has grown rapidly in Moorea [START_REF] Clua | Business partner or simple catch? The economic value of the sicklefin lemon shark in French Polynesia[END_REF][START_REF] Leenhardt | Exploring social-ecological dynamics of a coral reef resource system using participatory modeling and empirical data[END_REF], and has offered activities including interaction with sharks and rays in the lagoon since the 1980's [START_REF] Gaspar | Feeding sites frequentation by the pink whipray Himanturafai in Moorea (French Polynesia) as determined by acoustic telemetry[END_REF][START_REF] Kiszka | Using unmanned aerial vehicles (UAVs) to investigate shark and ray densities in a shallow coral lagoon[END_REF] as well as shark-feeding dives on the outer reef of the North coast [START_REF] Clua | Behavioural response of sicklefin lemon sharks Negaprion acutidens to underwater feeding for ecotourism purposes[END_REF]. In 2004, Moorea authorities implemented a Management Plan for the Marine Environment (Plan de Gestion de l'Espace Maritime) that restricted these shark provisioning activities to two sites on the outer reef in the North (Fig. 1) which primarily sought to attract sicklefin lemon sharks (Negaprion acutidens) and which banned shark provisioning in the lagoon but authorized ray-feeding (Himantura fai) that also ended up attracting large numbers of blacktip reef sharks (C. melanopterus)(site A2 in Fig. 1; [START_REF] Kiszka | Using unmanned aerial vehicles (UAVs) to investigate shark and ray densities in a shallow coral lagoon[END_REF]. Despite local regulations, another former provisioning site is also used, but on a more sporadic basis (site A1 in Fig. 1). Provisioning location A2 is used during the day by both professional operators who can bring up to 50 tourists per boat and individual users who feed pink whiprays and blacktip reef sharks with fish discards and frozen squid in less than 1.5 m depth on sandy banks. Provisioning has been present along the entire northern coast of Moorea since the 1980's and here we considered it as representative of an area potentially impacted by provisioning.

The blacktip reef shark Carcharhinus melanopterus (Quoy & Gaimard, 1824), the target species for the artificial feeding, is one of the most abundant and a common shark species in the coral reefs of the Western Indo-Pacific [START_REF] Vignaud | Blacktip reef sharks, Carcharhinus melanopterus, have high genetic structure and varying demographic histories in their Indo-Pacific range[END_REF], inhabiting shallow reefs and sand-flats of both atolls and high islands [START_REF] Papastamatiou | Scale-dependent effects of habitat on movements and path structure of reef sharks at a predator-dominated atoll[END_REF][START_REF] Speed | Spatial and temporal movement patterns of a multi-species coastal reef shark aggregation[END_REF][START_REF] Mourier | Population structure, spatial distribution and life-history traits of blacktip reef sharks Carcharhinus melanopterus[END_REF][START_REF] Chin | Population organisation in reef sharks: new variations in coastal habitat use by mobile marine predators[END_REF]. Blacktip reef sharks demonstrate a high degree of site attachment and individual spatial overlap [START_REF] Papastamatiou | Scale-dependent effects of habitat on movements and path structure of reef sharks at a predator-dominated atoll[END_REF][START_REF] Mourier | Evidence of social communities in a spatially structured network of a free-ranging shark species[END_REF], with individuals displaying larger range of movements only during the reproductive period (Mourier & Planes, 2013;[START_REF] Speed | Reef shark movements relative to a coastal marine protected area[END_REF]. These sharks have limited home ranges and are resident to specific reefs for periods of at least several years (up to 10 years or more in Moorea; [START_REF] Mourier | Evidence of social communities in a spatially structured network of a free-ranging shark species[END_REF]. Genetic studies demonstrated low connectivity at global [START_REF] Vignaud | Blacktip reef sharks, Carcharhinus melanopterus, have high genetic structure and varying demographic histories in their Indo-Pacific range[END_REF] and even local [START_REF] Vignaud | Microsatellite Analyses of Blacktip Reef Sharks ( Carcharhinus melanopterus ) in a Fragmented Environment Show Structured Clusters[END_REF] scales, highlighting limited and likely very rare large-scale migrations. Blacktip reef sharks have recently been found to demonstrate complex social interactions and assortment by sex at the group level [START_REF] Mourier | Evidence of social communities in a spatially structured network of a free-ranging shark species[END_REF] as well as some patterns of spatial sexual segregation [START_REF] Mourier | Population structure, spatial distribution and life-history traits of blacktip reef sharks Carcharhinus melanopterus[END_REF].

Control-Impact design

We used an array of 6 VR2W acoustic receivers (VEMCO Ltd., Halifax, Canada) deployed from June 2008 to November 2010 to build a Control-Impact design; for each of the three coasts, a pair of receivers was selected at a pass with one receiver in the lagoon and the other on the fore reef (Fig. 1). Receivers were anchored to the substratum using cement-filled car tires with a single metal bar (150 cm tall) through the centre, to which the receiver was attached. This design included two 'Impact' sites on the north coast (i.e., receivers A1 and A2) where provisioning was developed, and two groups of 'Control' sites where no provisioning was developed (i.e., Control 1 on the west coast with receivers B1 and B2; and Control 2 on the east coast with receivers C1 and C2) (Fig. 1). This design allowed us to compare the movements of sharks between the lagoon and the outer-reef for sites affected by provisioning activities and control sites where feeding does not occur. Two controls were chosen to account for influence of spatial effects in individual variability of movements. Range testing of the acoustic receiver array was conducted to determine the distance by which most tag emissions were detected by the receivers. Detection probabilities were found to drop at 400 meters from the receiver. For each detection, the receiver recorded the time, date and transmitter number. We retrieved and downloaded receivers every 6 months, from June 2008 to November 2010, in order to back-up the data, clean them from biofouling and change batteries.

Sampling and tagging

Blacktip reef sharks were caught from a boat using a fishing rod with barbless hooks at multiple locations within and outside the lagoon in each area where acoustic receivers were deployed (Fig. 1; TableS1; Mourier et al., 2013). They were brought alongside the boat where they were inverted and controlled in tonic immobility. Each individual was sexed and total length (TL) was measured from the snout tip to the end of the upper caudal lobe. Sharks were then equipped with VEMCO V16-4H transmitters (68 mm length x 13 mm width; VEMCO, Halifax, Nova Scotia, Canada). These tags continually transmitted an individual coded pulse, randomly within 40 and 80s intervals, producing a nominal transmission around 60 s apart at a frequency of 69 kHz, with an expected battery life of more than 800 days. The tags were externally attached and secured to the dorsal fins. We analysed data from the 35 blacktip reef sharks equipped with acoustic transmitters, ranging from 102 to 157 cm TL (mean ± SD = 132.7 ± 13.0) including 21 males and 14 females (Table S1).

Data analysis

To determine the amount of connectivity between provisioning sites and control sites and to confirm independency of "impact" and "control" sites, we inferred individual daily rate of movements between receivers defined as the number of total movements between each receiver divided by the total number of days of monitoring. Mean daily movement was used to examine the level of connectivity between sites and coasts. A connectivity plot was therefore constructed based on average daily movements between paired receivers indicating the magnitude of incoming and outgoing movements at receivers from the different coasts [START_REF] Heupel | Conservation challenges of sharks with continental scale migrations[END_REF]. The connectivity map was constructed using R package "circlize" package [START_REF] Gu | circlize implements and enhances circular visualization in R[END_REF]) in R v.3.3.0 (R Development Core Team, 2019).

We then assessed the degree of residency of sharks by calculating the proportion of hours each shark spent within the detection range of each receiver. For this, the number of monitoring days was inferred for each shark (i.e., the number of days a shark spent at liberty from tagging date to last recorded detection, with both extreme dates excluded). The degree of residency was calculated for each individual shark where diurnal and nocturnal periods were separated for each day at each receiver. Daily sunrise and sunset hours were used to determine the number of hours of diurnal and nocturnal periods for each day (Source -U.S. Naval Observatory: http://aa.usno.navy.mil/data/index.php).

We then used this measure of residency as the response variable in our Control-Impact analysis.

For this, the residency of sharks from impacted sites represented by a receiver on the ocean side and one on the inside of the lagoon (A1 and A2, respectively; Fig. 1) on the north coast was compared with the residency of sharks from inside and outside the lagoon at two sets of control sites on the west and east coasts (Control 1: B1 on the ocean side and B2 in the lagoon; Control 2: C1 in the lagoon and C2 on the ocean side; Fig. 1).

We used effect sizes to model the differences between impacted and control sites. This approach has proven to be powerful in Control-Impact frameworks (see [START_REF] Claudet | Marine reserves: Fish life history and ecological traits matter[END_REF][START_REF] Claudet | Marine reserves: size and age do matter[END_REF] as examples).

As we were interested in relative changes in residency between provisioning and non-provisioning sites, we used Rijk log-ratios as effect size [START_REF] Hedges | The Meta-Analysis of Response Ratios in Experimental Ecology[END_REF], calculated as follows:

𝑅 𝑖𝑗𝑘 = 𝑙𝑛 ( 𝑋 ̅ 𝐼,𝑖𝑗𝑘 𝑋 ̅ 𝐶,𝑖𝑗𝑘 )
where 𝑋 ̅ 𝐼,𝑖𝑗𝑘 and 𝑋 ̅ 𝐶,𝑖𝑗𝑘 are the mean shark residency indices in habitat i (i.e., lagoon or ocean), sex j (i.e., male or female) and diel phase k (i.e., day or night), in impact (I) or control (C) sites, respectively.

The variance 𝑣 𝑖𝑗𝑘 associated to each effect size was calculated as follows;

𝑣 𝑖𝑗𝑘 = 𝑠𝑑² 𝐼,𝑖𝑗𝑘 𝑛 𝐼,𝑖𝑗𝑘 𝑋 ̅ 𝐼,𝑖𝑗𝑘 + 𝑠𝑑² 𝐶,𝑖𝑗𝑘 𝑛 𝐶,𝑖𝑗𝑘 𝑋 ̅ 𝐶,𝑖𝑗𝑘 where 𝑠𝑑 𝐼,𝑖𝑗𝑘 and 𝑠𝑑 𝐶,𝑖𝑗𝑘 , and 𝑛 𝐼,𝑖𝑗𝑘 and 𝑛 𝐶,𝑖𝑗𝑘 are the standard deviations and sample sizes (i.e., number of sharks detected by the receivers in habitat i, of sex j, during diel phase k), associated to 𝑋 ̅ 𝐼,𝑖𝑗𝑘 and 𝑋 ̅ 𝐶,𝑖𝑗𝑘 , respectively.

We then weighted the effect sizes in our analyses by the inverse of their variance, as follows:

𝑤 𝑖𝑗𝑘 = 1 𝑣 𝑖𝑗𝑘 ⁄
Weighted analyses increase the precision of the combined estimates and increase the power of tests (Gurevitch & Hedges, 1999;[START_REF] Osenberg | Resolving Ecological Questions Through Meta-Analysis: Goals, Metrics, and Models[END_REF].

For each moderator of interest (i.e., habitat, sex and diurnal phase), we calculated weighted averaged effect sizes. All analyses were done using R v.3.3.0 (R Development Core Team, 2019).

Results

The maximum hourly residency, corresponding to the proportion of hours present at the receiver where the individual had the highest detection rate, ranged from 0.05 to 83.45% of hours for a single receiver, with a mean of 15.73 (Table S1). Inter-coast connectivity (i.e., from one side of the island to the other) was low compared to intra-coast connectivity (i.e., from lagoon to open ocean) (Fig. 2) as

shown by the rate of movements between receivers (maximum of 0.015 movements per day). There was no connection between the east and north coasts and only limited North-West and West-East connections (31 recorded N-W movements and only 1 W-E movement).

When considering both habitats together (lagoon or open ocean), time of the day (day or night), and sex (male and female), on average, there was no significant effect of provisioning on mean residency (i.e., the proportion of hour spent at a site) (R = -0.05 ± 0.44, 95% CI). However, there was some heterogeneity and effect of provisioning varied by habitat, time of the day, or sex.

On provisioning sites, both male and female sharks were disproportionately more resident on the fore reef than in the lagoon, when compared to control sites. There were significant 5 to 13-fold increases in residency in the fore reef, in comparison to control sites, not only during the day (R = 2.32 ± 1.78, 95% CI for females; R = 2.29 ± 1.04, 95% CI for males; Fig. 3) but also at night (R = 2.62 ± 2.42, 95% CI for females; R = 1.60 ± 1.10, 95% CI for males; Fig. 3).

Impact of provisioning varied between day and night. For females, loss of lagoon use in provisioning sites, when compared to control sites, was greater at night (R = -4.04 ± 2.12, 95% CI; Fig. 3) than during the day (R = -2.73 ± 1.90, 95% CI; Fig. 3). A decrease in lagoon use by males at provisioning sites followed a similar pattern during the night (R = -2.79 ± 1.04, 95% CI; Fig. 3) but no differences in habitat use was observed between provisioning and control sites during the day (R = 0.34 ± 0.98, 95% CI; Fig. 3).

Discussion

Here, we showed that blacktip reef sharks have high levels of attachment to small areas and low extent of movements at an island scale, and that this pattern of habitat use can be significantly modified by provisioning.

Effects of provisioning were materialized by a shift in habitat use and spatio-temporal dynamics of the monitored sharks. Overall, shark residency was higher during the day at provisioning sites than at control sites, but was lower at night. This is not surprising as provisioning is believed to increase residency of sharks when the activity is conducted [START_REF] Brena | Shark and ray provisioning: functional insights into behavioral, ecological and physiological responses across multiple scales[END_REF]. More importantly, we show that in those provisioning sites, sharks' habitat use was modified, with higher residency on the fore reef than in the lagoon, when compared to control sites. The presence of three main provisioning sites may have increase rate of movement between them and therefore between fore reef and lagoon (Figure 2). Even if one of the provisioning sites is located in the lagoon, the other site outside the lagoon may have attracted to the fore reef at least a portion of the population that would have remained in the lagoon under natural conditions.

Changes in habitat use induced by provisioning were different for males and females. Females lost some residency time in the lagoon at provisioning sites when compared to control sites, not only during the day (when provisioning occurs) but also at night. These patterns likely reflect a shift in habitat use via a translocation of their activity space from lagoon to fore reef habitats when under provisioning influence. While previous studies on the effects of provisioning on sharks have shown an increased residency at provisioning sites [START_REF] Clua | Behavioural response of sicklefin lemon sharks Negaprion acutidens to underwater feeding for ecotourism purposes[END_REF], reduced activity space [START_REF] Huveneers | The effects of cage-diving activities on the fine-scale swimming behaviour and space use of white sharks[END_REF] or change in depth niche during provisioning [START_REF] Fitzpatrick | Variation in depth of whitetip reef sharks: does provisioning ecotourism change their behaviour?[END_REF], no studies have shown evidence of a clear long-term shift in habitat use in any shark species. The presence of three official provisioning sites within a ~10km portion of reef both in the lagoon and fore reef in Moorea may partially explain the increase in time spent by sharks on the fore reef and a decrease in the level of residency in the lagoon. In fact, sharks likely increased movements and exchanges between these sites, simultaneously dividing their spatial use between different habitats.

On the north coast of Moorea, the effects of provisioning on habitat use may induce even larger impacts on natural population dynamics and social structure. Catch-data previously revealed that sex ratio was in favor of males in the fore reef and in favor of females in the lagoon at the island scale.

While the use of the lagoon by males is only weakly affected by provisioning, female habitat preferences shifted from lagoon to fore reef in the north as a consequence of the presence of provisioning sites on the fore reef. These changes may cause patterns of competition between the sexes with a potential increase in sexual harassment by males as previously documented [START_REF] Jacoby | Sex and social networking: the influence of male presence on social structure of female shark groups[END_REF], or may affect natural patterns of spatial sexual segregation found at the island scale [START_REF] Mourier | Population structure, spatial distribution and life-history traits of blacktip reef sharks Carcharhinus melanopterus[END_REF]. In addition, data from northern locations did not show any significant differences in sex-ratios. On the north coast, this lack of spatial sexual segregation was confirmed at the group scale with mixed-sex communities, although there were some assortment preferences for same sex at the level of associations [START_REF] Mourier | Evidence of social communities in a spatially structured network of a free-ranging shark species[END_REF]. Increase in use of the fore reef by females may favor a stronger pattern of spatial overlap between the sexes on the north coast.

A high level of male harassment can drive spatial segregation of the sexes. In many species of elasmobranch, sexual dimorphism in body size or differential activity budget between the sexes may be sufficient to cause very different behavioural strategies and movement patterns [START_REF] Jacoby | Social behaviour in sharks and rays: analysis, patterns and implications for conservation[END_REF]). These differences may select for male sharks to invest more time in pursuit of mates than female sharks that may themselves allocate a higher percentage of time in search of suitable environmental conditions to aid gestation [START_REF] Speed | Heat-seeking sharks: support for behavioural thermoregulation in reef sharks[END_REF]. As a consequence of mating, females often sustain bite marks and serious abrasions to the body and pectoral fins [START_REF] Chin | Blacktip reef sharks (Carcharhinus melanopterus) show high capacity for wound healing and recovery following injury[END_REF] which may in turn favour sexual segregation. Females need to avoid energetically expensive and potentially damaging multiple mating events, a mechanism which leads to sexual segregation. More research is needed to understand whether this change in sexual segregation patterns due to provisioning negatively affects reproduction in blacktip reef shark, and in turn fitness and population dynamics.

Female mobility seems to be increased by provisioning as demonstrated by the observed decrease in residency. Despite increased shark densities during provisioning [START_REF] Mourier | Evidence of social communities in a spatially structured network of a free-ranging shark species[END_REF][START_REF] Kiszka | Using unmanned aerial vehicles (UAVs) to investigate shark and ray densities in a shallow coral lagoon[END_REF], female sharks appear to spend shorter amount of time at provisioning sites in the lagoon, potentially moving across the seascape at higher rates. Previous findings have demonstrated that provisioning can modify activity rate and mobility of sharks [START_REF] Fitzpatrick | Variation in depth of whitetip reef sharks: does provisioning ecotourism change their behaviour?[END_REF][START_REF] Barnett | Ecotourism increases the field metabolic rate of whitetip reef sharks[END_REF][START_REF] Huveneers | Interacting with wildlife tourism increases activity of white sharks[END_REF] and can enhance competitive exclusion with increase intra-and interspecific densities (Juerg M. Brunnschweiler, Abrantes, & [START_REF] Brunnschweiler | Long-Term Changes in Species Composition and Relative Abundances of Sharks at a Provisioning Site[END_REF], but these changes tend to persist at night even if the provisioning activity has ended.

This is surprising and difficult to explain, although it can have important consequences for female sharks. In fact, female reef sharks use shallow warm waters of the lagoons to optimize their gestation [START_REF] Speed | Heat-seeking sharks: support for behavioural thermoregulation in reef sharks[END_REF]. If they are reducing the amount of time spent in the lagoon or increasing their activity, it can affect the time and energy allocated to gestation and could in turn have some implications for the fitness and survival of their pups, especially for a species with a relatively high turnover (1-year reproductive circle; 10-11 year gestation period; Mourier & Planes, 2013). Moreover, food from provisioning activities can induce trophic shifts in fed individuals and, in certain cases, can impact the health and the condition of an animal's body [START_REF] Semeniuk | Using Fatty-Acid Profile Analysis as an Ecologic Indicator in the Management of Tourist Impacts on Marine Wildlife: A Case of Stingray-Feeding in the Caribbean[END_REF][START_REF] Semeniuk | Hematological differences between stingrays at tourist and non-visited sites suggest physiological costs of wildlife tourism[END_REF][START_REF] Maljković | Effects of tourism-related provisioning on the trophic signatures and movement patterns of an apex predator, the Caribbean reef shark[END_REF][START_REF] Brunnschweiler | Hand feeding can periodically fuel a major portion of bull shark energy requirements at a provisioning site in Fiji[END_REF]). Yet, no negative effects have been reported so far on the reproduction of females using provisioning areas (Mourier & Planes, 2013), although further detailed investigations are needed.

Sharks from provisioning sites may allocate less time and energy in search of food and food should therefore be easier to acquire than under natural conditions. While food quantity delivered at the provisioning site in the lagoon is non-negligible [START_REF] Gaspar | Feeding sites frequentation by the pink whipray Himanturafai in Moorea (French Polynesia) as determined by acoustic telemetry[END_REF], the quantity of food that blacktip reef sharks may acquire at the fore reef provisioning locations may be much lower as blacktip reef sharks also have to compete for food with larger sharks such as the sicklefin lemon shark [START_REF] Clua | Behavioural response of sicklefin lemon sharks Negaprion acutidens to underwater feeding for ecotourism purposes[END_REF] and with higher numbers of conspecifics due to increased female's residency.

Nothing is currently known about the quantity of food that sharks may acquire from provisioning activities and what proportion of their diet it represents. Previous work on bull sharks in Fiji found that the amount of provisioned food was unlikely a significant contributor to the daily food budget [START_REF] Abrantes | You are what you eat: Examining the effects of provisioning tourism on shark diets[END_REF], although it can depends on the individual degree of residency (J. M. [START_REF] Brunnschweiler | Hand feeding can periodically fuel a major portion of bull shark energy requirements at a provisioning site in Fiji[END_REF]. Still, the observed changes in activity space may explain observed intraspecific variation in trophic interactions documented in blacktip reef sharks in Moorea [START_REF] Matich | Inter-individual differences in ontogenetic trophic shifts among three marine predators[END_REF]. Therefore, future research is required to better understand the influence of provisioning activities on diet, trophic ecology, metabolism and the health of sharks in Moorea. In particular, it may help us to understand whether food acquired easily from provisioning can offset the loss of energy allocation for reproduction, specifically in females.

The blacktip reef shark is a highly resident shark, like most reef-associated shark species [START_REF] Papastamatiou | Long-term movement patterns and trophic ecology of blacktip reef sharks (Carcharhinus melanopterus) at Palmyra Atoll[END_REF][START_REF] Barnett | Residency and Spatial Use by Reef Sharks of an Isolated Seamount and Its Implications for Conservation[END_REF][START_REF] Bond | Reef Sharks Exhibit Site-Fidelity and Higher Relative Abundance in Marine Reserves on the Mesoamerican Barrier Reef[END_REF][START_REF] Mourier | Extreme Inverted Trophic Pyramid of Reef Sharks Supported by Spawning Groupers[END_REF]. Degree of residency may be exacerbated in remote, isolated islands such as those of the Pacific due to low connectivity and low habitat availability in small islands and atolls [START_REF] Vignaud | Blacktip reef sharks, Carcharhinus melanopterus, have high genetic structure and varying demographic histories in their Indo-Pacific range[END_REF]Mourier et al., 2016). Sharks can remain within small 'sub-habitats' with very limited movements of individuals between sub-habitats, possibly as a result of intra-specific competition [START_REF] Brena | Concede or clash? Solitary sharks competing for food assess rivals to decide[END_REF][START_REF] Papastamatiou | Spatial separation without territoriality in shark communities[END_REF]. The rare movements from one coast to another mostly occurred during the mating period, likely representing reproductive migrations (Mourier & Planes, 2013). The independence between impact and control sites also reinforces the robustness of inferences we made about the influence of provisioning on the degree of residency. All together, these findings also suggest that blacktip reef sharks can be highly susceptible to local changes, including provisioning, as they tend to have restricted movements and hence low abilities to avoid these external pressures.

Wildlife tourism, including shark provisioning, has the potential to contribute significantly to the conservation of animals. However, both tourists and management agencies have the obligation to carefully consider the potential negative effects on health, fitness and long-term behaviour of targeted species and to ensure a sustainable activity by guaranteeing best practice/least impact tourism. If mobility and habitat use is modified for one sex (in our case for female sharks), this activity could have cryptic long-term detrimental impact on the local population, potentially affecting reproduction and population dynamics of resident, long-living fish like reef sharks. Further research is required to determine whether the results of the present study represent an isolated case and to what extend they apply to other localities for blacktip reef sharks and other taxa, as shark and ray provisioning is a popular activity throughout French Polynesia involving many species. In particular, it will be important to address provisioning management from an ecosystem perspective [START_REF] Vignon | Impact of shark-feeding tourism on surrounding fish populations off Moorea Island (French Polynesia)[END_REF][START_REF] Drew | Shark-based tourism presents opportunities for facultative dietary shift in coral reef fish[END_REF][START_REF] Meyer | The effects of wildlife tourism provisioning on non-target species[END_REF] as non-targeted species can also be affected, potentially leading to cascading effects within the ecosystem.

FiguresFig. 1 .

 1 Figures

Fig. 2 .

 2 Fig. 2. Connectivity map indicating rate of movement of individual blacktip reef sharks between

Fig. 3 .

 3 Fig. 3. Effects of provisioning on shark residency, split by habitat type, diurnal activity or sex. Effect

Acknowledgements

This study was implemented thanks to the financial support of the Direction à l'Environnement (DIREN) of French Polynesia, the Coordination Unit of the Coral Reef Initiatives for the Pacific (CRISP Program) and Proscience in French Polynesia. We are grateful to the Centre de Recherche Insulaire et Observatoire de l'Environnement (CRIOBE) staff for their technical support, as well as the volunteers and students who assisted in field data collection.