M Dudley

Nail Akhmediev

Adrian Ankiewicz

J M S O T O -C R E S P O

Universal triangular spectra in parametrically-driven systems

Universal triangular spectra in parametrically-driven systems

Introduction

A defining feature of parametric systems in physics is the generation of multiple frequency components as a result of nonlinear energy transfer from an initial single-mode spectral excitation. In many systems, the spectra of the evolving field are associated with an exponentially decaying energy transfer from the pump frequency, leading to a characteristic universal triangular spectral (UTS) form when, as is usual, the energy content is analyzed or measured logarithmically.

Examples of such spectra abound in physics, and include spin waves in ferrites (yttrium iron garnet) (see Fig. 3 of [1]) and parametrically-generated second sound waves in superfluid He 4 (see Fig. 3 of [2]). On the other hand, the explanations that have been given for the specific spectral broadening in these works are complicated and not completely satisfactory. A more recent result comes from the physics of capillary waves [3], where results have been obtained through careful experiments and are presented in a significantly clearer way. This has allowed the authors to make a specific link between the parametrically-driven system and the nonlinear dynamics of modulation instability and rogue waves that leads to the appearance of triangular spectra associated with capillary waves on the excited surface of water.

The most recent example of the excitation of triangular spectra is from optics. In fact, the initial stage of supercontinuum generation in optical fibers, starting with a continuous wave, has this specific triangular spectrum, with increasing width [4,[START_REF] Dudley | Supercontinuum Generation in Optical Fibers[END_REF].U ntil recently, this phenomenon has been explained using ideas of multiple four-wave mixing [4] or soliton fission [START_REF] Husakou | [END_REF],b u tas i m p l e r and more general explanation, based on modulation instability and parametric frequency conversion, can also be given for this case.

Clearly, parametric processes are the major cause of the transfer of energy between spectral components. In most cases (including surface waves [3] and optical fibres [7]), four-wave mixing (FWM) [8] processes play the major role, and there is no doubt that a set of coupled equations relating a multiplicity of frequency components can describe the above phenomena [4]. However, when proceeding in such a manner, the analysis can become cumbersome and difficult to understand in clear physical terms. Fortunately, there is a way to reduce these complications by using a single (master) equation that intrinsically describes all the frequency transformations and their mutual interactions. For fourwave mixing, this "master equation" is the well-known nonlinear Schrödinger equation (NLSE).

For the NLSE, it is possible to show that the self-phase modulation term adequately describes the interaction between a potentially infinite number of frequency components. This would be extremely difficult to treat analytically via a coupled mode formalism. While, in some cases, a truncated mode approach, considering only a limited number of components, can be used [9,10], it has recently been appreciated that there exist (and have indeed existed in the literature for some years) a number of analytic techniques that allow the NLSE to be solved more easily, allowing a more complete description of the underlying physics than that provided by truncated models. This is particularly the case for the phenomenon of modulation instability (MI), where a single spectral component in an NLSE system, together with noise or low amplitude sidebands, undergoes "instability" during propagation and the energy "drifts" from the central frequency component to the sidebands. MI and four- 7), (8), profiles of the AB solution at the point of maximal compression, Z = 0. On a log scale, the spectrum has a distinctive triangular shape.

wave mixing are, in a sense, different views of the same process, in which there is a dynamic energy transfer between adjacent sidebands. A defining feature of the many previous numerical studies of the process has been the observation of a (logarithmic) triangular spectrum. Although MI and four wave mixing processes within the NLSE are often studied via numerical methods, the underlying dynamics can, in fact, be described by a series of exact solutions of the NLSE, now often referred to as Akhmediev breathers (AB) [11][12][START_REF] Ten | Reports of RIAM Symposium No.17SP1[END_REF][START_REF] Clamond | [END_REF][15]. Significantly, the spectra of the AB solutions are logarithmically triangular, and, in optics they reflect the evolution dynamics seen in experiments, possessing the tendency to increase in width from a single frequency to a wide spectrum covering almost the whole optical range [7].

In the present work, we concentrate on the evolution of spectra that can be described by the exact and simple analytical expressions described by the AB theory. Although we focus on particular systems of optics and hydrodynamics, our analysis should be well-suited to describing results in other fields [1][2][3][4] and should motivate further studies in systems such as cold atoms and condensates where NLSE-like dynamics are also observed.

Modulation instability theory

The MI theory from the 1986 paper [17] predicts general results for the distance-dependent evolution of the amplitudes of the AB spectral sidebands. The results here show this spectral evolution more clearly than in Ref. [7], and confirm the theoretical predictions explicitly. We start the analysis with the NLSE in dimensional form [START_REF]Nonlinear Fiber Optics[END_REF] i

∂ A ∂ Z + β 2 2 ∂ 2 A ∂ T 2 + γ |A| 2 A = 0, (1)
where the coefficients are defined in [START_REF]Nonlinear Fiber Optics[END_REF]. Using a simple change of variables,

Z = zL NL , T = t |β 2 |L NL , A = u P 0 , (2)
where L NL = (γ P 0 ) -1 ,E q .(1) can be converted to its dimension- less form:

iu z + 1 2 u tt +|u| 2 u = 0. (3)
Thus, any solution of (3) can be easily transformed into a solution of (1), and vice versa. The parameter P 0 in the transformation is the peak power if we are dealing with pulses and is the average power if we deal with continuous waves. Eq. ( 1) has a solution in the form of the AB given by: A(Z , T )

= P 0 exp i Z L NL × (1 -4a) cosh( bZ L NL ) + ib sinh( bZ L NL ) + √ 2a cos(ω mod T ) √ 2a cos(ω mod T ) -cosh( bZ L NL )
. (4) This solution shows growth-return evolution over -∞ < Z < ∞ and represents a family of solutions with a variable independent parameter ω mod , corresponding to a single frequency modulation on the initial continuous wave (CW) solution. The coefficients a and b depend on ω mod and are defined by:

a = 1 2 1 - ω mod ω c 2 , ( 5
) b = 8a(1 -2a) (6)
where Initially, at z →-∞, the solution is a continuous wave with a constant amplitude which subsequently develops a small modulation at ω mod when z becomes finite. The modulation increases when z → 0 and has maximal compression of the pulses in a train when Z = 0. The amplitude profile of the solution at Z = 0f o r maximal gain is shown in Fig. 1(a).

In the frequency domain, the AB spectrum consists of discrete frequency sideband modes with separation ω mod , and amplitudes that evolve in z according to:

f 0 (z) = e iz × i √ 2sinh(bz) + √ 2cosh(bz) - √ cosh(2bz) √ cosh(2bz) , (7)
f n (z) = e iz i √ 2sinh(bz) + √ 2cosh(bz) √ cosh(2bz) × √ 2cosh(bz) -cosh(2bz) |n| . ( 8
)
Here f 0 (z) describes the pump amplitude and the f n (z), where n =±1, ±2, ±3,..., are the sideband amplitudes. We have used the relation cosh(2bz) = 2 cosh 2 (bz) -1 to simplify the formulae given in [17].U s i n gE q .(2), the spectrum can be converted to dimensional units.

A unique feature of the spectrum is that, for |n| > 0, log

f n (0) √ 2 =-|n| log(1 + √ 2) ≈-0.88|n|. (9)
Hence, the points plotted in Fig. 1(b) (|n| > 0) form two straight lines. At z →-∞ , the spectrum consists of a single central nonzero component that contains all the energy. At finite z,t h e nearest sidebands appear at the expense of the central component.

The sidebands grow and create additional sidebands, thus filling the whole spectrum with equidistant components. The spectrum given by Eqs. (7), ( 8) at the point Z = 0 is shown in Fig. 1(b). It is clearly triangular and represents the result of successive energy transfer between the spectral components from the central mode to sidebands due to the four-wave mixing process. On further evolution, when z > 0, the tendency is reversed and energy flows from the sidebands back to the central mode.

Initially, only the pair of sidebands which is closest to the central component is unstable due to MI (see Fig. 2(a)). As soon as this pair grows, FWM coupling transfers energy to higher frequencies, generating additional sidebands even at frequencies where amplitudes initially were zero. Describing these dynamics in terms of a set of coupled equations for this potentially infinite number of spectral components is difficult [4], but in fact the recent realisation that this process is described by the AB theory means that such a multi-mode description is in fact not necessary: the exact AB solution intrinsically and exactly describes the ensemble of FWM processes that occur in this evolution.

AB solutions are heteroclinic orbits (see Fig. 1 of [17]) that start from a CW which is a saddle point of the dynamical system. Thus, a single growth-return evolution will be observed only if the trajectory starts along this heteroclinic one. If we choose a starting position at another nearby point, we generally will obtain periodic motion (see Eq. ( 18) of [17]). For example, the dynamics around the heteroclinic AB solutions can be studied by numerically solv- ing the NLSE for the evolution of a modulated CW input field of the form A(Z = Z 0 , T ) = P 0 1 + α mod cos(ω mod T ) .

In this case, we obtain the periodic motion shown in Fig. 2.T h e period along the z-direction in this motion depends on how much the initial condition deviates from the heteroclinic orbit, i.e. it depends on the value of α mod . In order to obtain an infinite period, α mod has to be a complex number.

For waves in an optical fiber, we can estimate the real values of the magnitudes involved. We use the following values of the fiber parameters at 1550 nm: β 2 = 20 ps 2 km -1 , γ = 1.1W -1 km -1 , and P 0 = 30 W. Then MI gain is observed for modulation fre- An alternative and visually more attractive way to view the spectral evolution is as a surface plot. This is presented in Fig. 3. Two growth-return cycles of spectral evolution are clearly seen. The central component f 0 is restored in between the cycles (point B). The initial spectrum is completely restored, i.e. all higher-order components disappear and the spectrum consists of only the central component and the pair of nearest sidebands.

To compare with the AB theory for the spectral amplitudes, we follow the evolution along the "return phase" of the trajectory from the first point of spectral expansion where there is maximum conversion into the sidebands to the next point in the evolution where the energy is transferred back to the pump. These two points are identified as A and B, respectively, in Fig. 3. Point A is the location of the maximally-compressed AB pulse that we can associate with z = 0 in Eqs. ( 7) and ( 8).F r o mt h e numerical simulations, we can readily determine the intensity of each sideband and then compare with the predictions of Eqs. (7) and (8).

The results are shown in Fig. 4. Here the origin, z = 0, corresponds to point A above, the markers represent the numericallydetermined sideband amplitudes and the solid lines are the predictions of Eqs. ( 7) and ( 8). There is excellent agreement. 7) and ( 8).T h e range of z here is chosen between the points A and B in Fig. 3.

More accurate spectra for a higher-valued α mod should be based on exact solution in terms of elliptic Jacobi functions g i v e ni nt h ep a p e r [17] (see Eq. ( 18) therein when κ → 1).

These may have direct application in analytic approaches to understanding beat signal reshaping and frequency comb generation.

Multi-frequency initial modulation

When the initial modulation consists of only a single pair of sidebands, the analysis is reasonably simple, as it was shown in the previous section. The case of two or more sidebands can also be analyzed using the nonlinear superposition of AB solutions. There is a well-developed technique for constructing higher-order solutions that start with modulation instability [16,[START_REF] Akhmediev | [END_REF].T h er esults are cumbersome and will not be reproduced here. Another difficulty is calculating Fourier spectra for higher-order solutions as there is no technique for deriving analytic spectra in those cases.

Despite the absence of analytical solutions, we can still use the fact that the system described by the NLSE is integrable to present a qualitative analysis for triangular spectra in this case. Indeed, each AB component in the nonlinear superposition defined by the NLSE is independent of the others, in the sense that it can be used to construct a nonlinear superposition. Then, for every additional pair of sidebands, the spectra are triangular and evolve periodically, just as in Fig. 3. An important point here is that we do not have to deal with the frequency that provides us with the highest growth rate. As the growth rate depends on the frequency within the instability band, each component grows at its own rate. Their corresponding points A for the components with small growth rate may move infinitely far from the beginning. The periods of evolution for different sidebands also vary, depending on the amplitude of each component. Thus, for two or more sidebands, at every z there will be a triangular spectrum related to one of the components. Thus, dealing with nonlinear superposition, we cannot expect a complete recurrence to the original mode. Moreover, if the initial modulation is chaotic, i.e. it contains all frequency components within the gain bandwidth, there will be a continuous spectrum at any z.

The qualitative behaviour described above is confirmed using numerical simulations starting from a CW with multiple frequency modulation. One example is shown in Fig. 5. Ten initial small spectral components within the instability band used for this simula-

Its spectrum consists of the central component containing almost all energy and ten small equal sidebands at each side. Each sideband pair grows according to its particular growth rate and, after the initial stage, the spectrum assumes a steady triangular shape.

If we ignore small deviations, the spectrum hardly changes during further propagation. The different curves shown in Fig. 5 for various distances of propagation demonstrate the convergence to a steady triangular spectrum. Convergence takes place for random initial modulation as well. Suppose, we take all spectral components with small average amplitudes, namely we choose as the initial condition:

A(0, T ) = √ 30 1 + 0.0001a(T ) , (11)
where a(T ) is a stochastic Gaussian-distributed variable with unit variance and Gaussian correlation function. The resulting spectra for this case, after the initial stage of evolution has passed, are shown in Fig. 6. The function a(T ) has been chosen as purely real (red solid curve) or imaginary (blue dashed curve). Although initially the spectra may develop with different speeds, the final result is always the same, as we can see from Fig. 6. After reaching the "steady" state, around Z = 20L NL , the width of the triangular spectrum, on average, does not change.

Comparing the spectra in Fig. 6 with those in Fig. 5 shows that random and regular perturbations lead to triangular spectra with the same width. Moreover, the width coincides with the width of the triangular spectra given by the analytical expressions of Eqs. ( 7), (8), and shown in Fig. 1. This is a result of independent evolution of a multiplicity of ABs in the complex field that have their points of maximum compression (B) distributed homogeneously along the z-axis. This means that Eqs. ( 7), (8) provide analytic results for universal triangular spectra in systems governed by the NLSE.

Clearly, any additional term in the NLSE that lifts integrability results in deviations from the triangular spectra. This can be seen in Fig. 3 of [7]. The spectra that develop at the initial stages of MI are triangular. However, the influence of Raman and higher-order dispersion terms create asymmetry and result in a wide supercontinuum. One more example is shown in Fig. 7(a). Here we use the initial conditions given by Eq. ( 11), but the NLSE has been modified with an additional third-order dispersion (TOD) term, β 3

∂ 3 A ∂ T 3 ,
with a small coefficient, β 3 = 1.557 ps 3 /km. At small z, triangular spectra are generated (not shown). However, by Z = 20L NL , the asymmetry has developed, and later on during the propagation ( Z = 30L NL ), it creates a wide supercontinuum spectrum, as shown in Fig. 7(a). The widths of the spectral curve with, and without, TOD evolve differently, as can be seen from Fig. 7(b). With the pure NLSE, the width reaches a steady value around 0.2 THz, while, with TOD, it increases approximately as √ z.

Conclusion

Ubiquitous triangular spectra in parametrically-driven systems described by the NLSE can find a simple explanation in terms of exact AB solutions of the nonlinear Schrödinger equation. In particular, the spectra of the AB solutions that start from noise-driven modulation instability demonstrate the triangular shape with varying width. The maximum spectral width is reached at the point of maximal compression of the AB. Modulation instability that starts from a noisy background, with all the frequencies being present, also develops a triangular spectrum whose width is also approximately the same as for the spectrum of a single AB at its widest part.

It is clear that analytic AB solutions of the NLSE provide an exact description of the cascaded parametric interactions that underlie the dynamics of MI, and their use and interpretation eliminates the necessity of using multiple coupled equations. Certainly for some cases, numerical simulations must be used to obtain exact quantitative predictions, but the qualitative dynamics predicted by the AB theory that yield logarithmic triangular spectra appear universal.

Fig. 1 .

 1 Fig. 1. (a) Temporal, Eq. (4), and (b) spectral, Eqs. (7), (8), profiles of the AB solution at the point of maximal compression, Z = 0. On a log scale, the spectrum has a distinctive

ω 2 c

 2 = 4γ P 0 /|β 2 | and P 0 the power of the CW field at large |Z |.T h es o l u t i o ni s valid over the range of modulation frequencies that experience MI gain, namely 0 < ω mod < ω c . The coefficient a then varies in the interval 0 < a < 1/2, while the parameter b(> 0) governs the MI growth. The maximum gain condition, b = 1, occurs for a = 1/4, i.e. ω mod = ω c / √ 2.

Fig. 2 .

 2 Fig. 2. Spectral and temporal profiles of the AB solution over two periods of evolution.

  quencies f m = ω mod /2π in the range [0, f c ] where f c = ω c /2π = 408.88 GHz. The maximum gain at a = 0.25 corresponds to a modulation frequency f m = f c / √ 2 = 289.12 GHz. The spectral and temporal evolution for this case with α mod = 0.01 are shown in gray-scale in Fig. 2. We show the temporal intensity |A(Z , T )| 2 and spectral intensity | Ã(Z , f )| 2 , where A(Z , T ) ↔ Ã(Z , f ) are a Fourier transform pair. The propagation distance z is scaled relative to the characteristic nonlinear length L NL = (γ P 0 ) -1 = 30.3m .

Fig. 3 .

 3 Fig. 3. Spectrum of periodic AB solution in 3-D format. The arrow A shows the point of maximal compression, while the arrow B is the point of first return (recurrence) to the starting conditions.

Fig. 4 .

 4 Fig. 4. Comparison of spectral components, f n , for the pump n = 0 and sidebands n =±1 ...± 6 obtained numerically (dots), and from exact solutions (7) and (8).T h e

Fig. 5 .

 5 Fig. 5. (Color online.) The spectra of the field with initial modulation given by Eq. (10) at distances of propagation Z = 10L NL (red solid line), 20L NL (blue dotted line), and 30L NL (green dashed line). Inset: MI gain curve with vertical bars that show the spectral components within the gain curve relating to the initial conditions (10).

Fig. 6 .

 6 Fig. 6. (Color online.) The spectra of the field that starts with random initial modulation, Eq. (10), at distances of propagation (a) Z = 50L NL and (b) Z = 200L NL .R e d( s o l i d )curves are calculated for purely real initial perturbation while blue (dashed) curves are for purely imaginary one.

Fig. 7 .

 7 Fig. 7. (Color online.) (a) The spectra of the field with random initial modulation (10), taking into account TOD, at a distance of propagation Z = 20L NL (red solid line) and Z = 30L NL (dashed blue line). (b) Evolution of spectral width in z with (blue dashed curve) and without (red solid curve) TOD.

Acknowledgements N.A. is grateful to Alexander von Humboldt Foundation, as this work was prepared while he was visiting Germany. The authors N.A. and A.A. acknowledge the support of the Australian Research Council (Discovery Project DP0985394). J.M.S.C. acknowledges support from the Spanish Ministerio de Ciencia e Innovación under contracts FIS2006-03376 and FIS2009-09895. J.M.D. thanks the Institut Universitaire de France and the French Agence Nationale de la Recherche projects MANUREVA ANR-08-SYSC-019 and IMFINI ANR-09-BLAN-0065, for support.


Asma Bourougaaoui

Mohamed L Ben Jamaa

Christelle Robinet
email: christelle.robinet@inrae.fr

Has North Africa turned too warm for a Mediterranean forest pest because of climate change?

Keywords: pine processionary moth, Thaumetopoea pityocampa, climate warming, range retraction, Pinus, Tunisia

Climate warming is inducing dramatic changes in species distribution. While many studies report the poleward range expansion of some species, some others report the range retraction and extinction risk of other species. Here we explore how climate warming affects the southern edge in North Africa of the pine processionary moth, Thaumetopoea pityocampa, which is a model insect currently expanding northwards and toward higher elevation in Europe. This Mediterranean forest insect was found in southern Tunisia until 2003. Field surveys were conducted to map the current southern edge of the species in Tunisia. Pheromone traps were installed on a north-south gradient, a translocation experiment of egg-masses was conducted on this gradient and local temperature change was analysed. We thus proved that the pine processionary moth has disappeared from southern Tunisia, and that no more adult males were actually flying there. We also found a decrease of egg-hatching and of the proportion of individuals able to reach larval stages along this gradient, while daily minimal and maximal temperatures globally increased. Furthermore, we showed that daily maximal and minimal temperatures as well as indices of extremely high temperatures have substantially increased during the study period . This study reveals the retraction of the pine processionary moth from southern Tunisia due to higher mortality rates that could be attributed to a significant local warming. The role of other factors (mainly the response of host trees and natural enemies to climate change) may amplify this direct effect and should be further explored.

Introduction

Climate change is currently evidencing an unequivocally increase in temperature with an average warming of 0.85°C during the period 1880[START_REF] Imbert | Expansion d'un ravageur forestier sous l'effet du réchauffement climatique : la processionnaire du pin affecte-t-elle la biodiversité entomologique dans les zones nouvellement colonisées ?[END_REF][START_REF] Ipcc | Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, eds Core Writing Team[END_REF]. Extreme weather events, notably heatwaves, have also increased in both intensity and frequency in most parts of the world since 1950, and in particular in North Africa. For instance, in Tunisia, the mean annual temperature has increased by 1.4°C in the 20 th century and the increase has even reached +3°C locally in Tunis (northern Tunisia) [START_REF] Verner | Tunisia in a Changing Climate: Assessment and Actions for Increased Resilience and Development[END_REF]). Climate warming is significant throughout Tunisia, even in southern regions where average maximum temperatures can already exceed 40°C [START_REF] Verner | Tunisia in a Changing Climate: Assessment and Actions for Increased Resilience and Development[END_REF]. Climate change scenarios suggest that maximal temperature could further rise, for instance by 2.5-3.3°C in Tunis by the 2080s [START_REF] Verner | Tunisia in a Changing Climate: Assessment and Actions for Increased Resilience and Development[END_REF].

Higher temperatures and more frequent and intense heatwaves threaten lives as well as ecosystems [START_REF] Parmesan | Impacts of Extreme Weather and Climate on Terrestrial Biota[END_REF]. Indeed, a consequence of the ongoing global warming is a set of dramatic impacts on a broad range of organisms, particularly plants and invertebrates [START_REF] Parmesan | A globally coherent fingerprint of climate change impacts across natural systems[END_REF][START_REF] Robinet | Direct impacts of recent climate warming on insect populations[END_REF][START_REF] Walther | Ecological responses to recent climate change[END_REF]. One of these impacts is the distributional change for species whose geographic range is mainly limited by temperature. Until recently, most of the studies have reported that many species were expanding their distribution toward the poles and higher elevations due to climate change [START_REF] Andersen | Identification of winter moth (Operophtera brumata) refugia in North Africa and the Italian Peninsula during the last glacial maximum[END_REF][START_REF] Parmesan | Poleward shifts in geographical ranges of butterfly species associated with regional warming[END_REF][START_REF] Thomas | Birds extend their ranges northwards[END_REF]. Few species-level extinctions were reported so far, namely species of anurans including the golden toad (Bufo periglenes) (e.g., [START_REF] Pounds | Biological response to climate change on a tropical mountain[END_REF]) but, following species distribution models, [START_REF] Thomas | Extinction risk from climate change[END_REF] estimated that 15-37% of the studied species (Taxa of mammals, birds, frogs, reptiles, butterflies, other invertebrates and plants) could become extinct by 2050. The extinction risk is actually predicted to accelerate with the temperature increase, and up to one species out of six could become extinct in the worst climate change scenario (RCP8.5; 4.3°C rise compared to pre-industrial temperatures) [START_REF] Urban | Accelerating extinction risk from climate change[END_REF]. The overall positive or negative impacts of climate change on species distribution are thus very difficult to predict. It is notably unclear how the southern range is impacted for species that move their northern range limit northwards (e.g. southern pine beetle, Dendroctonus frontalis, in the United-States, [START_REF] Lesk | Threats to North American forests from southern pine beetle with warming winters[END_REF], or the pine processionnary moth, Thaumetopoea pityocampa, in Europe, [START_REF] Battisti | Expansion of geographic range in the pine processionary moth caused by increased winter temperatures[END_REF]. Modeling techniques such as species distribution models or phenological models can be used to predict insect distributions and thus assess the potential shifts of northern and southern limits (e.g., for the spruce budworm, Choristoneura fumiferana, in North America, [START_REF] Régnière | Predicting insect distributions under climate change from physiological responses: spruce budworm as an example[END_REF]). However, the southern limit is often not studied in the field and we have so far poor biological evidence to say whether some species face a real range expansion (e.g., the northern limit shifts northwards and the southern limit remains the same) or instead a range shift (i.e., both northern and southern limits shift northwards at about the same rate).

The pine processionary moth (PPM), Thaumetopoea pityocampa, is a well-known Mediterranean species that has been recognized by the Intergovernmental Panel on Climate Change (IPCC) to be one of the few species for which the underlying mechanism of response to climate warming was clearly demonstrated [START_REF] Rosenzweig | Assessment of observed changes and responses in natural and managed systems[END_REF]. Indeed, PPM is expanding northwards in Europe as an indirect effect of the rise of minimal temperatures on insect survival during the cold season [START_REF] Battisti | Expansion of geographic range in the pine processionary moth caused by increased winter temperatures[END_REF][START_REF] Robinet | Modelling the effects of climate change on the potential feeding activity of Thaumetopoea pityocampa (Den. & Schiff.) (Lep., Notodontidae) in France[END_REF]). The northern edge of its distribution has then been considered as an indicator of global warming in France (French Ministry of Ecology, ONERC 2010).

PPM is distributed not only in Southern Europe but also in Asia Minor and North Africa [START_REF] Roques | Climate Warming and Past and Present Distribution of the Processionary Moths (Thaumetopoea spp.) in Europe, Asia Minor and North Africa[END_REF], where summer temperatures can reach extremely high values. PPM occurrence has been documented at its southern edge in North Africa [START_REF] Roques | Climate Warming and Past and Present Distribution of the Processionary Moths (Thaumetopoea spp.) in Europe, Asia Minor and North Africa[END_REF] but the effect of climate change on the PPM distribution edge is unclear. PPM has expanded its distribution in the south of Algeria because of the massive pine planting program, The Green Dam, to combat desertification.

Climate conditions were likely suitable for the PPM establishment before these pine plantations. In Morocco, PPM has also spread due to large afforestation activities and introduction of exotic pine species. In addition, its southern limit is blurred by its altitudinal limit in the High Atlas. In Tunisia, no large afforestation was done since the 1990s and no mountainous barrier may limit the southern range of PPM. Thus, it offers a unique place to investigate the effects of climate change on the species' southern distribution. PPM had been seen in southern Tunisia until 2003 [START_REF] Azizi | Eclaircissement de la situation de la chenille processionnaire du pin dans le sud de sa répartition (le sud Tunisien)[END_REF].

However, since 2014, PPM has not been seen in this region during occasional field travels (Ben Jamâa, 2017), and it is unclear whether PPM is still present. Since monthly maximal temperatures above 25°C are detrimental for PPM populations [START_REF] Démolin | La processionnaire du pin en Tunisie[END_REF][START_REF] Robinet | Looking Beyond the Large Scale Effects of Global Change: Local Phenologies Can Result in Critical Heterogeneity in the Pine Processionary Moth[END_REF], and maximal temperatures usually exceed 40°C in July and August in southern Tunisia, the effects of their rise on these southern PPM populations are questioned.

In this context, two major objectives for this study were: 1) to update our knowledge about the PPM spatial distribution at its southern edge in Tunisia, and 2) to determine the effects of climate change on this possible distribution retraction. For this purpose, we first surveyed the current southern distribution of PPM. We installed pheromone traps to test the presence/absence of adult males, and we installed PPM colonies on a latitudinal gradient to measure their survival rate in relation with weather conditions. We also analyzed long time-series of temperatures along this gradient to characterize locally the intensity of observed climate change. This study on a model species known to expand northwards provides the first clear evidence about adverse effects of climate warming on its southern range.

Materials and methods

Study species

The pine processionary moth (PPM), T. pityocampa (Denis & Schiffermüller) (Lepidoptera, Notodontidae), is a univoltine forest insect species largely distributed in Mediterranean countries [START_REF] Démolin | Bioécologie de la processionnaire du pin, Thaumetopoea pityocampa Schiff. Incidences des facteurs climatiques[END_REF][START_REF] Roques | Processionary moths and climate change: an update[END_REF]. In Tunisia, adults generally fly in summer from 15 th July to 1 st August at high altitudes, and from 10 th August to 7 th September at mid-altitudes (Démolin et al.1968;[START_REF] Jamâa | Essai de lutte contre la processionnaire du pin : Thaumetopoea pityocampa Schiff. (Lep., Thaumetopoeïdae) à l'aide de Bacollus thuringiensis Kurstaki (ECOTECH-PRO)[END_REF]. Females lay eggs on pine needles, larvae hatch about one month later and develop in five instars (L1 to L5; generally from autumn to spring), they feed on pine needles and build collectively white silk nests on pine branches during the cold season. At the end of their development, larvae make a procession (occurring from the end of January to March depending on winter severity) to leave the tree and burry in the soil and pupate. In fact, pupation occurs in February in the mountains, whereas it occurs in September and January-February in the regions of coastline. After a flexible period of diapause, adults emerge. Each PPM life-stage has some thermal constraints given the time window at which they are exposed to. For instance, late instar larvae are vulnerable to cold in winter (< 0°C) while eggs and first instar larvae are vulnerable to heat in summer-autumn (> 25°C) across its range in Mediterranean countries [START_REF] Démolin | Bioécologie de la processionnaire du pin, Thaumetopoea pityocampa Schiff. Incidences des facteurs climatiques[END_REF][START_REF] Robinet | Looking Beyond the Large Scale Effects of Global Change: Local Phenologies Can Result in Critical Heterogeneity in the Pine Processionary Moth[END_REF]. In relation with these constraints, PPM phenology varies across bioclimate regions (i.e., adult emergence and egg hatching occur later in warmer areas, thus young larvae are less exposed to the highest summer temperatures ; [START_REF] Huchon | La bioécologie de la processionnaire du pin. Dispersion potentielle, dispersion actuelle[END_REF] and this phenology could furthermore be affected by climate change [START_REF] Robinet | Looking Beyond the Large Scale Effects of Global Change: Local Phenologies Can Result in Critical Heterogeneity in the Pine Processionary Moth[END_REF]. PPM, which is the most important insect defoliator of pine forests in Mediterranean countries, causes substantial impact on the tree growth [START_REF] Jacquet | Pine growth response to processionary moth defoliation across a 40-year chronosequence[END_REF]). In addition, PPM also causes important allergy risks to humans and animals because larvae can release urticating setae [START_REF] Battisti | Processionary Moths and Associated Urtication Risk: Global Change-Driven Effects[END_REF]). This species is thus not only a good model species to study the effects of climate change but also an important pest in everyday life in forestry and public health.

Historical distribution of the PPM in Tunisia

The PPM is one of the most serious pests of Aleppo pine forests in Tunisia [START_REF] Rive | La processionnaire du pin, notions de biologie et principes de lutte[END_REF]) and thus, it has been episodically surveyed in the past. This insect is historically present in all pine stands from north to south [START_REF] Démolin | La processionnaire du pin en Tunisie[END_REF][START_REF] Roques | Climate Warming and Past and Present Distribution of the Processionary Moths (Thaumetopoea spp.) in Europe, Asia Minor and North Africa[END_REF]) (Fig. 1B). First, [START_REF] Démolin | La processionnaire du pin en Tunisie[END_REF] conducted surveys to determine the potential dates of adult emergence across different localities. This dataset can be used to have PPM occurrence points in the 1960s (Fig. SM1-A). In the 1990s, PPM nests have been seen in three regions (Jebel Sidi Driss, Jebel Toual and Techine) around Gabes according to forestry technicians (Fig. SM1-B). Later, PPM occurrence was studied through egg-masses sampled by INRGREF in September-October 1992, 1993 and 1995 throughout different localities going from northern regions to the southernmost edge of natural Aleppo pine stands (Kasserine) (Fig. SM1-B).

Thereafter, PPM nests have been last seen in Tataouine in 2003 [START_REF] Azizi | Eclaircissement de la situation de la chenille processionnaire du pin dans le sud de sa répartition (le sud Tunisien)[END_REF] during field travels in Jebel Brourmet (Fig. SM1-C). In 2010, its presence can be inferred from egg-masses collected by INRGREF, and in 2013-2014 through field surveys conducted by [START_REF] Azizi | Eclaircissement de la situation de la chenille processionnaire du pin dans le sud de sa répartition (le sud Tunisien)[END_REF] (Fig. SM1-C). PPM absence was firstly assessed in 2013-2014 in Tataouine, Gabes and Medenine (Azizi 2015).

Mapping the southern edge of the PPM distribution in Tunisia

Due to characteristic signs of PPM presence (white silk nests, egg-masses, or signs of pine tree defoliation), the presence of PPM is relatively easy to detect on the field. However, recording locations of PPM presence was not enough to map the southern edge of the PPM distribution. We needed to record locations where PPM was absent as well in order to actually define the border. The mapping of PPM northern distribution in France has been standardized over a regular cell grid at a spatial resolution of 8 km × 8 km (see protocol for field sampling described by [START_REF] Rousselet | Assessing species distribution using Google Street View: a pilot study with the Pine Processionary Moth[END_REF].

We tried to derive a similar approach to map the southern edge of the PPM distribution in Tunisia.

However, it was difficult, if not impossible, to adapt this protocol since we had no baseline edge to rely on, and host trees could be very scarce and scattered in southern Tunisia. Therefore, we adopted an intermediate protocol: we reported the occurrence (presence/absence) of the PPM in potentially critical areas. In 2018, from the 2 nd to the 4 th of January, we visited the regions supposed to be at the southern edge or beyond (near Sidi Bouzid, Gabes and Tataouine; Fig. 12). We notably surveyed areas where pine stands were reported by the General Direction of Forests in Tunis (Fig. 2A), and the main roads along which some isolated pine trees could eventually be found. Trees outside forests could effectively play an important role in habitat connectivity and the presence of the PPM [START_REF] Rossi | Trees outside forests in agricultural landscapes: spatial distribution and impact on habitat connectivity for forest organisms[END_REF]. In addition to this extensive survey, we added observations done on October 10 th , 2017 near Gafsa during a preliminary survey trip.

Pheromone trapping

Sexual pheromone traps (BuxaTrap®) were installed at different dates in summer 2018 as PPM adults generally fly later in warmer areas [START_REF] Huchon | La bioécologie de la processionnaire du pin. Dispersion potentielle, dispersion actuelle[END_REF]. These traps were distributed on a latitudinal gradient going from north to south: Tunis, Testour, Kef, Sidi Bouzid, Gabes and Tataouine (Fig. 1). One to three traps per site were attached on the branches of Aleppo pine trees, Pinus halepensis, placed at 1.5 meter high, and checked weekly until the trap removal (end of September to end of October 2018). When several traps were installed on the same site, they were distant to 50 meters from each other.

Translocation of egg masses

General method. In order to assess the mortality rate under different bioclimatic zones, we chose to work on a thermal gradient consisting of three locations (with Aleppo pines): Zaghouan (core area), Sidi Bouzid (supposedly on the southern edge) and Tataouine (from where PPM may have disappeared) (Fig. 1). Since early larval stages are particularly sensitive to high temperatures [START_REF] Robinet | Looking Beyond the Large Scale Effects of Global Change: Local Phenologies Can Result in Critical Heterogeneity in the Pine Processionary Moth[END_REF], we aimed to measure how maximal temperatures recorded in Tunisia in autumn could impact the survival of eggs and early larval stages. Some PPM colonies were thus translocated and exposed to field weather conditions from September 2019 to January-February 2020 on this gradient. Then, they were retrieved before they left their nest (e.g., just before the end of their larval development and the processions). In this way, we avoid their potential establishment in the area where it has likely disappeared (even if the establishment risk was supposedly very low). This experiment was inspired from those done in France and Italy to explore the effects of low winter temperatures on the PPM survival at its northern and altitudinal edge [START_REF] Battisti | Expansion of geographic range in the pine processionary moth caused by increased winter temperatures[END_REF].

Field experiment in more details. Egg-masses were collected on P. halepensis near Sidi Bouzid on September 17 th , 2019 (i.e., as soon as egg-masses could be found in this site). The translocation of egg-masses was done on a latitudinal gradient, going from north (Zaghouan, within the PPM core distribution) to south (Tataouine, where PPM has likely disappeared), with an intermediate location at its southern limit (Sidi Bouzid) (Fig. 1). Our initial objective was to translocate 30 egg-masses at each location, but given the low population density, it was not possible and we decided to keep only 30 egg-masses for Tataouine. A total of 15 egg-masses have been translocated at Zaghouan on Sept 25 th , 2019 (three egg-masses per pine tree, on five trees), 15 egg-masses at Sidi Bouzid on Sept 17 th , 2019 (three egg-masses per pine tree, on five trees), and 30 egg-masses have been translocated at Tataouine on Sept 18 th , 2019 (five egg-masses per pine tree, on six trees). Since a very low number of pine trees were available for this experiment, we had to place more egg-masses per tree (five instead of three). Egg-masses and nests were then collected on Jan 7 th , 2020 from Tataouine and Sidi Bouzid, and on Feb 5 th , 2020 from Zaghouan.

Analysis of egg-masses and nests collected from field experiments. Egg-masses and nests retrieved from the previous experiments were brought to the laboratory at INRGREF (Tunis) and kept at ambient air temperature (25±2°C). For each egg-mass, scales covering the eggs were removed to count eggs and split them into four categories [START_REF] Tsankov | Structure and parasitism of egg-batches of a processionary moth population different from Thaumetopoea pityocampa (Den. & Schiff.) (Lep. Thaumetopoeidae) found in Bulgaria[END_REF][START_REF] Imbert | Expansion d'un ravageur forestier sous l'effet du réchauffement climatique : la processionnaire du pin affecte-t-elle la biodiversité entomologique dans les zones nouvellement colonisées ?[END_REF] following their observation under a binocular magnifier:

(1) hatched eggs, presenting a characteristic large and round larval emergence hole,

(2) non-hatched-eggs, having no larval nor parasitoid emergence holes,

(3) parasitized eggs, showing a characteristic small parasitoid emergence hole and/or having a dead parasitoid inside, and/or having a characteristic parasitoid meconium inside, and finally (4) other eggs, including dried, flat and non viable eggs, and dead PPM larva that failed to emerge and showed an irregular larval emergence hole.

For each nest, all the living and dead larvae were counted, and their larval instar was determined.

Temperature records. Small data loggers (iButtons ®) have been installed on the three locations from Dec 18 th , 2018 to Jan 7 th , 2020 in Tatatouine, from Jan 9 th , 2019 to Jan 7 th , 2020 in Sidi Bouzid and from Feb 2 nd , 2019 to Feb 5 th , 2020 in Zaghouan to record temperature every two hours. Each logger was attached in the shade of a branch of the selected pine tree among trees used for the translocation experiment. On each location, several temperature statistics were calculated to compare weather conditions: (1) monthly statistics: absolute minimal and maximal temperatures of each month, and mean temperature of each month, and (2) daily statistics: daily minimal and maximal temperatures, and daily mean temperature. Data analysis was done using the R software (R Core Team, 2018).

Statistical analysis related to the translocation experiments

As the percentages of the different egg categories and the living/dead larvae were not normally distributed, a Kruskal-Wallis H test (followed by a post-hoc Dunn's pairwise tests if the K-W test was significant) was performed to compare the results from the different studied locations. A pairwise Wilcoxon signed Rank test was applied to compare temperatures among the locations of the latitudinal gradient. 3 Results 10

Analysis of climate data

Mapping the PPM southern edge in Tunisia

In total, 186 locations were investigated. PPM was found in 25 locations (including one location visited in 2017 in north-east Gafsa) and absent in 161 locations (including one location visited in 2017 in south-east Gafsa) (Fig. 2B). All the southern sites visited around Gabes and Tataouine were free from PPM, whereas the most southern sites infested by PPM were the north-east of Gafsa and south-east of Sidi Bouzid.

Pheromone trapping

No PPM male adult was caught in the pheromone traps in Gabes and Tataouine, thus confirming that PPM is absent in the most southern sites. In contrast, 40-58 adults per trap were caught in the regions of Sidi Bouzid, Kef (wherein the highest number of captures occurred), Testour and Tunis (Table 2).

Translocation of egg masses

Comparison of hatched eggs. The proportion of hatched eggs clearly decreased from north to south (gradient Zaghouan -Sidi Bouzid -Tataouine) (Fig. 3; Tables SM2-SM4). In addition, in Tataouine, one whole egg-mass out of 30 (3%) has not hatched, whereas all egg-masses from Sidi Bouzid and Zaghouan provided some hatched eggs (Tables SM2SM4). There was a statistically significant difference in the percentage of hatched (Kruskal-Wallis X 2 = 7.06, df= 2, P=0.029) and unhatched eggs (Kruskal-Wallis X 2 = 12.19, df= 2, P= 0.002) between the different locations. There was a significant difference in hatched eggs between Tatatouine and Zaghouan (Dunn's pairwise test; P= 0.028, adjusted using the Bonferroni correction), and there was a significant difference in unhatched eggs, between Tatatouine and Zaghouan (P= 0.003, adjusted by the Bonferroni correction), and between Sidi Bouzid and Zaghouan (P=0.012, adjusted by the Bonferroni correction).

Comparison of parasitized eggs and other mortality factors. Surprisingly the rate of parasitized eggs was higher in Tataouine, including both the generalist egg parasitoids, Ooencyrtus pityocampae, and the supposedly specialist egg parasitoids, Baryscapus servadeii, and most of parasitoids survived and emerged because we found parasitoid meconium more than dead parasitoids inside eggs.

Nevertheless, this parasitism rate was very low (6% on average) and did not vary significantly from the other regions (X 2 = 1.67, df= 2, P=0.433). Other factors of egg mortality had little impacts (less than 11.5% on average) and did not show any clear pattern on the latitudinal gradient. The percentage of other egg mortality was not significantly different among locations (X 2 = 0.35, df= 2, P=0.836).

Comparison of larval survival. Unfortunately, when we retrieved PPM nests, none of them was found on the site of Zaghouan (they were likely removed by people). As a result, we could only compare larval survival between Sidi Bouzid and Tataouine. The proportion of living larvae was lower in Tataouine than in Sidi Bouzid (Fig. 3; Tables SM2-SM4). Surprisingly a few larvae succeeded to survive in Tataouine (46% of eggs succeeded their development in larvae). The percentage of larvae that succeeded to survive was not significantly different between Tataouine and Sidi Bouzid (X 2 = 0.83, df=1, P= 0.361). However, the percentage of dead larvae was significantly higher in Tataouine (12%) than in Sidi Bouzid (X 2 = 4.03, df=1, P=0.04). The highest larval mortality affected mainly young instars (L2 and L3) in Tataouine (table SM2). Besides, larval development was relatively faster in Tataouine compared to Sidi Bouzid. Percentages of L2, L3, L4 and L5 instars were respectively 11 versus 18%, 41 versus 54%, 41 versus 27% and 7 versus 1% in Tataouine compared to Sidi Bouzid (Tables SM2 andSM3).

Comparison of temperatures.

Temperatures globally increased from north to south (gradient Zaghouan -Sidi Bouzid -Tataouine) (Fig. 4). The average from Sept 1 st 2019 to Feb 5 th 2020 of daily minimal, mean and maximal temperatures was: 15.1°C, 20.1°C and 26.5°C in Tataouine; 14.3°C, 19.0°C and 26.0°C in Sidi Bouzid, and 11.7°C, 16.1°C and 24.0°C in Zaghouan. Daily minimal and mean temperatures were significantly different among the three locations (P< 10 -4 and P< 10 -9 ), while the maximal daily temperature was significantly different only between Tataouine and Zaghouan (P = 0.005), and Sidi Bouzid and Zaghouan (P=0.02) but not between Tataouine and Sidi Bouzid (P = 0.33) (pairwise Wilcoxon Signed Rank Test) (Fig. 5).

Analysis of climate data

Trends in annual maximal, minimal and mean temperatures. All temperature variables (Tmax, Tmin and Tmean) showed a significant rising trend for all sites (Fig. 6; Table 3). Within 39 years , maximal temperatures have thus increased by 1.2°C, 1.7°C, 1.2°C respectively in Tunis, Sidi Bouzid, Tataouine and they have increased by 0.8°C within 28 years in Zaghouan.

Minimal temperatures have increased by 2.6°C, 1.5°C and 1.7°C within 39 years respectively in Tunis, Sidi Bouzid, Tataouine, and they have increased by 2°C within 28 years in Zaghouan. The highest annual temperatures are recorded in Tataouine with Tmax = 28°C, Tmin = 16.5°C and Tmean = 22.5°C (Fig. 6).

Trends in extreme temperatures: warm indices. The number of warm nights (TN90p) significantly increased at all the stations (Table 4). Similarly, there was a positive trend on the number of warm days (TX90p) with significant increase observed only at Tunis and Sidi Bouzid. There was a positive trend on TXx at all stations (except Zaghouan, presented a negative insignificant trend) with a significant increase observed at Sidi Bouzid. Trend analysis for TNx showed a positive trend with significant increase throughout the study areas. The other indices (except SU40 and SU35, in Zaghouan) indicated an increasing trend. The number of days with a daily maximal temperature (TX) above 40°C (SU40) showed a significant increase at Sidi Bouzid and Tataouine. There was a highly significant increase in the number of days with a daily maximal temperature (TX) above 35°C (SU35)

at Tunis, Sidi Bouzid and Tataouine. The number of days with a daily maximal temperature (TX) above 25°C (SU25) presented a highly significant increase at Tunis, Sidi Bouzid and a significant increase at Tataouine (Fig. SM5-SM7).

Trends in extreme temperatures: temperature fluctuation indices. In diurnal temperature range (DTR), three sites (Tunis, Zaghouan and Tataouine) had significant downward trends (Table 4). The lowest DTR value was recorded in Tunis (in 2004), followed by Zaghouan (in 2018) and Tataouine (in 2017). Unlike the other regions, there was a positive but not significant trend for DTR at Sidi Bouzid.

Number of days with temperatures above 35°C and 40°C from June to October. Different patterns

were observed among months (SM8-SM13). Both Tataouine and Sidi Bouzid recorded extremely high temperatures from June to August, while Tataouine appears to be even warmer than Sidi Bouzid from September to October. Indeed, more than half of July and August months had temperatures higher or equal to 35°C in the southern regions (Sidi Bouzid and Tataouine) nearly each year (Fig. SM9 ; Fig. SM10). Similarly, the number of days with temperatures higher or equal to 40°C (even exceeded 45°C) was more frequent in the southern regions than the northern ones (SM9-SM10). In September, high temperatures were more frequent in Tataouine (exceeding 35°C for 8-20 days and 40°C for 1-8 days ; Fig. SM11) than in other sites. In October, temperatures exceeded frequently 35°C

(up to 10 days) in Tataouine (Fig. SM12). Although such phenomenon was observed also in the past (in 1990 and 1992), it was less frequent. Recently, several consecutive years reached such high temperatures (at least 4 days above 35°C were recorded each year from 2012 to 2016; Fig. SM12).

Discussion

Climate change has been demonstrated to be a major factor affecting populations of forest pests in numerous ways. Recent literature has mainly reported species expanding their distribution, such as T. pityocampa [START_REF] Roques | Climate Warming and Past and Present Distribution of the Processionary Moths (Thaumetopoea spp.) in Europe, Asia Minor and North Africa[END_REF]. In this species, range expansion in Europe has been explained by higher winter temperatures and thus more favorable weather conditions to survive in winter. On the other hand, temperature increase and drought may affect the southern edge of the species distribution (in North Africa). Thus, the present study explored the effects of climate warming on possible southern retraction based on experimental approaches conducted in the field.

Field evidence of a southern retraction

PPM was historically present in southern Tunisia (the presence was confirmed through the observation of PPM nests), notably in Tataouine [START_REF] Azizi | Eclaircissement de la situation de la chenille processionnaire du pin dans le sud de sa répartition (le sud Tunisien)[END_REF]. Following the survey of the PPM's southern edge, it appears that the PPM was not found neither around Gabes nor Tataouine, whereas it was present in the north-east of Gafsa and south-east of Sidi Bouzid. Pheromone traps confirmed that no PPM males were actually flying in these two areas. Former PPM nests (i.e., nests built the previous year(s)) were not found either. Following these results, we proved that PPM has effectively disappeared from south Tunisia and the reason for this southern retraction has then been explored.

Lower survival in warmer southern sites.

We observed a global decrease of egg hatching and proportion of living larvae (Fig. 3), and an increase of temperatures (Fig. 4) along the latitudinal gradient from north to south. Maximal temperatures as well as egg hatching rates were however not statistically different between the southern PPM edge (Sidi Bouzid) and the area where it has disappeared (Tataouine), but they were different from the northernmost site (Zaghouan) (Fig. 5). Minimal and mean temperatures were though statistically higher in Tataouine than in Sidi Bouzid during the study period (Fig. 5), and thus may not explain the same low level of hatching rates in both Sidi Bouzid and Tataouine. However, it could explain why larval development was substantially faster in Tataouine than in Sidi Bouzid (Tables SM2-SM3). Indeed these higher (minimal and mean) temperatures could trigger a more rapid development since insect development rate is usually linearly correlated to temperature, at least up to an optimal temperature threshold [START_REF] Rebaudo | Modeling temperature-dependent development rate and phenology in insects: review of major developments, challenges, and future directions[END_REF]. Maximal temperatures could be the main driver of the PPM mortality and disappearance in southern Tunisia. If the development rate of eggs increases with increasing mean temperatures, then young larvae hatching from these eggs appear earlier in autumn and are thus likely exposed to higher maximal and potentially lethal temperatures. We could thus expect a further retraction from the region of Sidi Bouzid in the coming years if minimal and mean temperatures continue increasing with climate warming. Difficulties in collecting egg-masses in this region provide evidence of the PPM population decline at its current southern edge.

Retrospective analysis of the role of local climate warming

Average annual maximal, mean and minimal temperatures significantly increased between 1980 and 2019 in all investigated sites. Except in Sidi Bouzid, minimal temperatures increased faster than maximal temperatures (Table 3). As already mentioned, increase of minimal temperature is likely beneficial to the development rate of the insect. However, if eggs could develop faster, young larvae would be exposed to higher summer temperatures (due to both their occurrence at an earlier timewindow in the season, and an increase -even if more moderate -of maximal temperatures) and this would potentially foster mortality rates. In Tataouine, given the observed warming trend (Table 3), the annual minimal and maximal temperatures have increased respectively by 0.7°C and 0.5°C over 16 years (since PPM was observed in 2003), which is a substantial warming. More importantly, extremely high temperature values are recorded more and more often in southern regions (Sidi Bouzid and Tataouine; SU40, SU35, TXx in Table 4). The increasing trend of these heat indices may have a detrimental effect on the early instars larvae and thus contribute to higher mortality rates, which may explain the disappearance of PPM from southern Tunisia. Although we found that some larvae could develop in Tataouine in 2019-2020, maximal temperatures were much lower that year than the previous years (Fig. 6) and heat extremes (e.g., TXX and SU35) were also lower (Fig. SM5).

Higher temperatures have mostly been recorded between 2009 and 2016 (Fig. 6).

Life-stages potentially impacted by extremely high temperatures

Looking at summer data more specifically, temperatures usually exceeded 35°C on more than 15 days per month, even 40°C for a few days, and reached sometimes 45°C in the southern regions in July and August (Fig. SM9-SM10). This period corresponds to the flight of adults, which may be subject to lethal temperatures (Ben Jamâa 2017). However, we should note that both Sidi Bouzid and Tatatouine had very similar maximal temperatures in summer (Fig. SM9-SM10). Regarding the autumn data, in September (period that coincides with the beginning of oviposition), temperatures above 35°C and 40°C were more frequent in Tataouine (Fig. SM11). These temperatures (in September) could negatively affect egg survival if they reach unusually high values. Indeed, [START_REF] Robinet | Are heat waves susceptible to mitigate the expansion of a species progressing with global warming?[END_REF] showed that egg-masses from a French PPM population were able to survive a daily maximal temperature of 40°C during 12 consecutive days. However, [START_REF] Rocha | Effect of heat waves on embryo mortality in the pine processionary moth[END_REF] showed that negative effects appeared on egg-masses from a Tunisian PPM population at 42°C after only three days of heatwaves, and none of them survived at 44°C. In October (period that corresponds to the development of young larvae), temperatures exceeded 35°C (up to 10 days) in Tataouine (Fig. SM12).

Young larvae are particularly sensitive to high temperatures [START_REF] Robinet | Looking Beyond the Large Scale Effects of Global Change: Local Phenologies Can Result in Critical Heterogeneity in the Pine Processionary Moth[END_REF]. First instar larvae from French and Portuguese winter PPM populations have a substantial mortality rate at 36°C [START_REF] Santos | Temperature niche shift observed in a Lepidoptera population under allochronic divergence[END_REF]. Thus, we suggest that high temperatures recorded especially in autumn in southern Tunisia had contributed to higher mortality rates and to the disappearance of the PPM from this area.

Other factors to explore

Even if the increase of temperatures could be largely involved in the retraction of the PPM southern range, other factors also need to be explored. Firstly, as already briefly mentioned, a change in PPM phenology could also potentially explain changes in its geographical distribution [START_REF] Robinet | Looking Beyond the Large Scale Effects of Global Change: Local Phenologies Can Result in Critical Heterogeneity in the Pine Processionary Moth[END_REF] and should be further explored in association with the direct effect of climate warming on PPM survival. Secondly, climate change can also affect the PPM host trees. The last IPCC report (2014) revealed that the frequency and intensity of droughts has likely increased in the Mediterranean areas especially in arid regions. Droughts could considerably reduce the needles' length [START_REF] Borghetti | Impact of long-term drought on xylem embolism and growth in Pinus halepensis Mill[END_REF]) and decrease the leaf life span [START_REF] Girard | Pinus halepensis Mill. crown development and fruiting declined with repeated drought in Mediterranean France[END_REF], and this may disturb the layout of eggmasses around the pine needles. In addition, it is currently unknown whether climate change affects the pine needles' composition and how this change could alter the larval survival. This indirect effect of climate change on host trees could be an additional mortality factor to consider. The effect of drought on host preference and performance was indeed demonstrated for the Mediterranean pine shoot beetle, Tomicus destruens, on maritime pine, Pinus pinaster [START_REF] Branco | Water stress affects Tomicus destruens host pine preference and performance during the shoot feeding phase[END_REF]. In addition, according to [START_REF] Henriksson | Effects of host shading on consumption and growth of the geometrid Epirrita autumnata: interactive roles of water, primary and secondary compounds[END_REF], leaf water content is one of the most important factors affecting the growth of defoliating larvae such as Epirrita autumnata. Hence, drought stress may have a detrimental effect on insect performance reflected as decrease in leaf water content, and/or increase in leaf toughness and in secondary metabolites concentration. However, there is little information available regarding the effect of drought on the PPM performance. Thirdly, increase of both droughts and temperatures in most arid conditions would increase P. halepensis stress, thus impacting their own growth and survival (Cherif et al. 2019;[START_REF] Lasanta | Aridification determines changes in forest growth in Pinus halepensis forests under semiarid Mediterranean climate conditions[END_REF]. In general, phenological plasticity of some tree species (such as Pinus sylvestris in Europe) may not allow the tree species to adapt to climate change at their distributional edge [START_REF] Duputié | Phenological plasticity will not help all species adapt to climate change[END_REF]. If the PPM host trees would no longer be able to survive in some regions, the PPM would also no longer be able to be there, unless shifting on other hosts. Lastly, another mortality factor to consider is parasitism.

Indeed, the rate of parasitized eggs was relatively high in Tataouine (Fig. 3). Although this parasitism rate increases mortality, it is quite low (6 % on average) compared to other parts of the PPM distribution (e.g. 24.4 % in France (Orléans) and 29 % in Bulgaria (Maglizh), [START_REF] Georgiev | Comparative studies of egg parasitoids of the pine processionary moth (Thaumetopoea pityocampa, Den. & Schiff.) in historic and expansion areas in France and Bulgaria[END_REF]) and may not be a key factor in the PPM retraction.

Contrasted effects of climate change on insects' range and outbreaks

Increase of outbreaks and spread of pests are often reported. Bark beetles, for example, can cause tree mortality due to severe droughts, which are amplified by higher temperatures, and this could trigger a positive feedback on the reproductive rate of the insect populations [START_REF] Raffa | Cross-scale Drivers of Natural Disturbances Prone to Anthropogenic Amplification: The Dynamics of Bark Beetle Eruptions[END_REF][START_REF] Raffa | Responses of tree-killing bark beetles to a changing climate[END_REF]. In contrast, range reduction and extinction risks are often reported for other species. In Britain, fine-resolution mapping revealed extinction gradients since 1970 for three of the four northern butterfly species: Aricia artaxerxes and Erebia aethiops retracted northwards, and Erebia epiphron retracted uphill [START_REF] Franco | Impacts of climate warming and habitat loss on extinctions at species' low-latitude range boundaries[END_REF]. Actual shifts of insect distribution are poorly documented so far, apart from predictions of bioclimate models. A good example documenting a shift in outbreaks is the larch bud moth, Zeiraphera griseana. Its cyclical outbreaks are actually shifting from medium to high elevations in the Alps, and this could be attributed to a better hostinsect synchronization with climate warming at those high elevations [START_REF] Johnson | Climatic warming disrupts recurrent Alpine insect outbreaks[END_REF][START_REF] Rozenberg | Direct and Indirect Analysis of the Elevational Shift of Larch Budmoth Outbreaks Along an Elevation Gradient[END_REF]). The study presented in this article, in association with the previous study of [START_REF] Battisti | Expansion of geographic range in the pine processionary moth caused by increased winter temperatures[END_REF] done on the PPM northern range expansion, thus provides the very first well-documented shift in species distribution due to climate warming, and it also reveals the causal mechanisms.

Conclusions

Increase of outbreaks and spread are often reported for some species while range reduction and extinction risks are often reported for other species. The scientific community is usually focusing on what is threatening (e.g., more damages due to pests, or disappearance of protected species), but the overall and comprehensive impact of climate change is not so often explored. Here we demonstrate that climate could become too warm for a Mediterranean forest pest in the south of its distribution. The next question to tackle is whether forest trees that could finally get rid of pests or pathogens could actually persist in such extreme climate conditions.

  Climate data source. Daily minimal and maximal temperatures were kindly provided by the National Institute of Meteorology (INM) in Tunis at the following stations: (1) Mograne (representative of the weather in Zaghouan, located 6.7 km away) during the period 1991 -2011, (2) Medenine (located at 47.6 km from Tataouine, characterized by desert climate as the latter) during the period 1980 -2011, (3) Sidi Bouzid during the period 1980 -2014, and (4) Tunis during the period 1980 -2014.Hereafter, for simplification, we will refer to the name of experimental sites: Zaghouan for the Mograne weather station, and Tataouine for the Medenine weather station. Since these datasets had missing data and do not provide recent weather conditions (between 2011 or 2014 to 2019)//fr.tutiempo.net/climat/tunisie.html). For the four stations, we checked that there was a good agreement between the two datasets (i.e., we compared daily temperatures in 2011 and they correctly matched). We thus considered the combination of these two datasets to have full timeseries without missing data until 2019.Calculation of annual maximal, minimal and mean temperatures. To explore the local warming trend during the period 1980-2019 in Tunis, Sidi Bouzid and Tataouine and during the period 1991-2019 in Zaghouan, we calculated the annual mean of daily maximal, minimal and mean temperatures (Tmax, Tmin, and Tmean, respectively). Then, we made a linear regression for each variable (Tmax, Tmin, Tmean) as a function of the years to elucidate the magnitude of warming on each site. These analyses were performed using R software (R Core Team 2018).Calculation of indices of extreme temperatures.Nine indices related to temperatures from those defined by the Expert Team on climate change Detection Monitoring and Indices (ETCCDMI) were considered (Table1;[START_REF] Peterson | Climate Change Indices[END_REF][START_REF] Zhang | RClimDex 1.0: User Manual[END_REF].The calculation of these indices as well as statistical analyses were performed using the RclimDex package of R software (https://github.com/ECCC-CDAS/RClimDex). Slopes of the annual trends of these extreme climate indices were calculated based on a least square linear fitting. A locally weighted linear regression using a loess smoother function in R[START_REF] Zhang | RClimDex 1.0: User Manual[END_REF] was also applied to each index. The statistical significance of trends was assessed based on the Student's t-test for the estimate of the slope at the 95% level (p value < 0.05). In addition to these annual trends, we explored on a monthly basis (from June to October) the fluctuation of the number of days at which daily maximal temperature exceeded 35 and 40°C.

FiguresFig. 1 :

 1 Figures

Fig. 2 :

 2 Fig. 2: (A) Map of forest pine stands on the 8 km x 8 km grid, as reported by the General Direction of Forests in Tunis, and (B) map of the PPM presence and absence in South Tunisia.

Fig. 3 :

 3 Fig. 3: Comparison of larval survival, egg hatching, parasitized eggs and other causes of egg mortality along the latitudinal gradient going from Zaghouan (within the PPM core distribution in North Tunisia), Sidi Bouzid (the southern PPM limit) and Tataouine (region of disappearance in South Tunisia). Error bars indicate the standard error.

Fig. 4 :

 4 Fig. 4: Comparison of monthly mean temperatures in Tataouine, Sidi Bouzid and Zaghouan between January 2019 (month 1) and February 2020 (month 14) following data loggers installed on field sites.

Fig. 5 :

 5 Fig. 5: Boxplot of daily minimal, mean and maximal temperatures recorded on sites from Sept 1 st 2019 to Feb 5 th 2020. Different letters on each panel indicate a statistical difference using a pairwise Wilcoxon Signed Rank Test.

Fig. 6 :

 6 Fig. 6: Annual mean of daily (a) minimal, (b) maximal, and (c) mean temperatures on a south-north gradient (going from Tataouine, Sidi Bouzid, Zaghouan to Tunis) between 1980 and 2019.

Acknowledgments

We are grateful to Jérôme ROUSSELET (INRAE-URZF, Orléans, France) for insightful advice on mapping the edge of species distribution and translocation experiment, to Olfa EZZINE and Samir DHAHRI (INRGREF, Tunis, Tunisia) for general discussions and support in the frame of the PhD of Asma BOUROUGAAOUI, and to Adel BEN ABADA (INRGREF, Tunis, Tunisia) for his valuable help in the field. Pheromone traps and protective suits were kindly provided by Jean Claude MARTIN (INRAE-UEFM, Avignon, France). The distribution of pine stands was kindly provided by the General Direction of Forests in Tunis, and a part of the temperature datasets by the National Institute of Meteorology (INM) in Tunis. This study was supported by funding from the University of Carthage (Tunisia): the University provided two grants to Asma BOUROUGAAOUI to spend her internships (two months/grant) at INRAE-URZF in France in 2018 and 2019 in the frame of her PhD. It was also supported by the ANR project PHENEC.

Funding

This study was supported by funding from the University of Carthage (Tunisia): the University provided two grants to Asma BOUROUGAAOUI to spend her internships (two months/grant) at INRAE-URZF in France in 2018 and 2019 in the frame of her PhD. It was also supported by the ANR project PHENEC (grant number 19-CE32-0007-01).

Availability of data and material

Most of the datasets are provided in the article and supplementary material. Other datasets are available from the corresponding author on request.

Code availability

The R code of analyses is available from the corresponding author on request.

Declarations

Conflicts of interest/Competing interests

The authors declare no conflicts of interest.

Authors' contributions

CR and MBJ designed the study. AB conducted the field experiments, the laboratory analysis and the statistical analyses ; CR and MBJ supervised her work. AB, MBJ and CR wrote and revised the article.

Ethics approval

This study was conducted on an insect species (the pine processionary moth, Thaumetopoea pityocampa), which is a pest species commonly controlled. In addition, this is a model species for science which has been exempted from any declaration regarding the Nagoya protocole.

Consent to participate

Not applicable.

Consent for publication

The authors of this article consent to publish the study.