Jean-Daniel Bontemps
email: jean-daniel.bontemps@engref.agroparistech.fr

Jean-Christophe Hervé

Jean-Michel Leban

Jean-François Dhôte
email: jean-francois.dhote@onf.fr

Jean-Christophe Herve

Nitrogen footprint in a long-term observation of forest growth over the twentieth century

Keywords: forest growth, long-term trends, nitrogen, climate, carbon dioxide, Fagus sylvatica

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

INTRODUCTION

Increases in temperate and boreal forest productivity have been evidenced over the past several decades [START_REF] Boisvenue | Impacts of climate change on natural forest productivityevidence since the middle of the 20th century[END_REF]Spiecker et al. 1996). Because there are several potential environmental drivers, including increased atmospheric CO 2 concentration, climate change, and nitrogen deposition, the role of each has been much debated [START_REF] Hyvönen | The likely impact of elevated [CO 2 ], nitrogen deposition, increased temperature and management on carbon sequestration in temperate and boreal forest ecosystems: a literature review[END_REF]). Contributing to the complexity of the issue are potentially positive and negative interactions among these factors, dependencies of responses on environmental limitations of the ecosystems under consideration, and delays in ecosystem response [START_REF] Ollinger | Interactive effects of nitrogen deposition, tropospheric ozone, elevated CO 2 and land use history on the carbon dynamics of northern hardwood forests[END_REF].

Existing contributions on environmental factor effects on increased productivity are consistently indecisive. Retrospective growth studies have demonstrated a global fingerprint of climate change [START_REF] Boisvenue | Impacts of climate change on natural forest productivityevidence since the middle of the 20th century[END_REF][START_REF] Esper | Low-frequency signals in long tree-ring chronologies for reconstructing past temperature variability[END_REF]. The enhancing effect of CO 2 increase on productivity has been putatively inferred by observations (e.g. [START_REF] Lamarche | Increasing atmospheric carbon dioxide: tree ring evidence for growth enhancement in natural vegetation[END_REF] and ecosystem experiments [START_REF] Norby | Forest response to elevated CO 2 is conserved across a broad range of productivity[END_REF], but this effect has also been questioned [START_REF] Jacoby | Tree rings, carbon dioxide, and climatic change[END_REF][START_REF] Körner | Carbon flux and growth in mature deciduous forest trees exposed to elevated CO 2[END_REF]. There is growing awareness of the fact that impacts of CO 2 increase may depend on ecosystem nitrogen availability [START_REF] Nowak | Functional responses of plants to elevated atmospheric CO2 -do photosynthetic and productivity data from FACE experiments support early predictions?[END_REF]. Reports indicate that increasing nitrogen supply throughout Europe has stimulated past growth (Kahle et al. 2008a) and also has effects on recent growth rates (e.g. [START_REF] Solberg | Analyses of the impact of changes in atmospheric deposition and climate on forest growth in European monitoring plots: a stand growth approach[END_REF]). However, positive effects of nitrogen deposition have been challenged [START_REF] Nadelhoffer | Nitrogen deposition makes a minor contribution to carbon sequestration in temperate forests[END_REF]. Negative impacts of nitrogen deposition are also supported by evidences of nitrate leaching and soil acidification [START_REF] Aber | Is nitrogen deposition altering the nitrogen status of northeastern forests?[END_REF][START_REF] Macdonald | Nitrogen input together with ecosystem nitrogen enrichment predict nitrate leaching from European forests[END_REF], which stem from excess nitrogen, called nitrogen saturation [START_REF] Aber | Nitrogen saturation in temperate forest ecosystems -Hypotheses revisited[END_REF]. The role of nitrogen availability thus remains controversial.

There have been 3 approaches to determining the contribution of environmental factors to productivity changes: ecosystem experiments, process-based modelling analyses, and observations of past growth. Ecosystem experiments have included nitrogen fertilization/suppression manipulations [START_REF] Wright | Introdution to the NITREX and EXMAN projects[END_REF] and free air CO 2 enrichment (FACE, [START_REF] Hendrey | A free-air enrichment system for exposing tall forest vegetation to elevated atmospheric CO 2[END_REF]. Interactions between these factors have also been explored experimentally [START_REF] Nowak | Functional responses of plants to elevated atmospheric CO2 -do photosynthetic and productivity data from FACE experiments support early predictions?[END_REF]). These experiments have limitations as they may reveal (or miss) transient (or delayed) responses due to their short duration. They also diverge from real conditions when factors are manipulated (e.g., [START_REF] Högberg | Tree growth and soil acidification in response to 30 years of experimental nitrogen loading on boreal forest[END_REF]. Process-based models encompass fundamental ecosystem processes that have explicit dependence on the environment [START_REF] Constable | Suitability of process-based tree growth models for addressing tree response to climate change[END_REF][START_REF] Landsberg | Physiology in forest models: history and the future[END_REF][START_REF] Luckai | Challenges in the application of existing process-based models to predict the effect of climate change on C pools in forest ecosystems[END_REF]. They have thus contributed to quantifying ecosystem responses to environmental changes [START_REF] Van Oijen | Evaluation of past and future changes in European forest growth by means of four process-based models[END_REF] and to the identification of key control processes (e.g. [START_REF] Medlyn | Soil processes dominate the long-term response of forest net primary productivity to increased temperature and atmospheric CO 2 concentration[END_REF]. However, their predictive ability is challenged by the frequently incomplete description of processes and uncertainties in parameter estimates [START_REF] Luckai | Challenges in the application of existing process-based models to predict the effect of climate change on C pools in forest ecosystems[END_REF][START_REF] Landsberg | Physiology in forest models: history and the future[END_REF][START_REF] Mäkelä | Process-based models for forest ecosystem management : current state of the art and challenges for practical implementation[END_REF]. Growth records offer a unique procedure for tracking past long-term dynamics of ecosystem productivity as an integrated response to environmental changes. They are thus crucial for testing the outcomes of short-term experiments and model explorations. In the growth record approach, environmental causality is investigated by comparing temporal and spatial patterns of environmental factors and productivity [START_REF] Jacoby | Tree rings, carbon dioxide, and climatic change[END_REF].

Empirical productivity-environment relationships can also aid in the determination of limiting factors [START_REF] Seynave | Is the spatial distribution of European beech (Fagus sylvatica L.) limited by its potential height growth?[END_REF].

However, the success of this observation-as-a-baseline strategy has not been comprehensively successful to date. Dendrochronology has contributed mainly to documenting growth trends from tree core samplings [START_REF] Boisvenue | Impacts of climate change on natural forest productivityevidence since the middle of the 20th century[END_REF][START_REF] Jacoby | Tree rings, carbon dioxide, and climatic change[END_REF]). Yet, radial growth indicators are sensitive to between-tree competition and may reflect past disturbances or silvicultural shifts. Whether they accurately depict ecosystem productivity is also an issue for debate [START_REF] Bouriaud | Is ring width a reliable proxy for stembiomass increment? A case study in European beech[END_REF]. Productivity-oriented indicators have also been investigated in Europe (Spiecker et al. 1996), but growth chronologies have seldom been determined, and the heterogeneity of approaches has hindered regional comparisons. Hence, it has not been possible to attribute causality to temporal shifts in environmental factors. To our knowledge, only 2 observational studies have addressed growth trends and their causality so far. [START_REF] Nellemann | Long-term changes in forest growth: potential effects of nitrogen deposition and acidification[END_REF] compared radial growth chronologies of Norway spruce over a gradient of nitrogen deposition in Norway. However, stand density was not controlled, and other drivers were not addressed. Kahle et al. (2008b) compared tree height growth data with climate and foliar nitrogen trends in 3 northern European species. Confounding of factors, however, prevented unambiguous interpretation and attribution of causality. In addition, both analyses were restricted to medium-term growth variations.

We previously developed a method for quantifying long-term growth trends in even-aged stands [START_REF] Bontemps | Long-term changes in forest productivity: a consistent assessment in even-aged stands[END_REF]) using: (i) dominant height growth reconstructed from stem analyses to serve as a proxy for aboveground productivity [START_REF] Assmann | The principles of forest yield study[END_REF]), (ii) the paired-plots sampling method, which allows for comparisons between young and old stands located in the same site conditions, and (iii) a statistical modelling approach for estimating growth chronology, factoring out the effects of ontogeny and site. We demonstrated a productivity increase of 50% throughout the twentieth century for common beech (Fagus sylvatica L.) in north-eastern France, mostly occurring in the 1950-1980s. The pattern we observed suggested a potential role for nitrogen deposition, in the frame of the nitrogen limitation hypothesis [START_REF] Johnson | Progressive N limitation in forests: review and implications for long-term responses to elevated CO 2[END_REF].

Accordingly, we repeated the study in an area with sites of contrasting nutritional status and nitrogen deposition level. Our objectives were: (i) to determine the regional variation in past productivity changes, (ii) to draw inferences on the environmental drivers of these changes, based on their documented historical variations and regional growth-environment relationships, and (iii) to test the relevance of the nitrogen driver hypothesis.

MATERIALS

Study area and forest species

The present study was based on temporary plots. We focused on even-aged stands of common beech in the northern plains sector of the species distribution [START_REF] Seynave | Is the spatial distribution of European beech (Fagus sylvatica L.) limited by its potential height growth?[END_REF]. To ensure management continuity through time, we restricted the sampling to State Forests. Beech has not been subjected to genetic selection programmes, and natural regeneration predominates, indicating that any historical changes are highly unlikely to be the result of human technical intervention. In addition, past disturbances caused by silvicultural practice were unlikely to have impacted the present sample, as beech stands are attested to have been managed as high forests from establishment onwards (Hüffel 1926). The historical intensity of forest management is moreover considered moderate for this species [START_REF] Polge | The influence of thinnings on the growth constraints of Beech (in French)[END_REF].

We selected the north-western Atlantic area (sample 2) of common beech distribution in France to complement the initial north-eastern sampling area (sample 1), because soils of the Atlantic area have lower nutritional status (European Commission 2005) and because current rates of total atmospheric nitrogen deposition (4-7 kg.ha -1 .year -1 ) in western France are lower than in eastern areas (8-13 kg.ha -1 .year -1 ) [START_REF] Croisé | Two independent methods for mapping bulk deposition in France[END_REF].

Dominant height as a productivity indicator

Dominant height corresponds to the height of the most vigorous trees in a pure and even-aged stand, and was defined as the mean height of the 100 thickest trees in 1 ha. Dominant height and aboveground net primary productivity (wood ANPP) are closely related, regardless of site fertility and stand density conditions [START_REF] Assmann | The principles of forest yield study[END_REF]. The reason for this lies in the insensitivity of both dominant height [START_REF] Lanner | On the insensitivity of height growth to spacing[END_REF] and ANPP (law of constant yield, [START_REF] Bormann | Stand density effects in young Red Alder plantations: productivity, photosynthate partitioning and nitrogen fixation[END_REF] to density within the usual range of closed-canopy stands. Hence, dominant height growth rate was used as a proxy for forest productivity. Stem analysis allows individual tree height growth to be reconstructed [START_REF] Curtis | A stem analysis approach to site index curves[END_REF]. We applied the procedure to samples of dominant trees. For each tree, the first 2 or 3 disc sections were sawn every 4 m along logs beginning 0.30 m above ground level, and then every 2 m in the upper part of the stem to the terminal bud.

Sampling design: controlling for tree age and site conditions

We partitioned out effects on growth of historical changes, ontogeny and site fertility by using a paired-plots sampling design, in which plots were sampled in pairs of young and old stands located within the same sites. This sampling also ensures some homogeneity of historical management and genetic similarity among trees. To encompass site variation within each area, we sampled around 15 stand pairs/area. Candidate pairs were first identified using forest management and geological maps. Old stands were selected around the standard rotation age (150-180 years) and young stands were half this age (70-80 years). Pairs were sampled after controlling for stand structure, species homogeneity, permanent environmental characteristics (topography, parent rock and soil) and current nutritional status (humus forms and summer surveys of understory vegetation).

Growth data

Information on the 14 pairs in sample 1 is provided in [START_REF] Bontemps | Long-term changes in forest productivity: a consistent assessment in even-aged stands[END_REF]. Sample 2 comprised 15 pairs whose location and age are indicated in the Online Resource. Both samples are mapped in Figure 1. In each stand, dominant trees (defined above) were sampled following the protocol of [START_REF] Duplat | Modelling of dominant height growth in Sessile oak (Quercus petraea Liebl) in France. Regional variability and effect of the recent period (1959-1993)[END_REF]; first, the five thickest trees in a circular plot of 0.06 ha were selected to avoid local sampling bias (Pierrat et al. 1995). Next, the first, third and fifth thickest trees were selected for stem analysis to reduce the sampling effort. Individual growth curves of these trees were averaged to give a plot mean curve. In total, 174 stem analyses were performed. From samples 1 and 2, respectively, we obtained 356 and 390 growth increments, covering on average 7.2 and 8.3 years (SD: 3.9 and 4.7 years), respectively. Height curves plotted for each sample are compared between generations in Figure 2.

Environmental data

In all plots, current local nutrition and water status were documented from soil analyses and understory vegetation-based indicators. Indicators of mean climate ) at 1 km spatial resolution were extracted for each stand pair (mean coordinates of the 2 plots) from the AURELHY model [START_REF] Bénichou | Prise en compte de la topographie pour la cartographie des champs pluviométriques statistiques[END_REF].

Soil analyses and vegetation-based indicators

Soil layer textures were analysed to estimate the 1 m depth soil water capacity (SWC100). Cation

(cobaltihexamin extraction at soil pH), pH, total nitrogen (N), organic carbon (C) and phosphorus (Phos, sodium bicarbonate extraction at pH 8.5 or Olsen's method) concentrations were also measured in each soil layer to a depth of 30 cm. C:N ratio, base cation concentration (S=Ca+K+Mg+Na), cation exchange capacity (CEC=S+Al+H+Fe+Mn) and base saturation rate (S:CEC) were also computed. Selected vegetation-based Ellenberg indicators [START_REF] Ellenberg | Zeigerwerten von pflanzen in Mitteleuropa[END_REF] included nitrogen (Nel), basicity (Rel), humidity (Fel) and ambient light (Lel) values. According to this system, each plant species is rated on an increasing relative 9-point scale for each factor.

Averages for each indicator were computed over each vegetation survey. Soil C:N ratio and the nitrogen Ellenberg indicator have been correlated with soil nitrogen mineralization [START_REF] Janssen | Nitrogen mineralization in relation to C:N ratio and decomposability of organic materials[END_REF][START_REF] Schaffers | Reliability of Ellenberg indicator values for moisture, nitrogen and soil reaction: a comparison with field measurements[END_REF] and nitrification [START_REF] Andrianarisoa | Comparing indicators of N status of 50 beech stands (Fagus sylvatica L.) in northeastern France[END_REF][START_REF] Schaffers | Reliability of Ellenberg indicator values for moisture, nitrogen and soil reaction: a comparison with field measurements[END_REF] and were consequently selected as measures of soil nitrogen availability.

Mean climate data

Mean climate indicators included monthly precipitation (P), minimum, maximum and mean monthly temperatures (Tn, Tx, Tm) and annual number of freezing days (NFD). Monthly potential evapotranspiration (PET) values were calculated from temperature and solar radiation (Rad; [START_REF] Piedallu | Multiscale computation of solar radiation for predictive vegetation modelling[END_REF]. Monthly water balances (WB) were calculated as P-PET. Monthly soil water budgets (SWB) and soil water deficits (SWD) were calculated from SWC100, P and PET.

Monthly indicators were indexed by corresponding calendar month numbers (e.g. Tx7 for maximum temperatures in July).

Regional chronologies of environmental factors

We assembled data for past climate, atmospheric CO 2 concentration and nitrogen deposition. As is general in retrospective sampling studies performed in temporary plots, no environmental monitoring has been conducted in the past, and information on the historical nutritional status of these forest stands is absent. The 1901-2000 chronologies for precipitation and mean temperature were extracted for each stand pair from the CRU TS 1.2 monthly dataset for Europe at a 10′ angular (approximately 15 km) spatial resolution [START_REF] Mitchell | A comprehensive set of highresolution grids of monthly climate for Europe and the globe : the observed records (1901-2000) and 16 scenarios[END_REF]) and averaged over each regional sample. Geographic variation of CO 2 concentration is reportedly marginal [START_REF] Chahine | On the determination of atmospheric minor gases by the method of vanishing partial derivatives with application to CO 2[END_REF] airs.jpl.nasa.gov/AIRS_CO2_Data/). A single chronology was thus computed from CDIAC datasets (cdiac.ornl.gov/). Although spatial patterns of recent deposition are known from monitoring networks [START_REF] Holland | Nitrogen deposition onto the United States and Western Europe: synthesis of observations and models[END_REF] [START_REF] Croisé | Two independent methods for mapping bulk deposition in France[END_REF][START_REF] Ulrich | Estimation of nitrogen deposition on 27 RENECOFOR plots (France) from 1993 to 1996[END_REF]). Bulk deposition is defined as the total deposition of material (wetonly deposition, dry particles, gases and fine aerosols) to a continuously-open sampler, and it is generally lower than the total deposition intercepted by vegetated surfaces [START_REF] Dämmgen | Practical considerations for addressing uncertainties in monitoring bulk deposition[END_REF].

Because it is independent from forest canopy structure and composition, it can be modelled and spatially extrapolated using abiotic predictors. Bulk depositions of NO 3 -and NH 4 + monitored in the CATANEAT network over the period 1993-1998 have been modelled as a function of elevation, monthly precipitation and period of the year. These models have been validated against deposition data for the period 1999-2001 and they have shown a good predictive ability [START_REF] Croisé | Two independent methods for mapping bulk deposition in France[END_REF].

We applied these models to the 1993-1998 CRU climate data extracted for each stand pair sampled, and plot elevations. Predictions from these models were averaged over each regional sample. EMEP data for nitrogen deposition were not used, as they apply at a coarse spatial scale (50km x 50 km grid) and have been found to partially capture the spatial distribution of deposition in France [START_REF] Croisé | Two independent methods for mapping bulk deposition in France[END_REF].

Posterior control of stand pairing and variations in site conditions

The accuracies of pairings were assessed by tracking inappropriately-matched individual pairings and systematic between-generation differences in soil-related indicators. Systematic differences were tracked using two-sided paired t tests. Accuracies of pairings were confirmed in sample 1 [START_REF] Bontemps | Long-term changes in forest productivity: a consistent assessment in even-aged stands[END_REF]) and in sample 2 for both criteria.

A selection of important environmental indicators is presented in

    b 2 b 0 1 0 S t t f H f = dt dH  (1)
The function f 2 is such that f 2 (t b ) = 1, with t b taken as 1900. Equation ( 1) was integrated with respect to time, to be fitted to the pluri-annual increment data:

      ) ( ) ( S ) ( ) ( 1 2 2 b 1 0 1 1 1 0     t F t F + t H F F = t H (2) where  ) ( 1 1 u f du = F and  dv v f = F ) ( 2 2
(u and v are dummy variables of integration)

The pattern of growth remained sigmoid in old and young stand generations (Figure 2). Therefore, sigmoid equations were used to model the effect of developmental stage. As an upper limit to height growth was undetected in these data, we selected three 3-parameter asymptotic growth equations showing slow convergence: the Richards, the Korf, and the Hossfeld equations [START_REF] Zeide | Analysis of growth equations[END_REF]. Their basic expression is given by:

  m K, , R 0 0 H f = dt dH (3)
They were parameterized so that R (m.yr -1 ) is the maximal growth rate, K (m) is the asymptotic height and m is a shape parameter (dimensionless). Parameter S b was assimilated to R, and f 1 to f. The asymptotic height parameter K in the growth equation may also depend on site conditions. Hence, its pair-wise variation was tested. As the uppermost sections of the height curves are lacking for the young stand generation (Figure 2), it was however impossible to accurately test whether this parameter also changed over time. The effect of date (f 2 ) was represented by a polynomial cubic spline function for assessing growth fluctuations on a decennial scale. Cubic splines were chosen to prevent excessive oscillations of ordinary polynomials. Mathematical expressions for the selected growth equations (f 1 and F 1 -1 ) and the effects of calendar date (f 2 ) are given in the Online Resource.

An in-depth justification of the model structure is provided in [START_REF] Bontemps | Long-term changes in forest productivity: a consistent assessment in even-aged stands[END_REF].

Data are longitudinal and have 2 hierarchical levels (pairs as level 1, stands within pairs as level 2) at which model parameters can vary. The model was fitted using a non-linear hierarchical mixedeffects model approach [START_REF] Lindström | Nonlinear mixed effect models for repeated measures data[END_REF]. Level 1 and 2 variations in maximal growth rate S b should correspond to the effect of site and residual within-pair site differences, and were tested. Level 1 variation of K and its correlation with S b were also tested. The random parameters were assumed to follow Gaussian distributions [START_REF] Lindström | Nonlinear mixed effect models for repeated measures data[END_REF]:

S b ~  Sb,2 Sb,1 b0 σ σ , S  N , K 1 ~   K,1 0 σ , K N ,   K , S cor ρ b 1 = (4)
where,  Sb,1 and  Sb,2 are random standard-deviations of S b at levels 1 and 2,  K,1 that of K at level 1,   is their correlation and S b0 and K 0 are fixed components of S b and K. Equation ( 2) was fitted to increments (forward difference data) using the maximum likelihood criterion of the nlme procedure in S-PLUS software [START_REF] Pinheiro | Deposition of fixed atmospheric nitrogen and foliar nitrogen content of Bryophytes and Calluna vulgaris (L.) Hull[END_REF]. Errors were assumed Gaussian and their standard-deviation (residual standard error) was assumed to be proportional to a power () function of the expected height increment. Selection of a growth equation and adequate parameterization for random effects and cubic splines were performed using the AIC criterion and the Khi-2 likelihood ratio test (LRT) for nested models.

Growth-environment relationships

To identify environmental factor constraints on growth, relationships between maximal height growth rate (S b ) and soil and mean climate indicators relative to each plot were investigated in each regional sample (28 or 30 estimates of parameter S b and environmental indicators, depending on the sample) and over the combined dataset. To make the analysis comparable between the 2 regions, estimates of S b were obtained using a common growth equation (see Results). Growth-environment relationships were established in 2 steps. First, correlations between maximal growth rate and environmental indicators were systematically computed and tested using Pearson's correlation t tests at a minimum significance level of 0.05. Significant correlations of the same sign over successive months were interpreted as seasonal patterns. To reduce inflation in variable numbers, sums (precipitations) or means (temperatures) of monthly climatic indicators were calculated over the corresponding periods whenever seasonal patterns were observed. Correlations between these seasonal indicators and maximal growth rates were again computed. Secondly, we built a multiple linear regression model of maximal growth rate, starting with indicators with most significant correlations. Indicators with weaker correlations were incorporated into the regression model until a minimum level of significance of 0.05 was exceeded. When both monthly and seasonal indicators of the same climatic variable were significant in the regressions, we systematically retained the model including seasonal indicators (Table 4), considered to be biologically more meaningful.

RESULTS

Historical changes in dominant height growth rate

Pair-wise comparisons of dominant height at the final ages of the younger stands indicated that younger stands were systematically taller than the older ones (Figure 2). Mean differences were highly significant, and were up to 5.7 m in sample 1 (p < 10 -4 , mean comparison age 75.4 years) and 2.2 m in sample 2 (p < 10 -3 , mean comparison age 73.5 years). The difference was significantly smaller in sample 2 than in sample 1 (p = 0.002).

Main characteristics of the growth models are reported in Table 2. The Hossfeld and Korf equations proved to be the most accurate for samples 1 and 2, respectively. In sample 2, the Hossfeld equation provided a lower but nevertheless very close fitting accuracy (+4.5 points in AIC, comparable RSE, Table 2). The random structure selected was the same for both samples: a level 1 random variation of parameters S b and K (the correlation between S b and K was not significant), and a level 2 variation of S b (Table 2). Level 1 standard-deviations of random effects were below 10% for K and ca. 15% for S b . Cubic splines for the effect of date were best parameterized with 15-and 20-year internodes in samples 1 and 2, respectively.

Regional growth chronologies are plotted in Figure 3. Consistent with the mean height differences between stand generations (Figure 2), increases in dominant height growth rate were identified and appeared to be irregular over time. Whereas no clear pattern was evident in the first half century, a noticeable increase occurred in following decades. Growth chronologies also revealed synchronous growth declines in the 2 samples, centred in the 1940s and 1990s, although the first was undetected by the sample 2 spline estimate, but appeared in the growth increments. The decline in the 1990s was also shown to be more acute in sample 2, in which the average growth rate was similar to the a indexes denote the level of random variation in standard deviation estimates, b in the absence of a variance function, c The best fit was obtained with the Korf equation for sample 2. The Hossfeld equation fit is also provided for comparison with sample 1, as it was used for elaboration of growthenvironment relationships (Tables 3 and4). Standard-errors of parameter estimates are in parentheses. level in 1900. The biggest difference between chronologies was apparent in the intensity of regional changes in growth rate. While the increase in growth rate reached +60% in sample 1 around 1980, it barely reached +30% in sample 2 before flattening off. Hence, the magnitude of change in sample 1 was double that in sample 2. The 95% confidence intervals for the estimated trends stopped overlapping between samples from the 1960s onwards.

Growth-environment relationships

Although the Korf equation provided the best fit for sample 2, the fit for the Hossfeld equation was very close (only +4.5 point difference in AIC). Plot estimates for maximal growth rate (S b ) were thus extracted from a Hossfeld-based model, as in sample 1 (the difference between these estimates was tiny; see Online Resource). Significant correlations between maximal growth rate and environmental indicators are reported in Table 3.

Sample 1

We identified negative correlations with soil acidity level (pH, aluminium, Ellenberg indicator for basicity) and positive correlations with soil nutrition indicators, including nitrogen (C:N ratio, Ellenberg indicator for nitrogen), phosphorus and base cation richness (S:CEC ratio and calcium concentration). A single correlation with a water indicator was demonstrated (precipitation in September), suggesting the absence of particular water constraints in the area. There were positive correlations with minimum temperatures in the early growing season (February to June), and strong negative correlations with maximum temperatures in the cold season (October to February). Finally, there was a negative correlation with ambient light (Lel) and summer radiation (June to August).

Due to the seasonal patterns detected for temperature and radiation, we computed the following mean indicators (correlations coefficients are given in parentheses): Tn4-5 (+0.40), Tn4-6 (+0.38) and Tn2-6 (+0.40), Tx12-1 (-0.63), Tx12-2 (-0.64), Tx10-1 and Tx10-2 (-0.61), and Rad6-8 (-0.49).

We incorporated these indicators into a multiple linear regression model of maximal growth rate following the procedure outlined in the Modelling Approach section. The final regression model is presented in Table 4.

Sample 2

We also identified positive correlations with soil nutrition indicators for sample 2, including nitrogen (C:N ratio, Ellenberg indicator for nitrogen) and cation richness (CEC, S:CEC ratio, concentration of base cations and calcium). A negative correlation was found with proton indicator for humidity, and soil water deficits in June and July). As in sample 1, we identified a negative effect of maximal temperatures in the cold season (October to January). Lastly, there was a positive effect of spring solar radiation (May to July). We then computed the following temperature averages: Tx11-12 (-0.60), Tx10-12 (-0.50) and Tx10-1 (-0.50). The regression model is presented in Table 4.

Samples 1+2

When samples 1 and 2 were combined, positive correlations were confirmed for nitrogen (C:N ratio, Ellenberg indicator for nitrogen), phosphorus concentration, and S:CEC ratio (Table 3).

Correlations with minimum temperatures were significant for all months, and we retained the mean annual minimum temperature. Negative correlations were found with maximum temperatures in September and October. The final regression model is reported in Table 4.

DISCUSSION

Regional chronologies and intensity of productivity changes

A major finding of this study was the strong regional difference in the magnitude of historical evolutions of dominant height growth rate, while temporal patterns were highly similar. The regional variation detected was thus a matter of signal intensity (Figure 3). Most of the changes occurred in the second half of the twentieth century and differed by a factor of 2 at their maximum level in the 1980s (+60% and +30% for samples 1 and 2, respectively). The growth decline in the 1990s was stronger in sample 2. The growth increase over time has been progressive reaching +27% and +6% on average in samples 1 and 2 over the twentieth century, respectively. The relationship between dominant height and total stand production [START_REF] Assmann | The principles of forest yield study[END_REF] suggests similar increases in productivity.

Interestingly, these magnitudes of increase are similar to those of retrospective process-based simulations. For sites spanning a latitudinal gradient in Europe, [START_REF] Van Oijen | Evaluation of past and future changes in European forest growth by means of four process-based models[END_REF] reported NPP increases ranging from +5 to +20% over the period 1920-2000 as a combined effect of changes in atmospheric CO 2 , nitrogen deposition and climate. A simulation of climate change and CO 2 impacts on common beech in northern France for the period 1960-2000 [START_REF] Loustau | Modeling climate change effects on the potential production of French plains forests at the sub-regional level[END_REF] estimated an NPP increase of 13%, a value that is between our 2 estimates.

Environmental interpretation of changes in productivity

The following necessary conditions should be met when attempting to attribute changes in productivity to environmental factors: (i) each putative factor must have experienced variation over time; (ii) temporal patterns of factor(s) and growth should match, assuming that the growth response is not delayed and that growth saturation is not attained (temporal argument); (iii) regional variations in growth changes must follow those of factor(s) (spatial argument), and (iv) the same factor(s) must be identified as limiting for growth, which further requires that no other factor is sufficiently constraining to totally suppress the productivity response (limiting factor argument).

The role of each putative environmental factor is now discussed with respect to these conditions.

The resulting deductions are summarized in Table 6.

Atmospheric carbon dioxide

Atmospheric CO 2 concentration increased from about 295 ppm in 1900 to 370 ppm in 2000, with acceleration in the 1950s, reaching a mean annual progression of 1.2 ppm yr -1 during the period 1950-2000 compared to 0.3 ppm yr -1 during the period 1900-1950 (CDIAC data). Hence, CO 2 chronology is not contradictory to trends reported. Because there is no spatial variation in CO 2 concentration at these geographic scales [START_REF] Chahine | On the determination of atmospheric minor gases by the method of vanishing partial derivatives with application to CO 2[END_REF], its limiting effect on productivity was not testable. To account for regional differences emphasized in growth chronologies, a combined role of either permanent constraints or synergies with other drivers is required (Table 6).

Temperature and precipitation

The mean temperature chronology for France [START_REF] Moisselin | Climatic changes in France during the 20th century. Analysis of long homogenised time series of temperature and precipitation[END_REF]) matches northern hemisphere chronology [START_REF] Jones | Hemispheric and Large-Scale Surface Air Temperature Variations: An Extensive Revision and an Update to 2001[END_REF]. There was a centennial increase of +1°C, with stability in the period in the 1940-1980s, after which temperature has been steadily rising. The stationary midperiod has been referred to repeatedly since [START_REF] Jones | Variations in surface air temperatures : Part 1. Northern hemisphere, 1881-1980[END_REF]. No significant trend in precipitation has been reported over the last century [START_REF] Moisselin | Climatic changes in France during the 20th century. Analysis of long homogenised time series of temperature and precipitation[END_REF]. Annual anomalies for mean temperature and precipitation (Figure 4) are comparable between our samples. Annual fluctuations were higher in sample 1, especially for precipitation.

The main objection to the role of temperatures in measured growth accelerations is the fact that temperatures remained stationary during the period 1950-1980 (Figure 4). Growth-environment relationships further suggest antagonistic effects of temperature (Table 4): positive effects of minimum annual temperatures (Tn) and of minimum spring temperatures (Tn4-6) were identified throughout the whole study area and in sample 1 (colder climate than in sample 2). Conversely, a negative effect of maximum cold season temperatures was detected in both samples. Both effects are of comparable magnitude and might have counteracted each other. When combined, these observations invalidate the role of temperature in the growth increases reported (Table 6).

Precipitation anomalies (Figure 4) followed no long-term trend paralleling productivity increases.

In sample 1, no water constraint was identified. In sample 2, where precipitation was lower (Table 1), SWC100 was significant in the model of maximal growth rate (Table 4). Thus, precipitation was not identified as a low-frequency driver of productivity changes, but may have acted as a permanent constraint in that area.

Nevertheless, singular climatic years may have played a significant role in decennial fluctuations in growth (Figures 3 and4). The positive precipitation anomalies observed for in the 1980s were synchronous with maximum levels of the growth rate. Growth declines were also synchronous with exceptional successions of warm years during the periods 1943-1949 and 1989-2000. Exceptional tree mortality following droughts in the period 1947-49 has also been reported for common beech in France [START_REF] Schaeffer | Decline in common beech in the administrative unit of Doubs (in French)[END_REF]. Lower water availability in sample 2 (Tables 1 and3) may also have resulted in a greater impact of dry years during the 1990s and may explain the higher intensity of the observed decline (Table 6).

Nitrogen availability and nitrogen deposition

A comprehensive record of past nitrogen deposition is lacking. However, there is a consensus that the mid-twentieth century increase tracked the course of emissions resulting from increasing fossil fuel combustion and fertilizer use in Europe [START_REF] Asman | Modelled historical concentrations and depositions of ammonia and ammonium in Europe[END_REF][START_REF] Holland | Nitrogen deposition onto the United States and Western Europe: synthesis of observations and models[END_REF]. This is corroborated by existing measurements [START_REF] Blake | Changes in soil chemistry accompanying acidification over more than 100 years under woodland and grass at Rothamstead Experimental Station, UK[END_REF], ice-core records [START_REF] Hastings | Anthropogenic impacts on nitrogen isotopes of ice-core nitrate[END_REF] and vegetation reconstitutions (Pitcairn et al. 1995). In France, the Alpine ice-core reconstitution of atmospheric nitrate [START_REF] Preunkert | A seasonally resolved alpine ice core record of nitrate : Comparison with anthropogenic inventories and estimation of preindustrial emissions of NO in Europe[END_REF]) also agrees with the records. Our reported evolutions in growth rate would thus be temporally consistent with the course of nitrogen deposition (temporal argument, Table 6).

Despite regional differences in nitrogen status (Table 1), we determined a significant positive effect of current nitrogen availability on growth (using both vegetation-related and soil indicators) in each sample and in the combined dataset (Table 4). Hence, nitrogen is identified as a limiting factor, even in the better nutritional conditions of sample 1, which is consistent with the findings of [START_REF] Seynave | Is the spatial distribution of European beech (Fagus sylvatica L.) limited by its potential height growth?[END_REF] for common beech (limiting factor argument, Table 6).

Sample 1 Sample 2 Difference (%) (sample 1)

Elevation (m a.s.l.) 371 ( 48) 125 ( 51)

Annual precipitation 1993-1998 (mm) 825 ( 60) 670 ( 52)

NO 3 --N (kg ha -1 yr -1 ) 3.3 (0.1) 2.6 (0.2) +28%

NH 4 + -N (kg ha -1 yr -1 ) 5.4 (0.2) 3.8 (0.4) +44%

Total N (kg ha -1 yr -1 ) 8.9 (0.3) 6.5 (0.6) +37%

Table 5. Modelled average level of bulk nitrate, ammonium and total nitrogen deposition for the period 1993-1998 for growth samples. Deposition levels were computed from models of NH 4 -N, NO 3 -N and total N bulk deposition calibrated from the RENECOFOR monitoring network (ICP forests) as a function of elevation, monthly precipitation and period of the year [START_REF] Croisé | Two independent methods for mapping bulk deposition in France[END_REF]. Standard deviations are given in parentheses.

Bulk nitrogen deposition levels predicted over the two samples remain moderate, as they do not exceed 10 kg.ha -1 .yr -1 (Table 5; [START_REF] Ulrich | Estimation of nitrogen deposition on 27 RENECOFOR plots (France) from 1993 to 1996[END_REF]. However, bulk deposition in the open field generally underestimates total deposition to the ecosystem [START_REF] Dämmgen | Practical considerations for addressing uncertainties in monitoring bulk deposition[END_REF][START_REF] Ulrich | Estimation of nitrogen deposition on 27 RENECOFOR plots (France) from 1993 to 1996[END_REF], in general by 20-30% [START_REF] Croisé | Two independent methods for mapping bulk deposition in France[END_REF]. Average 1993-1998 bulk deposition levels of nitrate and ammonium were found higher in sample 1 (Table 5), which had a very significant +37% total nitrogen deposition excess over sample 2 (p = 7 10 -11 ). This regional pattern has been found consistent with levels of precipitation higher in north-eastern than in north-western France (Table 1; [START_REF] Croisé | Two independent methods for mapping bulk deposition in France[END_REF], as observed on the CATANEAT monitoring network [START_REF] Ulrich | Estimation of nitrogen deposition on 27 RENECOFOR plots (France) from 1993 to 1996[END_REF]. Given the known determinants of nitrogen deposition (emission hotspots, topography, climatic conditions), this predicted difference may have been in effect over a longer time scale. Thus, it is in accordance with the regional differentiation of productivity changes (spatial argument, Table 6), and it may have increased the regional difference in soil richness and nitrogen status (see Materials and Table 1). As the study was carried out on temporary plots, we used modelled data of bulk nitrogen deposition.

Despite these models have shown a good predictive ability [START_REF] Croisé | Two independent methods for mapping bulk deposition in France[END_REF], these predictions remain subject to uncertainty. In addition, total bulk deposition in the open field does not take canopy interception and exchanges into account (Draajers et al. 1997). For these reasons, predicted nitrogen deposition levels and their regional difference may only deliver a rough picture of real deposition to the forest stands sampled.

Environmental change Temporal argument

Consistency of historical variations with growth chronologies on the long-term

Spatial argument

Consistency with regional differentiation of growth chronologies

Limiting factor argument

Constraint from growthenvironment relationships

Conclusion

Atmospheric CO 2 YES (CDIAC datasets)

No regional variation [START_REF] Chahine | On the determination of atmospheric minor gases by the method of vanishing partial derivatives with application to CO 2[END_REF] Impossible to determine by observation

Compatible

Cannot have acted alone Temperature NO (Figure 4)

Warm events may account for 1940s and 1990s growth declines NO (Figure 4) Antagonistic effects evidenced (Tables 3 and4)

Rejected

Role in short-term events

Precipitation

No long-term variation (Figure 4)

Wet events may account for growth increase in the 1980s NO (Figure 4)

YES in sample 2 (Table 3 and4

YES

(Table 5)

YES in both samples (Tables 3 and4) may be more pronounced in sample 2

Compatible

Table 6. Compatibility of environmental changes with long-term growth changes in common beech. Tables, figures or growth sources supporting the analysis are in parentheses.

Because no significant negative effect of acidity indicators was detected (though some trends were apparent sample 1), the eutrophicating role of nitrogen deposition may have overruled its acidifying side effect [START_REF] Macdonald | Nitrogen input together with ecosystem nitrogen enrichment predict nitrate leaching from European forests[END_REF]. Less acidic soil conditions in sample 1 (Table 1) may have also prevented such an impact. In addition, deposition levels were below the threshold of 10 kg ha -1 yr -1 , above which nitrogen leaching and saturation may triggered (Dise and Right 1995). Our results thus support the existence of nitrogen limitation during the past century [START_REF] Johnson | Progressive N limitation in forests: review and implications for long-term responses to elevated CO 2[END_REF]) and a growth enhancement by nitrogen deposition (Table 6).

Additional roles of regional environmental constraints and local events

Nutritional and water resources were less favourable in sample 2 (Table 1), and both were shown to influence growth (Table 4). The SWC100 was the single significant water indicator in the regression model of maximal growth rate, but spring/summer precipitation and water balances were positively correlated with maximal growth rate (water deficits in summer negatively correlated), which was not the case in sample 1 (Table 3). Therefore, as a permanent constraint, lower water availability may have hindered the productivity response in the area. However, other constraints were found to be specific to sample 1, such as phosphorus limitation (Phos) and minimum temperatures in the early growing season (Tn4-6). While soil concentrations were equivalent in samples 1 and 2, phosphorus limitation in sample 1 may be interpreted as the indication of a shift in growth control following increased nitrogen availability. This finding is consistent with the progressive phosphorus limitation evidenced in recent nitrogen fertilization experiments [START_REF] Braun | Does nitrogen deposition increase forest production? The role of phosphorus[END_REF]).

In the Atlantic beech forest domain, unusual tree mortality has been reported in the period 1965-1980, as ensuing from a simultaneous attack of the insect Cryptococcus fagisuga and the fungus Nectria coccinea. Mortality was exacerbated in the years following the drought event in 1976 [START_REF] Mormiche | Managing forest decline of common beech in Normandy 1961-1988 (in French)[END_REF]. While forest declines occurred in 3 of the forests we sampled (Lyons, Eawy and Eu, see Table 1 in Online Resource), increased height growth over this period (Figure 3 and pairwise height comparisons in these forests) suggests, either that these stands were not impacted at all by these disturbances, or that any resulting tree mortality may not have lowered between-tree competition to the point that height growth of the surviving trees would have been reduced.

CONCLUSIONS

We report increases in common beech productivity after the mid-twentieth century and their parallel course with strong differentiation in 2 regions of contrasting nutritional status. Growth and temperature chronologies did not match over the long term. Growth chronologies were consistent with CO 2 chronology, but CO 2 alone cannot account for the regional difference emphasized. The clearest environmental signal was related to nitrogen. This factor met the requirements of temporal and spatial concordance with growth chronologies, and was also identified as limiting for growth.

Nitrogen/CO 2 interaction was not contradicted by our observations, but could not be tested. Shortterm anomalies in precipitation and temperature, including recent climate warming, may have played roles in productivity declines in the 1940s and 1990s and in the productivity peak in the 1980s. Water constraint evidenced in north-western France may have also contributed to this regional differentiation.

Our results thus extend observations on the role of nitrogen over the medium- (Kahle et al. 2008a;[START_REF] Nellemann | Long-term changes in forest growth: potential effects of nitrogen deposition and acidification[END_REF]) and short-terms [START_REF] Solberg | Analyses of the impact of changes in atmospheric deposition and climate on forest growth in European monitoring plots: a stand growth approach[END_REF]) and the inferences of FACE experiments including nitrogen addition [START_REF] Nowak | Functional responses of plants to elevated atmospheric CO2 -do photosynthetic and productivity data from FACE experiments support early predictions?[END_REF]. They are also compatible with retrospective process-based modelling explorations of the role of nitrogen in the productivity of evergreen species in Europe [START_REF] Van Oijen | Evaluation of past and future changes in European forest growth by means of four process-based models[END_REF]. It is remarkable that past increases in NPP simulated for common beech in France [START_REF] Loustau | Modeling climate change effects on the potential production of French plains forests at the sub-regional level[END_REF]) showed no regional differentiation when nitrogen deposition was not taken into account, and were of a smaller magnitude than reported for the north-eastern area.

A comparison of retrospective process-based simulations and observations of past productivity changes remain crucial to further improve our understanding of the environmental origins of these changes. This comparison would also allow testing of the predictive accuracy of process-based models over the long-term [START_REF] Pussinen | Modelling long-term impacts of environmental change on mid-and high-latitude European forests and options for adaptive forest management[END_REF]. However, retrospective modelling explorations of environmental change impacts on productivity are still insufficient [START_REF] Hasenauer | Forest growth response to changing climate between 1961 and 1990 in Austria[END_REF][START_REF] Van Oijen | Evaluation of past and future changes in European forest growth by means of four process-based models[END_REF]. Furthermore, combinations of simulation and observation remain rare (Kahle et al 2008a). Because the parameterisation of process-based models most often requires intensive measurements of simulated forest plots, rapid progress seems unlikely. One option may lie in the use of more parameter-sparse hybrid modelling approaches that combine process-based and statistical models of forest dynamics [START_REF] Landsberg | Physiology in forest models: history and the future[END_REF][START_REF] Matala | Modelling the response of tree growth to temperature and CO2 elevation as related to the fertility and current temperature sum of a site[END_REF].

Growth equations tested

Differential expressions of growth equations (f 1 )

Equations parameterised with K as the height asymptote (metres) and S b as the maximal growth rate (metres/year):

 Richards: All equations admit close-form solutions for F 1 and F 1 -1 (integration on time interval [t -1 , t]):

-Richards:

  m 1 1 2 2 m b m 1 0 0 ) ( ) ( K C m R exp 1 K ) ( 1 K ) (                                    t F t F t H t H -Hossfeld:   ) ( K ) ( ) ( ) ( K C m R 1 K ) ( m 1 m 1 0 1 0 1 2 2 m b 0                          t H t H t F t F t H -Korf:                                   m 1 m 1 0 1 2 2 m b 0 ) ( K ln ) ( ) ( K C m R exp K ) ( t H t F t F t H

Fig. 1 Fig. 2

 12 Fig. 1 Location of samples 1 and 2. Each dot indicates a stand-pair located in either north-eastern or north-western France.

Fig. 3

 3 Fig. 3 Regional chronologies of dominant height growth rate. Thick lines: growth chronologies (cubic spline estimates), thin lines: 95% bilateral confidence intervals, circle symbols: growth increments plotted against their median dates. Dashed lines and empty circle symbols: sample 1 (Bontemps et al. 2009), full lines and solid circles: sample 2.

  ), S (+), S/T (+), Ca (+)

Fig. 4

 4 Fig. 4 Anomalies of mean annual temperature and annual precipitation (1901-2000). a: mean annual temperature, b: annual precipitation. Data extracted from Tyndall dataset CRU TS 1.2 for both samples (see Materials). Anomalies were smoothed using spline interpolation (thick line).

Table 2 . Characteristics and parameter estimates for models of dominant height growth.

 2

Table 3 . Correlations between maximal growth rate in dominant height and environmental indicators.

 3 Pearson correlation coefficients are in parentheses. p < 0.01: numerical correlations in bold, p < 0.05: numerical correlations, p < 0.1: +/-. Significant correlations were most often established with soil indicators in the first soil layer, otherwise(2) refers to soil layer 2. See Table1and text for explanation of abbreviations in Table3cells. concentration. No correlation was found with phosphorus concentration. Positive correlations with indicators of water availability suggested the existence of a water constraint in the region (i.e. precipitation from April to July, water balance in May and June, soil water capacity and Ellenberg

Mean geographic coordinates of stand pairs (ED 50 system),

mean elevation of stand pairs (a.s.l.),

stand age in 1998

ACKNOWLEDGEMENTS

The authors gratefully thank the French Ministry for Agriculture and Fisheries (MAP) and the French Forest Service (ONF) for providing funding and support to the present study. They also wish to thank Antoine Cazin (INRA), Jérôme Piat (ONF) and several ONF services for their helpful assistance in the sampling phase, Daniel Rittié (INRA) for conducting stem analyses, Vincent Perez (AgroParisTech) for climate data extraction, Erwin Ulrich and Manuel Nicolas (ONF) for useful discussions on nitrogen deposition measurements and associated models, and Jean-Claude Pierrat (INRA) and the two anonymous reviewers that greatly helped in clarifying the manuscript.

Parameter Estimation estimate p-value

Proportion of variation explained a

Sample 1

Intercept (m yr -1 ) 0.543 0.015 Tx10-1 (°C) b -0.044 0.020 25.9% Phoslayer 1 (g kg where

Cubic spline function of date (f 2 )

Cubic spline functions are piecewise continuous polynomials defined on successive time intervals. These time intervalsor internodeswere set equal along calendar date. Different internodes were tested to assess decennial-order fluctuations (20, 15 and 10 years). The general expression is:

where u = tt b , d 1 , d 2 , d 3 , p k and pm k are spline parameters to be estimated, [0, n] is the base internode of spline (n = 20, 15 or 10), and n, k 1 and k 2 externally specify the width and number of spline intervals necessary to describe the entire period range covered. For instance with n = 15, nodes are located at dates 1870, 1885, 1900, 1915, …, 1990.

Growth-environment relationships -S b estimates for sample 2

In sample 2, the best fit was obtained with the Korf equation. However, the difference was only -4.5 AIC units compared with the Hossfeld equation, also the most accurate in sample 1. For a clear comparison of regional growth-environment relationships, S b estimates were extracted from a Hossfeld-based model in sample 2.

Differences in S b estimates from Korf-and Hossfeld-based models were investigated (see Figure below). The correlation between both estimates was 1 -(p = 0). When Hossfeld S b estimates were predicted from the Korf ones using simple linear regression, the intercept was not significant: S b,Hossfeld = 0.9816; S b,Korf (RSE = 2.5 10 -3 m/year). Furthermore, the growth-environment relationship in sample 2 was insensitive to the set of S b estimates considered.