Achim Wirth
email: achim.wirth@legi.cnrs.fr

Bertrand Chapron

Empirical evidence of a fluctuation theorem for the wind mechanical power input into the ocean

The ocean dynamics is predominantly driven by the shear-stress between the atmospheric winds and ocean currents.

The mechanical power input to the ocean is fluctuating in space and time and the atmospheric wind sometimes decelerates the ocean currents. Building on 24-years of global satellite observations, the input of mechanical power to the ocean is analysed.

A Fluctuation Theorem (FT) holds when the logarithm of the ratio between the occurrence of positive and negative events, of a certain magnitude of the power input, is a linear function of this magnitude and the averaging period. The flux of mechanical power to the ocean shows evidence of a FT, for regions within the recirculation area of the subtropical gyre, but not over extensions of western boundary currents. A FT puts a strong constraint on the temporal distribution of fluctuations of power input, connects variables obtained with different length of temporal averaging, guides the temporal down-and up-scaling and constrains the episodes of improbable events.

Introduction

The exchange of heat, momentum and matter between the atmosphere and the ocean has a strong influence on our climate [START_REF] Stocker | Climate change 2013: the physical science basis[END_REF]). Recent advances in satellite and in-situ based global Earth Observation (EO) systems and platforms, have significantly improved our ability to monitor ocean-atmosphere interactions. In the present work the exchange of momentum is considered. More precisely, we investigate the flux of mechanical power into the ocean mixed-layer at the ocean surface.

It is caused by the shear-stress at the surface due to the difference between the atmospheric winds and the ocean currents near the surface, in the corresponding planetary boundary layers. Various physical processes occurring at the air-sea interface on a large range of scales in space and time are important for the momentum transfer. For a general discussion on air-sea interaction we refer to [START_REF] Csanady | Air-sea interaction: laws and mechanisms[END_REF] and [START_REF] Veron | Ocean spray[END_REF]), for ocean energetics to [START_REF] Ferrari | Ocean circulation kinetic energy: Reservoirs, sources, and sinks[END_REF] and for wind work to [START_REF] Wunsch | The work done by the wind on the oceanic general circulation[END_REF]. The atmospheric winds are usually stronger than the ocean currents and therefore the atmosphere mostly loses energy at the interface by friction and the ocean mostly gains energy. As a feedback mechanism, the presence of surface currents will then modulate the air-sea transfer of momentum [START_REF] Bye | Large-scale momentum exchange in the coupled atmosphere-ocean[END_REF], [START_REF] Renault | Satellite Observations of Imprint of Oceanic Current on Wind Stress by Air-Sea Coupling[END_REF]). The energy exchange is not conservative 1 and most of the mechanical energy is dissipated [START_REF] Duhaut | Wind stress dependence on ocean surface velocity: Implications for mechanical energy input to ocean circulation[END_REF], [START_REF] Wirth | A Fluctuation-Dissipation Relation for the Ocean Subject to Turbulent Atmospheric Forcing[END_REF], [START_REF] Wirth | On fluctuating momentum exchange in idealised models of air sea interaction[END_REF]). In the present work we are not concerned with the details of the exchange in the respective boundary layers (see e.g. [START_REF] Veron | Ocean spray[END_REF]) but suppose that it is well represented through bulk formulas of air-sea interaction [START_REF] Fairall | Bulk parameterization of air-sea fluxes for tropical ocean-global atmosphere coupled-ocean atmosphere response experiment[END_REF]). In those models the power input is estimated based on the shear-stress at the surface and the ocean current near the surface and also depends on the sea state and the density stratification in the atmosphere and the ocean.

More precisely, we consider the mechanical energy exchange between the atmosphere and the ocean at a time, t, over a fixed surface area A of the ocean (the area which spans 10 o in the longitudinal and the latitudinal direction) . For this area we evaluate the mechanical power the ocean gains at the interface P(t). Due to the turbulent dynamics in the atmosphere and the ocean the quantities are fluctuating over a large range of scales in time and space.

We focus on two properties of the mechanical power input to the ocean at the surface: (i) on average the ocean gains energy at the interface P(t) > 0 (where . represents an average over the observation period and several surface areas A i ) and (ii) the power input is fluctuating, in time and space, due to the turbulent motion in the atmosphere and the ocean and negative events, with P(t) < 0, occur.

Today, fluctuations are the focus of research in statistical mechanics, which was traditionally concerned with averages.

Fluctuations in a thermodynamic system usually appear at spatial scales which are small enough so that thermal, molecular, motion leaves an imprint on the dynamics as first noted by [START_REF] Einstein | Zur theorie der Brownschen Bewegung[END_REF] (see also [START_REF] Einstein | Investigations on the Theory of the Brownian Movement[END_REF] and [START_REF] Perrin | Atomes (Les)[END_REF]). The importance of fluctuations is, however, not restricted to small systems. Fluctuations can leave their imprint on the dynamics at all scales when (not necessarily thermal) fluctuations are strong enough.

Turbulent fluid motions are typical examples (i.e [START_REF] Frisch | Turbulence: the legacy of AN Kolmogorov[END_REF]), for which average motions can not be understood or modelled without some knowledge about the turbulent fluctuations, and vice versa . Turbulent fluctuations can be especially pronounced in geophysical flows, which are highly anisotropic due to the influence of gravity and rotation. This leads to a quasi two-dimensional dynamics and an energy cascade from small to large scales and strong fluctuations (see i.e. [START_REF] Boffetta | Two-dimensional turbulence[END_REF] for a review on 2D turbulence). Likewise, the description of air-sea interactions on large time scales may not be understood without some knowledge of the fluctuations at smaller and faster scales. Furthermore, the research interest in many natural systems lies also in the fluctuations not only in an average state, weather and climate dynamics are examples where we focus on the fluctuations of the same system on different time scales. For the weather the time scale of interest is from roughly an hour to a week, for the climate the focus is from tenths to thousands of years and beyond.

A recent concept which is presently subject of growing attention in non-equilibrium statistical mechanics are Fluctuation Theorems (FT) (see e.g. [START_REF] Evans | Probability of second law violations in shearing steady states[END_REF], Gallavotti and Cohen (1995a), Gallavotti and Cohen (1995b), [START_REF] Ciliberto | Experimental test of the Gallavotti-Cohen fluctuation theorem in turbulent flows[END_REF], [START_REF] Shang | Test of steady-state fluctuation theorem in turbulent Rayleigh-Bénard convection[END_REF] and Seifert (2012)). Not only the average values of quantities like entropy, work, heat or other, are studied, but their fluctuating properties are scrutinised. There are different forms of FTs, reviewed in detail by [START_REF] Seifert | Stochastic thermodynamics, fluctuation theorems and molecular machines[END_REF].

In the present paper we focus on the FT put forward in Gallavotti and Cohen (1995a), Gallavotti and Cohen (1995b) and [START_REF] Gallavotti | Equivalence of non-equilibrium ensembles and representation of friction in turbulent flows: the Lorenz 96 model[END_REF], corresponding to the detailed fluctuation theorem in the limit of large averaging times. When the FT applies to a fluctuating quantity, as i.e. P(t) in the present study, it relates the probability to have a negative event, i.e. the ocean loses energy, to the probability of a positive event, i.e. the ocean gains energy, of the same magnitude. The FT is not concerned with instantaneous values but considers the fluctuations of temporal averages over varying averaging time. The FT, which is stated precisely in the next section, thus puts a strong constraint on the fluctuations of the quantity considered and its temporal averages of varying length.

FTs have been established analytically for Langevin type problems with thermal fluctuations [START_REF] Seifert | Stochastic thermodynamics, fluctuation theorems and molecular machines[END_REF]). Most experimental data comes also from micro systems subject to thermal fluctuations. The thermodynamic frame of the quantities considered, as entropy, heat and work is not necessary to establish FTs. Examples of non-thermal fluctuations are the experimental data of the drag-force exerted by a turbulent flow [START_REF] Ciliberto | Experimental test of the Gallavotti-Cohen fluctuation theorem in turbulent flows[END_REF]) and the local entropy production in Rayleigh-Bénard convection [START_REF] Shang | Test of steady-state fluctuation theorem in turbulent Rayleigh-Bénard convection[END_REF]). For these non-Gaussian quantities the existence of a FT was suggested empirically. Our work is strongly inspired by these investigations of the FT in data from laboratory experiments of turbulent flows.

In [START_REF] Wirth | On fluctuating momentum exchange in idealised models of air sea interaction[END_REF] the FT was investigated for three parameterizations of air-sea interaction and we refer the reader to this work for the theory and analytical solutions on fluctuating air-sea interaction in these idealised models. In that publication the concept of FT is also placed in a broader context of fluctuating dynamics and the relation to the fluctuation-dissipation-relation and the fluctuation-dissipation-theorem is given (see also [START_REF] Seifert | Stochastic thermodynamics, fluctuation theorems and molecular machines[END_REF] for a general discussion). Here we extrapolate the research of [START_REF] Wirth | On fluctuating momentum exchange in idealised models of air sea interaction[END_REF] by applying the concept of FTs to data derived from satellite measurements and discuss their relevance. It is important to notice that even in the case of the idealised models the FT was not established by analytical calculation, but it was confirmed numerically that the FT is obtained asymptotically, in the long-time limit, when the averaging time is larger than the time-scale of the slow ocean-dynamics.

The Fluctuation Theorem

We are interested in the mechanical power, P(t), absorbed by the ocean over a given surface area, A, of the ocean surface and an observation period t obs . We suppose that P(t) is a statistically stationary random variable, meaning that its statistical properties (mean value, moments and temporal correlations) do not change when shifted in time. Its statistical properties, at every instance of time, are completely described by its probability density function (pdf), p(z), which gives the probability that P(t) takes values between z 1 and z 2 by integration: Pr[z 1 < P(t) < z 2 ] = z2 z1 p(z)dz. The symmetry function is:

S(z) = ln p(z) p(-z) . (1)
It compares the occurrence of events when the ocean receives power of magnitude z to the occurrence when the ocean loses power of the same magnitude. We further denote the normalised energy received during an interval τ starting at time t 0 , by:

E(t 0 ) τ = t0+τ t0 P(τ )dτ t obs -τ 0 t0+τ t0 P(τ )dτ dt 0 /(t obs -τ ) , ( 2
)
where t obs is the total length of the available data record. The corresponding pdf is denoted by p(z, τ ) and the symmetry function by S(z, τ ). Note that the averaging starts at time t 0 and extends over the interval τ .

The Galavotti-Cohen fluctuation theorem (called FT in the sequel for brevity) holds for P if, for averaging times larger than a convergence time-scale of the FT (τ τ 0 ), two conditions are satisfied: (i) the symmetry function depends linearly on the variable z, and (ii) on τ :

S(z, τ ) = στ z, (3)
where σ is called the contraction rate. The contraction rate σ > 0 (see Gallavotti and Cohen (1995a), Gallavotti andCohen (1995b), Ciliberto et al. (2004) and [START_REF] Shang | Test of steady-state fluctuation theorem in turbulent Rayleigh-Bénard convection[END_REF]) depends on the problem considered. We did not manage to determine it from observed quantities. In systems where the fluctuation is due to thermal motion its value is related to the thermal energy, that is the product of the Boltzmann constant and temperature, k B T . When fluctuations arise from turbulent motion the temperature has (almost) no influence on the fluctuations and the contraction rate σ depends on the turbulence. Indeed, in the incompressible Navier-Stokes equations temperature does not appear explicitly and only the kinematic viscosity has a slight dependence on temperature. There is, therefore no reason why k B T is a governing parameter of the problem.

If the FT holds it is sufficient to know the probability for either z > 0 or z < 0 to obtain the whole pdf, when σ is known.

The FT therefore constraints "half" of the pdf, a strong constraint in the absence of an equivalent of the Boltzmann distribution.

This property also allows to calculate the probability of the rare events of z < 0 from frequent events z > 0.

For a dynamical system the FT may or may not hold and it might only be valid for a range of values. It was already noted in Gallavotti and Cohen (1995a) and Gallavotti and Cohen (1995b) that the FT might only be valid for values z < z * , when the large deviation function (see i.e. Touchette ( 2009)) diverges outside the interval [-z * , z * ]. More recently it was recognised that boundary conditions, that is the value P(t) at t = t 0 and t = t 0 + τ , can leave their signature in the symmetry function S(z, τ ), even when the limit of τ → ∞ is taken, whenever the pdf p(z) has tails which are exponential or less steep than exponential (see [START_REF] Farago | Injected power fluctuations in Langevin equation[END_REF], [START_REF] Van Zon | Extended heat-fluctuation theorems for a system with deterministic and stochastic forces[END_REF] and [START_REF] Rákos | On the range of validity of the fluctuation theorem for stochastic Markovian dynamics[END_REF]). In such case an extended FT (EFT) should be expected, which shows a linear scaling of the symmetry function near the origin with a transition to a flatter curve for larger

values. An analytic expression of the symmetry function, or the value of z * is obtained only for very idealised cases and the results presented here are empirical.

Power Input

The calculations of the power input to the ocean are based on the shear-stress at the surface and the ocean velocity. The shearstress is usually evaluated, based on the difference between the horizontal wind velocity u s a , usually taken at 10m above the ocean surface and the horizontal ocean surface-current u s o , using the quadratic drag law (see i.e. [START_REF] Renault | Satellite Observations of Imprint of Oceanic Current on Wind Stress by Air-Sea Coupling[END_REF]):

F = C d (u s a -u s o ) 2 )(u s a -u s o ). (4)
The drag coefficient C d depends on the sea-state and the stratification in the atmosphere and the ocean, it is obtained using bulk formulas [START_REF] Fairall | Bulk parameterization of air-sea fluxes for tropical ocean-global atmosphere coupled-ocean atmosphere response experiment[END_REF]). Variations of the drag coefficient are not considered and all the results are independent of a constant C d .

To obtain the power input, the vector product between the shear-stress and the ocean current-velocity is taken:

P(t) = F • u o . (5)
For the work done on the large-scale geostrophic-circulation, [START_REF] Wunsch | The work done by the wind on the oceanic general circulation[END_REF] and [START_REF] Zhai | On the wind power input to the ocean general circulation[END_REF] used the surface geostrophic velocity estimates from altimetry for u o . Using model data, [START_REF] Rimac | The total energy flux leaving the ocean's mixed layer[END_REF] chose the velocity at the surface to calculate the total power input, to then evaluate that only a fraction of this power is transmitted to the interior ocean at the base of the mixed layer. In the present work, largely building on 15-m drogued drifter velocities [START_REF] Rio | Beyond GOCE for the ocean circulation estimate: Synergetic use of altimetry, gravimetry, and in situ data provides new insight into geostrophic and Ekman currents[END_REF]), we use for u o the estimation of the current velocity at 15m depth. By taking the velocity at 15m rather than at the surface we exclude the power that is promptly dissipated by viscosity in the Ekman-spiral, but include the power that is supplied to the ageostorpic and geostropic dynamics at the core of the mixed layer.

Data

In this study, we build on the newly released GlobCurrent products, now available via the Copernicus Marine Environment

Monitoring Service (CMEMS, http://marine.copernicus.eu/services-portfolio). Essentially building on the quantitative estimation of ocean surface currents from satellite sensor synergy, the production has been performed of a 25 years reanalysis of global, 1/4 o maps of ocean currents at two levels, the surface and 15m depth. The data is obtained from the combination of altimetry, GOCE, wind and in-situ data (largely building on 15-m drogued drifter velocities) [START_REF] Rio | Beyond GOCE for the ocean circulation estimate: Synergetic use of altimetry, gravimetry, and in situ data provides new insight into geostrophic and Ekman currents[END_REF]).

Strongly based on altimeter data, this global ocean surface current product, and also similar global observation-based products [START_REF] Bonjean | Diagnostic model and analysis of the surface currents in the tropical Pacific Ocean[END_REF] Satellite winds are from the Copernicus project (http://marine.copernicus.eu/services-portfolio/access-to-products/). They are from scatterometer and radiometer wind observations. It is a blended product based on the different missions (ERS-1, ERS-2, QuikSCAT, and ASCAT) available at 1/4 o spatial resolution and every 6 hours and is described in [START_REF] Bentamy | Homogenization of scatterometer wind retrievals[END_REF] and [START_REF] Desbiolles | Two decades [1992-2012] of surface wind analyses based on satellite scatterometer observations[END_REF]. The data record for which wind and current data is available extends over 24 years, 1993-2016, at the same resolution in space and time.

The FT is a property that concerns the tails of a pdf, and it is necessary to consider a large amount of data, as provided by the GlobCurrent products. Still, a time record of 24 years of data coverage at a single location is too small for empirically suggesting or refuting the existence of a FT. To increase the amount of data, we use different tiles A i that obey similar statistical properties. The tiles represent an effective area of 0.5 o in the longitudinal and latitudinal direction. The pdfs of E(t) τ are calculated for an interval that spans at least twice the standard deviation of each pdf from the origin.

Note, that the average is unity by definition. The pdfs are calculated with three different resolutions (bin sizes). The interval is separated into 21, 31 and 41 bins of equal size and the pdfs are obtained by counting the number of occurrences for each bin.

The symmetry function is only calculated when probabilities are lager than 10 -3 per bin, this led to an omission of bins in case ASG, only.

Results

The pdfs p(z, τ ) for the four domains and for different values of averaging times τ are presented in the left panels of figs. 1, 2, 3 and 4. All clearly display non-Gaussianity.

With increasing averaging period, the pdfs become more centred around unity, which is the average value, see eq. ( 2). This can be verified by comparing the graphs for τ = 125, 250, 500, and 1000 days. It is a consequence of the central limit theorem and occurrences of negative values become less likely for larger τ . In the right panels of figs. 1, 2, 3 and 4 the symmetry function divided by the averaging period is plotted. These plots are similar to those in [START_REF] Gallavotti | Equivalence of non-equilibrium ensembles and representation of friction in turbulent flows: the Lorenz 96 model[END_REF], who verified the FT in an idealised numerical model.

The verification of the FT, that is of eq. ( 3), is two-fold. First, we verify the linear dependence of the symmetry function on z for different averaging periods τ , and determine the slope. Second, we verify that the slope is a linear function of τ for times larger than the convergence time-scale of the FT, τ > τ 0 , of the system. That is, we first have to confirm that the lines in the right panels of figs. 1, 2, 3 and 4 converge towards straight lines for increasing averaging periods and second we see if the lines superpose when increasing averaging periods. This is demanding, and a large amount of data is necessary. For the first point, we have to consider the pdf for an extended range in z, including the tails, asking for ensemble sizes (number of intervals of length τ ) large enough so that we can observe a clear scaling behaviour. For the second point, we have to increase τ to verify convergence. Furthermore, for larger τ the pdfs are more and more peaked around unity and negative events become less and less likely.

For the four domains, we observed a convergence of the normalised symmetry function with increasing averaging time. This indicates the existence of a large deviation principle (see i.e. [START_REF] Touchette | The large deviation approach to statistical mechanics[END_REF]). For the domains within the recirculation area (ASG and PSG) of the subtropical gyre a convergence towards a linear variation with z is observed in less than τ 0 ≈ 1 year.

This points towards the existence of a FT, as both points put forward at the beginning of the previous paragraph are observed.

For the extensions of the western boundary currents (GSE and KUE) , the convergence does not achieve a linear behaviour of the normalised symmetry function. This shows that a FT does not hold, as the first point put forward at the beginning of the previous paragraph is not satisfied. The contraction rate σ is the slope of the curves in the right panels of figs. 1, 2, 3 and 4. To estimate the alignment of the points for τ = 1250 days, we constructed an index γ: the slope of the normalised symmetry function from the origin to the first bin divided by the slope from the origin to the last bin. A value γ = 1 indicates a perfect alignment of the first bin with the last.

The index is presented in table 1 for the four different domains and three different resolutions of the pdf. For the recirculation area of the subtropical gyre cases, the index varies around unity for the different bin sizes. It is significantly greater than unity in the Gulf Stream and the Kuroshio extension for all bin sizes considered.

We did not attempt to present error-bars in the figures and numbers in the tables, as uncertainties depend on the number of statistically independent events, that is the correlation time. In the case of air-sea interaction there are correlations due to the atmospheric dynamics (mostly synoptic), the ocean dynamics, the annual cycle, interannual variability and a climatic trend.

How these processes contribute to the tails of the pdf's, to improbable events, is currently a hot topic in climate science (see i.e. [START_REF] Ragone | Computation of extreme heat waves in climate models using a large deviation algorithm[END_REF]).

Discussion

We obtain evidence that a FT applies to data within the recirculation area of the subtropical gyre in the Atlantic and the Pacific Ocean. In these cases the FT can be used to estimate the occurrence of rare negative events from frequent positive events of the same magnitude for all averaging periods τ (measured in days). If the FT applies, the probability of the rare improbable negative events can be calculated from frequent positive events. The rare events when the ocean looses energy have recently been the focus of dedicated research (see e.g., [START_REF] Zhai | On the wind power input to the ocean general circulation[END_REF], [START_REF] Wirth | Determining the dependence of the power supply to the ocean on the length and time scales of the dynamics between the mesoscale and the synoptic scale, from satellite data[END_REF]). Improbable negative events also lead to strong transfer of energy to small-scale turbulence in the atmospheric and oceanic boundary layers, potentially causing strong mixing in the atmosphere and ocean. Improbable events are often key for the system in a variety of applications and are the focus of recent research in climate science [START_REF] Ragone | Computation of extreme heat waves in climate models using a large deviation algorithm[END_REF], [START_REF] Seneviratne | Changes in climate extremes and their impacts on the natural physical environment[END_REF]. As an example: in the Atlantic subtropical gyre (ASG) case the slope of the symmetry function is S(z, τ ) = 2.•10 -2 τ z, this means that an event of the magnitude z = -1 is p = exp(-2. • 10 -2 τ ) less likely than an event having the average value (z = 1) and an event of the magnitude z = -2 is p = exp(-4. • 10 -2 τ ) less likely than an event having twice the average value (z = 2). The variable τ gives the length of the averaging interval in terms of observations done every 6 hours, that is τ = 400 corresponds to a period of 100 days. A FT represents a tool to obtain the rare negative events from frequent positive events for all averaging times τ > τ 0 ≈ 1 year and demonstrates that, to leading order, the probability of negative events vanishes exponentially with the averaging time.

The FT does not seem to apply in the highly non-linear Gulf Stream extension for z 0.3 and Kuroshio extension z 0.5.

For these regions, the symmetry function follows a FT for small values of z, before the curve flattens. This resembles the behaviour found in the EFT (see section 2). Indeed, in these two cases (GSE & KUE) the tails of the pdf of P show pronounced super exponential tails and boundary values might be important leading to a behaviour predicted by an EFT. Nevertheless, a similar change of slope was also found using highly idealised models of air-sea interactions (discussed in Wirth ( 2019)), in which a friction term was added to the ocean. This suggests that other processes than air-sea interaction dominate in the extension of the boundary currents leading to a departure from a FT symmetry. When the scaling of the symmetry function flattens for higher power-input, the manifestation of a negative improbable event, versus a positive event of the same magnitude, becomes more likely.

During data analysis, we also found that a FT does not apply when islands or coastlines are present (not shown here).

Departure from a FT for the power input to the ocean is found where horizontal dynamics dominates over the vertical oceanatmosphere momentum exchanges. The influence of the horizontal transport of energy with respect to the injection of energy through the surface decreases with domain size considered, as the circumference of a domain grows linearly, whereas its surface growth is quadratic. Yet, determining the existence of a FT for larger ocean domains asks for more data, which is currently not available. Our results are purely empirical, a theory explaining why the power input follows a FT in some cases and not in others, is still missing.

Finally, we put the theory of FTs in the more general context of climate dynamics. A measurement, especially when coming from satellites always contains some averaging in space and time. A FT, when it applies, will help to relate averages over varying periods and is a powerful tool to guide the up and down-scaling of observational data in time and obtain the statistical information on shorter and longer time scales, which are not explicitly observed. More precisely, when the pdf of the power supply, and therefore also the symmetry function is known form observations for given averaging times the symmetry function can be calculated for shorter and larger averaging times and constrains "half" of the pdf. This is useful in down-scaling and the construction of statistical parameterizations of not directly observed dynamics over shorter time scales. On the other hand, the information can be useful for developing models for the persistence of events over large time-scales not yet observed. A FT can help to decide if the persistence in time of a phenomena is within the likeliness of the statistically stationary dynamics or due to external influences. Furthermore, when data from observations follow (or not) a FT, model data should do likewise. As such, the FT becomes a tool of investigating the fidelity of models.

Recently there is an increasing interest in the variability of atmosphere and ocean dynamics and in the exchange of the two, also at high frequency and their contribution toward the lower-frequency air-sea momentum and energy fluxes (e.g., [START_REF] Zhai | On the wind power input to the ocean general circulation[END_REF], [START_REF] Wirth | Determining the dependence of the power supply to the ocean on the length and time scales of the dynamics between the mesoscale and the synoptic scale, from satellite data[END_REF]). It was found that higher frequency wind forcing increases the mixed layer depth [START_REF] Zhou | The impact of high-frequency weather systems on SST and surface mixed layer in the central Arabian Sea[END_REF]).

There is evidence that long term variability of the atmosphere ocean system as the Madden-Julian Oscillation and El Niño-Southern Oscillation, are influenced by higher frequency wind forcing (e.g., [START_REF] Bernie | Impact of resolving the diurnal cycle in an ocean-atmosphere GCM. Part 1: A diurnally forced OGCM[END_REF], [START_REF] Terray | The role of the intra-daily SST variability in the Indian monsoon variability and monsoon-ENSO-IOD relationships in a global coupled model[END_REF]). FTs when they apply give a connection between events averaged over different lengths of time and can help to evaluate the impact of not explicitly resolved scales. When they do not apply they might help to identify a specific, non stochastic, mechanism responsible.

We conclude by looking at our results from the stand point of dynamical systems. Statistical mechanics of systems in equilibrium are described by the Boltzmann distribution, which is completely determined by the temperature. In non-equilibrium statistical mechanics no such universal distribution is known (see i.e. [START_REF] Derrida | Non-equilibrium steady states: fluctuations and large deviations of the density and of the current[END_REF], [START_REF] Touchette | The large deviation approach to statistical mechanics[END_REF] and [START_REF] Frisch | Turbulence: the legacy of AN Kolmogorov[END_REF]), but some quantities in some processes seem to follow a FT which constraints the pdf and might indicate some universality. The mechanical power-input to the ocean by air-sea interactions, as a forced and dissipative dynamical system, may thus belong to a class of particular non-equilibrium systems exhibiting a FT symmetry property and offer guidance for climate studies.

  , Sudre et al. (2013)), suffer from well-known limitations. The full spatio-temporal ocean dynamics is certainly not well captured, possibly missing part of the geostrophic component and a number of dominant ageostrophic signals (e.g. inertial oscillations). Also, accuracy is strongly reduced in the Equatorial Band where the geostrophic approximation fails, in coastal areas where altimetry accuracy decreases and where ageostrophic currents often dominate, and in the seasonally ice-covered Polar Seas. Nevertheless, this global ocean surface current product provides a consistent data set covering the last 25 years.

  For a trade-off between ensemble size and similar statistical properties, we choose to consider domains extending 10 o in the longitudinal and latitudinal direction, composed of 20 × 20 non-overlapping tiles each. Four domains are considered, the first is in the recirculation area of the subtropical gyre (20 o -30 o N , 20 o -30 o W ) of the North Atlantic (case: ASG), the second in the Gulf Stream extension (35 o -45 o N , 35 o -45 o W )(case: GSE). The third is in the recirculation area of the subtropical gyre of the North Pacific (15 o -25 o N , 150 o -160 o E)(case: PSG) and the fourth in the Kuroshio extension (30 o -40 o N , 150 o -160 o E) (case: KUE). The data record from which wind and current data is available extends over the 24 years, 1993-2016, at a resolution of 6h. At four occasions in time data was missing. The gaps were filled by linear interpolation.

Figure 1 .Figure 2 .Figure 3 .Figure 4 .

 1234 Figure 1. The pdf p(z, τ ) (left) and the symmetry function normalised by the averaging time S(z, τ )/τ (right) as a function of z for different τ (see caption). The variable τ gives the length of the averaging interval in terms of observations done every 6 hours. Data are for case ASG, res 0.5 o , 1993-2016, res 6h

  Author contributions. AW has performed the coding, writing was shared by both authorsCompeting interests. No competing interest

Acknowledgements. This work was funded by Labex OASUG@2020 (Investissement d'avenir -ANR10 LABX56). These data were provided by the Centre de Recherche et d Exploitation Satellitaire (CERSAT), at IFREMER, Plouzane (France) and CMEMS. Part of this work was performed when AW visited LOPS, Brest. We are grateful to Abderrahim Bentamy for explanations concerning the data and Mickael Accensi and Jean-Fancois Piolle for help with the data analysis.

Data availability. Data is available under: http://marine.copernicus.eu/services-portfolio/access-to-products/ and http://marine.copernicus.eu/servicesportfolio) Data availability. Data is available under: http://marine.copernicus.eu/services-portfolio/access-to-products/ and http://marine.copernicus.eu/servicesportfolio)


Nina Luisa Santostasi
email: n.santostasi@gmail.com

Olivier Gimenez

Assessing the dynamics of hybridization through a matrix modelling approach

Keywords: Genomic extinction, hybrids, hybridization, matrix models, Mediterranean common dolphin, population projections, prevalence, sensitivity analysis, striped dolphin, wolf 1

absence of reproductive isolation  Rare or depleted species are particularly vulnerable to genomic extinction  Genomic extinction depends mainly on demographic parameters of parental species  Maintaining healthy and abundant populations prevents genomic extinction

Introduction

Hybridization, defined as the interbreeding of individuals from genetically distinct populations, regardless of their taxonomic status [START_REF] Allendorf | The problems with hybrids: setting conservation guidelines[END_REF] is recognized as a relatively common phenomenon both in plants and animals [START_REF] Hewitt | Hybrid zones: natural laboratories for evolutionary studies[END_REF][START_REF] Olden | Ecological and evolutionary consequences of biotic homogenization[END_REF][START_REF] Grabenstein | Breaking barriers: causes, consequences, and experimental utility of human-mediated hybridization[END_REF]. Hybridization is most commonly observed between otherwise allopatric taxa that come into contact due to natural (natural hybridization) or anthoropogenic causes (anthropogenic hybridization, e.g., human encroachment or the release of non-native taxa). The widespread occurrence of natural hybridization is raising attention due to its implications in evolutionary biology [START_REF] Abbott | Genomics of hybridization and its evolutionary consequences[END_REF]. Additionally, the increasing occurrence of anthropogenic hybridization is considered a significant threat to biodiversity [START_REF] Rhymer | Extinction by hybridization and introgression[END_REF]Seehausen et al., 2008;[START_REF] Crispo | Broken barriers: human-induced changes to gene flow and introgression in animals[END_REF].

Hybridization can have different consequences for the evolution and conservation of species. If the fitness of hybrids is lower than that of parentals, hybridization can reinforce reproductive isolation between incompletely isolated species (Barton and Hewitt, 1985), but it can also cause extinction through demographic swamping [START_REF] Allendorf | The problems with hybrids: setting conservation guidelines[END_REF][START_REF] Wolf | Predicting the risk of extinction through hybridization[END_REF]. If the fitness of the hybrids is greater than or equal to that of parental individuals, hybridization can cause fusion of species (Seehousen et al., 1997;[START_REF] Allendorf | The problems with hybrids: setting conservation guidelines[END_REF], genetic swamping [START_REF] Allendorf | The problems with hybrids: setting conservation guidelines[END_REF], transfer of genetic material between species (potentially facilitating their adaptive evolution; [START_REF] Grant | Hybridization of bird species[END_REF][START_REF] Verhoeven | Population admixture, biological invasions and the balance between local adaptation and inbreeding depression[END_REF], and the origin of new species [START_REF] Demarais | Origin of Gila seminuda (Teleostei: Cyprinidae) through introgressive hybridization: implications for evolution and conservation[END_REF].

Understanding the potential consequences of hybrization is important to unveil evolutionary mechanisms such as how species integrity is maintained in the face of interspecific (and often intergeneric) gene flow [START_REF] Crossman | Hybridization in the Cetacea: widespread occurrence and associated morphological, behavioral, and ecological factors[END_REF] and how new species can arise from the introduction of new allelic combinations generated by hybridization. Furthermore, understanding anthropogenic hybridization dynamics can help identify effective and timely management actions for threatened species. To this end there has been an effort to calibrate management actions to the severity of the threat by classifiying hybridization in different types [START_REF] Allendorf | The problems with hybrids: setting conservation guidelines[END_REF][START_REF] Bohling | Strategies to address the conservation threats posed by hybridization and genetic introgression[END_REF] based inter alia on hybrids fitness and relative abundance of admixed individuals (i.e., hybrids prevalence, [START_REF] Santostasi | Use of hidden Markov capture-recapture models to estimate abundance in the presence of uncertainty: application to the estimation of prevalence of hybrids in animal populations[END_REF]. Depending on the hybridization type, a variety of management actions can be more effective or feasible to avoid genomic extinction: from hybrids removal and/or sterilization to the management of the human disturbances that cause hybridization in the first place [START_REF] Allendorf | The problems with hybrids: setting conservation guidelines[END_REF][START_REF] Bohling | Strategies to address the conservation threats posed by hybridization and genetic introgression[END_REF].

It is often difficult to define hybridization types, inter alia because of lack of information about admixed individuals prevalence and fitness. Projection models can tackle this uncertainty by simulating hybridization dynamics under different biological/evolutionary scenarios [START_REF] Wolf | Predicting the risk of extinction through hybridization[END_REF][START_REF] Fredrickson | Dynamics of hybridization and introgression in red wolves and coyotes[END_REF]. Sensitivity analysis and/or the simulation of different possible management actions can also be used to provide management recommmendatios [START_REF] Crouse | A stage-based population model for loggerhead sea turtles and implications for conservation[END_REF][START_REF] Cross | Using logistic regression to analyze the sensitivity of PVA models: a comparison of methods based on African wild dog models[END_REF].

The first attempts to model hybridization-extinction dynamics had a genetic focus and were based on changes in allelic frequencies at one or more loci [START_REF] Huxel | Rapid displacement of native species by invasive species: effects of hybridization[END_REF][START_REF] Ferdy | Extinction and introgression in a community of partially cross-fertile plant species[END_REF].

Compared to genetic models, ecological models place a greater emphasis on life-history traits, by explicitly examining the effects of fitness parameters (e.g., survival and reproductive rates) on the hybridization outcome. Within the ecological approach, two types of model have been used to model hybridization dynamics [START_REF] Hall | What can mathematical modeling tell us about hybrid invasions?[END_REF]: 1) individual-based models that simulate the contribution of each individual to the hybridization dynamics of the entire population (e.g., [START_REF] Thompson | Model-based analysis of the likelihood of gene introgression from genetically modified crops into wild relatives[END_REF][START_REF] Hooftman | Modelling the long-term consequences of crop-wild relative hybridization: a case study using four generations of hybrids[END_REF], and 2) population-based models that can be used when the mean fitness parameters of the main demographic stages are available (e.g., [START_REF] Wolf | Predicting the risk of extinction through hybridization[END_REF][START_REF] Campbell | Predicting patterns of mating and potential hybridization from pollinator behavior[END_REF]. Both modelling approaches were applied to hybridization in plant species [START_REF] Hall | What can mathematical modeling tell us about hybrid invasions?[END_REF][START_REF] Todesco | Hybridization and extinction[END_REF]. However, to our knowledge, few studies used individual-based models [START_REF] Fredrickson | Dynamics of hybridization and introgression in red wolves and coyotes[END_REF][START_REF] Nathan | A spatially-explicit, individual-based demogenetic simulation framework for evaluating hybridization dynamics[END_REF] and none used population-based models to simulate hybridization dynamics in animal species.

Here, we develop a population-based approach to project the dynamics of animal hybridization. By using a stage-based matrix model, and grouping individuals into genealogical categories, each described by their mean fitness parameters (i.e., age-specific survival, per capita fertility rate), our approach overcomes the need of realistic and accurate data at the individual level. To provide a practical example of this modeling approach, we illustrate its application to two case studies. The first refers to two delphinid species belonging to different genera, that interbreed in the Gulf of Corinth, Greece [START_REF] Bearzi | Dolphins in a scaled-down Mediterranean: the Gulf of Corinth's odontocetes[END_REF][START_REF] Antoniou | Evidence of introgressive hybridization between Stenella coeruleoalba and Delphinus delphis in the Greek Seas[END_REF]: the striped dolphin Stenella coeruleoalba and the common dolphin Delphinus delphis. The Mediterranean subpopulation of the former species is classified as Vulnerable, and the latter as Endangered in the IUCN Red List [START_REF] Bearzi | Ecology, status and conservation of short-beaked common dolphins Delphinus delphis in the Mediterranean Sea[END_REF][START_REF] Aguilar | Stenella coeruleoalba (Mediterranean subpopulation). The IUCN Red List of Threatened Species[END_REF]. Within the Gulf of Corinth, however, the subpopulation of common dolphin would qualify as Critically Endangered due to its small size (point estimate 22 individuals; [START_REF] Santostasi | A robust design capturerecapture analysis of abundance, survival and temporary emigration of three odontocete species in the Gulf of Corinth, Greece[END_REF] and high (≥50%) probability of extinction in three generations [START_REF] Santostasi | Common dolphins in the Gulf of Corinth are Critically Endangered[END_REF]. Recently confirmed hybridization with the much more abundant striped dolphin population (point estimate 1331 individuals; [START_REF] Santostasi | A robust design capturerecapture analysis of abundance, survival and temporary emigration of three odontocete species in the Gulf of Corinth, Greece[END_REF] threatens the persistence of common dolphins in the Gulf of Corinth [START_REF] Antoniou | Evidence of introgressive hybridization between Stenella coeruleoalba and Delphinus delphis in the Greek Seas[END_REF][START_REF] Santostasi | Common dolphins in the Gulf of Corinth are Critically Endangered[END_REF]. We regard this as a case of anthropogenic hybridization, considering that the steep decline of common dolphins due to anthropogenic impacts [START_REF] Bearzi | Ecology, status and conservation of short-beaked common dolphins Delphinus delphis in the Mediterranean Sea[END_REF] locally leads to a lack of available mates and to mating with more abundant species [START_REF] Frantzis | Mixed-species associations of striped dolphins (Stenella coeruleoalba), short-beaked common dolphins (Delphinus delphis), and Risso's dolphins (Grampus griseus) in the Gulf of Corinth (Greece , Mediterranean Sea)[END_REF][START_REF] Antoniou | Evidence of introgressive hybridization between Stenella coeruleoalba and Delphinus delphis in the Greek Seas[END_REF]. The possible hybridization outcomes for the two isolated dolphin subpopulations in the Gulf of Corinth have not been previously explored [START_REF] Antoniou | Evidence of introgressive hybridization between Stenella coeruleoalba and Delphinus delphis in the Greek Seas[END_REF].

In the second case study we evaluate the possible outcomes of hybridization between the wolf (Canis lupus) and its domestic counterpart, the dog (Canis lupus familiaris). Hybridization between wolves and dogs has been documented in several European countries and represents a well-known threat to wolf conservation [START_REF] Boitani | Action plan for the conservation of Wolves in Europe (Canis lupus)[END_REF][START_REF] Boitani | Wolf conservation and recovery[END_REF]. Still, little is known about the possible outcomes of hybridization between wolves and dogs. Hybridization with dogs may represent a problem for recovering wolf populations expanding into human-dominated landscapes, where few potential wolf mates compete with free-ranging dogs that are abundant and widespread [START_REF] Randi | Detecting hybridization between wild species and their domesticated relatives[END_REF][START_REF] Galaverni | Disentangling timing of admixture, patterns of introgression, and phenotypic indicators in a hybridizing wolf population[END_REF]. Although reproductive isolation due to behavioral or physiological barriers has been often assumed to contrast introgressive hybridization [START_REF] Vilà | Hybridization between wolves and dogs[END_REF][START_REF] Randi | Detecting rare introgression of domestic dog genes into wild wolf (Canis lupus) populations by Bayesian admixture analyses of microsatellite variation[END_REF][START_REF] Galaverni | Disentangling timing of admixture, patterns of introgression, and phenotypic indicators in a hybridizing wolf population[END_REF], admixed wolf populations are increasingly being reported where wolves live in close contact with free-ranging dogs population (e.g., Italy: [START_REF] Caniglia | Black coats in an admixed wolf dog pack is melanism an indicator of hybridization in wolves?[END_REF][START_REF] Galaverni | Disentangling timing of admixture, patterns of introgression, and phenotypic indicators in a hybridizing wolf population[END_REF][START_REF] Salvatori | Record-high levels of wolf x dog introgressive hybridization in agricultural landscapes of central Italy[END_REF].

The matrix model presented here allows quantitatively assessing the possibles outcomes of hybridization (i.e., genomic extinction vs persistence) under different fitness scenarios. Providing management recommendations is beyond the scope of this study, our model represents a valuable tool to inform management once appropriately customized and parametrized. While our focus is on mammalian species, the analytical approach described here is valid for other taxa, and it could be adopted to project the dynamics of admixed populations for situations entailing both natural and anthropogenic hybridization.

Methods

General model

We consider a system in which there are two parental groups (T1 and T2) that interbreed and produce an admixed progeny (H). We regard the admixed progeny as an absorbing state encompassing all offspring produced by pairs of different parental groups, where at least one of the parents is an admixed individual [START_REF] Wolf | Predicting the risk of extinction through hybridization[END_REF]. The possible crosses considered and the produced progeny are listed in Table 1. Based on the assumed dynamics of interbreeding and production of offspring, we project the future abundance of the three mixing groups over time. We present the projections step by step with linear equations and we introduce the equivalent matrix formulation to calculate the asymptotic growth rate and perform sensitivity analyses. Finally, we discuss the behavior of the model by applying it to our case studies.

Model equations

We assume that parental and admixed groups have similar life cycles, with three age classes:

offspring (C individuals up to 1 year old), juveniles (J non-reproductive individuals up to 3 years old) and adults (A reproductive individuals ≥ 3 years old). The transitions among age classes are described by survival parameters (S) and the reproductive parameters are described by per capita fertility rates (f). At time t, the total number of individuals of the different groups (the parentals and the admixed) is:

N TOT = N 1 + N 2 + N H,
where N 1 is the total number of individuals in group 1, N 2 is the total number of individuals in group 2 and N H is the total number of admixed individuals.

We model the dynamics of females and males separately. Below, we show the equations for females. Each group at time t is composed by females (f) and males (m) belonging to the three different age classes:

N 1f (t) = C 1f (t) + J 1f (t) + A 1f (t), N 2f (t) = C 2f (t) + J 2f (t) + A 2f (t), N Hf (t) = C Hf (t) + J Hf (t) + A Hf (t).
We assume that the probability of reproduction between adult females belonging to one group and adult males belonging to one of the three other groups, is proportional to the relative abundances of adult males α(t), β(t) and γ(t) which are given by the ratio between adult males (A m ) of one group divided by the total number of adult males TOT m (t):

α(t) = A 1m (t) / TOT m (t), β(t) = A 2m (t) / TOT m (t), γ(t) = A Hm (t) / TOT m (t).
Therefore, the number of offspring belonging to the parental group 1 produced at time t+1 is given by the number of females surviving to time t+1, A 1f (t)Sa 1 , multiplied by their per capita fertility rate (f 1 ) multiplied by the relative abundance of adult males of group 1 α(t). The number of offspring belonging to the parental group 2 produced at time t+1 is therefore:

C 2 (t+1) = A 2f (t)f 2 Sa 2 β(t).
The number of admixed offspring produced by e.g., the crossing between females of the group 1 and males of the group 2 is calculated as A 1f (t)f 1 Sa 1 β(t) and the total number of admixed offspring at time t +1 is given by the sum of the contribution of all the possible crosses (Table 1):

C H (t+1) = A 1 (t)f 1 Sa 1 β(t) + A 2 (t)f 2 Sa 2 α(t) + A 1 (t)f 1 Sa 1 γ(t) + A 2 (t)f 2 Sa 2 γ(t) + A H (t)f h Sa h α(t) + A H (t)f h Sa h β(t) + A H (t)f h Sa h γ(t),
where A 2f (t) is the number of adult females of group 2 at time t, Sa 1 and Sa 2 are group-specific adult survival values, f 1 and f 2 are the group-specific per capita fertility rates. Because we model separately males and females, the number of offspring of each sex produced every year by each group is obtained by multiplying the total number of offspring by 0.5, assuming a 50:50 sex ratio at birth:

C 1f (t+1) = C 1m (t+1) = C 1 (t+1)0.5, C 2f (t+1) = C 2m (t+1) = C 2 (t+1)0.5, C Hf (t+1) = C Hm (t+1) = C H (t+1)0.5.
Depending on the species mating system, hybridization can be modelled in different ways. For example, for species in which only the dominant individuals reproduce, it can be convenient to model hybridization at the level of the formation of the reproductive pairs. We consider this situation when dealing with the wolf x dog case study (section 3.2).

The number of female adults and juveniles at time t + 1 for the three groups are obtained as follows (the equations are showed only for group 1):

J 1f (t+1) = C 1f (t)Sc 1 A 1f (t+1) = J 1f (t)Sj 1 + A 1f (t)Sa 1
where Sc, Sj and Sa are respectively survival probabilities for offspring, juveniles and adults. The total number of females at time t + 1 is therefore:

N 1f (t+1) = A 1f (t)f 1 Sa 1 α(t)0.5 + C 1f (t)Sc 1 + J 1f (t)Sj 1 + A 1f (t)Sa 1 N 2f (t+1) = A 2f (t)f 2 Sa 2 β(t)0.5 + C 2f (t)Sc 2 + J 2f (t)Sj 2 + A 2f (t)Sa 2 N Hf (t+1) = [A 1f (t)f 1 Sa 1 β(t) + A 2f (t)f 2 Sa 2 α(t) + A 1f (t)f 1 Sa 1 γ(t) + A 2f (t)f 2 Sa 2 γ(t) + A Hf (t)f h Sa h α(t) + A Hf (t)f h Sa h β(t) + A Hf (t)f h Sa h γ(t)]0.5 + C Hf Sc h + J Hf (t)Sj h + A Hf (t)Sa h
The number of male individuals in the three groups is obtained with the following equations:

N 1m (t+1) = A 1f (t)f 1 Sa 1 α(t)0.5 + C 1m Sc 1 + J 1m (t)Sj 1 + A 1m (t)Sa 1 N 2m (t+1) = A 2f (t)f 2 Sa 2 β(t)0.5 + C 2m Sc 2 + J 2m (t)Sj 2 + A 2m (t)Sa 2 N Hm (t+1) = [A 1f (t)f 1 Sa 1 β(t) + A 2f (t)f 2 Sa 2 α(t) + A 1f (t)f 1 Sa 1 γ(t) + A 2f (t)f 2 Sa 2 γ(t) + A Hf (t)f h Sa h α(t) + A Hf (t)f h Sa h β(t) + A Hf (t)f h Sa h γ(t)]0.5 + C Hm Sc h + J Hm (t)Sj h + A Hm (t)Sa h

Matrix formulation

The model above can be conveniently formulated as N(t+1) = AN(t) where: and

In the following section, we consider two real-world scenarios illustrating the relevance of the model described above. All the analyses are performed with the software R (R core team, 2018).

Applications

Dolphin model

We built stage-based matrices (Fig. 1; [START_REF] Taylor | Generation length and percent mature estimates for IUCN assessments of cetaceans[END_REF] for the two parental species and the admixed individuals with the following stages: calf (individuals up to 1 year old), juveniles (nonreproductive individuals up to species-specific age of first reproduction reported by [START_REF] Taylor | Generation length and percent mature estimates for IUCN assessments of cetaceans[END_REF] and adults (individuals that reached the age of first reproduction). We used the available stage-specific demographic parameters for the target populations in the Gulf of Corinth, Greece [START_REF] Santostasi | A robust design capturerecapture analysis of abundance, survival and temporary emigration of three odontocete species in the Gulf of Corinth, Greece[END_REF] and we used parameters estimated for other populations of the same species for the non available stage-specific demographic parameters (see Table 2 and Appendix A for details about model parametrization). We built three fitness scenarios for admixed individuals (Table 2): i) in the Null Model we assumed that admixed individuals had intermediate demographic traits between the two species, ii) in the Hybrid Vigour scenario we assumed that hybrids have higher survival (the upper 95% confidence limit estimated for the population by [START_REF] Santostasi | A robust design capturerecapture analysis of abundance, survival and temporary emigration of three odontocete species in the Gulf of Corinth, Greece[END_REF] and annual per capita fertility rate (the highest annual pregnancy rate reported for Atlantic common dolphin subpopulations reported in [START_REF] Murphy | Importance of biological parameters in assessing the status of Delphinus delphis[END_REF] and that they become reproductively mature earlier (having the youngest age of first reproduction between the two parental species), iii) in the Outbreeding Depression scenario we assumed that admixed individuals have lower survival (the lower 95% confidence limit estimated for this population by Santostasi et al., 2016) and annual per capita fertility rate (the lowest annual pregnancy rate reported for Atlantic common dolphin subpopulations reported in [START_REF] Murphy | Importance of biological parameters in assessing the status of Delphinus delphis[END_REF], and that they become reproductivley mature later (having the oldest age of first reproduction between the two parental species). Transition probability from the juvenile to the adult stage (Tja) 0.94 (10) 0.94 (8) 0.94 (9) 0.92 (10) 0.96 (8) Adult survival (Sa) 0.94 1 0.94 1 0.94 0.92 1 0.96 1 1 Gulf of Corinth, Greece [START_REF] Santostasi | A robust design capturerecapture analysis of abundance, survival and temporary emigration of three odontocete species in the Gulf of Corinth, Greece[END_REF]; 2 Western Mediterranean (Calzada et al., 1995); 3 Atlantic [START_REF] Murphy | Importance of biological parameters in assessing the status of Delphinus delphis[END_REF]); 4 Iberian [START_REF] Murphy | Importance of biological parameters in assessing the status of Delphinus delphis[END_REF]; 5 control group [START_REF] Murphy | Importance of biological parameters in assessing the status of Delphinus delphis[END_REF]; 6 Eastern North Atlantic [START_REF] Mannocci | Assessing the impact of bycatch on dolphin populations: the case of the common dolphin in the Eastern North Atlantic[END_REF]; 7 theoretical calculation by [START_REF] Taylor | Generation length and percent mature estimates for IUCN assessments of cetaceans[END_REF].

The mating systems of odontocete cetaceans (toothed whales) have been reported as either polygynous (some males with multiple partners) or polygynandrous (both males and females with multiple partners; Murphy et al., 2007). We assumed that mature females belonging to one species would reproduce with mature males belonging to the same species, the other species or the admixed individuals, proportionally to their relative abundance. We used deterministic projections (i.e., with constant parameters) to compare the predicted time of extinction for the two species with and without the hybridization effect. The matrix formulation corresponding to the dolphin model is presented below. The subscript Sc refers to striped dolphin, the subscript Dd refers to common dolphins the subscript h refers to admixed individuals. In the Appendix A we show the corresponding R code, the linear equations, and we describe how to incorporate demographic stochasticity.

Wolf x dog model

To model wolf life cycle we used the pre-breeding stage-structured model described in [START_REF] Marescot | Reducing matrix population models with application to social animal species[END_REF], that simplifies the complex wolf social structure in four age-and social-stages: yearlings, individuals in dispersal, subordinates, and breeders (Fig. 2). Only the dominant couple breeds, producing a single litter of pups every year. The pups that survived their first year (yearlings), may survive to their second year (probability Sj) and disperse (probability Pdi), or they may remain in the natal pack (probability 1 -Pdi) as subordinates (sexually mature individuals 1 to 2 years old, inhibited from breeding by the presence of breeders; [START_REF] Zimen | Social dynamics of the wolf pack[END_REF]. The model makes the following assumptions [START_REF] Marescot | Reducing matrix population models with application to social animal species[END_REF]: 1) if subordinates survive (with probability Ss) one year in the natal pack they will then leave and become individuals in dispersal by their 3rd year of age; 2) the subordinates never directly transition to being breeders; 3) individuals in dispersal can either die or survive (with probability Sd), and gain access to reproduction (with probability Pes) or remain dispersers (with probability 1-Pes); 4) breeders never lose their status remaining in the breeder stage and surviving with probability Sa. To run the projections we used the demographic parameters that were estimated for the expanding wolf population in the Italian Alps [START_REF] Marucco | Wolf survival and population trend using non-invasive capture-recapture techniques in the Western Alps[END_REF]2010) in which hybridization has not been detected yet [START_REF] Fabbri | From the Apennines to the Alps: colonization genetics of the naturally expanding Italian wolf (Canis lupus) population[END_REF]. The details about model parametrization are shown in Appendix B. We modeled hybridization at the level of the formation of the reproductive pairs. Hybridization was documented in the almost totality of cases to occur between female wolves and male dogs [START_REF] Randi | Detecting hybridization between wild species and their domesticated relatives[END_REF][START_REF] Godinho | Genetic evidence for multiple events of hybridization between wolves and domestic dogs in the Iberian Peninsula[END_REF][START_REF] Pacheco | Spatial assessment of wolf-dog hybridization in a single breeding period[END_REF], while the opposite case of hybridization between male wolves and female dogs appears to be rare [START_REF] Hindrickson | Bucking the trend in wolf-dog hybridization between female dogs and male wolves[END_REF]. Therefore, in the model we did not consider reproductive pairs formed by male wolves and female dogs. For simplicity, we did not model the occurrence of backcrossing to dogs (the reproduction of admixed individuals and dogs; Table 3). We assumed that a constant number of reproductive events happen every year between female wolves and male dogs (parameter frd). The probability of the formation of a reproductive pair between a female wolf and male wolf at time t +1 is therefore expressed as the product between the abundance of wolf females in dispersal (Df w ) by their survival (Sd w ) by the probability that a female wolf establishes a pack (Pes w ) by the relative abundance of wolf males in dispersal α(t). The probability of formation of wolf-dog reproductive pairs, and of reproductive pairs in which at least one individual is admixed, is formulated in Table 3. For wolves, we used the probability of establishing a pack calculated by [START_REF] Marescot | Reducing matrix population models with application to social animal species[END_REF].

Such annual probability varies uniformly between 0.3 to 0.7. For the Null Model, we assigned the the same average value (Pes=0.5) to wolves and admixed individuals. We produced two alternative scenarios of reproductive isolation by simulating a reduced probability of establishing a pack for admixed females. In one scenario (Reproductive Isolation 1), we attributed the average value (0.5) to wolves and the minimum value (0.3) to the admixed individuals. In the second scenario (Reproductive Isolation 2) we assigned the maximum probability of establishing a pack to wolves (0.7) and the minimum probability (0.3) to the admixed individuals (Table 4). We used the prevalence of hybrids as a measure to define the type of hybridization reached at the end of the projection time-frame. The matrix formulation corresponding to the wolf x dog model is presented below.

In the Appendix B we show the corresponding R code and the linear equations.

Sensitivity analysis

We performed a sensitiviy analysis to explore the relationship between the hybridization outcome (genomic extinction vs persistence of parental and admixed populations after 50 years) and the demographic parameters of parental species and admixed individuals (Mc Carthy et al., 1995;[START_REF] Cross | Using logistic regression to analyze the sensitivity of PVA models: a comparison of methods based on African wild dog models[END_REF]. Here, we use "genomic extinction" to refer to the disappearence of parental individuals from the population [START_REF] Allendorf | The problems with hybrids: setting conservation guidelines[END_REF]. We generated 15,000 parameter sets by drawing them from uniform distributions (Table 5) to emphasize the effects of variability in vital rates on model sensitivity [START_REF] Cross | Using logistic regression to analyze the sensitivity of PVA models: a comparison of methods based on African wild dog models[END_REF]. We projected population abundance with each parameter set and cheked if the population went extinct or not using a quasi-extinction threshold of 5 parental mature females. We conducted logistic regressions to explore the relationship between the probability of genomic extinction of the two parental species as response variables and demographic parameters used in the projections as independent explanatory variables.

For each regression we built full models including all the parameters and used a backward stepwise selection procedures to identify the sets of parameters of potential importance [START_REF] Fredrickson | Dynamics of hybridization and introgression in red wolves and coyotes[END_REF]. From those sets of parameters we ranked their relative importance in affecting the probability of genomic extinction of the parental species based on their standardized regression coefficients which are the regression coefficient divided by their standard error [START_REF] Cross | Using logistic regression to analyze the sensitivity of PVA models: a comparison of methods based on African wild dog models[END_REF]). In the Appendix A and B we show the code to peform the logistic regression sensitivity analysis and how to perform a sensitivity analysis of the asymptotic growth rate using the matrix formulation for the two case studies (R codes are also reported there).

Table 5. Demographic parameters range used for drawing the demographic parameters from uniform distributions to project the abundance of parental and admixed individuals. These parameters are also used as explanatory variables in the logistic regression sensitivity analysis.

Dog

Frequency of mating with dogs 0-1 1 Gulf of Corinth, Greece [START_REF] Santostasi | A robust design capturerecapture analysis of abundance, survival and temporary emigration of three odontocete species in the Gulf of Corinth, Greece[END_REF]; 2 Western Mediterranean (Calzada et al., 1995); 3 control group [START_REF] Murphy | Importance of biological parameters in assessing the status of Delphinus delphis[END_REF]; 4 Eastern North Atlantic [START_REF] Mannocci | Assessing the impact of bycatch on dolphin populations: the case of the common dolphin in the Eastern North Atlantic[END_REF]; 5 Mech and Boitani, 2003; 6 Marescot et al., 2012 Results

Striped x common dolphin model

Comparing the projections with and without hybridization, we see that when hybridization is ignored, the model predicts an exponential population growth (Fig. 3). When the impact of hybridization is taken into account, the least abundant species (the common dolphin) reaches the quasi-extinction threshold after 16 years, no matter the hybrid fitness scenario. In the Null Model and in the Hybrid Vigour scenario, the growth rate of the most abundant species (the striped dolphin) is also affected, becoming negative after about 40 years in the first scenario and after about 30 years in the second. In the absence of mechanisms that counter hybridization (e.g., assortative mating), the population will eventually be composed of a continuum of hybrid classes (the "complete admixture" hybridization type described by [START_REF] Allendorf | The problems with hybrids: setting conservation guidelines[END_REF]Fig. 4). Conversely, under the Outbreeding Depression scenario, striped dolphin population abundance is expected to increase in the next 100 years, predicting a "widespread introgression" hybridization type (i.e., the coexistence of admixed and parental individuals; Fig. 4).

Wolf x dog model

Depending on the scenario, hybridization had different final outcomes. In the Null Model all the scenarios reached complete admixture (prevalence of hybrids = 1) after 50 years and prevalence increased more rapidly at the increasing of the recurrent gene flow from dogs (Fig. 5, upper-left panel). The reproductive isolation scenarios showed different outcomes. In the scenario with the weakest reproductive isolation (Reproductive Isolation 1) the prevalence increased less rapidly compared to the Null Model (Fig. 5, upper-right panel) but did not reach an asymptote, heading towards a final outcome of complete admixture [START_REF] Allendorf | The problems with hybrids: setting conservation guidelines[END_REF]. In the scenario of strongest reproductive isolation (Reproductive Isolation 2), prevalence reached an asymptote whose final value increased at the increasing of the intensity of the recurrent gene flow from dogs (Fig. 5, lower-left panel). In this last case the final outcome is the co-existence of admixed and parental individuals with a constant prevalence of hybrids (assuming that all the demographic parameters remain constant over time). This last scenario falls into the definition of "widespread introgression" [START_REF] Allendorf | The problems with hybrids: setting conservation guidelines[END_REF].

Sensitivity analysis

For the dolphin case study the probability of genomic extinction of both striped and common dolphins was most affected by the survival of mature individuals followed by juvenile survival with almost the same relative importance (Table 6). However, for common dolphins, the initial abundance was also important (Table 6).

Table 6. Logistic regression sensitivity analysis results for striped and common dolphin population projections in the presence of hybridization. The parameters are sorted in descending order of relative importance, evaluated by looking at the standardized regression coefficients. Sa = adult survival, Sj = juvenile survival, Sc = calf survival, f = annual per captita fertility rate, ia = initial abundance. The subscript Sc refers to striped dolphin parameters, the subscript Dd refers to common dolphin parameters, the subscript h refers to admixed individuals parameters. For the wolf x dog case study, our sensitivity analysis showed that wolf genomic extinction probability was mostly affected by wolf breeder survival (Sa w ), followed by annual frequency of mating with dogs (frd), wolf access to reproduction (i.e., probability of establishing a pack, pes w )

and annual wolf per capita fertility rate f w , with the same relative importance (Table 7).

Discussion

Identifying the demographic factors affecting the outcome of hybridization helps both understanding evolutionary mechanisms and developing meaningful management and conservation measures when hybridization is a threat [START_REF] Bohling | Strategies to address the conservation threats posed by hybridization and genetic introgression[END_REF]. Compared to the genetic approach [START_REF] Huxel | Rapid displacement of native species by invasive species: effects of hybridization[END_REF][START_REF] Ferdy | Extinction and introgression in a community of partially cross-fertile plant species[END_REF], our model has the advantage of making the link between hybridization and both demographic (e.g., demographic stochasticity) and ecological factors (e.g., environmental stochasticity). Moreover, our approach allows to relate the viability of the population to the contribution of different fitness components (such as survival and reproductive rates), ultimately providing the basis for sensitivity analyses.

On the other hand, ignoring genetic processes can be an important limitation [START_REF] Hall | What can mathematical modeling tell us about hybrid invasions?[END_REF]. In its present formulation, our model does not make a distinction between ancient (later generations backcrosses) and recent hybrids (first and second generation hybrids), which raises two issues. First, because later generations backcrosses have a limited mixed genomic content, in several cases they are indistiguishable from parentals with the current diagnostic techniques [START_REF] Vähä | Efficiency of model-based Bayesian methods for detecting hybrid individuals under different hybridization scenarios and with different numbers of loci[END_REF]. Second, later generations backcrosses are often considered as parentals from a management perspective. For example, in the management of the anthropogenic hybridization between bontebok Damaliscus pygargus pygargus and blesbok D. p. phillipsi in South Africa, van Wick et al., (2016) considered that backrosses to parental bonteboks with ≥ 90% of bontebok genes are not to be removed from the population. In the management of the red wolf and coyote, canids with > 87.5 % of red wolf genes are classified as red wolves and not sterilized or culled [START_REF] Gese | Managing hybridization of a recovering endangered species: the red wolf Canis rufus as a case study[END_REF]. To overcome this limitation, it may be desirable to split the "admixed" group into recent and ancient hybrids, and pool the ancient hybrids together with the parental populations. The threshold between the two groups could be set according to the defined acceptable threshold of admixture in the parentals, while also considering the power of the diagnostic system [START_REF] Vähä | Efficiency of model-based Bayesian methods for detecting hybrid individuals under different hybridization scenarios and with different numbers of loci[END_REF]van Wick et al., 2016). Despite these limitations, our modelling framework provides a clear way to make explicit the assumptions underlying the definition of hybrids [START_REF] Thompson | Model-based analysis of the likelihood of gene introgression from genetically modified crops into wild relatives[END_REF].

Even though population-based approaches like the one described here require less-detailed information compared to individual-based approaches [START_REF] Hall | What can mathematical modeling tell us about hybrid invasions?[END_REF], one still needs to know the survival and reproductive rates for all the population stages, and these may be unavailable for the target population. In this study, we overcome such limitation by resorting to parameters estimated for other populations, which may introduce an unknown level of uncertainty. Moreover, estimates of demographic rates of admixed individuals were not available and are likely to be unavailable in many studies. In such cases, one would need to make assumptions regarding admixed individuals survival and fertility rates. However, sensitivity analysis showed that genomic extinction probability was largely affected by the demographic parameters of parentals, whereas the parameters of admixed individuals had a lower relative importance. Nonetheless, the sensitivity of model predictions to unknown parameters should be tested through sensitivity analyses, and should be taken into account by considering different scenarios. Moreover, the comparison alternative simulated scenarios with empirical data (i.e the observed prevalence trends) can provide insight on the evolutionary mechanism shaping the observed patterns.

Another important limitation is that, for simplicity and for the absence of species-specific information, we assumed random mating and did not model more complex mating choice scenarios that could lead to assortative mating. If information about assortative mate choice is available, it can be modelled by changing the reproduction coefficients (α, β and γ in paragraph 2.1.1) from being proportional to the relative abundance of mature males to reflecting the observed percentage of conspecific and heterospecific matings.

At the moment, our projections cannot be validated due to the absence of long-term data on the quantitative evaluation of hybridization for the two case-studies. However, the predicted outcomes are consistent with previous studies showing that hybridization can lead to rapid genomic extinction in the absence of reproductive isolation [START_REF] Huxel | Rapid displacement of native species by invasive species: effects of hybridization[END_REF] or other forms of segregation [START_REF] Wolf | Predicting the risk of extinction through hybridization[END_REF][START_REF] Fredrickson | Dynamics of hybridization and introgression in red wolves and coyotes[END_REF]. Two studies on wolf x dog hybridization [START_REF] Santostasi | Use of hidden Markov capture-recapture models to estimate abundance in the presence of uncertainty: application to the estimation of prevalence of hybrids in animal populations[END_REF]Salvatori et al., 2019) showed that the prevalence of hybrids can locally reach high levels (around 50%) where reproductive isolation between wolves and dogs is disrupted by anthropogenic disturbance (e.g., presence of free ranging dogs, food provisioning, high wolf anthropogenic mortality), supporting the results of our projections. We suggest that, given the importance of reproductive isolation mechanisms in determining the wolf x dog hybridization outcome, more research should be addressed to understand to what extent and in which circumstances those mehcanisms are effective in contrasting wolf x dog hybridization. In absence of such information, the absence of strong reproductve isolation should be assumed as a precautionary measure.

Studies on hybridization dynamics in cetaceans that could be used to validate our results are absent.

However, our finding that the genomic extinction risk for the parental population increases as their initial frequency decreases, making the least abundant species particularly vulnerable, is consistent with previous studies [START_REF] Allendorf | The problems with hybrids: setting conservation guidelines[END_REF]Epifanio and Philipp, 2001).

We stress that the interest of our model is not the production of absolute predictions of population fate. Being long term hybridization dynamics otherwise difficult to test empirically, the practical value of our model is to illustrate the expected relative outcomes of alternative biological and managament scenarios. The comparison of those outcomes can be useful to hihglight future research priorities and to inform decision-making in a context of uncertainty (Gervasi and Ciucci 2018).

In alternative, not acknowledging hybridization in the projections may lead to underestimating the risk of genomic extinction. The common dolphin (the least abundant species) had a 100% probability of going extinct after a relatively short time (16 years) regardless the fitness scenario. A previous count-based projection [START_REF] Santostasi | Common dolphins in the Gulf of Corinth are Critically Endangered[END_REF] did not include the effect of hybridization (although it included the effect of demographic stochasticity) and estimated a considerably lower (50%) probability of extinction after 15 years. Although the two predictions are not entirely comparable due to the different methodologies, our results suggest that hybridization may pose an additional and important threat that must be considered when evaluating the conservation status of common dolphins. This is particulatly relevant at larger scales, considering that, in the Mediterranean Sea common dolphins occurr in simpatry with the more abundant striped dolphins and mixed-species groups are not rare [START_REF] Frantzis | Mixed-species associations of striped dolphins (Stenella coeruleoalba), short-beaked common dolphins (Delphinus delphis), and Risso's dolphins (Grampus griseus) in the Gulf of Corinth (Greece , Mediterranean Sea)[END_REF][START_REF] Giménez | Intra-and interspecific niche partitioning in striped and common dolphins inhabiting the southwestern Mediterranean Sea[END_REF][START_REF] Espada | Hybridization in the wild between Tursiops truncatus (Montagu 1821) and Delphinus delphis[END_REF].

Interestingly, the population growth rate of the most abundant species, the striped dolphin, was also affected and hybridization could eventually lead this population to genomic extinction, under scenarios in which the fitness of admixed individuals is equal to or greater than that of parentals (i.e., the Null Model and the Outbreeding Depression scenario). However, striped dolphins are predicted to persist and even increase in the next 100 years, if the fitness of admixed individuals is lower than that of parentals (i.e., the Outbreeding Depression scenario). More research on the fitness of admixed individuals would help understand which scenario represents the most likely alternative. Research on the genetic composition and ecology of admixed individuals would help understand if the predicted "hybrid swarm" may: 1) lead to a new species (Lamichhaneey et al., 2017;[START_REF] Larsen | Natural hybridization generates mammalian lineage with species characteristics[END_REF]; 2) preserve the integrity of the striped dolphin species by backcrossing and dilution of the common dolphin genetic contribution; or 3) allow for the persistence of adaptive common dolphin alleles in the striped dolphin, as suggested by other cases of adaptive introgression (Figueirò et al., 2017).

For the wolf x dog case study, our results confirm that hybridization with dogs should be expected to be a serious threat to the wolf genomic integrity, at least under weak reproductive isolation and/or frequent breeding with dogs. The sensitivity analysis pointed out that the most influential parameters decreasing the chances of genomic extinction are linked to social and reproductive integrity (i.e., survival of wolf breeders and the per capita fertility rate) and the annual frequency of mating with dogs. These results are in agreement with observations of hybridization in eastern wolves (Canis lycaon) and coyotes in Canada (Rutledge et al., 2011), and of red wolves and coyotes in North Carolina [START_REF] Bohling | Factors influencing red wolf-coyote hybridization in eastern North Carolina, USA[END_REF]. In both cases, high levels of mortality coupled with a large availability of coyotes have been identified as the main causes of hybridization. Our results offer further evidence that human-related factors contributing to hybridization (i.e., increased wolf mortality and dog presence) must be managed to avoid the risk of genomic extinction of wild wolves (Rutledge et al., 2011;[START_REF] Bohling | Factors influencing red wolf-coyote hybridization in eastern North Carolina, USA[END_REF].

In conclusion, our model provides estimates of genomic extinction risk in presence of hybridization by using data obtained during demographic monitoring programs of threatened populations. This model allows to link demographic parameters and environmental variables, therefore predicting hybridization dynamics under changing environments. Population projection models clearly represent a valuable tool to predict the outcome of hybridization, therefore contributing to management decision [START_REF] Kelly | The Arctic melting pot[END_REF].

Figure 1 .

 1 Figure 1. Life cycle used for the dolphin case study. Parameters are: Sc = survival of calves, Tja = transition probability from the juvenile to the adult stage modelled as the survival of juveniles (Sj) raised power of age of first reproduction -1, Sa = survival of adults, f = per capita fertility rate, approximated as the annual pregnancy rate.

Figure 2 .

 2 Figure 2. Life cycle for the wolf case study. Parameters are: Sj = survival of yearlings, Ss = survival of subordinates (yearlings and adults), Sd = survival of dispersers (yearlings and adults), Sa = survival of breeders, f = annual per capita fertility rate approximated as the litter size, Pdi = probability of going in dispersal, Pes = probability of establishing a pack.

Figure 3 .

 3 Figure 3. Population trajectories of striped, common and admixed dolphins under different parental and admixed individuals fitness scenarios (the continuos line represent the Null Model while different type of dashed lines represent the model in absence of hybridization, the Outbreeding Depression Scenario and the Hybrid Vigour Scenario).

Figure 4 .

 4 Figure 4. Projected prevalence of striped x common dolphin hybrids in the mixed dolphin population of the Gulf of Corinth, Greece under different parental and admixed individuals fitness scenarios (the continuos line represent the Null Model while different type of dashed lines represent the Outbreeding Depression Scenario and the Hybrid Vigour Scenario).

Figure 5 .

 5 Figure 5. Projected prevalence of wolf x dog hybrids under increasing levels of recurrent gene flow from dogs (parameter frd) and under different reproductive isolation scenarios: the Null Model (upper-left panel), the Reproductive Isolation 1 scenario (upper-right panel). and the Reproductive Isolation 2 scenario (lower-left panel).

Table 1 .

 1 Possible crosses in the system formed by two parental taxa (T1 and T2) and admixed individuals (H).

	Female parent	Male parent	Offspring
	T1	T1	T1
	T2	T2	T2
	T1	T2	H
	T2	T1	H
	H	T1	H
	H	T2	H
	T1	H	H
	T2	H	H
	H	H	H

Table 2 .

 2 Demographic parameters used for projecting the abundance of striped and common dolphins in the presence of hybridization using alternative parental and admixed individual fitness scenarios.

	Scenario	All scenarios	Null	Outbreeding	Hybrid Vigour
	Individuals	Striped dolphin	Common dolphin	Admixed dolphin	Admixed dolphin	Admixed dolphin
	Initial abundance	1331 1	22 1	55 1	55 1	55 1
	Per capita fertility rate (f)	0.25 2	0.26 3	0.26	0.19 4	0.33 5
	Age of first reproduction	11 2	9 3,6	10	11	9
	Calf survival (Sc)	0.80 7	0.80 7	0.80	0.80 7	0.80 7
	Juvenile survival (Sj)	0.94 1	0.94 1	0.94	0.92 1	0.96 1

Table 3 .

 3 Probability of formation of the reproductive pairs considered in the model. Df w = abundance of wolf females in dispersal, Pes w = probability that a female wolf establishes a pack, Pes

	Female parent	Male parent	Probability
	W	W	Df w sd w pes w α(t)
	W	D	Df w sd w pes h frd(t)
	W	H	Df w sd w pes h γ(t)
	H	W	Df w Sd h pes h γ(t)
	H	H	Df w Sd h pes h γ(t)

h = probability that an admixed female establishes a pack, frd = constant annual frequency of reproductive events between female wolves and male dogs, α(t) = relative abundance of male wolves in dispersal, γ(t) = relative abundance of admixed males in dispersal.

Table 4 .

 4 Demographic parameters used for projecting the abundance of wolves and admixed individuals in the presence of hybridization.Marucco et al., 2010; 2 Blanco and Cortés 2007; 3 Marucco et al., 2009; 4 Mech and Boitani, 2003; 5 Marescot et al., 2012.

	Scenario	Null Model	Repro. Isolation 1	Repro. Isolation 2
	Individual	Wolf	Admixed	Wolf	Admixed	Wolf	Admixed
	Initial abundance
	(reproductive	6 1	0	6 1	0	6 1	0
	pairs)
	Average litter size (f)	3.387 1	3.387 1	3.387 1	3.387 1	3.387 1	3.387 1
	Juvenile survival (Sj)	0.551 1	0.551 1	0.551 1	0.551	0.551 1	0.551
	Subordinate
	survival	0.82 1	0.82 1	0.82 1	0.82 1	0.82 1	0.82 1
	(Ss)
	Disperser
	Survival	0.69 2	0.69 2	0.69 2	0.69 2	0.69 2	0.69 2
	(Sd)
	Breeder survival (Sa)	0.82 3	0.82 3	0.82 3	0.82 3	0.82 3	0.82 3
	Probability of
	dispersing	0.25 4	0.25 4	0.25 4	0.25 4	0.25 4	0.25 4
	(Pdi)
	Probability of
	establishing a	0.5 4,5	0.5 4,5	0.5 4,5	0.3 4,5	0.7 4,5	0.3 4,5
	pack
	(Pes)

Table 7 .

 7 Logistic regression sensitivity analysis results for the wolf population projections in the

	presence of hybridization with dogs. The parameters are sorted in descending order of relative
	importance, evaluated by looking at the standardized regression coefficients. Sa = annual breeder
	survival, f = fertility (approximated as litter size), Pes (probability of establishing a pack), frd =
	annual frequency of mating with dogs, Sj = annual juvenile survival, Sd = annual disperser survival,
	Ss = annual subordinate survival, pdi = annual probability of dispersal. The subscript w indicates
	wolf demographic parameters and the subscript h indicates admixed demographic parameters.
	Coefficient	Standard Error	Standardized Coefficient	p value

Acknowledgements

OG and NLS were supported by the French National Research Agency with a Grant ANR-16-CE02-0007. NLS was supported by a PhD grant by Sapienza University of Rome. We thank the Editor and two anonimous reviewers for their insightful suggestions during the editorial process. We thank Dr. Federico Manna for his help.