Javier R Goicoechea

Olivier Berné

Maryvonne Gerin

Christine Joblin

David Teyssier

STAR FORMATION NEAR PHOTODISSOCIATION REGIONS: DETECTION OF A PECULIAR PROTOSTAR NEAR CED 201 1,2,3

Keywords: ISM : individual (IRAS 22129+7000, Ced 201, B175/L1219) -ISM : jets and outflows -stars: formation -stars: pre-main sequence

We present the detection and characterization of a peculiar low-mass protostar (IRAS 22129+7000) located ∼0.4 pc from Cederblad 201 Photodissociation Region (PDR) and ∼0.2 pc from the HH450 jet. The cold circumstellar envelope surrounding the object has been mapped through its 1.2 mm dust continuum emission with IRAM-30m/MAMBO. The deeply embedded protostar is clearly detected with Spitzer /MIPS (70 µm), IRS (20-35 µm) and IRAC (4.5, 5.8, and 8 µm) but also in the K s band (2.15 µm). Given the large near-and mid-IR excess in its spectral energy distribution, but large submillimeter-to-bolometric luminosity ratio (≃2%), IRAS 22129+7000 must be a transition Class 0/I source and/or a multiple stellar system. Targeted observations of several molecular lines from CO, 13 CO, C 18 O, HCO + and DCO + have been obtained. The presence of a collimated molecular outflow mapped with the CSO telescope in the CO J=3-2 line suggests that the protostar/disk system is still accreting material from its natal envelope. Indeed, optically thick line profiles from high density tracers such as HCO + J=1-0 show a red-shifted-absorption asymmetry reminiscent of inward motions. We construct a preliminary physical model of the circumstellar envelope (including radial density and temperature gradients, velocity field and turbulence) that reproduces the observed line profiles and estimates the ionization fraction. The presence of both mechanical and (non-ionizing) FUV-radiative input makes the region an interesting case to study triggered star formation.

INTRODUCTION

Low-mass stars (sunlike; M ⋆ <8 M ⊙ ) form from the gravitational collapse of dense molecular cores within giant molecular clouds. During the earliest stages of evolution, embryonic protostars are deeply embedded in cold and dusty envelopes of infalling material. In recent years, large millimeter (mm) and submm continuum surveys have been carried out to detect the cold dust emission associated with these early birth-sites of star formation (André et al. 1993;[START_REF] Johnstone | [END_REF]. During collapse, conservation of angular momentum combined with infall along the magnetic field lines leads to the formation of a rotating protoplanetary disk that drives the accretion process. At the same time, both mass and angular momentum are removed from the system by the onset of jets, collimated flows and the magnetic braking action (Bontemps et al. 1996;Cabrit et al. 2007). The resulting molecular outflow starts to erode and sweep up part of the natal envelope, contributing for the clearing of the circumstellar material and the termination of the infall phase [START_REF] André | Protostars and Planets IV[END_REF]. The so-called Class 0 sources are deeply embedded protostars at their main accretion phase and represent the earliest type of young stellar objects (YSOs). Observationally they should show: (i) Evidence of a central YSO (as revealed by the presence of a collimated outflow or an internal heating source); (ii) Centrally peaked but extended massive envelope (as traced by mm and submm dust continuum emission); (iii) High submm (λ >350 µm) to bolometric luminosity ratio (L smm /L bol >0.5%), and a spectral energy distribution (SED) similar to a single graybody at 15-30 K [START_REF] André | Protostars and Planets IV[END_REF]. On the other hand, Class I sources correspond to the late accretion phase (roughly half of the initial envelope mass is accreted in the Class 0 stage), they don't satisfy (iii) and can easily be detected at near-IR wavelengths [START_REF] Lada | The Origin of Stars and Planetary Systems[END_REF].

The detection of YSOs that do not exactly match the above canonical categories can provide the missing pieces of our current understanding of early disk evolution and early dissipation of the circumstellar envelope. These sources can thus potentially help to address, for example, how and when the velocity field around YSOs changes from infall to rotation dominated, which fraction of the initial mass envelope is dissipated and not accreted onto the star/disk, as well as the timescales for such dissipation (Brinch et al. 2007;Jørgensen et al. 2007) In this paper we present the detection of a lowmass protostar in the Bok globule B175 (L1219) within the Cepheus Flare, a nearby (∼400 pc) star forming region ∼85 pc above the galactic plane (Kun 1998;Nikolić & Kun 2004).

This cloud hosts the enigmatic Cederblad 201 reflection nebula (Wolf 1908; White contours are displayed from 22 (2σ) to 55 mJy/15 ′′ -beam in steps of 11 mJy/15 ′′ -beam (1σ). The CO J=3-2 integrated intensity (from 5 to 9 K km s -1 in steps of 1 K km s -1 ) from the outflow is shown in red (from v=-4 to -2 km s -1 ) and blue (from v=-6 to -8 km s -1 ). The integrated CO J=3-2 line intensity from Ced 201 PDR (v=-5.5 to -4.5 km s -1 ) is shown in black contours (from 6 to 20 K km s -1 in steps of 1 K km s -1 ). Left panels: 12 CO J=2-1 and 4-3, HCO + J=1-0 and DCO + J=3-2 lines towards IRAS 22129+7000 position. Right panels: Declination raster in the HCO + J=1-0 line. The vertical dashed line marks the rest velocity of the cloud (-5 km s -1 ). Cesarsky et al. 2000), probably the result of a chance encounter with a B9.5V runaway star (BD+69 • 1231, T ef f ≃10,000 K), which is moving through the cloud at ∼12 km s -1 (Witt et al. 1987). The star is thus not a only a source of dissociating photons (the FUV radiation field is ∼200 times the mean interstellar radiation field) but also a source of shocks and turbulence in the region. This interaction enables us to probe matter in a different situation than in classical PDRs in which the distance from the exciting star to the illuminated regions is large and/or in which the star was formed locally. In addition, the extended cloud hosts a Herbig-Haro jet (HH450), and several arsec-scale filaments of the supernova remnant (G110.3+11.3) that are rapidly approaching to the region [START_REF] Bally | [END_REF]. The presence of both mechanical (from shocks and outflows) and radiative input (from FUV photons) makes the region an interesting dynamics laboratory to study how star formation is proceeding through the cloud.

OBSERVATIONS

A 1.2 mm continuum emission map of the region was obtained at the IRAM-30m (Pico Veleta, Spain) in March 2005 using the 117-channel MAMBO bolometer array.

The angular resolution is ∼11 ′′ . A fast-mapping mode was used to map the region [START_REF] Teyssier | A Fast-Mapping Method for Bolometer Arrary Observations[END_REF]. The total integration time was 1 hour, achieving a rms noise of 7 mJy/11 ′′ -beam. Sky noise substraction and data analysis were carried out with the MOPSIC software [START_REF] Zylka | Pocket Cookbook for MOPSIC Software[END_REF]. A 31 ′′ ×22 ′′ fullwidth at half-maximum (FWHM) condensation with a peak intensity of 60±7 mJy/11 ′′ -beam is detected ∼225 ′′ East of Ced 201 PDR. The integrated flux is 249±50 mJy, which reflects the estimated absolute calibration uncertainty. The resulting 1.2 mm emission map is shown in Figure 1 (smoothed to a beam resolution of 15 ′′ ).

The HCO + J=1-0 line was observed also at the -30m telescope during September 2005 using the B100 single sideband receiver with a channel resolution of 40 kHz. The main beam efficiency at 89 GHz is 0.82 and the and2) was made "on-the-fly". The DCO + J=3-2 line was mapped in a ∼40 ′′ ×30 ′′ region around the dust condensation, whereas pointed observations in the CO J=2-1 and 4-3 lines were taken towards the dust emission peak (Fig. 1).

The spectra were analyzed with a 1024 channel acoustooptic spectrometer with a total bandwidth of 50 MHz and a resolution of ∼100 kHz. The main beam efficiencies of the telescope were 0.70, 0.75 and 0.53 at 230, 345, and 461 GHz respectively. The angular resolution is ∼30 ′′ at CO J=2-1, ∼20 ′′ at CO J=3-2 and ∼15 ′′ at CO J=4-3. Data was processed with GILDAS. The bright reflexion nebula (the PDR) is clearly seen at the west of the region. A well collimated outflow (as traced by CO J=3-2 high velocity line wing emission) emerges from the protostar position. Contours show the integrated intensity (from 5 to 9 K km s -1 in steps of 1 K km s -1 ) in the red (from v=-4 to -2 km s -1 ) and blue (from v=-6 to -8 km s -1 ) line wings respectively.

The Ced 201 region was also observed with Spitzer as part of our SPECPDR program4 [START_REF] Joblin | Astrochemistry: Recent Successes and Current Challenges[END_REF] with the Infrared Array Camera (IRAC) at 4.5, 5.8, and 8 µm and with the Infrared Spectrograph (IRS). IRS data were obtained in the spectral mapping mode with the LL1 module (spectral range: 20 to 35 µm). The spectral cube of the source was assembled using the CUBISM software (Smith et al. 2007) from the Basic Calibrated Data (BCD) files. The 20 to 35 µm data was obtained by integrating the cube spatially (see Berné et al. 2007). The 70 µm Multiband Imaging Photometer for Spitzer (MIPS) data that are part of the PID 30790 program (Fazio et al.) were retrieved from the Spitzer archive. MIPS and IRAC post-BCD calibrated images were used directly without additional processing.

RESULTS: DETECTION OF A YOUNG PROTOSTAR

Figure 2 shows the detection of a compact continuum source at ∼8 µm (warm dust) located ∼0.4 pc East of Ced 201 PDR (Berné et al. 2007). The compact object is also detected at 4.5, 5.8, 20-35 and 70 µm with Spitzer. The near-IR source position agrees with the continuum peak of the 1.2 mm extended emission (cold dust) detected in a larger field-of-view (Fig. 1). The best coordinates inferred from IRAC images are α 2000 = 22 h 14 m 08.3 s , δ 2000 = +70 • 15 ′ 06.7 ′′ , which coincides with the IRAS 22129+7000 source, clearly a YSO. A well collimated molecular outflow emanating from the object is seen along a SE-NW axis in the CO J=3-2 line (Figs. 2 and1). This outflow was previously detected in CO J=1-0 at lower angular resolution by Nikolić & Kun (2004), who first suggested that IRAS 22129+7000 could be the driving source. We now detect the nearand mid-IR scattered light from the embedded protostar and the mm thermal emission from the surrounding dusty cocoon. In the following we examine the nature of this source and whether or not the presence of a PDR plays a role in the star formation process in the region.

Peculiar SED Characteristics

In order to investigate the emitted power in the different wavelength domains, the SED of IRAS 22129+7000 has been built with our MAMBO and Spitzer observations and complemented with previous 2MASS

(2.15 µm), IRAS (12, 25, 60 and 100 µm) and KAO (200 µm; Casey 1991) detections (Table 1). We fitted the SED with two graybodies, B λ (T ) (1 -e -τ λ ) Ω, where Ω is the solid angle subtended by the emitting region. In the far-IR and mm domain τ λ ∝ κ 1200 (1200/λ) β and T can be associated with the dust opacity and color temperature (≃T d ) of the extended envelope. Here we take a dust opacity (per gas+dust mass column density) of κ 1200 =0.01 cm 2 g -1 at 1200 µm, the usual value adopted for the dust in protostellar envelopes [START_REF] André | Protostars and Planets IV[END_REF]. Because of the implicitly assumed grain growth and ice mantle formation (Ossenkopf & Henning 1994), this value is a factor ∼5 larger than κ 1200 for standard interstellar grains. Ω is fixed as the FWHM size of the envelope inferred from the 1.2 mm continuum map, while the column density of material is varied between A V ≃5 and 30. Satisfactory fits are obtained for T d =24-19 K and dust spectral indexes β=1.4-1.8 respectively (see the best fit as continuous blue curve in Fig. 3). A luminosity of L bol =5.5±0.5 L ⊙ and a "bolometric temperature" of T bol ≃175 K (the temperature of a blackbody with the same SED mean frequency; Myers & Ladd 1993) are inferred from the full fit.

In a first inspection, the resulting SED (Fig. 3) does not look like a typical Class 0 source in the near-and mid-IR (see e.g., Whitney et al. 2003). The improved sensitivity of the IR instrumentation compared to that available in the 90's has allowed the detection of scattered light from several Class 0 sources below ∼10 µm [START_REF] Tobin | [END_REF]). In addition, the K s -band emission due to scattered light observed through a cavity opened by the outflow (with a favorable inclination angle) and/or H 2 rovibrational emission from the outflow base, has also been detected in a few Class 0 sources (e.g., Tachihara et al. 2007). IRAS 22129+7000 shows these observable characteristics but it also displays an intriguing strong, and rather flat, near-and mid-IR emission (as shown by the relatively high T bol ). In fact, Spitzer fluxes towards IRAS 22129+7000 are consistent with the typical values observed in Class I sources with the same luminosity (Furlan et al. 2008), but an order of magnitude higher than in Class 0 sources [START_REF] Tobin | [END_REF]). On the other hand, IRAS 22129+7000 shows a very large fraction of submm (λ >350 µm) to bolometric luminosity (≃2%; a signature of younger Class 0 sources still surrounded by massive dusty envelopes). 1991) and 1200 µm (MAMBO). The Spitzer/IRS fluxes towards the source agree with those measured by IRAS at 25µm probing the compact nature of the mid-IR emission. The dashed curve is a simple fit using two graybodies with different temperatures. The continuous blue curve represents the circumstellar envelope with T d =22 K and a spectral index of 1.6. The magenta pointed-and pointed-dashed curves show the best 2D radiative transfer models obtained with the Robitaille et al. (2007) fitter for two different inclinations (i = 0 corresponds to "pole-on" and i = 90 to "edgeon" configurations).

In an attempt to better understand the peculiar IR shape we tried to fit the SED with the 2D radiative transfer tool by Robitaille et al. (2007). The best models fitting the data are shown in Fig. 3 for two different inclination angles. While the long-wavelength emission from the circumstellar envelope is well reproduced by the models (and is almost not affected by the viewing angle), the short-wavelength emission from the protostellar/disk system is clearly underestimated and strongly depends on inclination. Therefore, a small viewing angle favors the detection of near-and mid-IR scattered light from the embedded protostar through the cavities cleared by the outflow. However, the blue-and red-shifted lobes of the CO outflow (Fig. 2) appear relatively well separated, and thus the source can not be completely "polein" (Cabrit & Bertout 1986). In summary, the observed SED is peculiar, with a near-and mid-IR excess that can not be fully explained in terms of inclination effects.

Circumstellar Envelope: Radial Analysis

The cold dust surrounding the YSO has been mapped through its 1.2 mm emission (Fig. 1). As mm continuum emission is generally optically thin, the measured fluxes are proportional to the temperature-weighted mass content assuming that dust grains properties are known. In order to trace the mass distribution in the envelope, the radial intensity profile I(r) has been analyzed by averaging the 1.2 mm fluxes within equidistant circular annulus (Fig. 4). An envelope radius of R out ≃15,000 AU is inferred. Assuming spherical symmetry, the intensity profile I(r) ∝ r -m contains information on the density distribution, ρ(r) ∝ r -p , and temperature distribution, T d (r) ∝ r -q , with p = m + 1 -q (e.g., Adams 1991).

Before convolving the expected intensity profile with the bolometer angular resolution, we divide the envelope in two regions (inner and outer ) and try to fit the

L bol T bol Menv Lsmm/L bol Menv /L bol (L ⊙ ) (K) (M ⊙ ) (%) c (M ⊙ /L ⊙ )
5.5±0.5 ∼175 0.2 a -0.7 b ∼2 0.04 a -0.13 b a Assuming T d =20 K, κ1.2=0.01 cm 2 g -1 and r 5,000 AU.

b Mass within r=Rout 15,000 AU. c Lsmm is the "submm" luminosity radiated at λ >350 µm.

observed 1.2 mm emission profile with two power-laws. Fig. 4 shows that the dust emission in the outer envelope (r 5,000 AU one beam) follows a very steep profile, with m ≃1.9, connecting the envelope with the lower density ambient cloud at r ≥ R out . Since the outer envelope angular radii are larger than the beam's FWHM, the convolution does not modify the slope of the intensity profile (e.g., Adams 1991;Motte & André 2001). For the same reason, our observations can not resolve accurately the slope of the inner envelope (r < R in ≃5,000 AU). We take m 0.5 as the expected value in a free-falling envelope.

The line radiative transfer models shown in Sect. 3.4 are not inconsistent with this assumption.

Taking again a dust opacity of κ 1200 =0.01 cm 2 g -1 and T d =20 K (consistent with the SED analysis) we compute a total circumstellar mass of M env =0.7 M ⊙ within r ≤ R out . However, it is unlikely that all this mass can be finally accreted by the protostar/disk. First because the outflow will disperse a large fraction of it, and second, because differential rotation in the envelope is expected to decouple the inner rapidly infalling and rotating layers from the slowly collapsing and rotating outer regions (Belloche et al. 2002). Hence, we also computed M in env =0.2 M ⊙ , the mass within one 11 ′′ -beam (r R in ) as representative of the mass inside the inner envelope. 4.-1.2 mm radial intensity profile around IRAS 22129+7000. The intensity (in mJy/11 ′′ -beam) has been averaged over circular annulus in steps of 3.5 ′′ . The error bars show the standard deviation within each annulus. The dashed line represent a Gaussian beam with a FWHM of 11 ′′ . The right axis shows an estimate of the H 2 column density assuming κ 1.2 =0.01 cm 2 g -1 and T d =20 K. The red curve shows a fit to the data including 2 intensity power-law profiles (convolved with the beam), with I(r) ∝ r -1.9 for r 5,000 AU and I(r) ∝ r -0.5 for r 5,000 AU.

The steep slope of the 1.2 mm intensity profile suggests that a temperature gradient is also present. In particular, the dust thermal profile through the inner envelope should be dominated by accretion heating from the protostar/disk system with (Terebey et al. 1993):

T d (r) ≈ 38 K r 100 AU -q L ⋆ L ⊙ q/2 (1)
where q depends on the dust grain properties through q = 2/(4 + β). Given the proximity of Ced 201 PDR, FUV photons from BD+69 • 1231 star also heat the outer layers of IRAS 22129+7000 envelope. Penetration of this external FUV radiation field into the envelope greatly depends on the scattering properties of the dust particles, with larger penetration depths as grains grow toward bigger grains (Goicoechea & Le Bourlot 2007). In the following we adopt a temperature profile that follows Eq. ( 1), reaches a minimum of ∼10 K at r ∼10,000 AU and then increases up to ≃15 K at R out due to external FUV heating (this temperature is consistent with our analysis of CO lines in the ambient cloud). Taking q=0.36 (using β=1.6 from the SED) and assuming L ⋆ ≃ L bol , the estimated mass of the envelope increases by a factor ∼2 respect to the isothermal case because T d (r) drops below 20 K in a large fraction of the envelope. This difference adds to the dominant source of uncertainty in the mass determination, the dust composition and associated opacity.

Taking into account these caveats, the observed 1.2 mm intensity profile in the outer envelope is equivalent to a steep density profile ρ(r > R in ) ∝ r -2.54 with a number density of n(H 2 )≃(1-2)×105 cm -3 at R in . Under these assumptions, the (beam averaged) extinction towards the peak (e.g., Motte & André 2001) is A V ≃20-30, a factor ∼10 larger than the value implied by our 1.2 mm continuum observations in the ambient cloud.

Gas Kinematics: Outflow and Infall?

Figure 2 shows the outflow emanating from IRAS 22129+7000 and extending ∼0.3 pc in the NW direction (redshifted lobe) and ∼0.2 pc in the SE direction (blueshifted lobe). In addition, IRAS 22129+7000 is the driving source of the Herbig Haro HH450 jet detected ∼0.2 pc SE from the object [START_REF] Bally | [END_REF]. The CO outflow is well collimated, with an aspect ratio of ∼3-4 (the projected FWHM flow length/width ratio). CO J=2-1 to 4-3 line profiles show high velocity wing emission from accelerated gas. Typical full line widths are consistent with an expansion velocity of v exp ≃7 km s -1 (Fig. 1). We thus estimate that the dynamical age of the outflow (∼length/v exp ) is ∼35,000 yr (up to ∼10 5 yr if the inclination angle is i=18 o ). The presence of a molecular outflow suggests that the protostar/disk system is still accreting material from the natal envelope. The high density tracer HCO + J=1-0 line, with a critical density of n 1-0 cr ≃2 × 10 5 cm -3 , shows a deep red-shifted-absorption asymmetry in its optically thick line profile (when n(H 2 ) > n J-J ′ cr the associated J level is predominantly populated by collisions with H 2 molecules and T ex , the transition excitation temperature, tends to the kinetic temperature). The depth and velocity position of the absorption dip are due to a decrease of T ex at the outer edges of a (likely) infalling envelope (left insets in Fig. 1). As low-J CO lines easily thermalize at lower densities (n cr ≃10 3-4 cm -3 ), they may not be reliable tracers of the presence of infall motions but of the varying physical conditions in the beam. Among the observed CO lines, only the J=4-3 line starts to trace denser (n 4-3 cr ≃5 × 10 4 cm -3 ), possibly infalling, gas. In anycase, the HCO + J=1-0 self-absorption dip is redshifted 0.3 km s -1 from the rest velocity of the cloud (-5 km s -1 ) and is not observed at nearby positions of the cloud where the line shows a pure Gaussian emission profile (right insets in Fig. 1). Although higher-J optically thick lines from abundant high dipole moment molecules should be mapped to draw definitive conclusions, in the following we assume that the observed deep HCO + self-absorption is produced by inward motions.

The densest gas in the envelope has been traced with the DCO + J=3-2 line (n 3-2 cr ≃2 × 10 6 cm -3 ) which shows compact ( 40 ′′ ) emission around IRAS 22129+7000. Compared to the broad line wings from the CO outflow, the DCO + J=3-2 line shows the narrowest observed linewidth (∆v FHWM =0.60±0.04 km s -1 ). On the other hand, we used the C 18 O J=1-0 line to trace lower density gas (e.g., the outer envelope). This line displays a broader linewidth (∆v FHWM =0.96±0.02 km s -1 ). The radiative transfer calculations presented in the next section confirm that both lines are optically thin and thus they are not broaden by line saturation. Therefore, the different velocity dispersions reflect the distinct regions where both line profiles are formed. Assuming T k =10 to 20 K, the observed DCO + linewidths are 5 to 3.5 times larger than the expected thermal broadening. The required non-thermal velocity dispersion to fit the observed DCO + J=3-2 linewidth (from turbulence, infall or rotation) is σ nth ≃0.25 km s -1 (with ∆v FHWM = 2.355 × σ).

Line Radiative Transfer Analysis

In order to reproduce the observed molecular lines with a realistic model of the envelope, we have used a nonlocal and non-LTE excitation and radiative transfer code (Goicoechea et al. 2006). The code models an infalling spherical envelope and accounts for line trapping, collisional excitation 5 , and radiative excitation. Arbitrary velocity fields as well as density, temperature, and abundance radial gradients are included as input parameters. This is done by numerical discretization of the envelope in spherical shells. The nonlocal radiative transfer problem is then simulated by the emission of a determined number of model photons (cosmic mm background, continuum and line photons) using the Monte Carlo approach (Bernes 1979). The steady state statistical equilibrium equations are then solved iteratively and the non-LTE level populations are determined in the shells. The emergent line intensities along each line of sight are finally convolved with the telescope angular resolution.

A First Model for the Circumstellar Envelope

As our first attempt to understand IRAS 22129+7000 physical structure, we have modelled a spherical envelope with a radius R out =15,000 AU. We take the density and temperature profiles (ρ(r) ∝ r -2.54 and T (r) ∝ r -0.34 plus external FUV heating) derived in the previous section for the outer envelope (r > R in ≃5,000 AU). Given the high densities expected in most of the envelope, we assume that gas and dust are thermalized (i.e., T k =T d ).

Using the above density and temperature structure we run the model globally for HCO + , DCO + and C 18 O and try to find iteratively the best velocity fields and abundances fitting the observed line profiles. In order to reproduce the observed velocity of the CO J=4-3 and HCO + J=1-0 redshifted absorption dips (0.2 and 0.3 km s -1 ) we adopt a uniform infalling velocity of v inf =0.2 km s -1 in the outer envelope. For r ≤ R in (the inner envelope) we assume ρ(r) ∝ r -1.5 , and v inf (r) = 0.2 (r/R in ) -0.5 i.e., the density and velocity profiles in a freefall inside-out collapsing envelope (Shu et al. 1987). Numerical simulations of collapse induced by external compression waves also predict those approximate ρ(r) and v inf (r) profiles behind the expanding collapsewavefront (e.g., Hennebelle et al. 2003). These two regimes of infall (freefall and uniform) have been inferred in other YSOs with clear evidences of collapse (Belloche et al. 2002). Under the above conditions, DCO + optically thin lines constrain the remaining nonthermal broadening (turbulence in the following). Satisfactory fits to the DCO + J=3-2 line (with [DCO + /H 2 ]=(1-2)×10 -10 ) are obtained for σ turb =0.15 km s -1 (the turbulence on the size scale of the beam). Therefore, σ turb in the densest parts of the envelope (as traced by the DCO + ) is subsonic, and probably a vestige of the low degree of turbulence found in prestellar cores (Bergin & Tafalla 2007). Hereafter we adopt a uniform σ turb value of 0.15 km s -1 in the models of the envelope.

In order to reproduce the HCO + J=1-0 line profile we have added the contribution from the lower density (more turbulent) ambient cloud surrounding IRAS 22129+7000 and producing a line-of-sight extinction of A V =2-3 (from 1.2 mm observations). Given the high HCO + dipole moment and large abundance (producing subthermal excitation and large opacities), low-J line photons from the dense envelope will be greatly attenuated and scattered over large areas by this low density halo (Cernicharo & Guélin 1987). Note that the negligible DCO + and C 18 O column densities in the ambient cloud do not alter the DCO + and C 18 O line emission from the envelope itself. The density and temperature in the ambient cloud can be constrained with the observed CO lines. They are reproduced with n(H 2 )≃1000 cm -3 and T k ≃15 K, in agreement with estimations by Witt et al. (1987). The halo ambient cloud in the model is at rest (v inf =0 km s -1 ) but the broader observed linewidths require σ turb =0.30 km s -1 (supersonic), at least a factor 2 larger than in the circumstellar envelope. Figure 5 shows the best fits and the density and temperaure gradients used in the model. In spite of the good agreement, the model slightly underestimates the observed HCO + and C 18 O J=1-0 linewidths, likely due to a minor contribution from the high velocity outflow (not modelled here).

The resulting large DCO + /HCO + abundance ratio (∼0.02) found in IRAS 22129+7000 is similar to that derived in cold and dense prestellar condensations near other PDRs (e.g., Pety et al. 2007). This is consistent for a species that shows enhanced abundance in the cold FUV-shielded gas (low ionization degree) affected by ongoing deuterium fractionation (where dissociative recombinations with electrons are its main destruction route). As a consequence, the DCO + /HCO + abundance ratio is roughly inversely proportional to the electrons abundance and we can estimate its value. Taking the same gas phase chemical network used by us in Pety et al. (2007), and assuming a cosmic ionization rate of ζ=3×10 -17 s -1 (Dalgarno 2006), we derive an ionization fraction of ∼10 -8 (with respect to H nuclei) in the densest and coldest regions of the envelope where the DCO + abundance peaks. This value is, within a factor ∼2, similar to the ionization fraction inferred in prestellar cores such as Barnard 68 [START_REF] Maret | [END_REF].

CO outflow parameters

In the previous section we have estimated the amount of material towards IRAS 22129+7000 and also evaluated the CO line emission contribution from the ambient cloud and from the circumstellar envelope (through the optically thin C 18 O line). After subtracting the CO contribution from the above line-of-sight components we now estimate a CO column density of ∼10 16 cm -2 in the outflow. Taking into account the observed dimensions of the blue-and red-shifted lobes we compute the mass (∼0.14 M ⊙ ), momentum (∼1.0 M ⊙ km s -1 ) and kinetic energy (∼6.8×10 43 erg) carried out by the outflow. Note that the mechanical luminosity of the flow (∼0.02 L ⊙ ) is still a significant fraction (∼0.4%) of the radiant luminosity of the driving source. The estimated outflow momentum flux is F CO ∼3×10 -5 M ⊙ km s -1 yr -1 , although the correction factor for inclination effects (sin i/cos 2 i) will decrease the actual F CO if i <40 o . Hence, F CO is intermediate between the observed momentum fluxes found in early Class 0 sources and those found in late Class I sources (e.g, Bontemps et al. 1996, Bachiller 1996).

DISCUSSION AND CONCLUSIONS

The Nature of IRAS 22129+7000

As noted earlier IRAS 22129+7000 shows indirect evidence of a protostar/disk system (a collimated outflow) and it shows extended dust emission from the envelope with a large L smm /L bol ratio (≃2%). These characteristics are consistent with the standard observational definition of Class 0 sources [START_REF] André | Protostars and Planets IV[END_REF]. Comparison of the derived L bol and M in env values with theoretical evolutionary tracks (Saraceno et al. 1996;[START_REF] André | Protostars and Planets IV[END_REF] shows that the source is likely in the interesting Class 0 to I transition. Approximately 50% of the initial mass envelope must have been accreted at this stage. According to these diagrams, infall must have started in the parental core ∼(0.5-1)×10 5 yr ago. This age is consistent with the dynamical age of the outflow, specially if inclination effects are taken into account.

The peculiar SED of IRAS 22129+7000 shows however a near-and mid-IR excess that produces a high T bol of 175 K, consistent with more evolved Class I sources (∼70-650 K; Myers et al. 1998). Given the strong dependence of the short-wavelength IR emission on the viewing angle (Fig. 3), the inferred T bol is certainly affected by orientation effects. Up-to-date 2D SED models show however that inclination is not enough to reproduce the observed excess. Nevertheless, as mentioned by Robitaille et al. (2006), their large grid of modelled SEDs have some caveats. The most important for our analysis is that none of them includes multiple central sources. An interesting possibility is thus that the source in IRAS 22129+7000 is a multiple (e.g., binary) stellar system. In fact, prestellar cores can fragment a second time during protostellar collapse and lead to the formation of binary systems (e.g. Duchêne et al. 2007). However, it is not yet clear if the initial conditions in these cores (e.g., the cloud environment) play a major role in such dynamical fragmentation. Simulations predict that external compression waves arriving at the central object at different times or instabilities in the outer accretion disk could be the seed of multiple protostellar systems (e.g., Hennebelle et al. 2003Hennebelle et al. , 2004)). In this picture, the presence of both a young and a more evolved YSO sharing a common circumstellar envelope would explain the observed near-and mid-IR excess if the scattered light from the evolved YSO escapes from the outflow cavities.

Higher angular resolution IR observations are required to confirm this scenario. In addition, a larger set of molecular line observations and maps are needed to tightly constrain the physics (infall and rotation velocity fields) and chemistry (molecular content) of this protostar and its environment. ). Taking into account geometrical dilution and selective dust extinction of FUV dissociating radiation from BD+69 • 1231, Casey (1991) computed that a region with a radius of at least ∼0.6 pc (∼300 ′′ ) around the exciting star is permeated by FUV photons. In particular the FUV radiation field near IRAS 22129+7000 should be presently at least ∼2 times the mean interstellar radiation field. Aside from dynamical effects associated with the presence of a nearby PDR illuminated by a runaway star (e.g., shock waves), heating and radiation pressure in the region are larger than in other FUV radiation-free environments. Interestingly enough, the 1.2 mm dust emission (very weak in the PDR itself) peaks again in two dust condensations ∼0.25 pc North of Ced 201, roughly in the same direction of outflow red-lobe axis. These cores/clumps have a radius of ∼5,000 AU and coincide with an enhancement of the HCO + J=1-0 line emission (see Fig. 1). It is very likely that these condensations are being externally heated and compressed by the combined action of the outflow and the PDR associated FUV-driven shock wave. In fact, a similar scenario has presumably led to the formation of IRAS 22129+7000 itself. We thus suspect that star formation in the B175 globule is triggered by the shocks and radiative feedback induced by BD+69 • 1231 star and the new generation of low-mass YSOs.

Observations of cores and YSOs provide new insights into how the star formation process is affected by its cloud environment. The vicinity of late type B stars has interesting intermediate properties between isolated dark cores (where FUV-radiation does not play a major role) and cores compressed by the expansion of H ii regions (near more massive stars). Cool PDRs such as Ced 201 are nearly devoid of ionizing photons and thus they provide a very different environment compared to other star forming regions which are photoionized by brighter OB stars. In the special case of Ced 201, the exciting star is moving fast through the region, suggesting that its arrival to the cloud, and thus the triggering, is recent.

Fig. 1 .

 1 Fig. 1.-Middle panel: IRAM-30m/MAMBO 1.2 mm-continuum emission of the region. The image is centered at IRAS 22129+7000.White contours are displayed from 22 (2σ) to 55 mJy/15 ′′ -beam in steps of 11 mJy/15 ′′ -beam (1σ). The CO J=3-2 integrated intensity (from 5 to 9 K km s -1 in steps of 1 K km s -1 ) from the outflow is shown in red (from v=-4 to -2 km s -1 ) and blue (from v=-6 to -8 km s -1 ). The integrated CO J=3-2 line intensity from Ced 201 PDR (v=-5.5 to -4.5 km s -1 ) is shown in black contours (from 6 to 20 K km s -1 in steps of 1 K km s -1 ). Left panels: 12 CO J=2-1 and 4-3, HCO + J=1-0 and DCO + J=3-2 lines towards IRAS 22129+7000 position. Right panels: Declination raster in the HCO + J=1-0 line. The vertical dashed line marks the rest velocity of the cloud (-5 km s -1 ).

Fig. 2 .

 2 Fig. 2.-Spitzer /IRAC 8 µm image of Ced201 region including the detection of IRAS 22129+7000 protostar. The bright reflexion nebula (the PDR) is clearly seen at the west of the region. A well collimated outflow (as traced by CO J=3-2 high velocity line wing emission)emerges from the protostar position. Contours show the integrated intensity (from 5 to 9 K km s -1 in steps of 1 K km s -1 ) in the red (from v=-4 to -2 km s -1 ) and blue (from v=-6 to -8 km s -1 ) line wings respectively.

Fig. 3 .

 3 Fig. 3.-SED of IRAS 22129+7000 combining data in the H and Ks bands (2MASS; Skrutskie et al. 2006); 4.5, 5.8 and 8 µm (Spitzer/IRAC); 20-35 µm (Spitzer/IRS, smaller points), 12, 25, 60 and 100 µm (IRAS), 70 µm (Spitzer/MIPS), 200 µm (KAO: Casey1991) and 1200 µm (MAMBO). The Spitzer/IRS fluxes towards the source agree with those measured by IRAS at 25µm probing the compact nature of the mid-IR emission. The dashed curve is a simple fit using two graybodies with different temperatures. The continuous blue curve represents the circumstellar envelope with T d =22 K and a spectral index of 1.6. The magenta pointed-and pointed-dashed curves show the best 2D radiative transfer models obtained with theRobitaille et al. (2007) fitter for two different inclinations (i = 0 corresponds to "pole-on" and i = 90 to "edgeon" configurations).

Fig.

  Fig.4.-1.2 mm radial intensity profile around IRAS 22129+7000. The intensity (in mJy/11 ′′ -beam) has been averaged over circular annulus in steps of 3.5 ′′ . The error bars show the standard deviation within each annulus. The dashed line represent a Gaussian beam with a FWHM of 11 ′′ . The right axis shows an estimate of the H 2 column density assuming κ 1.2 =0.01 cm 2 g -1 and T d =20 K. The red curve shows a fit to the data including 2 intensity power-law profiles (convolved with the beam), with I(r) ∝ r -1.9 for r 5,000 AU and I(r) ∝ r -0.5 for r 5,000 AU.

Fig. 5 .

 5 Fig. 5.-(a) HCO + J=1-0, (b) DCO + J=3-2 and (c) C 18 O J=1-0 observations (histograms) towards IRAS 22129+7000. Best radiative transfer models for the circumstellar envelope are also shown (red curves). Predicted line profiles have been convolved with the telescope angular resolution at each frequency. Intensity scale is in main beam temperature. (d) Density and temperature profiles used in the model.

4. 2 .

 2 Star Formation Near Ced 201 PDR IRAS 22129+7000 is located near Ced 201 PDR (d ∼0.4 pc and A V ∼1

  ⊙ ) (M ⊙ km s -1 ) (10 43 erg) (10 -3 L ⊙ ) To calculate the energy and momentum we use an expansion velocity of 7 km s -1 and a dynamical age of 35,000 yr (see text).

			TABLE 3
	IRAS 22129+700: CO outflow parameters
				Kinetic	Mechanical
	Outflow	Mass	Momentum	Energy	Luminosity
	Component a (M Blueshifted... 0.06	0.4	2.9	6.9
	Redshifted....	0.08	0.6	3.9	9.3
	Total............	0.14	1.0	6.8	16.2

a

http://www.cesr.fr/∼joblin/SPECPDR public/Home

CO (and isotopologues) rotationally inelastic collisional rates have been derived from those ofFlower (2001). HCO + and DCO + collisional rates have been derived from those ofFlower (1999).

We thank P. Hennebelle for useful advice, J. Pety for his help with CLASS90 and the SPECPDR team for their contribution to the project. We also thank the referee for his/her constructive criticisms. JRG was supported by a Marie Curie intra-European Individual Fellowship, contract MEIF-CT-2005-515340. We made use of data products from the 2MASS Catalog, which is a joint project of the U. of Massachusetts and the IPAC/Caltech, funded by NASA and the NSF."

1 Based on observations obtained with the CSO telescope. The CSO telescope is operated by the California Institute of Technology under funding from the NSF, Grant No. AST-0540882.

2 Observations also obtained with IRAM telescopes, supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain).

3 This work is based on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA.


Aure ´lie Chaalali

Gre ´gory Beaugrand

Virginie Raybaud

Eric Goberville

Vale ´rie David

Philippe Boe

Benoit Sautour

Climatic Facilitation of the Colonization of an Estuary by Acartia tonsa

Global change has become a major driving force of both terrestrial and marine systems. Located at the interface between these two realms, estuarine ecosystems are probably the place where both direct and indirect effects of human activities conspire together to affect biodiversity from phytoplankton to top predators. Among European estuarine systems, the Gironde is the largest estuary of Western Europe and many studies have provided evidence that it has been affected by a variety of anthropogenic stressors such as thermal and chemical pollution, physical alterations and exploitation, especially for maritime traffic. In such a context, species introduction is also a current major issue with the establishment of strong competitive species that could lead to ecosystem reorganization with potential decrease or even disappearance of native species. In the Gironde estuary, this hypothesis was proposed for the invasive shrimp species Palaemon macrodactylus as a decrease in the native species abundance was observed at the same time. Although species introduction often takes place via ballast water, the influence of climate-driven changes on the establishment of new species remains a key issue. The calanoid copepod Acartia tonsa, observed in the Gironde estuary for the first time in 1983, have since colonized most part of the estuary, reaching a level of abundance comparable to the dominant native species Eurytemora affinis. In this study, using both the concept of the ecological niche sensu Hutchinson (fundamental and realized niches) and statistical models, we reveal that the dynamics of the colonization of A. tonsa was facilitated by environmental conditions that have become closer to its environmental optimum with respect to temperature and salinity.

Introduction

Marine biodiversity and ecosystems are being altered by many human-induced factors including overexploitation of marine resources [START_REF] Jennings | The effects of fishing on marine ecosystems[END_REF][START_REF] Omori | Impact Of Human Activities On Pelagic Biogeography[END_REF], chemical pollution and physical alterations [START_REF] Norse | Global marine biological diversity. A strategy for building conservation into decision making[END_REF], eutrophication and invasion of exotic species [START_REF] Carlton | Ecological Roulette -The Global Transport Of Nonindigenous Marine Organisms[END_REF][START_REF] Edwards | Case history and persistence of the non-indigenous diatom Coscinodiscus wailesii in the north-east Atlantic[END_REF][START_REF] Staneck | Human influences on coastal ecosystems: does overfishing create trophic cascades?[END_REF]. Fisheries have affected the marine environment through both direct and indirect effects [START_REF] Jennings | The effects of fishing on marine ecosystems[END_REF][START_REF] Omori | Impact Of Human Activities On Pelagic Biogeography[END_REF][START_REF] Omori | Anthropogenic short-term effects on distribution patterns and food webs of pelagic organisms[END_REF]. After habitat degradation, the introduction of nonnative species in an environment is a major cause of extinction [START_REF] Macdonald | Global change and alien invasions: implications for biodiversity and protected area management[END_REF][START_REF] Vitousek | Human domination of Earth's ecosystems[END_REF]. So far, marine invasions have been less investigated although their magnitude and frequency may lead to profound changes in ecosystem functioning and biological community structure [START_REF] Carlton | Ecological Roulette -The Global Transport Of Nonindigenous Marine Organisms[END_REF][START_REF] Jousson | Invasive alga reaches California -The alga has been identified that threatens to smother Californian coastal ecosystems[END_REF][START_REF] Ricciardi | Recent mass invasion of the North American Great Lakes by Ponto-Caspian species[END_REF].

Estuaries, located at the interface between terrestrial and marine realms, constitute habitat for many species or even nursery, refuge and growth areas [START_REF] Beck | The identification, conservation, and management of estuarine and marine nurseries for fish and invertebrates[END_REF][START_REF] Blaber | Factors Affecting The Distribution Of Juvenile Estuarine And Inshore Fish[END_REF][START_REF] Haedrich | Estuarine fishes[END_REF]. However, these ecosystems are being affected by human activities such as fishing, polluting, and maritime traffic [START_REF] Beaugrand | Spatial changes in the sensitivity of Atlantic cod to climate-driven effects in the plankton[END_REF][START_REF] Pauly | Ecosystem impacts of the world's marine fisheries[END_REF]. For example, maritime traffic increases the number of invasive species (e.g. copepods [START_REF] David | Successful colonization of the calanoid copepod Acartia tonsa in the oligo-mesohaline area of the Gironde estuary (SW France) -Natural or anthropogenic forcing?[END_REF]; shrimps [START_REF] Beguer | Presence of Palaemon macrodactylus in a European estuary: evidence for a successful invasion of the Gironde (SW France)[END_REF]) potentially transported via ballast waters. Survival of these species depends on environmental conditions within their new ecosystem and climate change may influence the establishment of alien species [START_REF] Gucu | Can overfishing be responsible for the successful establishment of Mnemiopsis leidyi in the Black Sea?[END_REF]. For example, Raitsos et al. [START_REF] Raitsos | Global climate change amplifies the entry of tropical species into the Eastern Mediterranean Sea[END_REF] showed that the appearance of benthic and pelagic tropical species in the eastern part of the Mediterranean Sea was highly positively correlated with sea surface warming.

When they manage to survive, invasive species may impact their new ecosystem by affecting both its structure and functioning with implications for the structure and function of the entire food web [START_REF] David | Successful colonization of the calanoid copepod Acartia tonsa in the oligo-mesohaline area of the Gironde estuary (SW France) -Natural or anthropogenic forcing?[END_REF][START_REF] Beguer | Presence of Palaemon macrodactylus in a European estuary: evidence for a successful invasion of the Gironde (SW France)[END_REF][START_REF] Aravena | Response of Acartia populations to environmental variability and effects of invasive congenerics in the estuary of Bilbao, Bay of Biscay[END_REF]. For example, the establishment of the invasive copepod species Acartia tonsa in the oligo-mesohaline zone of the Gironde estuary led to a phenological shift in the native copepod Acartia bifilosa production period being advanced one month earlier in the year [START_REF] David | Successful colonization of the calanoid copepod Acartia tonsa in the oligo-mesohaline area of the Gironde estuary (SW France) -Natural or anthropogenic forcing?[END_REF]. Acartia tonsa also invaded the estuary of Bilbao and altered the spatial distribution of the native copepod Acartia clausi [START_REF] Aravena | Response of Acartia populations to environmental variability and effects of invasive congenerics in the estuary of Bilbao, Bay of Biscay[END_REF].

In this study, we used the concept of the ecological niche sensu Hutchinson (fundamental and realized niches) to investigate the drivers involved in the colonization of A. tonsa. We considered two climate-driven parameters, temperature and salinity, for their wellknown influence on biological processes. Indeed, temperature controls the kinetics of many metabolic reactions (e.g. speed of enzymatic reaction, increase in metabolism; Arrhenius [START_REF] Arrhenius | On the reaction velocity of the inversion of cane sugar by acids[END_REF]) and influences species reproduction, locomotion, feeding rates and interaction between species [START_REF] Beaugrand | Long-term changes in phytoplankton, zooplankton and salmon related to climate[END_REF][START_REF] Magnuson | Temperature as an ecological resource[END_REF][START_REF] Mauchline | The Biology of Calanoid Copepods[END_REF]. Salinity may induce potential osmotic stress, increasing the risk of species mortality [START_REF] Vilas | Inter-and intra-specific differences in euryhalinity determine the spatial distribution of mysids in a temperate European estuary[END_REF]. Considering these two parameters, we evaluated the importance of climate change in the colonization of the Gironde estuary by the copepod Acartia tonsa. We applied a double modeling approach, based on both realized and fundamental niches, to provide evidence that the successful establishment of Acartia tonsa was the result of changes in both temperature and salinity.

Materials and Methods

Study area

The Gironde estuary (latitude 45u20'N, longitude 0u459W), a 70-km long estuary formed by the junction of the Dordogne and Garonne rivers [START_REF] Allen | Study on circulation of water masses at mouth of Gironde[END_REF], is the largest south-western European estuary. With a mean suspended mater concentration higher than 500 mg.L 21 [START_REF] Abril | Oxic/anoxic oscillations and organic carbon mineralization in an estuarine maximum turbidity zone (The Gironde, France)[END_REF][START_REF] Sautour | Comparative spring distribution of zooplankton in 3 macrotidal European estuaries[END_REF], it is one of the most turbid European estuaries. The Gironde is characterized by a strong spatial and temporal variability in both physical and chemical properties [START_REF] David | Long-term changes of the zooplankton variability in a turbid environment: The Gironde estuary (France)[END_REF]. In addition to the natural variability, this ecosystem is also subjected to anthropogenic pressures including climate changes [START_REF] David | Successful colonization of the calanoid copepod Acartia tonsa in the oligo-mesohaline area of the Gironde estuary (SW France) -Natural or anthropogenic forcing?[END_REF][START_REF] David | Long-term changes of the zooplankton variability in a turbid environment: The Gironde estuary (France)[END_REF][START_REF] Chaalali | Climate-caused abrupt shifts in a European macrotidal estuary[END_REF][START_REF] Coupry | E ´valuation des impacts du changement climatique sur l'estuaire de la Gironde et prospective a `moyen terme[END_REF][START_REF] Goberville | Climate-driven changes in coastal marine systems of western Europe[END_REF]. Among the other anthropogenic pressures, the maritime traffic (Bordeaux harbor; hull fouling, semi-dry ballast, and ballast water) represents an important cause of species introduction [START_REF] David | Successful colonization of the calanoid copepod Acartia tonsa in the oligo-mesohaline area of the Gironde estuary (SW France) -Natural or anthropogenic forcing?[END_REF][START_REF] Beguer | Presence of Palaemon macrodactylus in a European estuary: evidence for a successful invasion of the Gironde (SW France)[END_REF]. Therefore, the Gironde estuary can be considered as a hotspot for species introduction. The calanoid copepod species Acartia tonsa, first introduced in 1983 by ballast waters, is an example of species that found suitable conditions to establish in the Gironde estuary.

Environmental and zooplankton data

Environmental data were provided by both the ''Blayais'' nuclear power plant and the SOMLIT (Service d'Observation en Milieu LITtoral, a French coastal monitoring network: http:// somlit.epoc.u-bordeaux1.fr; see [START_REF] Goberville | Climate-driven changes in coastal marine systems of western Europe[END_REF] or [START_REF] Savoye | Origin and composition of particulate organic matter in a macrotidal turbid estuary: The Gironde Estuary, France[END_REF]) joint monitoring. Twelve environmental parameters and the five zooplankton species abundances (including 3 copepods and 2 mysids) were monitored at monthly intervals from March 1978 to November 2009. Here, we only focused on water temperature, salinity and abundance of Acartia tonsa per m 3 . Samples were collected nine times a year (except for site F monitored eight times a year) at three given sites (Figure 1): (1) site E, located at 52 km from Bordeaux city, since 1978; (2) site K, located at 30 km from Bordeaux city, since 1984; (3) site F, located at 67 km from Bordeaux city, since 1992 [START_REF] David | Long-term changes of the zooplankton variability in a turbid environment: The Gironde estuary (France)[END_REF]. Sampling was carried out at 1 m below the surface and 1 m above the bottom at 3 h intervals during a tidal cycle (high and low tide, flood and ebb tide).

We combined all triplets of temperature/salinity/Acartia tonsa abundance available in all sampling sites and obtained the data matrix: 5178 observations 63 variables (Matrix 1). To only consider the period since species establishment, from 1999 to 2010, (the year 2011 being incomplete), we extracted data from Matrix 1. Only considering data from site E, we obtained a second matrix: 920 observations 63 variables (Matrix 2).

Data analyses

Based on these two matrices, we determined the ecological niche of Acartia tonsa as a function of both monthly water temperature and salinity. Analysis 1 determined the realized niche, i.e. the niche including the effect of species interaction and dispersal [START_REF] Hutchinson | An introduction to population ecology[END_REF][START_REF] Pulliam | On the relationship between niche and distribution[END_REF] by using categories of temperature and salinity in a discrete model. Analysis 2, applying a mixed Gaussian-linear model, also estimated species realized niche. Analysis 3 determined the fundamental niche, i.e. the niche without the effect of dispersal or species interaction [START_REF] Hutchinson | An introduction to population ecology[END_REF]. Physiological thresholds obtained from literature were subsequently implemented in the new ecological niche model NPPEN (Non-Parametric Probabilistic Ecological Niche model; [START_REF] Beaugrand | A new model to assess the probability of occurrence of a species, based on presence-only data[END_REF]). This model presents several advantages: it requires presence-only data and is based on a non-parametric procedure using a simplification of the Multiple Response Permutation Procedure (MRPP) based on the Generalized Mahalanobis distance [START_REF] Beaugrand | A new model to assess the probability of occurrence of a species, based on presence-only data[END_REF]. NPPEN was already applied to many marine organisms, for instance as fishes [START_REF] Beaugrand | A new model to assess the probability of occurrence of a species, based on presence-only data[END_REF][START_REF] Lenoir | Modelled spatial distribution of marine fish and projected modifications in the North Atlantic Ocean[END_REF] or benthic macrofauna [START_REF] Rombouts | Potential changes in benthic macrofaunal distributions from the English Channel simulated under climate change scenarios[END_REF].

The three models (discrete, Gaussian-linear and NPPEN) were then compared using Taylor diagrams (Analysis 4; [START_REF] Taylor | Summarizing multiple aspects of model performance in a single diagram[END_REF]).

Analysis 1: Characterization of the realized niche by discrete model (DM). Using an approach similar to Beaugrand et al. [START_REF] Beaugrand | Macroecological study of Centropages typicus in the North Atlantic Ocean[END_REF], we used Matrix 2 (1999-2010 for site E) to compute the realized niche of A. tonsa (abundances transformed in log 10 (x+1)) as a function of monthly water temperature (from 0 to 35uC every 1uC) and salinity (from 0 to 35 psu every 1 psu). This estimation of A. tonsa realized niche for the period 1999-2010 was compared to the realized niche computed from all available data (for all sampling sites, i.e. Matrix 1). Using only data from site E (the longest time series) and the reference period 1999-2010, we removed any problem of circularity as the test of the model (Analysis 4) was performed on data from site E and other sites (F and K) for the period 1984-1998. From the knowledge of monthly water temperature and salinity of all sites, we retrieved the expected abundance of A. tonsa in all sites from 1984 to 2011 (see Figure 2B).

Analysis 2: Characterization of the realized niche by a Mixed Gaussian-linear (MGL) model. For this second analysis, we used monthly water temperature and salinity data of sampling site E for the period 1999-2010 (Matrix 2) to estimate the realized niche of A. tonsa. We graphically examined the relationships between the abundance of A. tonsa and the two environmental variables (not shown). The model was Gaussian for temperature because the examination of the data showed a maximum of abundance (as log 10 (x+1)) flanked by two decreasing slopes. Such a function is often used to model the ecological niche [START_REF] Fort | Biologging, Remotely-Sensed Oceanography and the Continuous Plankton Recorder Reveal the Environmental Determinants of a Seabird Wintering Hotspot[END_REF][START_REF] Helaouet | Physiology, ecological niches and species distribution[END_REF][START_REF] Nicolas | Impact of global warming on European tidal estuaries: some evidence of northward migration of estuarine fish species[END_REF]. As the relationship between salinity and abundance (as log 10 (x+1)) was linear between 0 and 17 psu and leveled off after, we used a linear function between 0 and 17 psu. For salinities under 17 psu, the resulting model was the product of a linear relationship for salinity 'S' and a Gaussian relationship for water temperature 't' (Eq. 1). For salinities over 17 psu, the model only considered the effect of temperature (Eq. 2):

For S,17 psu:

A~S 17 M Aobs e { t{u ð Þ 2 2s 2 ð1Þ
For S$17 psu:

A~M Aobs e { t{u ð Þ 2 2s 2 ð2Þ
with 91/179 the constant value of the linear relationship between A. tonsa and salinity S, u and s the thermal optimum and thermal amplitude of A. tonsa, respectively.

M Aobs is the maximum observed abundance for A. tonsa. The coefficients M Aobs , a, u, and s were estimated by minimizing the sum of squares of the residuals assessed by calculating the differences between observed and predicted abundance (as log 10 (x+1)) of A. tonsa [START_REF] Fort | Biologging, Remotely-Sensed Oceanography and the Continuous Plankton Recorder Reveal the Environmental Determinants of a Seabird Wintering Hotspot[END_REF]:

F (M Aobs ,u,s)~X n i~1 ½L i { f (M Aobs ,u,s; S i ,t i ) 2 ð3Þ
with n the number of observed data (920 observations). From the knowledge of the monthly water temperature and salinity of all sampling sites, we calculated the expected abundance of A. tonsa for each site from 1984 to 2011.

Analysis 3: Characterization of the fundamental niche by the model NPPEN. To characterize the fundamental niche of A. tonsa, we incorporated physiological thresholds to the ecological niche model NPPEN. The model was parameterized using species physiological survival limits. According to the literature [START_REF] Cervetto | Influence of salinity on the distribution of Acartia tonsa (Copepoda, Calanoida)[END_REF][START_REF] Gonzales | Critical thermal maxima and upper temperature for the calanoid copepods Acartia tonsa and A. clausi[END_REF],

we determined a range of salinity between 0 and 72 psu and a range of temperature between 22 and 42uC. The eurybiontic property of the species with respect to these two parameters is exceptional. From the knowledge of the monthly water temperature and salinity of all sites, we subsequently calculated the probability of occurrence of A. tonsa in all sites from 1984 to 2011.

For each site, the probabilities of occurrence of A. tonsa were then compared to the standardized abundances (between 0 and 1) of the species. Standardization was performed for each site by using the formula:

x à ~x-min ð Þ= max-min ð Þ ð 4Þ
with x the abundance of A. tonsa, min and max the minimum and maximum values of vector x, respectively.

Analysis 4: Models skill assessment using Taylor's diagrams. Validation and quantification of errors are important steps in any modeling study. Firstly, they provide the degree of confidence that may be given to the model for the prediction of phenomena. In addition, they can help to improve current models by revealing the mechanisms most poorly represented. There are many statistical methods to compare models and measurements and to define the level of performance of a model [START_REF] Allen | Quantifying uncertainty in highresolution coupled hydrodynamic-ecosystem models[END_REF]. Each of them gives different but complementary information. Each can be used separately to test one aspect of the model, but the simultaneous use of several methods can be performed to get an overall assessment of the adjustment given by the simulation to the observed variations. This is the solution chosen by Taylor [START_REF] Taylor | Summarizing multiple aspects of model performance in a single diagram[END_REF], which proposed to assess the performance of models by representing on a single diagram the standard deviation of the data and those of the model, the Root Mean Square Deviation (RMSD) and the correlation coefficient r between the model and the data. A detailed description of the Taylor diagram and formulas involved are given in ref 40 [START_REF] Taylor | Summarizing multiple aspects of model performance in a single diagram[END_REF]. This useful tool have first been developed for atmospheric general circulation models [START_REF] Taylor | Summarizing multiple aspects of model performance in a single diagram[END_REF] but has also been applied to physical-biogeochemical coupled models [START_REF] Raick | Model complexity and performance: How far can we simplify?[END_REF][START_REF] Raybaud | Short time-scale analysis of the NW Mediterranean ecosystem during summer-autumn transition: A 1D modelling approach[END_REF] and generalized additive models [START_REF] Loots | What controls the spatial distribution of the North Sea plaice spawning population? Confronting ecological hypotheses through a model selection framework[END_REF]. In our case, we used Taylor's diagrams to compare the results of our three models, in each sampling site: DM' for the Discrete Model, 'MGL' for the Mixed Gaussian-Linear model and 'NPPEN' for the ecological niche model. For the comparison, all data having been previously normalized between 0 and 1 using Eq. 4.

All methods were computed using the MATLAB language.

Results

Estimation of both realized and fundamental niches of A. tonsa

The discrete and mixed Gaussian-linear models, estimating the realized niche of A. tonsa (Figures 2B and2C), revealed patterns similar to the observed one (Figure 2A) despite the fact that we considered less data on Figure 2B (only data from period 1999-2010). Regarding species maximal abundances, the mixed Gaussian-linear model suggested a thermal optimum around 21.5uC (64.5uC according to the function minimizing the sum of squares of the residuals) and salinity suitable conditions over ,15 psu. Low species abundances were estimated for temperatures below 10uC and salinities below 3-5 psu (Figures 2A-C).

The fundamental niche of A. tonsa estimated by the NPPEN model was more extended than its realized niche (Figure 2D). Species environmental optima (i.e. probabilities .0.8) were obtained for a temperature ranging from 13 to 24uC and a salinity exceeding 15 psu.

Estimation of A. tonsa abundance and occurrence probabilities times series

When both observed and estimated abundances of A. tonsa were compared for each sampling site, we found that both the discrete and mixed Gaussian-linear models reproduced well the seasonality of A. tonsa. However, the two models systematically underestimated species abundances at site K where the species is less abundant (Figures 3C and4C). Regardless of the model used, both modeled and observed time series suggested a spatial saline gradient with maximal abundances observed at site F (the more downstream). Finally, both models reproduced the long-term trends observed at each site for the species (r varying between 0.35 and 0.81; Figure 5). Estimates from the two models were similar as suggested by their close positions on the Taylor diagrams (Figure 5). However, the mixed Gaussian-linear model (MGL) was more accurate than the two other models. Indeed, although the standard deviation of the discrete model (DM) was close to the observations (only for sampling site F, s DM = 0.32 and s OBS = 0.35 respectively; Figure 5A), the Pearson correlation coefficients were systematically higher for the mixed Gaussian-linear model than for the discrete model while the RMSD were low.

Year-to-year changes in the probabilities of occurrence of A. tonsa, assessed by the NPPEN model, showed a phase of increasing probabilities from 1983 (first observation of the species) to ,1992 in both sites E and K. Although the seasonality was detected, the pattern was less clear than the pattern retrieved from both the discrete and mixed Gaussian-linear models. However the NPPEN model estimated a spatial distribution of A. tonsa in accordance with the salinity gradient: higher probabilities were estimated for the site F (Figure 6A), intermediate for site E (Figure 6B) and low probabilities for site K (Figure 6C). The RMSD and correlation coefficients between observed (i.e. abundances having been standardized) and estimated times series from the NPPEN model were similar to those of the empirical models ('DM' and 'MGL'; Figure 5). Its position away from the observations on Taylor diagrams was due to the low standard deviations (s NPPEN #0.13) of the NPPEN model. One explanation could be that the occurrence probability range was already between 0 and 1 before normalization.

Relation between hydro-climatic evolution and A. tonsa abundance

The monthly temperature and salinity associated with each sampling of A. tonsa at site E for the period 1978-1983 were characterized by maximal values of ,26uC and , 14 psu (and 23.5uC and 13.5 psu for August; Figure 7 C). We recall here that we only disposed of data before 1983 for sampling site E. Changes in both salinity and temperature in the Gironde estuary (blue dots, Figure 7) between the two periods, before (1978-1983, Figure 7A) and after the apparition of A. tonsa (1984-2011, Figure 7B) highlighted the influence of warming and marinisation processes; the changes in environmental conditions being associated with higher estimated abundance of A. tonsa (Figure 7). Comparable results with the NPPEN model (results not shown) were observed. Although the estimated probabilities of occurrence of A. tonsa were low for the period 1978-1983 (between 0 and 0.4), high probabilities (between 0.7 and 1) were observed from 1984 onwards.

Discussion

Brylinski [START_REF] Brylinski | Report on the presence of Acartia tonsa Dana (Copepoda) in the harbour of Dunkirk (France) and its geographical distribution in Europe[END_REF] first presented a map of the spatial distribution of the copepod A. tonsa in Europe. Originating from the North American coast, A. tonsa could have been transported by maritime traffic to Europe. According to David et al. [START_REF] David | Successful colonization of the calanoid copepod Acartia tonsa in the oligo-mesohaline area of the Gironde estuary (SW France) -Natural or anthropogenic forcing?[END_REF], the copepod was observed in the Gironde estuary for the first time in 1983. Many studies have been conducted on the ecology of the species to assess its physiological traits [START_REF] Cervetto | Influence of salinity on the distribution of Acartia tonsa (Copepoda, Calanoida)[END_REF][START_REF] Gaudy | Comparison of the metabolism of Acartia clausi and A. tonsa: influence of temperature and salinity[END_REF][START_REF] Hojgaard | Salinity-induced quiescence in eggs of the calanoid copepod Acartia tonsa (Dana): a simple method for egg storage[END_REF][START_REF] Holste | The effects of temperature and salinity on egg production and hatching success of Baltic Acartia tonsa (Copepoda: Calanoida): a laboratory investigation[END_REF][START_REF] Pastorinho | Distribution, production, histology and histochemistry in Acartia tonsa (Copepoda: Calanoida) as means for life history determination in a temperate estuary (Mondego estuary, Portugal)[END_REF][START_REF] Peck | Effects of salinity, photoperiod and adult stocking density on egg production and egg hatching success in Acartia tonsa (Calanoida: Copepoda): Optimizing intensive[END_REF][START_REF] Svetlichny | Vulnerability of copepod eggs in salinity and temperature gradients of the Marmara and Black Seas[END_REF][START_REF] White | Egg-Production By The Calanoid Copepod Acartia tonsa In The Mesohaline Chesapeake Bay -The Importance Of Food Resources And Temperature[END_REF]. All the authors defined A. tonsa as a euryhaline and eurytherm species with a very large range of tolerance measured by both in situ and laboratory experiments [START_REF] Cervetto | Influence of salinity on the distribution of Acartia tonsa (Copepoda, Calanoida)[END_REF]. These two properties explain its ubiquity in many Atlantic Ocean estuaries [START_REF] Cervetto | Influence of salinity on the distribution of Acartia tonsa (Copepoda, Calanoida)[END_REF][START_REF] Gonzales | Critical thermal maxima and upper temperature for the calanoid copepods Acartia tonsa and A. clausi[END_REF][START_REF] Brylinski | Report on the presence of Acartia tonsa Dana (Copepoda) in the harbour of Dunkirk (France) and its geographical distribution in Europe[END_REF]. Gonzales [START_REF] Gonzales | Critical thermal maxima and upper temperature for the calanoid copepods Acartia tonsa and A. clausi[END_REF] evaluated species critical thermal maximum and lethal temperature. By gradually increasing temperature until ''uncoordinated'' movements of the species, a critical thermal maximal value around 37uC (50% of individuals being able to tolerate a temperature of 37uC for 4 h) was evaluated and a lethal temperature near 41-46.3uC was assessed. The species was able to tolerate a wide range of temperatures between 21uC and 32uC, Gonzales [START_REF] Gonzales | Critical thermal maxima and upper temperature for the calanoid copepods Acartia tonsa and A. clausi[END_REF] collected several individuals in the Narragansett Bay (Rhode Island, USA) at temperature of about 21uC. Cervetto et al. [START_REF] Cervetto | Influence of salinity on the distribution of Acartia tonsa (Copepoda, Calanoida)[END_REF] experimentally characterized both the maximum and minimum values of salinity in which A. tonsa could survive. They estimated an impressive range between ,1 and 72 psu with an optimal adaptation of the species between 15 and 22 psu for copepodites. These results were in accordance with previous observations [START_REF] Tester | Why is Acartia tonsa restricted to estuarine habitats?[END_REF], which found a maximal salinity threshold of 25 psu for nauplii stages. This upper value would limit the colonization of marine waters by the early stages of development of A. tonsa, these stages being more sensitive than the adults to changes in salinity [START_REF] Tester | Why is Acartia tonsa restricted to estuarine habitats?[END_REF].

Our estimations of both fundamental and realized niches are therefore in accordance with these experimental measurements. Estimated optimal salinity conditions were observed for values over 15-17 psu and a maximal value of 34 psu was observed in the estuary for the period 1999-2011. Moreover, our results estimated low abundances of A. tonsa for a salinity value under 3-5 psu. These results are supported by previous studies on species eggs quiescence and hatching success. Hojgaard et al. [START_REF] Hojgaard | Salinity-induced quiescence in eggs of the calanoid copepod Acartia tonsa (Dana): a simple method for egg storage[END_REF] showed that a strong quiescence process was induced at 0 salinity (and at both 17 and 25uC) and Holste and Peck ( [START_REF] Holste | The effects of temperature and salinity on egg production and hatching success of Baltic Acartia tonsa (Copepoda: Calanoida): a laboratory investigation[END_REF]; see also Peck and Holste [START_REF] Peck | Effects of salinity, photoperiod and adult stocking density on egg production and egg hatching success in Acartia tonsa (Calanoida: Copepoda): Optimizing intensive[END_REF]) observed an asymptotical increase in eggs hatching success for salinity between 0 and 17 psu, a maximum being observed between 17 and 25 psu.

Regarding temperature, we estimated an optimal value around 17-25uC with all modeling approaches and an abundance reduction being estimated under a temperature value of 10uC for empirical models and 13uC for the NPPEN model. However this species is also known to potentially produce diapause eggs (or low quality eggs; Zillioux and Gonzalez [START_REF] Zillioux | Egg dormancy in a neritic calanoid copepod and its implications to overwintering in boreal waters[END_REF]), at a thermal threshold experimentally estimated at about 10-13uC. Under this minimal temperature value, species may produce diapause eggs and no hatching was observed for eggs incubated at low temperatures (,12uC) by Holste and Peck [START_REF] Holste | The effects of temperature and salinity on egg production and hatching success of Baltic Acartia tonsa (Copepoda: Calanoida): a laboratory investigation[END_REF] and interestingly a maximum was observed between 22-23uC.

When the realized niche of A. tonsa was compared to its fundamental niche (Figures 2A,B, C versus 2D), we noticed some differences, especially for species occurrence at low temperatures (,13uC). These temperature values, mainly corresponding to winter conditions, are associated to species absence from samplings (see [START_REF] David | Long-term changes of the zooplankton variability in a turbid environment: The Gironde estuary (France)[END_REF][START_REF] Selleslagh | Seasonal succession of estuarine fish, shrimps, macrozoobenthos and plankton: Physico-chemical and trophic influence. The Gironde estuary as a case study[END_REF]) even if the NPPEN model did not estimate null probabilities of species occurrence (with all salinity values associated; from 0 to ,30 psu). However, the success of an invader population over time can not only be linked to abiotic parameters such as temperature or salinity. Biotic resistance exerted by native species can also possibly play a role and induce some differences between invader's realized and fundamental niches [START_REF] Dzialowski | Experimental effect of consumer identity on the invasion success of a non-native cladoceran[END_REF][START_REF] Parkes | Are zooplankton invasions in constructed waters facilitated by simple communities?[END_REF][START_REF] Taylor | Can biotic resistance be utilized to reduce establishment rates of non-indigenous species in constructed waters?[END_REF]. Indeed, the native copepod species Eurytemora affinis has a production period centered around March [START_REF] David | Long-term changes of the zooplankton variability in a turbid environment: The Gironde estuary (France)[END_REF]. E. affinis would be favored at this period when the Maximum Turbidity Zone ('MTZ', located upstream during summer and more downstream during both autumn and winter; Sottolichio & Castaing [START_REF] Sottolichio | A synthesis on seasonal dynamics of highlyconcentrated structures in the Gironde estuary[END_REF]) is particularly important because this species is known (i) to be able to feed on detritus and detrital vegetal matter in contrast with A. tonsa [START_REF] David | Successful colonization of the calanoid copepod Acartia tonsa in the oligo-mesohaline area of the Gironde estuary (SW France) -Natural or anthropogenic forcing?[END_REF], and (ii) to have a higher ability to select its prey among inorganic matter compared to Acartia spp. [START_REF] Gasparini | Autotrophic and heterotrophic nanoplankton in the diet of the estuarine copepods Eurytemora affinis and Acartia bifilosa[END_REF]. Thus, the NPPEN model would have not taken into account the presence of E. affinis. This hypothesis may be especially true for salinities between 0 and 18 psu (conditions encountered in the oligo-mesohaline estuarine regions). Regarding higher salinities, other species are present -some neritic copepod species -, that also have not been considered in the NPPEN model. Moreover, some of these species have recently been proved to significantly progress in the estuary in response to the estuary marinisation [START_REF] Chaalali | Changes in the distribution of copepods in the Gironde estuary: A warming and marinisation consequence? Estuarine Coastal And Shelf Science[END_REF]. In conclusion, the physiological niche of Acartia tonsa estimated by the NPPEN model only reflected what would be the fundamental niche if no other parameter was acting, for example, with no predation, resource or space competition. This may also explain the difference observed between species realized and fundamental niches.

Considering species temporal variability with only 2 climatedriven parameters, we estimated values of A. tonsa abundances at each site. Although our models underestimated species abundance especially at the more upstream site, the seasonality of A. tonsa was well reproduced with a production peak around August. This result was in accordance with species seasonal pattern characterized by an autumnal peak of abundance regardless the time series period or site (see [START_REF] David | Successful colonization of the calanoid copepod Acartia tonsa in the oligo-mesohaline area of the Gironde estuary (SW France) -Natural or anthropogenic forcing?[END_REF]) and with other studies in many North European estuaries (the Scheldt, [START_REF] Bakker | Comparison of plankton assemblages of identical salinity ranges in estuarine tidal, and stagnant environments. II. Zooplankton[END_REF]; the Ems, [START_REF] Baretta | Distribution and abundance of the zooplankton of the Ems estuary (North Sea)[END_REF]). Although a temporal lag is sometimes noticed between times series of modeled and observed data, no lag was detected in the present study. That explains why some models include a temporal lag to account for a development (generation time) of a species. Our model also estimated a spatial distribution responding to the estuary haline gradient with maximal abundances estimated at the more downstream site. This result was not surprising as A. tonsa was observed in the polyhaline area of the Gironde estuary a long time before it was first recorded in the oligo-mesohaline area in 1983 [START_REF] David | Successful colonization of the calanoid copepod Acartia tonsa in the oligo-mesohaline area of the Gironde estuary (SW France) -Natural or anthropogenic forcing?[END_REF]. This result could be more investigated by a study on the distribution of A. tonsa as a longitudinal shift in the distribution of its congeneric species, A. bifilosa, was already noticed by Chaalali et al. [START_REF] Chaalali | Changes in the distribution of copepods in the Gironde estuary: A warming and marinisation consequence? Estuarine Coastal And Shelf Science[END_REF] during the same period (from 1975 to 2003).

Since the first record of A. tonsa in the Gironde estuary in 1983, species abundance substantially increased [START_REF] David | Successful colonization of the calanoid copepod Acartia tonsa in the oligo-mesohaline area of the Gironde estuary (SW France) -Natural or anthropogenic forcing?[END_REF]. After 1999, this copepod became more abundant than its autochthonous congeneric species, A. bifilosa and even reached abundances comparable with the dominant native species E. affinis [START_REF] David | Successful colonization of the calanoid copepod Acartia tonsa in the oligo-mesohaline area of the Gironde estuary (SW France) -Natural or anthropogenic forcing?[END_REF][START_REF] David | Long-term changes of the zooplankton variability in a turbid environment: The Gironde estuary (France)[END_REF]. The colonization process of the species took place in different phases: from 1978 to 1982, A. tonsa was absent; from 1983 to 1987, species was sporadically observed; from 1988 to 1992, the calanoid was present in low abundance; from 1993 to 1998, the species exponentially increased; and since 1999 high abundances of A. tonsa are regularly observed [START_REF] David | Successful colonization of the calanoid copepod Acartia tonsa in the oligo-mesohaline area of the Gironde estuary (SW France) -Natural or anthropogenic forcing?[END_REF]. Regarding the NPPEN estimations, similar long-term trends were noticeable as we observed an increasing trend in estimated probabilities of occurrence; higher probabilities being estimated toward the end of the time series, reflecting the establishment of the species. However, the NPPEN model did not reproduce/estimate the seasonality of A. tonsa as well as the two empirical models. We suspect that one of the possible explanations could be linked/due to the parameterization of species physiological niche. The limit ranges of temperatures and salinities reported in the literature for species survival were different from the environmental conditions experienced by A. tonsa in the estuary. Indeed, the upper and lower limit values of temperature and salinity measured in the estuary during the study period ranged from 4 to 28uC for temperature (vs. 22 and 42uC reported in the literature as thermal limits; [START_REF] Gonzales | Critical thermal maxima and upper temperature for the calanoid copepods Acartia tonsa and A. clausi[END_REF]) and from 0 to 33 psu for salinity (vs. 0 and 72 psu [START_REF] Cervetto | Influence of salinity on the distribution of Acartia tonsa (Copepoda, Calanoida)[END_REF]).

Finally our results suggest that both warming and marinisation of the Gironde estuary between the beginning (1978)(1979)(1980)(1981)(1982)(1983) and the end (1984-2011) of the studied period, facilitated the establishment of an invasive species, in accordance with Raitsos et al. [START_REF] Raitsos | Global climate change amplifies the entry of tropical species into the Eastern Mediterranean Sea[END_REF] observations. Changes at a seasonal scale provided similar results (results not shown), comforting previous studies already pointing out strong links between temperature or salinity and A. tonsa abundances [START_REF] David | Successful colonization of the calanoid copepod Acartia tonsa in the oligo-mesohaline area of the Gironde estuary (SW France) -Natural or anthropogenic forcing?[END_REF][START_REF] Aravena | Response of Acartia populations to environmental variability and effects of invasive congenerics in the estuary of Bilbao, Bay of Biscay[END_REF]. Our study also highlights it is important to consider coupled/mixed approaches to deal with both species fundamental (i.e. physiological) and realized niches. Nevertheless previous studies (dealing with ecological niches modeling) only focused on realized niches of species (see Hirzel et al. [START_REF] Hirzel | Ecological-niche factor analysis: How to compute habitat-suitability maps without absence data?[END_REF] or Helaoue ¨t & Beaugrand [START_REF] Helaouet | Macroecology of Calanus finmarchicus and C. helgolandicus in the North Atlantic Ocean and adjacent seas[END_REF]), regardeless their fundamental niches. Recent works also pointed out the importance of considering both niches [START_REF] Helaouet | Physiology, ecological niches and species distribution[END_REF]. As an example, Helaoue ¨t & Beaugrand [START_REF] Helaouet | Physiology, ecological niches and species distribution[END_REF] showed a close relationship between Calanus finmarchicus fundamental and realized niches that remained constant at both biogeographical and decadal scales. They observed that changes in environmental forcing propagated from the physiological to the macroecological level and, projected a potential poleward shift in species spatial distribution in the North Atlantic over the 21 st century. Such approaches are relevant to examine the habitat suitability of species (Diekmann et al. [START_REF] Diekmann | Environmental cues and constraints affecting the seasonality of dominant calanoid copepods in brackish, coastal waters: a case study of Acartia, Temora and Eurytemora species in the south-west Baltic[END_REF]) and to improve our knowledge about the ecology of species [START_REF] Helaouet | Physiology, ecological niches and species distribution[END_REF].

In the pursue of our study on a species providing a new food resource for fish and shrimps juveniles [START_REF] Sorbe | La macrofaune vagile de l'estuaire de la Gironde -Distribution et migration des espe `ces: modes de reproduction, re ´gimes alimentaires[END_REF], there is a necessity of further investigations on species future evolution in the Gironde system in response to different climate scenarios [START_REF] Beaugrand | A new model to assess the probability of occurrence of a species, based on presence-only data[END_REF]. Regarding the use of the NPPEN model for invasive species, such studies would need to integrate co-occurring species as invasions could dramatically affect other species [START_REF] Williams | Two degrees of separation in complex food webs[END_REF] and some possible consequences were already observed [START_REF] David | Successful colonization of the calanoid copepod Acartia tonsa in the oligo-mesohaline area of the Gironde estuary (SW France) -Natural or anthropogenic forcing?[END_REF][START_REF] Chaalali | Changes in the distribution of copepods in the Gironde estuary: A warming and marinisation consequence? Estuarine Coastal And Shelf Science[END_REF]. Such predictions have indeed to be carefully considered as non-linear evolutions may occur as some abrupt shifts were previously documented (concerning other species) and make predictions difficult [START_REF] Chaalali | Climate-caused abrupt shifts in a European macrotidal estuary[END_REF][START_REF] Chaalali | Changes in the distribution of copepods in the Gironde estuary: A warming and marinisation consequence? Estuarine Coastal And Shelf Science[END_REF]. Moreover, while species distribution models assume equilibriumor at least pseudo-equilibrium -between the environment and observed occurrence, a non-equilibrium phenomenon is more realistic in ecology. However, such a drawback is compensated here by both our knowledge of the physiology of the species and a large-scale approach [START_REF] Guisan | Predictive habitat distribution models in ecology[END_REF] which allows us to cautiously extrapolate beyond baseline conditions to project future species' range.

Another key issue would be to study the consequences of species colonization on the structure and the function of the trophic food web, a progression of marine fish species such as the sprat (Sprattus sprattus) or the European anchovy (Engraulis encrasicolus) being already documented [START_REF] Pasquaud | Increase of marine juvenile fish abundances in the middle Gironde estuary related to warmer and more saline waters, due to global changes[END_REF]. Indeed, previous studies revealed an increase in the abundance of small pelagics in relation to global climate change (marinisation and water warming; [START_REF] Pasquaud | Increase of marine juvenile fish abundances in the middle Gironde estuary related to warmer and more saline waters, due to global changes[END_REF]). These species are known to predate on copepods and in particular on Acartidae [START_REF] Lobry | Diversity and stability of an estuarine trophic network[END_REF][START_REF] Pasquaud | A preliminary investigation of the fish food web in the Gironde estuary, France, using dietary and stable isotope analyses[END_REF][START_REF] Pasquaud | Determination of fish trophic levels in an estuarine system[END_REF]. For example, the abundance of the European anchovy in the Gironde showed a confounding trend close from the exponential trend noticed on A. tonsa (Delpech personal communication). While our study focused on the importance of climatic facilitation, investigations on the biotic factors that influence both the spatial and temporal colonization of A. tonsa would also be required.

Figure 1 .

 1 Figure 1. Map of the Gironde estuary with the three sampling sites (F, E, and K). doi:10.1371/journal.pone.0074531.g001

Figure 2 .

 2 Figure 2. Observed and modeled ecological niche of Acartia tonsa in the Gironde estuary. A. Observed niche of A. tonsa (log 10 (x+1)) as a function of monthly water temperature and salinity (considering all available data). B. Realized niche of A. tonsa (log 10 (x+1)) as a function of monthly water temperature and salinity. C. Realized niche of A. tonsa (log 10 (x+1)) as a function of monthly water temperature and salinity modeled by a mixed Gausian-linear model. D. Fundamental niche of A. tonsa (occurrence probability) modeled by the NPPEN model. doi:10.1371/journal.pone.0074531.g002

Figure 3 .

 3 Figure 3. Comparison between observed (in red) and estimated (in blue) abundance of Acartia tonsa (log 10 (x+1)) at three sampling stations in the Gironde estuary. A. Station F. B. Station E. C. Station K. Estimated data originated from the realized niche assessed by discretization from monthly water temperature and salinity (see Fig. 2B). doi:10.1371/journal.pone.0074531.g003

Figure 4 .Figure 5 .

 45 Figure 4. Comparison between observed (in red) and modeled (in blue) abundance of Acartia tonsa (log 10 (x+1)) at three sampling stations in the Gironde estuary. A. Station F. B. Station E. C. Station K. Modeled data originated from the realized niche assessed from monthly water temperature and salinity and using a mixed Gausian-linear model (see Fig. 2C). doi:10.1371/journal.pone.0074531.g004

Figure 6 .Figure 7 .

 67 Figure 6. Comparison between observed (in red) and modelled (in blue) abundance of Acartia tonsa (log 10 (x+1)) at three sampling stations in the Gironde estuary. A. Station F. B. Station E. C. Station K. Modeled data originated from the NPPEN model (see Fig. 2D). doi:10.1371/journal.pone.0074531.g006

PLOS ONE | www.plosone.org

September 2013 | Volume 8 | Issue 9 | e74531

Acknowledgments

We express our thanks to staff of the National Institute for the Sciences of the Universe (OASU), the Irstea Bordeaux, and the entire SOMLIT team

The data the authors used for their research article come from monitoring programs supported by Electricite ´de France (EDF), the Institut Franc ¸ais de Recherche pour l'Exploitation de la Mer (IFREMER), and the Centre National de la Recherche Scientifique (CNRS). This study was funded by the Conseil Re ´gional d'Aquitaine. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. The data the authors used for their research article come from monitoring programs supported by Electricite ´de France (EDF), the Institut Franc ¸ais de Recherche pour l'Exploitation de la Mer (IFREMER), and the Centre National de la Recherche Scientifique (CNRS). This study was funded by the Conseil Re ´gional d'Aquitaine. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Author Contributions

Conceived and designed the experiments: BS. Performed the experiments: BS VD AC. Analyzed the data: AC GB VR. Contributed reagents/ materials/analysis tools: BS. Wrote the paper: AC GB VR EG VD PB BS.