Nicolas Huybrechts
email: nicolas.huybrechts@cerema.fr

Hassan Smaoui
email: hassan.smaoui@cerema.fr

Sylvain Orseau
email: sylvain.orseau@cerema.fr

Pablo Tassi
email: pablo.tassi@edf.fr

Fabrice Klein
email: f-klein@bordeaux-port.fr

Automatic calibration of the bed friction coefficients to reduce the influence of their seasonal variation: the case of the Gironde estuary

An automatic procedure to identify the bed friction coefficient is tested on a 2D hydrodynamic model of the Gironde estuary (France). The proposed procedure involves an optimization algorithm based on evolution strategy, namely CMA-ES (Covariance Matrix Adaptation Evolution Strategy). Without optimization, application of the same friction distribution to different hydrological conditions leads to significant relative error in water level prediction up to 20-30%. For the tested configuration, 300 runs seemed to be sufficient to reach an optimal value whereas additional 200 runs would help to gain an accuracy of few millimetres (or 0.3%).

In order to reach the same level of accuracy for the different hydrological configurations, it is necessary to adapt for each configuration the bed friction coefficient. Such behaviour tends to confirm a seasonal variation of the friction coefficient and this particularly in the central part of the estuary. Different relationships of the friction coefficient according to the flowrate have been incorporated inside the 2D hydrodynamic model. These relationships effectively allow to maintain an accurate prediction of the water levels close to 10% for a wide range of hydrological configurations.

INTRODUCTION

The tide propagation inside estuary is mainly affected by the modification of the flow section and by energy loses due to bed friction [START_REF] Floch | Propagation de la marée dynamique dans l'estuaire de la Seine et la Seine-Maritime[END_REF]. Converging sections tend to increase the tidal amplitude whereas bottom friction rather decreases this amplitude [START_REF] Floch | Propagation de la marée dynamique dans l'estuaire de la Seine et la Seine-Maritime[END_REF].

In the estuarine uppest part, the interactions between tide and river discharge also impact the tidal propagation [START_REF] Moldwin | Tidal river dynamics[END_REF]. For instance, the flowrate magnitude influences the location of turbidity maximum (TM) and associated mud deposition [START_REF] Sottolichio | Modelling mechanisms for the turbidity maximum stability in the Gironde Estuary, France[END_REF].

The presence of fresh mud deposit induces a modification of the bottom friction (van Rijn 2007) and thus the tide attenuation. In contrast, harmonic analysis is generally used by harbours to predict the water level [START_REF] Moldwin | Tidal river dynamics[END_REF]. Prediction based on harmonic analysis is valid for harbours located near the shore but it becomes less accurate for ports located inside the estuary where interaction between river and tide becomes significant. A typical example is the Port of Bordeaux located 100 km upstream the mouth of the Gironde estuary. The macrotidal Gironde Estuary is located in South-West France covering a surface of 635 km 2 from the Bay of Biscay to 170 km landward (Fig. 1). The estuary is characterized by a complex geomorphology, high turbidity levels up to 20 g.l -1 and a heterogeneous bed composition (Allen 1972[START_REF] Castaing | Le transfert à l'océan des suspensions estuariennes -Cas de la Gironde[END_REF]. Over the years, a large number of hydrodynamic models with different complexity levels have already been developed. These models generally aimed at tracking the turbidity maximum zone [START_REF] Sottolichio | Modelling mechanisms for the turbidity maximum stability in the Gironde Estuary, France[END_REF], Jalon-Rojas et al., 2015) with twodimensional vertical (2DV) or three-dimensional approaches (3D) to compute the hydrodynamics, sediment transport, and salt intrusion. Alternatively, [START_REF] Huybrechts | Optimized predictive 2D hydrodynamic model of the Gironde Estuary (France)[END_REF] proposed a 2D depth-averaged horizontal model (2DH) that showed to be a good compromise between computational cost and accurate solution to efficiently capture the main hydrodynamic processes. Fast and robust models are indeed required in operational tools applied to various alert control systems, including flood control application [START_REF] Laborie | Impact of climate change on Gironde estuary[END_REF] and transport processes, such as sediment matter (Huybrechts andVillaret 2013, Orseau et al. 2020a), or pollutants in the environment. The model developed by [START_REF] Huybrechts | Optimized predictive 2D hydrodynamic model of the Gironde Estuary (France)[END_REF] has been further applied to forecast the ship welcoming capacity inside the Gironde estuary for an interval of 36 hours [START_REF] Orseau | Nav-Tel: an open-source system for ship routing and underkeel clearance management in estuarine channels[END_REF]. [START_REF] Huybrechts | Optimized predictive 2D hydrodynamic model of the Gironde Estuary (France)[END_REF] calibrated the bed friction coefficients by a trial and error procedure in order to reach water level differences lower than 15 cm at the estuary mouth (Verdon, Fig. 1) and at the central part of the estuary (Pauillac, Fig. 1). The calibration and the validation of this model have been performed with field measurements acquired in August 2006 and October-November 2009. These two events are characterized by low river discharges and calm weather conditions. The hydrodynamic model included river and tidal forcing whereas storm surges were not considered. The update of the Huybrechts et al.'s model (2012) to recent bathymetric information coming from up-today bathymetric surveys makes therefore necessary to assess the validity of the previous friction calibration procedure. In contrast to a flood control application where a robust calibration is needed especially for high water levels or storm conditions, a ship route plan requires a robust calibration for a wider range of hydrological condition. Therefore, the accuracy of the model needs to be evaluated under different flow scenarios and weather conditions. Since the trial and error methodology is not suitable to build a friction calibration procedure valid for different hydrological conditions, it is rather proposed to couple the hydrodynamic model with an optimization tool. As discussed by [START_REF] Dung | Multi-objective automatic calibration of hydrodynamic models utilizing inundation maps and gauge data[END_REF], automatic calibration is becoming popular for water-related applications mainly for groundwater, watershed applications. Application of the proposed methodology to large scale and unsteady hydrodynamic model, as observed for estuaries, is still rare [START_REF] Dung | Multi-objective automatic calibration of hydrodynamic models utilizing inundation maps and gauge data[END_REF] due to the required computational resources. The automatic calibration of physical coefficients looks for solution of an inverse problem. This solution corresponds to the minimization of the error between the experimental results (field data) and the results estimated by a numerical model (called direct model). To solve this inverse problem, two different methods have been proposed [START_REF] Fletcher | Practical methods of optimization: Volume 1 constrained optimization[END_REF][START_REF] Fletcher | Practical methods of optimization: Volume 2 unconstrained optimization[END_REF][START_REF] Holland | Adaptation in natural and artificial systems: An introductory analysis with applications to biology, control, and artificial intelligence[END_REF]: gradient-based and meta-heuristic methods.

The first category uses the objective function gradient to search for the optimum, while the second randomly searches for the optimum in a set of solutions (called the population of individuals). The gradient-based methods require that the objective function satisfies regularity conditions (differentiability, convexity). In addition, if the function has several local optima, these methods will be more likely to converge towards a local optimum than a global one. Meta-heuristic methods have been introduced to circumvent the disadvantages of the gradient-based methods. These methods will not use the calculation of the gradient of the objective function, but will explore the global research space based on stochastic processes on a population of individuals rather than on a single individual (solution). Meta-heuristic methods have the advantage of: (i) they are based on a random search and are therefore able to explore the whole space of the solution; (ii) the objective function does not have to be continuous allowing an efficient search for discrete problems, and (iii) they are robust, offering the guarantee of convergence towards the global optimum. However, these methods have the disadvantage of been computational costly at reaching the optimum since they are based on an iterative procedure with slow convergence [START_REF] Rudolph | Convergence analysis of canonical genetic algorithms[END_REF][START_REF] Smaoui | Combination of FEM and CMA-ES algorithm for tranmissivity field identification in aquifer systems[END_REF][START_REF] Smaoui | Transmissivity Identification by Combination of CVFEM and Genetic Algorithm: Application to the Coastal Aquifer[END_REF]. In Geosciences, several meta-heuristic methods have been proposed, e.g. instance genetic algorithms (GA, [START_REF] Goldberg | Genetic Algorithm in Search. Optimization and Machine Learning[END_REF]; simulated annealing (SA, [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF]; particle swarm optimization (PSO, [START_REF] Eberhart | A new optimizer using particle swarm theory[END_REF]; ant colony optimization (ACO, [START_REF] Dorigo | Ant colony system: a cooperative learning approach to the traveling salesman problem[END_REF]; cat swarm optimization (CSO, [START_REF] Ch | Computational intelligence based on the behavior of cats[END_REF], differential evolution (DE, Storn and Price, 1997) and evolution strategy (ES, Baeck et al. 2000aBaeck et al. , 2000b) ) In the present study, the Covariance Matrix Adaptation Evolution Strategy (CMA-ES, Hansen and Ostermeier 1996) is applied to the 2DH hydrodynamic model of the Gironde estuary. This meta-heuristic algorithm is first performed on six hydrological events selected between April and August 2015 with different flowrate values. The period was selected based on the availability of storm surges information provided by Météo-France. From these tests, results are analysed in term of accuracy and friction distribution. Finally, the robustness of the methodology is assessed by considering the effect of the mesh discretization and the number of friction zones.

STUDY AREA

The Gironde Estuary's width reaches 20 km at the mouth and decreases to 3 km downstream the confluence of the Dordogne and the Garonne Rivers. The tidal range varies from 1.5 m during neap tides to 5.5 m during spring tides at the mouth. Both Dordogne and Garonne contributions to the freshwater discharge are estimated to 35% and 65%, respectively [START_REF] Sottolichio | Modélisation de la dynamique des structurs turbides (bouchon vaseux et crème de vase) dans l'estuaire de la Gironde[END_REF]. Based on the bed composition, the estuary can be decomposed in 3 different zones comprising (i) a sandy facies in the estuary mouth; (ii) a mixed facies dominated by mud along the central part and (iii) a fluvial estuary, in the most upstream parts, characterised by the presence of sand, pebbles and gravels (Allen 1972). Fine suspendedsediments observed in the Gironde Estuary compose a pronounced Turbidity Maximum Zone (TMZ) with concentrations ranging between 1 and 20 g/l [START_REF] Sottolichio | A synthesis on seasonal dynamics of highlyconcentrated structures in the Gironde Estuary[END_REF]. Its location along the estuary depends on hydrological conditions [START_REF] Castaing | Le transfert à l'océan des suspensions estuariennes -Cas de la Gironde[END_REF][START_REF] Jalón-Rojas | Turbidity in the fluvial Gironde Estuary (southwest France) based on 10-year continuous monitoring: sensitivity to hydrological conditions[END_REF].

For the year 2015, a harmonic analysis [START_REF] Pawlowicz | Classical tidal harmonic analysis including error estimates in MATLAB using T_TIDE[END_REF] on measured water levels is performed month by month at Verdon (mouth), Pauillac (central part) and Bordeaux (Port) tidal gauge stations (Fig. 1). Variations of the M2 amplitude according to the monthly averaged discharge at previous stations are illustrated on Fig. 2.

At the mouth, a slight increase of the M2 amplitude from 1.44 to 1.53 m with the flowrate is observed. In the central part of the estuary, the M2 amplitude increases progressively until a relatively constant value. Conversely, the M2 amplitude tends to decrease at Bordeaux when flowrate is increasing from 400 to 1200 m 3 /s. A maximum value is reached around 300 m 3 /s.

For the lowest flowrate values, M2 amplitude is then also decreasing probably due to a migration of the turbidity maximum further upstream Bordeaux in the Garonne River (Jalón-Rojas et al. 2018).

MATERIALS AND METHODS

Hydrodynamic model

The hydrodynamics is computed by a two-dimensional formulation based on the solution for the depth-averaged shallow water equations (Eq. 1), with appropriate initial and boundary conditions: The module TELEMAC-2D of the TELEMAC-MASCARET modelling system [START_REF] Hervouet | Hydrodynamic of free surface flows modelling with the finite element method[END_REF]) is applied in this study to solve the shallow water equations (Eq. 1), with the finite element method. The computational domain is comprised from 30 km offshore the estuary mouth to 180 km landward up to the limit of the tidal dynamic and extends to 20 km from the North to the South (Fig. 1). The mesh is unstructured and composed of triangular elements. Two different meshes with different element size resolutions are used in this work: the mesh 1 containing 28000 nodes and the mesh 2 containing 76000 nodes (Fig. 3). The distance between nodes of mesh 1 ranges within 1000 -2000 m offshore 300 m in the central part (Fig. 3a), and within 75-200 m in the tributaries (Fig. 3c). Mesh 2 features an enhanced resolution along the navigation channel: within 300-2000 m offshore, within 60-300 m in the central part (Fig. 3b) and within 33-100m upstream the confluence of both tributaries (Fig. 3d).

𝜕ℎ 𝜕𝑡 + 𝜕(ℎ𝑈) 𝜕𝑥 + 𝜕(ℎ𝑉) 𝜕𝑦 = 0 𝜕𝑈 𝜕𝑡 + 𝜕(𝑈𝑈) 𝜕𝑥 + 𝜕(𝑈𝑉) 𝜕𝑦 = -𝑔 𝜕𝑍 𝜕𝑥 + 1 ℎ 𝜕 𝜕𝑥 (ℎ𝜈 𝜕𝑈 𝜕𝑥 ) + 1 ℎ 𝜕 𝜕𝑦 (ℎ𝜈 𝑡 𝜕𝑈 𝜕𝑦 ) - 𝑔 ℎ 1 𝐾 2 ℎ 1 3 ‖𝑈 ⃗ ⃗ ‖𝑈 + 𝑆 𝑥 𝜕𝑉 𝜕𝑡 + 𝜕(𝑈𝑉) 𝜕𝑥 + 𝜕(𝑉𝑉) 𝜕𝑦 = -𝑔 𝜕𝑍 𝜕𝑦 + 1 ℎ 𝜕 𝜕𝑥 (ℎ𝜈 𝜕𝑉 𝜕𝑥 ) + 1 ℎ 𝜕 𝜕𝑦 (ℎ𝜈 𝑡 𝜕𝑉 𝜕𝑦 ) - 𝑔 ℎ 1 𝐾 2 ℎ 1 3 ‖𝑈 ⃗ ⃗ ‖𝑉 + 𝑆 𝑦
Measured river discharge is imposed at fluvial boundaries for both Gironde tributaries (Fig. 1).

At the maritime boundary, astronomic tide elevation and tidal currents are reconstructed mouth, central part and tributaries. In the present work the number of zones is firstly increased up to 7 zones Ki, (Fig. 1b) to better characterize the bed roughness of the estuary's tributaries. Finally, a configuration accounting for two additional friction zones located at the central part of the estuary is considered (Fig. 1c). The delineation between the mouth and the central part of the estuary (respectively zone 1 and zone 2, Fig. 1b) corresponds to a change in the bed material from sand to mud, respectively. Other remaining delineations are arbitrarily defined mainly based on the geometrical features of the water body.

Optimization algorithm

The optimization algorithm implemented in this work is based on the evolutionary strategy algorithm (ESA, Baeck et al. 2000a, 2000b[START_REF] Dréo | Metaheuristics for hard optimization[END_REF]. According to the Darwin's theory, evolution will produce at the long-term organisms more adapted to their environment [START_REF] Dréo | Metaheuristics for hard optimization[END_REF]. Thus, in order to achieve better results, ESAs evolve in a set (called population) of solutions (called individuals) and a searching root on a random population instead of an individual. Research on a population increases the probability to find the optimum among individuals. During the iterative process (called generation) leading to the optimal solution, the populations evolve according to selection and variation cycles. From the ESA family, we have adopted the CMA-ES algorithm. This algorithm, due to [START_REF] Hansen | Adapting arbitrary normal mutation distributions in evolution strategies: The covariance matrix adaption[END_REF], has been proposed to improve several aspects of the others ESA but specially to overcome the main issues of the optimization solvers based on genetic algorithms [START_REF] Espana | Evolutionary strategies as applied to shear strain effects in reinforced concrete beams[END_REF].

CMA-ES offers good performance in optimizing functions that are not regular enough or even undefined explicitly. The CMA-ES research space has the advantage of evolving real numbers set, thus avoiding the coding/decoding steps that characterize the genetic algorithms (GA).

However, a complete description of the CMA-ES is out of the present scope. It is worth noting that metaheuristic optimization methods such CMA-ES can be effectively coupled with other numerical models to identify some parameters model not accessible from measurements [START_REF] Bayer | Evolutionary algorithms for the optimization of advective control of contaminated aquifer zones[END_REF][START_REF] Elshall | Parallel Inverse Modeling and Uncertainty Quantification for Computationally Demanding Groundwater-Flow Models Using Covariance Matrix Adaptation[END_REF][START_REF] Smaoui | Combination of FEM and CMA-ES algorithm for tranmissivity field identification in aquifer systems[END_REF][START_REF] Smaoui | Transmissivity Identification by Combination of CVFEM and Genetic Algorithm: Application to the Coastal Aquifer[END_REF].

Additional details are provided in the Appendix whereas full descriptions of the algorithm are available in [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF]; Hansen et al. (2003), [START_REF] Dréo | Metaheuristics for hard optimization[END_REF] or Hansen (2006 &2016).

Coupling between the hydrodynamics module and the optimization algorithm

The coupling interface between the optimization algorithm (CMA-ES) and the hydrodynamic module (TELEMAC-2D) is performed with the multi-paradigm numerical computing environment and proprietary programming language Matlab©, developed by MathWorks [START_REF] Moler | A history of MATLAB[END_REF]. The specificity of each application relies on the way of building the objective function. In our application involving 2D hydrodynamic modelling, the unknowns are the values of the different bed friction coefficients and the variable to optimize is the difference between measured and computed water levels. The coupling flowchart between the hydrodynamics module and the optimization algorithm is illustrated on Fig. 4.

An initial distribution of the bed friction coefficient is provided. A steering subroutine is implemented to build the objective function. This subroutine calls the module TELEMAC-2D

for launching the numerical simulations, it post-processes the numerical results and it evaluates the RMSRE (Root Mean Square Relative Error, Eq .4) between the computed water level depending on the friction distribution (Zc, Fig. 4) and the measurement (Zm, Fig. 4). The RMSRE is estimated at 8 tidal gauge stations (Fig. 1): Verdon, Laména, Pauillac, Medoc, Ambes, Bordeaux, Cadillac, Libourne. The first six stations are located along the navigation channel.

Cadillac station located more upstream in the Garonne River, while Libourne station located in the Dordogne River.

The CMA-ES algorithm searches for minimizing the mean value of the RMSRE of the 8 stations.

The minimized value is referred as RMSREm (m for mean between the 8 stations). For each station, the RMSRE is computed by (Eq. 3):

𝑅𝑀𝑆𝑅𝐸 = √ 1 𝑛 ∑ (𝑍 𝑚𝑖 -𝑍 𝑐𝑖 )² 𝑛 𝑖=1 𝑍 𝑚 Eq. 3
Where Zmi is the measured water levels [m] at a gage station and Zci is the computed values, n the number of "i" observations and Zm the mean measured value.

RESULTS

Convergence of the algorithm depending hydrological conditions

Six hydrological events are selected from April to August 2015 with flowrate varying from 150 to 1300 m 3 /s in the Garonne River and from 200 to 1500 m 3 /s in the central part considering the Dordogne contribution. Each event is simulated with 500 runs covering a period of 6 days.

The same initial friction distribution defined from [START_REF] Huybrechts | Optimized predictive 2D hydrodynamic model of the Gironde Estuary (France)[END_REF] is applied to all configurations. The configuration with 7 friction zones is firstly tested. The evolution of the mean error for the 8 stations during the optimization procedure is illustrated on Fig. 5a for the four events from April to June. Similar evolutions are also plotted on Fig. 5b for Pauillac station.

As shown in Fig 5a most of the gain is reached within the first 250 runs. The accuracy gain is more evident with the evolution of the relative error at Pauillac station (Fig. 5b). RMSRE starts around 0.18, then it is decreasing down to lower than 0.1 and it may even reach 0.06 (Fig. 5b).

At Pauillac for a mean tidal range of 1.6m, a decrease of 12% in relative error coincides with an absolute gain of 0.19 m. The results with 7 zones (Fig 1b) are summarized in Table 1.

Relative error is within the range [0.14 -0.18] before optimization and [0.09 -0.12] after optimization (Table 1, Fig. 5a). The global gain is thus within 3 and 8 % and the accuracy gain is increasing with the flowrate. For the 6 hydrological configurations, the mean ratio between final and initial error is 0.7. Along the different station, the mean ratio is almost equal to 1 at Verdon, between 0.51 and 0.59 at Laména, Pauillac and Medoc and within 0.73-0.84 at Ambes, Bordeaux, Cadillac and Libourne. At Verdon, no accuracy gain is observed. It might suggest that improving the accuracy at Verdon through bottom friction coefficient leads to deteriorate the accuracies of the other upstream stations. Improvement at the mouth may probably require enhanced offshore boundary conditions which is a combination of tidal atlases [START_REF] Huybrechts | Optimized predictive 2D hydrodynamic model of the Gironde Estuary (France)[END_REF] and prediction of the storm surges. Less accuracy gain could also be expected at the uppest estuarine part due to a sparse bathymetry dataset. However, for an application related to ship route and underkeel clearance management inside the estuarine configuration [START_REF] Orseau | Nav-Tel: an open-source system for ship routing and underkeel clearance management in estuarine channels[END_REF], it is crucial to attain an efficient prediction of water levels at the central part where navigable depths are more restricted.

Variation of the bed friction distribution related to the flowrate

The algorithm allows to reach a mean error relatively constant for the different hydrological conditions. Nonetheless, it requires for each case an adaptation of the values for bed friction coefficient. As suggested in Fig. 2, the flowrate variation might be responsible of the TM migration of the fluid mud deposits, and it thus has an influence on the bottom roughness.

Friction coefficients are plotted as a function of the flowrate to find a relationship that could be used to set an operational model. It is performed on Figs. 6 for the distribution with 7 friction zones. The zones are gathered as downstream part of the estuary for K1 and K2 (Fig. 6a), as Dordogne river for K3 and K4 (Fig. 6b) and as Garonne river for K5-K6 and K7 (Fig. 6c).

The evolution of the K1 coefficient (mouth, Fig. 6a) is in agreement with the evolution of M2 amplitude in Verdon: slight linear increase according to the flowrate. Similar patterns are also observed between evolution of K2 and M2 amplitude at Pauillac.

For the Dordogne river (Fig. 6b), a parabolic distribution is obtained for K3 and a third polynomial curve accurately describes the K4 evolution. For the Garonne river (Fig. 6c), K5 and K6 describe a second order decreasing curve according to the flowrate. It means that the friction increases due to the seaward migration of the TM which is in line with M2 evolution at Bordeaux. Similarly, for the lowest flowrates, the friction coefficient seems to reach a maximum value as observed with M2 amplitude. The most upstream coefficient K7 describes an inverse behaviour with maximum value around 700 m 3 /s. The plotted regression curves show the general tendencies of the friction evolution. Nonetheless, extrema values for K5 and K6 are not well captured by the simple second order equations. For the operational model, it would be rather suggested to use piecewise linear equations between the 6 optimized values.

Applicability of the methodology to finer mesh discretization and friction distribution

The proposed methodology requires between 250-500 TELEMAC-2D runs to reach the optimized friction distribution. For a 12 cores of 2.4 GHz RAM 48 Go workstation, 500 runs are performed in approximately one-day wall-clock time. Computational efficiently can be gained by avoiding some repetitive steps inside each individual run, by reducing the number of runs or by increasing the computing resources. Nevertheless, alternative way can be suggested to avoid a rough application of the methodology to a finer mesh.

The optimized distribution of bed friction obtained on a 28000 nodes is assessed on a finer grid resolution (76000 nodes). The finer mesh is characterized by a better resolution along the navigation channel and upstream the confluence. Table 2 summarises the averaged RMSRE at each station along the navigation channel obtained with the 6 hydrological events after an optimization of the coarser mesh. It reaches values ranging from 7.7% (Pauillac) to 13%

(Bordeaux) with a mean value of 9.8 %. Direct application of the optimized values to the finer grid leads to RMSE within 6.9 to 11% with a mean value of 8.7%. It means that the optimized values and associated abacuses are also valid on this finer mesh. As an alternative, the optimized value obtained from the coarse mesh could be used as initial solution for a second optimization with a finer mesh and a shorter number of runs.

To address the sensitivity to the number of friction zone, the optimization methodology is applied to the same 6 hydrological events, but with 9 zones (Fig. 1c) on Mesh 1. As detailed in Table 1, no significant differences can be noticed in term of accuracy. However, the values of the friction distributions are different in the central part. To distinguish the friction zone between the methodologies accounting for 7 or 9 zones, friction coefficients for the 9 zone distributions are noted as KK1 to KK9. In fact, K2 extension covers the area sum of KK2, KK3

and KK4 whereas KK1 is the same zone as K1. The evolution of the bed friction coefficients in the central part is shown in Fig. 7.

KK1 describes a linear relationship whereas KK2 and KK3 describe a parabolic relationship with a maximum value around 800 -900 m 3 /s. KK4 rather describes a parabolic relationship with a minimum value. It should be noted that KK1 is smaller than K1 and KK2 higher than K2. It results in a more abrupt transition occurring between the two zones which may affect the numerical results if the model is coupled to a sediment transport and bed evolution module.

Application of time-varying friction coefficients to a medium-term simulation (6 months).

From Fig. 6, a relationship can be built between the friction coefficient and the flowrate for all For the simulation 2 referred as "mean", the values of friction coefficients are not suited for low flowrates. The RMSRE can increase up to more than 20% (> 30 cm) at Pauillac. A similar behaviour is observed at Bordeaux. In contrast, for the simulation 1 referred as "Low", the accuracy tends to decrease at Pauillac once values are higher than 500 m 3 /s, while, at Bordeaux, the accuracy is more variable. However even if the prediction is correct at Bordeaux, the accuracy is not sufficient in the central part. The advantage of the PWL simulation is to maintain a constant accuracy for a wider range of flowrate since it combines the advantage of the two previous simulations.

  h is the water depth [m], 𝑈 ⃗ ⃗ is the depth-averaged flow velocity vector [m/s] , with eastwest, north-south components U and V, respectively, ‖𝑈 ⃗ ⃗ ‖|| is the velocity norm, g is the gravity acceleration [m²/s], Ζ is the free surface elevation [m], νt is the momentum diffusion coefficient [m²/s], ρ is the water density [m 3 /kg], K is the Strickler-Manning coefficient [m 1/3 /s], Sx and Sy are additional source terms. The Strickler coefficient used for the bed friction is just the inverse of the Manning coefficient. The mathematical system is therefore composed of 3 equations and 5 unknowns (h, U, V, K and ν). Bed friction and diffusion coefficients (K and ν) are provided by additional closure relationships or imposed values. A constant value equal to 1 m²/s is imposed for the diffusion coefficient over the whole numerical domain. In the shallow water equations, the bed friction term is included in the source term of the momentum equation.

  using NEA tidal atlases(North East Atlantic, Pairaud et al 2008[START_REF] Huybrechts | Optimized predictive 2D hydrodynamic model of the Gironde Estuary (France)[END_REF] as a superposition of harmonic waves[START_REF] Schureman | Manual of harmonic analysis and prediction of tides[END_REF] for each of the nodes of the offshore boundary (Eq. 2).𝐻 𝑡𝑖𝑑𝑒 = 𝐻 0 + 𝑑𝐻 0 (𝑡) + ∑ 𝐻 𝑛 𝑓 𝑛 𝑐𝑜𝑠(𝜎 𝑛 𝑡 -𝑔 𝑛 + 𝑉 𝑛 -𝑢 𝑛 ) 𝑛 Eq. 2 where Htide = the tidal height; Ho = the mean height of the water level; n = the harmonics number; Hn = the mean amplitude of the n-wave; f n = the nodal correction for the amplitude; σn = the frequency; t = the time; gn = the phase lag of the equilibrium tide; Vn = the astronomic argument; and un = the nodal correction for the phase lag. dH0 is the storm surge contribution.Sea levels variation due to storm surges are applied to the tidal signal to improve water level predictions. Storm surge data are provided by a Météo-France model and computed every 10 minutes at 12 nodes located along the maritime boundary. Linear interpolation is then performed to incorporate surge values for each boundary node. In the previous study,[START_REF] Huybrechts | Optimized predictive 2D hydrodynamic model of the Gironde Estuary (France)[END_REF] decomposed the bed friction into 4 different zones delineated as:

  the 7 friction zones. Concerning the regression, piecewise linear relationships are selected to interpolate the values for all flowrate values. Three simulations are conducted from April 1 st to the end of October 2015. The first two simulations are based on steady friction coefficients extracted from the optimization step. The first one corresponds to friction coefficient representative of low flowrate configuration (200 m 3 /s) and the second one rather corresponds to a configuration representative of mean discharge configuration (800 m 3 /s). The last and third simulation tests the piecewise linear relationships (PWL). Time series of storm surges and flowrate are imposed at the boundary conditions. The time step of flowrates is equal to 2 hours whereas it is equal to 10 min for the storm surge. Flowrate ranges within 130 and 2700 m 3 /s during this period. At Verdon, Pauillac and Bordeaux, the RMSRE are evaluated every two tidal cycles (25 hours) to provide an averaged estimation of the accuracy. Values of RMSRE at Bordeaux and Pauillac stations are plotted in regards to flowrate also averaged every 25 hours. At the Verdon station, the accuracy for each simulation is equivalent (not showed here).

Figure 1 Figure 2 Figure 3

 123 Figure 1 Model extension (1a) and bottom friction zones with 7 zones (1b) or 9 zones (1c)

Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 .

 45678 Figure 4 Flowchart of the coupling between the hydrodynamics module and the optimization algorithm

Table 1 RMSREm for the different configurations for different flowrate (Garonne and Dordogne contribution)

 1

	Total	RMSREm -7 zones	RMSREm -9 zones
	flowrate
	(m3/s)
		Before	After	Before	After
		optimization	Optimization	optimization	optimization
	187	0.149	0.119	0.152	0.121
	205	0.142	0.123	0.142	0.122
	450	0.157	0.117	0.153	0.116
	790	0.176	0.122	0.167	0.118
	883	0.183	0.107	0.173	0.111
	1490	0.173	0.091	0.163	0.089

Table 2 Reached RMSE for two different mesh sizes

 2

	Mesh	Mean	Verdon	Laména	Pauillac	Medoc	Ambes	Bordeaux
		RMSE
	M1 -28000	0.098	0.089	0.094	0.077	0.106	0.092	0.13
	nodes
	M2 -76000	0.087	0.090	0.070	0.069	0.100	0.085	0.11
	nodes

DATA AVAILABILITY

Data, models, and code scripts used for coupling CMA-ES and TELEMAC-2D developed in this study are available from the corresponding author upon request. Water levels data are available at www.vigicrues.gouv.fr.

CONCLUSIONS AND PERSPECTIVES

The CMA-ES algorithm has been coupled to the hydrodynamic module TELEMAC-2D applied to optimize the distribution of the bed friction coefficient inside the Gironde estuary. For the tested configuration, 300 runs seemed to be sufficient to reach an optimal value. Additional 200 runs would help to gain an accuracy of few millimetres (or 0.3%). For simulations performed on a 12 core workstation), 500 runs are completed in approximately one-day wallclock time for 12 tidal cycles. The application of the proposed methodology shows that it is necessary to modify the bed friction coefficient in order to reach the same level of accuracy for the different hydrological configurations. It also confirms a seasonal variation of the friction coefficient and this particularly in the central part of the estuary. Different relationships of the friction coefficient according to the flowrate have been incorporated inside the operational model. These relationships effectively allow to maintain an accurate prediction of the water levels for a wide range of hydrological configurations. However, further investigations on more extreme events, such as flood, storm and long dry periods, are still needed to provide more robust bed friction relationships.

For operational models, it would be interesting to further apply the methodology with several flow configurations in order to build a surrogate model providing the friction distribution according to hydro-meteorological forcing (flowrate, tidal range, storm surge) and to compare such variation to data related to the bed texture or water column as bed sample or satellite images of suspended matters.

ACKNOWLEDGMENTS

The research leading to these results has received funding from the Connecting Europe Facility Like all meta-heuristic algorithms, CMA-ES starts from an initial population randomly chosen.

To build a new generation of individuals, the CMA-ES algorithm follows on from the selection step in which the new candidate solutions are sampled using a multivariate normal distribution. Then individuals of this generation are evaluated via the objective function and selected according to their fitness (or objective function value) to be part of the next generation. Then comes the recombination stage to select a new mean value for the distribution. The penultimate step of the CMA-ES algorithm is the mutation which consists in adding a random vector acting as a perturbation with zero mean. The adaptation step terminates the algorithm by updating the various parameters involved in the construction of the covariance matrix. From this brief description of CMA-ES algorithm, we conclude that it is the mutation and adaptation stages that make this algorithm a robust and powerful tool for complex numerical optimization. In order to not burden the text and given their importance, we will briefly describe these two stages. 𝑔+1) is the ith individual of the population of the generation (𝑔 + 1), 〈𝑝〉 𝜇 (𝑔) is the mean value of 𝑝 (𝑔) at generation (𝑔) computed by

Mutation

𝜎 (𝑔) is the standard deviation at generation (𝑔) (but for CMA-ES, it is also called step size), 𝑁(0, 𝐶 (𝑔) ) note the normal distribution with center 0 and covariance 𝐶 (𝑔) at generation (𝑔).

It should be noted that the covariance matrix is a symmetric definite positive matrix, therefore diagonalizable. In this case the covariance matrix 𝐶 (𝑔) can be written as:

𝐶 (𝑔) = 𝐵 (𝑔) 𝐷 (𝑔) (𝐵 (𝑔) 𝐷 (𝑔) ) 𝑇 Eq. 5 Where the columns of the matrix 𝐵 (𝑔) are exactly the eigenvectors of 𝐶 (𝑔) and 𝐷 (𝑔) is a diagonal matrix whose diagonal elements are the square root of the eigenvalues of 𝐶 (𝑔) . The combination of expressions ( 1) and ( 2) allows to rewrite (1) in the new form as:

𝑝 𝑖 (𝑔+1) = 〈𝑝〉 𝜇 (𝑔) + 𝜎 (𝑔) 𝐵 (𝑔) 𝐷 (𝑔) 𝑧 𝑖 , 𝑖 = 1, … . . 𝜆

Eq. 6

With 𝑧 𝑖 = (𝐵 (𝑔) 𝐷 (𝑔) ) 𝑇 𝑁(0, 𝐼) , 𝑖 = 1, … . . 𝜆 Finally, the calculation of the covariance matrix at generation (𝑔 + 1) is based on the calculation of the evolution of the path 𝑝 𝑡 at generation (𝑔 + 1) according to the following scheme: 𝑔+1) -〈𝑝〉 𝜇 (𝑔) )

Eq. 7 𝐶 (𝑔+1) = (1 -𝑐 𝑐𝑜𝑣 ). 𝐶 (𝑔) + 𝑐 𝑐𝑜𝑣 𝑝 𝑡 (𝑔+1) . (𝑝 𝑡 (𝑔+1) ) 𝑇 Eq. 8

Where

is the cumulative time of the evolution path. The parameter 𝑐 can be interpreted a weight allowing the smoothing of 𝑝 𝑡 and can be normalized by 𝑐 𝑢 = √𝑐(2 -𝑐).

1 𝑐 𝑐𝑜𝑣 denotes the average time for the covariance matrix. In the other word, 𝑐 𝑐𝑜𝑣 allows the updating of the covariance matrix and it can be considered as the learning rate.

BIBLIOGRAPHY

Allen, G.P., (1972). Etude des processus sédimentaires dans l'estuaire de la Gironde. PhD Thesis, University of Bordeaux I, 310.


Dantas Natanael

Marcelo Farias

Gabriel Machado Falchetti

Paulina Matos

Cairé Schmitt

Nicolas Barreto

Jean-Luc Argenta

Evelyne Rolland

Bachère

Maria Luciane

Rafael Diego Perazzolo

Rosa

Natanael Dantas Farias

Marcelo Falchetti

Gabriel Machado Matos

Paulina Schmitt

Cairé Barreto

Nicolas Argenta

Jean-Luc Rolland

Evelyne Bachère

Luciane Maria Perazzolo

Rafael Diego Rosa
email: rafael.d.rosa@ufsc.br

Litopenaeus vannamei stylicins are constitutively produced by hemocytes

Keywords:

come

Introduction

Antimicrobial host defense peptides (AMPs) are important components of the innate immune system of both vertebrates and invertebrates. They are usually described as gene-encoded peptides (less than 10 kDa) with cationic and amphipathic properties which selectively target the negatively charged membranes of microbes [START_REF] Brogden | Antimicrobial peptides: pore formers or metabolic inhibitors in bacteria?[END_REF]. In addition to those classical cationic antimicrobial peptides (also known as CAMPs), the current classification of AMPs also includes polypeptides/proteins larger than 10 kDa, AMPs generated by the processing of precursor molecules and anionic peptides [START_REF] Yount | Advances in antimicrobial peptide immunobiology[END_REF]. Anionic antimicrobial peptides (AAMPs) comprise a non-phylogenetic group of either geneencoded or nonribosomally synthesized molecules with a high proportion of anionic amino acid residues (aspartate and glutamate). AAMPs are widely distributed in living organisms and play an important role in host defense against bacteria, fungi and viruses [START_REF] Harris | Anionic antimicrobial peptides from eukaryotic organisms[END_REF]. Like their cationic counterparts, AAMPs are multifunctional molecules engaged in different biological and immunological processes beyond antimicrobial functions [START_REF] Harris | Anionic antimicrobial peptides from eukaryotic organisms[END_REF].

Shrimp farming is an important economic activity for many developing countries in Asia and Latin America, which has been repeatedly threatened by infections caused by viruses and pathogenic bacteria from the genus Vibrio. Consequently, infectious disease outbreaks are clearly a major concern in aquaculture that has encouraged extensive research efforts. The scientific findings in the last decade have provided valuable information on the role of AMPs in shrimp defenses. More than natural antibiotics, shrimp AMPs are also involved in the control of the natural microbiota, wound healing, bacterial clearance and other immunomodulatory functions [START_REF] Destoumieux-Garzón | Antimicrobial peptides in marine invertebrate health and disease[END_REF]. To date, four gene-encoded AMP families have been identified in the hemocytes of penaeid shrimp: penaeidins, crustins, antilipopolysaccharide factors (ALFs) and stylicins [START_REF] Destoumieux-Garzón | Antimicrobial peptides in marine invertebrate health and disease[END_REF][START_REF] Rosa | Antimicrobial peptides in crustaceans[END_REF].

Stylicins were initially identified as transcripts associated to shrimp survival to pathogenic Vibrio infections [START_REF] De Lorgeril | Identification of genes that are differentially expressed in hemocytes of the Pacific blue shrimp (Litopenaeus stylirostris ) surviving an infection with Vibrio penaeicida[END_REF]. Characterization of stylicins revealed that they are anionic peptides of 8.9 kDa composed of an Nterminal proline-rich region followed by a C-terminal region containing 13 cysteine residues [START_REF] Rolland | Stylicins, a new family of antimicrobial peptides from the Pacific blue shrimp Litopenaeus stylirostris[END_REF]. In Litopenaeus stylirostris shrimp, stylicins form a diverse AAMP family composed of two members, named Lsty-Stylicin1 and Lsty-Stylicin2. Although the antimicrobial activity of the recombinant rLsty-Stylicin1 has been shown to be restricted to filamentous fungi, the rLsty-Stylicin1 displayed a strong lipopolysaccharide (LPS)-binding activity and the ability to agglutinate Gram-negative bacteria [START_REF] Rolland | Stylicins, a new family of antimicrobial peptides from the Pacific blue shrimp Litopenaeus stylirostris[END_REF]. In the kuruma prawn Marsupenaeus japonicus, the gene expression of its single stylicin (Mjap-Stylicin) showed to be modulated in gills and hepatopancreas in response to the White spot syndrome virus (WSSV) [START_REF] Liu | Identification and expression analysis of a novel stylicin antimicrobial peptide from Kuruma shrimp (Marsupenaeus japonicus)[END_REF]. Curiously, apart from these two reports [START_REF] Rolland | Stylicins, a new family of antimicrobial peptides from the Pacific blue shrimp Litopenaeus stylirostris[END_REF][START_REF] Liu | Identification and expression analysis of a novel stylicin antimicrobial peptide from Kuruma shrimp (Marsupenaeus japonicus)[END_REF], no other stylicin members have been characterized thus far.

In order to fill this research gap, novel members of the stylicin AMP family were identified and characterized in the most important cultured penaeid species, Litopenaeus vannamei (Lvan-Stylicin1 and Lvan-Stylicin2). We showed that Lvan-Stylicins (also known as Vibrio penaeicida-induced cysteine and proline-rich peptides or LvVICPs [START_REF] Wang | Analysis of expression, cellular localization, and function of three inhibitors of apoptosis (IAPs) from Litopenaeus vannamei during WSSV infection and in regulation of antimicrobial peptide genes (AMPs)[END_REF]) are highly anionic peptides encoded by two distinct genomic loci that follow different patterns of gene regulation during shrimp development and after microbial infections. Interestingly, while both genes responded to danger/damage-associated molecular patterns (shrimp muscle tissues), only Lvan-Stylicin2 showed to be up-regulated in circulating hemocytes in response to a Vibrio infection. Moreover, by combining immunohistochemistry and whole-mount immunofluorescence assays, we showed that Lvan-Stylicins are also constitutively produced by the midgut columnar epithelial cells. To our knowledge, this is the first evidence for the expression of a shrimp geneencoded AMP in other tissues than the immune cells from the hemolymph.

Materials and methods

Animals, immune challenge and tissue collection

Juvenile Pacific white shrimp (Litopenaeus vannamei) (10 ± 2 g) were obtained from the Laboratory of Marine Shrimp of the Federal University of Santa Catarina (Southern Brazil). After an acclimation period of seven days, animals (n=5) were stimulated by the injection of 5×107 colonyforming units (CFU)/animal of heat killed (70 °C for 20 min) Vibrio harveyi ATCC 14126 in 100 μL sterile seawater (SSW). Naïve (unchallenged) animals (n=5) were used as control. At 48 h post-stimulation, hemolymph was withdrawn into modified Alsever solution (MAS: 27mM sodium citrate, 336mM NaCl, 115mM glucose, 9mM EDTA, pH 7.0) and hemocytes were separated from plasma by centrifugation (800×g for 10 min at 4 °C). Shrimp were subsequently anesthetized (ice bath for 10 min) and sacrificed for the collection of the following tissues: gills, muscle, nerve cord, hepatopancreas, foregut, midgut and hindgut. Tissue samples were washed in Tris-saline solution (10mM Tris, 330mM NaCl, pH 7.4), homogenized in TRIzol reagent (Thermo Scientific) and immediately processed for total RNA isolation and tissue distribution analysis. Pleopods from naïve shrimp were collected and processed for genomic DNA (gDNA) extraction.

Experimental infections

Two unrelated shrimp pathogens were chosen for experimental infections, the Gram-negative Vibrio harveyi and the White spot syndrome virus (WSSV). For the bacterial infection, 6×107 CFU/animal of live V. harveyi ATCC 14126 in 100 μL SSW (median lethal dose within 2 days, LD50/2) or 100 μL SSW (injury control) were injected. For the viral infection, shrimp were injected with 100 μL of a WSSV inoculum containing 3×102 viral particles (median lethal dose within 15 days, LD50/15). The WSSV inoculum was prepared from muscle tissues of WSSV-infected shrimp as previously described [START_REF] Goncalves | Molecular signatures at imminent death: hemocyte gene expression profiling of shrimp succumbing to viral and fungal infections[END_REF]. Control animals for the viral infection were injected with 100 μL of a muscle tissue homogenate prepared from WSSV-free shrimp. At 48 h post-infections, circulating hemocytes and midguts were collected, pooled (3 pools of 5 animals per condition) and processed for total RNA extraction and quantitative PCR analysis of gene expression. Unchallenged animals (naïve shrimp at time 0 h) were used as control for all experiments.

In a second experiment, individual shrimp (n=3) were challenged by the oral administration of 7.5×105 CFU/animal of live V. harveyi ATCC 14126 in 50 μL SSW. Shrimp (n=3) that received SSW served as controls. The bacterial per os challenge was performed as previously described [START_REF] Pilotto | Exploring the impact of the biofloc rearing system and an oral WSSV challenge on the intestinal bacteriome of Litopenaeus vannamei[END_REF]. No mortalities were recorded during the course of the experiment. At 21 h post-challenge, midguts from both Vibriochallenged and unchallenged shrimp were collected and flushed with cold Tris-saline solution and then cut into three equal portions (anterior, middle and posterior). Each individual sample was homogenized in TRIzol reagent (Thermo Scientific) and processed for total RNA extraction and quantitative PCR analysis of gene expression.

Genomic DNA and total RNA extraction and cDNA synthesis

For gDNA extraction, individual pleopods were homogenized and incubated at 55 °C for 1 h in 500 μL of lysis buffer (100mM Tris-HCl pH 8.5, 100mM NaCl, 50mM EDTA pH 8.0, 1% SDS, 0.25 μg/μL proteinase K). After addition of 3M potassium acetate (1:2; v:v), samples were incubated at 4 °C for 30 min and centrifuged at 14,000×g for 10 min. Following precipitation with isopropanol, gDNA samples were washed in 70% ethanol and treated with 50 μg/mL RNase A (Fermentas) at 37 °C for 30 min.

Quantification and quality of gDNA samples were assessed by spectrophotometry and 0.8% agarose gel electrophoresis,respectively.

Total RNA was extracted using TRIzol reagent (Thermo Scientific) according to the manufacturer's instructions. RNA samples were treated with DNase I (Thermo Scientific) at 37 °C for 15 min and precipitated with 0.3M sodium acetate (pH 5.2) and isopropanol (1:1; v:v). RNA amount and quality were assessed by spectrophotometric analysis and the integrity of total RNA was analyzed by 0.8% agarose gel electrophoresis. Following heat denaturation (70 °C for 5 min), reverse transcription was performed using 1 μg of purified total RNA with 50 ng/μL oligo(dT)12-18 in a 20-μL reaction volume containing the RevertAid Reverse Transcriptase (Thermo Scientific), according to the manufacturer's instructions.

Molecular cloning

PCR amplifications for molecular cloning were conducted using primers based on the nucleotide sequence of two stylicin homologues (contigs: DN31608_c0_g1_i1 and DN31608_c0_g1_i2) identified in midgut transcriptomes of L. vannamei (unpublished data). PCR reactions were carried out in a 15-μL reaction volume containing 50-100 ng of gDNA, 2mM MgCl2, 0.4mM dNTP Mix, 0.4 μM of each primer (Table 1) and 1 U Taq DNA Polymerase (Sinapse). PCR conditions were as follows: 1 cycle of denaturation at 95 °C for 10 min followed by 35 cycles of 95 °C for 45 s, 55 °C for 45 s and 72 °C for 2 min, and a final extension step of 72 °C for 10 min. The amplification products were analyzed by electrophoresis (1.5% agarose gel) and cloned into a pCR2.1-TOPO vector (Thermo Scientific). The positive clones were identified by colony PCR and plasmid sequencing.

Sequence data analysis and phylogeny

Stylicin sequences from the penaeid species L. vannamei (contigs: DN31608_c0_g1_i1 and DN31608_c0_g1_i2), L. stylirostris (Lsty-Stylicin1: EU177435; Lsty-Stylicin2: EU177437) and M. japonicus (Mjap-Stylicin: KR063277) were used for the search of homologous sequences in the following publicly accessible databases: Expressed Sequence Tags (EST), Transcriptome Shotgun Assembly (TSA) and Whole-Genome Shotgun Contigs (WGS). Only full-length coding sequences were included. Homology searches were performed using tBLASTX at the NCBI web servers (http://www.ncbi.nlm.nih.gov/BLAST). All nucleotide sequences were manually inspected and translated using the ExPASy Translate Tool (http://web.expasy.org/translate/). Prediction of signal peptide was performed with the SignalP 4.1 program (http://www.cbs.dtu.dk/services/SignalP/) and the theoretical isoelectric point (pI) and molecular weight (MW) of the mature peptides were predicted using the ExPASy ProtParam Tool (http://web.expasy.org/protparam/). Phylogenetic analysis based on both nucleotide and predicted amino acid sequences were conducted in MEGA X [START_REF] Kumar | MEGA X: molecular evolutionary genetics analysis across computing platforms[END_REF] using the Maximum Likelihood method. Bootstrap sampling was reiterated 1000 times using a 50% bootstrap cutoff.

Semiquantitative RT-PCR analysis

The extraction of total RNA and cDNA synthesis were performed using the method described above. PCR reactions were carried out in a 15-μL reaction volume containing 1 μL of cDNA, 2mM MgCl2, 0.4mM dNTP Mix, 0.4 μM of each primer (Table 1) and 1 U Taq DNA Polymerase (Sinapse). PCR conditions were as follows: 1 cycle of denaturation at 95 °C for 5 min followed by 30 cycles of 95 °C for 30 s, 57 °C for 30 s and 72 °C for 30 s, and a final extension step of 72 °C for 10 min. PCR products were analyzed by electrophoresis (1.5% agarose gel) and stained by ethidium bromide. The expression of the LvActin gene (PCR conditions: 40 cycles of 95 °C for 45 s, 50 °C for 45 s and 72 °C for 1 min) was used as endogenous control to normalize the RT-PCR data for comparison.

Immunodetection of stylicins in shrimp tissues

Whole juvenile shrimp (n=3) were fixed in Davidson's fixative solution (22% formalin, 31.5% ethanol and 11.5% glacial acetic acid) for 24 h at room temperature, embedded in paraffin and cut into 5 μm thick sections. Histological sections of shrimp tissues were deparaffinized and hydrated through xylene-ethanol-water series, and washed in Tris-buffered saline (TBS: 50mM Tris-HCl, 200mM NaCl, pH 7.2). Then, sections were permeabilized (1×TBS, 0.1% Triton X-100) for 30 min and blocked in TBS-T solution (1×TBS, 1% BSA, 0.05% Tween 20) for 2 h followed by 16 h incubation at 4 °C with mouse antirLsty-Stylicin1 polyclonal antibodies (2.3 μg/mL) [START_REF] Rolland | Stylicins, a new family of antimicrobial peptides from the Pacific blue shrimp Litopenaeus stylirostris[END_REF]. After three washes in 1 × TBS+0.05% Tween 20 buffer, sections were incubated for 3 h at room temperature with alkaline phosphatase-labeled rabbit antimouse IgG (1:1000) (Thermo Scientific), followed by a 1 h incubation at room temperature in the dark in a solution of 100 mM Tris-HCl, 100 mM NaCl, 50 mM MgCl2 (pH 9.3) containing 0.175 mg/mL 5-bromo-4-chloro-3-indolyl phosphate (BCIP) (Sigma), 0.375 mg/mL nitro blue tetrazolium (NBT) (Sigma) and 0.24 mg/mL levamisole (Sigma). Negative controls consisted in replacing anti-rLsty-Stylicin1 antibodies with pre-immune mouse serum or TBS-T solution.

Whole-mount immunofluorescence assays were conducted in midgut samples from juvenile shrimp (n=3). Midguts were harvested by dissection, washed in ice-cold Tris-saline solution and immediately fixed in 4% paraformaldehyde. Just after removal of the intestinal content, midguts were longitudinally opened, washed in phosphatebuffered saline (PBS: 137mM NaCl, 2.7mM KCl, 10mM Na2HPO4, 2mM KH2PO4, pH 7.2) and blocked in PBS-T solution (1×PBS, 1% BSA, 0.1% Triton X-100) for 3 h followed by 16 h incubation at 4 °C with mouse anti-rLsty-Stylicin1 polyclonal antibodies (2.3 μg/mL) [START_REF] Rolland | Stylicins, a new family of antimicrobial peptides from the Pacific blue shrimp Litopenaeus stylirostris[END_REF]. Midguts were then washed 3 times with PBS-T solution for 20 min and Incubated for 4 h at room temperature with 4′,6-diamidino-2-phenylindole (DAPI) and FITC-conjugated anti-mouse secondary antibodies diluted at 1:500 (Thermo Scientific). Negative controls consisted in replacing the primary antibodies with pre-immune mouse serum or PBS-T solution. The experiments were repeated three times, and representative images were taken by confocal microscope (Leica DMI6000 B Microscope).

Fluorescence-based reverse transcription real-time quantitative PCR (RT-qPCR)

RT-qPCR reactions were performed in a final volume of 15 μL containing 0.2 μM of each primer (Table 1), 7.5 μL of Maxima SYBR Green/ ROX qPCR Master Mix (Thermo Scientific) and 1 μL of cDNA. The RTqPCR program was 95 °C for 10 min, followed by 40 cycles of 95 °C for 15 s and 60 °C for 60 s. Melt curve analysis (60-95 °C at a temperature transition rate of 0.05 °C/s) was performed to evaluate primer specificity. Primer pair efficiencies (E) were calculated from 2-fold dilution series of pooled cDNA for each primer pair. Primer pair efficiencies were calculated from the given slopes in the StepOne software v2.3, according to the equation: E=10(-1/slope).

The eukaryotic translation elongation factor 1-alpha (LvEF1α) and the ribosomal proteins LvL40, LvRpS3A and LvRpS6 (Table 1) were used as reference genes of expression data in circulating hemocytes. In midgut, LvEF1α, LvL40 and LvActin were used as reference genes for data normalization [START_REF] Silveira | An immune-related gene expression atlas of the shrimp digestive system in response to two major pathogens brings insights into the invoLvement of hemocytes in gut immunity[END_REF]. The relative expression levels of Lvan-Stylicins were calibrated with the expression profile of circulating hemocytes or midguts from naïve (unchallenged) shrimp, according to the 2-ΔΔCq method [START_REF] Livak | Analysis of relative gene expression data using realtime quantitative PCR and the 2(-Delta Delta C(T)) method[END_REF]. Differences in gene expression were considered statistically significant at P < 0.05 (cutoff of 1.5fold change in expression levels) using one-way ANOVA and Tukey's multiple comparison test.

The transcript abundance of Lvan-Stylicins in three midgut portions (anterior, middle and posterior) from Vibrio-challenged and naïve (unchallenged) shrimp was quantified by RT-qPCR and normalized with the gene expression of LvEF1α, LvRpS3A, LvRpS6 and LvActin (Table 1). The relative expression levels were calibrated with the gene expression of each midgut portion from unchallenged shrimp and differences were considered statistically significant at P < 0.05 (cutoff of 1.5-fold change in expression levels) using Student's t-test.

Gene copy number estimation and basal mRNA levels

The number of Lvan-STY gene copies in L. vannamei genome was estimated in five individual shrimp by absolute quantification through the qPCR technique. qPCR reactions were performed using specific primers (Table 1) and 80 ng of gDNA as template. The absolute quantification of the target genes (Lvan-STY1 and Lvan-STY2) was calculated using a standard curve derived from 10-fold dilution series of plasmids containing the DNA target sequences (107 to 103 plasmids/μL; R2=0.998 and 0.999 for Lvan-STY1 and Lvan-STY2, respectively). The transcript abundance of Lvan-Stylicin1 and Lvan-Stylicin2 in the circulating hemocytes of the same five shrimp individuals was assessed by absolute quantification using 1 ng of reverse-transcribed total RNA. The absolute mRNA quantification was performed using the same method described for the gene copy number estimation.

The comparison of the transcript abundance of Lvan-Stylicins between hemocytes and midgut was assessed by RT-qPCR (3 pools of 4 animals per tissue) and normalized with the gene expression of LvEF1α, LvRpS3A and LvRpS6 (Table 1), according to the 2-ΔΔCq method [START_REF] Livak | Analysis of relative gene expression data using realtime quantitative PCR and the 2(-Delta Delta C(T)) method[END_REF]. Differences in gene expression between tissues were considered statistically significant at P < 0.05 (cutoff of 1.5-fold change in expression levels) using Student's t-test.

Quantitative gene expression analysis during shrimp development

Three biological replicates of different development stages of L. vannamei were collected: fertilized eggs at 0-4 h (EI) and at 7-11 h postspawning (EII), nauplius I and V (NI and NV), protozoea I and III (ZI and ZIII), mysis I and III (MI and MIII) and postlarvae aged of 2, 9 and 17 days (PL2, PL9 and PL17). Transcript levels of L. vannamei stylicins during shrimp development were quantified by RT-qPCR and normalized with the gene expression of LvRpS6 and LvActin (Table 1), as previously described [START_REF] Quispe | Transcriptional profiling of immune-related genes in Pacific white shrimp (Litopenaeus vannamei) during ontogenesis[END_REF]. Hemocyte samples from juveniles (3 pools of 5 animals) were used as control for calibrating gene expression data. Genes were considered as "not expressed" in a specific development stage when PCR amplification yielded no product (no dissociation curves) whereas RT-qPCR reactions showing Cq values higher than the limit of quantification (but that generated an expected dissociation curve profile) were considered as "unquantifiable". Statistical significance was considered at P < 0.05 by oneway ANOVA followed by Tukey's multiple comparison test.

Results

L. vannamei stylicins comprise a diverse AAMP family

The transcriptomic analysis of the L. vannamei midgut (unpublished) has revealed the presence of two nucleotide sequences (contigs: DN31608_c0_g1_i1 and DN31608_c0_g1_i2) homologous to the stylicin antimicrobial peptides from the blue shrimp L. stylirostris (Lsty-Stylicin1 and Lsty-Stylicin2). Both sequences correspond to full-length transcripts that encode for precursors composed of a signal peptide followed by a mature peptide containing 13 cysteine residues. The full-length stylicin sequences from L. vannamei were cloned by PCR amplification from gDNA samples and re-sequenced for confirmation. The nucleotide sequence corresponding to the contig DN31608_c0_g1_i2 shared 86% identity with Lsty-Stylicin1 (GenBank: EU177435) whereas the contig DN31608_c0_g1_i1 shared 93% identity with Lsty-Stylicin2 (GenBank: EU177437). Thus, those sequences from L. vannamei were designated as Lvan-Stylicin1 (GenBank: MH108957) and Lvan-Stylicin2 (GenBank: MH108958), respectively (Fig. 1A).

The deduced amino acid sequences of Lvan-Stylicins start with a predicted 22-residue signal peptide followed by an anionic mature peptide of 82 amino acid residues (Fig. 1A; Table S1). Although Lvan-Stylicins contain no obvious protein domains, both mature peptides are characterized by the presence of a proline-rich N-terminal region and a C-terminal region holding the 13 conserved cysteine residues found in L. stylirostris stylicins and in the Mjap-Stylicin from M. japonicus. Lvan-Stylicins were quite similar to each other and to stylicins from L. stylirostris (> 80% identity), but less than 65% amino acid identity was observed between Lvan-Stylicins and Mjap-Stylicin.

In silico analysis of L. vannamei non-annotated databases publicly available on GenBank led to the identification of different isoforms for each L. vannamei stylicin: five for Lvan-Stylicin1 (GenBank: FE179060, FE156583, FE124653, FE125173 and GETD01027084) and three for Lvan-Stylicin2 (GenBank: GETD01027083, GETZ01053100 and GFRP01011277). Those sequences were identified in EST and TSA libraries from multiple shrimp tissues, such as hemocytes, lymphoid organ, nerve cord and hepatopancreas. Additionally, we included in this analysis the Vibrio penaeicida-induced cysteine and proline-rich peptides (LvVICP1 and LvVICP2), stylicin homologues identified in L. vannamei by Wang and colleagues [START_REF] Wang | Analysis of expression, cellular localization, and function of three inhibitors of apoptosis (IAPs) from Litopenaeus vannamei during WSSV infection and in regulation of antimicrobial peptide genes (AMPs)[END_REF]. While Lvan-Stylicin2 sequences differed only by synonymous and nonsynonymous substitutions, Lvan-Stylicin1 sequences showed two distinct lengths: 82 or 83 amino acid residues (Fig. 1B; Table S1). Besides from the synonymous and nonsynonymous substitutions, Lvan-Stylicin1 peptides can be distinguished from each other by the presence/absence of a tyrosine (Y) or histidine (H) residue at the position 51 of the precursor peptide (Fig. 1B). The 82-residue Lvan-Stylicin1 peptides showed a molecular weight of 8.83 kDa and a calculated pI of 4.98 while the 83-residue Lvan-Stylicin1 peptides showed a molecular weight ranging from 8.97 to 9.02 kDa and a calculated pI of 4.98 or 5.19 (Table S1). Comparatively, the mature Lvan-Stylicin2 peptides had a molecular weight of about 9 kDa and a calculated pI ranging from 4.47 to 4.69 (Table S1).

L. vannamei stylicins are encoded by distinct genomic loci

The genomic organization of L. vannamei stylicin genes (Lvan-STY1 and Lvan-STY2) was investigated in two individual shrimp using PCRbased and cloning strategies. Sanger DNA sequencing results showed that Lvan-Stylicins are encoded by distinct genomic sequences, but share a similar structural gene organization (Fig. 2A). Lvan-STY1 and Lvan-STY2 genes are composed by two exons interrupted by a single intron with a length of 217 bp and 205 bp, respectively (Fig. 2A). Both genomic DNA sequences followed the canonical GT/AG splicing recognition rule at the exon/intron boundaries (Fig. S1). In both genes, the first exon covers the 5′-untranslated region (UTR), the signal peptide and the first three residues of the mature peptide while the second exon encodes the remainder of the mature peptide sequence and the 3′-UTR.

Interestingly, the Lvan-STY1 genes from the two sampled shrimp were identical to each other while the nucleotide sequence of their Lvan-STY2 genes differed in ten nucleotides: one in the first exon, five in the intron and four in the second exon (Fig. S1). The five nucleotide substitutions found in the coding sequence resulted in the change of two amino acid residues. The obtained genomic sequences were deposited in GenBank under the accession numbers MH108959 to MH108962.

Then, we asked whether Lvan-STY1 and Lvan-STY2 genes have the same number of copies in L. vannamei genome. The relative gene copy number ratio of Lvan-STY2/Lvan-STY1 was estimated in five individual animals by quantitative PCR. The number of Lvan-STY2 gene copies was 2.37 ± 0.15-fold higher than the number of copies of the Lvan-STY1 gene (Fig. 2B). Finally, a phylogenetic analysis showed that the Lvan-STY genes were placed in the same clade with the stylicin genes from L. stylirostris (Lsty-STY1: EU177436 and Lsty-STY2: EU177437) (Fig. 2C). In this clade, Lvan-STY1 and Lsty-STY1 clustered together in a single group distinct to the STY2 genes (Fig. 2C). Stylicin genes from P. monodon (Pmon-STY: NIUS012084699) and M. japonicus (Mjap-STY: NIUR011088360) clustered in separate groups.

Stylicins cluster into three distinct phylogenetic groups

In silico mining of publicly accessible databases (EST, TSA and WGS) resulted in the identification of novel members of the Stylicin family in different penaeid species: Fenneropenaeus penicillatus (Fpen-Stylicin: GFRT01005742), M. japonicus (Mjap-Stylicin: NIUR011088360) and P. monodon (Pmon-STYlicin: JZ892895, DW678047, DW678039, DT366712, DW042940, GW996588, GEEP01015864, GEME01013089 and NIUS012084699). From our In silico analysis, stylicin sequences were only identified in penaeid shrimp species. All obtained sequences hold the 13 conserved cysteine residues at the C-terminal region (Fig. 3A).

Phylogenetic analysis revealed that the Stylicin family is a monophyletic group that evoLved from a common ancestor gene. Phylogenetic trees constructed with nucleotide and predicted amino acid sequences shared similar topological structures. As shown in Fig. 3B, stylicins clustered in three main clades. The first clade included only stylicin sequences from penaeid species from the genus Litopenaeus. In this clade, Stylicin1 and Stylicin2 were split into two distinct groups (Fig. 3B). Finally, while Pmon-STYlicins and Fpen-Stylicin clustered in a second phylogenetic group, Mjap-Stylicins formed a separate clade from all stylicin sequences (Fig. 3B).

Stylicins are constitutively produced by hemocytes and intestinal cells

The gene expression distribution of L. vannamei stylicins was first evaluated in eight different tissues of shrimp stimulated or not by the injection of heat-killed V. harveyi. Results from the semiquantitative RTPCR analysis evidenced the presence of Lvan-Stylicin1 transcripts in circulating hemocytes, foregut, midgut and gills (Fig. 4). In Vibrio-stimulated animals, Lvan-Stylicin1 gene expression was also detected in hindgut and nerve cord. Comparatively, transcripts of Lvan-Stylicin2 were detected in circulating hemocytes, foregut, midgut, hindgut, gills and nerve cord of both stimulated and non-stimulated animals (Fig. 4). For both Lvan-Stylicins, no signals were observed in muscle and hepatopancreas (Fig. 4).

To characterize the peptide localization of Lvan-Stylicins, immunohistochemistry analysis was subsequently performed on sections of different shrimp tissues using polyclonal antibodies raised against the rLsty-Stylicin1 from L. stylirostris [START_REF] Rolland | Stylicins, a new family of antimicrobial peptides from the Pacific blue shrimp Litopenaeus stylirostris[END_REF]. Due to the high degree of sequence conservation, anti-Lsty-Stylicin1 antibodies probably recognized both Lvan-Stylicins. Stylicin immunoreactivity was found in individual cells heterogeneously distributed across the shrimp tissues. In the cephalothorax, positive immunoreactivity was especially pronounced in some cells present in the connective tissue of the anterior midgut caecum (Fig. 5A andB). Based on morphological features, those stylicin-positive cells are likely tissue-infiltrating hemocytes.

Besides those cells, stylicin immunoreactivity was also observed in the apical region of the columnar epithelial cells lining the midgut and its anterior caecum (Fig. 5A andB). Whole-mount immunofluorescence assays were further performed to confirm the presence of stylicin peptides in those epithelial cells. Confocal images clearly evidenced the presence of stylicin-containing granules located at the apical region of the midgut columnar epithelial cells (Fig. 5C-E). No signals were observed in other cell types of the shrimp body, thus the results of the semiquantitative RT-PCR analysis are probably the consequence of the infiltration of stylicin-expressing hemocytes in shrimp tissues. Altogether, results from both immunohistochemistry and whole-mount immunofluorescence assays showed that Lvan-Stylicins are constitutively produced by the hemocytes and by the columnar epithelial cells of the midgut.

Lvan-STY1 and Lvan-STY2 genes are simultaneously transcribed in a single shrimp at different basal levels

Since Lvan-Stylicins are produced by both hemocytes and midgut cells, we focused on determining the main site of stylicin expression in penaeid shrimp. Results showed that the expression of Lvan-Stylicin1 and Lvan-Stylicin2 was, respectively, 52.5-fold and 62.2-fold higher in circulating hemocytes than in the midgut (Fig. 6A). Next, the basal mRNA expression levels of the Lvan-STY1 and Lvan-STY2 genes were analyzed in the circulating hemocytes of five individual shrimp by absolute quantification. Interestingly, both Lvan-STY genes showed to be constitutively and simultaneously transcribed in an individual shrimp, but at different transcriptional levels. The basal expression of Lvan-STY2 was 3.69-fold higher than Lvan-STY1 (Fig. 6B). Besides, the basal mRNA levels of each gene showed to be also variable among the individuals (Fig. 6B). For the Lvan-STY1 gene, differences in gene expression reached up to 3.78-fold whereas variations up to 2.36-fold were found in Lvan-STY2 gene expression.

The gene expression of L. vannamei stylicins is differentially modulated in response to infections

The gene expression profile of Lvan-Stylicins was further quantified by fluorescence-based quantitative PCR (RT-qPCR) in shrimp hemocytes and midgut at 48 h after infections with two unrelated pathogens, the Gram-negative V. harveyi and the WSSV. This time point was chosen on the basis of previous studies from our group [START_REF] Goncalves | Molecular signatures at imminent death: hemocyte gene expression profiling of shrimp succumbing to viral and fungal infections[END_REF][START_REF] Silveira | An immune-related gene expression atlas of the shrimp digestive system in response to two major pathogens brings insights into the invoLvement of hemocytes in gut immunity[END_REF]. Interestingly, while the expression of Lvan-Stylicin2 was induced in circulating hemocytes in response to the Vibrio infection (2.4-fold change), the expression of Lvan-Stylicin1 was not modulated (Fig. 7). No increase in Lvan-Stylicin2 expression was observed following the injection of sterile seawater (aseptic injury control). Besides, both genes were not modulated after the viral infection. In contrast, the expression of both genes was up-regulated in circulating hemocytes (2.5-fold change for Lvan-Stylicin1 and 2.9-fold change for Lvan-Stylicin2) after the injection of a muscle tissue homogenate prepared from WSSV-free shrimp (Fig. 7). In the midgut, the expression of Lvan-Stylicins was not regulated by the bacterial or by the viral infection (Fig. 7).

The lack of gene expression response in the midgut encouraged us to conduct an alternative experimental infection method, mimicking a more natural route of bacterial infection (per os challenge). Variations in gene expression of Lvan-Stylicins were assessed by RT-qPCR in three midgut portions (anterior, middle and posterior). The expression of both Lvan-Stylicins was up-regulated (2.89-fold change for Lvan-Stylicin1 and 4.31-fold change for Lvan-Stylicin2) only in posterior portion of the midgut (Fig. 8). By contrast, no increase in Lvan-Stylicin gene expression was observed in the two first midgut portions (anterior and middle) in response to the oral Vibrio challenge (Fig. 8).

L. vannamei stylicins show a different pattern of gene expression during shrimp development

We finally investigated the presence and the levels of stylicin transcripts in tweLve developmental stages of L. vannamei, from fertilized eggs to larval (nauplius, protozoea and mysis) and postlarval stages, and also in circulating hemocytes from juveniles. Transcript levels of Lvan-Stylicin1 were detected at very low levels in late protozoea stages (ZIII), but its expression was only quantified from mysis III (MIII). Then, Lvan-Stylicin1 expression increased gradually in the following developmental stages (Fig. 9). On the other hand, Lvan-Stylicin2 transcripts were found to be present early in shrimp development (fertilized eggs at 7-11 h post-spawning). However, Lvan-Stylicin2 expression was only quantified from the protozoea III (ZIII) stage (Fig. 9). For both genes, the highest mRNA abundance was quantified in hemocytes from juvenile shrimp (Fig. 9).

Discussion

We showed here that L. vannamei stylicins comprise a diverse family of anionic antimicrobial peptides (AAMPs) whose genes are differentially regulated in hemocytes and midgut cells in response to infections. From the four gene-encoded AMPs described in penaeid shrimp, only stylicins have not been fully characterized in the Pacific white shrimp (L. vannamei) and this is the first study exploring the diversity of the stylicin family in terms of sequence and gene expression distribution and regulation. By taking advantage of RNA-Seq technology, we have identified two stylicin homologues (Lvan-Stylicin1 and Lvan-Stylicin2) in midgut transcriptomes of L. vannamei showing high similarities to stylicins from the blue shrimp L. stylirostris [START_REF] Rolland | Stylicins, a new family of antimicrobial peptides from the Pacific blue shrimp Litopenaeus stylirostris[END_REF]. From our In silico analysis, stylicins form a diverse gene family in shrimp species of the genus Litopenaeus, but not in other taxa of the family Penaeidae, such as Fenneropenaeus, Penaeus and Marsupenaeus. Essentially, our results revealed that all known Litopenaeus gene-encoded AMPs are present as diverse multigene families composed of different members. For instance, while Litopenaeus penaeidins are composed of three members (Litvan PEN1/2, Litvan PEN3 and Litvan PEN4), at least four and seven members were identified in the crustin (Crustin Lv, Crustin-like Lv, LvSWD and LvSPLI) and in the ALF families (Litvan ALF-A to -G), respectively [START_REF] Destoumieux-Garzón | Antimicrobial peptides in marine invertebrate health and disease[END_REF][START_REF] Rosa | Antimicrobial peptides in crustaceans[END_REF][START_REF] Matos | Massive gene expansion and sequence diversification is associated with diverse tissue distribution, regulation and antimicrobial properties of anti-lipopolysaccharide factors in shrimp[END_REF]. However, unlike other gene-encoded AMPs from marine invertebrates, stylicins are exclusively composed of anionic peptides. Interestingly, the spectrum of activity of these anionic antimicrobials is restricted to filamentous fungi [START_REF] Rolland | Stylicins, a new family of antimicrobial peptides from the Pacific blue shrimp Litopenaeus stylirostris[END_REF], even if their gene expression has shown to be associated to shrimp survival to pathogenic Vibrio infections [START_REF] De Lorgeril | Identification of genes that are differentially expressed in hemocytes of the Pacific blue shrimp (Litopenaeus stylirostris ) surviving an infection with Vibrio penaeicida[END_REF]. Likewise, anionic peptides derived from the Cterminus of the shrimp respiratory protein hemocyanin are also specific against fungi [START_REF] Destoumieux-Garzón | Crustacean immunity[END_REF]. These hemocyaninderived peptides, named PvHCt, are able to bind and permeabilize fungal membranes [START_REF] Petit | A hemocyaninderived antimicrobial peptide from the penaeid shrimp adopts an alpha-helical structure that specifically permeabilizes fungal membranes[END_REF]. To date, the mechanism of action of stylicins and of other AAMPs is completely unknown. Actually, in comparison to CAMPs, few scientific groups have attempted to study the subject of AAMPs [START_REF] Harris | Anionic antimicrobial peptides from eukaryotic organisms[END_REF] and more functional studies are needed to achieve a more in-depth understanding of these unconventional AMPs.

We showed here that Lvan-Stylicin1 and Lvan-Stylicin2 are encoded by different genomic loci. Nonetheless, the presence of distinct stylicin genes was only observed in Litopenaeus species. Thus, it is likely that STY1 and STY2 are paralogous genes that arose from a single gene duplication event before the speciation of the genus Litopenaeus. Besides, although both Lvan-STY genes share a similar structural gene organization, their copy numbers in L. vannamei genome showed to be variable, suggesting that they have followed independent duplication events after STY1-STY2 divergence. Apart from the two main stylicins members (Lvan-Stylicin1 and Lvan-Stylicin2), we have identified Lvan-Stylicin1 sequences with distinct lengths (82 and 83 amino acids in length). No evidences for the presence of alternative splicing were found in the Lvan-STY1 gene, and thus it is most likely that different alleles may occur for L. vannamei stylicin genes. Taken together, our results suggest that L. vannamei stylicins belong to a diverse multigenic and multiallelic family of AAMPs. Moreover, Lvan-STY1 and Lvan-STY2 genes showed to be simultaneously transcribed in a single shrimp. The different penaeidin and ALF members are also simultaneously expressed in an individual shrimp [START_REF] Matos | Massive gene expansion and sequence diversification is associated with diverse tissue distribution, regulation and antimicrobial properties of anti-lipopolysaccharide factors in shrimp[END_REF][START_REF] Cuthbertson | Diversity of the penaeidin antimicrobial peptides in two shrimp species[END_REF], and it would be of great interest to colocalize the four gene-encoded AMP families in shrimp hemocytes and examine whether they function synergistically to enhance their antimicrobial spectrum of activity. Synergic activities have already been observed for AMPs from both vertebrate [START_REF] Yan | Synergistic interactions between mammalian antimicrobial defense peptides[END_REF] and invertebrate species [START_REF] Schmitt | Expression, tissue localization and synergy of antimicrobial peptides and proteins in the immune response of the oyster Crassostrea gigas[END_REF][START_REF] Marxer | Insect antimicrobial peptides act synergistically to inhibit a trypanosome parasite[END_REF], but unfortunately, no studies of this type have been addressed in crustaceans.

One of the important findings of this study is that the expression Lvan-Stylicins is not limited to the immune cells (hemocytes) as observed for the other three shrimp gene-encoded AMPs [START_REF] Destoumieux-Garzón | Antimicrobial peptides in marine invertebrate health and disease[END_REF]. Both immunostaining and gene expression analysis revealed that Lvan-Stylicins are constitutively produced by the midgut columnar epithelial cells and that their expression is induced in response to Vibrio infections. Notably, the expression of Lvan-Stylicin was pronounced in the anterior midgut caecum, an intestinal region primarily invoLved in the production, secretion and activation of digestive enzymes [START_REF] Mcgaw | A review of gastric processing in decapod crustaceans[END_REF]. Thus, it is plausive to hypothesize that Lvan-Stylicins are secreted into the midgut lumen. This finding brings new insights into the role of crustacean AMPs in the control of the gut microbiota and in shrimp intestinal defenses. Indeed, the shrimp midgut lacks the cuticular lining found in the other portions of the intestine (foregut and hindgut), representing a potential route of entry for many pathogens into the hemocel. In a previous study, we showed that gut is an important source for the expression of immune-related genes in penaeid shrimp [START_REF] Silveira | An immune-related gene expression atlas of the shrimp digestive system in response to two major pathogens brings insights into the invoLvement of hemocytes in gut immunity[END_REF], thus the presence of stylicins in the midgut cells suggests the participation of these antimicrobial effectors in the first intestinal line of defense. Interestingly, the expression of some ALF members showed to be invoLved in the maintenance of the microbiota residing in the shrimp hemolymph [START_REF] Wang | A shrimp C-type lectin inhibits proliferation of the hemolymph microbiota by maintaining the expression of antimicrobial peptides[END_REF]. Taken all together, our results highlight the importance of stylicins in both hemolymph-based and gut-based immunities.

Another relevant conclusion that can be drawn from this study is that Lvan-STY genes are differentially regulated after infections. Such diversity in gene expression regulation was previously reported for the different members of the L. vannamei ALF family in response to fungal infections [START_REF] Rosa | Functional divergence in shrimp anti-lipopolysaccharide factors (ALFs): from recognition of cell wall components to antimicrobial activity[END_REF]. Interestingly, while some AMPs such as Lvan-Stylicin1, penaeidins [START_REF] Destoumieux | Penaeidins, antimicrobial peptides with chitin-binding activity, are produced and stored in shrimp granulocytes and released after microbial challenge[END_REF], Litvan ALF-A [START_REF] Rosa | Functional divergence in shrimp anti-lipopolysaccharide factors (ALFs): from recognition of cell wall components to antimicrobial activity[END_REF] and Type II crustins [START_REF] Barreto | Specific molecular signatures for type II crustins in penaeid shrimp uncovered by the identification of crustin-like antimicrobial peptides in Litopenaeus vannamei[END_REF] are notregulated in response to infections, the expression of other L. vannamei AMPs (Lvan-Stylicin2 and other ALF members) is induced in immune cells to improve host antimicrobial responses against pathogens. As elegantly shown by Wang et al. [START_REF] Wang | Analysis of expression, cellular localization, and function of three inhibitors of apoptosis (IAPs) from Litopenaeus vannamei during WSSV infection and in regulation of antimicrobial peptide genes (AMPs)[END_REF], the gene expression of most shrimp AMPs, including penaeidins, Type II crustins and stylicins (LvVICPs), is controlled by the IMD/NF-κB pathway, an evolutionarily conserved signaling cascade invoLved in the regulation of the antimicrobial responses of arthropods. On the other hand, the expression of almost all shrimp AMPs, including Lvan-Stylicins, can be drastically affected in circulating hemocytes during lethal infections [START_REF] Goncalves | Molecular signatures at imminent death: hemocyte gene expression profiling of shrimp succumbing to viral and fungal infections[END_REF]. Unlike the Mjap-Stylicin from the kuruma prawn M. japonicus [START_REF] Liu | Identification and expression analysis of a novel stylicin antimicrobial peptide from Kuruma shrimp (Marsupenaeus japonicus)[END_REF], Lvan-Stylicins were not modulated neither in circulating hemocytes nor in midgut by the WSSV. This result could be probably due to the time course response of Lvan-Stylicins in those shrimp tissues. On the other hand, both Lvan-STY genes have been shown to be responsive to a tissue homogenate prepared from shrimp muscle (injury control for the WSSV infection). This finding strongly suggests that Lvan-Stylicins are induced in response to danger/damage-associated molecular patterns (DAMPs) and that these AAMPs could be invoLved in early inflammation and in wound healing processes as proposed for penaeidins [START_REF] Destoumieux | Penaeidins, antimicrobial peptides with chitin-binding activity, are produced and stored in shrimp granulocytes and released after microbial challenge[END_REF].

L. vannamei stylicin genes showed to be differentially regulated not only in response to infections, but also during shrimp development. Interestingly, Lvan-Stylicin2 and the M. japonicus stylicin showed a very similar pattern of gene expression during shrimp development [START_REF] Liu | Identification and expression analysis of a novel stylicin antimicrobial peptide from Kuruma shrimp (Marsupenaeus japonicus)[END_REF]. Comparatively, Lvan-Stylicin1 expression was quite similar to that observed for other L. vannamei AMPs, such as Litvan PEN1/2, Litvan PEN4 and Litvan ALF-D [START_REF] Quispe | Transcriptional profiling of immune-related genes in Pacific white shrimp (Litopenaeus vannamei) during ontogenesis[END_REF]. Apart from that, we cannot discard the possibility of these changes in gene expression could be the result of the differences in the mRNA basal levels observed between the Lvan-STY genes. Indeed, the detection of stylicin transcripts, as well as of other gene-encoded AMPs, highlights the importance of these antimicrobial effectors during shrimp development.

Conclusions

In conclusion, we showed that the stylicin family from shrimp species of the genus Litopenaeus is composed of two members encoded by distinct genomic loci that exhibit different patterns of gene expression distribution and regulation. According to the best of our knowledge, this is the first evidence for the expression of a shrimp gene-encoded AMP in other tissues than the hemocytes. Even though the expression of stylicins has been shown to be a marker of shrimp survival to pathogenic Vibrio infections, the of these effectors in shrimp immune defense is still largely unknown. The application of RNAibased methods could make significant contributions to understanding the significance of the molecular and transcriptional diversity of Litopenaeus stylicins in host-microbe interactions.

FIGURE 9

Gene expression profile of Litopenaeus vannamei stylicins during shrimp development. EI: fertilized eggs at 0-4 h post-spawning; EII: fertilized eggs at 7-11 h post-spawning; NI: nauplius I; NV: nauplius V; ZI: protozoea I; ZIII: protozoea III; MI: mysis I; MIII: mysis III; PL2: postlarva 2; PL9: postlarva 9; PL17: postlarva 17. The red dotted line indicates the basal expression level in hemocytes from juvenile shrimp while the solid blue underline highlights the stages at which the gene expression was detected (valid dissociation curve profile) but not quantified (Cq values higher than the limit of quantification). Results are presented as mean ± standard deviation (three biological replicates). Different letters indicate significant differences among the development stages while asterisks (*) shows significant differences between each developmental stage and hemocytes from juveniles (oneway ANOVA/Tukey, P < 0.05). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

FIGURE 6 (

 6 FIGURE 6 (A) Quantitative comparison of the relative abundance of Lvan-Stylicin1 (white bars) and Lvan-Stylicin2 (black bars) transcripts in circulating hemocytes (HE) and midgut (MG). Results are presented as mean ± standard deviation of relative expressions (three biological replicates) and statistical differences are indicated by (*) (Student's t-test, P < 0.05). (B) Transcript abundance of Lvan-Stylicins in five individual shrimp. The absolute quantification was assessed by qPCR using a standard curve derived from 10-fold dilution series of plasmids containing each target gene. Results are presented as the ratio of the abundance of Lvan-Stylicin2 transcripts per ng of total RNA to that of Lvan-Stylicin1. The β-actin gene (LvActin) was used as endogenous expression control for each individual shrimp (lower panel).

FIGURE 7

 7 FIGURE 7Relative gene expression profile of Lvan-Stylicin1 (white bars) and Lvan-Stylicin2 (black bars) in circulating hemocytes and midgut of shrimp at 48 h after experimental infections with the Gramnegative Vibrio harveyi ATCC 14126 (6×107 CFU/animal) or the White spot syndrome virus (WSSV: 3×102 viral particles/ animal). Results are presented as mean ± standard deviation of relative expressions (three biological replicates) and statistical differences are indicated by asterisks (*) (oneway ANOVA/Tukey, P < 0.05). W-free: tissue homogenate inoculum prepared from WSSV-free shrimp.

FIGURE 8

 8 FIGURE 8 Relative gene expression profile of Lvan-Stylicin1 and Lvan-Stylicin2 in three portions of shrimp midgut (anterior, middle and posterior) at 21 h after the oral administration of 7.5×105 CFU/animal of live V. harveyi ATCC 14126. Results are presented as mean ± standard deviation of relative expressions (three biological replicates) and statistical differences are indicated by asterisks (*) (Student's t-test, P < 0.05). The figure (not-to-scale) shown at the top of the figure indicates the anatomic location of the three midgut portions from naïve (white bars) and Vibrio-challenged (black bars) shrimp.

Acknowledgments

We are grateful to the Laboratory of Marine Shrimp (Federal University of Santa Catarina -UFSC) for providing the shrimp used in this study and to Fabio S Ribeiro for his technical assistance in the whole-mount immunofluorescence experiments. The authors are also thankful to the Electronic Microscopy Central Laboratory (LCME-UFSC) and the Multiuser Laboratory of Biology Studies (LAMEB-UFSC). This work was supported by the Brazilian funding agencies CNPq (MEC/MCTI/CAPES/CNPq/FAPs PVE 401191/2014-1 and MCTI/CNPq Universal 406530/2016-5) and CAPES (CIMAR 1974/2014). P Schmitt was funded by FONDECYT grant 11150009. ND Farias and GM Matos were supported by scholarships provided by FAPESC while N Argenta received a scholarship provided by CNPq. M Falchetti and C Barreto were supported by fellowships provided by CAPES.

TABLE 1

Nucleotide sequences of primers used in this study.