Assioua Yasmine
email: yasmine.assioua@renault.com

Ameur-Boulifa Rabea
email: rabea.ameur-boulifa@telecom-paris.fr

Guitton-Ouhamou Patricia
email: patricia.guitton-ouhamou@renault.com

Requirements to Models of Automotive Software: Application to the Automatic Park Assist function

Keywords: Requirements analysis, reliable systems, modelbased design

In the software development lifecycle, errors and flaws can be introduced in the different phases and lead to failures. Establishing a set of functional requirements helps producing safe software. However, ensuring that the (being) developed software is compliant with those requirements is a challenging task due to the lack of automatic and formal means to lead this verification. In this paper, we present our approach that aims at analysing a collection of automotive requirements by using formal methods. The proposed approach for formal verification is evaluated by the application to the Automatic Park Assist (APA) function.

I. INTRODUCTION

Usually, software project involves a huge set of requirements, that can be written by different stakeholders (software and requirements engineers, etc). These requirements are usually expressed in natural language which arises a problem with the expressiveness, the completeness, and the accuracy. The design problem in system engineering is to ensure that the relevant functional and safety requirements are respected at all stages during the development phase. Within that context, many efforts have been done in the processes of software development and in particular in the validation phase. However, although methodologies and tools have been proposed for that purpose, very few works have been devoted to supporting rigorous design, in particular to provide a systematic stepwise design approach for building a validated design model. Such an approach would be valuable for software developer, it should provide confidence that the requirements are consistent and realizable.

The idea behind this paper is to propose a rigorous modelbased process leading from requirements to correct implementations. The advocated approach is applied for designing software for reliable autonomous vehicles, will provide assistance and guidance to the developer in the development phase to assess the correctness of his design. Indeed, our methodology allows to early specify and validate a system design from a collection of functional requirements. In order to highlight and assess our methodology, we applied it on a real industrial case study, that is the the APA (Automatic Park Assists) function. This function is implemented into the vehicle so it can park itself. Inside this function, there is many sub-functions which have different autonomous levels. From the HFP (Hand Free Parking) function, the system only controls the steering wheel to the RPK (Remote Park) function where the driver is outside the vehicle and the system controls all the elements involved in parking the car. The systems requirements are specified in an excel file called a System Technical Requirement Component (STRComp, for short). This file involves all the requirements (functional and safety) and the definitions of a given system. In this particular instance, the APA's STRComp consists of 327 requirements among them functional requirements and safety requirements.

II. FRAMEWORK -FROM REQUIREMENTS TO FORMAL MODELS

Our goal is to use formal methods both to increase the quality of the systems developed at Renault through enhancing the verification activity, and to prevent unnecessary tests. We propose a general model-based approach for the software development activity, that aims to the systematic development of a design solution for a set of system requirements. The big picture of the proposed approach is shown in Figure 1, highlighting the relevant steps towards fulfilling the transformation of requirements written in natural language into exploitable design models that can be automatically verified for errors and validated against behavioural properties.

Step1. Requirements Analysis: Aiming to provide assistance and guidance to the engineers to ensure the quality of their developed software with respect to functional and safety requirements. Our methodology requires proceeding with the specification of the identified requirements. Actually, the analysis is naturally focused on the requirements with clear relevance to the system under design (software under development). A system is defined by a collection of states and the set of relevant requirements are those that can affect states of the collection. Starting from those requirements described in an informal manner, we proceed by analysing them for gathering knowledge and extracting the key elements/keyconcepts from their textual description. Parsing requirements is not a trivial task, as those requirements, are often written by different engineers in various style, and use natural language that is inherently ambiguous, as it is not tied to a formal semantics. This activity should be conducted together with domain experts (requirements engineers). At this final parsing stage unambiguous specification of requirements are derived. Step2. Model Construction: Once all the requirements of a system (software or architecture) have been analysed and the patterns generated, an initial design model is automatically built from all the patterns. Actually, patterns bring a set of states and transitions between them. So the design model in its initial form is a state machine built from all the set of obtained states and transitions. However, as requirements may be incorrect, missing and conflicting, so the set of states may be incomplete and not sound. For instance, some requirement patterns specify only the target state (source state is lacking), so the state machine should be completed, in a first step so that this state is potentially accessible from all other states of the state machine. In a second step, the plausibility of the additional information is checked by using to the priority table -this table is provided by requirement engineers that provides the forbidden transitions by specifying the final states. If the extra-transitions are not coherent with the priority table, we correct the design model. After that, if inconsistencies remain the system software engineer role is required, he can choose among the transitions those that make sense for creating design model that satisfy a set of properties. So the design model is progressively revised and completed until a coherent and complete design model has been obtained.

Step3. Verification: This step aims at automating the systematic verification of compliance with requirements. The model checking method is applied to determine the error-free design of the generated state machine model and automatically find the logical errors. The strength of this method [START_REF] Baier | it Principles of model checking[END_REF] is returning a feasible sequence through a counter-example and simultaneously verifying its correctness using properties. To facilitate the production of tool chains the popular model checker UPPAAL [START_REF] Glarsen | UPPAAL in a nutshell[END_REF] has been used to verify the generated design models. If assumed properties are satisfied, then so is the whole set of requirements; otherwise, the design model should be refined or certain unsatisfied requirements have to be revised.

III. CASE STUDY

The APA system is related to the ADAS (Advanced Driver Assistance System) ECU system. This system assists the driver during the parking maneuver. The automatic park assist contains various sub-functions with different level of autonomy. We have developed an early version of the parser, it has been implemented in Python language.

1) Requirement analysis. An initial step in the parsing is to extract information about states and conditions are extracted from the set of analysed requirements. to Safe State 3 are not plausible, so the model is modified accordingly. In this way, the plausibility rules are incrementally applied to the design model. These rules are those defined in the priority table or provided by software engineer. 5) Model Generation. Where there is more rule to apply, the final design model is generated. We use the UPPAAL model checker tool to display the obtained design model, and especially to verify safety properties that can be expressed as logical formulas [START_REF] Clarke | Automatic verification of finite-state concurrent systems using temporal logic specifications[END_REF] or simply as observer automata [START_REF] Halbwachs | Synchronous observers and the verification of reactive systems[END_REF]. Fig. 2 shows the resulting automaton of our use case, with 6 reachable states and guarded transitions.

To improve the readability of the generated models, we code the names of states and the names of the variables used over transitions in a concise manner. Indeed in practice, during the parsing phase for each identified variable and associated values, a new variable name is created and inserted in a table together with the domain from where it gets its value. Such that for each variable encountered the table is initially consulted to check whether the variable exists or not. More specifically, we simplify the case study encoding in Fig. 2, as shown by denoting the state Safe State with SS and the variable pgen PLC ParkingSequenceStatus with Pstatus.

Additionally, as the pgen PLC ParkingSequenceStatus variable can take 7 different values vgen PLCIdle, vgen PLCRamp, vgen PCLStanby, vgen MovePause, etc, we encode these values using numerical numbers 0, 1, 3, etc.

It is also noted that, still for the sake of readability we encode some conditions by these practical meaning, for instance the condition pgen APA Failure = vgen Fail corresponds to a failed situation, so we denote it by Fail predicate over the associated transition. In the last decade, companies, in particular in the area of aeronautics and aerospace have made great efforts to reduce their development times and to improve the quality of their products. Indeed, for developing complex systems, they are increasingly adopting model-based design tools in system development, e.g. Simulink [START_REF] Anonymous | The MathWorks Home Page[END_REF] and SCADE suite [START_REF] Anonymous | SCADE Suite -Esterel Technologies Home Page[END_REF]. However, the various tools have been developed with the aim of improving the specifications development are mostly focused on requirements management and trace-ability. Unfortunately, even if there is a great interest about the use of formal methods for early validation of system requirements, e.g., [START_REF] Bozzano | Spacecraft early design validation using formal methods[END_REF] and [START_REF] Miller | Proving the shalls: Early validation of requirements through formal methods[END_REF], there are very few requirements validation tools available for checking their correctness and functional consistency before any design or coding fulfills. Among them Argosim [START_REF] Anonymous | The Argosim Home Page[END_REF] that provides practical tool for debugging the requirements, and modelling and simulation capabilities for the validation of systems, from functional requirements engineering to automatic test-case generation. Similarly to the Argosim approach our work aims at a design process that allows for the early validation of system requirements. Important differences from the existing approaches are (i) our solution is automotive-specific. The approach focuses on functional and safety requirements coming from the automotive domain. This means that the validated models can be optimally adapted to automotive projects; (ii) and furthermore, the adoption of formal methods in the analysis process. Indeed, the model checking technique is used as a validation technique in order to verify the correction of system requirements, otherwise to guide their revision.

V. CONCLUSION This paper introduces a systematic process for building design model from functional and safety requirements with the aim of reducing the validation testing during the late stages of the software development lifecycle. The design model provides an early assurance that the requirements specification are complete, correct and realizable. This first effort led to a formalisation of a set of requirements the APA function, and we showed a first proof-of-concept regarding the feasibility of our approach that aims at the exhaustive verification of the generated model; the correctness of the model is proved by using the model checking technique. Naturally, our next goals is to propose a requirements specific language for the expression of requirements that can fit the specific automotive domain. This language with constrained grammar and welldefined semantics will lead to express system requirements in a language close to the natural language, and will make possible to prove formally their correctness.

Fig. 1 .

 1 Fig. 1. Framework for the formal analysis of requirements

Fig. 2 .

 2 Fig. 2. The final model for the APA case study

  Completion. The requirement REQ2 is defined in a complete manner, it provides an initial state, a final state, and a condition enabling the transition. However, some requirements are incomplete, the source state or the destination state are not specified. As for the following requirement: REQ3: IF ADAS receives: -pgen APA Failure= vgen Fail THEN the ADAS ECU shall switch to "APA ADAS ECU Safe State 3". This requirement describes that the function should switch to Safe State 3, so specifying only the target state and making no mention of source state; So, this incomplete requirement generates a state on the initial model that should be completed in the next step. Indeed, in the model construction phase a transition will be added from each state of the predefined list. 4) Check the plausibility. However, the transitions added are not all semantically correct. The model construction relies on the priority table to cut out the meaningfulness and undesirable transitions. For the example we are considering, the priority table reveals that the transitions from Safe State 2 to Safe State 3, and Safe State 4

	elsewhere in
	the STRComp document.
	2) Requirements selection. When the process of gathering
	knowledge has been completed, the set of collected
	states is submitted to an expert. He will select ac-
	cordingly the set of most relevant requirements for the
	analysis. Actually, the set of relevant requirements are
	all the requirements that can affect the list of provided
	states. For the analysed system, only 15 functional
	and safety requirements are kept. Let us consider the
	following two example requirements:
	REQ1: IF (ADAS ECU is in "APA ADAS ECU
	Safe State 1" AND Standstill is Activated AND
	pgen PLC ParkingSequenceRequest=vgen stop2) THEN
	(ADAS ECU shall continuously request Standstill AND set
	pgen PLC SequenceStatus to vgen PLC Move Ready).
	REQ2: IF (ADAS ECU is in "APA ADAS ECU Safe state 1"
	AND pgen PLC ParkingSequenceRequest=vgen PLCRamp)
	THEN ADAS ECU shall switch to Ramp state.

For our case study we identified six states: Safe State 1, Safe State 2, Safe State 3, Safe State 4, Idle State and Ramp State. The last two states they are also called out of maneuver state and in maneuver state resp. The two requirements concern the same function APA ADAS ECU within the function ADAS ECU. Both REQ1 and REQ2 refer to selected states Safe State 1 and Ramp state, but as one can notice REQ1 causes no change of states whereas REQ2 do. So, only REQ2 is kept for the rest of the analysis. 3)


Nicolas Dubos

Isabelle Le Viol

Alexandre Robert

Céline Teplitsky

Manon Ghislain

Olivier Dehorter

Romain Julliard

Pierre Yves Henry
email: henry@mnhn.fr

Disentangling the effects of spring anomalies in climate and net primary production on body size of temperate songbirds Running title: Climate and Body size of temperate songbirds

Body size is implicated in individual fitness and population dynamics. Mounting interest is being given to the effects of environmental change on body size, but the underlying mechanisms are poorly understood. We tested whether body size and body condition are related to ambient temperature (heat maintenance hypothesis), or/and explained by variations in primary production (food availability hypothesis) during the period of body growth in songbirds. We also explored whether annual population-level variations of mean body size are due to changes of juvenile growth and/or size-dependent mortality during the first year. For 41 species, from 257 sites across France, we tested for relationships between wing length (n = 107,193) or body condition (n = 82,022) and local anomalies in temperature, precipitation and net primary production (NDVI) during the breeding period, for juveniles and adults separately. Juvenile body size was best explained by primary production: wings were longer in years with locally high NDVI, but not shorter in years with low NDVI. Temperature showed a slightly positive effect. Body condition and adult wing length did not covary with any of the other tested variables. We found no evidence of climate-driven size-dependent mortality for the breeding season. In our temperate system, local climatic anomalies explained little of the body size variation. A large part of wing length variance was site-specific, suggesting that avian size was more dependent on local drivers than global ones. Net primary production influenced juvenile size the most through effects on body growth. We suggest that, during the breeding season in temperate systems, thermoregulatory mechanisms are less involved in juvenile growth than food assimilation.

INTRODUCTION

Body size is a key determinant of individual fitness and population dynamics, affecting reproductive performance and survival [START_REF] Ozgul | Coupled dynamics of body mass and population growth in response to environmental change[END_REF][START_REF] Yom-Tov | Recent spatial and temporal changes in body size of terrestrial vertebrates: probable causes and pitfalls[END_REF], Gardner et al. 2014b). With the increasing body of evidence of temporal changes in body size, there is an emerging interest in the impact of environmental factors, and in particular the influence of climate change (Gardner et al. 2011, Sheridan andBickford 2011).

Rapid body size changes in a population can arise from changes in growth conditions (Gardner et al. 2014b). Body size is largely determined by environmental conditions during the period of growth (i.e., from egg laying to the post-fledging period in birds; Yom-Tov and Geffen 2011). In cool climates, warmer temperatures can increase body growth as a result of a change in the cost of heat maintenance [START_REF] Kendeigh | Tolerance of Cold and Bergmann's Rule[END_REF]. Thus, in warmer years, juveniles can allocate more energy towards body growth, resulting in larger individuals [START_REF] Gillooly | Effects of size and temperature on metabolic rate[END_REF]).

In addition to external temperature, metabolic allocation to growth is also dependent on the amount of protein intake [START_REF] Dawson | The importance of microclimate variation in determining size, growth and survival of avian offspring : experimental evidence from a cavity nesting passerine[END_REF]. Between-year changes in body size may be driven by fluctuations in food availability [START_REF] Yom-Tov | Recent spatial and temporal changes in body size of terrestrial vertebrates: probable causes and pitfalls[END_REF]. Food limitation is particularly expected when the temperature is high and precipitation is low (or when precipitation is low, independently of temperature in arid systems, e.g., Gardner et al. 2014b), i.e. conditions that reduce primary production, and ultimately result in low prey availability for secondary consumers [START_REF] Aber | A generalized, lumped-parameter model of photosynthesis, evapotranspiration and net primary production in temperate and boreal forest ecosystems[END_REF]. Those climatic conditions are particularly constraining in arid ecosystems [START_REF] Holmgren | Extreme climatic events shape arid and semiarid ecosystems[END_REF]. For instance, in honeyeaters Ptilotula penicillatus, individuals are smaller in drier years, which are presumably the years with the lowest food availability (Gardner et al. 2014b). Hence, temperature can have both direct effects through thermoregulation, and indirect effects through its influence on food availability. The relative dependence of temporal change in body size on climate and primary production has not been assessed for any taxa in temperate systems yet (but see Gardner et al. 2014b for a case study in a semi-arid system).

In addition to body growth effects, body size composition in a given population can also be driven by size-dependent mortality (Gardner et al. 2014b). In temperate climates, temperatures rarely reach lethal or sub-lethal levels [START_REF] Tewksbury | Putting the Heat on Tropical Animals The Scale of Prediction[END_REF][START_REF] Khaliq | Global variation in thermal tolerances and vulnerability of endotherms to climate change[END_REF]). Hence, selective pressure related to thermoregulatory mechanisms is unlikely to drive significant body size change in a temperate system. Size-dependent mortality may also be related to food availability. Small individuals suffer higher mortality when food becomes scarce [START_REF] Ozgul | Coupled dynamics of body mass and population growth in response to environmental change[END_REF]. Studies that have investigated size-dependent mortality in terrestrial vertebrates were performed in arid, or semi-arid systems (McKechnie & Wolf 2010;du Plessis et al. 2012;Gardner et al. 2014b). In fact, little is known about the mechanisms underlying changes in body size induced by environmental changes in temperate climates. Specifically, there is a need to identify whether temperature and food availability operate through effects on juvenile growth or size-dependent mortality.

Differences in the relative importance of body growth effects and size-dependent mortality among species and populations could explain the lack of consensus about the direction of body size change. The few existing studies focusing on temporal change in body size showed contrasting trends between species (Gardner et al. 2014a[START_REF] Salewski | Morphological change to birds over 120 years is not explained by thermal adaptation to climate change[END_REF], but see [START_REF] Gardner | Shifting latitudinal clines in avian body size correlate with global warming in Australian passerines[END_REF][START_REF] Van Buskirk | Declining body sizes in North American birds associated with climate change[END_REF]) and/or between sites [START_REF] Meiri | Global change and carnivore body size: data are stasis[END_REF][START_REF] Collins | Heterogeneous changes in avian body size across and within species[END_REF]). This disagreement is reinforced by a probable publication bias towards cases exhibiting significant changes [START_REF] Meiri | Global change and carnivore body size: data are stasis[END_REF]). Among these studies, only a few tested the effect of interannual variation in temperature on body size, and even fewer assessed the relative importance of temperature versus net primary production (but see Gardner et al. 2014b). Moreover, some studies were based on Museum data (e.g., [START_REF] Salewski | Morphological change to birds over 120 years is not explained by thermal adaptation to climate change[END_REF] and may have been prone to temporal collection and curation biases. Some were based on one or two localities [START_REF] Van Buskirk | Declining body sizes in North American birds associated with climate change[END_REF][START_REF] Collins | Heterogeneous changes in avian body size across and within species[END_REF], which precludes drawing macroecological conclusions. For this reason, there is a need to assess the influence of temperature and food availability on body size at larger taxonomic and geographic scales.

When relying on wing length measurements, effects on growth or mortality apply to juveniles (through ontogeny), but also to adults as they undergo a moult of their flight feathers after breeding (Jenni and Winkler 1994). Therefore, the distribution of adult wing length in a population results from both the quality of feather growth during the period of moult and from size-dependent mortality. Here we assessed whether juvenile body growth and adult feather growth are influenced by interannual variation in local environmental conditions, which would account for the possibility that population-level body size variation is driven by changes in population composition through size-dependent mortality.

Body constitution is characterized by two independent dimensions: body size and body condition [START_REF] Canale | Differential plasticity of size and mass to environmental change in a hibernating mammal[END_REF]. Unfortunately, many studies infer temporal patterns of body size from body mass data, and typically confuse changes of these two dimensions. Body mass conveys information on both body size and individual body condition, i.e. the body fat and protein content [START_REF] Labocha | Morphometric indices of body condition in birds: A review[END_REF]. We are rarely able to access independent measurements of both dimensions. For birds, the best data available at a large scale are wing length, used as an index of body size [START_REF] Gosler | The field determination of body size and condition in passerines: a report to the British Ringing Committee[END_REF], and wing length-adjusted body mass, used as an index of body condition [START_REF] Labocha | Morphometric indices of body condition in birds: A review[END_REF]. Because the response of body size and body condition to climate variation can differ (e.g., [START_REF] Gardner | Individual and demographic consequences of reduced body condition following repeated exposure to high temperatures[END_REF], we assessed the relative importance of climatic conditions and primary production on both of these traits.

In the present study, for the 41 commonest songbird species of continental France, we explored between-year changes in body size and condition at the population-level by analysing the interannual variation in wing length and wing length-adjusted body mass, and their dependence on interannual fluctuations in climate and primary production during the breeding period, over the past 15 years. We tested the effect of local, interannual environmental fluctuations only during the breeding season, which corresponds to the period of body growth, as body size has been shown to be mainly driven by the conditions during this period (e.g., Gardner et al. 2014b), and because we had no information on bird locations and conditions experienced during the rest of their annual cycle. We analysed first-year birds and adults separately, to distinguish potential effects due to changes in body growth from changes in body size distributions in the population. As the distribution ranges of our study species exceed the bounds of our study area, this study represents a case in a temperate climate, regardless of potential 'edge effects' [START_REF] Jenni | Bird population trends are linearly affected by climate change along species thermal ranges[END_REF]. We addressed the following questions: (1) Are birds larger during or after warmer years, as expected under the heat maintenance hypothesis? (2) Are birds smaller or have lower body condition during or following poorly productive years? (3) Which driver (temperature or food availability) is the most important for body size, during the period of juvenile growth? (4) If they have any effect, do these drivers operate through effects on body growth or size-dependant mortality?

MATERIALS AND METHODS

Bird survey

We used individual records of juveniles for the 41 most captured songbird species extracted from the French Constant bird ringing Effort Sites (CES) scheme from 2000 to 2014 [START_REF] Robinson | Constant effort: Studying avian population processes using standardised ringing[END_REF], more information at http://crbpo.mnhn.fr, see 'STOC Capture').

Biometric data were collected by 382 volunteer bird ringers at 257 sites spread across France (Fig. S1 in Supporting Information), each site being monitored during 5.1 (± 3.9 SD) years [START_REF] Dehorter | Bird ringing and movement database for France[END_REF]. Overall, study sites were evenly distributed across the years (Fig. S2), and there was no bias in site-specific average temperature or latitude (Table S1).

Captures occurred 2.5 ± 1.3 SD times per breeding season per site, from May/early June (3 June ± 12.7 SD), until late June/early July (26 June ± 11.8 SD). For a given site, the number and date of capture sessions and the number and location of mist-nets were kept constant throughout the years. Each individual captured was individually marked, its species identified and aged (juvenile for birds born during the ongoing breeding season, or adult if born in previous years; [START_REF] Svensson | Identification Guide to European Passerines[END_REF]. Mist-netting of birds is most efficient in habitats with a low canopy (3-to-4 m high), so most CES sites are settled in shrublands, woodlands with dense understory, or reedbeds [START_REF] Eglington | Latitudinal gradients in the productivity of European migrant warblers have not shifted northwards during a period of climate change[END_REF].

Biometric data

We used wing length as a proxy for body size [START_REF] Gosler | The field determination of body size and condition in passerines: a report to the British Ringing Committee[END_REF], and body mass adjusted to wing length as a proxy for body condition (hereafter 'body condition index', [START_REF] Labocha | Morphometric indices of body condition in birds: A review[END_REF]. Wing length was measured with a butt-ended ruler as the length of the flattened wing chord from the carpal joint to the tip of the longest primary, at an average precision of 0.8 mm (± 0.1 across observers; i.e. most observers rounded the value to the nearest integer).

Body mass was measured with spring or electronic balances, to the nearest 0.5 g. We used fully grown juveniles measured during the breeding period (15 May -15 July, i.e. some days to weeks after fledging) and adults measured during their species-specific breeding period (Appendix S1 in supporting information). Measurements were performed by 2.4 (± 2.3 SD) trained measurers per site. Within-observer repeatability of wing length and of body mass measurements was 0.92 (± 0.11 SD) and 0.91 (± 0.10 SD), respectively (Appendix S2). For a given site, the observer is usually the same over the years. Differences in measurements between observers and sites were accounted for in our models (see statistical analysis below for model descriptions). For each species, we accounted for potential measurement errors by removing measurements that were beyond the upper and lower 0.5% limits of a Gaussian distribution fitted to the data. Only one measurement per individual was used in the analyses (with random choice for the 9.1% of data that were repeated measurements on the same individual). Hereafter, wing length measurements were from 40,071 juveniles and 67,122 adults, and body mass measurements from 30,783 juveniles and 51,239 adults.

Environmental variables

We used environmental variables that are already known to explain interannual fluctuations in avian body size, and that were relevant for the study period and area [START_REF] Keller | Effects of local environmental conditions on nestling growth in the great tit Parus major[END_REF]Van Noordwijk 1994, Gardner et al. 2014b). For each site and each year between 2000 and 2014, we computed average environmental conditions during the breeding period (1 st April-1 st July, i.e. the main period of reproduction, including parental allocation to egg laying, incubation, and nestling and post-fledging growth). These metrics were computed from daily records of mean temperature and total precipitation, and monthly records of NDVI images.

Populations are supposedly adapted to local thermal and trophic conditions [START_REF] Both | Climate change and population declines in a long-distance migratory bird[END_REF]).

To reveal the influence of between-year fluctuations in the local environment at each study site, independently from average local conditions, each raw variable was transformed into site-specific yearly anomalies by subtracting the average value of the corresponding site for the 2000-2014 period. This allowed us to control for confounding spatial effects (e.g., latitudinal size gradients). Site-specific yearly anomalies in mean temperature were expected to document fluctuations of the thermal constraint [START_REF] Kendeigh | Tolerance of Cold and Bergmann's Rule[END_REF], whereas total precipitation (alone or in synergy with mean temperature) or net primary production anomalies would record fluctuations of trophic resources (Yom-Tov and Geffen 2011) during the breeding period (Appendix S3). We used these three complementary trophic-related surrogates as they may characterise resource availability differently (e.g., Gardner et al. 2014b). We extracted daily mean temperature and daily total precipitation from the E-OBS meteorological dataset [START_REF] Haylock | A European daily high-resolution gridded data set of surface temperature and precipitation for[END_REF], with a 0.25° pixel resolution using climateExtract R package (https://github.com/RetoSchmucki). We used monthly averaged raster images of remotely sensed Normalized Different Vegetation Index (NDVI; Copernicus Service information 2016) as a proxy for net primary production [START_REF] Pettorelli | Using the satellite-derived NDVI to assess ecological responses to environmental change[END_REF]. The spectral reflectance covered by these images were the red and near infrared wavebands (0.61 to 0.68 µm and 0.78 to 0.89 µm, respectively), generally used for vegetal characterisation [START_REF] Kerr | From space to species: Ecological applications for remote sensing[END_REF]. We computed mean NDVI inside a buffer zone with a 5 km radius around each site from the raster images. The 5 km radius was chosen because it coincides with the level of precision of CES site geolocations. All variables were included in the analyses since their variations were largely uncorrelated (r < 0.35; Table S2). They were centred and scaled so that relative effect sizes could be compared between variables, regardless of their order of magnitude of mean and variance.

Adjustment variables

All statistical null models presented hereafter included the effects of (i) species (fixed term), and random variation between (ii) observers, (iii) sites and (iv) year. Wing length of juveniles increases slightly throughout the breeding season; this was accounted for by adding (v-a), a fixed effect of log-transformed date of the year, with a species interaction term. Body mass increases during the morning until it reaches a plateau at about noon; this was accounted for by adding (v-b) log-transformed time of the day (hour) as a fixed term, with a species interaction term (see justifications and model details in Appendix S4).

Statistical analysis

As population responses to environmental parameters can vary with their relative position in the species' thermal niche [START_REF] Jenni | Bird population trends are linearly affected by climate change along species thermal ranges[END_REF], we verified that the majority of the populations studied were not located at the edge of species distributions. Less than 4% of the data were located in the upper 10% of species' thermal ranges, and less than 0.4% were located in the lower 10% (See details in Appendix S5).

Interannual variation in environmental and biometric variables

The first step was to identify whether environmental conditions and biometric variables varied between years, and whether those variations were consistent.

For environmental variables, we examined temporal fluctuations using a spline function to estimate parsimonious, smoothed patterns of interannual variation. This was implemented with generalised additive mixed models (GAMMs, gamm4 R package; Wood and Scheipl 2014), with 'year' as a smoothed term and 'site' as a random effect.

For body constitution features, we assessed the proportion of interannual variation that was common to all sites (i.e. nation-wide variation) and/or common to all species (i.e. independent of species life history traits). This was assessed with a variance-partitioning method [START_REF] Grosbois | Modeling survival at multi-population scales using mark-recapture data[END_REF]. For wing length and body condition index, we used Linear Mixed Models (LMM; lme4 package version 1.1.7; [START_REF] Bates | lme4: Linear mixed-effects models using Eigen and S4. R package version 1[END_REF] to estimate betweenyear (byear), between-site (byear) and (bi) between-species i variances using random terms for the corresponding effects (Appendix S6). With these variance estimates, we identified the proportion of temporal variance in body size and condition (byear + byear,i + byear:site + byear:site i) that was common to all species and sites (byear), common to all sites but species-specific (byear + byear,i), or common to all species but site-specific (byear + byear:site).

Effect of environmental anomalies on wing length and body condition index

We assessed the dependence of wing length and body condition index on local environmental anomalies during the breeding period for both adults and juveniles. For adults, we also assessed the dependence on environmental anomalies during the post-breeding moult period of the preceding year, of species for which a moult period could be identified at the capture site (n = 13; Appendix S1; S7), and also during the breeding period of the year prior to capture (to account for potential lags between change in environmental conditions and the response of body size and condition; e.g., Gardner et al. 2014b). The post-breeding moult occurred at the breeding site for most of the study species [START_REF] Morrison | Flexibility in the timing of post-breeding moult in passerines in the UK[END_REF], so we tested the effect of environmental anomalies at the same location as during breeding periods. Models included mean temperature, total precipitation, the interaction between mean temperature and total precipitation, and mean NDVI for the breeding period (see equations in supporting information; Appendix S6).

Attempting to infer climate-driven size-dependent mortality in first-year birds

The risk of mortality before the first breeding attempt (as a yearling) may depend on the environmental conditions experienced throughout the year and varies among individuals, depending on their size: larger individuals would die more frequently in hotter years, and/or smaller individuals would die more often when born during poorly productive years. Ideally, the link between environmental conditions, size and survival should be investigated using mark-recapture models. However, the small sampling area (2-4 ha) of our study sites meant that our mark-recapture data were unsuitable for this analysis (high natal dispersal, transiency, and female breeding dispersal). Hence, size-dependent mortality was inferred from differences in average wing length (hereafter ΔWL) between juveniles in breeding season t, and yearling birds in breeding season t+1, after adjusting for feather abrasion and differences in capture probability between sexes (Appendix S9). ΔWL was computed when at least 10 measurements were available per site, year and species. This resulted in 138 data points for six species, from 46 sites (one data point representing one ΔWL for one species, at one site for two consecutive years) obtained from 2020 individual measurements in total. We then assessed whether ΔWL depended on local environmental anomalies of breeding season t using LMMs accounting for random variation between sites and years.

Model selection process

The dependence of wing length, body condition index and ΔWL on environmental anomalies were inferred using a multi-model selection based on Akaike's Information Criteria (AIC; [START_REF] Burnham | Model selection and multimodel inference: a practical information-theoretic approach[END_REF]; adjusted for small sample size for ΔWL, i.e. AICc). Models containing only effects with the highest statistical support have the lowest AIC values. The support for a model m relative to all other models considered was quantified by its AIC weight (wm). The relative importance of an explanatory variable i (Σwm,i) was quantified as the sum of wm of models containing this variable. To account for model selection uncertainty, model-averaged estimates of variable coefficients were computed using the 'best model set', defined as the set of models for which the cumulative sum of wm ≤ 95%. We also showed averaged estimates, once uninformative models (Arnold 2010) had been removed. Model averaging was performed only if the best model set did not include the linear effect of a given variable together with its quadratic effect, or an interaction [START_REF] Banner | Considerations for assessing model averaging of regression coefficients[END_REF]. Model selection and averaging were implemented using MuMIn R package version 1.9.13 [START_REF] Barton | MuMIn: Multi-model inference[END_REF]. The MuMIn function builds all possible combinations of the aforementioned effects (with each combination corresponding to a single model). Full models corresponded to null models (i.e. with adjustment variables), to which were added the additive fixed effects of environmental variables, an interaction between mean temperature and total precipitation anomalies (to allow for synergistic or antagonistic effects), and interactions between each environmental variable and species identity to allow species-specific responses. To accommodate for potential non-linear effects of environmental variables, a quadratic effect was also allowed for all environmental variables. Full model equations are described in Appendix S6. When a model included a quadratic term or an interaction term, the linear or additive effect, respectively, was systematically maintained in the model. For significant relationships, we verified the robustness of the linearity and quadratic assumptions using smoothed estimates obtained with a spline function of a GAMM version of the corresponding LMM.

Finally, we quantified the proportion of temporal variation in body size that was explained by each influential environmental variable (i.e. similar to a R²; [START_REF] Grosbois | Modeling survival at multi-population scales using mark-recapture data[END_REF]). This proportion was computed as the ratio of interannual variances estimated respectively with the model including the environmental variable and the null model. All analyses were performed using R version 3.3.0 (R core team 2016).

RESULTS

Temporal variations in environmental conditions and body constitution

Between 2000 and 2014, during the breeding periods, climate variables showed significant temporal variations (degrees of freedom > 8 for all variables, all p values for smooth terms < 0.001) but no temporal trends. NDVI showed a positive temporal trend over the study period (Fig. 1).

Interannual variations in wing length and body condition index were extremely heterogeneous across sites and species (Year:Site:Species interaction, Table 1). Juvenile wing length and body condition index varied between years (17% and 51% of the total variance estimated by random effects, respectively). Most interannual variation occurred at the site level and was species-specific (88% and 96%; Table 1). These patterns were similar for juveniles and adults.

Influence of local environmental anomalies on wing length and body condition index

We found strong statistical support for a quadratic effect of NDVI on juvenile wing length.

We did not find any effect of the tested environmental variable on adult wing length, nor on adult and juvenile body condition (Table 2; Appendices S7-S8).

Temperature

The effect of temperature anomalies on juvenile wing length received weak statistical support, an absence of effect being the most likely (Σwm= 0.34). In case of an effect, it would be linear and positive (GAMM: estimated degrees of freedom = 1, p < 0.0001; Fig. 2). Temperature anomaly explained 10% of the interannual variance in wing length common to all species (i.e. byear + byear:site; Fig. 3). Coefficients were robust to model averaging pitfalls, as temperature was only included as an additive, linear effect. After rescaling, wing length would increase by 0.09 mm ± 0.03 SE per degree Celsius (+0.31 mm ± 0.11 SE in years with the highest anomalies compared to years with average conditions).

There was no support for an effect of temperature in any of the other analyses (i.e., juvenile body condition index, adult wing length and adult body condition index; Table 3).

Net primary production

Juveniles were larger in years with positive NDVI anomalies (Σwi = 1, Table 3) with a quadratic relationship (see GAMM on Fig. 2; df = 2.34, p (smooth term) = 0.02). This effect of NDVI would be largely common to all species as models allowing for species-specific responses were not statistically supported (ΔAIC > 100). The squared-effect of NDVI explained 13% of the temporal variance in wing length common to all species (i.e. byear + byear:site; Fig. 3). Coefficients were robust to model averaging pitfalls, as NDVI was systematically included with its squared-effect in the best model set. These were similar for the average model and the best model (Table 3). Wing length increased by 1.1 mm ± 0.21 SE in years with the highest anomalies, compared to years with average conditions. Alternative variables related to primary production (precipitation, with interactive effects with temperature) did not receive any statistical support (Table 2). An interaction between temperature and precipitation was included in the best model set, but this effect was uninformative, as its statistical support was very weak (Σwm= 0.09) and it only captured 2% of the interannual variance common to all species. When removing this interaction from the best model set, our results remained qualitatively unchanged: the coefficient for temperature decreased from 0.019 to 0.014 (NDVI coefficients remained unchanged).

Climate-driven size-dependent mortality

Differences in average wing length between juveniles in year t and yearling birds in year t+1 (i.e. ΔWL) were not related to any of the environmental variables (Appendix S9). Hence, there is no indication of climate-driven size-dependent mortality.

DISCUSSION

In our dataset, body size fluctuations of songbirds were best explained by interannual variations in food availability (NDVI), although much of the temporal variance remained unexplained. We showed that the effect of temperature anomaly is fairly positive as expected under the heat maintenance hypothesis, but the effect is weak and relatively less important than net primary production in determining juvenile size. As we found no evidence for dependence of average adult body size, and yearling size-dependent mortality on NDVI or climatic anomalies, the relationship is probably driven by effects on body growth.

In contrast with former studies, we did not find a negative effect of temperature, a result most often found in species inhabiting arid regions (e.g., [START_REF] Yom-Tov | Global warming and body mass decline in Israeli passerine birds[END_REF]. In tropical and arid systems, species are more exposed to hyperthermia as they live closer to their upper thermal limit [START_REF] Tewksbury | Putting the Heat on Tropical Animals The Scale of Prediction[END_REF]. In a temperate climate such as that of France, the thermal envelope of species is much wider than at lower latitudes and climatic fluctuations rarely expose temperate animals to lethal or sub-lethal temperatures [START_REF] Addo-Bediako | Thermal tolerance, climatic variability and latitude[END_REF][START_REF] Deutsch | Impacts of climate warming on terrestrial ectotherms across latitude[END_REF][START_REF] Khaliq | Global variation in thermal tolerances and vulnerability of endotherms to climate change[END_REF]. If temperature was to have an effect on body size (Table 3), it would be positive, which is contradictory to the expectations of Bergmann's rule over time. A similar effect was found for adults in Australia during the breeding period, and this effect was attributed to size-dependent mortality (Gardner et al. 2014b). This may not be the case in France, as the 2003 heatwave increased avian productivity [START_REF] Julliard | Evidence for the impact of global warming on the long-term population dynamics of common birds[END_REF]).

Besides, avian mortality was not particularly high after the 2003 heatwave [START_REF] Ghislain | Réponses démographiques des oiseaux aux changements globaux dans l'espace : des tendances nationales à la gestion locale[END_REF].

Moreover, preliminary analyses did not provide support for a potential effect of the number of days with maximum temperature > 35°C (Gardner et al. 2014b) on juvenile body size (Appendix S10). These extremely hot conditions still remain rare in the French temperate climate. As temperatures rarely exceed near-lethal points in temperate systems, even at the hot edge of species distributions [START_REF] Khaliq | Global variation in thermal tolerances and vulnerability of endotherms to climate change[END_REF]), a positive effect of temperature would probably be mediated through body growth. This is consistent with the conclusions of a recent review [START_REF] Teplitsky | Climate warming and Bergmann's rule through time: is there any evidence?[END_REF] which suggested that body size decline may be caused mainly by changes in body growth as a result of non-adaptive plasticity. Under temperate climates, increasing temperature reduces the cost of juvenile heat maintenance [START_REF] Kendeigh | Tolerance of Cold and Bergmann's Rule[END_REF]. This results in higher metabolic allocation for growth, which could explain the positive effect of high temperatures on juvenile size [START_REF] Gillooly | Effects of size and temperature on metabolic rate[END_REF]. The direction of the relationship is thus consistent with previous studies performed in cool or cold regions (e.g., [START_REF] Collins | Heterogeneous changes in avian body size across and within species[END_REF][START_REF] Pérez | Nestling growth rates in relation to food abundance and weather in the Arctic[END_REF], or during a cool period of the year (e.g., Gardner et al. 2014b), and supports the fact that warming can benefit juvenile growth through changes in metabolic allocation in temperate climates. In France, species are probably more coldconstrained than hot-constrained, and hot anomalies are probably less detrimental for reproduction and juvenile growth [START_REF] Julliard | Evidence for the impact of global warming on the long-term population dynamics of common birds[END_REF]) than at lower latitudes. It is also possible that species sensitivity to climatic anomaly depends on the location of a given population with reference to the cold edge of its geographic distribution [START_REF] Jenni | Bird population trends are linearly affected by climate change along species thermal ranges[END_REF].

Warmer years would benefit species that are located near the cold edge, while it could harm those located at the hot edge. However, France is closer to the hot edge for most of our study species (i.e., typically Palearctic), and thus it is far from the cold edge of their distribution.

Hence, a positive effect of temperature is expected in temperate climates, even in the core of species distribution ranges.

In cool regions, warming induces an increase in primary production, provided that precipitation is not limiting, which contributes to improved food supply and results in larger individuals [START_REF] Searcy | Effects of early nutrition on growth rate and adult size in song sparrows Melospiza melodia[END_REF], so in contrast with arid systems where temperature and food availability are decoupled (Gardner et al. 2014b), it is hard to conclude whether the likely positive effect of temperature is related to reductions in the cost of body heat maintenance, or to increases in food availability in the system.

Higher NDVI values were associated with longer wings, presumably due to a positive relationship between vegetation production and invertebrate abundance [START_REF] Wimp | Increased primary production shifts the structure and composition of a terrestrial arthropod community[END_REF].

During juvenile growth, most songbirds are insectivorous, so higher invertebrate abundance improves juvenile protein intake, and may result in larger individuals (Lindström 1999). We did not detect any effect of NDVI on body condition, presumably because body mass is highly labile, varying over time-scales of hours, days and weeks [START_REF] Canale | Differential plasticity of size and mass to environmental change in a hibernating mammal[END_REF], and therefore within-breeding period fluctuations may not be correctly documented when analysing a single value per individual. In accordance with [START_REF] Keller | Effects of local environmental conditions on nestling growth in the great tit Parus major[END_REF], we found no effect of NDVI in years with low or average primary production. In years with poor resource availability, species may adjust the number of offspring to maintain a fair body size (i.e., size-number trade-off; [START_REF] Lack | Ecological Adaptations for Breeding in Birds[END_REF]). The absence of effect of negative NDVI anomalies might be explained by the adjustment of brood size according to climatic conditions and expected food abundance [START_REF] Parker | Optimal egg size and clutch sizeeffects of environment and maternal phenotype[END_REF]. In years with lower food availability, birds may produce less juveniles, thus enabling higher parental investment in individual offspring [START_REF] Smith | The trade-off between offspring number and quality in the great tit Parus major[END_REF], resulting in unchanged body size in years with scarce resources.

As expected, given the complexity of biological and ecological processes, the predictive power of our variables was relatively low. Primary production and climatic anomalies do influence juvenile size. However, the true proportion of body size variation that can be formally attributed to fluctuations in primary production and climate remains unknown, and is probably under-estimated given the simple, averaged variables used in the present study.

Body size is largely heritable, and most of the interannual variation should be captured by heritability (e.g., 75% in great tits; [START_REF] Garnett | Body size, its heritability and influence on juvenile survival among great tits[END_REF]. For this reason, only a small part of body size variance can be captured by environment. Yet, our variables were still able to capture some interannual variation in juvenile size, suggesting that their effect is not negligible. The influence of temperature fluctuations on wing length was weak, with a maximum increase estimated at 0.31 mm for the highest temperature anomaly (+3.6°C) compared to years with average conditions. High primary production had a stronger effect on wing length, resulting in an increase estimated at +1.1 mm in the year with the highest NDVI value. As climate warming is expected to increase the frequency of positive anomalies in primary production [START_REF] Melillo | Global climate change and terrestrial net primary production[END_REF], we predict that climate change will induce a body size increase in temperate songbirds. Temperature and precipitation showed important fluctuations over the study period. The extent of interannual fluctuations in our variables was similar to other studies (e.g., Gardner et al. 2014b, Teplitsky and[START_REF] Teplitsky | Climate warming and Bergmann's rule through time: is there any evidence?[END_REF], treating each year as independent category, with cold and hot years, and dry and wet years. This suggests that the relatively low explanatory power of the tested variables does not result from a lack of temporal variance. Our statistical models were based on extensive long-term data obtained from a large-scale monitoring program and took into account most of the possible sources of bias and noise: such as observer effect, spatial and temporal variability. With the high statistical support attributed to NDVI anomalies, we can safely conclude that net primary production is a better predictor of juvenile growth than temperature in our system.

Interannual variation was largely heterogeneous between sites, suggesting that variation in body size was largely related to local, rather than global factors. This is consistent with other studies performed on multiple species at different sites [START_REF] Meiri | Global change and carnivore body size: data are stasis[END_REF][START_REF] Collins | Heterogeneous changes in avian body size across and within species[END_REF]. This emphasises the necessity of considering alternative environmental variables or variables at finer spatial resolution to understand the very local determinants of size. The sensitivity to climate change is known to vary within species range [START_REF] Jenni | Bird population trends are linearly affected by climate change along species thermal ranges[END_REF][START_REF] Pearce-Higgins | Geographical variation in species' population responses to changes in temperature and precipitation[END_REF]. In our study area, mean temperature of the breeding season ranged from 6°C to 20°C, depending on the site. The effect of temperature anomaly may then differ between the hottest and the coldest sites. Yet, our results rely on the assumption of a uniform response to climatic variation across species ranges, thus ignoring possible spatial heterogeneity in the response to climatic anomalies. Variation in body size may also be driven by land use changes (e.g., [START_REF] Schmidt | Concomitant patterns in avian and mammalian body length changes in Denmark[END_REF]Jensen 2005, Desrochers 2010). Human activities such as agriculture, logging, garbage and gardening may affect food availability, and in turn could impact body size. However, volunteer bird ringers generally settle CES in places that are protected from deleterious anthropogenic activities. Food availability at the study sites should be closely related to primary production, with a limited confounding effect of human activities.

Another factor that may contribute to limitations of the explanatory power of NDVI is the temporal resolution (e.g., [START_REF] Kruuk | Contrasting effects of climate on juvenile body size in a Southern Hemisphere passerine bird[END_REF]. Bones and feathers are grown over a short period (ca. 1 month for juvenile growth, and ca. 2 weeks for primary moult), and the environmental conditions during this growing period are expected to be the most influential. This period varies between sites (latitude, altitude and habitat), species and individuals. Therefore, despite significant correlations with invertebrate abundance [START_REF] Wimp | Increased primary production shifts the structure and composition of a terrestrial arthropod community[END_REF], proxies such as NDVI, which are averaged for the whole breeding season, inevitably document only a limited fraction of the dependence of primary production on growth. The limited explanatory power of NDVI could also be explained by a temporal mismatch between species phenology and prey dynamics which are known to affect body size [START_REF] Husby | Testing mechanisms of Bergmann's rule: phenotypic decline but no genetic change in body size in three passerine bird populations[END_REF]). Nonetheless, we can hardly improve the temporal relevance of NDVI proxies, as this would require knowledge about each species, the sites and the yearly basis of the actual periods of growth.

As we did not perform a mark-recapture analysis, strictly speaking, we could not formally disentangle the contributions of tissue growth versus size-dependant mortality. Yet, only juvenile wing length depended on NDVI, and it did not explain differences in wing length between juvenile and 1 st -year birds in the next year. We can thus reasonably presume that the proximate mechanism is mainly a direct influence of food availability on nestling and post-fledgling growth. This direct effect on growth may also have been reinforced by sizedependent mortality in the nest, but we could not document mortality before the first capture.

After the juvenile stage, we did not obtain evidence of climate-driven size-dependant mortality or impaired feather growth during moulting. To our knowledge, only one study revealed a contribution of both growth and size-dependent mortality on population body size variation (Gardner et al. 2014b). A key challenge for future studies documenting the influence of environmental variability on body size, is to convincingly disentangle the respective contributions of growth and mortality. We did not consider winter conditions here because individual bird locations were unknown during the winter. However further work to investigate winter conditions and mortality is needed, since mortality in temperate climates is the most prevalent during the winter ( [START_REF] Balen | Population fluctuations of the great tit and feeding conditions in winter[END_REF], and winter conditions may drive sizedependent mortality [START_REF] Van Buskirk | Declining body sizes in North American birds associated with climate change[END_REF][START_REF] Brown | Fluctuating viability selection on morphology of cliff swallows is driven by climate[END_REF], Björklund et al. 2015[START_REF] Danner | A critical season approach to Allen's rule: bill size declines with winter temperature in a cold temperate environment (L Manne[END_REF].

CONCLUSIONS

This study shows that, in a temperate system, temporal variation in body size is better predicted by net primary production than climatic variables. Our results support the role of food availability during the breeding period on juvenile body size. In turn, body size variation was largely asynchronous between species and sites, and the predictive power of climate and net primary production was limited, emphasising the need to account for finer-grained local factors.

Table 1 Interannual variance partitioning of wing length and body condition index, for juveniles and adults of 41 songbird species. Variances were estimated from random effects of Linear Mixed Models. Models were adjusted for feather growth/abrasion, within-day body mass variation and adult sexual dimorphism. Percentages (in parenthesis) correspond to the proportion of the total temporal variance (i.e. the sum of all temporal variances).

Figure legends Fig. 1

 legends1 Figure legendsFig. 1 Partial residual climate anomalies during the breeding season in France for the 2000-

Fig. 2

 2 Fig. 2 Relationship between partial residual wing length and local anomalies of (a) NDVI and

Fig. 3

 3 Fig. 3 Variance partitioning Juvenile wing length for 41 songbird species in France between

Figure 2 .

 2 Figure 2.

Figure 3

 3 Figure 3

ACKNOWLEDGEMENTS

We thank the 382 bird ringers who voluntarily collected biometric data (Appendix S11) and made the present study possible, along with the continuous support of MNHN and CNRS. We are also grateful to Morgane Barbet-Massin, Léo Bacon, Eric Edeline, Vincent Devictor, Pierre Fiquet, Aliénor Jéliaszkov, Frédéric Jiguet, Boris Leroy, Karine Princé and Victor Saito for providing useful material and wise advice. N.D. was funded by LabEx BCDiv.; CT and AR were funded by ANR Bioadapt ANR-12-ADAP-0006. We acknowledge Copernicus for providing NDVI data, the E-OBS dataset from the EU-FP6 project ENSEMBLES (http://ensembles-eu.metoffice.com) and data providers in the ECAandD project (http://www.ecad.eu). We finally thank Liam Bailey, and six anonymous reviewers and Editors for their useful comments on previous versions of the manuscript. The final dataset used in this study is available online as supporting information.


Dr J Royes

Dr C Tribet

Pasteur

Dr Mrs O Ilioaia

Dr F Angius

B Miroux
email: bruno.miroux@ibpc.fr

Dr M Bally

Dr G V Dubacheva

Bacterial-based production of Thiol Clickable Genetically Encoded Lipid Nanovesicles

Jorge Royes, [a,b] Oana Ilioaia, [b] Quentin Lubart, [c] Federica Angius, [b] Galina V. Dubacheva, [c] Marta Bally, [c] Bruno Miroux,* [b] Christophe Tribet* [a] Abstract: Despite growing research efforts on the preparation of (bio)functional liposomes, synthetic capsules cannot reach the densities of protein loading and the control over peptide display that is achieved by natural vesicles. Here we present a microbial platform for high yield production of lipidic nanovesicles, with clickable thiol moieties in their outer corona. These nanovesicles show low size dispersity, are decorated with a dense, perfectly oriented and customizable corona of transmembrane polypeptides. In addition, this approach enables encapsulation of soluble proteins into the nanovesicles. Due to the mild preparation and loading conditions (absence of organic solvents, pH gradients or detergents) and their straightforward surface functionalization taking advantage of the diversity of commercially-available maleimide derivatives, engineering bacterial-based proteoliposomes are an attractive ecofriendly alternative that can outperform current liposome preparation methods.

Microbial biosynthesis offers significant advantages when in glassware synthesis is cost-intensive and/or provides poor yield and is particularly relevant for the preparation of fragile biomolecules. [1] Nowadays it is no longer limited to small molecules. With a growing demand for renewable sources and sustainable industrial processes, bioproduction can be an ecofriendly response for the preparation more advanced materials, such as polymers, [2] metal oxide nanoparticles [3,4] and supramolecular assemblies of proteins. [5,6] To our knowledge, bioproduction has hardly been considered as an alternative to in glassware synthesis for lipidic nanovesicles, though membranes and compartments are ubiquitous in cells. This contrasts with the variety of applications of synthetic liposomes, used as encapsulation systems in consumer products [7] and nano-medicine to protect fragile active compounds. [8] Synthetic liposomes are advantageously displaying a variety of molecules on their outer coronae, including polymers, peptides, fluorophores, etc, which is key to achieve optimal shelf-life, targeting, imaging functions. Accordingly, post-modifications in the corona of biological proteoliposomes (by conjugation with functional groups) is needed to approach the state of art of artificial lipid capsules. A recent renewal in this field came from nanovesicles secreted by cells in various physiological and pathological conditions. [9] Such extracellular vesicles hold promises on innovative therapies and diagnostics, but their recovery from cell extracts and batch-tobatch reproducibility are poor. [10,11] The interest for extracellular vesicles is an additional motivation for a versatile bioproduction of proteoliposomes. To obtain artificial proteoliposomes and/or mimics of membrane cells, the common (in vitro) route is to reconstitute purified membrane proteins into liposomes form detergent solutions. [12,13] The quality of insertion of membrane proteins (native folding and control over in/out orientation) and protein densities are however far from approaching the natural, cell-secreted vesicles. [11] Due their easy manipulation, bacteria provide the most affordable and handy source for nanovesicle bioproduction. [14] Diverting the lipid synthesis and protein insertion cellular machineries may unlock the access towards more controlled proteoliposomes, providing a versatile, ecofriendly source of functional lipid:protein capsules.

We present here a platform for the production of bacterialsourced proteoliposomes. We obtained well-defined, proteinloaded nanovesicles displaying a high density surface clickable groups (thiols from cysteine-containing tags). At variance with previously described approaches, it does not rely on outer membrane vesicle secretion suffering from low rate of vesiculation. [15] In contrast, we triggered the proliferation of internal membranes to increase the vesicle yield, [16,17] producing at the same time a dense array of engineered proteins on the nanovesicle surface. Internal membrane proliferation was achieved by overproducing a engineered membrane protein derived from AtpF, the b subunit of F0F1 ATP synthase (ATP-b), in C43(DE3) Ecoli strain. This strain specifically evolved and trained to produce and fold ATP-b at high levels. [17,18] Up to now, they were never considered as a tool to prepare and functionalize proteoliposomes. To this aim, we constructed a chimeric membrane protein, derived from AtpF, the b subunit of F0F1 ATP synthase (ATP-b). [17] We characterized the composition, lipid:protein ratio and explored the surface functionalization (via thiol-maleimide click chemistry) of the bioproduced vesicles (Figure 1).

Nanovesicle production, isolation and composition analysis

The membrane proliferation induced by ATP-b overproduction in C43 (DE3) E.coli strain was the chassis for the preparation of lipidic nanovesicles containing custom polypeptides (Figure 1a). Overproduction of the unmodified wild type ATP-b was compared to ATP-b functionalized with a FlAsH tag and a Tobacco Etch Virus protease (TEV) digestion sequences (quoted ATP-b/F/T) both inserted at position 35, [17] immediately after the membrane interaction domain (Figure 1c and Figure S1 in Supporting Info. N.B.: Attempts to produce ATP-b C-ter and N-ter fusion proteins failed to induce membrane proliferation and conduced to protein inclusion bodies). [19] FlAsH tag enables in situ specific fluorescence labeling, [20] allowing to identify and characterize the recombinant protein (Figure 1b). Additionally, its tetra-cysteine motif offers a mean for in vitro post-modifications via thiol-maleimide click chemistry. Production yield and qualitative compositions of lipid:protein membranes were determined on isolated vesicles obtained after bacterial lysis and differential centrifugation(Figure S2). Applying this straightforward and fully scalable procedure, we typically obtained 50 mg of proteins and 10 mg of lipids per liter of culture (Table S1). These amounts well above the yield of outer membrane vesicles production, which is generally around 0.5 mg of vesicles per liter of culture containing c.a. 100 µg of proteins. [21,22] Determination of lipid composition of purified membranes indicated a mixture of phosphatidylglycerol, cardiolipin and phosphatidylethanolamine in molar percentages of 18%, 12% and 70% respectively. Purified nanovesicles exhibited an unexpectedly low lipid:protein ratio (ca. 0.2 wt/wt), showing that the overproduced proteins are densely packed in the membranes. For comparison, common formulation methods of membrane protein reconstitution into artificial liposomes in a large excess of lipids, at lipid:protein ratios of about 10 to 800 wt/wt and require the use of detergents. [23] The present in vivo bacterial bioproduction system achieves higher protein density than any reported method with no need for detergent, neither organic solvents that may denature proteins.

The isolated vesicles contained predominantly the overproduced ATP-b scaffold (Figure 1b). The bands of ATP-b (17 kDa) and ATP-b/F/T (20 kDa) were identified in SDS-PAGE by both Western Blot anti His tag and FlAsH tag in gel fluorescence (Figure 1b, using FlAsH-EDT2 pro-fluorophore). In overloaded gels, Coomasie blue staining showed two additional major protein bands (Figure 1b). They were identified as OmpA (37 kDa) and OmpF (39 kDa) by MS. The three major protein components (ATP-b scaffold, OmpA and OmpF) were individually quantified using the fluorescence of tryptophan (Table S2). [24] Target peptides (ATP-b and ATP-b/F/T) represent ca. 20% wt/wt., whereas the sum of OmpA and OmpF accounts for ca. 40% wt/wt. of all proteins. Of note, this weight percentage corresponds to a clearly dominant molar fraction of the ATP-b peptide because of its low molecular weight (Mw of ATP-b is about a tenth of OmpA and OmpF ones). Based on an estimate of total mass of all proteins (determined by bicinchoninic acid assay), the sum of other proteins (each one being present in trace amount, see Figure 1b overloaded lanes) can be estimated to account for 40% wt/wt. The presence of OmpA and OmpF (major outer membrane proteins forming water-filled unspecific pores) in cytosolic membranes is not fully understood. It cannot be excluded that membrane proliferation may have facilitated unusual insertion of Omp or that contamination occurred during cell disruption, a point that will deserve future optimizations.

The orientation of the overproduced polypeptides in the lipidic bilayer was evaluated by assessing the accessibility of ATP-b Cterminal soluble domain to water soluble proteases. [17,[START_REF] Rigaud | Liposomes, Part B[END_REF] As trypsin do not spontaneously cross lipid bilayers, the solutionaccessible domains (outer corona of proteoliposomes) are degraded much faster than inner components. Both ATP-b and ATP-b/F/T were significantly degraded after 12h of incubation at 30 °C in the presence of trypsin (< 5% of initial intensity remains) (Figure 2). When the proteoliposomes were solubilized by addition of a surfactant, further progression of trypsinolysis was observed (complete degradation of ATP-b, ATP-b/F/T, OmpA and of proteins from 40 to 100 kDa, see Figure S4

Size distribution and thiol-Michael conjugation of the nanovesicles

The size and dispersity of vesicles were first measured by DLS to explore the impact of mechanical shear used during cell lysis (Figure S3 in SI). Increasing the extrusion pressure (from 0.25 to 2 kBar) gradually decreased the vesicle mean diameter from ~200 nm to 125 nm. Polydispersity index (PDI) was also slightly decreased by application of higher pressures. 2.0 kBar was selected as optimal lysis pressure to obtain vesicles with the narrowest possible size distribution (Table S3). Under these conditions, the diameter of nanovesicles ranged from 100 to 125 nm and PDI was < 0.2. This size dispersity is similar to what is typically achieved by a few extrusion cycles of artificially formulated liposomes (PDI of about 0.25 compared to 0.18 in our case). [START_REF] Frisken | [END_REF] Negative staining TEM images show circular objects in both ATP-b and ATP-b/F/T samples, confirming the presence of lipidic nanovesicles (Figure S6 in SI). Next, the population of nanovesicles was characterized by single particle counting using Nanoparticle Tracking Analysis (NTA) and the microfluidic device described by Friedrich et al. [27] The lipidic bilayer was nonspecifically marked using the lipophilic green fluorescent dye spDiO whereas the peptide was stained with Alexa-647, emitting at near infrared to avoid channel cross-talk. Alexa-647 was selectively conjugated to ATP-b/F/T FlAsH moieties using thiol-maleimide click chemistry (Figure S7). A single and slightly shifted band compared to unstained sample was observed in SDS-PAGE gels for ATPb/F/T after conjugation (Figure S8). These results suggest that i) the coupling reaction can be detected and ii) the vesicles could be selectively marked on ATP-b/F/T polypeptides. The fluorescence distribution curves (Figure 3a-b) allow to compare protein and lipids distributions in the vesicle population. Similarity between NTA and spDiO fluorescence intensity distributions (Figure 3a) suggests that the nonspecific spDiO dye is evenly distributed among the nanovesicle population and is a fair reporter of the lipid surface. The agreement between sizedistributions measured by DLS and NTA reinforces this statement (Figure S9). ATP-b/F/T fluorescence was homogeneously distributed among the nanovesicle population in a lognormal distribution, slightly broader than the lipid one (Figure 3b). A possible explanation of broadening is the presence of a minor amount of polypeptides with low or no lipids. To assess whether lipids and proteins are colocalized or form distinct assemblies we imaged by fluorescence microscopy vesicles immobilized into Koliphor ® hydrogels. Measurements of pixel intensity correlation [28] (Pearsons coefficient = 0.77±0.08 and Manders coefficients, Mgreen = 0.87±0.03 and Mred = 0.82±0.02) confirmed that most proteins are embedded in liposomes. The presence of a few particles stained in protein channel but unstained in lipid channel supports the aforementioned hypothesis of a minor fraction of proteinenriched clusters (Figure 3c).

In cellulo loading

Finally, we applied this bacterial expression platform to encapsulate soluble proteins (Figure 4). Usual in vitro encapsulations of water soluble molecules in the interior of liposomes require conditions that may denature proteins, such as freeze-drying, the presence of detergents, or interfaces with organic solvents. [29][30][31] For this reason, encapsulation of denaturation-prone enzymes or proteins remains challenging. The soluble fluorescent protein mCherry was co-expressed with ATP-b under the same T7 promotor using a single plasmid (Figure S1). Nanovesicles were isolated using the same procedure as described above. Fluorescence microscopy images of mCherry loaded nanovesicles immobilized in a Koliphor ® hydrogel showed fluorescent dots dispersed in a dark background (Figure 4b, red channel). Colocalization with nonspecific lipophilic spDiO staining revealed that most of the nanovesicles are loaded with mCherry (Figure 4b green channel and Figure S9). In addition, Western-Blot confirmed the presence of mCherry into nanovesicles (Figure S10c). Fluorescence spectra of free and nanovesicle entrapped mCherry are similar (Figure S11). Altogether, these results point to an encapsulation of native soluble mCherry inside the inner cavity of the nanovesicles. Protective properties of the lipidic nanovesicles against external environment were tested by protease digestion experiments. Fluorescence of soluble mCherry steadily decreased when non-encapsulated mCherry was mixed with proteases. In contrast, the fluorescence of mCherry-loaded nanovesicles was preserved under the same experimental conditions (Figure 4c).

Conclusions and perspectives

In summary, a robust and high yield method is proposed for the bioproduction of lipidic nanovesicles decorated with a genetically-engineered, custom polypeptide corona that can be easily (post)-modified. The density of the polypeptide corona largely exceeds (>20 fold) values achieved in artificial formulations of proteo-liposomes and enables to reach >95% orientation of C-ter end toward the outer solution. This C-ter domain remains accessible to post-modifications including in vitro sequence-specific cleavage and/or selective thiolmaleimide addition click conjugation. Furthermore, co-expressed soluble proteins can be encapsulated in the interior of the nanovesicles, as exemplified with mCherry.

The dense and almost fully oriented array of clickable peptides displayed on the top of vesicles could be particularly advantageous for applications requiring a high surface density functionalization (e.g. introduction of photodynamic therapy agents or macromolecules for enhanced shelf-life or vesicleinducing fusions). [32,33] As multiple tags can be introduced, they could allow multiple-functionalizations; for instance to combine specific targeting ligands with bio-repellent molecules. [START_REF] Curk | Proc. Natl. Acad. Sci[END_REF] Alternatively, proteoliposomes are specially adapted as signal amplification systems in biosensors. [START_REF] Liu | [END_REF] The present method provides a mild in vivo encapsulation of fluorescent proteins and possibilities for a dense surface display of fluorophores that should be useful for developing probes with high brightness. Finally, reconstitution of artificial cell mimics and liposome-based sensors is often struggling for co-encapsulating soluble enzymes and membrane proteins in capsules. [36][37][38] The present bioproduced lipidic vesicles offers a simple mean to combine all these constituents into well-defined assemblies.

Figure 1 .

 1 Figure 1. a) Schematic representation of the preparation of lipid:protein nanovesicles b) SDS-PAGE analysis of isolated vesicles and c) cartoon of ATP-b and ATP-b/F/T recombinant proteins. FlAsH tag contain four clickable thiol sites and TEV digestion site enables for post-cleavage of ATP-b cytosolic domain.

  in SI). Regular progression of the proteolysis was observed in a timecourse experiment (Figure S5 in SI). ATP-b/F/T were also treated with TEV protease. Specific cleavage of ATP-b/F/T soluble domain occurred at the expected position, generating a small polypeptide, which retained the fluorescence intensity of FlAsH (Figure 2). These results indicate that both the tags and cytosolic domain of ATP-b and ATP-b/F/T point towards the exterior of the nanovesicles. It should be noted that this almost 100% oriented display of custom, clickable polypeptides is hardly achieved in artificial formulation of liposomes.

Figure 2 .

 2 Figure 2. In vesicle enzymatic digestion of ATP-b and ATP-b/F/T using nonspecific (trypsin) or sequence specific (TEV) proteases.

Figure 3 .

 3 Figure 3. Size and fluorescence intensity distribution of ATP-b/F/T nanovesicles a) nonspecific lipid bilayer spDiO staining and b) Alexa-647 selective staining of ATP-b/F/T peptide. c) Fluorescence microscopy image of nanovesicles stained with spDiO (green), Alexa-647 (red) and merged images (yellow).

Figure 4 .

 4 Figure 4. a) Schematic representation of the in vivo production of mCherry loaded nanovesicles. b) Fluorescence microscopy image of mCherry loaded nanovesicles (mCherry in red channel and spDiO in green channel) and c)

Acknowledgements

This work was supported by CNRS, INSERM and Ecole Normale Supérieure. We acknowledge funding to support JR from the 'Initiative d'Excellence' program from the French State (Grant 'DYNAMO', ANR-11-LABEX-0011-01) and from the ANR GeneCap (ANR-17-CE09-0007).

We are grateful to Sandrine Masscheleyn for the training on proteomics (MS identification of OmpA and OmpF). We thank Dr. Audry Solgadi for lipid analysis by MS and the Région Ile de France for co-funding the SAMM MS Facility at IPSIT (Chatenay-Malabry, France)