Rodrigo Vargas
email: rvargas@udel.edu

Oliver Sonnentag

Gab Abramowitz

Arnaud Carrara

Jing Ming Chen

Philippe Ciais

Alexandra Correia

Trevor F Keenan

Hideki Kobayashi

Jean-Marc Ourcival

Ourcival

Dario Papale

David Pearson

Joao S Pereira

Shilong Piao

Serge Rambal

Dennis D Baldocchi

R V Wrote

Drought Influences the Accuracy of Simulated Ecosystem Fluxes: A Model-Data Meta-analysis for Mediterranean Oak Woodlands

Keywords: ecosystem models, eddy covariance, FLUXNET, forest survival, model evaluation, 75

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

INTRODUCTION 77

Water availability is the dominant control of roughly 40% of global terrestrial primary 78 productivity (Beer and others, 2010) with concurrent effects on ecosystem respiration 79 (Reichstein and others, 2002a;Schwalm and others, 2010). Recent increases in the global 80 incidence of drought (Dai and others, 2004), and decreased global evapotranspiration (Jung and 81 others, 2010) suggest that current terrestrial ecosystems are finding themselves in more water-82 limited environments [START_REF] Hartmann | Will a 385 million year-struggle for light become a struggle for water and for 577 carbon? -How trees may cope with more frequent climate change-type drought events[END_REF]. Such observations are in line with projections of climate 83 change, with higher temperatures and reduced soil moisture, leading to longer and stronger 84 seasonal droughts in the coming decades [START_REF] Giorgi | Climate change projections for the Mediterranean region[END_REF].

85

Process-oriented ecosystem models are critical tools for studying land-atmosphere fluxes 86 and for up-scaling this information to regional and global scales (Raupach and others, 2005).

87

Drought has a large impact on ecosystem state and function, and previous studies have 88 highlighted data-model disagreement for CO 2 fluxes in ecosystems suffering regular drought 89 such as those found in Mediterranean-type climate zones (e.g., Krinner and others, 2005; 90 Morales and others, 2005;Jung and others, 2007). The limited understanding of the hierarchy 91 and the temporal/spatial variation of drivers controlling the ecophysiology/biogeochemistry of 92 Mediterranean ecosystems hamper the accurate representation of H 2 O and CO 2 fluxes by models.

93 (Reichstein and others, 2002b;Keenan and others, 2010a;Misson and others, 2010). Thus, it is 94 critical to understand how drought stress influences these fluxes, and how process-oriented 95 ecosystem models predict them across different ecosystems.

96

Mediterranean ecosystems constitute a well-suited test-bed for the study and simulation 97 of drought effects, due to the combination of cool, wet winters with hot, dry summers (Joffre and 98 others, 2007). The asynchrony between the supply and demand for water causes annual potential 99 evaporation to exceed annual precipitation [START_REF] Baldocchi | What limits evaporation from Mediterranean oak woodlands -The supply of 511 moisture in the soil, physiological control by plants or the demand by the atmosphere?[END_REF], leading to seasonal 100 droughts which cause large declines in both assimilation and ecosystem respiration (Pereira and 101 others, 2007; Misson and others, 2010;Schwalm and others, 2010). Soil water availability is the 102 main factor limiting the metabolism of Mediterranean ecosystems, making them climate change 103 hotspots [START_REF] Giorgi | Climate change hot-spots[END_REF]. We performed simulations for all study sites considering differences in leaf habit (i.e., evergreen 150 vs. deciduous). We used five process-oriented ecosystem models (i.e., BEPS, Biome-BGC, 151 CABLE, ORCHIDEE, JULES; In this study we evaluated the performance of these models using standard plant 157 functional type parameterizations chosen in accordance with site-specific plant characteristics 158 (i.e., evergreen vs. deciduous) and individual model requirements (Table 1). The model runs 159 were without optimization to avoid enhancing of local agreement between observed-and model-160 calculated fluxes. For all models we used the same data for each site-specific characteristics 161 (e.g., soil texture, nitrogen content, canopy height, maximum leaf area index, soil depth; year and each site. This modeling exercise provided a total of 105 site-year simulations that were 167 analyzed for this study at the daily time-step (Supplementary Figures 1234).

168

Analysis of Ecosystem Drought Conditions 170

For subsequent model-data analyzes we assumed ecosystem drought conditions at a daily time 171 step to occur when the relative extractable soil water (REW) dropped below a threshold of 0.4 as 172 assumed in previous studies (Bernier and others, 2002). Briefly, REW was calculated from soil 173 water content at a daily time step:

174 (1)
175
where θ is soil water content reported in the La Thuile FLUXNET dataset, θ w is soil water 176 content at permanent wilting point, and θ fc represents water content at field capacity. The values 177 for θ w and θ fc were obtained from site-specific sand and clay contents (Table 1) and equations 178 provided by [START_REF] Saxton | Estimating Generalized Soil-Water 712 Characteristics from Texture[END_REF]. 191 wavelet and one of the most-used for geophysical applications (Torrence and Compo, 1998). To 193 analyze the data we first normalized the time series of observations and model outputs by: X' = 194 (x -mean(x))/std(x), where x represents the flux values of the time series of fluxes (e.g., GPP 195 either from observations (X' obs ), models (X' mod ), or residuals (measurements minus model 196 outputs)). Second, we calculated the global power spectrum using the continuous wavelet 197 transform of each normalized time series of fluxes or their residuals (Torrence and Compo, 198 1998). The length of each time series analyzed for observations (X' obs ), models (X' mod ), or 199 residuals was of 3 years. All time-series were analyzed using a daily time step for all 3 years of 200 measurements at each study site (i.e., 4 evergreen and 3 deciduous sites; 21 site-years), and for 201 each model simulation for each site (i.e., 5 models and 5 sites; 105 site-years). For the meta-analysis, we first calculated the response ratio as the model outputs (either 214 for the spectral power or the daily X' mod ) divided by the site observations (either for the spectral 215 power or the daily X' obs ) (Hedges and others, 1999). Second, the mean effect size was calculated 216 as the natural log of the response ratio for each flux (e.g., GPP) and each group (i.e., evergreen applied a random effects model using the MetaWin software (Rosenberg and others, 2000).

222

Furthermore, bias-corrected bootstrap 95% confidence intervals (CIs) were calculated for each 223 mean effect size (Rosenberg and others, 2000). Thus, if the 95% CIs of mean effect size did not 224 overlap with 1, then differences were considered to be significant at P < 0.05. We complemented 225 the meta-analysis by analyzing the root mean square error (RMSE), standard deviation, and 226 correlation coefficient (r values) between simulations and measurements for each individual 227 model and summarize them using Taylor diagrams [START_REF] Taylor | Summarizing multiple aspects of model performance in a single diagram[END_REF]. The three-year averages of GPP, Reco and ET for measurements, model outputs, and the 233 residuals (i.e., measurements minus model outputs) were summarized for evergreen sites (Figure 234 1). Subsequent analyzes in this study were done using all available data for measurements (n=12 235 site-years) and simulations (n=60 site-years; Supplementary Figure 1) for evergreen sites. These 236 sites had an average measured annual sum for GPP of 1518±241 (gC m 2 yr -1 ), Reco of 1260±191 237 (gC m 2 yr -1 ), and ET of 390±43 (mm yr -1 ; Figure 1). Using wavelet analysis we calculated the 238 global wavelet power spectra of GPP, Reco and ET from measurements, model outputs, and 239 residuals (3 years of data for each site and model analyzed). For a visual representation of the 240 differences in the power spectra only the measurements and residuals were plotted (Figure 2). In the time series. Overall, models underestimated the spectral power of GPP at time-periods of 2-254 10 days (95% CI of mean effect size <1; Figure 3a). For Reco, as a generality models 255 underestimated the spectral power at time-periods of 10-30 and 30-200 days (95% CI of mean 256 effect size <1; Figure 3b), but overestimated at 2-10 and 365 days (95% CI of mean effect size 257 >1; Figure 3b). In contrast, we found a general overestimation of the spectral power of ET at 258 time-periods of 2-10 and 10-30 days (95% CI of mean effect size >1; Figure 3c). An Subsequent analyzes were done using all available data for measurements (n=9 site-years) and 271 simulations (n=45 site-years; Supplementary Figure 2) for deciduous sites. Overall, deciduous 272 sites had an average measured annual sum for GPP of 1972±257 (gC m 2 yr -1 ), Reco of 1499±177

273

(gC m 2 yr -1 ), and ET of 500±89 (mm yr -1 ). Similar to evergreen sites, the global wavelet power 274 spectra of GPP, Reco and ET measurements showed strong power at 365-days and a lesser 275 seasonal power at ~200-days for deciduous sites (Figure 5). This means that all the time series 276 showed systematic oscillations that represent year-to-year variability (i.e., 365-days) and intra-277 annual variability that may be linked to plant phenology (i.e., ~200-days). Similarly to evergreen 278 site, residuals had a comparable spectral signature for all fluxes showing that model errors were 279 associated with these important fluctuations in the measured fluxes, but also showed that there 280 was large variability in the representation of these fluxes at these time-scales.

281

power of GPP at time-periods of 2-10, 30-200, and 365-days (95% CI of mean effect size <1; 283 Figure 3a). For Reco, models underestimated the spectral power at time-periods of 10-30 and 30-284 200 days (95% CI of mean effect size <1; Figure 3b), but overestimated at 365-days (95% CI of 285 mean effect size >1; Figure 3b). In contrast, we found a general overestimation of the spectral 286 power of ET at time-periods of 2-10 and 10-30 days (95% CI of mean effect size >1; Figure 3c), 287 but an overestimation at 365-days (95% CI of mean effect size >1; Figure 3c). 1) characterized drought conditions at these sites.

293

These results contrast with higher precipitation (2.8 mm day -1 ), lower photosynthetic photon flux 294 density (28.9 mol m -2 day -1 ) and cool temperatures (12.8 o C) characteristic of no drought 295 conditions (Table 3).

296

Individual model results for daily values of ecosystem fluxes at evergreen sites were 297 summarized in Supplementary Figure 5. In general, for GPP (Supplementary Figure 5a,b) and 298 Reco (Supplementary Figure 5d,e) the correlation coefficient decreased, and the RMSE and 299 standard deviation increased under drought conditions. In contrast, the correlation coefficient, 300 RMSE and standard deviation for daily values of ET tend was more consistent between drought 301 and no drought conditions (Supplementary Figure 5e,f). Importantly, we did not find that any 302 single model that consistently over-or under-represented ecosystem fluxes across all sites 303 (Supplementary Figures 1,2,5).

304

To synthetize and find general patterns from these results we applied meta-analysis 305 techniques. For evergreen sites the meta-analysis showed that models agreed with observations 306 in representing mean daily GPP under no drought and under drought conditions (i.e., 95% CI of 307 mean effect size not different from 1; Figure 6a). Models overestimated mean daily Reco under 308 no drought conditions but had even higher overestimation under drought conditions (i.e., 95% CI 309 of mean effect size >1; Figure 6a). Finally, models agreed with observations in representing 310 mean daily ET both under no drought and under drought conditions (i.e., 95% CI of mean effect 311 size not different from 1; Figure 6a).

312

Drought conditions (i.e., REW <0.4) occurred for 41% of the measured days at deciduous 313 sites. Drought conditions at these sites were characterized by low precipitation (0.7 mm day -1 ), 314 high photosynthetic photon flux density (46.7 mol m -2 day -1 ) and warm temperatures (20.8 o C; 315 Table 1). These results contrast with higher precipitation (3.3 mm day -1 ), lower photosynthetic 316 photon flux density (24.3 mol m -2 day -1 ), and cool temperatures (11.3 o C) found under no 317 drought conditions (Table 3). Noteworthy, these patterns were similar to those observed for 318 evergreen sties, but deciduous sites present lower sensitivity to drought.

319

Individual model results for daily values of ecosystem fluxes at deciduous sites were 320 summarized in Supplementary Figure 4. In general, for GPP (Supplementary Figure 6a,b) and 321 Reco (Supplementary Figure 6d,e) the correlation coefficient decreased, and the RMSE and 322 standard deviation increased under drought conditions. In contrast, the correlation coefficient, 323 RMSE and standard deviation of daily values of ET tend were more consistent between drought 324 and no drought conditions (Supplementary Figure 6e,f). We also did not find that any single

327

For deciduous sites the meta-analysis showed that models underestimated the 328 observations in representing mean daily GPP under drought conditions (i.e., 95% CI of mean 329 effect size <1; Figure 6b). In contrast, models overestimated mean daily Reco under drought 330 conditions (i.e., 95% CI of mean effect size >1; Figure 6b). Finally, models agreed with 331 observations in representing mean daily ET both under no drought and under drought conditions 332 (i.e., 95% CI of mean effect size not different from 1; Figure 6b). Our results show that overall there were data-model disagreements at multiple temporal scales 336 for GPP, Reco, and ET at both evergreen and deciduous Mediterranean oak woodland sites.

337

These results support the hypothesis that overall data-model disagreements are substantial at 338 temporal scales between weeks to months (i.e., 10-200 days), but important disagreements were 339 also found at shorter (i.e., 2-10 days) and larger (i.e., 365-days) time scales. Furthermore, the 340 meta-analysis showed that simulated Reco was overestimated under drought conditions for both 341 plant functional types. Previous studies suggest that with appropriate parameterization ecosystem 342 process-oriented models accurately represent ecosystem fluxes in Mediterranean ecosystems 343 (Keenan and others, 2009;Keenan and others, 2010b), but the necessity remains to identify the challenge for models to simulate phenological cycles especially for deciduous forests 357 (Richardson and others, 2012). Importantly, the spectral power of ET was systematically 358 overrepresented at both vegetation types between 2 and 30 days. This has implications for water 359 balance budgets, and the representation of the temporal variability of water-use-efficiency in 360 these and other semi-arid ecosystems. Thus, there is increasing interest in evaluating this variable 361 in process-oriented ecosystem models (Lawrence and others, 2007;Wang and others, 2011). [START_REF] Trumbore | Age of soil organic matter and soil respiration: Radiocarbon constraints on 739 belowground C dynamics[END_REF]Vargas and others, 2009); and e) potential 373 confounded parameterization during drought stress (Reichstein and others, 2003), or that 374 parameters are not representative of specific time scales (Mahecha and others, 2010).

375

Mediterranean oak woodlands have complex heterogeneous canopies that represent a 376 challenge to accurately estimate leaf area index (Ryu and others, 2010) and understory net 377 ecosystem exchange (Misson and others, 2007), especially for open canopy sites. These sites 378 may have an understory of grasses and forbs with different photosynthesis rates and phenology 379 (Joffre and others, 2007) that were not explicitly considered in the generic parameterizations of 380 the process-oriented models used in this study. For example, the understory could be a large Drought conditions were present for 30% of the measured days for evergreen sites and 41% of 396 the days for deciduous sites. Thus, the challenge imposed by drought conditions to simulating 397 ecosystem fluxes may affect the annual estimates of productivity and ecosystem carbon balance 398 at water limited ecosystems (Richardson and others, 2010;van der Molen and others, 2011) (Reichstein and others, 2003). Importantly, daily GPP was underestimated under 404 drought conditions at deciduous sites, which could be a result of the role of stomatal 405 conductance, maximum carboxylation rate (V cmax ) and maximum rate of electron transport (J max ) 406 influence on leaf scale photosynthesis under drought [START_REF] Xu | Seasonal trends in photosynthetic parameters and stomatal conductance 765 of blue oak (Quercus douglasii) under prolonged summer drought and high temperature[END_REF]Limousin and 407 others, 2010). This may be especially important when taken into account the REW and its 408 influence on photosynthesis rates at the canopy scale (Rambal and others, 2003) and model 409 parameters (Keenan and others, 2009;Keenan and others, 2010b).

410

Across sites and years we observed a slight reduction on the mean daily values of GPP 411 and Reco under drought conditions but these were associated with large variability (Table 3).

412

Previous studies have shown an important reduction in assimilation [START_REF] Xu | Seasonal trends in photosynthetic parameters and stomatal conductance 765 of blue oak (Quercus douglasii) under prolonged summer drought and high temperature[END_REF] 413 Limousin and others, 2010) and respiration (Reichstein and others, 2003) under drought.

414

However, there are also lag effects at the seasonal scale where assimilation and respiration 415 increase at the time when soil water content starts to decrease (Baldocchi and others, 2010; 416 Vargas and others, 2010a). Thus, lag effects between wet and dry seasons and their influence on 417 calculation of these fluxes divided by drought conditions (Table 3).

419

The potential access to ground water by the trees through deep roots is an important 420 adaptation that influence plant function (David and others, 2007;Miller and others, 2010) Reco may be challenging to represent under drought conditions especially after light rain pulse 430 events that could rapidly increase heterotrophic respiration. These light rain pulse events may not 431 be sufficient to reduce drought conditions as water may not permeate and increase soil water 432 content at deeper soil depths (Thomey and others, 2011). Furthermore, heterotrophic respiration 433 may be substantially increased after these light rain pulses due to organic matter breakdown by 434 photodegradation during the dry and hot summers (Brandt and others, 2009). Thus, (Sonnentag and others, 2008), and has been identified as a key 483 challenge for the next generation of global climate models (Lawrence and others, 2007).

104

  Oak woodlands are an important forest type in Mediterranean regions of the world (Joffre 105 and others, 2007). They range from open to dense forest that allows a well developed to sparse 106 understory comprising mostly of grasses, and are predominantly composed of evergreen and 107 deciduous trees (Reichstein and others, 2002b; Pereira and others, 2007; Baldocchi and others, 108 2010). Oak trees in Mediterranean woodlands have developed structural and functional attributes 109 that enable them to survive on a limited annual water supply by: a) having a generally low to 110 medium leaf area index, an erectophile leaf inclination angle distribution to reduce thermal road, 111 or low to medium tree cover to reduce total water loss at the landscape scale (Joffre and others, 112 2007; Ryu and others, 2010); b) down-regulating photosynthesis, respiration, and stomatal 113 conductance with decreasing soil water availability (Baldocchi and others, 2010); and c) the 114 potential of extending their root systems to tap ground water (David and others, 2007; Miller and 115 others, 2010). 116 In this study we use data from multiple Mediterranean Oak woodland sites to address the 117 following question: How well do process-oriented ecosystem models simulate the temporal 118 variability of CO 2 (i.e., gross primary production and ecosystem respiration) and H 2 O vapor 119 fluxes (i.e., evapotranspiration) in Mediterranean oak woodlands across years and under drought 120 or no-drought conditions? Based on previous studies we hypothesize that data-model 121 disagreement may be largest under drought stress (Reichstein and others, 2002b; Schwalm and 122 study including meteorological drivers were onsite measurements and available from the La 146 Thuile 2007 FLUXNET 2.0v dataset (www.fluxdata.org).

  spectral properties of the time series of daily GPP, Reco, and ET from 182 measurements and model outputs using the continuous wavelet transform (Torrence and Compo, 183 1998). Previous studies have reviewed in detail the concepts of wavelet analysis (Torrence and 184 Compo, 1998), and used this technique for model performance evaluation (Vargas and others, 185 2010b; Dietze and others, 2011; Wang and others, 2011). Briefly, this technique provides 186 information about the periodicities of the time series, and allows us to test for differences in the 187 spectral properties to identify the frequencies at which there are (if any) substantial data-model 188 disagreements. In other words, we are able to evaluate data-model agreement in the frequency 189 domain to know when (i.e., at which time scales) there is a data-model agreement or 190 disagreement.

  study we applied meta-analysis techniques(Hedges and others, 1999; 205Gurevitch and others, 2001) to determine the generalities of data-model disagreement between 206 observations of fluxes (i.e., control) and model outputs (i.e., treatments). The meta-analysis was 207 performed for: a) mean spectral power of GPP, Reco, and ET at different time-periods (i.e., 2-10, 208 10-30, 30-200 and 365 days); and b) mean daily flux estimates of GPP, Reco, and ET under no 209 drought and under drought conditions (i.e., REW<0.4). The different time-periods represent 210 weekly, monthly, seasonally, and yearly oscillations in the time series, respectively. With the 211 first analysis we were able to identify the time-periods for data model disagreement, and with the 212 second the effect of drought conditions for daily estimates of ecosystem fluxes.

  213

217

  and deciduous)(Hedges and others, 1999). A mean effect size of 1 indicates that the model 218 outputs were not different from observations, whereas a mean effect size >1 and <1 indicate that 219 the model outputs overestimated and underestimated the observations, respectively. To 220 determine if model outputs were different from site observations (i.e., mean effect size≠1) we 221

  241general, the global wavelet power spectra of GPP, Reco and ET from measurements showed high 242 power at 365-days and a lesser seasonal power at ~200-days. Importantly, a similar spectral 243 signature was found in the residuals of these time series showing that model errors were 244 associated with these important fluctuations in the measurements. A peak in the global power 245 spectrum means that there were strong systematic oscillations at that specific time-scale. For 246 example, a strong power at 365-days suggests systematic oscillations in the time series at the 247 annual scale, and from Figure 2 one can see that it represents the strongest oscillation but with 248 larger variability across the time series analyzed. 249 We were interested in testing how well the model outputs represented the temporal 250 oscillations of the measurements of GPP, Reco, and ET. Thus, we did a meta-analysis on the 251 spectral power of model measurements and model outputs at different time-periods (i.e., 2-10, 252 10-30, 30-200 and 365 days) representing weekly, monthly, seasonally, and yearly oscillations in

  253

  259 underestimation of the spectral power implies a reduction in the magnitude the oscillations of the 260 time series of a model output in comparison with the observations of each flux. Therefore, an 261 overestimation of the spectral power represents an increase in the magnitude of the oscillations 262 of the time series of model outputs in comparison with the measurements of each flux. In other 263 words, an overestimation of the spectral power of Reco at 365-days suggests a potential 264 overestimation of this variable when looking at the integrated response of daily simulations at the 265 averages of GPP, Reco and ET for measurements, model outputs, and the residuals 269 (i.e., measurements minus model outputs) were summarized for deciduous sites (Figure 4).

  270

  Means of CO 2 Fluxes and ET to Drought Conditions 290 Drought conditions (i.e., REW <0.4) were occurred for 30% of the measured days at evergreen 291 sites. Low precipitation (0.5 mm day -1 ), high photosynthetic photon flux density (47.6 mol m -2 292 day -1 ), and warm temperatures (21 o C; Table

325

  Figure 3, 4, 6).

  Figure2) were those between 30-200 and at 365-days for GPP, Reco and ET at both evergreen

362

  Our results show systematic errors in Reco for evergreen and deciduous sites that 363 underestimate the oscillations of the time series of Reco at 30-200 days and overestimate the 364 oscillations at 365-days. This is likely a result of an underestimation of the magnitude of Reco in 365 response of water pulses and synoptic events at scales between 30-200 days, but an 366 overestimation of the magnitude of Reco at 365-days. The large underestimation of the spectral 367 power for Reco may be a result of the challenge in representing: a) the overall response of Reco 368 to rainfall pulses in Mediterranean ecosystems (Baldocchi and others, 2006); b) the dependence 369 of heterotrophic and autotrophic respiration to seasonal variation; which are highly controlled by 370 water availability in Mediterranean ecosystems (Carbone and others, 2011); c) the control of 371 recent photosynthesis on Reco (Vargas and others, 2011); d) changes of the carbon residence 372 time in soils or in plants

  381contributor to GPP, Reco or ET along the year or may be dependent to seasonality(Misson and 382 others, 2007) driving interannual variability in ecosystem fluxes(Ma and others, 2007). In 383 addition, most land surface and process-oriented models assume that a vegetated canopy can be 384 abstracted as a turbid medium, but they fail to simulate radiation environments in heterogeneous 385 landscapes (Mariscal and others, 2004). 386 Our results support previous observations where ecosystem process-oriented models fail 387 to represent ecosystem fluxes at multiple time scales, and bring attention on the study of the 388 spectral characteristics of the time series as an alternative way for model evaluation (Mahecha 389 and others, 2010; Vargas and others, 2010b; Dietze and others, 2011; Wang and others, 2011).390We further hypothesize that model-data disagreement at short time scales may be related to 391 factors controlled by light, the intermediate timescales to drought conditions, and the annual 392 scales potentially to antecedent conditions and inter-annual weather variability.

435

  understanding pulse dynamics in these ecosystems has been an important topic to understand 436 shifts in contribution of ecosystem fluxes(Baldocchi and others, 2006; Thomey and others,

  and Future Considerations 440 Historically, model improvements have focused on how drought stress influence photosynthesis, 441 but fewer efforts have been made for the role of drought stress of Reco (van der Molen and 442 others, 2011) and ET(Ichii and others, 2009; Soylu and others, 2011; Thompson and others, 443 2011). Our results show that systematic overestimation in Reco during drought conditions at both 444 evergreen and deciduous Mediterranean oak woodlands, and support previous observations that 445 models need improvements on how drought conditions control Reco (Reichstein and others, 446 2002a). Finally, evergreen sites have lower photosynthetic capacity but longer growing seasons 447 than deciduous sites, which results at similar annual sums of ecosystem fluxes (Baldocchi and 448 others, 2010; and this study). Our results show that deciduous sites are more resilient to drought, 449 whereas evergreen sites are more sensitive and decrease substantially ecosystem fluxes under 450 drought. An incorrect representation of ecosystem fluxes may have larger implications under the 451 context of climate change where droughts may become more common (Kharin and others, 2007). 452 For example, if plants will encounter more water stress in future years and model outputs 453 overestimate Reco, then plants could be depleted faster from carbon (i.e., carbon starvation) and 454 die faster than in reality. This scenario could be enhanced by the potential underestimation of 455 GPP under drought conditions for deciduous sites limiting the carbon input into the system and 456 generating a positive feedback for carbon depletion and therefore plant mortality. Thus, the 457 omission or misrepresentation of how models response to water deficit will influence predictions 458 about how ecosystems may evolve, adapt and respond to past and future environmental changes. 459 Potential modeling experiments and advancements could focus on: a) model sensitivity 460 for water holding capacity, soil drainage and percolation that could affect model design and 461 parameterizations; b) the effect of changes in rooting depth and the ability of roots for extracting 462 water across the soil-water compartment (e.g., root weighted soil moisture, access of water by and plant functional types; c) consistently evaluate the autotrophic and heterotrophic contribution 465 of Reco across a wide range of ecosystems and weather conditions including different ranges of 466 drought stress; d) consider the uncertainty for model parameters, initialization and random and 467 systematic errors within measurements under the proposed experiments; e) studies on radiative 468 transfer models and the light absorbed by tree leaves and the understory layer (Kobayashi and 469 others, 2012); and e) hierarchy model comparisons to better identify differences in model-data 470 agreement based on a systematic evaluation of model structures (Vogel and others, 1995). 471 This study represents a first attempt to evaluate the overall performance of multiple 472 models for whole ecosystem carbon and water cycling across deciduous and evergreen 473 Mediterranean oak woodlands within FLUXNET. Our results summarize generalities on model-474 data disagreement but these vary for any specific site and any specific model that could be 475 improved by site-specific model parameterizations (Keenan and others, 2009; Keenan and 476 others, 2010b). Therefore, the generality of these multi-model meta-analysis results needs to be 477 tested at more study sites, across different years, and using a larger set of ecosystem process-478 oriented models. It is clear that longer time series are needed in order to account for model-data 479 agreements especially to account for inter-annual variability. Finally, the difficulty in 480 representing responses of ecosystem fluxes in drought prone ecosystems is exacerbated by the 481 biogeochemistry of a vertically stratified canopy with multiple layers of vegetation (e.g., grasses 482 and trees), soil hydrology

  Table 2) based on volunteer participation by attendants at a 152 FLUXNET workshop during 2009 (Asilomar, California, USA). The models used site 153 meteorological forcing input variables along with site-specific ancillary information (Table 1) to 154 calculate GPP, Reco and ET. All data were drawn from the La Thuile 2007 FLUXNET 2.0v

	156

155 dataset (www.fluxdata.org) and represent onsite measurements.

Table 1 )

 1 162and same forcing variables (e.g., precipitation, longwave/shortwave radiation, air temperature)

	163
	164	typically available at the global scale for general model simulations. Specific details about model
	165	architecture, parameters, and calculation of CO 2 and ET fluxes are found in Table 2. Model
	166	outputs and site measurements were averaged into daily sums for GPP, Reco and ET for each

  . Our 399 results support the hypothesis of overall data-model disagreement by systematically 400 overestimating daily Reco under drought stress at both plant functional types. The results support 401 that drought has a strong influence on simulating Reco, and therefore it is advisable that model 402 parameters should follow a seasonal course in response to this condition in Mediterranean

	403
	ecosystems

  . A 421 deeper rooting depth may result in access to more soil, and a potential access to ground water 422 that may have a direct effect on stomata conductance and an indirect effect on Reco. Furthermore, oak trees that have access to ground water may perform hydraulic lift and 424 redistribute water to shallow soil layers stimulating microbial activity and potentially 425 heterotrophic soil respiration. Recently there have been efforts improving how to simulate water 426 table depth and how roots utilize the water as the water table fluctuates in order to improve 427 evapotranspiration estimates (Ichii and others, 2009; Soylu and others, 2011; Thompson and 428 others, 2011), and could have implications for representation of other ecosystem fluxes. Finally,

	429

423

Mediterranean evergreen oaks: surviving the summer drought. Tree Physiology 27, 793-803.

Dietze, M.C., Vargas, R.,Richardson, A.D., Stoy, P.C., Barr, A.G., Anderson, R.S., Arain, M.A., Baker,

models: Bridging the gap between conflicting reports on the relative roles of stomatal, mesophyll

conductance and biochemical limitations to photosynthesis. Agr ForestMeteorol 150,[443][444][445][446][447][448][449][450][451][452][453]

Kharin, V.V., Zwiers, F.W., Zhang, X.B., Hegerl, G.C., 2007. Changes in temperature and precipitation

extremes in the IPCC ensemble of global coupled model simulations. J.Clim. 20, 1419Clim. 20, -1444. .

Kobayashi, H., Baldocchi, D.D., Ryu, Y., Chen, Q., Ma, S.Y., Osuna, J.L., Ustin, S.L., 2012. Modeling

energy and carbon fluxes in a heterogeneous oak woodland: A three-dimensional approach. Agr

Mahecha, M.D., Reichstein, M., Jung, M., Seneviratne, S.I., Zaehle, S.,Beer, C., Braakhekke, M.C.,

Carvalhais, N., Lange, H., Le Maire, G., Moors, E., 2010. Comparing observations and process-

based simulations of biosphere-atmosphere exchanges on multiple timescales. Journal of

Geophysical Research-Biogeosciences 115, G02003.

Mariscal, M.J., Martens, S.N., Ustin, S.L., Chen, J.Q., Weiss, S.B.,Roberts, D.A., 2004. Light-

Ecosystem Exchange measured with eddy covariance technique: algorithms and uncertainty

estimation. Biogeosciences 3, 571-583.

Paco, T.A., Rodrigues, A., 2007. Net ecosystem carbon exchange in three contrasting

Mediterranean ecosystems -the effect of drought. Biogeosciences 4, 791-802.

applications. 1. Hydrologic balance, canopy gas-exchange and primary production processes.

Ecol. Model. 42,

Ryu, Y.,Sonnentag, O., Nilson, T., Vargas, R., Kobayashi, H., Rebecca, W., Baldocchi, D.D., 2010. How

Water Resour. Res. 47.

Tirone, G., Dore, S., Matteucci, G., Greco, S., Valentini, R., 2003. Evergreen Mediterranean forests.

carbon and water fluxes, balances, ecological and ecophysiological determinants. In: Valentini, R.

tropical forest. NewPhytol. 182,[710][711][712][713][714][715][716][717][718]

Vogel, C.A., Baldocchi, D.D., Luhar, A.K., Rao, S., 1995. A comparison of a hierarchy of models for

determining energy balance components over vegetation canopies. JApMe 34, 2182-2196.

Wang, Y.P.,Kowalczyk, E., Leuning, R., Abramowitz, G., Raupach, M.R., Pak, B., van Gorsel, E.,

Luhar, A., 2011. Diagnosing errors in a land surface model (CABLE) in the time and frequency

domains. Journal of Geophysical Research-Biogeosciences 116, -.

ACKNOWLEDGEMENTS

486


Dai-Hua Tsai

Michael Riediker

Antoine Berchet

Fred Paccaud

Gerard Waeber

Peter Vollenweider

Murielle Bochud
email: murielle.bochud@chuv.ch

Effects of short-and long-term exposure to particulate matter on inflammatory marker levels in the general population

Keywords: Particulate matter, Inflammation, Short-term, Long-term, modeling, air pollution

Background:

The effect of particulate matter (PM) on health increase with exposure duration but the change from short to longer term is not well studied.

Objectives: We examined exposure to PM smaller 10 μm (PM10) from short to longer duration and their associations with levels of inflammatory markers in the populationbased CoLaus cohort in Lausanne, Switzerland.

Methods: Baseline and follow-up CoLaus data were used to study the associations between PM10 exposure and inflammatory markers, including the high-sensitivity Creactive protein (CRP), as well as interleukin 1-beta (IL-1β), interleukin 6 (IL-6), and tumornecrosis-factor alpha (TNF-α) using mixed models. Exposure was determined for each participant's home address from hourly air-quality simulations at 5-meter resolution. Shortterm exposure intervals were one day, one week, and one month prior to the hospital visit (blood withdrawal); long-term exposure intervals were three and six months prior to the visit.

Results:

In most time-windows, IL-6, IL-1β and TNF-α were positively associated with PM10. No significant associations were identified for CRP. Adjusted associations with longterm exposures were stronger and more significant than those for short-term exposures.

In stratified models, gender, age, smoking status and hypertension only led to small modifications in effect estimates, though a few of the estimates for IL-6 and TNF-α became non-significant.

Conclusion:

In this general adult cohort exposed to relatively low average PM10 levels, clear associations with markers of systemic inflammation were observed. Longer duration of elevated exposure was associated with an exacerbated inflammatory response. This may partially explain the elevated disease risk observed with chronic PM10 exposure. It also suggests that reducing prolonged episodes of high PM exposure may be a strategy to reduce inflammatory risk.

Introduction

Air pollution is a widely-acknowledged environmental risk factor globally, estimated to cause 3.3 million premature deaths per year worldwide [START_REF] Lelieveld | The contribution of outdoor air pollution sources to premature mortality on a global scale[END_REF]. Several studies have shown that exposure to particulate matter (PM) is associated with adverse cardiovascular effects [START_REF] Brook | Particulate matter air pollution and cardiovascular disease: An update to the scientific statement from the american heart association[END_REF]. Moreover, PM-induced inflammation has been postulated to be one of the mechanisms behind cardiovascular disease [START_REF] Brook | Particulate matter air pollution and cardiovascular disease: An update to the scientific statement from the american heart association[END_REF][START_REF] Pope | Health Effects of Fine Particulate Air Pollution: Lines that Connect[END_REF]. Inflammatory markers are known to be associated with increased risk of developing cardiovascular disease [START_REF] Kaptoge | Inflammatory cytokines and risk of coronary heart disease: new prospective study and updated meta-analysis[END_REF][START_REF] Kofler | Role of cytokines in cardiovascular diseases: a focus on endothelial responses to inflammation[END_REF]. Levels of inflammatory markers can also be mediated by lifestyle factors, such as smoking, or physical activity [START_REF] Bermudez | Interrelationships among circulating interleukin-6, c-reactive protein, and traditional cardiovascular risk factors in women[END_REF][START_REF] Colbert | Physical activity, exercise, and inflammatory markers in older adults: findings from the Health, Aging and Body Composition Study[END_REF]).

A series of epidemiologic studies have documented associations between short-term exposure to PM and inflammatory markers [START_REF] Hassanvand | Short-term effects of particle size fractions on circulating biomarkers of inflammation in a panel of elderly subjects and healthy young adults[END_REF][START_REF] Hertel | Influence of short-term exposure to ultrafine and fine particles on systemic inflammation[END_REF][START_REF] Steinvil | Short-term exposure to air pollution and inflammation-sensitive biomarkers[END_REF][START_REF] Tsai | Effects of particulate matter on inflammatory markers in the general adult population[END_REF], while long-term exposure also appears to play an important role. In one study in the United Kingdom (UK), an association was detected between long-term exposure to PM having aerodynamic diameter of 10 µm or smaller (PM10) and inflammatory markers in the general population [START_REF] Forbes | Chronic exposure to outdoor air pollution and markers of systemic inflammation[END_REF].

Meanwhile, investigators in the United States (US) reported one year exposure to PM2.5 to be associated with high-sensitivity CRP among midlife women [START_REF] Green | Long-and short-term exposure to air pollution and inflammatory/hemostatic markers in midlife women[END_REF][START_REF] Ostro | Chronic pm2.5 exposure and inflammation: Determining sensitive subgroups in mid-life women[END_REF]. In a German study, an association was identified between 365-day PM10, PM2.5 (PM < 2.5 µm) and number of particles, and CRP in the general population [START_REF] Hennig | Association between source-specific particulate matter air pollution and hs-crp: Local traffic and industrial emissions[END_REF]. Fewer studies have evaluated the consequences of PM exposure on inflammation markers, after both short-and long-term exposure, in the same population.

The upper thresholds for PM10 recommended by the World Health Organization are 50 !g/m 3 for a 24-hour mean and 20 !g/m 3 for an annual mean (WHO, 2018). In Switzerland, the air quality is mostly below these values, which provides an opportunity for the health effects of relatively-low concentrations of PM to be studied in the general population.

Assessing exposure is a critical step in epidemiological studies. Many studies focusing on short-term and acute effects have relied on data extracted from air monitoring stations, taking advantage of the fact that day-to-day variations in pollutant concentrations generally are much smaller than spatial variations within a given region. For long-term studies, however, data on spatial distribution is necessary, and this requires the use of models. Our objective was to examine for associations between both short-term and longterm exposure to PM10, and several circulating inflammatory markers: C-reactive protein (CRP), interleukin 1-beta (IL-1β), interleukin 6 (IL-6), and tumor-necrosis-factor alpha (TNF-α).

Materials and Methods

Study population

The CoLaus study is a large, longitudinal, population-based cohort survey that includes residents, ages 35 to 75 years old, of Lausanne, Switzerland, a city of 117,161 inhabitants (data for 2003 from city of Lausanne population registry). The study's focus has been on cardiovascular and metabolic risk factors, as described elsewhere [START_REF] Firmann | The Colaus study: A population-based study to investigate the epidemiology and genetic determinants of cardiovascular risk factors and metabolic syndrome[END_REF].

The Institutional Ethics Committee of the Faculty of Medicine of the University of Lausanne approved the CoLaus study. The informed written consent was obtained from all the participants.

A non-stratified random sample of 35% of the overall population was drawn. Recruitment was done between 2003 and 2006. The sample of 8,121 subjects who agreed to participate represented 41% of the initially-sampled population ("CoLaus 1"). Participants were asked to attend a specific outpatient clinic in the morning after an overnight fast. Between 2009 and2012, all CoLaus participants (n=6,183) were invited for a first followup ("CoLaus 2"), amongst whom 4,679 (75.7%) accepted. Both CoLaus 1 and CoLaus 2 included standardized questionnaires (as well as questions on medications), a physical examination, and blood exams. The final sample size were 3860 subjects, who had modeling estimates of exposure levels for both waves.

Health data

During both examinations, study participants reported, for their health examination, to the outpatient clinic at University Hospital of Lausanne (CHUV), always between 7 am and noon. Fasting blood samples were analyzed for the high-sensitively inflammatory marker C-reactive protein (CRP), as well as for interleukin 1-beta (IL-1β), interleukin 6 (IL-6), and tumor-necrosis-factor alpha (TNF-α). The laboratory equipment, analysis methods used, the lower detection limit and the method to substitute values below the lower detection limit were identical in CoLaus 1 and CoLaus 2, as described earlier [START_REF] Tsai | Effects of particulate matter on inflammatory markers in the general adult population[END_REF]. On follow-up (CoLaus 2), 1139, 197, and 21 values for IL-1β, IL-6 and TNF-α were below the lower detection limit.

Extreme outliers were defined as those that were 10-fold more than the corresponding measurement's 99 th percentile (p99) value; all such outliers were excluded from further analysis. The extreme outliers in CoLaus 1were described earlier [START_REF] Tsai | Effects of particulate matter on inflammatory markers in the general adult population[END_REF]. In CoLaus2, no extreme outliers were detected for IL-1β (p99=74.8 pg/mL). For IL-6, there was one extreme outlier at 2281 pg/mL (p99=222 pg/mL); and for TNF-α, two outliers, of 3242 and 3793 pg/mL (p99=93.8 pg/mL).

Exposure data

Assigning exposure to participants: Hourly PM10 concentrations at the home addresses of individuals living in Lausanne were simulated with the GRAMM/GRAL model (see below). We assigned to each participant home exposure data for one day, seven days (1-week), 30 days (1-month), 90 days (3-months) and 180 days (6-months) before their health assessment. The model was validated by comparing its output with measurements conducted at a central monitoring station in the city of Lausanne (see below). The daily modeled and measured concentrations were correlated with an R 2 coefficient of 0.83 [START_REF] Berchet | A costeffective method for simulating city-wide air flow and pollutant dispersion at building resolving scale[END_REF]. We calculated central monitoring site-time windows that corresponded to the modeled time windows for each of the participants, to understand gains in the population's exposure gradient, by using the model instead of central-site monitoring. As not all subjects lived within the GRAMM/GRAL model coverage area at both examinations, modeling data could only be assigned to 3860 participants who participated in both CoLaus 1 and CoLaus 2.

GRAMM/GRAL exposure model:

A coupled mesoscale-microscale model system -GRAMM/GRAL with 5-meter grid-resolution -was set up to simulate PM10 concentrations in Lausanne. Briefly, computational fluid dynamic simulations were computed to simulate 1008 different typical weather situations that influence the threedimensional wind and concentration fields in and around the city. To generate long-term series from pre-computed reference simulations, the catalogue of simulations was "matched to observations" (i.e., situations were selected from the catalogue of 1008 situations for which the simulated winds best matched observed winds at five different sites around Lausanne). Once the situations were identified, the simulated source-specific concentrations were scaled by characteristic time functions describing diurnal, day-ofweek, and seasonal variability of the emissions from different source types (traffic, industry, residential heating, etc.) and added together to obtain the total concentration field of PM10 [START_REF] Berchet | A costeffective method for simulating city-wide air flow and pollutant dispersion at building resolving scale[END_REF]. This modeling did not require any local pollution measurements.

Study area and central monitoring:

Lausanne covers an area of 70 x 15 km 2 and is built on the ridges and valleys of the northern shore of Lake Geneva. The highest point of the city is 600 m above lake level. One of the Swiss National air monitoring network (NABEL) stations is located in central Lausanne at César-Roux. A detailed description of the station can be found from NABEL website. We downloaded and averaged PM10 data into 24hour means (0:00-23:59). NABEL had missing data only during CoLaus 1. There were 15% of hourly PM10 data missing, most of the gaps were one or two hours of missing data. This affected 27 (0.7%) of the subjects' daily averages. In case of one or two hours of missing data, we replaced the missing hourly data by the weighted average of the two hours before and after the gap. Longer gaps were replaced by the average of the corresponding three daytime hours of the day before and after the gap. Outdoor air temperature and air pressure data was obtained from MeteoSwiss, Switzerland's Federal Office for Meteorology and Climatology, which provides quality-assured weather and climate services for the protection and benefit of Switzerland.

Analytic framework: to investigate the short-term effects of PM exposure, we were able to use data from both the central monitoring station and the models. For assessing effects of exposures longer than one month, however, only modeled data provided the necessary degree of variation related to the spatial variation between homes.

Statistical analysis

Stata 14.1 (StataCorp, College Station, TX) was used for all statistical analyses. Logtransformation was performed to approximate a normal distribution of skewed variables (CRP, IL-1ß, 1L-6, and TNF-α). Descriptive results were expressed as number of participants, percentage for categorical variables and as mean (95% Confidence Interval) or median (Interquartile range) for continuous variables. Bivariate analyses were performed using chi-square test for qualitative variables and Student's t-test for quantitative variables between the two waves.

In a prior study reported by Tsai et al., a positive association between PM10 exposure and pulse pressure and systolic blood pressure was evident only when the outdoor temperature was above 5°C [START_REF] Tsai | Shortterm effects of particulate matter on pulse pressure in two general population studies[END_REF]. Therefore, we created a covariate 'coldday ratio', corresponding to the percentage of days within a given exposure period when the outdoor daily average temperature was below 5°C.

Linear mixed-effect regression models were used to examine the associations between the inflammatory markers and PM10. The outcome variables were CRP, IL-1ß, 1L-6, and TNF-α, and the exposure variable was the PM10 level from modeling and monitoring stations. Each regression model included fixed effects for age, sex, BMI, smoking status, diabetes, hypertension, cold-day ratio, outdoor air pressure, and CoLaus time period. The participants were fitted as random effects in the mixed-effects models.

Estimated effects for both long-and short-term exposure were presented as unit change per 1 µg/m 3 increase in PM10. Short-term exposure interval means were 1-day, 1-week, and 1-month averages of PM10 level prior to the hospital visit (blood withdraw); while longterm exposure interval means were calculated for three and six months prior to the visit.

Results

Characteristics of the CoLaus 1 and CoLaus 2 participants are summarized in Table 1.

Only participants with health and environmental data available for both data-collection waves were included in this analysis. Body mass index was lower at baseline than at follow-up. Smoking rate dropped at follow-up, while diabetes and hypertension prevalence rates increased. Except for IL-1ß, the inflammatory marker levels increased from baseline to follow-up. Table 1 also shows the distribution of PM10 and temperature. In general, concentrations in CoLaus 1 were higher than in CoLaus 2. The distributions of inflammatory markers were similar in the two waves (Supplementary Figure S1).

For short-term exposure, PM10 levels derived from the NABEL station, averaged over the 1-day, 1-week, and 1-month intervals, were positively (directly) associated with IL-1β, IL-6 and TNF-α (Figure 1). Similar results were obtained when using short-term exposure data extracted from the GRAMM-GRAL models. Long-term modeled PM10 exposure (3month and 6-month) exerted more profound effects on these inflammatory markers. For long-term analyses (3-month and 6-month), only modeled data were used, due to the small variations in long-term pollutant levels at the single, central site. No association was found for CRP with either short or long-term exposure. The PM10 levels also were associated with IL-1β, IL-6, and TNF-α after dichotomizing these cytokines into high (higher quartile) versus low (quartiles 1 to 3) levels, which suggests that the associations between PM10 with the continuous cytokine were very robust to any distributional assumption, and to potential outlier values (Supplementary Table S1).

Associations between 6-month PM10 exposure levels and inflammatory markers are provided for different subject subgroups in Table 2. Six-month exposure to PM10 was significantly associated with IL-1β and IL-6 in men, in younger participants, and in normotensive participants; but with IL-1β and TNF-α in women and older participants.

Associations with IL-1β, IL-6 and TNF-α were more significant in non-smokers. However, the interactions between PM10 exposure and the selected subgroups were not statistically significant for the various inflammatory cytokines.

In general, PM10 concentrations decreased significantly from 2003 to 2012, covering both examinations (Supplementary Figure S2). The variability in PM10 levels was greater for modeled than for NABEL data (Supplementary Figures S3 andS4).

Discussion

In the current study, short-term and long-term exposures to PM10 were associated with elevated levels of IL-1β, IL-6 and TNF-α, but not of CRP. The short-term effects obtained with modeled air pollution exposures were consistent with our previous findings that used only NABEL data [START_REF] Tsai | Effects of particulate matter on inflammatory markers in the general adult population[END_REF]). The long-term effects found in the current analysis also were much stronger than the short-term effects.

To date, only a few studies have assessed the association between long-term exposure to PM and CRP, IL-1β, IL-6 and TNF-α in general population. In a Swedish study, IL-6, TNF-a, and CRP were not significantly associated with one-year, five-year or 30-year exposure to PM10 [START_REF] Panasevich | Associations of long-and short-term air pollution exposure with markers of inflammation and coagulation in a population sample[END_REF]). In the US MESA study, long-term exposure (one-year) to PM2.5 was associated positively with IL-6, but not with CRP [START_REF] Hajat | Longterm exposure to air pollution and markers of inflammation, coagulation, and endothelial activation: A repeat-measures analysis in the multi-ethnic study of atherosclerosis (MESA)[END_REF]. Meanwhile, in the US SWAN study, one-year exposure to PM2.5 was associated with increased CRP in women, but IL-1β, IL-6 and TNF-α were not included in the investigators' analysis [START_REF] Green | Long-and short-term exposure to air pollution and inflammatory/hemostatic markers in midlife women[END_REF]).

The possible pathophysiological mechanism underlying the link between PM and systemic inflammation could be the increased production of the pro-inflammatory cytokines such as TNF-a and IL-6 by alveolar macrophages in response to PM [START_REF] Panasevich | Associations of long-and short-term air pollution exposure with markers of inflammation and coagulation in a population sample[END_REF]van Eeden et al. 2001). Prolonged exposure to PM will result in more oxidative stress and prolonged activation of cell-and cytokine-mediated pathways including apoptosis, which can progress into a chronic inflammatory response with insufficient clearance of apoptotic cells [START_REF] Grabiec | The role of airway macrophages in apoptotic cell clearance following acute and chronic lung inflammation[END_REF][START_REF] Robb | Key mechanisms governing resolution of lung inflammation[END_REF].

Interleukin 6 is an acute phase reactant released by T-cells and macrophages. In vitro studies have revealed increased production of the proinflammatory cytokines, TNF-α and IL-6, by alveolar macrophages in response to PM [START_REF] Becker | Regulation of cytokine production in human alveolar macrophages and airway epithelial cells in response to ambient air pollution particles: Further mechanistic studies[END_REF][START_REF] Imrich | Alveolar macrophage cytokine response to air pollution particles: Oxidant mechanisms[END_REF]van Eeden et al. 2001). The stronger effects that we observed with prolonged exposure levels suggest that there is no adaptation to such exposure over time; instead, there appears to be propagation of the inflammatory process by air pollutants. This may explain the increased rates of chronic inflammatory diseases associated with living in polluted areas [START_REF] Chen | Metabolic syndrome and inflammatory responses to longterm particulate air pollutants[END_REF]. In a cross-sectional study conducted in a representative sample of elderly Taiwan residents, long-term PM10, PM2.5 and NO2 concentrations were all linked to IL-6 levels [START_REF] Chuang | Long-term air pollution exposure and risk factors for cardiovascular diseases among the elderly in Taiwan[END_REF].

We identified no association between exposure to PM10 and circulating CRP levels, with either long-term or short-term exposure, which is consistent with most prior populationbased studies. In the US MESA study, circulating CRP levels were weakly, positively associated with exposure to 30-day and 60-day mean PM2.5 exposures, but not to shorter term (one-day, two-day) exposures (Diez [START_REF] Roux | Recent exposure to particulate matter and c-reactive protein concentration in the multi-ethnic study of atherosclerosis[END_REF]. No association between shortterm exposure to PM10 and CRP was found in the Tel-Aviv Sourasky Medical Center Inflammation Survey [START_REF] Steinvil | Short-term exposure to air pollution and inflammation-sensitive biomarkers[END_REF]. Similarly, no association between long-term exposure to PM10 and circulating CRP levels was detected in 17'000 adults who participated in the Health Survey for England [START_REF] Forbes | Chronic exposure to outdoor air pollution and markers of systemic inflammation[END_REF]. The absence of an association between PM10 and CRP does not imply the absence of an inflammatory response to PM. Rather, CRP, as a ubiquitous inflammatory marker, just might not be ideal to detect the low-grade, chronic, respiratory inflammation that can result from relatively-low levels of air pollution. One of the results identified in the European Study of Cohorts for Air Pollution Effects (ESCAPE) was that elevated PM10 levels were not associated with elevated CRP; however, living close to busy traffic was [START_REF] Lanki | Air pollution from road traffic and systemic inflammation in adults: a cross-sectional analysis in the European ESCAPE project[END_REF].

Assessment of subgroups showed no significant interaction terms and the effect estimates were also in a similar range for most explored strata. Estimates across the various strata lied mostly well within the confidence intervals. A significant positive association of PM10 exposure with TNF-a levels was observed in women, but not in men, whereas similar significant positive associations with IL-1b and IL-6 were observed in both men and women. Participants younger than 55 years tended to have a weaker association with IL-1β and a stronger one with IL-6 than older participants. A significant positive association of PM10 exposure with TNF-a levels was observed in non-smokers, but not in smokers, possibly the result of the pro-inflammatory action of smoking [START_REF] Mcevoy | Relationship of Cigarette Smoking With In ammation and Subclinical Vascular Disease The Multi-Ethnic Study of Atherosclerosis[END_REF].

Normotensive participants showed a stronger positive association of CRP with exposure to PM10 than hypertensive participants, which could reflect the anti-inflammatory property of statins used by many hypertensive people [START_REF] Miyata | Statins reduce ambient particulate matter-induced lung inflammation by promoting the clearance of particulate matter <10 μm from lung tissues[END_REF].

One important strength of the current study is that it is representative of the general adult population of a city. Many previously-published studies focused on specific population groups, like women [START_REF] Green | Long-and short-term exposure to air pollution and inflammatory/hemostatic markers in midlife women[END_REF][START_REF] Ostro | Chronic pm2.5 exposure and inflammation: Determining sensitive subgroups in mid-life women[END_REF], the elderly [START_REF] Rioux | Residential traffic exposure, pulse pressure, and c-reactive protein: Consistency and contrast among exposure characterization methods[END_REF] or patients. Our study also included several inflammatory cytokines, while previous general population studies only included CRP (Diez [START_REF] Roux | Recent exposure to particulate matter and c-reactive protein concentration in the multi-ethnic study of atherosclerosis[END_REF][START_REF] Hertel | Influence of short-term exposure to ultrafine and fine particles on systemic inflammation[END_REF][START_REF] Hoffmann | Chronic residential exposure to particulate matter air pollution and systemic inflammatory markers[END_REF][START_REF] Rioux | Residential traffic exposure, pulse pressure, and c-reactive protein: Consistency and contrast among exposure characterization methods[END_REF][START_REF] Steinvil | Short-term exposure to air pollution and inflammation-sensitive biomarkers[END_REF]. Our longitudinal design included repeated measurements at baseline and 5-year follow-up, which reinforces the robustness of the results. A final important strength of the current study are the accurate estimates of PM10 exposure [START_REF] Berchet | A costeffective method for simulating city-wide air flow and pollutant dispersion at building resolving scale[END_REF]), since we estimated exposures for 5meter grid-lengths at hourly intervals, using air dispersion models, also taking into account both meteorological data and emission inventories. Some limitations exist. First, for particles, only modeled and measured data on PM10 was available in Lausanne, thus this study does not provide granularity on how much of the inflammatory response may be attributable to the size fractions PM2.5 and ultrafine particles. Second, gaseous pollutants such as nitrogen dioxide and ozone are also known to influence inflammatory processes. Third, Lausanne is a comparably clean city, thus caution should be applied when translating these findings to other, more polluted regions of the world. Fourth, exposure estimates were based on home addresses. However, we did not have the information on when the participants moved to their current addresses.

While they are usually fairly well correlated with personal exposure [START_REF] Culyba | Comparing residence-based to actual path-based methods for defining adolescents' environmental exposures using granular spatial data[END_REF][START_REF] Shafran-Nathan | Exposure estimation errors to nitrogen oxides on a population scale due to daytime activity away from home[END_REF], they are associated with an error, which increases the uncertainty of the effect estimates. Finally, our analysis focused on PM10, but further work should extend the analysis to other air pollutants, such as NO2 or O3.

Conclusions

Our results indicate that long-term exposure to high PM10 levels is associated with elevated markers of systemic inflammation in the general adult population, when PM levels generally are low. Prolonged elevated exposure to PM10 seems to exacerbate systemic inflammatory response, which may contribute to the elevated disease risk generally observed with chronic PM10 exposure. A clear increase in effect estimates started after one month of high exposure, which suggests that reducing elevated exposures to less than one month may be helpful to reduce negative long-term effects.

GlaxoSmithKline, the Faculty of Biology and Medicine of Lausanne, and the Swiss National Science Foundation (grants 33CSCO-122661, 33CS30-139468 and 33CS30-148401). We also thank the funding from Nano-Tera.ch to support Dai-Hua Tsai for doing this work (grant RTD2013: 20NA20_150956).

van Eeden SF, Tan WC, Suwa T, Mukae H, Terashima T, Fujii T, et al. 2001

Figure1.

  Figure1. Association between exposure to PM10 on short-term and long-term averages from NABEL monitoring stations and Model data. Data are regression coefficients and 95% confidence intervals of mixed models. Cytokine units are µg/mL for CRP, pg/mL for the other cytokines. All estimates are adjusted for age, sex, BMI, smoking status, diabetes, hypertension, cold-day ratio, outdoor air pressure, and CoLaus time period.

Figure S2 .

 S2 Figure S2. Daily concentrations of PM10 at César-Roux Station from 2003 to 2012 (CoLaus1 : 2003-2006, CoLaus2 : 2009-2012).

Figure S4 .

 S4 Figure S4. Box plot of different averaging PM10 concentrations of César Roux station and modeled data.

Table 1 . Characteristics of participants, that had participated in both waves and that had the addresses within the pollution map (N=3860). The data indicates categorical data (% of total), and continuous data (mean, 95% CI or median, 25%-75%).

 1 . Cytokines involved in the systemic inflammatory response induced by exposure to particulate matter air pollutants (PM10). American Journal of Respiratory and Critical Care Medicine

	164:826-830.
	World Health Organization. 2018. Website: http://www.who.int/news-room/fact-
	sheets/detail/ambient-(outdoor)-air-quality-and-health

*** p-value <0.001. Bivariate analyses were performed using chi-square test for qualitative variables and Student's t-test for quantitative variables between the two waves.

Acknowledgments

We express our gratitude to participants in the CoLaus studies, and to the healthcare professionals and collaborators who have contributed to data collection. We thank MeteoSwiss for providing meteorological data, and the National Air Pollution Monitoring Network (NABEL) for providing air-pollution data in Lausanne. We thank the Direction Générale du Canton de Vaud for providing in-situ data and a pollution inventory for the city of Lausanne. The CoLaus study was and is supported by research grants from

a Adjusted factors: gender, age, BMI, smoking status, diabetes, hypertension, outdoor air temperature, outdoor air pressure and season. b Cytokine levels equal to or above the 75 th percentile were coded as 1, and as 0 otherwise.