M C Corona-Lozada

S Morin

Philippe Choler
email: philippe.choler@univ-grenoble-alpes.fr

Drought offsets the positive effect of summer heat waves on the canopy greenness of mountain grasslands

Keywords: Climatic water balance, European Alps, heat wave, MODIS, mountain grasslands, phenology

Heat waves are becoming more frequent in a warming climate and understanding the impacts of these extreme events on terrestrial ecosystems remains a major challenge. Seasonally snowcovered mountain grasslands are temperature-limited ecosystems and one may question whether the summer heat waves there have the same negative effect on primary productivity as they have on lowland ecosystems. Here, we examined the remotely-sensed phenology of mountain grasslands in the French Alps with a particular focus on the four unusually warm summers that occurred in 2003, 2012, 2015, and 2016. Our results showed that an early and rapid senescence and a lack of lateseason regrowth were the most significant responses in 2003 and 2012 when heat waves coincided with climatic water deficit, i.e. negative anomalies of precipitation minus reference evapotranspiration. Mountain grasslands located below 2350m were the most affected. In contrast, we found low senescence and moderate to high regrowth in 2015 and 2016 which were years unaffected by severe water deficit. When water was non-limiting, warm and even very warm summers have a positive effect on the canopy greenness of the highest grasslands, particularly by extending the green period in late summer and early fall. Our study highlights the pivotal role that summer water balance has in the phenological response of mountain grasslands exposed to heat waves.

Introduction

There has been a marked increase in air temperature in the European Alps from the 1980s onwards and this has caused warmer and longer snow-free periods for mountain vegetation (Durand et al., 2009a;[START_REF] Gobiet | 21st century climate change in the European Alps-A review[END_REF][START_REF] Beniston | The European mountain cryosphere: a review of its current state, trends, and future challenges[END_REF]. A growing body of evidence indicates that above treeline ecosystems have undergone significant changes over the last decades, including changes in community composition [START_REF] Gottfried | Continent-wide response of mountain vegetation to climate change[END_REF][START_REF] Cannone | Ecological responses of plant species and communities to climate warming: upward shift or range filling processes?[END_REF] and ecosystem functioning [START_REF] Choler | Growth response of temperate mountain grasslands to inter-annual variations in snow cover duration[END_REF]. Many of these documented changes, such as the upward shift of species [START_REF] Lenoir | A significant upward shift in plant species optimum elevation during the 20th century[END_REF][START_REF] Steinbauer | Accelerated increase in plant species richness on mountain summits is linked to warming[END_REF], the thermophilization of high summit flora [START_REF] Gottfried | Continent-wide response of mountain vegetation to climate change[END_REF], or the long-term greening of vegetation at high elevation [START_REF] Carlson | Observed long-term greening of alpine vegetation-a case study in the French Alps[END_REF], are consistent with the expected positive effects of rising temperatures in temperature-limited ecosystems. Nonetheless, recent reports indicate that earlier frost exposure due to the earlier snow melt-out date [START_REF] Wheeler | Increased spring freezing vulnerability for alpine shrubs under early snowmelt[END_REF][START_REF] Liu | Extension of the growing season increases vegetation exposure to frost[END_REF], more frequent heat waves [START_REF] Cremonese | Heat wave hinders green wave: The impact of climate extreme on the phenology of a mountain grassland[END_REF] and droughts [START_REF] Marchin | Productivity of an Australian mountain grassland is limited by temperature and dryness despite long growing seasons[END_REF] may be additional effects of on-going climate change which could attenuate this positive response. These studies highlighted the need to consider both trends and extreme events when assessing mountain vegetation response to climate change.

In a warming climate, the duration, intensity, and frequency of heat waves is expected to markedly increase [START_REF] Stocker | Climate change 2013: the physical science basis[END_REF]. The few reports that examine the response of mountain grasslands to summer heat waves have yielded contrasting results. Using both remote sensing and in situ measurements, [START_REF] Jolly | Divergent vegetation growth responses to the 2003 heat wave in the Swiss Alps[END_REF] found a positive effect of the 2003 heat wave on the summer greenness of high-elevation ecosystems in the Swiss Alps, and this was in stark contrast with what was happening in lowland ecosystems. [START_REF] Cremonese | Heat wave hinders green wave: The impact of climate extreme on the phenology of a mountain grassland[END_REF] used flux data and time-lapse cameras to characterize the effect of the 2015 heat wave on a subalpine Nardus-rich community and noted a reduction of canopy development in addition to advanced senescence. Experimental approaches using monolith transplantation have highlighted that the negative effect of an increased air temperature only occurred when it coincided with reduced water supply [START_REF] De Boeck | Simulated heat waves affected alpine grassland only in combination with drought[END_REF].

This calls for a comparative study that includes a broader range of temperature-precipitation combinations to better capture the interplay of heat waves and drought on mountain grassland phenology.

The remote sensing of vegetation phenology provides a unique opportunity to track ecosystem responses to climate change at a regional scale and over decades [START_REF] Pettorelli | Using the satellite-derived NDVI to assess ecological responses to environmental change[END_REF][START_REF] Richardson | Climate change, phenology, and phenological control of vegetation feedbacks to the climate system[END_REF]. In mountainous environments, several studies have successfully used high or medium resolution satellite images to examine the relationships between climate and the growth response of mountain grasslands [START_REF] Fontana | Alpine grassland phenology as seen in AVHRR, VEGETATION, and MODIS NDVI time series -a comparison with in situ measurements[END_REF][START_REF] Dunn | Land surface phenology of North American mountain environments using moderate resolution imaging spectroradiometer data[END_REF]Zhang et al., 2013b;[START_REF] Choler | Growth response of temperate mountain grasslands to inter-annual variations in snow cover duration[END_REF]. For seasonally snow-covered ecosystems, the use of snow and vegetation indices may provide several phenological metrics including the onset of growth, the date and magnitude of maximum canopy greenness, the dynamics of senescence, and the date of the first snowfalls. However, many comparative studies (either inter-annual or inter-spatial comparisons) have focused on the maximum or averaged summer greenness and, by doing so, have largely overlooked the late and end of season phenology [START_REF] Gallinat | Autumn, the neglected season in climate change research[END_REF][START_REF] Garonna | Shifting relative importance of climatic constraints on land surface phenology[END_REF]). Yet extreme events occurring after the date of peak standing biomass, such as heat wave or drought, may accelerate the senescence or diminish the possibility of regrowth and this has important consequences for the management of mountain pastures [START_REF] Nettier | Resilience as a framework for analyzing the adaptation of mountain summer pasture systems to climate change[END_REF]. We believe that remotely-sensed phenological metrics that capture the kinetics of both growth and senescence phases will provide a more comprehensive assessment of mountain grassland responses to heat waves.

In this paper, we used time series of the Normalized Difference Vegetation Index (NDVI) for the period 2000 -2016 to examine how the phenology of mountain grasslands has been impacted by the summer heat waves in 2003, 2012, 2015 and 2016. We performed our comparative analysis by using the NDVI from Moderate-Resolution Imaging Spectroradiometer (MODIS) and meteorological variables derived from the SAFRAN (Système d'Analyse Fournissant des Renseignements Atmosphériques à la Neige) reanalysis system (Durand et al., 2009a(Durand et al., , 2009b)). We focused on the French Alps for which the availability of high-resolution land cover data allowed us to select MODIS pixels with predominant coverage of mountain grasslands. More specifically, we addressed the following three sets of questions: (1) what was the intensity, duration, and timing of summer heat waves in the French Alps from 2000 to 2016 and to what extent did these events coincide with climatic water deficit, ( 2) what are the remotely-sensed indicators that best capture the impact of extreme meteorological events on the phenology of mountain grasslands, and (3) how do heat wave and water balance interfere on the growth and senescence of mountain grasslands?

Materials and Methods

The distribution of mountain grasslands in the French Alps was derived from a 10 m resolution land cover map of France produced by the Centre d'Expertise Scientifique Occupation des sols, downloadable at http://osr-cesbio.ups-tlse.fr/~oso/ [START_REF] Inglada | Operational High Resolution Land Cover Map Production at the Country Scale Using Satellite Image Time Series[END_REF]. To obtain a binary map of mountain grasslands, we aggregated the three land cover classes corresponding to meadows ("pelouses", code 34), heathlands ("landes ligneuses", code 36) and grasslands ("prairies", code 211). We noticed that the class "landes ligneuses" included a limited number of pixels and was not very reliable due to the difficulty in distinguishing mountain low shrubs from nearby grasslands and meadows. We calculated the percent cover of the aggregated class at a 250 m resolution to align with the MODIS data, selecting those pixels with >90% cover and a mean elevation ranging from 1350 m to 2550 m m a.s.l.. This yielded a total of 44,295 pixels which fall into four elevation bands (1350 m -1650 m, 1650 m -1950 m, 1950 m -2250 m, and 2250 m -2550 m) and 23 SAFRAN massifs (Fig. 1 and Table S1). Pixel-based estimates were averaged by elevation band and massif. This resulted in a table containing 87 combinations of elevation band and massif and 17 years (2000 -2016). Annual standardized anomalies were calculated from this table.

For meteorological data, we used the SAFRAN reanalysis that simulates surface conditions as a function of elevation in 300 m increments and for 23 massifs of French Alps. These massifs were defined for their climatological homogeneity in a previous spatial analysis of all weather data available for the French Alps (Durand et al., 2009a;2009b). The SAFRAN reanalysis was evaluated against in-situ observations by Durand et al. (2009a) and is used for many real-time and climatological applications in the French mountain areas (e.g. [START_REF] Verfaillie | Multi-component ensembles of future meteorological and natural snow conditions for 1500m altitude in the Chartreuse mountain range, Northern French Alps[END_REF]. The period from 1981 to 2010 was set as a reference period to estimate standardized anomalies. For each MODIS pixel, we extracted all the required meteorological data to compute a Heat Wave Magnitude index, hereafter HWI, and a Reference Crop Evapotranspiration, hereafter ET 0 as described below.. We simplified the methodology proposed by [START_REF] Russo | Magnitude of extreme heat waves in present climate and their projection in a warming world[END_REF] to calculate a HWI that accounts for both magnitude and duration of heat waves. First, we identified days exhibiting a maximum temperature above a threshold value. For each day of year (DOY) i, we assembled a dataset (DS i ) comprising the daily Tmax values of the n days before and after day i for the 30 yearlong reference period according to Eq. 1 124 125 where the symbol U denotes the union of values. The 80 th percentile of DS i defined the threshold value for day i. As in [START_REF] Russo | Magnitude of extreme heat waves in present climate and their projection in a warming world[END_REF], we chose n = 15 days to estimate the threshold.

We examined the sensitivity of HWI to varying values of n. The relative change in HWI was predominantly less than 5% when comparing n=15 with n=13, and less than 10% when comparing n=15 with n=11 (Fig. S1). In all cases, the ranking of years was the same and the main findings of the study did not change. Second, we selected all episodes of at least three consecutive days during for which the maximum air temperature exceeded the threshold and summed the daily maximum temperature of each episode to obtain an unscaled magnitude of the heat wave. Finally, we summed the unscaled magnitudes of all heat waves occurring during the three periods of interest (see below).. [START_REF] Russo | Magnitude of extreme heat waves in present climate and their projection in a warming world[END_REF], we did not consider sub-heat waves and did not rescale heat wave magnitude into probability values.

Compared to

⋃ ⋃

(1) ET 0 is a Penman-Monteith estimate of potential evapotranspiration of a reference crop surface that quantifies the evaporating demand of the atmosphere when water is not limiting. We followed the recommendations of the Food and Agricultural Organization [START_REF] Allen | Crop evapotranspiration-Guidelines for computing crop water requirements-FAO Irrigation and drainage paper 56[END_REF] to compute ET 0 (in mm day -1 ) from weather data (Eq. 2).

( ) ( ) ( ) ( ) (2)
where R n is the net radiation at crop surface (MJ m -2 day -1 ), G is the soil heat flux density (MJ m -2 day -1 ), Tm is mean daily air temperature at 2 m (°C), u 2 is the wind speed at 2 m (m s -1 ), e s is saturation vapour pressure (kPa), e a is the actual vapour pressure (kPa), Δ is the slope of vapour pressure curve (kPa °C-1 ), and γ is the psychrometric constant (kPa °C-1 ). The use of SAFRAN reanalysis to compute ET 0 according to Eq. 2 is described in [START_REF] Vannier | Calcul d'une évapotranspiration de référence spatialisée pour la modélisation hydrologique à partir des données de la réeanalyse SAFRAN de Météo-France[END_REF]. The cliamtic water balance corresponding to the difference between precipitation (P) and reference crop evapotranspiration, hereafter P-ET 0 , was taken as a proxy to assess the dry-wet conditions of a given period.

For each year, HWI and P -ET 0 were computed over three periods of interest: (i) the climatological summer (June, July and August), hereafter JJA, (ii) the growth phase defined as the 30 days preceding the date of maximum observed NDVI, hereafter B30, and (iii) the senescent phase defined as the 60 days following the date of maximum observed NDVI, hereafter A60 (Fig. 2). We examined the sensitivity of the results to the length of the latter two periods and did not find any significant effect if periods were increased or decreased by 10 to 15 days.

Our remote sensing analysis was based on the surface reflectance acquired by MODIS, a sensor onboard Terra satellite. We downloaded the collection 6 of the MOD09Q1 products, corresponding to tile h18.v4 (40-50• N, 0-15.6• E), from the Land Processes Distributed Active

Archive Center (https://e4ftl01.cr.usgs.gov/). The 8-day composite at 250 m resolution, covering the period from 18 February 2000 to 27 December 2016, was processed following [START_REF] Choler | Growth response of temperate mountain grasslands to inter-annual variations in snow cover duration[END_REF]. To summarize, re-projected surface reflectance values (ρ) in the red and near-infrared (NIR) bandwidths were used to calculate a NDVI according to Eq. 3, keeping only reflectance values of optimum quality (according to the MOD09Q1 Quality Control flag). Following [START_REF] Choler | A simple ecohydrological model captures essentials of seasonal leaf dynamics in semi-arid tropical grasslands[END_REF], we first removed outliers defined as NDVI values that were two times larger or smaller than the average value of the two preceding and the two following observations. Such sudden changes in NDVI are inconsistent with the growth or senescence dynamics of mountain grasslands and may be due to quality control issues, ephemeral snowfall during summer etc. Only 0.2% of the total NDVI values were discarded at this step. Missing values and removed outliers were gap-filled using a cubic spline interpolation. We then applied a Savitzky and Golay low pass filter in which the original data value at time t i was replaced by the value of a cubic polynomial fitted for 2n+1 points, with n being the number of time steps preceding or following t i [START_REF] Savitzky | Smoothing and Differentiation of Data by Simplified Least Squares Procedures[END_REF]. Consistent with our previous studies [START_REF] Choler | A simple ecohydrological model captures essentials of seasonal leaf dynamics in semi-arid tropical grasslands[END_REF][START_REF] Choler | Growth response of temperate mountain grasslands to inter-annual variations in snow cover duration[END_REF], we empirically found that a window size of seven days (i.e. n=3) was able to track ecologically significant changes in NDVI while still able to reduce noise in the NDVI time series efficiently. Finally, we used cubic splines to daily interpolate the NDVI time series and to derive the following phenological metrics (Fig. 2): (i) the growth onset, hereafter Onset, defined as the Julian day on which NDVI reached half its maximum, (ii) the NDVI maximum, hereafter NDVImax, as a proxy for the maximum canopy greenness, (iii) the Julian day of NDVImax, hereafter TNDVImax, as the date of the peak vegetation (iv) the absolute NDVI decrease during the 60 days following NDVImax, hereafter SEN (A60), see Eq. 4 as a proxy for the canopy senescence, and (v) the absolute NDVI increase during the 60 days following NDVImax, hereafter REG (A60), see Eq. 5 as a proxy for the end-of-season regrowth..

(3)

∑ ( ( )) (4) ∑ ( ) (5)
Pixel-based estimates of phenology were averaged by elevation band and massif in order to align with meteorological variables. For REG (A60) we transformed the raw values into binary data (zero or positive values) and we calculated the percent of pixels exhibiting regrowth for each combination of elevation band and massif.

Standard Major Axis (SMA) regression analysis was used to analyse the bivariate relationships between anomalies of phenology and to test the statistical significance of any deviation of the slope from unity and of the intercept from zero [START_REF] Warton | smatr: (Standardised) Major Axis Estimation and Testing Routines[END_REF]. We modelled anomalies of NDVImax, SEN (A60) and REG (A60) as a linear function of anomalies of HWI and P-ET 0 . The model was written as follows (Eq. 6)

where Y is the dependent variable (phenology),  are the regression coefficients, anom means standardized anomaly and  is an error term. The  parameters were solved separately for each elevation band (z) using ordinary least squares regression.

All statistical analyses were performed in R (R Core Team, 2017), using the signal library (signal developers, 2013) for the time series analyses and the raster library [START_REF] Hijmans | raster: Geographic Data Analysis and Modeling[END_REF] for the spatial analysis. All figures were produced using the ggplot2 [START_REF] Wickham | ggplot2: Elegant Graphics for Data Analysis[END_REF] and the rasterVis (Perpinan Lamigueiro and Hijmans, 2018) libraries.

( ) ( ) ( ) (6)

Results

Since the 1980s, warm to very warm summers with above normal HWI have become more frequent in the French Alps (Fig. 3A). The four summers with the highest HWI all occurred in the last 15 years, i.e. during the years 2003, 2012, 2015 and 2016 (Fig. 3A). The intensity of summer heat waves reached a maximum in 2003 with an HWI (JJA) median value of 184 for the 1950 m -2250 m elevation band (Fig. S2A). During the same period, we found a wide range of water balance conditions during the summer period, with years when high values of HWI and negative P-ET 0 values co-occurred and other years when they occurred independently (Fig. 3, Fig. S2B). Summer droughts coincided with the 2003 and 2012 heat waves, while summer water balance was close to normal in 2015 and 2016 (Fig. 3B). The French Alps also experienced a severe drought in 2009, but this was not accompanied by a significant heat wave (Fig. 3B). m elevation band (mean DOY of 157) (Table 1). The date of maximum canopy greenness followed around 1.5 months later, with the length of the growing period (TNDVImax-ONSET) decreasing with elevation (Table 1). Late season phenology demonstrated a much larger range of inter-annual variabilitywith a coefficient of variation ranging from 18% to 31% for senescence and regrowthas compared to those of Onset, TNDVImax and NDVImax (Table 1, Fig. 5), and this was observed across all elevation bands (Table 1). Overall, year-to-year variations in phenology were much smaller than those of HWI and P-ET 0 (Table 1).

The last four years (2013-2016) were among the five greenest summers of the study period, while the four strongest negative anomalies of NDVImax occurred before 2006 (Fig. 5D). The summers of 2003 and 2012 were characterized by an above normal senescence and a lack of regrowth, although the NDVImax was close to normal (Fig. 5D,E,F). A lack of regrowth was also observed in 2016, although senescence was close to normal and NDVImax was particularly high (Fig. 5D,E,F). By comparison, the year 2015 exhibited a below normal senescence and above normal regrowth and NDVImax (Fig. 5D,E,F).

Figure 6 shows the bivariate relationships among phenological anomalies and illustrates the particularity of the four years marked by heat wave occurrence. We found a strong positive relationship between anomalies of Onset and of TNDVImax (Table 2, Fig. 6A). The SMA slope was significantly below one indicating that maximum canopy greenness was reached in a shorter time in case of a delayed onset (Table 2). In addition, the two years with early season drought (e.g., 2003

and 2012) exhibited an earlier peak of NDVI than the one predicted by the SMA (Fig. 6A). The relationship between TNDVImax and NDVImax showed that a delayed phenology tend to be associated with higher values of NDVImax (Table 2, Fig. 6B). Remarkably, the two years with a pronounced early season heat wave (e.g., 2003 and 2015) exhibited a higher anomaly of NDVImax than predicted (Fig. 6B). Finally, we found a significant negative relationship between SEN (A60)

and REG (A60) (Table 2, Fig. 6C). The years 2003, 2012 and 2016 exhibited a late-season phenology that was marked by normal to rapid senescence and limited regrowth while the year 2015 clearly showed an opposite pattern (Fig. 6C).

Finally, we examined the effects of HWI and P-ET 0 on phenology using multiple regression models. For this quantitative analysis, we distinguished between the two periods, i.e. before and after the peak, and among different elevation bands. Irrespective of elevation, there was a significant effect of P-ET 0 (B30) anomalies on NDVImax while no effect was detected for HWI (B30) (Table 3, Fig. 7A,B). Wetter conditions in the early season significantly enhanced the growth of grassland while heat waves did not. P-ET 0 after the peak was also a significant driver of SEN (A60) and REG (A60) and the magnitude of this effect markedly differed among elevation bands. Wet summers significantly slowed down senescence and enhanced regrowth at low elevation, with up to 47% of variance explained in SEN (A60) and 24% of variance explained in REG (A60) in the 1350 m 1650 m elevation band (Table 3, Fig. 7C,E). These effects vanished higher up and became non-significant (SEN) or marginally significant (REG) above 2250 m (Table 3). Regarding HWI, we only found a significant and negative effect on senescence above 2250 m (Table 3). Noticeably, the effect of lateseason heat waves on senescence depended on elevation with an increased senescence at low elevation (although not significantly) and a delayed senescence at high elevation (Fig. 7D).

Overall, these results indicated that maximum canopy greenness, senescence and regrowth were much more responsive to summer water balance than to heat waves. This is consistent with our initial hypotheses concerning the combined effects of heat wave and climatic water deficit on plant phenology. An early heat wave that did not coincide with a severe drought -as in 2015 and 2016did not significantly impact NDVImax. A late heat wave coinciding with a severe drought -as in 2003, 2012 and 2016 -resulted in an increased senescence and a low regrowth. A late heat wave not accompanied by a pronounced droughtas in 2015 -induced a delayed senescence, especially at high elevation, and a substantial regrowth.

Discussion

Our study of mountain grassland phenology in the French Alps for the period from 2000 to 2016 has provided the framework to analyse more thoroughly how these ecosystems respond to heat waves. Using data from the past two decades, we could show that (i) warmer summers generally have a positive effect on mountain canopy greenness especially at high elevation, (ii) water deficit hinders the positive temperature effect of heat waves on canopy greenness, and (iii) faster senescence and weaker late-season regrowth are the main negative effects of heat waves that co-occur with water deficit. Our findings emphasize the need to consider both the growth and senescent phases to capture the effect of climate extremes on phenology. We also highlighted the pivotal role of the magnitude of summer water balance and its within-summer distribution for understanding whether mountain grasslands can benefit or not from higher summer temperatures.

We found a general increase in canopy greenness (NDVImax) in the French Alps during the study period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) which is consistent with our previous remote sensing study that showed a greening trend for high-elevation ecosystems in the Ecrins Massif [START_REF] Carlson | Observed long-term greening of alpine vegetation-a case study in the French Alps[END_REF]. This greening has also been documented in other temperate mountains [START_REF] Xu | High positive correlation between soil temperature and NDVI from 1982 to 2006 in alpine meadow of the Three-River Source Region on the Qinghai-Tibetan Plateau[END_REF] as well as in the arctic and sub-arctic [START_REF] Jia | Greening of arctic Alaska, 1981-2001[END_REF][START_REF] Forbes | Russian Arctic warming and 'greening' are closely tracked by tundra shrub willows[END_REF]. The combination of earlier snow melt-out date and higher summer temperatures has been put forward as the main driver of changes in primary productivity [START_REF] Jonas | How alpine plant growth is linked to snow cover and climate variability[END_REF]Zhang et al., 2013a;[START_REF] Choler | Growth response of temperate mountain grasslands to inter-annual variations in snow cover duration[END_REF][START_REF] Dedieu | On the Importance of High-Resolution Time Series of Optical Imagery for Quantifying the Effects of Snow Cover Duration on Alpine Plant Habitat[END_REF][START_REF] Choler | Winter soil temperature dependence of alpine plant distribution: Implications for anticipating vegetation changes under a warming climate[END_REF]. Besides temperature, other factors such as the increased deposition of atmospheric nitrogen possibly led to some eutrophication of upland ecosystems and further contributed to greening [START_REF] Bowman | Limited ecosystem recovery from simulated chronic nitrogen deposition[END_REF]. In alpine plants, nitrogen addition may induce stronger phenological changes than warming alone [START_REF] Smith | Phenological Changes in Alpine Plants in Response to Increased Snowpack, Temperature, and Nitrogen[END_REF]. Our study indicated that sparsely vegetated meadows in the French Alps are the most prone to greening, and we hypothesized that the prolonged growing season in high-elevation contexts has created new thermally suitable niches for early plant colonizers which could lead to an increase in primary productivity [START_REF] Carlson | Observed long-term greening of alpine vegetation-a case study in the French Alps[END_REF]. These conclusions are in agreement with ground-based observations of community dynamics that emphasise the densification of plant cover in open grasslands and/or the upward/northward shift of more productive functional groups [START_REF] Kullman | Alpine flora dynamics -a critical review of responses to climate change in the Swedish Scandes since the early 1950s[END_REF][START_REF] Myers-Smith | Shrub expansion in tundra ecosystems: dynamics, impacts and research priorities[END_REF][START_REF] Steinbauer | Accelerated increase in plant species richness on mountain summits is linked to warming[END_REF].

Our results also showed that the frequency of summer heat waves has increased over the last decades in the French Alps. This is consistent with reports highlighting the increased probability of heat waves in a warming climate [START_REF] Stocker | Climate change 2013: the physical science basis[END_REF][START_REF] Russo | Magnitude of extreme heat waves in present climate and their projection in a warming world[END_REF]. We pointed out that the four main heat waves that hit the Alps over the last 20 years did not have the same effect on mountain grassland phenology. Remarkably, the 2015 heat wave had opposite late-season effects to those of 2003, 2012, and 2016. This discrepancy is due to the combined effects of temperature and climatic water balance, both in magnitude and timing. Water deficit coincided with heat wave in 2003 and 2012, and to a lesser extent in 2016. We hypothesize that a threshold in the water balance was reached for these three years leading to the diminished canopy greenness during the late-season.

Other reports have suggested that the disproportionate effect of the 2003 summer heat wave may not only be due to the unusually warm temperatures, but also to a negative soil water anomaly at the beginning of the growth period [START_REF] Loew | The European heat wave 2003: Early indicators from multisensoral microwave remote sensing[END_REF]. Compared to the lowland ecosystems, mountain grasslands did not exhibit a strong reduction of NDVImax in 2003 and 2012, probably because the snowmelt provided sufficient water for the initial growth period, at least in the Northern Alps.

Subsequently, the soil water content has likely decreased during the summer due to a climatic water deficit causing a significant reduction of canopy greenness during the late-season. For the year 2015, we hypothesize that the above normal water balance in the early season and the limited water deficit in the late-season have favoured the positive effect of the late-season heat wave, especially for the highest, i.e. the most temperature-limited, mountain grasslands. In their study of the 2003 heat wave, [START_REF] Jolly | Divergent vegetation growth responses to the 2003 heat wave in the Swiss Alps[END_REF] reported marked regional differences in the response of high elevation ecosystems.

In particular, it was striking that mountain ecosystems of the Southern Alps showed a reduced canopy greenness compared to the Northern French or Swiss Alps, which was likely related to the interplay of heat wave and water availability. Similarly, [START_REF] Cremonese | Heat wave hinders green wave: The impact of climate extreme on the phenology of a mountain grassland[END_REF] pointed out the colimitation of high temperature and soil water content to impede the greening of a Nardus-rich grassland in 2015. Recent experimental studies conducted on temperate mountain grasslands showed that the direct effects of air temperature on photosynthesis and aboveground phytomass production strongly depend on soil water content [START_REF] De Boeck | Simulated heat waves affected alpine grassland only in combination with drought[END_REF]. These findings should prompt us to improve our understanding and regional modelling of the soil water balance of mountain grasslands.

Unfortunately, this remains notoriously difficult due to the many unknowns (e.g., the soil water capacity and lateral run-offs). However, with further temperature increases looming in the not so distant future, there is an urgent need to formulate eco-hydrological models of seasonal leaf dynamics and to identify the different thresholds that drive ecosystem response. This is particularly true for grasslands of the subalpine belt that are the most sensitive to drought and are a key resource for summer livestock farming [START_REF] Nettier | Resilience as a framework for analyzing the adaptation of mountain summer pasture systems to climate change[END_REF] Numerous studies on arctic and alpine ecosystems have underlined that the response to warming varies across functional groups [START_REF] Arft | Responses of tundra plants to experimental warming: Meta-analysis of the international tundra experiment[END_REF][START_REF] Jonas | How alpine plant growth is linked to snow cover and climate variability[END_REF][START_REF] Jung | Intraspecific trait variability mediates the response of subalpine grassland communities to extreme drought events[END_REF][START_REF] Livensperger | Earlier snowmelt and warming lead to earlier but not necessarily more plant growth[END_REF]. In relation to the recent observed greening it is plausible that the expansion of more productive habitats has changed the responsiveness of mountain grasslands to climate extremes. The coexistence of a richer set of plant life forms may arise, e.g., through the migration of low shrubs onto the grasslands, or of tall perennial grasses into otherwise low-stature habitats. There are several lines of evidence that an increased functional diversity may provide some ecosystem insurance by buffering the impact of environmental changes [START_REF] Diaz | Vive la difference: plant functional diversity matters to ecosystem processes[END_REF].

However, further studies are needed to examine a possible relationship between plant functional diversity and the phenological responses of mountain grasslands to heat waves and drought. For example, the combined use of high-resolution satellite imagery in combination with a regional-scale database of vegetation could provide new insights into this line of enquiry.

Our comparative analysis of mountain grassland phenology over the last two decades provides a knowledge base to help anticipate the long-term response of grasslands under future climate conditions. Regional climate models used to compute climate projections are operating at coarse resolutions (10-50 km) [START_REF] Jacob | EURO-CORDEX: new high-resolution climate change projections for European impact research[END_REF] and are insufficient to provide ecologicallyrelevant surface parameters in mountainous regions [START_REF] Ford | Spatial Heterogeneity in Ecologically Important Climate Variables at Coarse and Fine Scales in a High-Snow Mountain Landscape[END_REF]. Although downscaling methods can refine projections of meteorological variables and snow cover dynamics [START_REF] Verfaillie | Multi-component ensembles of future meteorological and natural snow conditions for 1500m altitude in the Chartreuse mountain range, Northern French Alps[END_REF], further work is needed to translate these model outputs into bioclimatic indicators and assess to what extent they are driving canopy greenness. In the longer term, efforts have to be conducted towards prognostic models of mountain grassland phenology that capture both temperature and soil water constraints. The dynamics of leaf growth and senescence remains poorly represented in land surface models [START_REF] Arora | A parameterization of leaf phenology for the terrestrial ecosystem component of climate models[END_REF][START_REF] Richardson | Climate change, phenology, and phenological control of vegetation feedbacks to the climate system[END_REF], and this holds particularly true for seasonally snow-covered mountain grasslands. Our study has focused on interannual fluctuations to provide a broad picture over the French Alps. Further work examining between-site phenological variation along elevational and bioclimatic gradients will provide a complementary view and will be helpful to improve our representation of the feedback between phenology and climate within process-oriented models.

Conclusions

Our study highlighted the combined effect of summer temperature and water balance on the phenology of temperate mountain grasslands. In the European Alps, the increased or prolonged canopy greenness of mountain grasslands that result from warmer conditions becomes less pronounced when water is limiting. While there has been a significant decadal-scale increase in both the summer temperature and the frequency of heat waves, summer precipitation and reference evapotranspiration did not show a similar trend and exhibited large inter-annual variability. In addition, there is a growing body of evidence to support long-term changes of mountain vegetation in relation to climate warming and land use change. The extent to which these changes modify the sensitivity of mountain grassland to extreme events will also be pivotal to make possible the prediction of canopy greenness under future climate. Grey boxplots represent the four main summer heat waves (2003, 2012, 2015, and 2016). and the first and third quartiles, with whiskers showing the 95% confidence intervals. Grey boxplots represent the four main summer heat waves (2003, 2012, 2015, and 2016). third quartiles, with whiskers showing the 95% confidence intervals. Grey boxplots represent the four main summer heat waves (2003, 2012, 2015, and 2016).

Figure 4

 4 Figure4showed that the 2012 and 2016 heat waves mainly concerned the late season while

Fig

  Fig. S2.

Figure 2 .

 2 Figure 2. An idealized yearly course of vegetation index showing the different variables and periods

Figure 3 .

 3 Figure 3. Annual standardized anomalies of (A) the Heat Wave magnitude Index (HWI (JJA)) and

Figure 4 .

 4 Figure 4. Bivariate plots of standardized anomalies of HWI and P-ET 0 during (A) the 30 days before

Figure 5 .

 5 Figure 5. Annual standardized anomalies of mountain grassland phenology: (A) onset of growth, (B)

Figure 6 .

 6 Figure 6. Bivariate relationships between annual standardized anomalies of TNDVImax, NDVImax,

Figure 7 .

 7 Figure 7. Effect of HWI and P-ET 0 on: NDVImax (A, B), SEN (A60) (C, D) and REG (A60) (E, F).

  Figure1

Table 1 .

 1 Annual means and inter-annual coefficients of variation (CV in %) of the meteorological 383 and phenological variables used in the study. 384

					Mean (CV)
				1350 m -	1650 m -	1950 m -	2250 m -
	Metric	Period	Abbreviation	1650 m	1950 m	2250 m	2550 m
	Meteorology
		Summer	HWI (JJA)	55 (79)	65 (71)	58 (79)	45 (89)
	Heat Wave magnitude Index	Preceding the peak	HWI (B30)	18 (76)	21 (77)	20 (107)	14 (127)
		Following the peak	HWI (A60)	35 (111)	43 (98)	34 (116)	22 (130)
	Precipitation -	Summer	P-ET 0 (JJA)	61 (79)	58 (110)	51 (24)	16 (22)
	Reference Crop	Preceding the peak	P-ET 0 (B30)	35 (150)	24 (135)	14 (98)	5 (-285)
	Evapotranspiration (mm)	Following the peak	P-ET 0 (A60)	32 (25)	36 (-41)	39 (11)	25 (103)
	Phenology
	Onset of growth (DOY)		ONSET	112 (8)	130 (8)	145 (8)	157 (7)
	Date of maximum NDVI (DOY)		TNDVImax	175 (7)	184 (5)	192 (5)	200 (6)
	Offset of growth (DOY)		OFFSET	327 (5)	315 (5)	306 (4)	300 (3)
	Maximum NDVI		NDVImax	0.79 (2)	0.78 (2)	0.75 (3)	0.69 (4)
	Senescence	Following the peak	SEN (A60)	0.13 (22) 0.14 (18) 0.16 (18) 0.19 (21)
	Regrowth	Following the peak	REG (A60)	69 (18)	62 (19)	54 (24)	48 (31)
	385

Table 2 .

 2 Summary statistics for Standardized Major Axis (SMA) regression analysis of phenological variables. Shown are correlation coefficients (r), estimates of slopes (a) and intercepts (b) with 95% confidence intervals, Mean Square Error (MSE) and Mean Absolute Error (MAE). For the intercept,

	(NS) indicates a non-significant offset. For the slope, (NS) indicates an estimate statistically
	indistinguishable from one.
				SMA	H 0 : b=0	SMA	H 0 : a=1
	Relation	n	r	intercept (b)	P value	slope (a)	P value	MSE MAE
	ONSET -TNDVImax	1479	0.69	-0.11	***	0.73	***	0.39	0.50
				(-0.13, -0.08)		(0.7, 0.76)
	TNDVImax -NDVImax	1479	0.31	0.04	*	0.97	NS	0.51	0.55
				(0.01, 0.08)		(0.93, 1.02)
	SEN (A60) -REG (A60)	1478 -0.45	-0.1	***	-1.82	***	1.70	1.06
				(-0.16, -0.05)		(-1.74, -1.91)

Table 3 .

 3 Regression coefficients and adjusted R 2 for the linear models testing the effect of P-ET 0 and HWI on NDVImax, SEN (A60) and REG (A60). (x) denotes the interaction between the two effects.Plain text is for highly significant effect (P < 0.001), italics is for significant effect (P < 0.01). Empty cells indicate non-significant effects. Bauges; 5. Beaufortin; 6. Haute-Tarentaise; 7. Chartreuse; 8. Belledonne; 9. Maurienne; 10. Vanoise; 11. Haute-Maurienne; 12. Grandes Rousses; 13. Thabor; 14. Vercors; 15. Oisans; 16. Pelvoux; 17. Queyras; 18. Devoluy; 19. Champsaur; 20. Embrunnais Parpaillon; 21. Ubaye; 22. Haut-Var Haut-Verdon; 23. Mercantour. Account of pixels by massif and elevation band is given in

	NDVImax

Acknowledgements

Thanks are due to Brad Carlson for helpful discussions and to Deborah Verfaillie, Matthieu Lafaysse and Matthieu Vernay (CNRM/CEN) for maintenance and provision of the SAFRAN reanalysis used in this study. Comments and suggestions from the two reviewers have greatly improved the manuscript. This research was conducted within the Long-Term Socio-Ecological Research (LTSER) Zone Atelier Alpes, a member of the eLTER-Europe network. It received funding from the LTSER Zone Atelier Alpes and from a scholarship of CONACYT. LECA and CNRM/CEN are part of Labex OSUG@2020 (ANR10 LABX56).


Rym Lalaoui

Ana Djukovic

Sofiane Bakour

Jaime Sanz

Eva M Gonzalez-Barbera

Miguel Salavert

Jose Luis López-Hontangas

Miguel A Sanz

Karina B Xavier

Bernhard Kuster

Laurent Debrauwer

Carles Ubeda

Jean-Marc Rolain
email: jean-marc.rolain@univ-amu.fr

Detection of plasmid-mediated colistin resistance, mcr-1 gene, in Escherichia coli isolated from high-risk patients with acute leukemia in Spain

Keywords: Escherichia coli, colistin resistance, mcr-1 gene, leukemic patients, Spain

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Acute leukemia patients are considered a high-risk group because of their weakened immune system. Infections with multidrug-resistant (MDR) bacteria are highly responsive in this type of patients, mainly due to the combination of several risk factors such as: hematological disease, intensive/or repeated chemotherapy, neutropenia, healthcare-associated infections, gastrointestinal mucositis and prolonged hospitalization, which promote their colonization by this type of bacteria [1,[START_REF] Alp | Antibacterial resistance in patients with hematopoietic stem cell transplantation[END_REF].

Antimicrobial therapy in hematology patients, such as leukemia patients, is often used for its important contribution to the survival of these patients, but the emergence of MDR bacteria due to selection pressure complicates the management of these bacterial infections [1].

Colistin, an antibiotic long abandoned for its neurological and renal toxicity, has been reintroduced for its effectiveness against MDR Gram-negative bacteria, especially against carbapenemase producers [START_REF] Olaitan | Emergence of polymyxin resistance in Gram-negative bacteria[END_REF][START_REF] Baron | Molecular mechanisms of polymyxin resistance: knowns and unknowns[END_REF]. Indeed, colistin used alone or in combination with other antibiotics, has shown its effectiveness in the treatment of certain bacteremia due to MDR bacteria in hematology patients [1]. Unfortunately, since its use, colistin resistance has increased considerably, represented mainly by chromosomal gene mutations involving a variety of lipopolysaccharide (LPS) modifications [START_REF] Baron | Molecular mechanisms of polymyxin resistance: knowns and unknowns[END_REF][START_REF] Olaitan | Mechanisms of polymyxin resistance: acquired and intrinsic resistance in bacteria[END_REF]. In 2016, Liu et al. reported for the first time colistin resistance mediated by mobile genetic elements identified in Enterobacteriaceae, called plasmid-mediated colistin resistance gene mcr-1 [START_REF] Liu | Emergence of plasmidmediated colistin resistance mechanism MCR-1 in animals and human beings in China: a microbiological and molecular biological study[END_REF]. Since it first detection, mcr-1 gene was widespread worldwide in both animals and humans, and several variants of this gene were detected in Enterobacteriaceae [START_REF] Wang | Emergence of a novel mobile colistin resistance gene, mcr-8, in NDM-producing Klebsiella pneumoniae[END_REF]. Very few data on PubMed are available on the occurrence of mcr genes in leukemic patients and no studies have been conducted in Spain. For this reason, we have sought to detect the presence of these genes in fecal samples collected from patients with acute leukemia in a single institution in Spain.

Materials and methods

Study design

The FloraStopMRE project (2015 Infect-ERA call) is a collaboration of a multidisciplinary consortium of scientists aiming to understand the role of the human gut microbiome in conferring protection against MDR Enterobacteriaceae infections. Between November 2013 and April 2015, a total of 802 fecal samples were collected from 133 patients with acute leukemia at the University Hospital La Fe (Valencia, Spain). All subjects gave their informed consent for inclusion before they participated in the study. The study was approved on the 1st of July 2013 by the Ethics Committee of CEIC Dirección General de Salud Pública y Centro Superior de Investigación en Salud Pública (20130515/08). Samples were collected every week during their hospitalization period. These samples were then screened for the presence of MDR Enterobacteriaceae (MRE) by plating them on Brilliance ESBL Agar (Oxoid) in order to quantify MRE levels and to study the impact of clinical factors and commensal bacteria on MRE intestinal colonization levels (results from this study will be published elsewhere). A subset of 56 patients was included in the present study. These patients had at some point been colonized by an MRE strain, and one or more consecutive samples were collected after the first MRE detection. Two hundred and seventeen samples representing the first positive MRE sample collected during a hospital admission period, plus all consecutive samples collected from this patient during this particular admission period until the MRE is no longer detected, were included in this study. Samples from 2 additional patients matching the criteria described above could not be included in this study since all the fecal material was used in a parallel study involving microbiome analysis.

Microbiological tests and molecular characterization

The 217 fecal samples were screened in September 2017 for the presence of the plasmid mediated colistin resistance mcr genes (including mcr-1, -2, -3, -4 and -5) by multiplex PCR [START_REF] Rebelo | Multiplex PCR for detection of plasmid-mediated colistin resistance determinants, mcr-1, mcr-2, mcr-3, mcr-4 and mcr-5 for surveillance purposes[END_REF]. PCR positive samples (N=5) were cultured to isolate the colistin-resistant strains harboring mcr genes by culture on LBJMR agar (containing 4 µg/ml colistin and 50 µg/ml vancomycin) and MacConkey agar supplemented with colistin (2 µg/ml) [START_REF] Bardet | LBJMR medium: a new polyvalent culture medium for isolating and selecting vancomycin and colistin-resistant bacteria[END_REF]. The isolated colonies were identified by matrix-assisted laser desorption ionization-time of flight mass spectrometry (MALDI-TOF MS) (Microflex, Bruker Daltonics, Bremen, Germany).

Subsequently, the antibiotics susceptibility of the isolates was determined by evaluating the minimum inhibitory concentration (MIC) using the broth microdilution method according to the Clinical and Laboratory Standard Institute (CLSI) guidelines for the colistin antibiotic, and using Etest method on Mueller Hinton agar for the other antibiotic families tested. In addition, the ESBL (blaTEM, blaSHV and blaCTX-M) encoding genes were screened in the colistin-resistant isolates.

Molecular epidemiology

The epidemiological relationship between the colistin resistant strains was determined by multilocus sequence typing (MLST). The seven housekeeping genes (adk, fumC, gyrB, icd, mdh, purA, recA) were amplified, sequenced and then blasted on the MLST database available on the Warwick web site (http://enterobase.warwick.ac.uk/species/ecoli/allele_st_search).

Conjugation experiments and plasmid analysis

Conjugative experiments were conducted using azide-resistant Escherichia coli J53 as a recipient, as described [START_REF] Bachiri | First report of the plasmid-mediated colistin resistance gene mcr-1 in Escherichia coli ST405 isolated from wildlife in Bejaia, Algeria[END_REF]. The transconjugants were selected on Luria Bertani (LB) agar (Beckton Dickinson, Le Pont de Claix, France) supplemented with sodium azide (120 µg/ml) and colistin (2 µg/ml). Transconjugant strains were screened for the presence of the colistin resistance genes (mcr genes) by PCR, and were subjected to antibiotic susceptibility testing as described above. Plasmid typing experiments were conducted on transconjugant strains using standard PCR [START_REF] Carattoli | Identification of plasmids by PCR-based replicon typing[END_REF].

Results

Microbiological, molecular and epidemiological characterization

Five of the 217 fecal samples tested were positive for the presence of the mcr-1 gene variant.

The PCR product sequence showed 100% identity to the published sequence [START_REF] Olaitan | Emergence of polymyxin resistance in Gram-negative bacteria[END_REF]. These samples corresponded to three patients with acute myeloid leukemia (AML) named Patient-1 (two positive samples), Patient-2 (one positive sample) and Patient-3 (two positive samples) aged 59, 48 and 63 years respectively (Table 1). All patients received an antimicrobial therapy prior sampling, including ciprofloxacin, piperacillin, meropenem, vancomycin and teicoplanin, but none had received colistin (Table 1). The patients also been subjected to coli-913). The E. coli strains isolated were resistant to at least five antibiotics among the sixteen tested, including colistin with MIC= 4 µg/ml (Table 1). The four colistin resistant E. coli carried the mcr-1 gene also harbored ESBL genes, including blaTEM-206 and blaTEM-98 (Table 1).

According to the MLST analysis, three different sequence types (STs) were assigned to the four E. coli isolates, including ST1196, ST140 and ST10. E. coli-643 and E. coli-648 strains retrieved from Patient-1 belonged to the same sequence type, ST1196.

Conjugation experiments and plasmid analysis

Conjugation experiment was conducted on the three E. coli harboring mcr-1, including E. coli-643, E. coli-866 and E. coli-913. The E. coli-648, considered a duplicate of Patient-1 E. coli strain, was not included in this experiment. Conjugative experiment allowed us to isolate 3 transconjugants (E. coli J53-643 Azide r , E. coli J53-866 Azide r and E. coli J53-913 Azide r ) resistant to colistin with MIC= 4 µg/ml. The antibiotic susceptibility profile of these transconjugant strains is presented in Table 2. All transconjugant strains were positive to the mcr-1 gene and the plasmid typing showed that this gene was located on IncP plasmid.

Discussion

Antibiotic resistance in immunocompromised patients was mainly represented by the Grampositive cocci group, such as methicillin-resistant Staphylococcus aureus (MRSA) and vancomycin-resistant Enterococcus (VRE). In the last decades, the emergence of infections due to Gram-negative bacteria, especially Enterobacteriaceae and Pseudomonas, has been noted [START_REF] Blennow | The challenge of antibiotic resistance in haematology patients[END_REF]. The Gram-negative bacilli identified as the predominant infectious agent in hematological patients is E. coli. In this group of patients, E. coli most often exhibits a high rate of resistance to quinolones, ceftazidime and beta-lactam antibiotics [START_REF] Blennow | The challenge of antibiotic resistance in haematology patients[END_REF]. These results have been observed in our study where the E. coli strains, isolated from our patients with acute leukemia, in addition to other antibiotics, also exhibited resistance to beta-lactam and quinolones.

To counter the emergence of MDR bacteria in hematological patients, the use of colistin as monotherapy or in combination with other antibiotics, such as beta-lactams, aminoglycosides, tigecycline or fosfomycin, has been suggested by the current American and European guidelines on febrile neutropenia [1]. Unfortunately, since its use in clinical settings, colistin resistance has increased and this is becoming very alarming, especially since the detection of a plasmid containing the mcr-1 colistin resistance gene [START_REF] Liu | Emergence of plasmidmediated colistin resistance mechanism MCR-1 in animals and human beings in China: a microbiological and molecular biological study[END_REF]. The result that emerged in our study and which worried us was the fact that we isolated from the leukemia patients E. coli strains resistant not only to the antibiotics mentioned above but also to the antibiotic of last resort, colistin, due to the presence of the mcr-1 gene.

The colistin resistance mcr-1 gene has been reported worldwide, mainly in animals [START_REF] Baron | Molecular mechanisms of polymyxin resistance: knowns and unknowns[END_REF]. In Spain, the mcr-1 gene has also been identified in animals (poultry, pigs and swine) [START_REF] Quesada | Detection of plasmid mediated colistin resistance (MCR-1) in Escherichia coli and Salmonella enterica isolated from poultry and swine in Spain[END_REF][START_REF] Poirel | MCR-2mediated plasmid-borne polymyxin resistance most likely originates from Moraxella pluranimalium[END_REF], in the environment (wastewater and sewage water) [START_REF] Lekunberri | Detection and quantification of the plasmidmediated mcr-1 gene conferring colistin resistance in wastewater[END_REF][START_REF] Ovejero | Spread of mcr-1-carrying Enterobacteriaceae in sewage water from Spain[END_REF] and a few studies reported the mcr-1 gene in clinical isolates (urine, blood, sputum) [START_REF] Baron | Molecular mechanisms of polymyxin resistance: knowns and unknowns[END_REF]. In patients with leukemia, the colistinresistance mcr-1 variant has been reported in five studies over the world including China, Austria and Italy. This gene was mostly detected in E. coli strains followed by Klebsiella pneumoniae and conferred to these strains a resistance to colistin with MIC ranging from 4 to 8 µg/ml [START_REF] Zhang | Evolution of tigecyclineand colistin-resistant CRKP (carbapenem-resistant Klebsiella pneumoniae) in vivo and its persistence in the GI tract[END_REF][START_REF] Tietgen | Impact of the colistin resistance gene mcr-1 on bacterial fitness[END_REF][START_REF] Tian | MCR-1-producing Klebsiella pneumoniae outbreak in China[END_REF][START_REF] Hartl | Detection of the mcr-1 Gene in a Multidrug-Resistant Escherichia coli Isolate from an Austrian Patient[END_REF][START_REF] Pilato | mcr-1.2, a New mcr Variant Carried on a Transferable Plasmid from a Colistin-Resistant KPC Carbapenemase-Producing Klebsiella pneumoniae Strain of Sequence Type 512[END_REF]. This was reported in our study, where the colistin resistance mcr-1 gene was detected in E. coli strains resistant to colistin with MIC= 4 µg/ml in high-risk patients with a weak immune system, who are leukemia patients, in Spain.

Concerning the genetic support of colistin resistance, the mcr-1 gene was generally identified in different plasmid types such as IncI2, IncHI2, IncX4 and IncP [START_REF] Hadjadj | Study of mcr-1 gene-mediated colistin resistance in Enterobacteriaceae isolated from humans and animals in different countries[END_REF]. In our strains, the mcr-1 gene was located on IncP transferable plasmid type. In Spain, the IncP plasmid type has been previously reported to be associated with different resistance genes, such as carbapenem resistant genes [START_REF] Yao | Insights into a novel blaKPC-2-encoding IncP-6 plasmid reveal carbapenem-resistance circulation in several Enterobacteriaceae species from wastewater and a hospital source in Spain[END_REF], but has never been associated with the mcr-1 gene.

Worldwide, different sequence types of E. coli harboring the mcr-1 gene have been detected in leukemia patients such as ST10, ST46, ST58, ST156, ST607 and ST3944 [START_REF] Tian | MCR-1-producing Klebsiella pneumoniae outbreak in China[END_REF][START_REF] Hartl | Detection of the mcr-1 Gene in a Multidrug-Resistant Escherichia coli Isolate from an Austrian Patient[END_REF][START_REF] Pilato | mcr-1.2, a New mcr Variant Carried on a Transferable Plasmid from a Colistin-Resistant KPC Carbapenemase-Producing Klebsiella pneumoniae Strain of Sequence Type 512[END_REF]. ST10 has made a significant contribution to the dissemination of the mcr-1 gene worldwide [START_REF] Hadjadj | Study of mcr-1 gene-mediated colistin resistance in Enterobacteriaceae isolated from humans and animals in different countries[END_REF]. In our study, the ST10 was reported in only one strain. ST1196 was the predominant sequence type detected in our strains (2 out of 4). This ST was reported in only one study in which the mcr-1 positive E. coli strain was recovered from a wastewater treatment plant at West China Hospital [START_REF] Zhao | Remarkable diversity of Escherichia coli carrying mcr-1 from hospital sewage with the identification of two new mcr-1 variants[END_REF]. In addition, our study reports for the first time the association of E. coli ST140 with the presence of the mcr-1 gene. Unfortunately, the lack of information on the history of patients; if they have been hospitalized, if they have travelled to a high-risk country where the mcr-1 gene is endemic or if the patient has been in contact with animals, does not allow us to determine the epidemiology and origin of the mcr-1 gene detected in our study.

In the present study, we detected the plasmid-mediated colistin resistance mcr-1 gene in patients not treated with colistin who were hospitalized between 2013 and 2015 at La Fe University Hospital. This suggests that the appearance of the mcr-1 gene in these patients was not due to selection pressure with this antibiotic. Despite the fact that our mcr-1 positives E. coli strains do not exhibit resistance to all the antibiotics tested, the detection of the mcr-1 gene on transferable plasmid is alarming. Indeed, the worrying scenarios emerging from this study are the acquisition of the plasmid encoding mcr-1 genes by MDR bacteria in these leukemia patients whose immune systems are already weakened, and the transmission of these resistant bacteria from one patient to another. This situation can lead the clinician into therapeutic impasse. This is exactly what was reported in a study conducted by Di Pilato et al.

where the mcr-1 gene was identified in an MDR KPC-producing K. pneumoniae ST512 [START_REF] Pilato | mcr-1.2, a New mcr Variant Carried on a Transferable Plasmid from a Colistin-Resistant KPC Carbapenemase-Producing Klebsiella pneumoniae Strain of Sequence Type 512[END_REF].

In these cases, what measures should be taken to prevent the spread of colistin resistance in this high-risk population? In our situation, the isolation of patients with mcr genes encoding colistin resistance could be an urgent solution to avoid any risk of transmission of such genes to other immunocompromised patients. In the long-term, the implementation of rapid tests to detect antibiotic resistance genes and the systematic screening of the carrying of antibiotic resistance genes as soon as they arrive at the hospital, would be a great solution for controlling the spread of antibiotic resistance genes and thus avoid any therapeutic impasse.

Patients

  either chemotherapy, transplant or both. Patient-1 received two chemotherapy and one bone marrow transplant prior sampling the two mcr-1 positive samples, unlike the positive mcr-1 sample from Patient-2 was collected while the patient showed clinical signs of infection and after undergoing bone marrow transplantation. During chemotherapy treatment of the Patient-3, the two mcr-1 positive samples were collected. The culture method allowed us to isolate four E. coli strains among the five mcr-1 positive samples; two strains from Patient-1 (E. coli-643 and E. coli-648), one strain from Patient-2 (E. coli-866) and one strain from Patient-3 (E.

Acknowledgements

The authors thank Linda Hadjadj for her technical assistance.

The authors also thank CookieTrad for proofreading the text.

Funding source

This work was supported by the ANR FloraStopInfectMRE project, InfecERA-ERANET-Acciones complementarias (PCIN-2015-094) to CU and by IHU Méditerranée Infection, Marseille, France and by the French Government under the «Investissements d'avenir» (Investments for the Future) proGram managed by the Agence Nationale de la Recherche (ANR, fr: National Agency for Research), (reference: Méditerranée Infection 10-IAHU-03). This work was supported by Région Provence Alpes Côte d'Azur and European funding FEDER PRIMI.

Conflict of interests

The authors declare they have no conflicts of interest