Emeline Talansier
email: emeline.talansier@univ-grenoble-alpes.fr

Audrey Bacconnier

François Caton

Carine Chastel

Lucy Costa

Deniz Z Gunes

Denis C D Roux
email: denis.roux@univ-grenoble-alpes.fr

Accurate methodology to determine slip velocity, yield stress and the constitutive law for molten chocolate

Keywords: Yield stress fluid, Slip velocity, Solid-liquid transition, Herschel-Bulkley model, Casson model, Molten chocolate

Flow of concentrated soft materials like molten chocolate often has the ambiguity of behaving like a solid or a liquid depending on the external constraints.

Determination of the constitutive law is often unreliable and consequently the yield stress value too. This difficulty of determining yield stress, particularly when using a standards protocol such as the IOCCC for chocolate, is a consequence of the presence of wall slip. Most attempts to assess slip velocity using Yoshimura and Prud'homme's method have resulted in failure. We propose to tackle this problem by using an easy and reliable methodology applied to a commercial chocolate in a melted state with, as equipment, a plate-plate rotational rheometer. This method allows precise determination of the yield stress value, the constitutive law and the slip velocity, as well as also providing an explanation to the apparent Newtonian plateau and thickening behavior of the raw flow curve.

Introduction

Molten chocolate, a concentrated suspension of large and solid particles within a continuous fat system, is known to have a non-Newtonian behavior.

Sedimentation, wall slip or solid to liquid transition occur in such complex systems [START_REF] Barnes | A review of the slip (wall depletion) of polymer solutions, emulsions and particle suspensions in viscometers: its cause, character, and cure[END_REF][START_REF] El Kissi | Glissement et Rhéométrie[END_REF] making their rheological characterization difficult although this is essential for the sizing of industrial processes.

However, proper flow characterization is of critical importance for industrial process design, as well as for determining the best adapted formulation window for a given factory process. In the food industry, this is particularly true for chocolate manufacture, since most industrial processes have been optimized since a long time ago. In this context, the International Office of Cocoa, Chocolate and Confectionary office (IOCCC) and the International Chocolate Awards (ICA) organizations have developed specific protocols for molten chocolate [START_REF] Afoakwa | Comparison of rheological models for determining dark chocolate viscosity[END_REF][START_REF] Aeschlimann | International inter-laboratory trials to determine the factors affecting the measurement of chocolate viscosity[END_REF]. This protocol, currently accepted as an accepted standard, consists of successive increasing and decreasing shear rate ramps set at a series of values between 2 and 50 s -1 (see Table 1a).

Using either a Casson or an Herschel-Bulkley model, the fit of the experimental data gives the parameters of a constitutive equation [START_REF] Steffe | Rheological Methods in Food Process Engineering[END_REF]. Among these parameters, the yield stress below which the material behaves as a solid, is an important parameter difficult to determine. From the IOCCC standard method the following question arises: "Why are only shear rates higher than 2s -1 used to characterize liquid chocolate?" This is quite surprising since it is well-known that lower shear rates, below 2s -1 , are today easily attainable with any commercial rheometer.

The answer is certainly hidden within the complex solid to liquid transition of concentrated particles systems for which a resurgence of interest has been seen in the scientific literature [START_REF] Balmforth | Yielding to Stress: Recent Developments in Viscoplastic Fluid Mechanics[END_REF][START_REF] Bonn | Yield Stress Materials in Soft Condensed Matter[END_REF][START_REF] Bonn | Yield stress fluids slowly yield to analysis[END_REF][START_REF] Damianou | Cessation of viscoplastic Poiseuille flow with wall slip[END_REF][START_REF] Divoux | From stress-induced fluidization processes to Herschel-Bulkley behaviour in simple yield stress fluids[END_REF][START_REF] Grenard | Timescales in creep and yielding of attractive gels[END_REF][START_REF] Herzhaft | Correlation between transient shear experiments and structure evolution of aqueous foams[END_REF][START_REF] Herzhaft | Measurement and modeling of the flow behavior of aqueous foams using a recirculating pipe rheometer[END_REF][START_REF] Piau | Carbopol gels: Elastoviscoplastic and slippery glasses made of individual swollen sponges: Meso-and macroscopic properties, constitutive equations and scaling laws[END_REF][START_REF] Salmon | Towards local rheology of emulsions under Couette flow using Dynamic Light Scattering[END_REF]; "Special Issue on Yield Stress Fluids: a 100 Years after Bingham's Landmark Paper," 2017). However, even if our understanding of the solid to liquid transition is gradually improving, few works have been reported on industrial soft condensed matter. Moreover, all the works which have shown interest in solid to liquid transition face the issue of material slip at the solid wall boundaries.

Wall slip velocities are usually quantified by means of capillary rheometry [START_REF] Mooney | Explicit Formulas for Slip and Fluidity[END_REF] applied to Rabinowitsch-corrected data. Using the same hypothesis as Mooney, Yoshimura and Prud'homme [START_REF] Yoshimura | Wall Slip Corrections for Couette and Parallel Disk Viscometers[END_REF] developed a fast method based on two sets of experiments to assess wall slip using rotational rheometry. Yilmazer and Kalyon [START_REF] Yilmazer | Slip Effects in Capillary and Parallel Disk Torsional Flows of Highly Filled Suspensions[END_REF] generalized this method by applying the Mooney plots procedure to parallel disk geometry and obtained a better value for wall slip velocity.

Nonetheless, the Yilmazer generalization on rotational rheometry has never been assessed on yield stress fluids with the exception of Bertola's work [START_REF] Bertola | Wall slip and yielding in pasty materials[END_REF] on aqueous foam and emulsion.

Nevertheless, we should also notice that in liquid chocolate, yield stress determination using the Vane or Couette geometries along with the IOCCC standard method does not take into account the slippage of the material at the wall boundary [START_REF] Aeschlimann | International inter-laboratory trials to determine the factors affecting the measurement of chocolate viscosity[END_REF][START_REF] Afoakwa | Comparison of rheological models for determining dark chocolate viscosity[END_REF][START_REF] Afoakwa | Effects of particle size distribution and composition on rheological properties of dark chocolate[END_REF][START_REF] Servais | Determination of Chocolate Viscosity[END_REF][START_REF] Sokmen | Influence of some bulk sweeteners on rheological properties of chocolate[END_REF].

In the present paper we present and discuss the applicability of conventional rotational rheometry methods for accurate quantification of wall slip velocity and yield stress determination in a commercial molten chocolate.

Materials and methods

Composition and samples preparation

"Couverture" chocolate containing 64.5% (w/w) cocoa, 27% sugar and 8% fat from the company Barry Callebaut was used. Emulsifier soya lecithin and natural vanilla flavor completed the 0.5% of the remaining components of the mix. This highly fluid chocolate is usually used for molded, coating or tablet use.

Molten chocolate samples were put into a sealed glass bottle which was placed onto a rotating mixer (RM5) inside a Memmert thermostatic chamber (Model UFE 600, 256L). The chamber was maintained at T = 50°C (+/-0.5°C) for a minimum period of two days. This protocol ensured that the chocolate was fully liquid and free of any remaining crystals. This preparation acted as the thermal and sample homogenization step.

Rheometry

Rheological measurements were performed using an imposed torque rheometer MCR501 (Anton Paar, GmbH) with minimum torque and rotational velocity measures of 0.1 µN.m and 10 -6 rpm respectively possible. A stainless steel, smooth plate-plate geometry of 25 mm diameter was employed on the rheometer with different gaps of 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.6, 1.7 and 1.8 mm.

In accordance with the continuum mechanics hypothesis, all gaps were much larger than the largest particles in suspension ( ) in this kind of molten chocolate [START_REF] Aeschlimann | International inter-laboratory trials to determine the factors affecting the measurement of chocolate viscosity[END_REF]. All measurements were performed at 40°C, using the rheometer's Peltier system. To avoid temporal behavior on creep flow at low stresses, but also aging or complex dynamics induced by residual stresses, thixotropy and retarded viscoelasticity [START_REF] Benmouffok-Benbelkacem | Nonlinear viscoelasticity and temporal behavior of typical yield stress fluids: Carbopol, Xanthan and Ketchup[END_REF][START_REF] Caton | Plastic behavior of some yield stress fluids: from creep to long-time yield[END_REF][START_REF] Lidon | Power-law creep and residual stresses in a carbopol gel[END_REF], measurements were conducted by imposing successive increasing and decreasing shear rate ramps as fully described by Table 1b and Fig. 1. Increasing shear rate ramps are known to be different from decreasing ones, as illustrated in Fig. 1 by the two first ramps of 0 and 5min. Indeed, when passing from solid to liquid state, yield stress fluids may have several different behaviors [START_REF] Cheddadi | Steady Couette flows of elastoviscoplastic fluids are nonunique[END_REF][START_REF] Saramito | A new constitutive equation for elastoviscoplastic fluid flows[END_REF] that make the exploitation of increasing ramps impossible. Increasing ramps were only used to fluidize material and the point settings were therefore fewer than those of the decreasing ramps. If compared with the standards presented in Table 1a, the proposed decreasing ramps attained lower shear rates, as low as 0.005s -1 , thus permitting an accurate determination of the constitutive law parameters. Whatever the starting time of the measured decreasing ramp (5min, 1h23 or 2h41min), all apparent flow curves were shown to be perfectly superimposed, providing that a steady state was reached. Thus, each measured point was accepted as valid if the comparison of five consecutive recordings were within a tolerance limit of 2% within a time interval of 50 seconds maximum (cf Table 1b). Apparent shear rate [s -1 ] 0min 5min 1h23min 2h41min

If the tolerance limit was over 2%, the last value measured (at 50s) was selected.

In addition, measurements were made at minimum in triplicate giving a maximum experimental error of 11%. The middle ramp was only selected for further exploitation in the case of perfect replicate superimposition implying that the surrounding ramps have been used as references.

By assuming either sugar and cocoa particles to be spheres of diameter dpart 20µm with a density of ~1500 kg/m 3 dispersed in a continuous fat cocoa butter phase at 40°C ( =900 kg/m 3 , =100Pa.s), a quick estimate of the sedimentation velocity, assuming no inter-particle interactions, is given in the order of 5.10 -8 m/s by :

(1)

The corresponding sedimentation time inside the plate-plate geometry is therefore around 50 hours for a 1mm gap width. Knowing that an increasing-decreasing shear rates ramps cycle takes about 38 minutes, we can deduce that no sedimentation took place in our experiments. This is confirmed by the perfect superimposition of all the curves in Fig. 1.

Finally, the normal force, i.e. perpendicular to the rotation plane, will not be displayed here as the entire measurement campaign values remained within the 15 10mN interval.

Results and discussion

1mm-gap raw flow curves of Apparent shear rate [s -1 ] 1.50 mm 1.00 mm 0.75 mm 0.50 mm instance, gels [START_REF] Clasen | Determining the true slip of a yield stress material with a sliding plate rheometer[END_REF]Meeker et al., 2004aMeeker et al., , 2004b)), mayonnaise [START_REF] Clasen | The flexure-based microgap rheometer (FMR)[END_REF], fiber pulp suspensions [START_REF] Derakhshandeh | Rheology of pulp suspensions using ultrasonic Doppler velocimetry[END_REF], toothpaste [START_REF] Ardakani | Thixotropic flow of toothpaste through extrusion dies[END_REF], shaving foam, emulsions [START_REF] Bertola | Wall slip and yielding in pasty materials[END_REF] or liquid foams [START_REF] Laporte | Rheological properties of food foams produced by SMX static mixers[END_REF]. This type of behavior, for which the measurements are well superimposed at high values, while remaining highly dependent on the characteristic gap length at low values is usually interpreted as the presence of a fluid wall slip [START_REF] Bertola | Wall slip and yielding in pasty materials[END_REF][START_REF] Cloitre | A review on wall slip in high solid dispersions[END_REF][START_REF] El Kissi | Glissement et Rhéométrie[END_REF][START_REF] Hatzikiriakos | Slip mechanisms in complex fluid flows[END_REF]Meeker et al., 2004aMeeker et al., , 2004b;;[START_REF] Poumaere | Unsteady laminar flows of a Carbopol® gel in the presence of wall slip[END_REF][START_REF] Sochi | Slip at Fluid-Solid Interface[END_REF].

Corrections of rheometrical raw data

Shear stress correction

The science of rheology is based on the experimental determination of the material constitutive equation that links shear rate and shear stress, respectively imposed and measured by the rheometer. In the simplest shear flow configuration, the fluid is contained and sheared between two parallel circular plates in relative rotation, at a velocity [rad.s -1 ] around the vertical axis (Oz)

and separated by a distance H (Fig. 3a).

In plate-plate geometry, a homogeneous fluid shows a velocity field reduced to its tangential component which depends solely on the distance of the rotational axis and of the vertical position in the gap z:

(2)

The corresponding apparent local shear rate is expressed as:

(3)

Inside the gap, the shear stress is related to the torque as follows:

(4)

Considering Eq.3, the torque may be rewritten as:

(5)

and at the rim of plate-plate geometry, the differentiation of Eq.5 finally gives [START_REF] Mooney | Explicit Formulas for Slip and Fluidity[END_REF][START_REF] Rabinowitsch | Über die viskosität und elastizität von solen[END_REF][START_REF] Yilmazer | Slip Effects in Capillary and Parallel Disk Torsional Flows of Highly Filled Suspensions[END_REF]:

(6)
In all these equations, the "a" and "R" indices are representative of the words "apparent" and "rim". In the case of Newtonian viscous fluid, the raw shear stress which may also be called apparent shear stress, is expressed as:

At this stage, we can report that a similar correction on a Non-Newtonian aqueous hydroxyethyl cellulose solution directly carried on the viscosity was published by [START_REF] Krieger | Direct Determination of the Flow Curves of Non-Newtonian Fluids. IV. Parallel-Plane Rotational Viscometer[END_REF][START_REF] Krieger | Direct Determination of the Flow Curves of Non-Newtonian Fluids. IV. Parallel-Plane Rotational Viscometer[END_REF]. This alternative yields corrected viscosity values obviously identical to those obtained via Eq.6. 1.8mm 1.5mm 1.0 mm 0.75 mm 0.50 mm stresses, the raw data is below the corresponding solid lines leading to an underestimation of the apparent shear stress without correction.

Shear rate correction

Once the shear stress had been corrected, attention was focused on the correction of the shear rate values. To do this, the sample inside the plate-plate tool illustrated in Fig. 3 was used again as follows. With a motionless lower plate, the upper plate was rotated around the (Oz) axis at a rotation speed . Where there was presence of slippage phenomena at the solid walls, i.e. at z=0 and z=H in Fig. 3b, the apparent shear rate was no longer homogeneous. This wall slippage can be understood macroscopically as the effect of a thin depleted layer at the solid walls sandwiching the bulk material [START_REF] Barnes | A review of the slip (wall depletion) of polymer solutions, emulsions and particle suspensions in viscometers: its cause, character, and cure[END_REF][START_REF] Cloitre | A review on wall slip in high solid dispersions[END_REF]Meeker et al., 2004aMeeker et al., , 2004a)). This depleted layer of thickness h is submitted to a high shear rate near the solid walls whereas on the bulk length the rest of the fluid remains less sheared.

Based on Mooney's pioneer work [START_REF] Mooney | Explicit Formulas for Slip and Fluidity[END_REF], the apparent slip velocity at the solid walls is expressed as:

where and are the shear rates in the bulk and in the depleted layer, respectively. A schematic drawing of the slip between the plates can be seen in Fig. 3b.

Assuming and , with is the total length between the two plates, the apparent shear rate at the rim , can be written as:

(9)

Using the fact that the constitutive equation links the shear stress and the shear rate, Eq.9 can also be written as:

(10)

where is the bulk shear rate at a given wall stress value, , and is the wall slip velocity at the same shear stress. In Eq.10, the bulk shear rate and the slip velocity are expressed at a fixed shear stress . In the rheological trials conducted here, all measurements were handled at fixed shear rates. As previously mentioned ( § 3.1.1), all shear stress data was corrected and interpolated with a "csaps" tool in order to use Eq.10 at a fixed shear stress . This representation of the apparent shear rate vs 1/H has been used for highly filled suspensions in the past [START_REF] Yilmazer | Slip Effects in Capillary and Parallel Disk Torsional Flows of Highly Filled Suspensions[END_REF]. As illustrated by Figs.5a, 5b and 5c, our data has been split into three different plots.

In Fig. 5a for low shear stresses, the apparent shear rates are shown to be strictly proportional to 1/H. According to Eq.10, the best fit of experimental data is reduced to an affine function where only wall slip occurs (Fig. 6, green region).

The apparent shear of the sample is then reduced to a solid-like slippage at the wall, whereas the "fluid" bulk remains unsheared but only elastically distorted. In Fig. 5c for high shear stresses, apparent shear rates are represented by the horizontal straight lines (Fig. 6, blue region). Thus, in Eq.10 the bulk shear rate contribution largely exceeds the terms of slippage implying equality between applied and measured shear rates.

In Fig. 5b for intermediate shear stresses, the apparent shear rates are no longer proportional to 1/H but, rather, composed of two straight lines joined at 1/H=1.33 mm -1 (or 0.75 mm equivalent gap), therefore implying two ranges for the use of Eq.10. At 1/H values below 1.33mm -1 (blue region, Fig. 6) the chocolate behaves as a liquid-like material without any slip velocity. This behavior is identical to that of Fig. 5c where the slip velocity term in Eq.10 is negligible.

Above 1.33 mm -1 (red region), the chocolate may take on a partially liquid and also a partially solid form so both the two terms of Eq.10 must be considered when data are adjusted. For instance, at 6.6 Pa, the fit of experimental data gives the slip velocity equal to 0.03m/s and a bulk shear rate of 0.05s -1 . The reciprocal length 1/H=1.33mm -1 is consequently the separating line between these two different behaviors or regimes and is surely a consequence of confinement of the particles which measure approximately 20 µm. This same transition 1/H=1.33 mm -1 corresponds to H=750 µm, indicating that the gap is only ~40 times the size of an average particle and therefore relatively small. This value seems to be the limit for distinguishing between an "equivalent granular" or alternatively a continuous view of the molten chocolate used here. The "equivalent granular" material flows more like a solid while the continuous material will flows more like a liquid [START_REF] Liu | Nonlinear dynamics: Jamming is not just cool any more[END_REF]. Without going further in the interpretation of such behavior, it can be argued that the transition driven by the 1/H value, i.e. by the gap size, is due to a localized "cooperative flow event". Such an event affects the fluidity which is highly dependent on the ratio of H to particle size [START_REF] Coussot | Yield stress fluid flows: A review of experimental data[END_REF][START_REF] Goyon | How does a soft glassy material flow: finite size effects, non local rheology , and flow cooperativity[END_REF].

A second important parameter that drives chocolate flow or not is the experiment duration. For experiments of short duration, the material will be elastically distorted but will not flow. For experiments which are longer, the material will first be distorted and then, depending on the stress applied, will flow after either a short or a substantially longer time as shown in the existing literature [START_REF] Auffret | Stress and strain controlled rheometry on a concentrated lyotropic lamellar phase of AOT/Water/Iso-octane[END_REF][START_REF] Benmouffok-Benbelkacem | Nonlinear viscoelasticity and temporal behavior of typical yield stress fluids: Carbopol, Xanthan and Ketchup[END_REF][START_REF] Caton | Plastic behavior of some yield stress fluids: from creep to long-time yield[END_REF][START_REF] Cheddadi | Steady Couette flows of elastoviscoplastic fluids are nonunique[END_REF][START_REF] Coussot | Yield stress fluid flows: A review of experimental data[END_REF][START_REF] Ovarlez | Rhéologie des pâtes granulaires[END_REF][START_REF] Ovarlez | Influence of shear stress applied during flow stoppage and rest period on the mechanical properties of thixotropic suspensions[END_REF][START_REF] Saramito | A new constitutive equation for elastoviscoplastic fluid flows[END_REF].

To conclude this section, Fig. 6 summarizes what must be borne in mind and what interpretation is to be given to the symbols on the right. For low shear stresses (in green), the flow chocolate behaves as a solid-like slippage at the wall (insert (a), Fig. 6). For high or intermediate shear stresses and small 1/H (in blue), the chocolate behaves as a liquid-like material without any slip velocity at the walls (insert (c), Fig. 6). For intermediate shear stresses (in red), the chocolate behaves partially as liquid and also partially as solid material showing simultaneously the behaviors of both bulk flow and slip at the walls (insert (b), Fig. 6). This latter case is interpreted to be a consequence of the non-homogeneous shear rate distribution within the gap induced by the progressive fluidization of the yield stress fluid. As a consequence, a compromise has to be found between the required level of characterization and the corresponding period available to carry out the assays. In fact, the cost necessary to accurately determine the yield stress value will be proportional to the measurement duration. From an industrial point of view, applying the present method to a few well selected gap values, e.g. two or three, may be sufficient, particularly if a well-defined programming procedure is used. By contrast, to distinguish between products possessing very similar behaviors, e.g. two molten chocolates of similar composition and particle size, a longer duration and a greater number of experiments will be required.

From a scientific point of view, the generalization of this method could help to provide some clues to a deeper understanding of flow types close to the yield stress. In addition, the tool presented here can also be used in numerical studies to set the boundary conditions at solid walls.

Conventional Yoshimura and Prudhomme correction

Based on the continuum hypothesis of shear stress within the gap (cf 3.1.1), the conventional Yoshimura and Prudhomme method [START_REF] Yoshimura | Wall Slip Corrections for Couette and Parallel Disk Viscometers[END_REF] proposed the use of two series of experimental data carried out at H1 and H2 gap sizes, which allows the calculation of the apparent shear rate and viscosity at a given shear stress :

(11) (12)
Slip velocity may also be directly deduced:

(13) Hereinafter, we will apply these formula to the interpolated and shear stress corrected data in §3.1.1.

Wall slip velocity

In Fig. 7, the solid triangles represent the wall slip velocity in relation to the shear stress (Pa). The wall slip velocity was obtained from the data represented in Fig. 5a to 5c and by using Eq.10. Below 5.4Pa in the green region, the slip velocity follows an affine function of slope =4.55×10 -3 ±4.11×10 -5 mm/(Pa.s). For these stresses, the chocolate behaves as a solid-like material. Above 8Pa in the blue region, the slip velocity symbolized by the thick solid line at the bottom is negligible. In this domain, the chocolate is fully sheared without any slip velocity at the walls. ). The open circle symbols represent the wall slip velocity evolution obtained using the Yoshimura and Prud'homme method with 0.50-0.75 mm; 0.50-1.50 and 1.00-1.50 mm gaps pair.

Between 5.4Pa and 8Pa, two slip velocities are obtained. For a constant shear stress at low 1/H values, there is no slip velocity as represented by the horizontal lines in Fig. 5b. These results are represented by the thick solid line in the blue region in Fig. 7. At high 1/H values, the fit of the linear shapes in Fig. 5b using the Eq.10, as previously explained, returns not only values of an effective shear rate but also a constant slip velocity of approximately 0.03mm/s whatever the applied stress, as represented by the triangles in the red region in In the existing literature, slip velocity is assumed to be either proportional [START_REF] Russel | Distinguishing between dynamic yielding and wall slip in a weakly flocculated colloidal dispersion[END_REF] or quadratically dependent on the shear stress [START_REF] Clasen | Determining the true slip of a yield stress material with a sliding plate rheometer[END_REF]Meeker et al., 2004b;[START_REF] Salmon | Towards local rheology of emulsions under Couette flow using Dynamic Light Scattering[END_REF][START_REF] Yeow | Non-Newtonian flow in parallel-disk viscometers in the presence of wall slip[END_REF]. However, at this stage there is no obvious explanation for these differences.

A path worthy of consideration could be seen in Fig. 7 by open circle symbols representing data obtained using the Yoshimura and Prud'homme method [START_REF] Yoshimura | Wall Slip Corrections for Couette and Parallel Disk Viscometers[END_REF] applied at 0.50-0.75 mm; 0.50-1.50 and 1.00-1.50 mm gaps pair, respectively. Clearly, the Yoshimura and Prud'homme method produces different results in relation to the gaps pair used.

Furthermore, the velocity curves show surprising concave shapes characterized by a maximum occurring in the intermediate shear stress region (red region). At low stresses (green region), the difficulty in reproducing the same slip velocity evaluation when the gaps pair is changed can be interpreted as the result of a combination of experimental errors and of the difference between the two gaps.

At intermediate stresses (red region), the decrease in the slip velocity is induced by the use of Eq.13 at two reciprocal gaps 1/H1 and 1/H2, within two different regions, for instance, the red region for 1/H1 and the blue region for 1/H2 in Fig. 6.

In this case, it becomes obvious that the necessity of following a constant shear stress line to join a shear rate couple, at 1/H1 and 1/H2 values respectively, will not be satisfied.

For these reasons, we do not recommend the use of the Yoshimura and Prud'homme method on all yield stress fluids in order to determine the slip vs stress evolution and even less so for conditions close to the yield stress.

Reconstitution of the steady state bulk flow curve

The effective viscosity vs the bulk shear rate figure can now be reconstituted (Fig. 8). The effective viscosity is computed by a simple division of the corrected shear stress ( §3.1.1) by the effective (or bulk) shear rate ( §3.1.2). In contrast to Fig. 1 the present work led to a yield stress fluid behavior depicted by a continuous decrease of the effective viscosity vs the bulk shear rate. The apparent Newtonian plateau followed by an apparent shear-thickening behavior in Fig. 1 completely disappeared when the corrections proposed in this article were applied. Nonetheless, this is not the case for our data (open circles, Fig. 8) to which we applied the Yoshimura and Prud'homme correction. Indeed, this last figure shows a Newtonian plateau for the different selected gaps pair. Moreover, the plateau value will depend on the selected gaps pair. Table 2 summarizes all obtained parameters from the fitting of corrected data from the IOCCC standard conditions (2-50s -1 ) and also for the extended conditions (0.005-50s -1 ) with Eq.14 and Eq15. From Fig. 9 and Table 2, it can be observed that the Herschel-Bulkley model is systematically above experimental data whereas the Casson model is systematically below them. Moreover, both models applied within the range of the IOCCC standard recommendation reveal a larger discrepancy. Indeed, as an example, the Herschel-Bulkley model gives a yield stress of 8.27±0.14Pa rather than the 5.4±0.2Pa directly determined in Fig. 5. It is slightly better for the Casson model which gives a yield stress of 4.23±0.10 Pa.

For both models, the viscosity, consistency and pseudo plastic index are summarized in Table 2 in addition to the yield stress value. All these estimations constitute useful data to characterize a yield stress fluid properly.

It should be noted that the evaluation of the yield stress by using both the Casson and Herschel-Bulkley models is of particular interest in the present case, as these models respectively under and overestimate it. This point could be a subject worthy of further study. Finally, we can also notice that the first stress value for which the molten chocolate begins to flow in Fig. 5b ( §3.1.2), represented by the horizontal line in Fig. 9 at 5.4Pa, is in agreement with the yield stress obtained by using the best fits.

Conclusion

In this article we have proposed a simple and reliable method particularly well-adapted for the characterization of food materials showing yield stress behavior. This method allows identification of the constitutive and slip velocity law, and consequently an accurate evaluation of rheological parameters, essential for material characterization such as viscosity, consistency, pseudo plastic index, or yield stress. The level of accuracy of such parameters increases if, and only if, (i) an extended range of reciprocal gap is used and (ii) if all collected data are in an apparently steady state. This last condition can take a long time to obtain close to the yield stress.

The originality of this method lies in its capacity to take into account the correction of the slip velocity simultaneously with the non-homogeneity of the flow field within the geometry gap. It also highlights the fact that the Yoshimura and Prud'homme method is inappropriate for yield stress fluid characterization.

Molten chocolate, given here as an example of an industrial multiphase product, was found to be a relatively simple material as it is not subject to aging, retarded viscoelasticity, thixotropy or sedimentation phenomena within the measurement duration. Nevertheless, the validity of our method remains to be demonstrated for more complex yield stress soft materials. Indeed, we can consider this method as applicable to model fluids as well as to more complex soft condensed shear thinning or viscoplastic materials or to industrial foodstuffs, which are not too heterogeneous at the scale of gap of the geometry. Provided, that a sufficient level of accuracy is reached, one can expect to be able to distinguish molten chocolates of similar formulations or other liquid foodstuffs. ). The open circle symbols represent the wall slip velocity evolution obtained using the Yoshimura and Prud'homme method with 0.50-0.75 mm; 0.50-1.50 and 1.00-1.50 mm gaps pair. Fig. 9: Evolution of corrected shear stress vs bulk shear rate (solid circle symbols). The lines correspond to the different data fittings using Casson (in blue) and Herschel-bulkley (in orange) as models. Broken lines refer to the IOCCC standard conditions, i.e. a 2s -1 -50s -1 fitting while solid lines refer to a fit over the entire range of data, i.e. a 0.005-50s -1 fitting.

List of Figure Captions

  Fig.1 showed an apparent Newtonian plateau followed by an apparent shear-thickening and a shear-thinning behavior when increasing shear rate values. If the shear stress vs shear rates are represented (Fig.2), one can see distinct but slightly curved lines in the lower part of the graph which change abruptly into a further much more accentuated curve. The pseudo Newtonian plateau and shear-thickening behavior in Fig.1 correspond to the lower part of the graph in Fig.2 while the shear-thinning behavior in Fig.1 is illustrated by the more accentuated curve in the second figure. At other gap values, similar behaviors were observed albeit with an abrupt change observed for the lower shear rate when increasing the gap size.

Fig. 2 :

 2 Fig.2:Examples of apparent shear stress vs shear rates curves obtained at 1.50mm, 1.00mm, 0.75mm and 0.50 mm gaps.

Fig. 3 :Fig. 4 :

 34 Fig.3: (a) Schematic drawing of parallel plate-plate geometry. (b) Focus on the velocity field inside the geometry gap in the presence of slip at the solid wall.

Fig.

  Fig.5a-c details the apparent shear rate vs the reciprocal gap at different

Fig. 6 :

 6 Fig.6: Apparent shear rate as a function of the reciprocal gap for different constant shear stresses represented by the solid lines. A summary of the different behaviors encountered on the left is symbolized by the three representations on the right: (a) at low shear rate (in green), a pure slip regime, (b) at intermediate shear rate (in red), a transition regime (c) at high shear rate (in blue), a homogeneous bulk shear regime.

Fig. 7 :

 7 Fig.7: Wall slip velocity as a function of the shear stress. The solid triangle symbols illustrate experimental data from the present work. The solid line represents the correspondingly best linear fit (). The open circle symbols represent the wall slip velocity evolution obtained using the Yoshimura and Prud'homme method with 0.50-0.75 mm; 0.50-1.50 and 1.00-1.50 mm gaps pair.

Fig

  Fig.7.

Fig. 8 :Fig. 9 :

 89 Fig.8: Evolution of the effective viscosity vs bulk shear rate. The solid line represents the present work's corrected data. Open circles represent the data obtained via the Yoshimura and Prud'homme method with 1.00-1.50mm; 0.50-1.50 and 0.50-0.75mm gaps pair.

Fig. 1 :

 1 Fig.1: Apparent viscosity vs apparent shear rate evolutions at different times for the same molten chocolate sample measured with a 1mm gap smooth plate-plate. Solid dots (0min) represent the first increasing ramp. Open symbols: triangles (5min), diamonds (1h23min) and squares (2h41), represent decreasing ramps. The symbols show respectively steps number 1,2,6 and 10 of the proposed protocol.

Fig. 2 :

 2 Fig.2:Examples of apparent shear stress vs shear rates curves obtained at 1.50mm, 1.00mm, 0.75mm and 0.50 mm gaps.

Fig. 3 :

 3 Fig.3: (a) Schematic drawing of parallel plate-plate geometry. (b) Focus on the velocity field inside the geometry gap in the presence of slip at the solid wall.

Fig. 4 :

 4 Fig.4: Corrected shear stress versus apparent shear rate curves. The solid symbols are the mean values of interpolated data. The solid lines are the mean values of both the corrected and interpolated data. All experiments were carried out with a smooth plateplate geometry at varying gap size: 1.8mm, 1.5mm, 1mm and 0.5 mm.

Fig. 5 :

 5 Fig.5: Apparent shear rate as a function of the reciprocal gap for different constant shear stresses. Details of each encountered behavior: (a) low shear rate: pure slip regime, (b) intermediate shear rate: transition regime (c) high shear rate: homogeneous bulk shear regime.

Fig. 6 :

 6 Fig.6: Apparent shear rate as a function of the reciprocal gap for different constant shear stresses represented by the solid lines. A summary of the different behaviors encountered on the left is symbolized by the three representations on the right: (a) at low shear rate (in green), a pure slip regime, (b) at intermediate shear rate (in red), a transition regime (c) at high shear rate (in blue), a homogeneous bulk shear regime.

Fig. 7 :

 7 Fig.7: Wall slip velocity as a function of the shear stress. The solid triangle symbols illustrate experimental data from the present work. The solid line represents the correspondingly best linear fit (). The open circle symbols represent the wall slip velocity evolution obtained using the Yoshimura and Prud'homme method with 0.50-0.75 mm; 0.50-1.50 and 1.00-1.50 mm gaps pair.

Fig. 8 :

 8 Fig.8: Evolution of the effective viscosity vs bulk shear rate. The solid line represents the present work's corrected data. Open circles represent the data obtained via the Yoshimura and Prud'homme method with 1.00-1.50mm; 0.50-1.50 and 0.50-0.75mm gaps pair.

Table 1

 1

	(a) IOCCC2000 and ICA standards.
	(b) Rheometrical protocol and corresponding time sequences.
		Pre-shear	Increasing	Constant	Decreasing
			shear rate	high shear	shear rate
				rate
	IOCCC 2000	5s -1 on 500s 2 to 50s -1	50s -1	50-2s -1
	ICA	5s -1 on 300s	on 180s	on 60s	on 180s
			(a)

b) Fig.1: Apparent viscosity vs apparent shear rate evolutions at different times for the same molten chocolate sample measured with a 1mm gap smooth plate-plate. Solid

	dots
	(0min) represent the first increasing ramp. Open symbols: triangles (5min), diamonds
	(1h23min) and squares (2h41), represent decreasing ramps. The symbols show
	respectively steps number 1,2,6 and 10 of the proposed protocol.

Table 2 :

 2 Summary of the different parameters obtained by fitting the molten chocolate data with the Herschel-Bulkley and Casson models. The symbol * indicates deviation from the bulk yield stress 5.4Pa determined directly in Fig.5.

	Shear rate range	Corrected data
	for fittings	Herschel-Bulkley	Casson
			0,HB=8.27±0.14Pa (1.7 %) (53%)*	0,c=4.23±0.10 Pa (2.3%) (22%)*
	Standard	conditions	=2.38±0.06Pa.s n (2.3%)	c=0.94±0.01 Pa.s (1.3%)
	(2-50s -1 )		nHB=0.876±0.006 (0.7%)
	Extended conditions	0,HB=6.05 ±0.08Pa (1.3%) (12%)*	0,c=5.10±0.08 Pa (1.7%) (5.6%)*
	(0.005-50s -1 )	=3.39±0.08Pa.s n (2.4%)	c=0.85±0.01 Pa.s (1.3%)
			nHB=0.780±0.008 (1.0%)

Acknowledgments

The Laboratoire Rhéologie et Procédés is part of the LabEx Tec 21 (Investissements d'Avenir -grant agreement n°ANR-11-LABX-0030) and of the PolyNat Carnot Institut (Investissements d'Avenir -grant agreement n°ANR-11-CARN-030-01). Authors would like to gratefully acknowledge Susan Sinclair for proofing the manuscript.


P Gentile
email: pascal.gentile@cea.fr

T David

F Dhalluin

D Buttard

N Pauc

M Den Hertog

P Ferret

T Baron

Growth of small-diameter silicon nanowires to nanotrees

Keywords: Chemical vapour deposition, nanowires, Silicon

In this work we have studied a way to control the growth of small diameter silicon nanowires by the Vapour Liquid Solid (VLS) mode. We have developed a method to deposit colloids with good density control, which is a key point for the control of the nanowires (NWs) diameter. We also show the high dependence of the allowed growth diameter with growth conditions opening the door to the realization of as grown 2 nm silicon NWs. Finally we developed a smart way to realise nanotrees in the same run, by tuning the growth conditions and making benefit of gold on the sidewall of nanowires, without two catalyst deposition steps.

Introduction

The growth of semiconducting nanowires is now of a great interest for nanoelectronics, biosensors, and solar cells technology, due to their particular electrical and optical properties arising from their one dimensional geometry [1]. In particular, small diameter (d) silicon nanowires (SiNWs) and nanotrees (SiNTs) are viewed as attracting building blocks for future devices operating in the quantum limit predicted for d<5 nm wires [2,3]. To fabricate devices with homogeneous properties, getting monodispersed diameter populations is also a prerequisite. In this communication, we present a method to grow silicon nanotrees (Si NTs) with small diameter branches (<10 nm) using a two step process. This method comes from our studies which showed that (i) in case of larger diameter nanowires (> 50 nm) we have observed gold droplets on Si NWs sidewalls, based on in situ gold diffusion, (ii) the growing condition could be adjusted to grow selectively small or large NWs. By combining the catalytic activity of the Au particles at the sidewalls with the selective growth condition for different diameter we can realise NTs with branches presenting diameter down to 2 nm and reduced size dispersion. To sum up the presented method, the first step process consists in growing large NWs (>50 nm) with Au droplets at the sidewalls and the second step consists in the selective growth of branches from these Au sidewalls droplets.

Experimental details

Silicon nanowires are grown by Low Pressure Chemical Vapour Deposition (LPCVD) via the VLS mechanism [4][5][6][7], in which NWs grow from gold catalyst. The silicon substrate orientation (111) is chosen to obtain vertical growth for large NWs as their growth direction is <111>.

Silane (SiH4) diluted in hydrogen is used as the reactive gas. In this method, key parameters that control the NWs characteristics are the catalyst size and the growth conditions (temperature, silane pressure). To achieve control of the catalyst size, we tested two deposition procedures, i.

e. gold film dewetting and direct deposition of size calibrated catalysts from a gold colloids solution (Ted Pella) onto a silicon surface. For both methods, silicon surface is first cleaned in alcohol and deoxidized into an ammonium fluoride mixture before catalyst deposition. The first family of samples is immediately coated with a thin gold layer (< 2nm) in a joule effect evaporator, whereas "colloidal" samples are placed in an electrochemical cell for colloids deposition. In this second approach, negatively charged gold colloids (2, 5, 10, 20, 50, 100 and 200 nm diameter) dispersed in water are deposited in a controlled fashion onto the sample surface by using an external electric field as a driving force for colloids motion. This technique provides a fine control of the catalysts density over large surfaces by adjusting electrical fields or deposition times as shown in table 1. Samples are then introduced in the LPCVD reactor, and annealed under H2 flow. To form efficient catalyst from a thin gold layer, in situ annealing (between 650°C and 850°C, 5 minutes) is required before crystal growth to induce dewetting of the films. It produces droplets with large size dispersion as compared to the colloidal way, making fine control of the diameter distribution harder. Note that for colloid deposited samples an annealing step in the 450°C -750°C range is performed, but this treatment was not found to be critical to achieve efficient crystal growth. Temperature is then decreased down to the growth condition and SiH4 reactive gas is introduced into the reactor, initiating the SiNWs growth.

Growth modes study

In this section, growth studies are carried out to obtain (i) large diameter (> 50 nm) nanowires with Au at the sidewalls, (ii) very small diameter NWs diameter (< 10 nm) and (iii) nanotrees using a combination of the former points.

Large diameter nanowires

Optimum growth conditions for regular and vertical growth of large diameter NWs (from d> 50 nm to d≈ 1 µm, see in figure 1) on gold coated Si(111) were obtained for temperature between 550°C to 650°C and H2 and SiH4 flow rates of 4 l/mn and 10 to 25 standard cube cm per minute (sccm) respectively. The total base pressure was held at 20 mbar. The same conditions are used for monodispersed colloidal catalyst and a high growth yield for d>50 nm is demonstrated. For these conditions, we noticed (i) an inhibited growth of small (less than 70 nm diameter) NWs' diameter (see figure 1b) and (ii) the formation of gold droplets (d≈5-10 nm) at the sidewalls NWs as shown by the bright spots on the high magnification SEM view in figure 1c, which is corroborated by the observation from J. B. Hannon et al [8]. Size homogeneity and chemical composition of these nanodroplets have been also studied by Scanning Transmission Electron Microscopy (STEM) and Energy Dispersive X-ray (EDX) analysis [9], which confirmed the gold content in these nanodroplets and their tight size distribution. We want also to point out the fact that the growth of some NWs stops as shown in figure 1b. We attribute this effect either to the gold depletion at the Au/Si catalyst during SiNWs growth for the used experimental conditions (temperature, silane partial pressure) as shown by J. B. Hannon et a l [8] or by a poisoning of the Au/Si catalyst. This effect sets the lowest NWs' diameter for fixed temperature and silane partial pressure independently on the size of the original Au catalyst. We found that this phenomenon is not dependant on the fabrication method of the catalyst (dewetted films or colloids). Increasing the growth temperature favours large diameter NW growth (in agreement with temperature activated diffusion of gold) and high growth rate.

Small diameter nanowires

To reduce the gold diffusion the growth temperature (TG) was decreased and in the same time the silane flux was enhanced to conserve an equivalent growth rate. Then w e found that decreasing the TG and enhancing the silane flux give the best results to grow small diameter NWs. To study the diameter dependence of SiNWs as a function of the growth conditions, we used various size gold colloids (2, 5, 10 20 and 50 nm diameter) and we varied the TG and silane flow. SEM side views of SiNWs synthesized with 50 nm colloids are reported in figure 2 for some typical process points. We observe inhibited growth and catalyst loss on the top of NWs at 600°C and 10 sccm silane flux (figure 2(a)). The growth efficiency increased progressively as temperature decreased and silane flux increased as indicated by the following series:

550°C/silane 20 sccm (figure 2b), 450°C/silane 20 sccm (figure 2c). For smaller diameter colloids (2 nm to 20 nm), we found optimum growth conditions for lower temperatures (in the 400°C to 500°C range) and silane flux in the 35 sccm to 200 sccm range as shown in table 2.

Low temperature and rich silane atmosphere are then necessary conditions to get small diameter nanowires. High resolution transmission electron microscopy (HRTEM) images are realised on 20 nm SiNWs (figure 3). These two images show crystalline nanowires with different growth direction <111>

and <110> [7].

Adjusting the temperature and the silane partial pressure leads to a controlled growth of ultra small NWs diameter as seen in the figure 4. The size dispersion observed is one to three times the nominal diameter of the gold colloids corresponding to the size dispersion of the gold colloids solution. This suggests that for the small diameter conditions no catalyst diffusion and coalescence occur onto the sample at the annealing step, making the starting diameter of the wires only dependent on the colloids size dispersion. Using small diameter calibrated colloids, NWs as small as 2 nm in diameter were routinely synthesised, which opens up the door to the realisation of devices where quantum effects play a significant role. It is worthwhile noting that our method offers a straightforward way to get ultra-small NWs without the need to reduce the diameter by thermal oxidation in an uncontrollable fashion for these nanostructures [10].

To conclude this study, we carried out similar experiments on gold coated samples and found the same NW diameters evolution with the temperature and the silane partial pressure (PSiH4). As an example, we grew SiNWs under small diameter conditions (400°C and 200 sccm silane) from a dewetted 0.7 nm gold film. We obtained a population of small NWs with diameters ranging from 5 nm to 50 nm (see figure 5). Note that low annealing temperatures are required to obtain small diameter NWs due to the limited ripening of gold droplets, which results in small diameter catalysts.

Discussion

To sum up the influence of the experimental conditions, a high TG and low PSiH4 leads to the growth of large NWs diameter, whereas small NWs diameter are only form ed at low TG and high PSiH4. In addition, we want to point out that growth of NW from large diameter catalysts under small diameter conditions produces not only large diameter wires but also small specimens at their base (see in figure 6), making thermodynamic conditions leading parameters for the growth diameter. We propose two mechanisms to explain this result based on the formation of small catalyst. First, during the annealing phase, the large colloids particles (50 nm) could separate and form also small particles. Second, under high PSiH4 and low TG, a local Si supersaturation in the large Au particle could be obtained and formed an efficient seed for the growth of small Si NWs. Indeed, by reducing TG, the Si diffusion length is reduced leading to the formation of Si rich Au/Si droplets at the surface of the Au particle. Then in these two cases, at high PSiH4 these droplets have a radius superior to the critical radius necessary for the occurence of Si NWs growth [11].

Nanotrees

From the combination of (i) the presence of small gold droplets on the sidewalls on the NWs for the larger diameter growth conditions and (ii) the selective control of growth conditions for small diameter, we propose a way to grow nanotrees with NWs branches presenting very small diameter (<10 nm) without post-growth catalyst deposition. Up to now [12,13], growth of nanotrees is a two steps process which combines initial catalyst deposition for the trunk growth and a post-growth deposition of seeding material on as grown trunk sidewalls to generate the catalyst for the branches growth.

In our method, at first, trunks are made under large diameter growth conditions (high temperature, low silane partial pressure) from 200 nm gold colloids to obtain 200 nm diameter SiNWs. The selectivity of the process we have developed allow us to grow a low dispersed NWs population in diameter for every kind of diameter. Growth conditions are then changed within the same run (low temperature, high silane partial pressure) to favour growth of small NWs diameter in order to get branches from the gold diffused droplets on the trunk sidewalls.

Figures 7a and7b show a typical view of a nanotree, with a 200 nm diameter trunk and branches with diameters lower than 10 nm. The branch diameter distribution is re markably restricted within the 5 nm to 10 nm range, which fits with the size of the gold droplets as shown in figure 1c and unambiguously proves that diffused droplets catalyse the branch growth.

By using colloidal particle or patterned metal as catalyst for the trunk growth, branches growth occurs only on nanowires and not on substrates, as it is the case for already published method where a two steps catalyst deposition leads to uncontrolled growth over the entire free surface [12,13]. This parasite growth must be eliminated by an ex-situ process [13] which uses a polymer mask at the trunk base to remove the redeposited catalyst from the free surface by liftoff. In addition, our method offers a single run process to produce large amounts of nanotrees with small SiNWs diameter branches with a low size dispersion, without supplemental oxidation and contamination problems resulting from air exposure and catalyst redeposition processes (in case of ex-situ processing). Similarly, we also obtained nanotrees on gold dewetted samples, but we noticed gold droplets formation on the whole sample surface presumably due to the large amount of deposited catalyst and here small diameter NWs grow on the free surface during the second growth step.

In figure 7a, we also observed a striking phenomenon: branch growth is accompanied by simultaneous regrowth at the trunk summit of a branch free section. Indeed, as shown in figure 6, under small diameter (<10 nm) growth conditions we also obtained a growth of larger NWs diameter. This demonstrates that no gold droplet forms on the regrowth section since no branch is observed at the tree summit. Moreover, we noticed a diameter reduction of the summit in comparison with the underlying trunk section. This remarkable size modification is probably due to the reduction of the catalyst diameter during the temperature ramp down at the end of the trunk growth (silane flow stopped, H2 atmosphere). Gold free NW sidewalls in the regrowth section obtained at lower process temperatures contrast with the rich gold trunk sidewalls seen at high temperatures. This trend suggests a thermally assisted diffusion mechanism of gold as the origin of the trunk section reduction and as the source of the catalyst droplets along the trunk sidewalls.

Conclusion

In summary, we have fabricated SiNWs down to 2 nm diameter and we have demonstrated different modes for SiNW growth in the gold assisted VLS method. By modulating and adjusting the growth conditions within the same run, we have developed a new and origin al method to synthesize nanotrees without the need for a second catalyst deposition to seed the branches. Our results indirectly give evidence of diffusion and formation of gold nanodroplets on the SiNWs sidewalls with a uniform and small (lower than 10 nm) diameter. Further studies and thermodynamic calculation are needed to explain these different growth modes [11] for improving the size distribution, density and length of branches.

Figure 1 Figure 2 Figure 3 Figure 4 Figure 5

 12345 Figure 1 SEM images of SiNWs grown at 650°C, 20 mbar, silane 15 sccm, (b) zoom of small diameter NWs with stopped growth during process, (c) zoom of NWs sidewall showing gold droplets.

Figure 6

 6 Figure6SEM image of SiNWs grown from 50 nm diameter gold colloid catalyst at 430°C and 50 sccm silane. We can observe small diameter NWs (from 5 nm to 10 nm) at the base of the larger NWs.

Figure 7

 7 Figure 7 SEM images of nanotrees grown from 200 nm diameter trunks. The branches are grown from gold droplets present on the trunk sidewall. In (a) we can see the overgrowth NW on the top with diameter reducing, in (b) zoom of nanotree.

Table 1

 1 (1) Density of 50nm gold colloids as function of applied voltage for 2nm.

	1	Density as function of Applied Voltage for 50 nm and 2mn colloids
	Voltage (V)	0	20	30	40	50	60
	Density (cm -3 )	0,00E+00	2,50E+08	5,00E+08	6,00E+08	8,00E+08	1,10E+09
	2	Density as a function of time for 20nm colloids and 40V
	Time (mn)	0,5	1	2	4	6	10
	Density (cm -3 )	1,00E+09	1,20E+10	2,00E+10	2,80E+10	3,30E+10	5,00E+10

Table 2

 2 Growth temperature and silane flux for the different nanowires diameters.

	diameter (nm)	2	5	10	20
	Temperature °C	400 -430	400 -430	430 -450	450 -500
	silane flux (sccm)	200	100 -200	100 -50	50 -35


Miriam Reverter
email: miriam.reverter@nuigalway.ie

Marie-Aude Tribalat
email: marieaude.tribalat@gmail.com

Thierry Pérez
email: thierry.perez@imbe.fr

Olivier P Thomas
email: olivier.thomas@nuigalway.ie

Metabolome variability for two Mediterranean sponge species of the genus Haliclona: specificity, time, and space

Keywords:

Introduction

The study of natural variation of metabolites brings valuable information on the physiological state of the organisms as well as their phenotypic traits. In marine organisms, metabolome variability has mostly been addressed through targeted studies on metabolites of ecological or pharmaceutical interest. 1 2 3 1,2✉ 1 2 3 8/25/2018 2/26

However, comparative metabolomics has demonstrated its potential to address the overall and complex metabolic variability of organisms.

Objectives

In this study, the intraspecific (temporal and spatial) variability of two Mediterranean Haliclona sponges (H. fulva and H. mucosa) was investigated through an untargeted and then targeted metabolomics approach and further compared to their interspecific variability.

Methods

Samples of both species were collected monthly during 1 year in the coralligenous habitat of the Northwestern Mediterranean sae at Marseille and Nice. Their metabolomic profiles were obtained by UHPLC-QqToF analyses.

Results

Marked variations were noticed in April and May for both species including a decrease in Shannon's diversity and concentration in specialized metabolites together with an increase in fatty acids and lyso-PAF like molecules. Spatial variations across different sampling sites could also be observed for both species, however in a lesser extent.

Conclusions

Synchronous metabolic changes possibly triggered by physiological factors like reproduction and/or environmental factors like an increase in the water temperature were highlighted for both Mediterranean Haliclona species inhabiting close habitats but displaying different biosynthetic pathways. Despite significative intraspecific variations, metabolomic variability remains minor when compared to interspecific variations for these congenerous species, therefore suggesting the predominance of genetic information of the holobiont in the observed metabolome.

Introduction

Metabolites are the products of metabolic reactions within an organism and their levels reflect the response of biological systems to physiological and environmental changes. The metabolome represents a snapshot of the metabolites present in an organism at a given time and will change according to different life stages and environmental conditions [START_REF] Kosmides | Metabolomic fingerprinting: Challenges and 8/25/2018 22/26 opportunities[END_REF][START_REF] Goulitquer | Mass spectrometry-based metabolomics to elucidate functions in marine organisms and ecosystems[END_REF]. The variability of the metabolome can bring very useful information in ecology but also in systems biology. For instance, studying the natural variation of metabolites proved to be useful in unravelling the evolution of phenotypic traits in terrestrial plants [START_REF] Li | Navigating natural variation in herbivory-induced secondary metabolism in coyote tobacco populations using MS/MS structural analysis[END_REF][START_REF] Routaboul | Metabolite profiling and quantitative genetics of natural variation for flavonoids in Arabidopsis[END_REF]. In marine organisms, study of metabolome variability mostly focused on so-called secondary or specialized metabolites [START_REF] Kuhlisch | Metabolomics in chemical ecology[END_REF], low molecular weight (< 1500 Da) molecules playing key roles in various ecological processes [START_REF] Agrawal | Insect herbivores drive real-time ecological and evolutionary change in plant populations[END_REF]Ferrer andZimmer 2013, 2012). These metabolites being involved in chemical mediation between organisms, their variation does not only affect the resilience of one species but it can have cascading effects on the networks of interactions [START_REF] Hay | Challenges and opportunities in marine chemical ecology[END_REF][START_REF] Bornancin | Chemical mediation as a structuring element in marine gastropod predator-prey interactions[END_REF]. Intraspecific variation of metabolites with time and space is crucial for understanding species ecology and ecosystem functioning [START_REF] Paul | Chemical defenses: From compounds to communities[END_REF][START_REF] Glassmire | Intraspecific phytochemical variation shapes community and population structure for specialist caterpillars[END_REF]).

In the marine environment, spatial and temporal variations of specialized metabolites have been studied in algae [START_REF] Payo | Variability of non-polar secondary metabolites in the red alga Portieria[END_REF][START_REF] Soares | Intraspecific variation of meroditerpenoids in the brown alga Stypopodium zonale guiding the isolation of new compounds[END_REF], cnidarians [START_REF] Kelman | Variation in secondary metabolite concentrations in yellow and grey morphs of the red sea soft coral Parerythropodium fulvum fulvum: Possible ecological implications[END_REF][START_REF] He | NMR-Based Metabolomic Analysis of Spatial Variation in Soft Corals[END_REF][START_REF] Costa-Lotufo | Chemical profiling of two congeneric sea mat corals along the Brazilian coast: Adaptive and functional patterns[END_REF], tunicates [START_REF] López-Legentil | Temporal variation in the production of four secondary metabolites in a colonial ascidian[END_REF], bryozoans [START_REF] Peters | Variation of brominated indoles and terpenoids within single and different colonies of the marine bryozoan Flustra foliacea[END_REF]) and sponges [START_REF] Reverter | Secondary metabolome variability and inducible chemical defenses in the Mediterranean sponge Aplysina cavernicola[END_REF][START_REF] Ternon | How environmental factors affect the production of guanidine alkaloids by the Mediterranean sponge Crambe crambe[END_REF][START_REF] Sacristan-Soriano | Relevant spatial scales of chemical variation in Aplysina aerophoba[END_REF]. Spatial variation of metabolites has been mostly related to local adaptation to the environment and fluctuating parameters such as temperature, light, or chlorophyll concentrations. In a recent study, [START_REF] Soares | Intraspecific variation of meroditerpenoids in the brown alga Stypopodium zonale guiding the isolation of new compounds[END_REF] found different chemotypes of the brown alga Stypopodium zonale in different geographic locations, probably indicative of genotypic sorting across habitats. Seasonal changes of specialized metabolites possibly related to changes of environmental parameters such as temperature have been observed in some marine invertebrates [START_REF] Page | Spatial and temporal variability of cytotoxic metabolites in populations of the New Zealand sponge Mycale hentscheli[END_REF][START_REF] Sacristán-Soriano | Temporal trends in the secondary metabolite production of the sponge Aplysina aerophoba[END_REF][START_REF] Ternon | How environmental factors affect the production of guanidine alkaloids by the Mediterranean sponge Crambe crambe[END_REF]. According to the optimal defence theory (ODT) [START_REF] Rhoades | Offensive-defensive interactions between herbivores and plants: Their relevance in herbivore population dynamics and ecological theory[END_REF], the energy allocated to the biosynthesis of specialized metabolites results from the continual trade-off with the energy invested in primary functions such as growth and reproduction. Several studies on marine invertebrates have indeed shown a decrease in the production of specialized metabolites during reproduction [START_REF] López-Legentil | Secondary metabolite and inorganic contents in Cystodytes sp. (Ascidiacea): Temporal patterns and association with reproduction and growth[END_REF][START_REF] Ivanisevic | Biochemical trade-offs: Evidence for ecologically linked secondary metabolism of the sponge Oscarella balibaloi[END_REF]. However, this phenomenon may not be universal in all marine invertebrates, and might be related to the reproduction strategy and associated metabolic cost [START_REF] Reverter | Secondary metabolome variability and inducible chemical defenses in the Mediterranean sponge Aplysina cavernicola[END_REF].

Sponges (phylum Porifera) are important components of benthic communities worldwide where they play major roles in some ecosystem functioning like the Caribbean Sea or the Central Indo-Pacific [START_REF] De Goeij | Surviving in a marine desert: The sponge loop retains resources within coral reefs[END_REF][START_REF] Haas | Global microbialization of coral reefs[END_REF]. Sponges are also well known to produce a large diversity of bioactive specialized metabolites that confer protection against predators and pathogens and advantage against competitors [START_REF] Proksch | Defensive roles for secondary metabolites from marine sponges and sponge-feeding nudibranchs[END_REF][START_REF] Hay | Marine chemical ecology: Chemical signals and cues structure marine populations, communities, and ecosystems[END_REF][START_REF] Blunt | Marine natural products[END_REF][START_REF] Loh | Chemical defenses and resource trade-offs structure sponge communities on Caribbean coral reefs[END_REF]. Within the class Demospongiae, several sponges of the order Haplosclerida are considered as prolific sources of bioactive marine natural products. In particular, species from the genus Haliclona (family Chalinidae) have attracted much attention due to the high diversity of species and the variety of specialized metabolites they contain [START_REF] Tribalat | Does the chemical diversity of the order Haplosclerida (Phylum Porifera: Class Demospongia) fit with current taxonomic classification?[END_REF]. Importantly, the systematics of this genus and more generally the Haplosclerida group is being revised deeply. Haliclona species are found throughout the Mediterranean, from intertidal to continental shelf habitats (2 to > 100 m) [START_REF] Santín | Sponge assemblages on the deep Mediterranean continental shelf and slope (Menorca Channel, Western Mediterranean Sea)[END_REF].

In this study we focused on the Mediterranean species Haliclona fulva and Haliclona mucosa, two sponge species found in shallow waters (10-50 m) in semi-dark conditions. Haliclona (Halichoclona) fulva is known to produce several acetylenic compounds (renierins, fulvinol, and fulvynes) and sesquiterpenes (fulvanins) [START_REF] Cimino | New acetylenic compounds from the sponge Reniera fulva[END_REF][START_REF] Casapullo | Paniceins and related sesquiterpenoids from the Mediterranean sponge Reniera fulva[END_REF][START_REF] Ortega | Fulvinol, a new long-chain diacetylenic metabolite from the sponge Reniera fulva[END_REF][START_REF] Nuzzo | Fulvynes, antimicrobial polyoxygenated acetylenes from the Mediterranean sponge Haliclona fulva[END_REF]. Recently, cladocroic acid and the diastereomer 3-epicladocroic acid were also isolated from H. fulva (Genta-Jouve and Thomas 2013). H. mucosa (subgenus Soestella) is known to produce sesquiterpenes hydroquinones from the panicein family [START_REF] Zubía | Sesquiterpene hydroquinones from the sponge Reniera mucosa[END_REF]) and the bicyclic eicosanoid mucosin [START_REF] Casapullo | Mucosin: A new bicyclic eicosanoid from the Mediterranean sponge Reniera mucosa[END_REF]). Thus, these two species, phylogenetically related and living in the same environment, produce totally different families of specialized metabolites. The objective of this study was (1) to assess the intra-specific variability (temporal and spatial) of these two species metabolomes, in order to investigate if despite their different biosynthetic pathways they display similar patterns of variation and (2) to evaluate which type of intra-specific variation is more important (spatial or temporal) and their relative contribution when compared to inter-specific variation. To reach these objectives, we first deployed a non-targeted UHPLC-HRMS metabolomics approach to investigate patterns of metabolome variability.

Materials and methods

Biological material

Specimens of Haliclona fulva and Haliclona mucosa were collected by SCUBA between 20 and 25 m depth in the Northwestern Mediterranean Sea. For the temporal variability, 6 individuals of each species were collected monthly at the Grotte du Lido in the bay of Villefranche-sur-Mer (N43°41′30.0′′, E7°19′11.9′′) from April 2014 until March 2015. For the spatial variability, H. mucosa was also sampled in Grotte de Jarre (N43°11′52.0′′, E5°21′53.8′′; n = 4) and Riou Caramasseigne (N43°10′19.3′′; E5°23′37.4′′; n = 3) near Marseille and H. fulva in the Grotte de Jarre (Marseille, n = 3) and Monaco (N43°43′34.9′′, E7°25′15.2′ ′; n = 4) in August 2014 (Supp. Fig. S1). Samples were immediately frozen at -20 °C, freeze-dried and stored dry at - 20 °C until extraction.

Chemical extraction

The freeze-dried samples were first ground to obtain a homogeneous powder. A mass of 250 mg of powder was extracted 3 times in 5 mL of a mixture methanol/dichloromethane (1:1) for 5 min in an ultrasonic bath. The combined organic extracts were centrifuged 10 min at 2000 rpm and the supernatant was collected carefully and dried using a Speedvac. The extracts were dissolved in 5 mL of methanol and stored at 4 °C for 24 h. A volume of 0.1 mL of the supernatant was then filtered through 0.2 µm polytetrafluoroethylene (PTFE) and diluted with 0.9 methanol for the UHPLC-MS analysis. The quality control (QC) samples were prepared by pooling small aliquots from all extracts.

UHPLC-HRMS analyses

UHPLC-HRMS analyses were performed with a UHPLC-DAD equipped with an Ultimate 3000 (ThermoScientific, MA, USA) connected to a qToF Impact II mass spectrometer (Bruker). Chromatographic separation was achieved on a Nucleodur® Phenyl-hexyl column (100 × 4.6 mm, 1.8 µm, Macherey-Nagel, Germany) using a linear elution gradient of H O/Acetonitrile both solvents containing 0.01 M ammonium formate and 0.1% formic acid, from 95:5 (v/v) to 0:100 (v/v) in 12 min followed by 5 min at 100% acetonitrile. Mass spectra were acquired in positive mode with an elution flow 0.5 mL min and a constant column temperature of 40 °C. Injection volume was set to 10 L. The mass spectrometer analyser parameters were set as follows: nebulizer sheath gas, N (44 psi); dry gas, N (12 L/min); capillary temperature, 200 °C; capillary voltage, 2500 V; collision energy, 4 eV and a full-scan mass window of 50-1200 m/z. Samples sequence was randomized to reduce systematic error associated with instrumental drift. QC (pool) and blank (MeOH) samples were injected at the beginning, the end and equidistantly throughout the sequence.

AQ1

Metabolomics data have been deposited to the EMBL-EBI MetaboLights database (https://doi.org/10.1093/nar/gks1004. PubMed PMID: 23109552) with the identifier MTBLS638.

The complete dataset can be accessed here https://www.ebi.ac.uk/metabolights/MTBLS638.

Data treatment and statistical analysis

LC-MS raw data files were converted to mzML files using the open-source msConvert tool from the ProteoWizard library [START_REF] Kessner | ProteoWizard: Open source software for rapid proteomics tools development[END_REF]. mzML files were processed using the R package XCMS (R version 3.3.1., XCMS version 1.50.0) to detect, deconvolute and align features (molecular entities with a unique m/z and a specific retention time). Parameter for XCMS processing were set as described in [START_REF] Patti | Meta-analysis of untargeted metabolomic data: Combining results from multiple profiling experiments[END_REF] for UPLC-QTOF (high resolution). XCMS analysis of these data provided a matrix containing the retention time, m/z value and integrated peak areas of the identified features. Finally, only features with a coefficient of variation less than 20% among the QC samples were kept for further data analyses. Data were normalized by log transformation prior to statistical analysis.

Metabolite diversity of each sample was determined by the Shannon's diversity index (H). ANOVA test and Tukey's HSD post-hoc test were used to detect temporal and spatial differences in metabolite diversities of both species. Principal Component Analysis (PCA) was used to visualize the spatial and temporal (monthly) variability of both sponge metabolomes. Specialized metabolites were annotated with an in-house library of metabolites from H. fulva and H. mucosa, and variation of their expression level tested with Kruskal-Wallis non parametric analysis and Kruskal post-hoc test (package pgirmess for R). The multivariate Partial Least Squares Discriminant Analysis (PLS-DA) was used to select the metabolites driving the spatial and temporal differences observed. PLS-DA model accuracy was calculated using a two-fold crossvalidation (2CV) and the resulting number of misclassifications (NMC) was compared to its permutated null distribution (999 permutations) to test for model significance (P value < 0.05). Variable Importance in Projection (VIP) were obtained (package RV AideMemoire for R), and features with VIP scores higher than 1 were retained. Further selection of VIPs was performed based on their intensities (threshold of 10.000 at least in two replicates). Kruskal-Wallis test and Kruskal post-hoc test were used to test for significant differences between the different groups (the two different temporal periods detected: 1: April and May and 2: the rest of the months). Putative identification of selected VIPs was performed by searching the exact mass in Metlin (https://metlin.scripps.edu/) with 5 ppm accuracy, comparing with standards present in our laboratory and comparing with published literature on sponge lipids.

Finally, a multiple regression tree (package mvpart in R) was used to study the relative importance of inter-specific variability, temporal and spatial intraspecific variability of the Haliclona metabolomes. MRT allows ranking of factors according to their contribution by producing a tree where at each split data is partitioned into two mutually exclusive groups, each of which is as homogeneous as possible, displaying this way the major factors determining the variability [START_REF] De'ath | Multivariate regression trees: A new technique for modeling species-environment relationships[END_REF]. Only factors with a contribution higher than 0.1% were plotted.

Results

Temporal variability

The metabolome diversity (Shannon index, H) of H. fulva decreased significantly in April and May (P < 0.001) (Fig. 1). A decreased metabolite diversity was also observed for H. mucosa in April (significantly lower diversity, P < 0.001, in April than the rest of months but February and May) and May (significantly lower diversity, P < 0.001, than January, March and September) (Fig. 1a).

PCA showed significant differences in the metabolome composition of both Haliclona species between April and May and the rest of the months (PERMANOVA, P < 0.001) (Fig. 2).

Fig. 2

Principal components analysis (PCA) score plots of H. fulva (a) and H. mucosa (b) on the first two components (PC1 and PC2). Shape represents seasons (winter: Dec, Jan, Feb; spring: Mar, Apr, May; summer: Jun, Jul, Aug; Autumn: Sept, Oct, Dec) and the colors indicate different months Our in-house library of specialized metabolites allowed identification of several features in both sponges based on the comparison between MS spectra and retention times. A total of 10 known metabolites from H. fulva were assigned: renierin 1 modified structure (acid 17-bromo-6-hydroxytricosa-16-ene-2,4dinoic [START_REF] Tribalat | Specialized metabolisms of Mediterranean sponges of genus Haliclona Grant[END_REF])), 18-dihydrorenierin 1, cladocroic acid and fulvynes A-G (found as dicharged ions, see supplementary information). Eight paniceins of our chemical library (A, A hydroquinone, B2, B3, C, D, F1 and F2) were identified in H. mucosa extracts (See supplementary information). All identified metabolites from both sponges significantly decreased in April and May (P < 0.001) (Fig. 3). We were also able to putatively annotate some families of other metabolites based on their isotopic profile and/or fragmentation patterns like lysophospholipids, peptides or fatty acids. However, a clear identification was not sought as our work is mainly focused on the variability of specialized metabolites. The PLS-DA models accurately predicted differences between the metabolomic profiles of April and May and the rest of the months for both sponge species (NMC = 0, P = 0.001; NMC = 0, P = 0.001). After cleaning the matrix (removal of isotopes and adducts), 133 features in H. fulva and 202 in H. mucosa were highlighted with a VIP score higher than 1 and minimum fold-change of 2, and significantly over or under-expressed in April-May (Kruskal-Wallis, P < 0.001). Putative annotation of metabolites using Metlin and the in-house built library allowed the assignment of the families of 75 metabolites in H. fulva and 30 in H. mucosa.

Fig. 3

A total of 38 metabolites were overexpressed in H. fulva in April and May, which included 11 brominated non polar compounds (assigned through their characteristic isotopic pattern), 11 fatty acids and one lysophospholipid (lyso-PAF like molecule). From the 95 metabolites that significantly decreased in April and May in H. fulva, 8 were specialized metabolites identified in our chemical library (6 fulvynes and two renierins), 42 minor doubly charged metabolites with similar masses and retention times than fulvynes (m/z ~ 400 Da, RT ~ 500 s), two fatty acids and four peptides (found as triply charged ions with monoisotopic masses 1600-2000 Da) (Fig. 4A). Over-expressed and under-expressed metabolites in each of the sponge species (a: H. fulva; b: H. mucosa) plotted according to their retention time (x axis) and exact mass (y axis)

In H. mucosa 35 metabolites were overexpressed, 12 of which were putatively annotated as lysophospholipids (lyso-PAF derivatives) and 6 as fatty acids. 167 metabolites significantly decreased in April and May, including 8 paniceins present in our chemical library, one fatty acid and the rest unknowns (Fig. 4b). Six metabolites were overexpressed in both Haliclona species, 2 of which were putatively annotated as fatty acids.

Spatial variability

Metabolomic profiles of H. fulva did not present significant differences between the studied sites (Marseille, Monaco, and Villefranche-sur-Mer) (Fig. 5a). H. mucosa samples from Riou Caramasseigne are significantly different from samples from Villefranche-sur-Mer and Grotte de Jarre (P < 0.001) (Fig. 5a). PCA showed significant differences in the metabolome composition of both Haliclona species sampled at different sites (PERMANOVA, P = 0.002, P .

= 0.02) (Fig. 5b,c). However, these differences were not affected by the geographic distance between the sites. H. fulva samples from Villefranche and Monaco (9 km distance) were not significantly more similar than samples from Grotte de Jarre (167 km away from Villefranche). Similarly, H. mucosa samples from Grotte de Jarre and Riou Caramasseigne (4 km distance) did not show more similar metabolome compositions than Villefranche.

The PLS-DA model was appropriate to predict metabolite differences in H. fulva between different sampling sites (NMC = 0.09, P = 0.001). However, it could not predict accurately the metabolite differences in H. mucosa in different sampling sites (NMC = 0.33, P = 0.032). In consequence, only the VIP from the H. fulva spatial study were obtained. None of the specialized metabolites identified in neither H. fulva nor H. mucosa displayed significant differences among the different sampling sites.

Intra-specific versus inter-specific variability

A MRT analysis was conducted to assess the relative importance of inter-specific variability and the two types of intra-specific variability studied (temporal and spatial) on the variations of the sponges metabolomes. Results show that 97.72% of the differences observed in the sponge metabolomes are explained by the species (H. mucosa and H. fulva). Temporal variability (differences between April and May and the rest of the months) was the second source of variability in both sponge species, accounting for 0.2% in H. mucosa and 0.3% in H. fulva (Fig. 6). Spatial variation accounted for less than 0.1% of the total metabolome variability and therefore it was not plotted in the MRT.

Fig. 6

Multiple regression tree of the H. fulva and H. mucosa metabolomes (both temporal and spatial samples) with the % of variance explained by each split

Discussion

The use of comparative metabolomics in the study of intra-specific variabilities provides numerous advantages compared to classical targeted studies [START_REF] Kuhlisch | Metabolomics in chemical ecology[END_REF]. Metabolic responses of organisms are highly complex and therefore a comprehensive study of the metabolome is required to capture the whole variability and its biologic implications. So far, most studies on metabolomic variabilities of marine invertebrates were performed using targeted approaches [START_REF] Sacristán-Soriano | Temporal trends in the secondary metabolite production of the sponge Aplysina aerophoba[END_REF][START_REF] Duckworth | Effects of water temperature and pH on growth and metabolite biosynthesis of coral reef sponges[END_REF][START_REF] Reverter | Secondary metabolome variability and inducible chemical defenses in the Mediterranean sponge Aplysina cavernicola[END_REF]. Although these studies have the advantage of being quantitative, they explore only a very limited part of the metabolism. Conversely, comparative metabolomics, has the potential to study and uncover previously unidentified metabolites with important biological and ecological roles [START_REF] Kuhlisch | Metabolomics in chemical ecology[END_REF]. This study on the spatio-temporal variability of two Mediterranean Haliclona metabolomes using a comparative LC-MS metabolomics approach showed significant changes in the metabolomic profiles of both species occurring in spring (April and May). As highlighted in previous similar researches on sponge metabolomes (e.g. [START_REF] Ivanisevic | Biochemical trade-offs: Evidence for ecologically linked secondary metabolism of the sponge Oscarella balibaloi[END_REF]), the life cycle may explain at least part of the detected variation. The life cycle of H. fulva has been recently characterised, and a shift from asexual to sexual reproduction has been reported in early April [START_REF] Ereskovsky | Asexual and puzzling sexual reproduction of the Mediterranean sponge Haliclona fulva (Demospongiae): Life cycle and cytological structures[END_REF]. We therefore hypothesize that this reproduction shift induces a change in the metabolite content as reflected for both sponges metabolomes in April and May. The significant decrease in the specialized metabolites of both species is in accordance with the ODT theory, indicating an investment towards the start of the sexual reproduction cycle. For both H. fulva and H. mucosa, this shift occurs only at the early stages of sexual reproduction when the animal needs to trigger the development of sexual gametes, and not throughout all the gametogenesis process, which in H. fulva lasts until September [START_REF] Ereskovsky | Asexual and puzzling sexual reproduction of the Mediterranean sponge Haliclona fulva (Demospongiae): Life cycle and cytological structures[END_REF]. Tradeoff between the production of specialized metabolites and reproduction has been previously reported in the homoscleromorph sponge Oscarella balibaloi, however in this case it was observed at the peak of the reproductive effort [START_REF] Ivanisevic | Biochemical trade-offs: Evidence for ecologically linked secondary metabolism of the sponge Oscarella balibaloi[END_REF]. The increase in fatty acids and lysophospholipids would also be concurrent with the start of reproduction. Fatty acids are important for cell metabolism and are involved in multiple signalling roles in reproduction and early developmental processes [START_REF] Vrablik | Polyunsaturated fatty acid derived signaling in reproduction and development: Insights from Caenorhabditis elegans and Drosophila melanogaster[END_REF][START_REF] Stanley | Eicosanoids in invertebrate signal transduction systems[END_REF]. Physiological roles of lysophospholipids (LPLs) from marine invertebrates remain mostly unknown [START_REF] Sugiura | Distribution of alkyl and alkenyl ether-linked phospholipids and platelet-activating factorlike lipid in various species of invertebrates[END_REF][START_REF] Nishikawa | Cytoprotective effects of lysophospholipids from sea cucumber Holothuria atra[END_REF]. Some studies suggest that LPLs might play a role in invertebrate defence and several phosphatidylcholines and lyso-PAF (Platelet Activating Factor) derivatives isolated from sponges display biological activities such as cytotoxicity, antibacterial, antifungal and antifouling [START_REF] Müller | Molecular/chemical ecology in sponges: Evidence for an adaptive antibacterial response in Suberites domuncula[END_REF][START_REF] Alam | Additional bioactive lyso-PAF congeners from the sponge Spirastrella abata[END_REF][START_REF] Shin | Lyso-PAF analogues and lysophosphatidylcholines from the marine sponge Spirastrella abata as inhibitors of cholesterol biosynthesis[END_REF][START_REF] Lin | Lysophospholipids from the Guangxi sponge Spirastrella purpurea[END_REF][START_REF] Zhao | New lysophosphatidylcholines and monoglycerides from the marine sponge Stelletta sp[END_REF][START_REF] Butler | Antifouling activity oflyso-platelet-activating factor extracted from australian sponge Crella incrustans[END_REF]. Interestingly, [START_REF] Ivanišević | Lysophospholipids in the Mediterranean sponge Oscarella tuberculata: Seasonal variability and putative biological role[END_REF] also observed that level of two bioactive LPLs increased in reproductive sponges suggesting that LPLs could play a role in sponge embryogenesis and morphogenesis.

Sponges are a prolific source of specialized metabolites including rare fatty acid derivatives such as brominated acetylenic acids [START_REF] Pham | Isolation of xestosterol esters of brominated acetylenic fatty acids from the marine sponge Xestospongia testudinaria[END_REF][START_REF] Aoki | Bromotheoynic acid, a brominated acetylenic acid from the marine sponge Theonella swinhoei[END_REF][START_REF] Ohta | Petroacetylene, a new polyacetylene from the marine sponge Petrosiasolida that inhibits blastulation of starfish embryos[END_REF][START_REF] Dembitsky | Lipid compounds of freshwater sponges: Family Spongillidae, class Demospongiae[END_REF]. Several brominated fatty acids were described from Haliclona sponges [START_REF] Alarif Walied | A new cytotoxic brominated acetylenic hydrocarbon from the marine sponge Haliclona sp. with a selective effect against human breast cancer[END_REF][START_REF] Aratake | A new polyunsaturated brominated fatty acid from a Haliclona sponge[END_REF] including renierins from H. fulva [START_REF] Cimino | New acetylenic compounds from the sponge Reniera fulva[END_REF]. In this study, eleven unknown brominated non polar compounds were found overexpressed in H. fulva in April and May, while the concentrations of most of the known renierins identified in our chemical library plummeted. Specialized metabolites can experiment enzymatic bioconversion to best serve the organism's physiological and ecological needs. For example, bioconversion of metabolites into more active ones is well reported in terrestrial and marine organisms following mechanical damage, which is known as activation of chemical defenses [START_REF] Pohnert | Chemical defense strategies of marine organisms[END_REF][START_REF] Paul | Activation of chemical defenses in the tropical green algae Halimeda spp[END_REF]. Bioconversion of metabolites in marine invertebrates related to reproduction has not been reported to date, although this could have passed unnoticed by targeted studies or untargeted approaches which failed to characterize the whole chemical diversity. The decrease of renierins and increase in unknown brominated non polar compounds with similar retention times could indicate the biotransformation of renierins. However, correlation between renierins and the brominated non polar compounds is not well supported. More research is needed to elucidate the structures of these brominated compounds and understand their biological role in the sponge, biosynthesis, and putative contribution to reproduction.

Although H. mucosa and H. fulva possess specialized metabolites from different biosynthethic pathways, we observed similar changes in the composition of the metabolome, characterized by an overall decrease in the concentration of specialized metabolites in April and May for both species. Furthermore, both species experienced an increase in fatty acids and lysophospholipids in this period of the year, and 6 common metabolites were putatively annotated, two of them as fatty acids. These similar variation patterns suggest a similar metabolic response in both species. H. mucosa life cycle has not yet been established, but it will be interesting to study whether there is a reproduction shift in April and May as observed in H. fulva. There is still a debate whether changing environmental conditions such as temperature produces variations in the production of specialized metabolites. Several studies have found higher concentrations of specialized metabolites in summer corresponding to maximum yearly temperatures [START_REF] Abdo | Temperature and spatiotemporal variability of salicylihalamide a in the sponge Haliclona sp[END_REF][START_REF] Ternon | How environmental factors affect the production of guanidine alkaloids by the Mediterranean sponge Crambe crambe[END_REF][START_REF] De Caralt | Sources of secondary metabolite variation in dysidea avara (Porifera: Demospongiae): The importance of having good neighbors[END_REF] while other studies failed to find seasonal changes in production of metabolites [START_REF] Duckworth | Effects of water temperature and pH on growth and metabolite biosynthesis of coral reef sponges[END_REF][START_REF] López-Legentil | Temporal variation in the production of four secondary metabolites in a colonial ascidian[END_REF]. In this study, a clear decrease of specialized metabolites occurred in April and May, when temperature in the Northwestern Mediterranean starts to increase after winter (http://somlit.mio.univ-amu.fr), but no further variation occurred during the warmest period (August-October), indicating the fluctuation in the metabolomes might be due to sponge physiological responses rather than environmental parameters.

Although spatial variabilities at a scale of 200 km for H. fulva and H. mucosa metabolomes were significative, they were less intense than temporal variability.

For the spatial study, different samplings were performed for both sponges at the same time (August 2014) to avoid any temporal bias. Metabolic changes due to spatial variation have been mostly associated to local adaptation to environmental conditions or genotypic selection across habitats [START_REF] Moore | Explaining intraspecific diversity in plant secondary metabolites in an ecological context[END_REF]. However, parameters driving spatial variability and spatial scales at which variation is the highest are still uncertain. In this study, we observed metabolomic differences between sampling sites, which were not dependent on the geographic scale studied (10 km compared to > 100 km). However, the spatial scale studied remains considerably small, and thus larger geographical scales should be studied to assess distance decay in chemical fingerprint similarities [START_REF] Noyer | Relationship between genetic, chemical, and bacterial diversity in the Atlanto-Mediterranean bath sponge Spongia lamella[END_REF][START_REF] Noyer | Patterns of chemical diversity in the Mediterranean sponge Spongia lamella[END_REF]. Other factors that could participate to spatial variations are biotic interactions [START_REF] De Caralt | Sources of secondary metabolite variation in dysidea avara (Porifera: Demospongiae): The importance of having good neighbors[END_REF] and the variable part of the sponge metabolome [START_REF] Borchert | Diversity of natural product biosynthetic genes in the microbiome of the deep sea sponges Inflatella pellicula, Poecillastra compressa, and Stelletta normani[END_REF][START_REF] Schmitt | Assessing the complex sponge microbiota: Core, variable and species-specific bacterial communities in marine sponges[END_REF]. In a recent study on two zoantharian species collected in the tropical and temperate Brazilian Atlantic Ocean, intense metabolomic variations were observed at a large geographical scale [START_REF] Costa-Lotufo | Chemical profiling of two congeneric sea mat corals along the Brazilian coast: Adaptive and functional patterns[END_REF]. However, no clear data were reported for the precise taxonomic identification of the sampled species therefore questioning the conclusions given. Finally, spatial variations in the metabolome of marine invertebrate could result from a complex set of factors linked to environmental parameters of the sampling sites and the ecosystem structure. We thus suggest that temporal variations observed in sponge metabolomes would reflect some physiological changes such as the change in life-stages or stress episodes that affect strongly their metabolism, whilst spatial variations would be the result of a complex set of factors linked to the environment and more difficult to unravel.

Haplosclerids are an important group of sponges in terms of diversity and production of bioactive compounds. However, their taxonomy remains challenging due to major discrepancies between morphological and molecular data [START_REF] Redmond | Phylogenetic relationships of the marine Haplosclerida (Phylum Porifera) employing ribosomal (28S rRNA) and mitochondrial (cox1, nad1) gene sequence data[END_REF][START_REF] Tribalat | Does the chemical diversity of the order Haplosclerida (Phylum Porifera: Class Demospongia) fit with current taxonomic classification?[END_REF]. Molecular data indicates that many very distant related species have been placed in the genus Haliclona, that may even be polyphyletic (organisms that do not share a common ancestor), and re-examination of the systematics of the group Haploscerids is urgently needed [START_REF] Tribalat | Does the chemical diversity of the order Haplosclerida (Phylum Porifera: Class Demospongia) fit with current taxonomic classification?[END_REF][START_REF] Redmond | Phylogenetic relationships of the marine Haplosclerida (Phylum Porifera) employing ribosomal (28S rRNA) and mitochondrial (cox1, nad1) gene sequence data[END_REF]. Integrative taxonomy is a concept that is growing in popularity among taxonomists and consists in a combination of several independent datasets for the classification of species: morphology, anatomy, cytology, spicule shapes, reproduction, molecular data, chemical information etc. [START_REF] Schlick-Steiner | Integrative taxonomy: A multisource approach to exploring biodiversity[END_REF][START_REF] Boury-Esnault | The integrative taxonomic approach applied to porifera: A case study of the homoscleromorpha[END_REF][START_REF] Ivanišević | Lysophospholipids in the Mediterranean sponge Oscarella tuberculata: Seasonal variability and putative biological role[END_REF]. Sponge chemotaxonomy was first proposed by Bergquist and co-workers in the late 60 s [START_REF] Bergquist | Free amino acid patterns and the classification of the demospongiae[END_REF] however its popularity decreased in the last decades due to several pitfalls [START_REF] Erpenbeck | Status and perspective of sponge chemosystematics[END_REF]. Variability of metabolites is considered one of the major pitfalls of chemotaxonomy, since suitable chemotaxonomic markers require stability in their occurrence [START_REF] Erpenbeck | Status and perspective of sponge chemosystematics[END_REF]. In our case and even if we observed very clear patterns of intraspecific variation in both Haliclona species, contribution of intraspecific variability is much lower than interspecific variation, confirming the suitability of comparative and targeted metabolomics as a tool for assisting in the classification of some well groups of sponges. The main limit of its application for the genus Haliclona could be more associated to the correct taxonomic classification of this extremely complex group rather than to wrong insights from metabolomic analyses.

In conclusion, we found similar metabolic variations in both sponge species in April and May, suggesting common physiological responses that could be triggered by changes in environmental factors (e.g. start of reproduction related to thermal cues). This study draws attention to the importance of knowledge on the biological cycle as a key to better assess the associated metabolic changes of marine invertebrates. It also cast lights in production dynamics of specialized metabolites, which is relevant for marine natural product research and drug discovery. Finally, we show for the first time that despite marked intraspecific variations in the metabolomes of both sponge species, these are minor when compared to interspecific variations, indicating the suitability of using comparative metabolomics to help in Haliclona classification.

Fig. 1

 1 Fig. 1Shannon's diversity indices of H. fulva and H. mucosa metabolomes sampled at different months

  Identified specialized metabolites in the metabolomes of H. fulva (a) and H. mucosa (b). Size of bubbles represents the peak area of the feature in the LC-MS chromatogram obtained by XCMS

Fig. 4

 4 Fig. 4

  Fig.5

Acknowledgements

This project (Grant-Aid Agreement No. PBA/MB/16/01) is carried out with the support of the Marine Institute and is funded under the Marine Research Programme by the Irish Government. M.-A.T. received a Ph.D. scholarship from the French Ministry for Higher education and Research. Metabolomic analyses were performed on the MALLABAR platform (Funded by the CNRS, the Provence Alpes Côte d'Azur Region and the Total Foundation). S. Greff (IMBE Marseille, France) is acknowledged for his help in recording and analysing the metabolomic data.

Funding

This project (Grant-Aid Agreement No. PBA/MB/16/01) is carried out with the support of the Marine Institute and is funded under the Marine Research Programme by the Irish Government. The Ph.D. scholarship of M.-A. Tribalat has been funded by the French "Ministère de l'Enseignement supérieur, de la Recherche et de l'Innovation".

Compliance with ethical standards

Author contribution

Methodology and Formal Analysis, M.-A.T., T.P., M.R.; Validation, O.P.T.; Writing-Original Draft Preparation, M.R.; Writing-Review & Editing, O.P.T.; Supervision, T.P., O.P.T.; Project Administration, O.P.T.; Funding Acquisition, T.P., O.P.T.

Conflict of interest

The authors declare that they have no conflict of interest.

8/25/2018

Ethical approval All applicable international, national, and/or institutional guidelines for the care and use of animals were followed. This article does not contain any studies with human participants performed by any of the authors.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary material 1 (DOCX 320 KB)