Remi R Pilon
email: picon@clermont.inra.fr

Catherine Picon-Cochard

Juliette M G Bloor

Sandrine S Revaillot

Estelle Kuhn Kuhn

Robert R Falcimagne

Philippe Balandier

J.-F Soussana

Angela R Hodge

Pilon

: R Falcimagne

Grassland root demography responses to multiple climate change drivers depend on root morphology

Keywords: Root dynamics, Warming, Summer precipitation reduction, Elevated CO 2, Minirhizotron, Root lifespan, Cox proportional hazards

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

Introduction

Increased global carbon dioxide concentration, air temperature, and in many regions reductions in precipitation during the growing season, are some of the most important aspects of global change predicted to affect ecosystems in coming years [START_REF] Ipcc | Climate Change 2007: the physical science basis[END_REF]. It is well established that warming stimulates plant growth and transpiration and affects soil environment through increased soil temperature and reduced soil water content (SWC) [START_REF] Luo | Terrestrial carbon cycle feedback to climate warming[END_REF]. A reduction in growing season precipitation should further reduce SWC, while elevated CO 2 may increase SWC through increased plant water-use efficiency [START_REF] Niklaus | Soil moisture dynamics of calcareous grassland under elevated CO 2[END_REF][START_REF] Wullschleger | Plant water relations at elevated CO 2 -implications for water-limited environments[END_REF]. In temperate regions, air warming and elevated CO 2 are expected to both stimulate plant growth, while variability in precipitation may counterbalance aboveground growth responses via changes in SWC and in belowground processes [START_REF] Norby | Evaluating ecosystem responses to rising atmospheric CO 2 and global warming in a multifactor world[END_REF][START_REF] Pendall | Below-ground process responses to elevated CO 2 and temperature: a discussion of observations, measurement methods, and models[END_REF]). Thus ecosystem responses to multiple climate drivers are complex, which indicates possible non-additive effects of simultaneouslyoccurring factors [START_REF] Shaw | Grassland responses to global environmental changes suppressed by elevated CO 2[END_REF][START_REF] Luo | Modeled interactive effects of precipitation, temperature, and [CO 2 ] on ecosystem carbon and water dynamics in different climatic zones[END_REF][START_REF] Zhou | Modeling patterns of nonlinearity in ecosystem responses to temperature, CO2 and precipitation changes[END_REF].

It is expected that the roots and the below-ground compartment of grassland will be particularly affected by climate change as, in this ecosystem, most of the C is stored below-ground [START_REF] Steinaker | Belowground litter contributions to nitrogen cycling at a northern grassland-forest boudary[END_REF]. However previous single-factor climate manipulation studies showed inconsistent patterns of root system response to warming and elevated CO 2 . For example, warming has been shown to stimulate root production and root mortality [START_REF] Fitter | Root production, turnover and respiration under two grassland types along an altitudinal gradient: influence of temperature and solar radiation[END_REF][START_REF] Fitter | Root production and turnover in an upland grassland subjected to artificial soil warming respond to radiation flux and nutrients, not temperature[END_REF][START_REF] Bai | Increased temperature and precipitation interact to affect root production, mortality, and turnover in a temperate steppe: implications for ecosystem C cycling[END_REF]) and reduce root biomass [START_REF] Boeck | Combined effects of climate warming and plant diversity loss on above-and below-ground grassland productivity[END_REF]). Under elevated CO 2 , root production may increase [START_REF] Fitter | Root production and mortality under elevated atmospheric carbon dioxide[END_REF][START_REF] Fitter | Root production and turnover and carbon budgets of two contrasting grasslands under ambient and elevated atmospheric carbon dioxide concentrations[END_REF][START_REF] Soussana | Long-term effects of CO 2 enrichment and temperature increase on a temperate grass sward. 2. Plant nitrogen budgets and root fraction[END_REF][START_REF] Higgins | Root production and demography in a California annual grassland under elevated atmospheric carbon dioxide[END_REF][START_REF] Sindhøj | Root biomass dynamics in a semi-natural grassland exposed to elevated atmospheric CO 2 for 5 years[END_REF]Milchunas et al. 2005a, b;[START_REF] Lecain | Root biomass of individual species, and root size characteristics after 5 years of CO 2 enrichment on native shortgrass steppe[END_REF][START_REF] Anderson | Root responses along a subambient to elevated CO 2 gradient in a C3-C4 grassland[END_REF], decrease [START_REF] Pendall | Soil carbon storage under simulated climate change is mediated by plant functional type[END_REF] or show no response [START_REF] Arnone | Dynamics of root systems in native grasslands: effects of elevated atmospheric CO 2[END_REF][START_REF] Dukes | Responses of grassland production to single and multiple global environmental changes[END_REF]. For root mortality, elevated CO 2 is also associated with an increase [START_REF] Fitter | Root production and mortality under elevated atmospheric carbon dioxide[END_REF][START_REF] Fitter | Root production and turnover and carbon budgets of two contrasting grasslands under ambient and elevated atmospheric carbon dioxide concentrations[END_REF]Milchunas et al. 2005a) or no change [START_REF] Arnone | Dynamics of root systems in native grasslands: effects of elevated atmospheric CO 2[END_REF][START_REF] Higgins | Root production and demography in a California annual grassland under elevated atmospheric carbon dioxide[END_REF]. Unlike root growth and mortality, root lifespan has faced little attention in climate change studies [START_REF] Zak | Elevated atmospheric CO 2 , fine roots and the response of soil microorganisms: a review and hypothesis[END_REF]: Milchunas et al. (2005a) observed higher lifespan under elevated CO 2 , but [START_REF] Anderson | Root responses along a subambient to elevated CO 2 gradient in a C3-C4 grassland[END_REF] reported opposite results. Responses in root lifespan could be mediated by root diameter responses because of the tight link observed between these two variables [START_REF] Eissenstat | The ecology of root lifespan[END_REF][START_REF] Eissenstat | Building roots in a changing environment: Implications for root longevity[END_REF][START_REF] Gill | Longevity and turnover of roots in the shortgrass steppe: Influence of diameter and depth[END_REF]. To date, relatively few studies have focused on the combined effects of warming, precipitation reduction and elevated CO 2 on the below-ground compartment (e.g. [START_REF] Pendall | Below-ground process responses to elevated CO 2 and temperature: a discussion of observations, measurement methods, and models[END_REF][START_REF] Wan | Responses of soil respiration to elevated CO 2 , air warming, and changing soil water availability in a model old-field grassland[END_REF]. Available data indicate no root mass response to combined warming and elevated CO 2 [START_REF] Soussana | Long-term effects of CO 2 enrichment and temperature increase on a temperate grass sward. 2. Plant nitrogen budgets and root fraction[END_REF][START_REF] Kandeler | The response of soil microorganisms and roots to elevated CO 2 and temperature in a terrestrial model ecosystem[END_REF][START_REF] Shaw | Grassland responses to global environmental changes suppressed by elevated CO 2[END_REF], while amounts of precipitation seem to modulate the positive effect of elevated CO 2 (Milchunas et al. 2005b) and warming [START_REF] Bai | Increased temperature and precipitation interact to affect root production, mortality, and turnover in a temperate steppe: implications for ecosystem C cycling[END_REF]) on root growth. Overall, these findings reflect the complexity of interactions, occurring between direct and indirect effects of climate drivers on above and belowground processes. Furthermore, methodology can substantially vary between studies, concerning root dynamics measurement [START_REF] Arnone | Dynamics of root systems in native grasslands: effects of elevated atmospheric CO 2[END_REF][START_REF] Luo | Terrestrial carbon cycle feedback to climate warming[END_REF][START_REF] Milchunas | Estimating root production: comparison of 11 methods in shortgrass steppe and review of biases[END_REF] and the intensity of the climate factors imposed in experiments, e.g. range of air warming between 1 to 4 °C, and a doubling of atmospheric CO 2 concentration.

Root responses to climate change may be complicated by interactions occurring between plant physiological and morphological processes. In perennial grassland, standing root mass is the result of growth, mortality and decomposition processes that occur simultaneously and vary through season [START_REF] Matthew | Relationship between root production and tiller appearance rates in perennial ryegrass (Lolium perenne L.)[END_REF][START_REF] Steinaker | Phenology of fine roots and leaves in forest and grassland[END_REF][START_REF] Garcia-Pausas | Seasonal patterns of belowground biomass and productivity in mountain grasslands in the Pyrenees[END_REF]. Moreover, root dynamics may be influenced by variability of root morphology within a root system, i.e. by the presence of different root classes. Although several authors have outlined the importance of separating fine and coarse roots [START_REF] Fitter | An architectural approach to the comparative ecology of plant-root systems[END_REF][START_REF] Jackson | Integrating resource heterogeneity and plant plasticity: modelling nitrate and phosphate uptake in a patchy soil environment[END_REF][START_REF] Craine | Functional traits, productivity and effects on nitrogen cycling of 33 grassland species[END_REF][START_REF] Waisel | Functional diversity of various constituents of a single root system[END_REF][START_REF] Bahn | Root respiration in temperate mountain grasslands differing in land use[END_REF][START_REF] Pierret | Root functional architecture: a framework for modeling the interplay between roots and soil[END_REF][START_REF] Picon-Cochard | Effects of species, root branching order and season on the root traits of 13 perennial grass species[END_REF]), this has not been considered in recent climate change studies. Variation in root order is likely to affect root lifespan and main root functions including uptake and respiration in response to climate change [START_REF] Eissenstat | Building roots in a changing environment: Implications for root longevity[END_REF][START_REF] Norby | Root dynamics and global change: seeking an ecosystem perspective[END_REF][START_REF] Wan | Responses of soil respiration to elevated CO 2 , air warming, and changing soil water availability in a model old-field grassland[END_REF]. Root diameter may be a useful proxy of root order, and a more easily measurable variable to determine how root turnover responds to climate change.

This paper reports root length growth and mortality (based on root disappearance, including decomposition) assessed with minirhizotron data obtained for the upper soil horizon of an upland perennial grassland as part of the Clermont Climate Change Experiment [START_REF] Bloor | Effects of warming, summer drought, and CO 2 enrichment on aboveground biomass production, flowering phenology, and community structure in an upland grassland ecosystem[END_REF]. In this experiment, grassland has been exposed to four climate treatments in a replicated additive design since 2005: control (C), elevated temperature (T), combined elevated temperature and summer precipitation reduction (TD), combined elevated temperature, summer precipitation reduction and elevated atmospheric CO 2 (TDCO 2 ). The TDCO 2 treatment simulates projected climate conditions for the end of the century (ACACIA A2 scenario for 2080 in central France, [START_REF] Ipcc | Climate Change 2007: the physical science basis[END_REF]. Previous results from this study system showed a short-term stimulation of above-ground production under air warming, associated with an advance in flowering time for the dominant grass species [START_REF] Bloor | Effects of warming, summer drought, and CO 2 enrichment on aboveground biomass production, flowering phenology, and community structure in an upland grassland ecosystem[END_REF]. This effect disappeared after 2 years, and a progressive decrease in aboveground production was observed with decline in grass abundance at the expense of legumes [START_REF] Cantarel | Four years of simulated climate change reduces aboveground productivity and alters functional diversity in grassland[END_REF]. Overall, the results obtained for the aboveground part suggested a main effect of air warming without marked impact of other climate drivers (precipitation reduction, elevated CO 2 ). In this context, we studied how the root system responded to the climate treatments. We particularly evaluated how single or combined warming, summer precipitation reduction and elevated CO 2 affect root dynamics by analysing abiotic (soil water content, air temperature, elevated CO 2 ) and biotic (diameter classes) factors. In addition, complementary methods were used to quantify effect of treatments on median root longevity (Kaplan-Meier) and on risk of mortality induced by co-variates (Cox). We hypothesized that: i) root growth is stimulated by air warming like above-ground biomass, ii) long-term warming effects on soil water content reduction may decrease root longevity, fine roots being more sensitive than coarse roots, iii) seasonal dynamics of the grassland root system masks root responses to climate change.

Materials and methods

Study site and experimental design

A full description of the Clermont Climate Change Experiment design and study site is given in [START_REF] Bloor | Effects of warming, summer drought, and CO 2 enrichment on aboveground biomass production, flowering phenology, and community structure in an upland grassland ecosystem[END_REF]. Briefly, the study system is a native upland grassland, extensively managed, located in the French Massif Central area (control site C, Theix, France, 45°43'N, 03°01'E, 850 ma.s.l.), subjected for 15 years to low-intensity sheep grazing without inorganic fertilisation [START_REF] Louault | Plant traits and functional types in response to reduced disturbance in a semi-natural grassland[END_REF]. The grassland had relatively low species diversity (11 species) and was dominated by perennial grasses (more than 90 % of phytomass) such as Agrostis capillaris, Elytrigia repens and Festuca arundinacea, but also included a small number of dicotyledon species (Stellaria graminea, Taraxacum officinale, Vicia sativa). The soil is a Cambisol developed on granite (C:N, 9.8; pH water , 6.2) and the climate is semi-continental with oceanic influences (mean annual temperature, 8.7 °C, mean annual precipitation, 780 mm on average of [1986][1987][1988][1989][1990][1991][1992][1993][1994][1995][1996][1997][1998][1999][2000][2001][2002].

In November 2004, 80 monoliths (0.5 m x 0.5 m, 0.4 m deep) were extracted from homogeneous, multispecies zones to avoid monoliths with a dominance of any one species. Botanical surveys at the start of the experiment indicated that the 0.5 x 0.5 m monoliths contained 8±0.2 species with the same composition on average as the field plot.

Monoliths were encased in stainless-steel boxes with drainage holes in the bottom. After a six-month on-site acclimatization period, the monoliths were randomised and assigned to experimental treatments: 20 monoliths were left at the control site (850 m) while the remainder was transferred to a lowland area 15 km away from the control site (lowland site, Clermont-Ferrand, France, 45°46'N, 03°08'E, 350 ma.s.l.; mean annual temperature, 12.4 °C; mean annual precipitation, 575 mm). Four monoliths were grouped together to form a single experimental unit (1 m x 1 m). Five experimental units (i.e. 20 monoliths) and 15 experimental units (i.e. 60 monoliths) were allocated to the control and lowland sites, respectively. Units were placed in holes in the ground lined with a 0.5 m-deep layer of volcanic gravel and a 0.4 m layer of subsoil (0.4 -0.8 m) collected from the original monolith extraction site. The experimental units were placed 7 m apart to ensure that there was no risk of atmospheric contamination between CO 2 enriched and ambient CO 2 units.

Four climate treatments were applied in a replicated additive design with five replicate experimental units per treatment: control (C, control site), elevated air temperature (T, lowland site, called warmed site thereafter), elevated temperature and summer precipitation reduction (TD, warmed site), elevated temperature and summer precipitation reduction plus elevated atmospheric CO 2 (TDCO 2 , warmed site). Climate manipulations were based on ACACIA A2 climate change scenario projections for the area to the 2080 time-horizon (IPCC 2007) and corresponded to a 3.5 °C temperature increase, a 20 % reduction in summer precipitation, and atmospheric CO 2 level of 600 μmol mol -1 . In our study, a mean annual warming of 3.5 °C (T) was achieved by transplanting monoliths along an altitudinal gradient (500 m). To ensure the same water supply (i.e. precipitation) for the C and T treatments, rainfall was monitored at both sites and treatment T was irrigated weekly when necessary. For TD and TDCO 2 treatments, reduction of summer (June, July and August) precipitation was achieved by intercepting rainfall with plastic roofs at 1.2 m above experimental units during rainfall events which lasted for 24-48 h (for details see [START_REF] Bloor | Effects of warming, summer drought, and CO 2 enrichment on aboveground biomass production, flowering phenology, and community structure in an upland grassland ecosystem[END_REF]. The average precipitation reduction between T and TD and TDCO 2 , respectively, was 57, 78 and 56 mm in 2006, 2007 and 2008, corresponding to a 20 % rainfall reduction during summer period. Elevated CO 2 was applied using a mini-FACE (Free Air CO 2 Enrichment) system following [START_REF] Miglietta | Free-air CO 2 enrichment (FACE) of a poplar plantation: the POPFACE fumigation system[END_REF].

Once all the experimental units were in place (on April 27, 2005), grass monoliths were cut to a height of 5 cm twice a year, in spring and autumn (end April and end October), to mimic the previous extensive management. No fertiliser was applied during the experimental period in line with previous management.

Measurements

Climate data (air temperature, °C; precipitation, mm; atmospheric CO 2 concentration, μmol mol -1 ; net global radiation, W m -2 ) were recorded every 30 min using two automatic dataloggers (Campbell Scientific, UK) installed at each site. Mean air temperature difference between the two sites was on average 3.3±0.1 °C over the three-year study period. Differences of monthly precipitation and cumulated potential Penman evapotranspiration (PET, mm month -1 ) were used as an index of water balance at both sites (P -PET). For the three year-study, cumulative P -PET was positive (656 mm) at the control site (C) and negative (-566 mm) at the warmed site. Volumetric soil water content (SWC, cm 3 cm -3 ) was measured every 30 min, in each experimental unit, with a vertically inserted 20 cm long Echo dielectric capacity probe connected to a datalogger (Ech 2 O, Decagon Devices, USA). Soil water potential was determined from gravimetric soil water content (SWCg) measurements at field capacity (0.03 MPa) and at wilting point (1.5 MPa) with a dewpoint meter (PotentiaMeter WP4, Decagon Devices, USA). The relationship between SWCg obtained by weighing 5-L pots and SWC given by the probes inserted in the pots made it possible to convert SWC to SWCg, as: SWCg0 0.0948 Ln (SWC)+0.4021, n0108, R 2 00.81.

Root data were obtained using a minirhizotron system (BTC-2, Bartz Technology, CA, USA). In April 2005, 50 cm long transparent acrylic tubes were inserted at an angle of 45°with one tube per monolith. However, only two tubes per experimental unit were analysed due to logistic constraints. Root images (1.35 cm x 1.8 cm) were recorded from March 14, 2006 to January 15, 2009; images from the same localizations were recorded once a month during the growing period in 2006 and 2007, and twice a month from 2008 onwards. Three successive images were analysed from the main root zone (5-10 cm depth) of each tube as described in [START_REF] Fitter | Root production and mortality under elevated atmospheric carbon dioxide[END_REF].

Length, diameter and status (live, dead) of root segments were manually recorded using WinRHIZO-TronMF software (V2005a, Regent Instruments, CA). For each date, length was modified according to growth event. Root diameter was determined after the first image and not modified thereafter; once grasses form roots, their diameter remains stable because there is no secondary growth in monocots (the majority of species in the community). Root length was determined per unit tube area (mm cm -2 ). Three diameter classes were defined: ≤ 0.1 mm, 0.1<D≤ 0.2 mm and>0.2 mm, representing 22, 55 and 23 % of total root length across treatments. Root mortality was based on root disappearance following [START_REF] Johnson | Advancing fine root research with minirhizotrons[END_REF], and included time of decomposition. Root length growth and mortality rates were calculated as the difference in length of live and dead roots, respectively, between two successive dates divided by the number of days (mm cm -2 day -1 ). Mean values by season are shown for each year. In addition, annual increment in new root length and root length mortality (mm cm -2 ) was calculated for 2006, 2007 and 2008 and for the three-year study period (2006)(2007)(2008). Standing root length (mm cm -2 ) was calculated as the difference between annual increment of new root length and root length mortality per year and for the three-year study period. All root variables were calculated as the average values of the two tubes present in the same experimental unit. A total of 10,360 root segments were drawn and an average of 2,590 root segments per treatment was used in the analysis.

Statistical analyses

Root length dynamics (growth, mortality rates, annual and total increment of new root length, mortality and standing root length) were analyzed using mixedmodel repeated-measures ANOVA. All interactions between treatment, diameter classes and season or year were included in the model. Treatments were compared using orthogonal contrasts and Tukey's honest significant difference (HSD) procedure. The effect of climate change drivers was analyzed as follows: warming effect, C vs T; summer precipitation reduction effect under warming, T vs TD; CO 2 effect under warming and summer precipitation reduction, TD vs TDCO 2 ; simultaneous application of warming, summer precipitation reduction and elevated CO 2 (2080 climate scenario), C vs TDCO 2 . When necessary, variables were square-root or log-transformed to normalize the data prior to ANOVA.

For estimation of root survival times, dates of root birth and root death were assumed to be the mid-point between recording dates before and after appearance and disappearance, respectively. Roots living beyond the last recording time were considered to be right censored, i.e. root survival time is known only to be at least as long as the time to the last recording. The proportion of roots surviving as a function of age was described by the Kaplan-Meier survival function, a non-parametric statistical technique which allows estimation of median longevity [START_REF] Allison | Survival analysis using SAS. A practical guide[END_REF]. Median root lifespan was estimated for 50 % of surviving roots. The comparison of the survival curves of groups is based on a formal non-parametric statistical test called the logrank test. All these statistical analyses were carried out using Statgraphics Plus 4.1 (Statistical Graphics Corp., Maryland, USA).

A Cox proportional hazards regression [START_REF] Cox | Regression models and life tables[END_REF][START_REF] Allison | Survival analysis using SAS. A practical guide[END_REF][START_REF] Wells | Marked differences in survivorship among apple roots of different diameters[END_REF]) was used to examine the effects of biotic and abiotic covariates on risk of root mortality. Cox analysis on the full dataset included the effects of all treatments (C, T, TD, TDCO 2 ), year of appearance (2006,2007,2008), season of appearance (growth season: spring+summer; dormant season: autumn+winter), root diameter classes, air temperature and SWC measured between the last date the roots were observed and date of disappearance. The model was fitted using a maximum likelihood method to estimate the regression parameters (parameter estimate) associated with each explanatory variable and calculating the standard error and a chi-square test to test the null hypothesis that each parameter estimate is equal to zero. The sign of the parameter estimate indicates whether the hazard of mortality decreases (negative sign) or increases (positive sign) as the covariate value increases. A conditional risk ratio (or hazard ratio) was also computed from the parameter estimates. This ratio has several interpretations depending on the nature of the covariate. For dichotomous covariates, such as season of appearance, which was assigned as growth season or dormant season, the risk ratio can be interpreted as the ratio of hazard in the dormant compared with the growth season, controlling for all other covariates. For quantitative covariates such as air temperature and SWC, subtracting one from the risk ratio and multiplying it by 100 gives the estimated percent change in the hazard of mortality associated with a one-unit change in the covariate, controlling for all other covariates [START_REF] Allison | Survival analysis using SAS. A practical guide[END_REF]. Cox analyses were run using the PROC PHREG procedures of the SAS software package (v8.1, SAS Institute Inc., Cary, NC, USA).

Results

Soil water content

During the course of the experiment, values of soil water content (SWC) ranged between 0.04 and 0.14 cm cm -3 in the warmed site (T, TD and TDCO 2 ) and between 0.08 and 0.15 cm cm -3 in the control site (C), leading to lower (-18 %, P<0.01) values in the warmed site (Fig. 1). However, due to large withintreatment variation, no between-treatment differences were found for T, TD and TDCO 2 during study period. Irrespective of treatment, SWC values were lower in the summer period (particularly in 2006, a drier year than 2007 and 2008) and higher in winter and most SWC values fell within the range of soil water potential measured between field capacity and wilting point (1.5 MPa).

Root length growth rate and cumulative new root length

Over the three-year study, seasonal root growth rate was higher (P < 0.01) in the warmed site (T, 0.31 mm cm -2 day -1 ) than in the control site (C, 0.23 mm cm -2 day -1 ) (Fig. 2, Table 1), but the effects were more pronounced in spring 2007 and summer 2008. Cumulative new root length ranged from 490.9 (C treatment) to 650.3 mm cm -2 (TDCO 2 ) with significant higher (P<0.01) values in the warmed site (C vs T, 490.9 vs 593.5 mm cm -2 ) and in the combined climate change drivers (C vs TDCO 2 ) (Fig. 3, Table 1). On an annual basis, a significant warming effect (C vs T) on growth rate was observed for 2 years (2007 and 2008), whereas cumulative new root length was only higher in 2008 (Table 1). Similarly in combined climate change drivers (C vs TDCO 2 ), root growth was higher (P<0.05) in 2007 and new root length increased in 2007 and 2008. In contrast, summer precipitation reduction (T vs TD) and elevated CO 2 (TD vs TDCO 2 ) had no significant effect on growth rate and on cumulative new root length (Tables 1,2, Figs. 2,3).

Over the course of the experiment, growth rate and cumulative new root length showed strong seasonal and annual effects across all climate treatments (Fig. 2, Table 2): growth rates were higher in spring compared to summer and autumn-plus-winter (+43 % and +140 %, respectively); and in summer growth rate was higher (+68 %) than in autumn and winter (Fig. 2).

Treatment effects on cumulative new root length varied depending on diameter class (significant treatment x class interaction, P<0.001, Table 2, Fig. 3). For the fine roots (≤ 0.1 mm), new root length increased from 70.1 to 135.9 mm cm -2 in response to air warming and from 7.1 to 114.1 mm cm -2 in response to combined climate change drivers. In contrast elevated CO 2 (TD vs TDCO 2 ) showed an opposite effect (-28 %, 158 vs 114.1 mm cm -2 ). In the case of the coarse roots (> 0.2 mm), elevated CO 2 and combined climate change drivers (C vs TDCO 2 ) increased new root length (219.2 vs 145.7 and 219.2 vs 136.5 mm cm -2 , respectively), whereas Fig. 2 Effects of climate treatment on seasonal dynamics in root growth rate. Treatment codes are as follows: C 0 control; T 0 elevated temperature (+3.5 °C), TD 0 elevated temperature and summer precipitation reduction (-20 %); TDCO 2 0 elevated temperature, summer precipitation reduction and CO 2 enrichment (+200 μmol mol -1 ). Vertical bars correspond to standard error (n05). sp 0 spring; su 0 summer; aut 0 autumn; wint 0 winter warming and summer precipitation reduction had no effect (Fig. 3). The intermediate diameter class (0.1< D≤0.2 mm) always exhibited higher root growth than fine (x 2.5) and coarse (x 1.8) roots, fine roots having the lowest new root length (Fig. 3).

Root length mortality based on disappearance

Over the three year-study, climate treatment had no overall effect on cumulative root length mortality, but varied with diameter class (significant treatment x class interaction, P<0.001, Table 2, Fig. 3). For fine roots (≤ 0.1 mm), root length mortality increased (P<0.001) in response to warming (75 %) and ranged from 56.8 (treatment C) to 99.2 mm cm -2 (treatment T) and also increased in response to summer precipitation reduction (T vs TD, 99.2 to 137.7 mm cm -2 ). In contrast, elevated CO 2 had an opposite effect (74.4 vs 137.7 mm cm -2 ), and combined climate change driver showed no effect. For coarse roots (> 0.2 mm), elevated CO 2 (TD vs TDCO 2 , 62.4 vs 103 mm cm -2 ) and combined climate change drivers (C vs TDCO 2 , 62.8 vs 103 mm cm -2 ) increased root mortality, whereas warming and summer precipitation reduction had no effect (Fig. 3).

Standing root length

Standing root length ranged between 201 and 305 mm cm -2 in the C and TDCO 2 treatments, respectively. Significant warming (C vs T, 201 vs 283.2 mm cm -2 , P<0.05) and combined climate change drivers (C vs TDCO 2 , 201 vs 305.3 mm cm -2 , P<0.01) effects were observed across diameter classes, while summer precipitation reduction (T vs TD) and elevated CO 2 (TD vs TDCO 2 ) had no effect (Tables 1, 2; Fig. 3). Furthermore, across treatments and diameter classes, standing root length, was significantly higher in 2006 than in 2007 and 2008 (+178 and +60 %, respectively) (Table 2). Only in 2008, elevated CO 2 and combined climate change 1).

Median root longevity based on root disappearance

Treatments and diameter classes significantly affected the median root longevity (MRL) given by Kaplan-Meier analysis (Table 3). MRL of all root types reached about 404 days in C and TDCO 2 treatments, which showed significant higher values than T and TD treatments, 369 and 355 days, respectively. This trend was also observed for the fine roots, but with lower MRL: average values of 343 days for C and TDCO 2 , and 285 days for T and TD treatments. In the case of coarse roots (> 0.2 mm), the effects of air warming on MRL was also observed in T and TD, but the future Fig. 3 Cumulative new root length, root length mortality and standing root length (mm cm -2 ) measured over 3 years for three diameter classes (≤ 0.1 mm, 0.1<D≤ 0.2 mm, > 0.2 mm) and for all roots subjected to four treatments C, T, TD, TDCO 2 . Treatment codes are as follows: C 0 control; T 0 elevated temperature (+3.5 °C), TD 0 elevated temperature and summer precipitation reduction (-20 %); TDCO 2 0 elevated temperature, summer precipitation reduction and CO 2 enrichment (+200 μmol mol -1 ). Vertical bars correspond to standard error (n05). Different letters correspond to statistical differences (P≤0.05). NS: P>0.05 climate treatment showed the lowest MRL (421 days). Furthermore, median root longevity was positively related with root diameter classes for three treatments (C, T, TD), while for the TDCO 2 treatment MRL did not differ between root diameters (Table 3).

Risk of root mortality based on root disappearance

Across all root cohorts, treatment, year of appearance, diameter classes and SWC were significantly related to the risk of root mortality (Table 4). The sign of the parameter estimate indicates whether the hazard of mortality decreases (negative sign) or increases (positive sign) as the covariate value increases. Risk of mortality was significantly (P<0.0001) lower in dry year (2006) than in wet years (2007,2008) and for coarse roots compared with fine roots. In addition, effect of treatments and of SWC was marginally significant (P<0.05): the risk of mortality decreased in TDCO 2 in comparison with C treatment and in case of high SWC values. However, season of appearance and air temperature had no effect on risk of root mortality.

Discussion

Grassland root responses to air warming and combined climate change drivers A better understanding on ecosystem responses to multiple climate change drivers are needed to anticipate future ecosystem functioning and provision of ecosystem services. Given that roots represent the plant interface with the soil, it is crucial to gain deeper insight into root responses to climate change. Here we report root dynamics of a grassland ecosystem subjected to combined effects of air warming, reduction of summer precipitation and elevated atmospheric CO 2 over 3 years. Our results show that root growth rate, measured in the most active root zone (5-10 cm), was sensitive to warming (C vs T) and to simultaneous application of warming, precipitation reduction and elevated CO 2 (C vs TDCO 2 ), leading to higher net root length production and consequently higher standing root length. This is in agreement with other studies showing that root growth is stimulated by warming [START_REF] Fitter | Root production, turnover and respiration under two grassland types along an altitudinal gradient: influence of temperature and solar radiation[END_REF][START_REF] Fitter | Root production and turnover in an upland grassland subjected to artificial soil warming respond to radiation flux and nutrients, not temperature[END_REF][START_REF] Pregitzer | Responses of tree fine roots to temperature[END_REF][START_REF] Wan | Responses of soil respiration to elevated CO 2 , air warming, and changing soil water availability in a model old-field grassland[END_REF][START_REF] Bai | Increased temperature and precipitation interact to affect root production, mortality, and turnover in a temperate steppe: implications for ecosystem C cycling[END_REF]. Notably, root growth mirrored the increase of springtime aboveground production (see [START_REF] Bloor | Effects of warming, summer drought, and CO 2 enrichment on aboveground biomass production, flowering phenology, and community structure in an upland grassland ecosystem[END_REF]), indicating that above and below-ground production increase simultaneously under warming, confirming our first hypothesis.

In our experiment, air warming also extended the growing season to early spring and late autumn which explains the higher root growth in the warmed treatments in these periods. However 2 years after the onset of treatment the positive effect of air warming on the above-ground compartment was not sustained [START_REF] Cantarel | Four years of simulated climate change reduces aboveground productivity and alters functional diversity in grassland[END_REF]). This decline is possibly linked to acclimation of photosynthesis [START_REF] Chapin | Responses of Arctic tundra to experimental and observed changes in climate[END_REF][START_REF] Shaw | Grassland responses to global environmental changes suppressed by elevated CO 2[END_REF], induced by negative water balance index (P -PET) and lower SWC (Fig. 1). Thus maintenance of root growth stimulation induced by air warming may correspond to a shift of C allocation to roots at the expense of the above-ground compartment, in warmer and drier soil conditions. Contrary to expectations, precipitation reduction had no effect on root growth and mortality, presumably due to relatively important soil moisture values, which indicated that no marked edaphic drought occurred during the study period (Fig. 1). Reduced soil moisture and subsequent limitation of nutrients under combined elevated CO 2 and warming may explain the lack of response under elevated CO 2 [START_REF] Shaw | Grassland responses to global environmental changes suppressed by elevated CO 2[END_REF][START_REF] Luo | Progressive nitrogen limitation of ecosystem responses to rising atmospheric carbon dioxide[END_REF], but see Hovenden et al. 2008). Positive effects of elevated CO 2 may be reduced when combined with air warming and reduction of summer precipitation. Warming was, thus, the main climate factor shaping root system response of our upland perennial grassland.

Linking morphology, median root longevity and covariates associated with risk of root mortality based on disappearance Root longevity is one of the most important variables for understanding belowground C cycling in ecosystems, as it determines how quickly C allocated to roots is transferred to microbial and soil C pools [START_REF] Gill | Global patterns of root turnover for terrestrial ecosystems[END_REF][START_REF] Norby | Evaluating ecosystem responses to rising atmospheric CO 2 and global warming in a multifactor world[END_REF]). However, multiple effects of climate change on plants and their soil environment makes simple prediction of root longevity a difficult task [START_REF] Eissenstat | Building roots in a changing environment: Implications for root longevity[END_REF]. This is further compounded by the numerous determinants of root longevity or root mortality [START_REF] Eissenstat | The ecology of root lifespan[END_REF][START_REF] Gill | Longevity and turnover of roots in the shortgrass steppe: Influence of diameter and depth[END_REF].

Consequently, identifying determinants of root mortality and root longevity and predicting source/ sink ecosystem activity under future climate scenarios remains a complex challenge. One of the most cited determinants of root longevity is diameter, with coarse roots demonstrating higher longevity than fine roots [START_REF] Eissenstat | The ecology of root lifespan[END_REF][START_REF] Gill | Longevity and turnover of roots in the shortgrass steppe: Influence of diameter and depth[END_REF]. Our data confirmed these previous studies for three of the four treatments (C, T, TD) (Tables 3,4). This may be , the hazard ratio can be interpreted as the estimated hazard for roots with a value of 1 vs the estimated hazard for roots with a value of 0 [START_REF] Allison | Survival analysis using SAS. A practical guide[END_REF] because thicker roots function as conduits and initiate new laterals, and are thus probably preferentially retained by the plant [START_REF] Wells | Marked differences in survivorship among apple roots of different diameters[END_REF]. They are also more suberized with more lignified cell walls in cortex or stele than fine roots, increasing their defence against herbivores or pathogens. Overall this contributes to increase their longevity. However the decrease in longevity of coarse roots and increase of fine root longevity observed under elevated CO 2 were unexpected. This is a novel finding that raises questions on the contrasted functioning under elevated CO 2 of these extreme root classes that together represent about 45 % of total root length. The faster turnover of the coarse roots can be due to higher foraging activity in depleted resource patches and to replacing them with more active roots in new soil sites. This would explain the higher rhizosphere activity for coarse roots grown under elevated CO 2 , as suggested in earlier study [START_REF] Hungate | The fate of carbon in grasslands under carbon dioxide enrichment[END_REF]. Milchunas et al. (2005a) also suggested that the higher diameter of roots grown under elevated CO 2 may be related to their water status. If these roots have higher water content, i.e. low dry matter content, we would also expect lower root tissue density and consequently lower root lifespan [START_REF] Ryser | Intra-and interspecific variation in root length, root turnover and the underlying parameters[END_REF]. In contrast, combined effects of warming and elevated CO 2 may have induced higher carbohydrate and nitrogen contents in the tissues of fine roots, which in turn may increase their lifespan (see [START_REF] Farrar | The Effects of increased atmospheric carbon dioxide and temperature on carbon partitioning, source-sink relations and respiration[END_REF]. Air warming (C vs T treatment) significantly reduced root longevity, irrespective of diameter class, thus confirming previous results obtained in grasslands [START_REF] Forbes | Temperature-induced alteration to root longevity in Lolium perenne[END_REF][START_REF] Fitter | Root production, turnover and respiration under two grassland types along an altitudinal gradient: influence of temperature and solar radiation[END_REF][START_REF] Fitter | Root production and turnover in an upland grassland subjected to artificial soil warming respond to radiation flux and nutrients, not temperature[END_REF]. Nevertheless, fine roots were more prone to mortality in warmed treatments than other diameter class roots, which confirm our hypothesis on the lifespan sensitivity of fine roots. The sensitivity could be explained by indirect effects of reduced soil water content due to higher plant evapotranspiration demand and lower water balance index. Indeed, our study revealed effects of soil water content on the risk of root mortality, depending on year, diameter class, and treatment (Table 4). These results suggest the importance of soil moisture in shaping the risk of root mortality in warmer seasons and drier years for upland grasslands, with fine roots being particularly sensitive. Root sensitivity to both average climate (comparison of C with T, TD and TDCO 2 ) and climatic variability (between-year and betweenseason differences) showed a strong difference depending on root diameter class, indicating that not all the grassland root system is similarly affected by climate change components. This may have consequences on plant and soil relationships, especially root and microbe activities that are critical for the decomposition of soil organic matter and its mineralization in future climates [START_REF] Zak | Elevated atmospheric CO 2 , fine roots and the response of soil microorganisms: a review and hypothesis[END_REF].

Seasonal dynamics of perennial grassland root system

Our study confirms the importance of seasonality in grassland ecosystem root dynamics [START_REF] Matthew | Relationship between root production and tiller appearance rates in perennial ryegrass (Lolium perenne L.)[END_REF][START_REF] Fitter | Root production, turnover and respiration under two grassland types along an altitudinal gradient: influence of temperature and solar radiation[END_REF][START_REF] Steinaker | Phenology of fine roots and leaves in forest and grassland[END_REF][START_REF] Garcia-Pausas | Seasonal patterns of belowground biomass and productivity in mountain grasslands in the Pyrenees[END_REF]. Global radiation and air temperature are the main abiotic determinants of root growth in perennial grasslands not limited by water inputs [START_REF] Edwards | Root production is determined by radiation flux in a temperate grassland community[END_REF]), because of the rapid transfer of photosynthates to roots in grassland ecosystem [START_REF] Kuzyakov | REVIEW: time lag between photosynthesis and carbon dioxide efflux from soil: a review of mechanisms and controls[END_REF]. This is confirmed by the fact that the year 2006, for which we measured higher standing root length, was warmer (+0.7 °C) and more sunny (+10 %) particularly during the growing season compared with 2007 and 2008. However disentangling the respective role of radiation flux and of temperature is difficult due to their tight relationship. In addition, our results may reflect plant ontogeny patterns (i.e. species-specific differences in phenology, [START_REF] Pontes | Seasonal productivity and nutritive value of temperate grasses found in semi-natural pastures in Europe: responses to cutting frequency and n supply[END_REF]), as it is likely that not all species responded at the same time and at the same rate to the climate treatments applied here. While species number was unaffected by the treatments, abundance of grasses declined from 98 % to 89 % in favour of legumes in the warmed site in comparison with the control site [START_REF] Cantarel | Four years of simulated climate change reduces aboveground productivity and alters functional diversity in grassland[END_REF]. These changes are expected to affect ecosystem functioning both at above and below-ground levels, but due to the impossibility to recognise species or functional group from minirhizotron images, it was not possible to assess directly its effect on root dynamics. Finally, this data set shows that grassland root system response to multiple climate change drivers is mainly driven by intra and inter-annual climate variability, which may in turn control the dynamic of plant functional diversity activity and confirms the third hypothesis.

Conclusions

Here, we report changes in root dynamics, longevity and morphology of grassland monoliths submitted to air warming and simultaneous application of warming, summer precipitation reduction and elevated CO 2 over 3 years. Our results suggest that under warmer climates and the future IPCC A2 scenario, upland grassland ecosystems may increase their total root system length but lose a proportion of their fine roots growing in the upper soil layer. This underlines the importance of studying root diameter class responses in climate change experiments. Our finding also show that climate variability, both at inter-and intra-annual scales, mediates grassland root responses under both current and future climate scenarios, and may mask clear patterns of response to climate drivers. Future studies should examine whether climate-induced impacts of root dynamics may increase soil C storage in the long term and for deeper soil layers.

Fig. 1

 1 Fig.1Effects of climate treatment on seasonal variations in volumetric soil water content (0-20 cm). Treatment codes are as follows: C 0 control; T 0 elevated temperature (+3.5 °C), TD 0 elevated temperature and summer precipitation reduction (-20 %); TDCO 2 0 elevated temperature, summer precipitation

Table 1

 1 Effects of warming, summer precipitation reduction, elevated CO 2 and combined climate change drivers (2080 climate scenario) on root length growth rate, cumulative new root length, mortality and standing root length measured during the three-year study period (all years). Treatment codes are as follows: C 0 control; T 0 elevated temperature (+3.5 °C), TD 0 elevated temperature and summer precipitation reduction (-20 %); TDCO 2 0 elevated temperature, summer precipitation reduction and CO 2 enrichment (+200 μmol mol -1 )

	Plant Soil (2013) 364:395-408
		Year	Warming Summer precipitation reduction Elevated CO 2 2080 climate scenario
			C vs T		T vs TD		TD vs TDCO 2 C vs TDCO 2
			F	P	F	P	F	P	F	P
	Root growth rate (mm cm -2 day -1 ) 2006	0.69 NS 0.01	NS	0.19	NS	1.43	NS
		2007	3.85 *	0.15	NS	0.09	NS	4.05	*
		2008	7.39 ** 0.13	NS	0.75	NS	3.11	NS
		All years 10.05 ** 0.05	NS	0.15	NS	7.19	**
	New root length (mm cm -2 )	2006	0.07 NS 0.47	NS	0	NS	0.37	NS
		2007	3.41 NS 0.01	NS	2.4	NS	10.2	*
		2008	5.35 *	0.33	NS	1.02	NS	7.22	*
		All years 5.25 *	0.02	NS	0.98	NS	7.45	*
	Mortality root length (mm cm -2 )	2006	0.51 NS 0.62	NS	0.49	NS	0.86	NS
		2007	0.12 NS 0.36	NS	0.85	NS	0.46	NS
		2008	0.03 NS 2.69	NS	0.13	NS	0.5	NS
		All years 0.13 NS 0.33	NS	0	NS	0.79	NS
	Standing root length (mm cm -2 )	2006	0.07 NS 0.07	NS	0.6	NS	0.03	NS
		2007	0.87 NS 0.41	NS	0.01	NS	4.34	NS
		2008	4.71 NS 3.18	NS	5.72	*	8.09	*
		All years 6.31 *	1.01	NS	3.9	NS	19.4	**
	F is F-ratio, NS: P>0.05; *: P≤0.05; **: P≤0.01

Table 2

 2 Effects of climate treatment, root diameter classes and year on cumulative new root length, root length mortality and standing root length (mm cm -2 ) measured during the three-year study period For each class, all roots and each treatment, different letters correspond to statistical differences (NS: P>0.05; *: P≤0.05; **: P≤0.01.). Uppercase and lowercase letters correspond to comparison of diameter class per treatment and comparison of treatments per diameter class, respectively

	Source	d.f.	F	New root length	F	Mortality root length	F	Standing root length
	Treatment	3	5.3	***	1.2	NS	3.0	*
	Diameter classes	2	123.9	***	43.2	***	54.4	***
	Year	2	7.0	**	31.5	***	30.5	***
	Treatment x Classes	6	3.0	***	4.3	***	0.7	NS
	Treatment x Year	6	0.7	NS	0.9	NS	1.7	NS
	Classes x Year	4	19.5	***	9.8	***	11.7	***
	Treatment x Classes x Year	12	0.5	NS	0.6	NS	1.0	NS
	Significances of treatment, diameter classes, year and their interactions are indicated for the different variables; NS: P>0.05; *: P≤0.05;
	***: P≤0.001. d.f.: degree of freedom, F: F-ratio

Table 4

 4 Proportional hazards regression performed on pooled root lifespan data collected from 2006 to 2008. The parameter estimate and chi-square statistic for each model covariate is shown, for the effects of all treatments Parameter estimates (β) can be interpreted as follows: for quantitative covariates such as air temperature or volumetric soil water content, a one-unit change in the covariate is associated with a 100[exp(β)-1] percent change in the hazard of mortality, controlling for other covariates. For dichotomous categorical covariates (e.g. coded 0 and 1, season of appearance)

	Variable	d.f.	Parameter estimate	Standard error	χ 2 value	P-value	Hazard ratio
	All treatments	1	-0.040	0.019	4.35	0.037	0.960
	Year of appearance	1	0.857	0.033	676.41	<.0001	2.356
	Season of appearance	1	0.038	0.039	0.92	0.336	1.039
	Diameter classes	1	-0.313	0.029	112.80	<.0001	0.731
	Air temperature	1	0.005	0.003	2.71	0.099	1.005
	Volumetric soil water content	1	-0.016	0.007	6.16	0.013	0.984
	Significant P-values are shown in bold (P<0.05).

Plant Soil (2013) 364:395-408

PlantSoil (2013) 364:395-408

Acknowledgements We thank Patrick Pichon, Nicolas Deltroy, Amélie Cantarel for their help with data collection and site management during the study period. Katja Klumpp is acknowledged for improving a previous version of the paper. This study received financial support via a doctoral fellowship from the Region Auvergne and the INRA (EFPA Department) awarded to RP, a CNRS postdoctoral fellowship awarded to JMGB, an IFB-GICC project grant, an ANR [French national research agency] project grant, and the EC FP6 'NitroEurope-IP' project.


Bouziane Moumen
email: bouziane.moumen@jouy.inra.fr

Christophe Nguen-The

Alexei Sorokin
email: alexei.sorokine@jouy.inra.fr

Sequence Analysis of Inducible Prophage phIS3501 Integrated into the Haemolysin II Gene of Bacillus thuringiensis var israelensis ATCC35646

Diarrheic food poisoning by bacteria of the Bacillus cereus group is mostly due to several toxins encoded in the genomes. One of them, cytotoxin K, was recently identified as responsible for severe necrotic syndromes. Cytotoxin K is similar to a class of proteins encoded by genes usually annotated as haemolysin II (hlyII) in the majority of genomes of the B. cereus group. The partially sequenced genome of Bacillus thuringiensis var israelensis ATCC35646 contains several potentially induced prophages, one of them integrated into the hlyII gene. We determined the complete sequence and established the genomic organization of this prophage-designated phIS3501. During induction of excision of this prophage with mitomycin C, intact hlyII gene is formed, thus providing to cells a genetic ability to synthesize the active toxin. Therefore, this prophage, upon its excision, can be implicated in the regulation of synthesis of the active toxin and thus in the virulence of bacterial host. A generality of selection for such systems in bacterial pathogens is indicated by the similarity of this genetic arrangement to that of Staphylococcus aureus β-haemolysin.

Introduction

Many bacterial strains of the B. cereus group are pathogenic to different eukaryotic organisms, including animals, insects, and nematodes [START_REF] Lapidus | Extending the Bacillus cereus group genomics to putative food-borne pathogens of different toxicity[END_REF][START_REF] Ehling-Schulz | Bacillus cereus, the causative agent of an emetic type of food-borne illness[END_REF][START_REF] Kotiranta | Epidemiology and pathogenesis of Bacillus cereus infections[END_REF][START_REF] Mock | Anthrax[END_REF][START_REF] Rasko | Genomics of the Bacillus cereus group of organisms[END_REF][START_REF] Schnepf | Bacillus thuringiensis and its pesticidal crystal proteins[END_REF]. The caused illnesses are mainly attributed to the synthesis of toxins and protective cellular structures, usually encoded by plasmids or, in the cases of diarrheic food intoxications, in the chromosome. In addition to the importance of plasmids, carrying the toxins, it was also suggested that the temperate phages can be involved in the adaptation of these bacteria to animal hosts [START_REF] Lapidus | Extending the Bacillus cereus group genomics to putative food-borne pathogens of different toxicity[END_REF][START_REF] Schuch | The secret life of the anthrax agent Bacillus anthracis: bacteriophage-mediated ecological adaptations[END_REF]. It is indicative in this respect that several sequenced genomes of these bacteria possess large extrachromosomal elements encoding plasmid-related and phage-related functions [START_REF] Lapidus | Extending the Bacillus cereus group genomics to putative food-borne pathogens of different toxicity[END_REF][START_REF] Rasko | Genomics of the Bacillus cereus group of organisms[END_REF][START_REF] Han | Pathogenomic sequence analysis of Bacillus cereus and Bacillus thuringiensis isolates closely related to Bacillus anthracis[END_REF][START_REF] Rasko | The genome sequence of Bacillus cereus ATCC 10987 reveals metabolic adaptations and a large plasmid related to Bacillus anthracis pXO1[END_REF][START_REF] Tourasse | The Bacillus cereus group: novel aspects of population structure and genome dynamics[END_REF][START_REF] Hoffmaster | Identification of anthrax toxin genes in a Bacillus cereus associated with an illness resembling inhalation anthrax[END_REF][START_REF] Challacombe | The complete genome sequence of Bacillus thuringiensis Al Hakam[END_REF]. The pathogen evolution can thus be regarded as a constant dynamic exchange of genes between plasmids and temperate phages integrated or not into the bacterial chromosome. A notable illustration of such prophageplasmid coevolution is the similarity between the genome of a large phage 0305f8-36, isolated from B. thuringiensis and a contig of genome of the strain B. weihenstephanensis KBAB4 [START_REF] Thomas | Complete genomic sequence and mass spectrometric analysis of highly diverse, atypical Bacillus thuringiensis phage 0305Φ8-36[END_REF][START_REF] Hardies | Comparative genomics of Bacillus thuringiensis phage 0305ϕ8-36: defining patterns of descent in a novel ancient phage lineage[END_REF]. In fact, this contig corresponds to the 417 kb extrachromosomal element pBWB401 and could therefore be regarded as either a plasmid or a nonintegrated prophage. Still no evidence exists for the association between phages and toxicity genes in the B. cereus group, although the presence in plasmids of genes for highly effective toxins, like that of anthrax, entomotoxic crystal protein, or emetic toxin, was well documented [4-6, 11, 15]. A recent report relates the evolution of temperate phages and the regulatory elements involved in adaptation of B. anthracis to the animal host [START_REF] Schuch | The secret life of the anthrax agent Bacillus anthracis: bacteriophage-mediated ecological adaptations[END_REF].

A partial sequence of the B. thuringiensis var. israelensis ATCC35646 genome was deposited earlier as a "permanent draft" to the NCBI Entrez database [START_REF] Anderson | Comparative genome analysis of Bacillus cereus group genomes with Bacillus subtilis[END_REF]. The future use of these data would rather envisage looking for the answers to defined biological questions, instead of completing of the whole genome. In particular, our interest was related to the prophages integrated into this genome. We found during analysis of the ATCC35646 contigs, corresponding to 3,500-3,800 kb of the 5.4 Mb B. cereus ATCC14579 genome sequenced earlier [START_REF] Ivanova | Genome sequence of Bacillus cereus and comparative analysis with Bacillus anthracis[END_REF], that several temperate phages were integrated in this area into the genome of ATCC35646. One of these prophages has appeared to be particularly interesting, since it was inserted into the gene hlyII, encoding haemolysin II, characterized in several B. cereus strains [START_REF] Andreeva | Purification and cytotoxic properties of Bacillus cereus hemolysin II[END_REF][START_REF] Baida | Complete nucleotide sequence and molecular characterization of hemolysin II gene from Bacillus cereus[END_REF][START_REF] Budarina | A new Bacillus cereus DNA-binding protein, HlyIIR, negatively regulates expression of B. cereus haemolysin II[END_REF][START_REF] Shadrin | Bacillus cereus pore-forming toxins hemolysin II and cytotoxin K: polymorphism and distribution of genes among representatives of the cereus group[END_REF]. This toxin belongs to the same family of proteins as cytotoxin K, responsible for important cases of food poisoning [START_REF] Lund | A new cytotoxin from Bacillus cereus that may cause necrotic enteritis[END_REF][START_REF] Fagerlund | Toxin production in a rare and genetically remote cluster of strains of the Bacillus cereus group[END_REF][START_REF] Fagerlund | Genetic and functional analysis of the cytK family of genes in Bacillus cereus[END_REF]. The excision of such a prophage can lead to the formation of active toxin gene and thus, the prophage can be implicated in the toxin synthesis regulation. We therefore decided to identify the contigs and complete the sequencing of the region corresponding to this prophage. The completed prophage sequence permits to derive the genomic organization of this prophage. We also tested if the excision of the corresponding phage DNA can take place in this strain.

Results

Detection of Multiple Prophages in the B. thuringiensis ATCC35646 Genome and Nucleotide Sequence of phIS3501

Integrated into the hlyII Gene. The assembly of B. thuringiensis var israelensis ATCC35646 genome, deposited in NCBI Entrez, consisted of 866 contigs with accession numbers NZ AAJM01000001 to NZ AAJM01000866 [START_REF] Anderson | Comparative genome analysis of Bacillus cereus group genomes with Bacillus subtilis[END_REF]. The sequencing coverage was estimated to be ∼6.2-fold and only contigs longer than 1500 bp were used for comparison to other genomes [START_REF] Anderson | Comparative genome analysis of Bacillus cereus group genomes with Bacillus subtilis[END_REF]. Like most bacteria of the B. cereus group, the genome of this strain encodes diarrheic toxin components: cytotoxin K (the gene cytK, locus RBTH0664, contig #1900), NheA (RBTH01882, contig #1388), NheB (RBTH01881, contig #1388), NheC (RBTH01258, contig #1255), HblC (RBTH03191, contig #1573), and HblA (RBTH03163, contig #1589). The gene encoding HblD must be located in the noncovered by sequencing regions. The paralogous to the cytK gene, hlyII, was found in this strain interrupted by a cluster of genes coding for several phagerelated proteins. Two parts of this gene, N-terminal (locus RBTH03378) and C-terminal (RBTH01357), were separated on two contigs (##1604 and 1648, resp.). We found this situation interesting since the excision of this prophage should lead to the formation of the entire functionally active hlyII gene. This arrangement resembles sporulation sigK gene of B. subtilis [START_REF] Stragier | Chromosomal rearrangement generating a composite gene for a developmental transcriptional factor[END_REF][START_REF] Kroos | Switch protein alters specificity of RNA polymerase containing a compartmentspecific sigma factor[END_REF][START_REF] Kunkel | Temporal and spatial control of the mother-cell regulatory gene spoIIID of Bacillus subtilis[END_REF][START_REF] Kunkel | The Bacillus subtilis gene for the developmental transcription factor σ(K) is generated by excision of a dispensable DNA element containing a sporulation recombinase gene[END_REF] and Clostridia [START_REF] Haraldsen | Efficient sporulation in Clostridium difficile requires disruption of the σ K gene[END_REF] and βhaemolysin gene hlb of S. aureus [START_REF] Baba | Genome sequence of Staphylococcus aureus strain newman and comparative analysis of staphylococcal genomes: polymorphism and evolution of two major pathogenicity islands[END_REF][START_REF] Bae | Prophages of Staphylococcus aureus Newman and their contribution to virulence[END_REF][START_REF] Goerke | Extensive phage dynamics in Staphylococcus aureus contributes to adaptation to the human host during infection[END_REF].

Finishing of the sequencing of this prophage was done as described in Methods. The completed phIS3501 genome has the size of 44,401 bp and G+C content of 34.9%. The Figure 1(a) represents experimental data that confirm the correctness of sequence assembly. For this purpose, we used LR PCR amplification with primers specific to several regions of the prophage. The coincidence of the lengths of amplified products with that predicted from nucleotide sequence indicated the correctness of the assembling. This verification was necessary since the assembly of the B. thuringiensis ATCC35646 genome, available from NCBI, contains other contigs with high similarity to phIS3501. This similarity reflects the existence of multiple similar prophages inserted into the host genome.

A total of 52 protein-coding genes, varying in size from 53 to 1344 amino acids, were identified as described in Methods. Figure 1(b) (see also Supplementary Table 1 available on line at doi: 10.1155/2012/543286) presents the functional map of phIS3501. Similarity search of the completed sequence of phIS3501 against the 866 contigs deposited in NCBI revealed 19 DNA-DNA hits with the scores less than e -13 , corresponding to complete contigs or their parts with identities on the nucleotide level ranging from 80 to 99%. Apart from eight contigs corresponding to the prophage phIS3501 and having the identity scores of 97 to 99%, the search revealed parts of 11 contigs scoring from 80 to 94% of identity (Figure 1(c)). These hits represent the regions of contigs corresponding to other similar prophages. In the same way, we tested similarities to other completely sequenced phages or prophages of the B. cereus group. These included the Gamma phage and four Lambda prophages of B. anthracis [START_REF] Sozhamannan | The Bacillus anthracis chromosome contains four conserved, excision-proficient, putative prophages[END_REF][START_REF] Fouts | Sequencing Bacillus anthracis typing phages gamma and cherry reveals a common ancestry[END_REF], phBC6A51, and phBC6A52 of B. cereus ATCC14579 T [START_REF] Ivanova | Genome sequence of Bacillus cereus and comparative analysis with Bacillus anthracis[END_REF] and also phBC391A1 and phBC391A2 of "B. cytotoxicus" NVH391-98 [START_REF] Auger | The genetically remote pathogenic strain NVH391-98 of the Bacillus cereus group is representative of a cluster of thermophilic strains[END_REF]. Although the set of phages used in this analysis was not exhaustive, it revealed 13 additional contigs corresponding to prophages (not shown). Six of the identified phage-related contigs, with scaffold numbers 1740, 1746, 1749, 1759, 1797, and 1824 in Entrez, contained integration-replication gene clusters. Several other contigs containing significant phage-related gene clusters were also detected. Thus four contigs (##1427, 1666, 1720 and a part of the mentioned above #1824) corresponded to the phage DNA-packaging gene clusters. Four others (##1436, 1576, 1601, and 1880) contained the lysis module genes. Although, because of sequence gaps in the available shotgun assembly, it was not possible to join together the different functional gene clusters, we concluded that at least six other similar prophages could be found in the genome of B. thuringiensis ATCC35646. An attempt to make such estimation using individual phage proteins did not appear to be successful due to ambiguity during interpretation (not shown). During the time that this paper was under preparation, a partial sequence of the strain B. thuringiensis IBL4222 isolated from a cat has appeared in NCBI Entrez database (accession identifier CM000759.1) in the form of 383 contigs. It appeared that this strain is in fact also B. thuringiensis var israelensis and five of the contigs corresponded to phIS3501 having higher than 99.9% identity to our prophage sequence (Supplementary Table 2).

For 24 of all proteins encoded in phIS3501 (46% of total number), we have found significant similarity with proteins of known biochemical functions (Figure 1(b)). For six of these proteins, the precise biological function in phage development cannot be assigned. These are four proteins involved in phage development regulation (p04, p07, p21 and p29), integrase p22, and FtsK family protein p50. In grey are the genes with hypothetical functions. Other colours correspond to the functional phage modules: red (p01, p04, p05 and p22) and yellow (p07, p09 and p21), lysogeny and lysogenic regulation, brown (p10 and p11), replication; dark green (p29-p32)-DNA packaging and maturation; green, blue and purple (p33-p36, p39, and p41-p43)-head and tail structural module; violet (p45 and p46), host lysis; and cyan (p50) indicates the FtsK family protein, probably involved in lysogenic recombination. Scale bar is in bp. (c) Contigs from Genbank entries NZ AAJM01000001-NZ AAJM01000866 revealed by BLASTN search against the completed phIS3501 sequence. The names of significant contigs with identities 80 to 94% are shown in italics. Other contigs, with the names shown in bold, have the identity of 97-99% and correspond to the phage. The entire phage sequence is represented on the top as the grey scale bar in bp.

For 28 (54%) of encoded proteins we, cannot predict a function, although 24 of them share similarity with other phage proteins and therefore can be regarded as conserved.

For some of these hypothetical proteins putative phagerelated biological functions can be postulated based on their locations in the phage genome (Figure 1(b)). No similarity was detected for four (8%) of predicted proteins. The genome of phIS3501 also contains tRNA-Met gene (not shown) located about 200 bp upstream of the gene p01 for the lysogenic integrase. In this region, no genes-encoding proteins were predicted. This location is interesting since it can suggest the involvement of this tRNA gene in the regulation of integrase expression and thus of the prophage DNA integration or excision. Five functional modules or gene clusters can be recognized in the genomic organization of phIS3501. These are lysogeny and lysogenic regulation module, replication, DNA packaging and maturation, and head and tail structural module and lysis (Figure 1(b)). The modular organization of these genes can prove to be useful in postulating of biological functions for the hypothetical proteins [START_REF] Canchaya | Prophage genomics[END_REF]. It is interesting to mention that most of the similar phages of the B. cereus group contain in their genome an additional module, which includes the gene encoding for an ATPdependent DNA transporter of FtsK family. The role of this protein in phage biology was not yet experimentally studied and is unclear at present. However, the vicinity of this gene, in the autonomously replicating phage DNA, to lysogenic integrase suggests that this module can be involved in phage lysogeny. The lysogeny and lysogenic regulation module of phIS3501 contains nine genes (p01-p09), the most evident functions are the integrase (p01) and the putative lysogenic repressor (p05). Integrase is the protein essential for the specific integration and, together with excisionase, which we cannot reliably recognize, excision of this phage. The abovementioned ftsK-like gene can also be involved in this process, facilitating the integrase-mediated reaction. The lysogenic regulation involves also the repressor gene (p05) whose disruption would result in the clear-plaque phenotype.

In addition, the lysogenic regulation module includes also a helix-turn-helix protein and antirepressor (p04 and p07, resp.). The precise role of these genes in the regulation of similar phages of the B. cereus group was not yet experimentally studied. Based on its location, the tRNA-Met gene can also be involved in the lysogeny regulation. Insufficiency of sequence data on experimentally characterized excisionases from Gram-positive bacteria does not allow detecting by similarity search of the potential excisionase. Its function is needed to direct the integrase reaction out of the bacterial chromosome. Nevertheless, the appropriate localization and the similarity of size (85 aminoacids) and of content of positively charged aminoacids (15%) suggest that this could be the protein encoded by the p03 gene.

Another two regulatory proteins of this phage are encoded by p09 and p21. These are the RNA-polymerase sigma factor and ArpU-family regulator, the former could be involved in the regulation of late phage genes, as it was demonstrated for the Fah phage of B. anthracis [START_REF] Minakhin | Genome sequence and gene expression of Bacillus anthracis bacteriophage Fah[END_REF], similar in its genomic organization to phIS3501. It is not clear, what is the exact role of the second integrase (p22), although it cannot be excluded that this protein is involved in the lysogenic conversion of the phage, participating in acquisition of functions useful for the bacterial host [START_REF] Schuch | The secret life of the anthrax agent Bacillus anthracis: bacteriophage-mediated ecological adaptations[END_REF]. The replication module of phIS3501 contains two genes encoding DNA replication initiator (p11) and primosome loader (p10). These proteins are very characteristic for the temperate phages of the B. cereus group and their counterparts are found in most of the prophages residing in the sequenced genomes [START_REF] Rasko | Genomics of the Bacillus cereus group of organisms[END_REF].

The phylogenic position of proteins encoded in the packaging module can indicate the type of DNA processing during the phage packaging. Large subunit of terminase in phIS3501 (p31) clusters with terminases characteristic to the phages this DNA maturation is driven by the recognition of cossites [START_REF] Casjens | The generalized transducing Salmonella bacteriophage ES18: complete genome sequence and DNA packaging strategy[END_REF]. The most similar protein from the phage genomes is that of Geobacillus phage E2. Highly similar are also the terminases of B. subtilis phage phi105 and of B. anthracis phages Cherry, Gamma, Wβ and Fah. The phIS3501 can therefore be positioned among the phages using the cos-type of DNA-packaging initiation.

The head and tail module genes of tailed phages define, in a great extent, the specificity of phages to their bacterial hosts. Although in the phage phIS3501, this module includes more than ten genes, only two of them, p42 and p43, are supposed to be directly involved in the bacterial host recognition. The orthologues of these two genes, encoding tail fiber protein and a minor structural protein, respectively, were shown to be involved in specific binding of similar phages from B. anthracis to the host [START_REF] Schuch | Detailed genomic analysis of the Wβ and γ phages infecting Bacillus anthracis: implications for evolution of environmental fitness and antibiotic resistance[END_REF]. Specific binding of the counterpart of p42 from the phage Wβ was experimentally demonstrated using fusion with the Green Fluorescent Protein (GFP). Involvement of the p43 counterpart was postulated due to finding of many mutations changing the phage specificity [START_REF] Schuch | Detailed genomic analysis of the Wβ and γ phages infecting Bacillus anthracis: implications for evolution of environmental fitness and antibiotic resistance[END_REF]. It was also shown that the host cell encoded protein GamR is one of essential components for binding of similar phages, presumably the phage receptor [START_REF] Davison | Identification of the Bacillus anthracis γ phage receptor[END_REF].

The lysis module genes, in particular the phage endolysin, were also shown to be highly related to the bacterial host specificity, suggesting their great usefulness for diagnostic and therapeutic purposes [START_REF] Schuch | A bacteriolytic agent that detects and kills Bacillus anthracis[END_REF][START_REF] Fischetti | Reinventing phage therapy: are the parts greater than the sum?[END_REF]. It appeared that the p46 lysine gene of phIS3501 contained a frame shift in the fifth codon, thus leading to encode a nonfunctional protein (not shown). Therefore, the synthesized phage phIS3501 is not able to autonomously finish the lytic phage cycle, due to its inability to provide the lysis of bacterial host. However, since the bacterial genome encodes six or more of other similar prophages, it is possible that this function can be accomplished by one of them.

We concluded from the genomic sequence that the region encoding phIS3501 is able to provide autonomous phage DNA replication, excision, and synthesis of structural components, but it cannot provide the autonomous release of the phage particles, if they are formed, from the bacterial cells.

Induction of Lytic Development of Temperate Phages in B. thuringiensis ATCC35646.

The induction experiments reported here have the goal to show if the genomic region that we designated phIS3501 is able to be induced for selfreplication and to be excised from the genome. Also, we have tried to test whether this sequence corresponds to an inducible prophage that can be detected by using a sensible bacterial host. The phage DNA was tested by PCR, using specific oligonucleotides, and the biological activity of the phage was tested using several strains of bacteria phylogenetically closely related to the strain B. thuringiensis ATCC35646. The mitomycin C induction experiment is represented on Figure 2(a). Total DNA was extracted from the supernatant of lysed bacterial culture and from the noncompletely lysed cells collected by centrifugation. PCR amplification experiments and sequencing of the corresponding products show that phIS3501 is able to induce the synthesis of the replicative form of its DNA, and to excise it from the host chromosome (Figure 2(b) and 2(c)). This form can be detected by PCR even in the cells nontreated with mitomycin C (Figure 2(c)). The excision of this prophage from the chromosome results in formation of the wholelength hlyII gene, thus enabling the cells to produce this toxin. Twenty strains, closely related to B. thuringiensis ATCC35646, from recently described collection [START_REF] Sorokin | Multiplelocus sequence typing analysis of Bacillus cereus and Bacillus thuringiensis reveals separate clustering and a distinct population structure of psychrotrophic strains[END_REF], were tested for the formation of phage-related plaques. However, we did not succeed in finding any sensitive host to be used as indicator (data not shown). This is surprising since, as reported above, the strain B. thuringiensis ATCC35646 contains at least seven different prophages in the chromosome and also possesses the nonintegrated linear phage-like element, presumably also corresponding to a phage [START_REF] Anderson | Comparative genome analysis of Bacillus cereus group genomes with Bacillus subtilis[END_REF]. We consider this result as not definitive and rather due to our ignorance of proper conditions for obtaining the plaques for these phages and hosts.

Analysis of sequences of the circularized form of phage DNA and its integration site in the bacterial chromosome permitted us to identify the 13 bp sequence common in the phage DNA and in the chromosome (Figure 2(d)). This site-specific phage integrase p01.

Discussion

Many bacterial sequencing projects do not result in establishing the completed sequences of studied genomes. The relative amount of such data will certainly increase with the advent of extrahigh throughput methods. The reason is that, even if the finished data represent much higher value than the "skimmed" ones, the so-called finishing is a laborious enterprise, containing a significant amount of manual work that can hardly be automated. Comparative utility of both approaches was discussed [START_REF] Fraser | The value of complete microbial genome sequencing (you get what you pay for)[END_REF][START_REF] Branscomb | On the high value of low standards[END_REF][START_REF] Gardner | Draft versus finished sequence data for DNA and protein diagnostic signature development[END_REF][START_REF] Selkov | Functional analysis of gapped microbial genomes: amino acid metabolism of Thiobacillus ferrooxidans[END_REF] but it is certain that the investment in precise experimental data acquisition requires case-to-case justification and cannot be generalized. The work that we present here uses previously assembled shotgun data deposited into public databases [START_REF] Anderson | Comparative genome analysis of Bacillus cereus group genomes with Bacillus subtilis[END_REF]. The initial complete assembly of the 40 kb long prophage was hampered by insufficiency of experimental data and also by the complexity that was due to existence of highly similar, but different, sequences in the same genome. We completed the sequence of this prophage that permitted us to understand the extent of its integrity. For this purpose, we used the similarity of organization of such prophages to those already entirely sequenced. This revealed the potential sequencing contigs covering this region and thus minimized the subsequent combinatorial PCR experiment. Since we detected also the existence of at least six other similar prophages in this genome, the independent verification of the final assembly of this area was applied. Finally, our complete prophage assembly was also confirmed by independent shotgun sequencing and assembly of a strain B. thuringiensis IBL4222, which appeared to be another sequenced israelensis strain (NCBI accession identifier CM000759.1).

This prophage was chosen to be entirely sequenced since it is integrated into an important haemolysin gene. The availability of completed prophage sequence permitted us to make preliminary conclusions of its functionality and to test it experimentally. We demonstrated the ability of this prophage to induce its DNA excision and ligation. However, since the endolysin gene of the prophage is interrupted by internal frameshift, we cannot expect formation of the mature phage particles without functional involvement of other induced prophages. The excision of this sequence from the bacterial chromosome leads to the formation of uninterrupted hlyII gene, encoding a potentially active haemolysin. We did not, however, find the experimental conditions to demonstrate the synthesis of this protein (not shown).

Two other intensively studied systems can be mentioned that represent a similar mode of gene regulation by creating of uninterrupted coding sequence upon a temperate phage excision. The first is the well-known skin element, inserted into the coding sequence of sigK gene of B. subtilis [START_REF] Stragier | Chromosomal rearrangement generating a composite gene for a developmental transcriptional factor[END_REF][START_REF] Kroos | Switch protein alters specificity of RNA polymerase containing a compartmentspecific sigma factor[END_REF][START_REF] Kunkel | Temporal and spatial control of the mother-cell regulatory gene spoIIID of Bacillus subtilis[END_REF][START_REF] Kunkel | The Bacillus subtilis gene for the developmental transcription factor σ(K) is generated by excision of a dispensable DNA element containing a sporulation recombinase gene[END_REF]. Genomic sequencing revealed that the skin element is in fact a prophage that lost many functions essential for formation of phage particles [START_REF] Takemaru | Complete nucleotide sequence of a skin element excised by DNA rearrangement during sporulation in Bacillus subtilis[END_REF]. Actually, the skindependent regulation of this sigma-factor synthesis is not essential for sporulation of B. subtilis and such insertions do not exist in other Bacillus species [START_REF] Kunkel | The Bacillus subtilis gene for the developmental transcription factor σ(K) is generated by excision of a dispensable DNA element containing a sporulation recombinase gene[END_REF][START_REF] Adams | Molecular cloning and characterization of two genes encoding sigma factors that direct transcription from a Bacillus thuringiensis crystal protein gene promoter[END_REF][START_REF] Takami | Complete genome sequence of the alkaliphilic bacterium Bacillus halodurans and genomic sequence comparison with Bacillus subtilis[END_REF]. However, such a construction is important for the correct sporulation timing in Clostridia [START_REF] Haraldsen | Efficient sporulation in Clostridium difficile requires disruption of the σ K gene[END_REF]. Another remarkable example of gene regulation by prophage excision is the haemolysin gene of Staphylococcus aureus. The phage φNM3 is integrated into the β-haemolysin gene (hlb) of many S. aureus strains [START_REF] Baba | Genome sequence of Staphylococcus aureus strain newman and comparative analysis of staphylococcal genomes: polymorphism and evolution of two major pathogenicity islands[END_REF][START_REF] Bae | Prophages of Staphylococcus aureus Newman and their contribution to virulence[END_REF][START_REF] Goerke | Extensive phage dynamics in Staphylococcus aureus contributes to adaptation to the human host during infection[END_REF]. Since φNM3 carries several virulence factors, its induction not only leads to the formation of active haemolysin gene but also increases its physiological effects by weakening the host immune system. Surprisingly, the endolysin gene in this prophage is also mutated [START_REF] Bae | Prophages of Staphylococcus aureus Newman and their contribution to virulence[END_REF][START_REF] Goerke | Ciprofloxacin and trimethoprim cause phage induction and virulence modulation in Staphylococcus aureus[END_REF]. In the case of the prophage phIS3501, described here, we did not detect any obvious virulence factor associated with the prophage. However, since six other similar prophages are integrated in the host genome, some of them, being simultaneously induced, may provide such factors.

At present, the involvement of hlyII gene of the B. cereus group bacteria in virulence for animals is not entirely clear. Several papers report that this system can be rather important. A pathogenicity-related protein of this family, cytotoxin K, was described in the remote strain "B. cytotoxicus" NVH391-98 and considered as the main factor responsible for toxicity in a case of collective food poisoning [START_REF] Lund | A new cytotoxin from Bacillus cereus that may cause necrotic enteritis[END_REF]. The toxicity was regarded as to be related to a particular allele of this gene, the others do not seem to be so hazardous due to the weakness of gene expression level and cell lysis ability of the protein [START_REF] Fagerlund | Toxin production in a rare and genetically remote cluster of strains of the Bacillus cereus group[END_REF][START_REF] Fagerlund | Genetic and functional analysis of the cytK family of genes in Bacillus cereus[END_REF][START_REF] Brillard | Comparison of cytotoxin cytK promoters from Bacillus cereus strain ATCC 14579 and from a B. cereus food-poisoning strain[END_REF]. Two other papers demonstrate the increase of virulence of bacterial cells due to synthesis of the hlyII gene product [START_REF] Sineva | Expression of Bacillus cereus hemolysin II in Bacillus subtilis renders the bacteria pathogenic for the crustacean daphnia magna[END_REF][START_REF] Tran | Haemolysin II is a Bacillus cereus virulence factor that induces apoptosis of macrophages[END_REF].

The similarity of the arrangement of prophage phIS3501 associated with the hlyII gene to the one of S. aureus phage φNM3 integrated into the β-haemolysin gene may indicate a general character of selection for such systems. The fact that a similar genetic arrangement is found for sigK of such remote organisms as B. subtilis and C. difficile and for haemolysins of S. aureus and B. thuringiensis indicates that common evolutionary advantages could exist. Therefore the functionality of phIS3501 merits further investigation especially with relatedness to the virulence, if any. An important factor would be development of an experimental system for measuring of pathogenicity of the strain B thuringiensis ATCC35646 or similar in relation to the haemolysin II gene functioning.

Conclusions

The complete nucleotide sequence and genomic organization of the prophage phIS3501 inserted into the hlyII gene of B. thuringiensis var israelensis ATCC35646 was established. The prophage encoded lysine gene contains a frameshift that can prevent formation of phage particles. Nevertheless, the concerted cell lysis that we observed applying mild mitomycin concentrations could be due to induction of one of other six prophages residing elsewhere in the host genome. Excision of phIS3501 from bacterial chromosome leads to formation of uninterrupted hlyII gene, encoding a potentially active haemolysin. A general character of selection for such systems in bacterial pathogens is indicated by the similarity of this gene arrangement to that of Staphylococcus aureus β-haemolysin.

Methods

Bacterial Strains, Growth

Conditions and DNA Manipulations. The strain Bacillus thuringiensis var israelensis ATCC35646, the same that was used for the genomic sequencing [START_REF] Anderson | Comparative genome analysis of Bacillus cereus group genomes with Bacillus subtilis[END_REF], was obtained from Dr. Alla Lapidus (Integrated Genomics). Bacterial growth, total DNA preparation, and PCR conditions were as described [START_REF] Sorokin | Multiplelocus sequence typing analysis of Bacillus cereus and Bacillus thuringiensis reveals separate clustering and a distinct population structure of psychrotrophic strains[END_REF][START_REF] Candelon | Two distinct types of rRNA operons in the Bacillus cereus group[END_REF].

Sequencing of the phIS3501 Genome.

To determine the entire sequence of the prophage phIS3501, we used the data of shotgun assembly, available from the ERGOlight database (http://www.ergolight.com/ERGO/), deposited to NCBI Entrez under the accession numbers NZ AAJM01000001-NZ AAJM01000866 (866 contigs). The contig #1648, containing a part of the gene hlyII, interrupted by the phageencoding DNA, was tested, using the Pinned Regions tool implemented in the ERGOlight database, for colinearity with other B. cereus group genomes containing similar prophages elsewhere in the chromosome. This procedure, interactively applied, permitted to detect five contigs (##1501, 1657, 1422, 754 and 1696 in ERGOlight) presumably corresponding to the phIS3501 prophage. Contiguity these regions in the bacterial chromosome, together with the contigs #1604 and 1648, containing parts of the hlyII gene, was verified by Long Range PCR (LR PCR, Figure 1). Thus, the sequencing substrates were generated which were used to complete the sequence of the entire phIS3501 by primer walking. For that 129 oligonucleotides and 416 additional sequencing runs of 700 bases of average length were produced. Xbap or Gap4 software [START_REF] Dear | A sequence assembly and editing program for efficient management of large projects[END_REF][START_REF] Bonfield | A new DNA sequence assembly program[END_REF] was used for the sequence assembly. The average coverage of the de novo sequenced regions was approximately 7-fold. The completed assembly of phIS3501, once again confirmed by LR PCR, revealed two additional contigs (#1660 and 1024) from the ERGOlight corresponding to this prophage. Supplementary Table 3 presents the oligonucleotides, used for sequence assembly confirmation and in phage induction experiments and their positions in relation to the phage. The de novo annotated sequence of phIS3501 was deposited to NCBI Entrez under the accession number JQ062992.

Phage Induction.

The overnight culture of B. thuringiensis var israelensis ATCC35646 grown at 37 • C in liquid medium under aeration (200 rpm/min) was diluted 100-fold and incubated in the same conditions. Mitomycin C was added after 1 h to the final concentration 0.2 μg/mL. This concentration was found in pilot tests (not shown) to be the minimal needed to induce the prophages without causing the inhibition of bacterial growth. Culture supernatant or total cellular DNA extracted from 1 mL of precipitated by centrifugation cells were used for PCR analysis of the presence of excised phage DNA.

Phage Genome Analysis and Gene

Predictions. The consensus sequence of phI3501 genome, generated by the assembly software, was used for gene prediction by GeneMark program [START_REF] Lukashin | GeneMark.hmm: new solutions for gene finding[END_REF] implemented at web site of the Georgia Institute of Technology (http://opal.biology.gatech.edu/GeneMark/). The start positions of predicted genes were manually scrutinized and corrected, if a better potential Ribosome Binding Site was detected, using the Sequin NCBI tool (http://www.ncbi.nlm.nih.gov/Sequin/). Search for homology of phage-encoded proteins to the NCBI database was done using BLAST tools implemented there or at the LIRMM (Le Laboratoire d'Informatique, de Robotique et de Microélectronique de Montpellier) web-site (http://www.phylogeny.fr/) [START_REF] Dereeper | Phylogeny.fr: robust phylogenetic analysis for the non-specialist[END_REF]. The ACLAME database (http://aclame.ulb.ac.be/) was used to detect the similarity with bacterial virus proteins [START_REF] Leplae | ACLAME: a CLAssification of mobile genetic elements, update 2010[END_REF]. The gene functions were assigned using similarity to known proteins and also taking into account the positions of genes in phage genomes [START_REF] Canchaya | Prophage genomics[END_REF]. tRNAscan-SE [START_REF] Lowe | tRNAscan-SE: a program for improved detection of transfer RNA genes in genomic sequence[END_REF] was used for tRNA prediction.

Figure 1 :

 1 Figure 1: Analysis of the chromosomal region containing phIS3501. (a) Bars on the top show the sequencing contigs from GenBank entries NZ AAJM01000001-NZ AAJM01000866, corresponding to the region encoding phIS3501. Unfilled parts correspond to bacterial sequences, grey parts correspond to the phage contigs detected using ERGOlight, in black are the GenBank contigs assigned to this phage during the sequencing proceeds. The unfilled bars on the bottom show the PCR products seen on the gel to the right, used for the sequence assembly verification. The numbers indicate gel lanes and the expected product sizes, corresponding to the following primer couples: (1) PHISB8 × PHISE8; (2) PHISC3 × PHISE5; (3) PHISB3 × PHISD7; (4) PHISD6 × PHISG8; (5) PHISA1 × PHISH2. A few molecular weight marker sizes (lane M) are indicated to the left of the gel. (b) Sequence-based genetic map of phIS3501. Similarity-or position-based identification of gene functions are shown by short descriptions of encoded proteins. Predicted protein coding genes correspond to the Entrez sequence acc. # JQ062992.In grey are the genes with hypothetical functions. Other colours correspond to the functional phage modules: red (p01, p04, p05 and p22) and yellow (p07, p09 and p21), lysogeny and lysogenic regulation, brown (p10 and p11), replication; dark green (p29-p32)-DNA packaging and maturation; green, blue and purple (p33-p36, p39, and p41-p43)-head and tail structural module; violet (p45 and p46), host lysis; and cyan (p50) indicates the FtsK family protein, probably involved in lysogenic recombination. Scale bar is in bp. (c) Contigs from Genbank entries NZ AAJM01000001-NZ AAJM01000866 revealed by BLASTN search against the completed phIS3501 sequence. The names of significant contigs with identities 80 to 94% are shown in italics. Other contigs, with the names shown in bold, have the identity of 97-99% and correspond to the phage. The entire phage sequence is represented on the top as the grey scale bar in bp.

Figure 2 :

 2 Figure 2: phIS3501 induction. (a) Effect of mitomycin C treatment on growth of B. thuringiensis var israelensis ATCC35646. Optical densities are shown of cultures treated (squares) or untreated (triangles) by mitomycin C (0.2 μg/mL) added at time point of 60 min. (b) Detection of replicative form of phIS3501 DNA. PCR reactions were done using the total DNA extracted from the cell pellet (lanes 2 and 3) or supernatant (lane 4) 20 min after mitomycin C addition. Primers used: lane 2, PHISI2 × PHISI3 (specific to phage DNA); lane 3, PHIST6 × PHIST7 (specific to host chromosomal DNA); lane 4, PHISK6 × PHISS8 (specific to replicative or mature form of phage DNA). A few molecular weight marker sizes (lane 1) are indicated in bp to the left of gel. Obtained product sizes correspond to the ones expected from the sequence. (c) Sequencing tracks corresponding to integrated and excised DNA of phIS3501 and of the host. DNA for PCR amplification was extracted from noninduced overnight culture of B. thuringiensis var israelensis ATCC35646. Cartoons under the sequencing tracks show their interpretation in relevance to the phage integration status. Nter and Cter are the chromosomal parts of the hlyII gene. Phage L and Phage R are the ends of integrated phIS3501. The phage sequences are shown inside of red bar above the sequencing tracks. PCR amplification and (sequencing) primers were, from top to the bottom: PHISC7 × PHISE8 (PHISA2); PHISA2 × PHISC8 (PHISA2); PHISK5 × PHISC5 (PHISC5); PHISK6 × PHISS8 (PHISK6). (d) The bacterial and phage attachment sites. The derived phage (attP) and bacterial (attB) attachment sites based on the data presented in (c). Thirteen bases common sequence, presumably recognized by the phage integrase, is underlined.

Acknowledgments

The authors thank Dr. Alla Lapidus for providing us the strain ATCC35646, Benoit Quinquis, and Nathalie Galleron for the technical help with sequencing reactions. B. Moumen held a thesis fellowship from franco-algerian intergovernmental program. The work was partially supported by the French National Research Agency (project ANR-05-PNRA-013).