Tamiris Roque

Jérôme Delettre

Nicolas Hardy

Alvin W Nienow

Frédéric Augier
email: frédéric.augier@ifpen.fr

Fadhel Ben Chaabane

Catherine Béal

The impact of fluid-dynamic stress in stirred tank bioreactors on the synthesis of cellulases by Trichoderma reesei at the intracellular and extracellular levels

Keywords: Trichoderma reesei, cellulases, fluid dynamic stress, scale-up/scale-down, proteomic analysis. Highlights

Cellulases for bioethanol production are mainly made by fed-batch fermentation using a filamentous fungus, Trichoderma reesei. Agitation at different scales impacts on morphology, rheology and growth rate and can be correlated by EDCF εmax . Typically, EDCF εmax is much smaller at commercial scale and fungal size, viscosity and growth rate are greater. Here, to increase understanding, continuous culture in 3 L bioreactors using two EDCF εmax values were undertaken. The higher EDCF εmax decreased the cellulase production (concentration, 21 %; specific production rate, 24 %; protein yield, 20 %) whilst proteomic analysis showed, at an intracellular level, a decrease of cellulase and hemicellulase synthesis. An increase of stress proteins also occurred, which may help cells to limit the impact of fluid dynamic stress. Also, cellulase production during continuous culture at the bench varied with EDCF εmax similarly to that between bench and commercial scale during fed-batch culture.

Introduction

Second-generation (2G) bioethanol from lignocellulosic materials is a good candidate to replace conventional fuels. Indeed, 2G bioethanol allows a reduction of about 85% of greenhouse gas emissions compared to fossil fuels [START_REF] Morales | Life cycle assessment of lignocellulosic bioethanol: Environmental impacts and energy balance[END_REF]. Also, 2G

bioethanol uses lignocellulosic waste from the agri-food and forestry and does not compete with food. However, production costs remains high compared to ethanol from starch. This increase is mainly due to the high price of cellulases, the enzymes which hydrolyze lignocellulosic biomass into simple sugars. The 2G bioethanol price can be reduced if cellulases production costs are lowered [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF]. Economy of scale is a way of doing so, requiring the use of bioreactors of several hundred cubic meters, a huge size compared with that at which the process is developed. With increasing scale, the medium in the bioreactor experiences increasing spatial heterogeneities [START_REF] Amanullah | Mixing in the fermentation and cell culture industries[END_REF][START_REF] Lara | Living with heterogeneities in bioreactors: Understanding the effects of environmental gradients on cells[END_REF] and the fluid dynamic stresses change [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF]; and it is important to know what their impact will be.

Cellulolytic enzymes are mainly produced at industrial scale by the filamentous fungus Trichoderma reesei because of its high cellulases secretion capacity [START_REF] Ferreira | Use of cellulases from Trichoderma reesei in the twenty-first century -Part I: Current industrial uses and future applications in the production of second ethanol generation[END_REF][START_REF] Soni | Cellulases: Role in lignocelullosic biomass utilization[END_REF]. As enzyme production is dissociated from growth, industrial cultures are often carried out in two steps: a batch phase for growth of the fungus followed by a fed-batch phase for the production of cellulases [START_REF] Jourdier | Cellulase activity mapping of Trichoderma reesei cultivated in sugar mixtures under fed-batch conditions[END_REF]. During growth,

T. reesei can use different carbon sources, but glucose is the most common [START_REF] Hardy | Identification des critères d'extrapolation du procédé de production de cellulases par Trichoderma reesei en utilisant l'approche « scale-down[END_REF].

As the T. reesei grows, it develops hyphae which form complex structures (Hardy et al., 2017b). The complex morphology and high concentrations at the end of growth lead to a high viscosity medium with shear thinning behaviour [START_REF] Gabelle | Impact of rheology on the mass transfer coefficient during the growth phase of Trichoderma reesei in stirred bioreactors[END_REF][START_REF] Hardy | Rhéologie des champignons filamenteux: Un outil pour la compréhension d'un procédé de production de biocatalyseurs utilisés pour la production de bioéthanol[END_REF][START_REF] Quintanilla | Fungal morphology in industrial enzyme production-Modelling and monitoring[END_REF].

The production phase starts when the growth stops by limiting the carbon source to induce the production of cellulases [START_REF] Dos Santos Castro | Comparative metabolism of cellulose, sophorose and glucose in Trichoderma reesei using highthroughput genomic and proteomic analyses[END_REF]. In this phase, soluble sugars such as sophorose, cellobiose and lactose are used to induce the secretion of cellulases. In order to avoid further growth, an inductive substrate is fed continuously at an optimal limiting rate [START_REF] Jourdier | Cellulase activity mapping of Trichoderma reesei cultivated in sugar mixtures under fed-batch conditions[END_REF][START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF]). As T. reesei is a strictly aerobic microorganism, the dissolved oxygen concentration (dO 2 ) must be above the critical value, 15 % dO 2 (1.2 mg of O 2 .L -1 ) [START_REF] Marten | Rheological, mass transfer, and mixing characterization of cellulase-producing Trichoderma reesei suspensions[END_REF] for the strain RUT-C30. However, the high viscosity at the end of the batch phase reduces mass transfer rates, and oxygen transfer [START_REF] Albaek | Evaluation of the energy efficiency of enzyme fermentation by mechanistic modeling[END_REF][START_REF] Gabelle | Impact of rheology on the mass transfer coefficient during the growth phase of Trichoderma reesei in stirred bioreactors[END_REF]. To ensure dO 2 > 15 % dO 2 at this time, the power input has to be increased [START_REF] Gabelle | Impact of rheology on the mass transfer coefficient during the growth phase of Trichoderma reesei in stirred bioreactors[END_REF], enhancing the fluid dynamic stress on the mycelia.

There are many ways of characterizing fluid dynamic stress with changes in agitation intensity and with scale [START_REF] Amanullah | Mixing in the fermentation and cell culture industries[END_REF]. Until the 1990s, the two most common were tip speed or specific power input, P/V (W.m -3 ). However, in 1996, it was shown with the fungus, P. chrysogenum, that whilst tip speed increased with scale at constant specific power, damage to the organism was reduced [START_REF] Jüsten | Dependence of mycelial morphology on impeller type and agitation intensity[END_REF]. However, a new function which reduced with increasing scale, the energy dissipation/circulation function, EDCF, was able to correlate fungal breakage and various growth parameters, better than P/V [START_REF] Jüsten | Dependence of mycelial morphology on impeller type and agitation intensity[END_REF]. In essence, the EDCF concept is that fungus break to an equilibrium size dependent on the maximum specific energy dissipation rate, ε max (W.m -3 ) in a region close to the impeller and the frequency with which the organism passes through that [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF] an improved EDCF function, EDCF εmax . We then used EDCF εmax to correlate the impact of the fluid dynamic stress associated with a variety of impellers on various process parameters during the initial batch growth phase of T. reesei fermentations at the bench scale and also with a multiple impeller combination at the industrial scale.

Other filamentous organisms to which the EDCF has been applied are T. harzianum [START_REF] Rocha-Valadez | The influence of circulation frequency on fungal morphology: dissipation rate cultures of Trichoderma harzianum[END_REF], A. oryzae [START_REF] Amanullah | Dependence of morphology on agitation intensity in fed-batch cultures of Aspergillus oryzae and its implications for recombinant protein production[END_REF] and P. ostreatus [START_REF] Fernández-Alejandre | Diffusional and transcriptional mechanisms involved in laccases production by Pleurotus ostreatus CP50[END_REF], where it was related to the transition of the organism from clumps to pellets. With T. reesei [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF], higher values of EDCF εmax reduced the mycelial size and specific growth rate (by up to 20 %) at the bench scale whilst at the commercial scale, as explained above, though the tip speed increased, EDCF εmax decreased leading to larger mycelia, higher viscosity and higher growth rates. However, as the bench runs were limited to the growth phase, cellulases were not produced; so the impact of fluid dynamic stresses on their production could not be investigated.

Previous work on the proteome of T. reesei has been studied in order to better identify the enzymes involved in lignocellulosic biomass degradation and the protein secretion profile under different environmental conditions [START_REF] Herpoël-Gimbert | Comparative secretome analyses of two Trichoderma reesei RUT-C30 and CL847 hypersecretory strains[END_REF].

Apart from the common proteins involved in general metabolism, it mostly consists of two cellobiohydrolases (CBHI and CBH2), five endoglucanases (EGI, EGII, EGIII, EGIV and EGV) and two β-glucosidases (BGLI and BGLII) that act in synergy to degrade lignocellulosic materials [START_REF] Herpoël-Gimbert | Comparative secretome analyses of two Trichoderma reesei RUT-C30 and CL847 hypersecretory strains[END_REF]. Proteomic studies have investigated the secretome and/or intracellular proteins synthesized by filamentous fungi [START_REF] Peterson | Trichoderma reesei RUT-C30 -thirty years of strain improvement[END_REF] by combining one-dimensional (1D) or two-dimensional (2D) gel electrophoresis with protein identification by liquid chromatography/ mass spectrometry. This technique allows the analysis of the multicomponent cocktail of cellulolytic enzymes secreted by T. reesei in culture media [START_REF] Herpoël-Gimbert | Comparative secretome analyses of two Trichoderma reesei RUT-C30 and CL847 hypersecretory strains[END_REF][START_REF] Jun | Insights into enzyme secretion by filamentous fungi: Comparative proteome analysis of Trichoderma reesei grown on different carbon sources[END_REF][START_REF] Kubicek | Systems biological approaches towards understanding cellulase production by Trichoderma reesei[END_REF], and to investigate the differential protein synthesis at intracellular level [START_REF] Bianco | Methodologies and perspectives of proteomics applied to filamentous fungi: From sample preparation to secretome analysis[END_REF]. The influence of carbon source [START_REF] Jun | Insights into enzyme secretion by filamentous fungi: Comparative proteome analysis of Trichoderma reesei grown on different carbon sources[END_REF][START_REF] Dos Santos Castro | Comparative metabolism of cellulose, sophorose and glucose in Trichoderma reesei using highthroughput genomic and proteomic analyses[END_REF][START_REF] Peciulyte | Morphology and enzyme production of Trichoderma reesei Rut C-30 are affected by the physical and structural characteristics of cellulosic substrates[END_REF], pH [START_REF] Adav | Proteomic analysis of pH and strains dependent protein secretion of Trichoderma reesei[END_REF], light intensity [START_REF] Stappler | Abundance of secreted proteins of Trichoderma reesei is regulated by light of different intensities[END_REF] and agitation intensity [START_REF] Mukataka | Variation in cellulase-constituting components from Trichoderma reesei with agitation intensity[END_REF][START_REF] Lejeune | Effect of agitation on growth and enzyme production of Trichoderma reesei in batch fermentation[END_REF] on extracellular enzyme secretion has also been studied. However, except for [START_REF] Jun | Insights into enzyme secretion by filamentous fungi: Comparative proteome analysis of Trichoderma reesei grown on different carbon sources[END_REF], which compared its intracellular proteome with various carbon sources; and Arvas et al., 2011 that studied the intracellular effects of growth rate and cell density in chemostat culture, the impact of culture conditions on intracellular protein synthesis by T. reesei remains poorly investigated.

To our knowledge, the impact of fluid dynamic stress (here expressed as EDCF εmax ) on the synthesis of intracellular proteins by T. reesei has never been reported.

The objective of this work therefore was to characterize the effects of fluid dynamic stress on the production of both intracellular proteins and extracellular cellulases by T. reesei.

Material and Methods

Strain, preculture and culture media

The T. reesei strain used and the composition of the preculture medium and of the batch-culture medium were fully described previously [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF]. For the production phase, the medium had the same composition as the batch-culture medium but it was supplemented with 0.83 g.L -1 yeast extract instead of powdered corn steep liquor.

A 150 mL preculture was prepared in a 2 L Fernbach flask by mixing 135 mL of preculture medium and 15 mL of a sterile glucose solution at 250 g L -1 . The flasks were seeded with 1 mL of conidia to give a seeding rate of 10% and held in a shaker incubator (Multiron II, Infors, Bottmingen, Switzerland) for 72 hours at 30°C and 180 rpm with an orbit of 50 mm before being added into the bioreactor culture.

Bioreactor cultures using two different fluid dynamic stress conditions

This work is a follow on to that reported in [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF]. In that study, the initial batch production of T. reesei was studied at the bench scale in a 3.5 L bioreactor with a working volume of 2.5 L and at commercial scale. At the bench scale, a Rayneri centripetal turbine impeller (VMI-mixing, Montaigu, France) (see Fig 1) was used along with 3 others. All of the changes in process performance brought about by fluid dynamic stress for each impeller were successfully related to changes in EDCF εmax . In order to move to continuous culture, a different bioreactor had to be selected and one of a similar size was available which used a Rayneri impeller. Since that impeller had fitted in well with the others in the earlier work, in order to have bench scale data in a similar sized bioreactor to aid comparison between the different aspects of the overall study, the present configuration was chosen.

A 3.0 L, dished bottomed, stirred bioreactor (ezControl BioBundle, Applikon Biotechnology, Foster, CA, USA) with a working volume of 1.5 L, diameter T = 126 mm with a broth height, H mm such that H/T = 1 and without baffles but with sufficient probes inserted to provide a significant baffling effect. The Rayneri centripetal turbine impeller was 0.08 m diameter, height, 0.045 m with Po = 2.0 and Fl = 1.3. The temperature was controlled at 27 °C, pH at 4.8 and dissolved oxygen concentration at 70 % by introducing a blend of air and nitrogen (thereby altering the driving force for mass transfer) through a sparger. The total flow rate was held constant at 0.8 L/min so that neither the dO 2 or flow rate changed during the experiment [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF]. However, as in earlier work [START_REF] Herpoël-Gimbert | Comparative secretome analyses of two Trichoderma reesei RUT-C30 and CL847 hypersecretory strains[END_REF], to control to this dO 2 at the bench scale even with gas blending required a high P/V [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF] of at least 6 kW. m -3 , which needed a minimum speed of 800 rpm. Such speeds led to very high EDCF εmax compared to the commercial scale, where the circulation time is significantly longer and the specific power was kept as low as possible for economic reasons.

Two complete runs, each consisting of 3 parts, were undertaken. In each, there was an initial batch phase during which the fungus was grown at an agitation speed of 800 rpm until a fungal biomass concentration of  7.9 ± 0.2 g.kg -1 was reached. Then, the continuous phase of the run was begun, thereby commencing cellulase production. That time is marked as 'start of continuous feed' on Fig 1 and from then, it was held at 800 rpm for another 150 hrs which allowed time for approximately 5 volume changes so that a steady state should have been reached [START_REF] Macauley-Patrick | Modes of Fermenter Operation[END_REF]. The concentration of biomass and the production of cellulases also indicated an approximate steady state for 50 hrs. At that point, the speed was increased to 1700 rpm to give a higher fluid dynamic stress and then kept constant with the same feed and discharge rate. This condition was held constant for another approximately 140 hrs, equivalent to approximately 4 volume changes, sufficient to approach another steady state and allow the impact of the change on biomass and intracellular and extracellular proteins to be assessed. The experiment was duplicated and samples were withdrawn at regular intervals for biomass and protein concentration measurements and those taken at the times noted on Fig 1 were used for proteomic analysis.

In the batch phase, the bioreactor contained 0.75 L 2N culture medium, 0.3 L water and 0.3 L glucose solution (150 g.L -1 ) and it was inoculated with 0.15 L preculture to give an initial fungal concentration in the medium of 0.8 g.kg -1 . During the whole of the continuous phase, the constant feed and discharge rates of two solutions were lactose (180 g.L -1 ) at 6 mL.h -1 and the enriched 2N medium at 39 mL.h -1 giving a dilution rate, D (0.03 h -1 ) which since D = , the growth rate, is also 0.03 h -1 .

Determination of concentrations in the culture medium

Biomass was quantified as described previously [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF] as was lactose following the method used for glucose. The concentration of the cellulases produced was quantified on filtered broth samples according to the Lowry method, by using the Protein Assay DC TM kit (Biorad, Hercules, CA, USA) and a range of bovine serum albumin (BSA) concentrations from 0 to 1.5 g.L -1 as standards [START_REF] Hardy | Identification des critères d'extrapolation du procédé de production de cellulases par Trichoderma reesei en utilisant l'approche « scale-down[END_REF]. β-glucosidase activity was measured as described by [START_REF] Jourdier | Cellulase activity mapping of Trichoderma reesei cultivated in sugar mixtures under fed-batch conditions[END_REF]. All analyses were performed in duplicate.

Proteomic analyses

The proteomic analyses of intracellular proteins were carried out in 3 steps: protein extraction, two-dimensional (2D) electrophoresis on sodium dodecyl sulfate polyacrylamide gel (SDS PAGE) and identification of proteins by mass spectrometry coupled to liquid chromatography.

The cytoplasmic proteins were first extracted by cell grinding from a biomass cake obtained by filtration as previously described. Initially, 0.2 g of biomass cake was washed and re-suspended as described previously [START_REF] Wang | Starvation induces physiological changes that act on the cryotolerance of Lactobacillus acidophilus RD758[END_REF]. The cells were lysed with 0.6 g of glass beads (0.3 g of diameter 0.1 mm and 0.3 g of diameter 0.5 mm, Biospec Products, Bartlesville, OK, USA) in a FastPrep FP 120 equipment (Bio 101, Savant Instruments, Holbrook, AZ, USA), for four times with cycles of 30 s at an intensity of 6.5. After grinding the cells, the samples containing solubilized cytoplasmic proteins were centrifuged, purified with 60 U endonuclease (Biorad, Hercules, CA, USA) and precipitated, and their protein concentration was assayed using the colorimetric method of Bradford (1976) using bovine serum albumin as a standard in order to ensure 300 µg of proteins were present [START_REF] Wang | Starvation induces physiological changes that act on the cryotolerance of Lactobacillus acidophilus RD758[END_REF].

The separation of the extracted proteins was achieved by successive separation according to their isoelectric point by isoelectric focusing (IEF) and to their molecular weight by SDS-page electrophoresis. All samples from the two conditions were analyzed in triplicate from three different extractions. The cell extracts containing 300 μg of proteins were diluted in 350 μL of extraction buffer [31]. To perform IEF, the diluted samples were loaded onto Immobiline DryStrip gels (pI 4-7, 17 cm, Biorad) that were submitted to electric focusing at 20°C on an IEF Cell (Biorad, Hercules, CA, USA) to reach 73,000 Wh for about 24 hours.

After equilibration of the DryStrip gels, the second dimension separation was performed in SDS-polyacrylamide gels as described by [START_REF] Wang | Starvation induces physiological changes that act on the cryotolerance of Lactobacillus acidophilus RD758[END_REF]. Proteome image analysis was carried out with PDQuest 2-D software (Biorad) to compare the proteins density on the gels obtained from different conditions. Spots only present on one gel condition or showing a difference in density greater than a factor of 2 were considered as different between two conditions and retained for identification.

The spots corresponding to the differentially synthesized proteins were collected, washed separately with 50 mM ammonium hydrogen carbonate and 50 % acetonitrile and then dried. After digestion with 0.5 µg of trypsin (Promega, Madison, WI) and 50 mM ammonium hydrogen carbonate for 12 h at 37 °C, the supernatant containing peptides were used for protein identification by tandem mass spectrometry coupled to liquid chromatography (LC-MS / MS) using a LTQ-XL mass spectrophotometer (ThermoFinnigan, Waltham, MA, USA) (PAPPSO, Gif-sur-Yvette, France). The sequences were identified using the Joint Genome Institute databases (JGI, Walnut Creek, CA USA).

Results and Discussion

Estimation of the fluid dynamic stress, EDCF εmax ,

In order to estimate 𝜀 (see Section 1), a value of 𝑥 is required. The values [START_REF] Grenville | Mixing: Impeller performance in stirred tanks: Characterizing mixer impellers on the basis of power, flow, shear and efficiency[END_REF] for different impellers range from 12 to 17 but one is not available for the Rayneri impeller, which is a very different shape to others (see Fig 1). However, it produces radial flow as does a Rushton turbine which has an 𝑥 value of 12 and this value is initially assumed here. Thus EDCF εmax was 1884 and 38,400 kW.m -3 .s -1 respectively, a little higher than used previously [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF] because though similar speeds were used, here the bioreactor volume was smaller.

At the commercial scale, multiple impellers were used but the Rushton turbine provided most of the power input [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF]. Thus, both the Po value and x of the Rushton turbine were chosen for the determination of EDCF εmax for the commercial scale, which as explained earlier was very much lower at ~ 6.65 kW.m -3 .s -1 [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF]. Overall, EDCF εmax was the best correlator of hyphal size, broth rheology and growth rate for all the impellers and across these scales.

Effect of fluid dynamic stress, EDCF εmax , on the production of cellulases

The details of the runs are given in Section 2.2. Figure 1 shows the biomass and protein concentration throughout the culture whilst 3.7 ± 0.1 2.9 ± 0.1 0.001 Specific protein (cellulases)production rate, q p (g.kg -1 .h -1 ) 13.5 ± 0.5 10.3 ± 0.8 0.0002 Specific substrate consumption rate, q s (g.kg -1 .h -1 ) 0.084 ± 0.002 0.077 ± 0.002 0.0001 β-Glucosidase activity (IU/mL) 129.0 ± 1.7 107.9 ± 1.2 0.001 𝑌 ⁄ (g.g -1 ) 0.35 ± 0.02 0.36 ± 0.04 0.471 𝑌 ⁄ (g.g -1 ) 0.17 ± 0.01 0.13± 0.01 0.005 𝑌 ⁄ (g.g -1 ) 0.48 ± 0.02 0.37 ± 0.02 0.002 The data presented corresponds to the mean values for each condition (800 and 1700 rpm) for the period shown on Figure 1 and their standard deviations. Legend as in Nomenclature ; na: not applicable.

The specific growth rate remained unchanged since a constant value was imposed by the constant dilution rate, D. However, constant biomass concentration is globally determined as cell dry weight. Ideally, the biomass concentration should also be subdivided into the percentage of active and inactive/damaged cells. However, reliable techniques for quantitatively measuring the viability of filamentous fungi are not available whether by fluorescent dyes, because of sampling issues, or by flow cytometry, because of the filamentous structure of the fungi which tends to block the capillaries. On the other hand, the constant biomass implies that it should all be active.

Table 1 and Figure 1 show that the specific production rate decreased by 24 % and the protein production yields in relation to sugar and to biomass concentration decreased by 20 % and 23 % respectively with p values again showing statistical significance. These results agree with previous studies on the effect of fluid dynamic stress on secondary metabolites production by filamentous fungus in stirred bioreactors. [START_REF] Jüsten | Dependence of Penicillium chrysogenum growth, morphology, vacuolation, and productivity in fed-batch fermentations on impeller type and agitation intensity[END_REF] observed increase in EDCF with different impellers decreased penicillin production during fed-batch fermentations of P. chrysogenum. [START_REF] Lejeune | Effect of agitation on growth and enzyme production of Trichoderma reesei in batch fermentation[END_REF] and [START_REF] Mukataka | Variation in cellulase-constituting components from Trichoderma reesei with agitation intensity[END_REF] found a lower production of cellulases at higher speeds with T. reesei cultures but did not analyze the impact of different impellers or quantify hydrodynamic stress beyond agitator speed. [START_REF] Reese | Shear inactivation of cellulase of Trichoderma reesei[END_REF] determined the effect of 'shear' on the stability of a crude cellulase preparation of T. reesei and found that it caused deactivation of the enzyme exoglucanase cellobiohydrolase (CBH), which slowed down cellulose digestion in a cellulosic biomass hydrolysis. [START_REF] Ganesh | Cellulase deactivation in a stirred reactor[END_REF] observed similar inactivation with increasing speed in stirred reactors. These studies confirm that higher fluid dynamic stresses lead to a decrease in cellulases production by T. reesei and to a deactivation of components of its cellulolytic cocktail.

Effect of higher fluid dynamic stress on the synthesis of intracellular proteins

Scanned pictures of the two SDS-PAGE gels from proteomic analysis of cells are given in Figure 1 (Supplementary Material). The gels indicate intracellular proteins that were differently synthesized at the two EDCF εmax values of 1884 and 38,400 kW.m -3 .s -1 respectively. Image analysis of these gels allowed these proteins to be collected and identified. In total, 30 spots, which corresponded to 24 different proteins, had their synthesis modified as EDCF εmax increased. They were identified by tandem mass spectrometry coupled with liquid chromatography (LC-MS/MS) and compared with in silico data from sequence libraries as described in Section 2. At the higher EDCF εmax , six proteins were newly-synthesized (NW1 to NW6 in Table 2); six proteins were increasingly-synthesized (up-regulated) by times 2 or more (UP1 to UP6, Table 3); and 18 proteins were under-synthesized (down-regulated) by times 2 or more (UN1 to UN18, Table 4). Overall, two groups of proteins can be highlighted: under the higher EDCF εmax , stress proteins were over-synthesized (up-regulated) and cellulases were under-synthesized (down -regulated).

An analysis of

Effect of higher fluid dynamic stress on the synthesis of stress proteins

It is known that environmental stresses cause protein denaturation by aggregation and misfolding leading to loss of biological functions and cell apoptosis [START_REF] Roller | The molecular chaperone GRP78/BiP in the development of chemoresistance: mechanism and possible treatment[END_REF][START_REF] Tiwari | Role of heat-shock proteins in cellular function and in the biology of fungi[END_REF]. Table 2 andFigure 1 (Supplementary Material) show the higher speed (EDCF εmax ) produces the chaperonin protein HSP60 and the ubiquitin-proteasome system which are known to reduce cell damage by acting on post-translational processes [START_REF] Tiwari | Role of heat-shock proteins in cellular function and in the biology of fungi[END_REF], correcting protein folding, helping refolding and stabilizing them under stress [START_REF] Pickart | Ubiquitin and the stress response[END_REF]. Specifically, HSP60 has a crucial role in synthesis, transportation, folding and degradation of proteins [START_REF] Chen | Human 60-kDa heat-shock protein: A danger signal to the innate immune system[END_REF] and its over-expression under heat [START_REF] Raggam | Molecular detection and characterisation of fungal heat shock protein 60[END_REF][START_REF] Tiwari | Role of heat-shock proteins in cellular function and in the biology of fungi[END_REF] and heavy metal [START_REF] Enjalbert | Role of the Hog1 stress-activated protein kinase in the global transcriptional response to stress in the fungal pathogen Candida albicans[END_REF] stress has been reported in fungal pathogens. [START_REF] Kashyap | Deciphering the salinity adaptation mechanism in Penicilliopsis clavariiformis AP, a rare salt tolerant fungus from mangrove[END_REF] observed its expression under saline stress in halotolerant fungus P. clavariiformis. The ubiquitin-proteasome system performs a similar function [START_REF] Shang | Ubiquitin-proteasome pathway and cellular responses to oxidative stress[END_REF][START_REF] Kashyap | Deciphering the salinity adaptation mechanism in Penicilliopsis clavariiformis AP, a rare salt tolerant fungus from mangrove[END_REF] with various fungi; for example, nitrogen deprivation [START_REF] Shang | Ubiquitin-proteasome pathway and cellular responses to oxidative stress[END_REF], heat shock [START_REF] Staszczak | The role of the ubiquitin-proteasome system in the response of the ligninolytic fungus Trametes versicolor to nitrogen deprivation[END_REF] and exposure to cadmium [START_REF] Goller | Trichoderma reesei prs12 encodes a stress and unfoldedprotein-response-inducible regulatory subunit of the fungal 26S proteasome[END_REF].

The present study clearly shows for the first time that the cells response to higher fluid dynamic stress is similar to that with other environmental stresses in order to protect themselves against irreversible damage.

Effect of fluid dynamic stress on the synthesis of cellulases

Figure 1 (Supplementary Material) and Table 4 show that at the higher EDCF εmax , the synthesis of four intracellular cellulases decreased: glucosidases (β-D-glucoside glucohydrolase I), endoglucanases (endo-1,4-β-glucanases), cellobiohydrolases (exo-1,4-βglucanases) and xyloglucanases. The first three make up the T. reesei cellulolytic cocktail [START_REF] Jourdier | Cellulase activity mapping of Trichoderma reesei cultivated in sugar mixtures under fed-batch conditions[END_REF]. Cellobiohydrolases attack the cellulose chains from the ends, releasing cellobiose molecules. Endoglucanases hydrolyze cellulose chains randomly, generate new ends accessible to cellobiohydrolases. β-glucosidase has the role of hydrolyzing into glucose the cellobiose released by the action of endoglucanases and cellobiohydrolases. This synergistic action is essential for hydrolysis of lignocellulosic biomass by the fungus [START_REF] Kumar | Bioconversion of lignocellulosic biomass: biochemical and molecular perspectives[END_REF] as are xyloglucanases [START_REF] Herpoël-Gimbert | Comparative secretome analyses of two Trichoderma reesei RUT-C30 and CL847 hypersecretory strains[END_REF] for improving hydrolysis efficiency as the hemicellulolytic activity of this enzyme increases the surface area [START_REF] Benkő | Evaluation of the role of xyloglucanase in the enzymatic hydrolysis of lignocellulosic substrates[END_REF].

The fact that the synthesis of cytoplasmic cellulases is negatively affected by higher fluid dynamic stress is strongly linked to the reduction of extracellular cellulases (Figure 1 and Table 1). This result indicates that the profile of extracellular proteins is directly linked to the synthesis of intracellular proteins, as demonstrated above from the proteomic analysis. This finding is original because, although some authors have reported the negative effects of 'shear' stress on the secretion of cellulases by T. reesei [START_REF] Mukataka | Variation in cellulase-constituting components from Trichoderma reesei with agitation intensity[END_REF][START_REF] Lejeune | Effect of agitation on growth and enzyme production of Trichoderma reesei in batch fermentation[END_REF], the literature has not made any connection with the effect of fluid dynamic stress on the intracellular synthesis of cellulases.

Effect of scale on the production of cellulases

For this study in chemostat culture, EDCF εmax was 1884 kW.m -3 .s -1 and 38400 kW.m - 3 .s -1 and the proteomic analysis showed that, as a result, a change took place at the intracellular level to protect the cells and partially maintain the production of cellulases. As a result of this protective action, though EDCF εmax increased substantially, the fall in cellulase production was relatively small, from 13.5 to 10.3 g.kg -1 . h -1 .

In the batch cultures previously reported [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF] (i.e. during the growth phase of the fermentation), EDCF εmax was a little smaller at the bench scale than here because a larger bench scale bioreactor was used. As usual, it was also much lower than at the commercial scale (6.65 kW.m -3 .s -1 ) for the reasons given earlier [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF] and in Section 3.1. The lower EDCF εmax was mirrored by increases in the growth rate, the viscosity of the broth and the size of the fungi.

At the commercial scale [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF], since the aim was to produce as much cellulase as possible, fed batch culture was used for the production phase to give a q p of 17 g.kg -1 .h -1 . In bench scale work [START_REF] Hardy | Identification des critères d'extrapolation du procédé de production de cellulases par Trichoderma reesei en utilisant l'approche « scale-down[END_REF] related to this continuous culture study and the earlier batch study [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF], two fed batch cultures were also undertaken. These fed-batch studies were done in the same 3.5 L bioreactor as used by [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF] with the Rayneri centripetal impeller at 1000 rpm to give an EDCF εmax of 2760 kW.m -3 .s -1 . The biomass concentration was controlled to 3.5 g/kg or 12 g/kg but the impact on the production of cellulases was not statistically significant (q p = 11.9 ± 3.3 g.kg -1 .h -1 ). Clearly, the scatter was large [START_REF] Hardy | Identification des critères d'extrapolation du procédé de production de cellulases par Trichoderma reesei en utilisant l'approche « scale-down[END_REF]. It is interesting to compare the small change in q p for cellulase production with the very different EDCF εmax values in the continuous culture conditions reported here and those arising in fed-batch at the two very different scales. It is also convenient to assess whether the choice of x = 12 for the Rayneri impeller is an appropriate one. [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF] showed that in batch culture, the growth rate µ is linked with EDCF εmax and under the conditions studied, the growth rate was well described by the empirical function, Eq 1 for all the impellers tested at the bench scale and for the data obtained at the commercial scale:

µ 𝐴 • 𝑙𝑛 𝐸𝐷𝐶𝐹 𝐵 ( 1 )
In order to investigate the effect of EDCF εmax on the specific production rate q p , a similar analysis is proposed. However, given the low number of experiments performed, it is not possible to validate the function proposed by [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF] or indeed any other function. Given this situation, it has been assumed as a working hypothesis, that a similar function of EDCF εmax is usable with other constants A and B. Thus,

𝑞 𝐴 • 𝑙𝑛 𝐸𝐷𝐶𝐹 𝐵 ( 2 )
Based on this hypothesis, the results of the present work are used to fit constants to Eq 2 and then used for scale-up consideration. Using the q p values obtained in the present work at the two EDCF εmax gives A = -1.05 and B = 21.4. Eq 2 can now applied to the fed-batch case at the bench (subscript fb1) and commercial scale (subscript fb2) to give the empirical Eq 3:

𝑞 𝑞 𝐴 • 𝑙𝑛 𝐸𝐷𝐶𝐹 𝐸𝐷𝐶𝐹 (3)
If it is assumed that though different fermentation methods (fed-batch and continuous) give different product yields, the impact of fluid dynamic stress is the same for the particular organism, then the specific production rate at commercial scale can be obtained by Eq 4: 𝑞 11.9 1.05 • 𝑙𝑛 2760 6.65 ⁄ 18.3

This predicted value of q p (18.3 g.kg -1 .h -1 ) at the commercial scale is similar to the actual one (17 g.kg -1 .h -1 ), especially considering the large difference in the scales used for these two fed-batch runs and that other factors can come into play at the commercial scale such as increased levels of inhomogeneity due to the greater mixing time and pCO 2 and dO 2 due to the static head. Nevertheless, it is clear that more data are required to confirm this relationship.

Table 5 summarizes the different culture conditions, scales, process parameters, EDCF εmax and results obtained in the whole of this wide ranging study. x value of 17 for the bench scale fed batch run, EDCF εmax becomes 3910 kW.m -3 .s -1 . If this value is used in Eq 3 then, since EDCF εmax for the large scale where the Rushton turbine dominates remains the same, the predicted value of q p for the commercial scale becomes 18.6 g.kg -1 .h -1 . In this case, the choice of x clearly is not a significant factor in the analysis of the impact of fluid dynamic stress on the production of cellulases. Therefore given the similarity in flow pattern between the Rayneri impeller and the Rushton turbine, x = 12 remains the value of choice.

Overall, it can be concluded that not only can the use of the EDCF εmax function enable bench scale data to predict hyphal size, broth rheology and growth rate at the commercial scale [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF], it also seems to be a promising tool in aiding scale-up in relation to predicting the production of cellulases. Indeed, the production of cellulases at the commercial scale should not be compromised by fluid dynamic stresses. More importantly, the agitation intensity chosen should be based on ensuring it is sufficient to meet the demands of oxygen mass transfer and carbon dioxide stripping along with an adequate blending strategy, including sub-surface addition of nutrients and pH control chemicals [START_REF] Amanullah | Mixing in the fermentation and cell culture industries[END_REF]). In addition, economic considerations must be taken into account.

Conclusion

T. reesei was cultivated in continuous cultures at two EDCF εmax values, the lower one (1884 kW.m -3 .s -1 ) until a steady state was reached and the higher (38400 kW.m -3 .s -1 ) until a near-steady state, sufficient to show a significant fall in protein production linked to measureable changes in intracellular proteins. This method of defining fluid dynamic stress has already been successful in correlating its impact in the batch growth phase of T. Reesei at the bench and commercial scale [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF]. Here, the use of EDCF εmax provided a link between this bench scale continuous culture study with that batch work [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF] and also related fed-batch culture to produce cellulases at the bench [START_REF] Hardy | Identification des critères d'extrapolation du procédé de production de cellulases par Trichoderma reesei en utilisant l'approche « scale-down[END_REF] and commercial scale ( 100 m 3 ) [START_REF] Hardy | Scale-up agitation criteria for Trichoderma reesei fermentation[END_REF]. It is worthy of note, that a recent extensive review of scale-up methods concluded that the EDCF function was the best way of defining agitation intensity in bioprocessing [START_REF] Böhm | Multiphase stirred tank bioreactors -New geometrical concepts and scale-up approaches[END_REF].

At the higher EDCF εmax at the bench scale, the concentration of the cellulases in the medium as well as the production yield and the specific production rate decreased by 21 %, 20 % and 24 % respectively. A proteomic analysis provided further insight that helped to explain these observations. Firstly, at the higher EDCF εmax , the intracellular synthesis of all enzymes involved in cellulosic activity was decreased. On the other hand, stress proteins such as the HSP60 chaperone protein and the ubiquitin-proteasome system whose role is to protect the cells from the effects of fluid dynamic stress appeared. These findings demonstrate that at higher fluid dynamic stress, the fungus adapts itself by triggering protection and damage repair mechanisms at the intracellular level and by favoring its central metabolism to the detriment of other less essential functions, the combination of the two explaining the small decrease of cellulase biosynthesis even at much higher EDCF εmax .

The higher extracellular cellulase production (q p ) at the lower EDCF εmax in continuous culture compared to that at the higher EDCF εmax is quantitatively similar to the increase in q p during fed-batch fermentation at the commercial scale compared to that at the bench scale.

Perhaps more importantly for the commercial aim of producing cellulases, the sensitivity to fluid dynamic stress is relatively small as the proteomic analysis explains.

  for the time periods indicated by the dotted lines on Fig 1. The variables describing the biomass (biomass concentration, specific growth rate, growth yield) were not influenced by the agitation rate (p-values higher than 0.05). The variables that characterize the protein production (protein concentration, specific production rate, production and conversion yields and β-glucosidase activity) were significantly affected by the EDCF εmax value (p-values lower than 0.05).

Figure 1 :

 1 Figure 1: Biomass and protein concentrations during continuous culture agitated by a Rayneri centripetal turbine impeller (embedded image) at initially 800 rpm increasing to 1700 rpm (EDCF εmax 1884 and 38,400 kW.m -3 .s -1 respectively). Dotted lines indicate time periods used for ANOVA tests. Error bars represent the standard deviation.

  , EDCF generally decreases with scale, sometimes markedly. There are a number of ways of assessing ε max and a recently developed definition by Grenville et al. (2017) is 𝜀 1.04 . 𝑥 . 𝜌 . 𝑃𝑜 . 𝑁 . 𝐷 . Here, x is the ratio of impeller diameter to trailing vortex diameter, a dimensionless characteristic of the impeller being used, as is Po, the power number; and 𝜌 (kg.m -3 ) is the density of the broth. In 2017, we used this definition of

	𝜀	to develop
				region, 1/t c
	(s) where the circulation time, 𝑡	𝑉 𝐹𝑙 • 𝑁 • 𝐷 ⁄	. Here, 𝑉 (m 3 ) is the volume of broth in

the bioreactor, 𝑁 (rev.s -1 ) the agitator speed, 𝐷 (m) its diameter and 𝐹𝑙 (dimensionless), its flow number. Because the circulation time increases with increasing scale even at constant 𝑃 𝑉 ⁄

  Table 1 gives the mean biomass and cellulase concentrations at the first steady state (the lower EDCF εmax ) and the near-steady state at the higher EDCF εmax . Table 1 also gives the value of a range of other process parameters obtained for these two conditions and p values from ANOVA tests for the results

Table 1 :

 1 Measured values of culture parameters and associated p-values at the two different speeds and EDCF εmax values

	Stirring speed	800 rpm	1700 rpm	p-value
	𝐸𝐷𝐶𝐹	(kW.m -3 .s -1 )	1884	38400	na
	Biomass concentration, X (g.kg -1 )	7.7 ± 0.2	7.6 ± 0.2	0.538
	Specific growth rate, µ (h -1 )	0.028 ± 0.01 0.028 ± 0.03	0.064
	Protein concentration, P (g.kg -1 )

  Table 2 shows the increase in EDCF εmax only produced five new proteins, as two of the 6 spots are isoforms of the same protein. Of the others, three assist

	in protein quality control (including one chaperonin and two ubiquitin dependent
	proteasomes) and are recognized as stress proteins. From Table 3, three up-synthesized
	proteins are involved in the central carbon metabolism as they are related to
	oxidoreductase activity (dihydrolipoyl dehydrogenase and FAD/NAD(P) domain-containing

protein) and Acetyl CoA biosynthesis (acetate kinase). The four others are involved in carbon metabolism (adenosylhomocysteinase), proteolysis (glutamate carboxylpeptidase) and GDPmannose biosynthesis (phosphomannomutase). Finally, of the under-synthetized proteins (Table

4

), there are a number of isoforms and only 12 different ones are identified from 18 spots of which four relate to cellulolytic (glucosidases, endoglucanase and cellobiohydrolases) and hemicellulolytic (xyloglucanases) activities. Of the eight others, one involved putrescine biosynthesis (putative agmatine deiminase), two to transferase activity (transketolase and uracil phosphoribosyl transferase) plus one isomerase (UDP-glucose-4epimerase), one kinase (diphosphomevalonate decarboxylase), one oxidoreductase (NAD(P)H-dependent-D-xylose reductase), one involved in Ca 2+ homeostasis (SGL domaincontaining protein) and one implicated in proteolysis (Zn-dependent exopeptidases).

Table 2 :

 2 Proteins newly-synthesized by the increase from 800 to 1700 rpm (EDCF εmax 1884 and 38,400 kW.m -3 .s -1 respectively)

	Spot	Protein ID* Protein description	MW (KDa) in the 2D gel	IP in the 2D gel	Coverage** (%)	Functional category
	NW1	135423	Chaperonin HSP60	63.1	5.25	75	Chaperone molecular family
	NW2	26797	20S proteasome subunit alpha type (homologs PSMA4 / PRE9)	26.2	6.12	73	Threonine endopeptidase
							activity -Ubiquitin-
	NW3	99350	20S proteasome alpha type (homologs PSMA6 / SCL1)	26.3	6.30	73	dependent protein
	NW4	80225	Uroporphyrinogen decarboxylase	39.2	5.99	71	Porphyrin biosynthesis
	NW5 & NW6***	39524	Malate dehydrogenase	35.5 35.5	5.92 6.22	37 49	Oxidoreductase activity -Carbohydrate metabolism

MW: molecular weight; IP: isoelectric pH; *Protein ID were obtained from JGI database; **Coverage corresponds to the percentage of the protein sequence covered by identified peptides; ***Spots that correspond to isoforms of the same protein

Table 3 :

 3 Proteins over-synthesized from the increase from 800to 1700 rpm (EDCF εmax 1884 and 38,400 kW.m -3 .s -1 respectively)

	Spot	Protein ID*	Protein description	MW (KDa) in the 2D gel	IP in the 2D gel	Coverage** (%)	Log2 fold change	Functional category
	UP1	78299	Dihydrolipoyl dehydrogenase	55.3	6.47	77	1.9	Oxidoreductase activity -Cell
	UP2	25190	FAD/NAD(P)-binding domain-containing protein	59.5	6.07	78	1.4	redox homeostasis
	UP3	97186	Phosphomannomutase	28.1	5.38	71	1.3	GDP-mannose biosynthetic process
	UP4	137356	Glutamate carboxypeptidase	57.8	5.69	82	2.5	Proteolysis
	UP5	99242	Acetate kinase	46.6	6.49	81	1.3	Acetyl-CoA biosynthesis
	UP6	142425	Adenosylhomocysteinase	46.0	6.22	68	1.5	Nucleic acid and amino-acid metabolism

MW: molecular weight; IP: isoelectric pH; *Protein ID were obtained from JGI database; **Coverage corresponds to the percentage of the protein sequence covered by identified peptides

Table 4 :

 4 Proteins under-synthesized by the increase from 800 rpm to 1700 rpm (EDCF εmax 1884 and 38,400 kW.m -3 .s -1 respectively)

	345
	Spot	Protein ID* Protein description	MW (KDa) in the 2D gel	IP in the 2D gel	Coverage** (%)	Log2 fold change	Functional category
	UN1	122470	Exoglucanase II (1,4-β-cellobiohydrolase)	63.5	5.22	67	-2.3
	UN2	124438	Endo-β -1,4-glucanase	22.8	5.81	71	-2.3
				94.8	5.98	66	-1.5
				91.8	6.15	65	-1.5	Hydrolase activity -Cellulases
	UN3 to UN7***	136547	β-D-glucoside glucohydrolase I	74.9	6.15	67	-2.8
				74.5	6.34	69	-2.6
				75.7	6.49	71	-2.8
	UN8 & UN9***	111943	Xyloglucanase	96.5 92.1	4.97 4.87	61 63	-2.4 -1.5	Hydrolase activity -Hemicellulase
	UN10	108605	Uracil phosphoribosyl transferase	25.4	5.82	85	-1.2	Transferase activity -Nucleoside metabolism
	UN11	111063	Peptide hydrolase	67.5	4.60	58	-2.3	Proteolysis
	UN12 & UN13***	102903	Putative agmatine deiminase	47.0 46.8	4.55 4.65	84 40	-1.5 -2.6	Putrescine biosynthesis

MW: molecular weight; IP: isoelectric pH; *Protein ID were obtained from JGI database; ** Coverage corresponds to the percentage of the protein sequence 346 covered by identified peptides; ***Spots that correspond to isoforms of the same protein 347

Table 5 .

 5 Comparison of the three different culture conditions and results obtained in this overall wide ranging study.The question regarding the impact of the value of x (the ratio of impeller diameter to trailing vortex diameter) chosen for the Rayneri impeller can also be addressed in relation to this type of analysis.[START_REF] Grenville | Mixing: Impeller performance in stirred tanks: Characterizing mixer impellers on the basis of power, flow, shear and efficiency[END_REF] give a range of x values from 12 to 17 with 12 being the value for the Rushton turbine and chosen for the Rayneri because, though very different from most impellers, they are both radial flow. If the other extreme, x = 17, is chosen to test the range of possible impacts on this analysis, then the EDCF εmax values of the continuous cultures reported in Table1are multiplied by 17/12, the B parameter of Eq.1 becomes 21.9 instead of 21.4, and the A parameter remains unchanged. Applying the same

	Source	Culture Conditions	Bioreact or Scale	Process Parameters	EDCFεmax (kW.m -3 .s -1 )		Results
	Hardy (2016) Fed batch	3.5 L and 100 m 3	𝑞	2760 and 6.65	𝑞 : 0.012 and 0.017 g.kg -1 .h -1
				Rheology,
	Hardy et al., (2017)	Batch	3.5 L and 100 m 3	biomass and morphology	24000 to 6.6	𝑌	𝐴 • 𝑙𝑛 𝐸𝐷𝐶𝐹	𝐵
				Y
	This work	Continuous	3.0 L	𝑞	38400 and 1880	𝑞 : 0.010 and 0.014 g.kg -1 .h -1

Acknowledgements

The authors acknowledge the Plateforme d'Analyse Protéomique de Paris Sud Ouest (PAPPSO, Gif-sur-Yvette, France) for help in identifying the proteins. Funding: This work was supported by the Agence de l'Environnement et de la Maîtrise de l'Energie (ADEME, Angers, France) by providing financial support for the PhD study of Nicolas Hardy [number 2016SACLA017].

Nomenclature

Supplementary material

Figure 1 : Proteomic analysis of T. reesei cells: A) 800 rpm (EDCF εmax 1884 kW.m -3 .s -1 ); B) 1700 rpm (EDCF εmax 38400 kW.m -3 .s -1 ); UN1-UN19: under-synthesized proteins; UP1-UP7: over-synthesized proteins; NW1-NW6: newly-synthesized proteins.


Sylvain Caurla
email: sylvain.caurla@inra.fr

Antonello Lobianco

Estimating climate service value in forestry: The case of climate information on drought for maritime pine in Southwestern France

Keywords:

In Western Europe, future climate changes go hand-in-hand with increasing risks of droughts and heat waves during summer. For forest ecosystems, a drought may result in both an increase in tree mortality and a reduction in tree growth. These impacts are delayed over time, i.e., there is a time gap between the drought and its impacts on the forest stand, which makes it possible to adjust forest management practices and, in particular, to prematurely harvest the impacted stand and to replant a new one if it is economically profitable to do so. Consequently, we define Climate Services (CS) as the information that supports forest owners in their decision to prematurely harvest or not after a drought. Our paper aims at developing a method to estimate the economic value of these CS in the case of a maritime pine stand in Southwestern France. Using a comparison of Land Expected Values (LEV) over an infinite period of forest rotations, our analysis suggests that the Climate Service Value (CSV) is highly dependent on three characteristics: (1) the age of the forest exposed to the drought; (2) the intensity of the drought in terms of both mortality and growth impacts; and (3) the discount rate value used. Overall, for a 2% discount rate for a young stand (less than 15-20 years old), the CSV is rather low and ranges from 0 to €50/ha depending on the intensity of the drought. However, for a mature stand, the CSV rapidly increases, up to a maximum the year before the optimal harvest date. In this latter case, the CSV may be as high as €4900/ha for intense droughts.

Overall, for a 2% discount rate, we show that, for a young stand (less than 15-20 years old), the CSV is rather low and ranges from 0 to €50/ha depending on the intensity of the drought.

Practical implications

A review of drought and heat-induced tree mortality by [START_REF] Allen | A global overview of drought and heat-induced tree mortality reveals emerging climate change risks for forests[END_REF] suggests that no forest type or climate zone around the world is invulnerable to such risks, even zones that are not considered as water-limited. However, forest owners cannot easily predict the time and the intensity of a drought since it occurs at random points in time and causes random-sized damages. For forest ecosystems, a drought results in both an increase in tree mortality and a reduction in tree growth. These impacts are delayed over time, usually begin to appear the year after the drought, and last for several years. The time gap between the observation of the climate event and its consequences makes it possible to adjust forest management practices and, in particular, to prematurely harvest the impacted stand and to replant a new stand if it is economically profitable to do so. In this study, we define a Climate Service (CS) as the information about the intensity of a drought event and its impacts on tree growth and mortality provided to forest owners just after a drought in order to help them decide to prematurely harvest (or not). For example, this CS could take the form of a combination of an eco-physiological model with an economic model linking climate and forest dynamics with economic outputs.

Estimating the economic value of this type of CS provides public institutions and decision-makers with an instrument to calibrate their investments in anticipation of a CS. In addition, as more and more private CS providers emerge, comparing the theoretical CS value with the willingness-to-pay of CS consumers makes it possible to better design the services exchanged on this market.

Using a comparison of Land Expected Values over an infinite period of forest rotations, our analysis suggests that the Climate Service Value (CSV) for a maritime pine stand undergoing a drought event is highly dependent on three characteristics:

(1) the age of the forest subject to the drought;

(2) the intensity of the drought both in terms of mortality and growth impacts;

(3) the value of the discount rate used.

However, for a mature stand, the CSV rapidly increases, up to a maximum the year before the optimal harvest date. In this latter case, the CSV can reach a maximum of €4900/ha for an intense drought.

These results suggest that a potential market for CS provision in the forest sector exists. In the case of drought, the main challenge for CS providers is to better assess the impacts of drought in terms of additional mortality and tree growth reduction. To tackle this challenge, mechanistic eco-physiological models are being developed [START_REF] Bréda | Temperate forest trees and stands under severe drought: a review of ecophysiological responses, adaptation processes and long-term consequences[END_REF][START_REF] Davi | Assessing drought-driven mortality trees with physiological process-based models[END_REF]. One potentially important challenge for CS providers in the forest sector is to combine information on the impacts of climate events with economic information that can support forest managers in their decisions.

Introduction

In Western Europe, even under conservative scenarios, future climate changes are likely to come with larger winter rains, an increasing risk of winter storms and more severe precipitation deficits and heat waves during summer [START_REF] Collins | Long-term climate change: projections, commitments and irreversibility[END_REF]. Extreme heat waves and drought episodes like those experienced during summer 2003 in France, Germany and Spain are therefore expected to occur at increasing frequency. For forest ecosystems, a drought may result in both an increase in tree mortality and a reduction in tree growth. First, by limiting water availability, drought may induce large-scale tree growth decline [START_REF] Bréda | Temperate forest trees and stands under severe drought: a review of ecophysiological responses, adaptation processes and long-term consequences[END_REF]. Second, at some point, when the intensity of drought increases, this may result in the premature mortality of roots or twigs, and could ultimately lead to tree death [START_REF] Battaglia | Prediction of leaf area index in eucalypt plantations: effect of water stress and temperature[END_REF][START_REF] Dantec | Interannual and spatial variation in maximum leaf area index of temperate deciduous stands[END_REF]. The intensity of these impacts depends on the characteristics of the drought (duration, quantity and temporal distribution of precipitation), as well as to the combination of the drought with other climatic parameters (sunshine duration, relative humidity, mean and extreme temperatures) or with other climate-related disturbances such as pathogen invasion [START_REF] Bréda | Temperate forest trees and stands under severe drought: a review of ecophysiological responses, adaptation processes and long-term consequences[END_REF][START_REF] Caminero | Growth responses to climate and drought at the southernmost European limit of Mediterranean Pinus pinaster forests[END_REF][START_REF] Mina | Forward modeling of treering width improves simulation of forest growth responses to drought[END_REF].

Importantly, these impacts are not immediate but are delayed over time, with an increase in mortality observed during the very first years after the drought [START_REF] Renaud | Results from the European Network For Damage Monitoring, Département de la santé des forêts[END_REF][START_REF] Vennetier | Impact of climate change on pines forest productivity and on the shift of a bioclimatic limit in Mediterranean area[END_REF], whereas tree ring width and leaf area are frequently smaller for several years following a severe drought [START_REF] Caminero | Growth responses to climate and drought at the southernmost European limit of Mediterranean Pinus pinaster forests[END_REF][START_REF] Mina | Forward modeling of treering width improves simulation of forest growth responses to drought[END_REF]. This time gap between the observation of the climate event and its consequences makes it possible to adjust forest management practices and, in particular, to prematurely harvest the impacted stand and to replant a new stand if it is economically profitable to do so.

In forestry, climate services (CS) can support this kind of decision by taking the form of recommendations concerning initial plantation densities, choices of new species or species mixes, changes in rotation lengths or landscape planning for the purpose of minimizing fire and insect damage [START_REF] Easterling | Food, fibre and forest products[END_REF]. In particular, the existence and the quality of information to support the decision to change the rotation length in order to reduce economic losses are precious for forest managers. However, the economic value of this type of CS is rarely estimated and, to our knowledge, there is no example for forest activities in the literature. Yet, estimating this value could provide public institutions and decision-makers with an instrument that would allow them to better calibrate their investments in the supply of CS for the forest sector. Meanwhile, as more and more private CS providers emerge, comparing the theoretical CS value with the willingness-to-pay of CS consumers makes it possible to better design the services exchanged on this market.

In this study, we propose to estimate the value of CS in the case of a drought for a monospecific maritime pine forest in the Landes region in Southwestern France, as shown on Fig. 1. In this case, the CS consists of the information provided to forest owners about the timing and the intensity of the drought, both in terms of increases in mortality and growth decline. We compare two scenarios: one where CS are available and the other where they are not. The numerical results suggest that for an average 2% discount rate, the CS ranges between €0-50/ha for a young forest stand and between €0-4900/ha for a mature forest stand close to its optimal year of harvesting.

The article is structured as follows. In Section 1, we provide a literature review of climate services in the agricultural sector. Since no study dealing with CS value assessment in forestry was found, this literature review investigates the possibility of adapting methods developed in agriculture to the forest sector. We then describe the model developed in our framework in Section 2. In Section 3, we present the results and the analysis, discuss them in Section 4, and provide a conclusion in Section 5.

Review: the economic value of climate services in agriculture and forestry

In the agriculture sector, CS may take the form of raw climate observations and projections, climate indices (days of frost, length of a heatwave), and information such as mean temperatures [START_REF] Vaughan | Climate services for society: origins, institutional arrangements, and design elements for an evaluation framework[END_REF]. They can also encompass adaptation solutions or advice regarding the best practices to adopt [START_REF] Street | Towards a leading role on climate services in Europe: a research and innovation roadmap[END_REF]. The economic literature extensively explores the concept of the expected value of information. Assuming that the information is perfect (i.e., predicted climatic conditions before decision-making will occur), the expected value of the information and, therefore, the value of CS are defined as the difference between the expected value of an optimal action with information and the expected value of an optimal action without information. In this vein, [START_REF] Mjelde | Valuing forecast characteristics in a dynamic agricultural production system[END_REF] focused on the value of seasonal climate forecasts in corn production, defined as the difference between the expected net returns when the information is used optimally and when no additional information is used. They showed that a less accurate forecast received earlier in the production process might be more valuable than a more accurate forecast received later. However, as [START_REF] Runge | Which uncertainty? Using expert elicitation and expected value of information to design an adaptive program[END_REF] point out, information is rarely perfect. Several sources of uncertainty may arise and the relative importance of the different sources of uncertainty can be tested in order to decide which factor is more important to focus on. [START_REF] Hilton | The determinants of information value: synthesizing some general results[END_REF] proposed to take these uncertainties into account by defining the expected value of the information as the difference between the expected utility of the agent using the weather forecast distribution and the expected utility of the agent according to the climatic conditions of the past. However, this method assumes that users have complete knowledge of historical climate conditions, which is highly unlikely.

To assess the economic consequences of climate arising from various El Niño-Southern Oscillation (ENSO) phases, [START_REF] Adams | The benefits to Mexican agriculture of an El Niño-southern oscillation (ENSO) early warning system[END_REF] developed a stochastic price-endogenous economic model for Mexico that accounts for changes in cropping patterns, production and consumption arising from the yield changes under each ENSO phase forecast. They show that the benefits of an ENSO early warning system for Mexico are approximately US$ 10 million annually, based on a 51-year time period of ENSO frequencies when a forecast skill of 70% is assumed.

Similarly, [START_REF] Roudier | An ex-ante evaluation of the use of seasonal climate forecasts for millet growers in SW Niger[END_REF] computed the economic value of seasonal climate forecasts for millet producers in Nigeria according to the simulated yields derived from the different management strategies under different types of forecasts (imperfect, perfect and perfect with rainy season onset and offset dates). More recently, [START_REF] Lechthaler | The climate challenge for agriculture and the value of climate services: application to coffee-farming in Peru[END_REF] estimated the value of CS for the coffee sector in Peru using a stochastic life-cycle model of a rural household that faces uncertainty with respect to the timing and the size of an adverse weather shock. In their paper, the one-year time horizon is divided into two sub-periods. In the first period, farmers receive an income that is either consumed or saved for precautionary measures against future coffee rust. In the second period, the savings are invested in the precautionary measure as the climate event occurs. The farmer's profit is then impacted according to the timing of the occurrence of the weather shock, its intensity and the amount of the savings. The value of climate service is defined as the error made by the decision-maker between his discounted welfare with climate services (no uncertainty) and the one expected under uncertainty.

One missing feature of these studies is the treatment of multiannual systems such as forests. One important strategy to adapt to climate change in multiannual systems is to modify the rotation length. For example, by harvesting the trees just before the consequences of a climate event occur, forest owners can reduce their economic losses. For a climate event such as a storm, the event and its impacts on the forest stand are concomitant. In that case, it is impossible to adjust the forest rotation length. However, for other events such as a drought or a pathogen invasion, there is a time gap between the onset of the event and its physical consequences on the forest stand. In this latter case, the forest owner can deliberately choose to prematurely harvest and to replant or to leave the forest standing. This choice is likely to depend on the age of the forest and on the intensity of the climate event. In the subsequent sections, we develop a methodology to explore this issue.

Material and methods

We define the CS as the climate information related to the date of occurrence and the intensity of the climate event. We define the intensity as the potential growth reduction and rise in mortality in the forest stand. Following [START_REF] Lechthaler | The climate challenge for agriculture and the value of climate services: application to coffee-farming in Peru[END_REF] and [START_REF] Adams | The benefits to Mexican agriculture of an El Niño-southern oscillation (ENSO) early warning system[END_REF], we developed a methodology based on the comparison between a situation with CS and a situation without CS. In order to deal with the multiannual time horizon, we used the forest's Land Expectation Value (LEV), i.e., the Faustmann (1849) criterion, as an economic value index between scenarios. The LEV is the present value, per unit area, of the projected costs and revenues from an infinite series of identical even-aged forest rotations, initially starting from bare land. It assumes that each rotation is of equal length, that the sequence of events within each rotation is the same, and that the net revenue associated with each event within a rotation is similar for all rotations. Contrary to other economic criteria such as the Net Present Value or the Discounted Profit, using the LEV makes it possible to compare the values of projects of different lengths, which is precisely what we want to do here.

The model

We use the standard LEV model in the case of even-age stands with constant discounted rate and prices over time. To simplify, we assume there are only two silvicultural operations during the rotation: plantation and final harvest. 1We define the Climate Service Value (CSV) as the difference between the LEV with CS and the LEV without CS.

= CSV LEV LEV CS noCS

(1)

where:

= LEV R T e D e ( ) (0) 1 CS CS rT rT CS CS
(2) and:

= LEV R T e D e ( ) (0) 1 no rT rT CS (3)
These equations are the so-called Faustmann formula [START_REF] Reed | The effects of the risk of fire on the optimal rotation of a forest[END_REF][START_REF] Clark | Mathematical Bioeconomics[END_REF]. The left component of the numerator in (2) and ( 3) is the revenue derived from the final harvest, whereas the right component consists of the initial plantation costs. In detail:

• T is the optimal harvesting age without climate service and T CS is the optimal harvesting age with climate service.

• R T ( )

CS

is the revenue from selling trees at age T ,

CS where = × R T V T p v T ( ) ( ) ( ( )) CS CS CS o V T ( )

CS

is the total volume per ha: . The Climate Service (CS) considered in this paper consists of the economic information provided by a model (or a set of models) to help forest owners to decide whether or not it is better to harvest the stand before the consequences of the drought occur. The CS is not about the date of occurrence of a drought, which we consider to be unpredictable, but about the intensity of its impacts in terms of an increase in mortality and a reduction in growth and, subsequently, about the value of the forest stand itself.

= × V T v T q T ( ) ( ) ( )
the price dependent on the unitary tree volume assumes that the type of wood use depends on the size of the trees. Schematically, the biggest trees are expected to be used in the roundwood transformation industries and sold at higher prices, whereas the thinnest are expected to be used as industrial or energy wood and sold at lower prices.

• D (0) is the initial plantation cost and is the discount rate.

Optimal rotation ages with and without climate services

The forest owner's problem is to determine the age that maximizes the forest's LEV, i.e., to find the date t such that:

= = + LEV t R t t r LEV R t 0 ( ) ( ( )) (4)
This equation can be interpreted as equating the marginal benefit of waiting for an additional year before harvesting ( )

R t t

( ) with the sum of the interest costs of postponing the next harvest ( × r R t ( )) plus the interest costs of postponing all the subsequent harvests ( × r LEV ) [START_REF] Reed | The effects of the risk of fire on the optimal rotation of a forest[END_REF]. Now, let us assume that a drought (designated CC for "Climate Change") occurs and impacts both tree growth through the unitary volume v CC and tree mortality through the number of trees q CC . The forest owner's problem with CS is then expressed as:

= = + LEV T R T T r LEV R T 0 ( ) ( ( )) CS CC CS CC CC CS CC CS CC CS CC CC CS CC _ _ _ _ _ _ (5) where = × × R T v T q T p v T ( ) ( ) ( ) ( ( )) CC CS CC CC CS CC CC CS CC CC CS CC _ _ _ _
Whereas the same problem without CS is expressed as:

= = + LEV T R T T r LEV R T 0 ( ) ( ( )) CC CC CC CC (6) = × × R T v T q T p v T ( ) ( ) ( ) ( ( )) CC CC CC CC
Eq. ( 5) implies that if the forest owner receives the climate service, she/he will adjust the optimal age of rotation T CS CC _ depending on the impacts of the drought CC over the unitary volume v T . In the case where the forest owner does not receive a climate service (Eq. ( 6)), the forest stand still undergoes the same drought CC, (i.e., v CC and q CC are the same as in Eq. ( 5)) but the forest manager is shortsighted and acts as if there were no drought by conserving the optimal harvest age without the climate service T .

In both cases, we assume that the drought occurs at the same date and with the same intensity. However, the optimal harvesting age is likely to be different since the forest owner with CS will adapt the rotation length accordingly.

To summarize, CSV is the difference in the LEV of the following two scenarios:

Scenario CS_CC:

Scenario CC:

• Forest owner has no information about the incidence of a drought at the beginning of the rotation

• Drought incidence takes place • Forest owner has information about the impacts of the drought, just after the incidence

• Optimal rotation age TCS CC _ is modified

• Forest owner has no information about the incidence at the beginning of the rotation • Drought incidence takes place • Forest owner has no information about the impacts of the drought • Optimal rotation age T remains the same as without a drought incidence

Modeling the impacts of a drought

Modeling plan

To solve Eqs. ( 5) and ( 6), we modelv t ( )

CC

, q t ( )

CC and p v t ( ( ))

CC for a maritime pine stand in Southwestern France. We present the modeling and calibration processes in Appendix 1.

We then choose to estimate the impacts of the drought CC on the forest stand by varying three parameters: is the year of occurrence of the impacts of the drought on the forest stand, s is the additional mortality rate (in %), and h is the rate of growth decline. If a drought occurs at year

, = v v v t t t 1 is modified such that = × t v t v h ( ) (1 )
CC t and q t is modified such that as

= × t q t q s ( ) (1 )
CC t

. We compute the optimal rotation ages (i.e., we solve Eqs. ( 5) and ( 6)) using Julia software. 2 We then compute the CSV by subtracting the LEV with and without CS, as mentioned in Eq. ( 1).

We calculate the CSV for a range of these three parameters and by varying the discount rate (Table 1) . This results in 4000 simulations.

Results: optimal age of rotation

Without a climate event, (i.e., s = 0, h = 0, φ=Ø), the optimal year of harvesting derived from Eq. ( 4) is T = 32 years for r = 2%. Now, if a drought takes place and if the forest owner receives the climate service, the optimal age of rotation changes according to the values of s, h and φ, as presented in Fig. 2.

Fig. 2 shows that the optimal harvesting age with climate services

T CS CC _

depends on the date of occurrence of the climate event . We distinguish two ranges of values for with different implications and we define a threshold value for that we refer to as .

For < , we have > T T

CS CC _ (solid blue line above red dotted line on Fig. 2), i.e., when the drought occurs at the beginning of the forest stand life, knowing the consequences of the drought on the stand dynamics tends to postpone the optimal harvesting age. This is because when the drought occurs early in the forest stand life, the cumulated standing capital is low and the profit made by harvesting and selling it would not compensate for new plantation costs. Therefore, the optimal economic option is to keep the impacted forest standing and conduct it up to age T CS CC _ when the forest volume will be big enough to generate sufficient revenue. We can observe on Fig. 1 that in this case, T CS CC _ decreases when increases: the earlier the climate event occurs, the later the age for which the forest volume will be big enough to be harvested. This is due to two factors. First, we assume that the growth reduction is permanent as of the date of the drought, which implies that the earlier it happens, the higher the cumulated growth loss would be. Second, the younger the trees are, the higher their growth rate is (see Fig. 7 in the Appendix) and, therefore, the greater the impact of a growth reduction. For the very same reasons, all other things being equal, increasing h or s postponesT

CS CC

_ . For T , we have < T T

CS CC _

. At that point, the forest owner's optimal strategy completely changes since it is now economically profitable to harvest the stand just after the drought but before the impacts of the drought on the stand occur in order to start a new forest stand from bare land once again. This threshold value appears quite early, between 15 and 20 years of age in our simulations, depending on the values of h and s. Increasing either h or s (or both) moves forward the threshold value and therefore reducesT

Results: climate service value

Once we obtainT CS CC _ , we can compute the LEV CS (Eq. ( 3)) and LEV noCS (Eq. ( 4)) to obtain the climate service value CSV (Eq. ( 1)).

Fig. 3 shows the variations in the CSV for r = 3%, s = 0.3 and h = 0.3.

We can once again distinguish two trends of results regarding the date of occurrence :

For < , the CSV is low and decreases with . In that case, as ex- plained in the previous section, the forest owner adapts by adjusting the optimal harvest age to For T , the CSV increases with since the value of the stand and, consequently, the value of the expected loss due to the drought increases as the stand ages. 3 In other words, the later the climate event occurs in the life of the stand, the higher the CSV will be.

Fig. 4 shows the variations in the CSV for s, h and the discount rate r. Three main results can be derived from this figure.

First, increasing h or s, all other things being equal, increases the CSV. This is rather intuitive and means that the value of climate services increases with the intensity of the drought. However, the relative importance of s and h differs according to the date of occurrence of the drought φ. A decline in the growth rate relatively affects the CSV value more than a rise in the mortality rate for early φ, whereas the opposite holds for late φ. This is because an additional mortality has the same impact regardless of the age of the stand. However, since the marginal growth is higher for young trees and tends to an asymptote as the stand matures (Fig. 9 in Appendix 2), an early drought will have a relatively greater impact on growth than a later one.

Second, the CSV is highly sensitive to the discount rate value. This is visible on Fig. 5 for a 1% increase in r and for several values of s and h. The discount rate r affects both the optimal harvesting age and the absolute CSV value. In addition, a higher discount rate gives a relatively greater importance to an early drought and less importance to a late one (see Fig. 10 in Appendix 2).

Third, when testing "reasonable" values for s, h and r, the CSV ranges from 0 to €4900/ha depending on the value of φ, s, h and r. However, we observe a large discrepancy of the CSV regarding the date of occurrence of the drought. When the drought occurs before , the CSV remains low and ranges between €0-50/ha. However, when drought occurs after ,therefore affecting an almost mature stand, the CSV is much higher and ranges between €0-4900/ha in our simulations. Consequently, the forest owner's willingness-to-pay for CS is likely to be non-linear regarding the age of the impacted stand, i.e., while it remains low for young stands, it rapidly increases as stands mature. One direct implication is that investments for CS should focus on mature forest areas.

Discussion

Our study is a first important step in defining a methodology to estimate the CSV in forestry. We propose four possible future research paths here to improve and expand our methodology.

First, relying on the Faustmann model leads us to assume an infinite rotation sequence with a drought occurring at the same age during the future rotations. In addition, using the Faustmann framework requires constant parameter values over time. One possible improvement consists in making these parameters dynamic. This point is extensively discussed in the literature. For instance, [START_REF] Buongiorno | Further generalization of Faustmann's formula for stochastic interest rates[END_REF] worked on a Faustmann model with stochastic prices and discount rates. [START_REF] Price | Optimal rotation with declining discount rate[END_REF][START_REF] Price | Optimal rotation with differently-discounted benefit streams[END_REF]) tested different discount rates patterns over time, whereas [START_REF] Chang | Application of the generalized Faustmann model to uneven-aged forest management[END_REF] extended the Faustmann model to uneven-age stands. In particular, it is likely that a risk-averse forest owner suffering from a drought during the first rotation will change her/his behavior for the second rotation, even without a climate service. For this purpose, [START_REF] Chang | A generalized Faustmann model for the determination of optimal harvest age[END_REF] proposes a generalization of the Faustmann formula that allows the harvest age to vary from timber crop to timber crop by varying the stumpage price, timber yield, regeneration cost and interest rate. However, empirical studies show that forest owners are usually risk-prone or risk-neutral [START_REF] Andersson | Assessing non-industrial private forest owners' attitudes to risk: do owner and property characteristics matter?[END_REF]. A recent study by [START_REF] Brunette | Attitude towards Risk and Production decision : an empirical analysis on French private forest owners[END_REF] shows that French forest owners are characterized by a relative risk aversion coefficient close to 1, i.e., they are risk-neutral. In addition, the length of rotation makes the "learning-by-doing" process limited in the forest sector since one individual forest owner might experience only one forest rotation during his/her active life.

Second, we assumed that growth decline is permanent as of the date of drought. In reality, trees undergo a time of recovery that depends on several factors: climate and carbon cycle dynamics, with biodiversity and CO 2 fertilization as secondary factors [START_REF] Schwalm | Global patterns of drought recovery[END_REF]. Data presented in [START_REF] Schwalm | Global patterns of drought recovery[END_REF] suggest that the current average time of recovery for the Landes forest ranges between one and two years. Yet, their model shows that given the changes in temperature as well as the increases in drought frequency and severity projected for the 21st century, terrestrial ecosystems will take longer to recover after droughts.

Third, we modeled the impacts of a single drought on growth and mortality during one rotation without considering the implications of multiple droughts during one forest rotation. In particular, do multiple droughts lead to a permanent growth decline with/without additional mortality or do they increase both growth decline and mortality? [START_REF] Bréda | Temperate forest trees and stands under severe drought: a review of ecophysiological responses, adaptation processes and long-term consequences[END_REF] show that long-term tree mortality is a complex function of the interactions among depleted carbohydrate stores in trees, decreased water transport efficiency and pest and disease outbreaks that affect weakened trees. Multiple droughts could therefore both lead to additional growth decline and a rise in mortality (possibly in a non-linear way since trees already weakened by a prior drought could die with a less intense drought than healthy trees). To consider these impacts, further work would consist of linking our modeling framework to a mechanistic eco-physiological model that estimates the mortality risk of forest trees under global change. One example is [START_REF] Davi | Assessing drought-driven mortality trees with physiological process-based models[END_REF] who use the CASTANEA model to simulate the development of tree functioning over time with different ontogenetic and phenotypic characteristics (age, diameter, Leaf Area Index, leaf traits) and that grow in different site conditions (elevation, soil water content).

Fourth, the combination of several climate events could make the CSV estimation even more challenging. [START_REF] Bréda | Temperate forest trees and stands under severe drought: a review of ecophysiological responses, adaptation processes and long-term consequences[END_REF] suggest that the consequences of a drought are more severe for trees already weakened by previous incidents such as storms, and the effects of such climatic accidents can occur more than 10 years later. In this vein, [START_REF] Manion | Tree Disease Concepts[END_REF] assumes that there is a three stage decline over many years: a long-term stress, a short term one like drought, and then death, occurring via a vector such as pathogens since climate change induces demographic changes that increase the number of pathogens [START_REF] Mcdowell | Mechanisms of plant survival and mortality during drought: why do some plants survive while others succumb to drought?[END_REF][START_REF] Rouault | Effects of drought and heat on forest insect populations in relation to the 2003 drought in Western Europe[END_REF].

Conclusion

Using a comparison of Land Expected Values, our analysis suggests that the Climate Service Values (CSV) for a maritime pine stand undergoing a drought is highly dependent on the age of the stand, on the intensity of the drought in terms of mortality and growth impacts and on the discount rate used. Overall, for a 2% discount rate, we show that for a young stand (less than 15-20 years old), the CSV is rather low, ranging from 0 to €50/ha. However, for a mature stand, the CSV rapidly increases, up to a maximum the year before the optimal harvest date. In the latter case, in the event of an intense drought, the CSV can reach €4900/ha.

The agricultural economics literature suggests that the CSV for coffee plantations in Peru ranges between $17-28/ha [START_REF] Lechthaler | The climate challenge for agriculture and the value of climate services: application to coffee-farming in Peru[END_REF]. For maize, [START_REF] Bert | Climatic information and decision-making in maize crop production systems of the Argentinean Pampas[END_REF] estimate the economic value of climate information at between $20-30/ha, with peaks at $320/ha corresponding to the occurrence of ENSO phenomena. In addition, when reviewing surveys about the economic value of seasonal climate forecasts for agriculture, [START_REF] Meza | Economic value of seasonal climate forecasts for agriculture: review of ex-ante assessments and recommendations for future research[END_REF] reported that the CSV considerably increases with risk aversion up to a certain point at which the agent is so risk-averse that he does not trust the weather forecasts, which are too uncertain. Our results do not directly compare with these findings since we compute the CSV for an infinite number of forest rotations, whereas agriculture studies calculate the values for one year. The annuity corresponding to a LEV of €4900 using a discount rate of 2% is 4900 * 0.02 = €98/ha/year. It appears that this value is still higher than the values given in the agricultural economics literature. One interpretation is that in the case of a forest, the drought destroys a capital that was built up over several years, whereas in agriculture, a climate disaster "only" ruins one year of production.

Moreover, in our case study, the values for s and h remain theoretical and disconnected from observations. In order to compute the CSV for a specific drought, it is necessary to rely on dendrochronological and inventory observations, which make it possible to estimate h and s. [START_REF] Allen | A global overview of drought and heat-induced tree mortality reveals emerging climate change risks for forests[END_REF] reviewed several studies that estimate tree mortality after the 2003 drought. It appears that mortality rates are highly dependent on the species and the location, and range between 0.8% and 80%. For France specifically, [START_REF] Renaud | Results from the European Network For Damage Monitoring, Département de la santé des forêts[END_REF] estimate the tree mortality following the 2003 summer drought at 1.2% for conifers at the national scale. The impact of drought on growth reduction is even more heterogeneous between forests and within a specific forest. Exploring the impacts of the 1994-1995 drought on Aleppo pine in Spain, [START_REF] Gazol | Aleppo pine forests from across Spain show drought-induced growth decline and partial recovery[END_REF] show that on 27 observation sites, 88% of the trees showed a growth reduction of about 60% the year after the drought. The values presented in Fig. 3 are therefore rather conservative and show an intermediate level of impact on growth reduction.

Eventually, our method could apply in many other situations. Indeed, the significance of the Faustmann formula and its conditions is not restricted to forestry. As pointed out by [START_REF] Gaffney | Keeping land in capital theory[END_REF], the Faustmann formula ostensibly deals with timber growth but can be adapted to deal with all capital assets and with any time patterns of inputs and outputs whatsoever. When restricted to the forestry field, our methodology is likely to apply to any events (climate-related or not) that have time-delayed impacts, making it possible to adjust the forest rotation length. This could be the case for several forest pathogens (Heterobasidion annosum or Thaumetopoea pityocampa, for example, in the case of maritime pine).

Declaration of Competing Interest

None.

. To do that, we start from an observed silvicultural standard defined by the Société Forestière de la Caisse des Dépôts et Consignations (Table 2). This standard is made for a monoculture of maritime pine in the Landes de Gascogne Forest. It refers to a plantation of 1000 seedlings/ha with two thinnings performed at years 13 and 23, and a final harvest at year 35. We assume the two thinnings to be cost-neutral, so that the only expenditure is the initial plantation cost at year 0 with D(0) = €700/ha.

In practice, two thinnings are planned in order to both "concentrate" the timber volumes on a smaller number of individuals and to allow a dense initial plantation, which, in turn, favor slim logs with few intermediate branches and knots. On the basis of this standard, our "data crunching" objective is to obtain the dynamics of the forest in terms of the number of trees and volumes, as if the two thinnings are not performed and when the final harvesting age is not chosen ex ante. In other words, we aim at defining a "natural" dynamics of the planted forest.

For this, two main assumptions are made: (a) the evolution of the stand, measured through its mass (proxied by the inventoried volumes) in the absence of thinning would follow a logistic curve; and (b) thinning does not have any effect on the total growth and mortality at the stand level.

To do that, we start computing the survival rate on the basis of the number of trees in the silvicultural standard:

= survRate qstd qstd t t t 1 (7)
We then interpolate the series (excluding the two years of thinning) to obtain the "natural" survival rate natSurvRate (intercept: 0.99987; slope: -0.000185 -Fig. 6) that we used, in turn, to (a) compute the dynamic of trees q t that would have survived in the stand if thinnings were not performed; and (b) to decompose the "naturally dead trees", distinct from the removed (thinned) trees in the observed variation in the number of trees in the years of thinning.

Removed "thinned" trees are computed as:

= qr qstd qstd natSurvRate qstd t t t t t 1 1
Modeling V t and v t On the basis of the number of trees thinned, we computed the volume removed = Vr qr vstd c t t t t , where c t is an adjustment coefficient, in order to take the possibility that thinned trees may have lower-than-average volumes into account (we used the coefficients 0.65 and 0.7 for the two thinnings).

The stand volumes that would exist without thinning are computed as the volumes given by the silvicultural standard with thinning, plus all the removed volumes:

=

Vnothinnings qstd vstd Vr t t t j

Once we obtained Vnothinnings , t we then used the Julia Package LsqFit 4 to estimate a generalized logistic model:

obtaining the parameters K, B, M and G (519.7, 0.0926, -6.44 and 0.0874, respectively).

We then computed v t as = v V q / t t t .

Modelingp v ( ) t Timber prices are highly dependent on log size since the array of possible uses increases with size. We obtained the standing timber market price (stumpage) for each log's volume class (Table 3) from the Société Forestière de la Caisse des Dépôts et Consignations.

Using the LsqFit[1] package once again, we fitted the prices according to a logarithmic model (Fig. 8) with parameters 29.939 and 6.631 for the fixed and logarithmic terms, respectively.

Appendix 2: Relative importance of s, h and r on the CSV

Fig

  Fig.1. . The Climate Service (CS) considered in this paper consists of the economic information provided by a model (or a set of models) to help forest owners to decide whether or not it is better to harvest the stand before the consequences of the drought occur. The CS is not about the date of occurrence of a drought, which we consider to be unpredictable, but about the intensity of its impacts in terms of an increase in mortality and a reduction in growth and, subsequently, about the value of the forest stand itself.

.

  As increases and moves towards , the harvest is less and less postponed and the T

Fig. 2 .

 2 Fig. 2. . Optimal harvesting ages when the forest owner obtains the climate service. For each chart, the x-axis indicates the date of the occurrence φ (the green dotted line refers to φ = 17 years as a benchmark for all the scenarios), and the y-axis indicates the optimal harvesting age with climate services T CS_CC (solid blue line), where the red dotted line refers to T = 32 years. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .Fig. 5 .

 45 Fig. 4. Climate service values (CSV) for different r, φ s and h. For each chart, the x-axis indicates h and the y-axis indicates s. Colored contour lines indicate the climate service values in €/ha.

Fig. 6 .

 6 Fig. 6. Observed and estimated tree survival rate.

Fig. 7 .

 7 Fig. 7. Comparison between forest volumes (with thinning) from the standard and computed volumes without thinning.

Fig. 8 .Fig. 9 .Fig. 10 .

 8910 Fig. 8. Observed and fitted timber prices per log volume.

Table 1

 1 Values of the parameters tested.

	CS CC _ , contrary to what

Table 2

 2 Sylvicultural standard defining forest operations and dynamics as given by the Société Forestière de la Caisse des Dépôts et Consignations.

	age	qstd	vstd	age	qstd	vstd
			(m 3 /tree)			(m 3 /tree)
	1	1000		19	487	0.31
	2	1000		20	485	0.34
	3	1000		21	483	0.38
	4	1000		22	480	0.42
	5	1000		23	321	0.5
	6	1000		24	320	0.55
	7	1000		25	319	0.6
	8	1000		26	318	0.64
	9	1000		27	317	0.69
	10	1000		28	315	0.74
	11	997	0.06	29	313	0.79
	12	994	0.07	30	311	0.83
	13	496	0.11	31	309	0.88
	14	495	0.14	32	307	0.93
	15	494	0.17	33	305	0.97
	16	493	0.2	34	303	1.02
	17	491	0.24	35	301	1.07
	18	489	0.27

4 John Myles White, LsqFit Julia package, v 0.3.0, https://github.com/JuliaNLSolvers/LsqFit.jl.

Table 3

 3 Stumpage prices per log volume class as computed with the fitted logarithmic model.

	v	p	v	p	v	p
	m 3 /tree	€/m 3	m 3 /tree	€/m 3	m 3 /tree	€/m 3
	0.05	10	1	30.5	2	34.74
	0.1	13.63	1.1	31.2	2.1	34.96
	0.2	18.71	1.2	31.84	2.2	35.16
	0.3	21.68	1.3	32.42	2.3	35.33
	0.4	23.79	1.4	32.97	2.4	35.47
	0.5	25.42	1.5	33.17	2.5	35.6
	0.6	26.76	1.6	33.56	2.6	35.7
	0.7	27.89	1.7	33.91	2.7	35.79
	0.8	28.86	1.8	34.22	2.8	35.86
	0.9	29.73	1.9	34.5	2.9	35.92
					3	35.96

In particular, we assume there is no other forest operation such as thinning.

https://julialang.org/.

Acknowledgments

This work was funded by the European project MARCO-730272 (MArket Research for a Climate Services Observatory) of the Horizon 2020 program. The BETA contributes to the French project Labex ARBRE ANR-11-LABX-0002-01. The authors thank Jemma Lemarchand for her valuable work.

Appendix 1:. Model calibration

To solve Eqs. ( 5) and ( 6), we need to model v t ( )

CC

, q t ( )

CC

and p v t ( ( ))

Appendix 3:. Supplementary material description

This article is accompanied by a supplementary material archive containing the following items:

• input_data.ods is the input data for the simulations, including q(age), v(age) and p(v) for the studied maritime pine stand, in OpenDocument format;

• model_code.{jl|pdf} is the modeling code using to produce the simulations, written in the Julia language;

• sim_results.ods includes the full output of our simulations, i.e., (LEV CS ,T CS ,CSV) as a function of (r, ,h,s), also in OpenDocument format.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.cliser.2019.100106.