# Towards a More Realistic Simulation of Plant Species with a Dynamic Vegetation Model Using Field-Measured Traits: The Atlas Cedar, a Case Study

Alain Hambuckers, Franck Trolliet, Marie Dury, Alexandra-Jane Henrot, Kristof Porteman, Yassine El Hasnaoui, Jan van den Bulcke, Tom de Mil, Cécile Remy, Rachid Cheddadi, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Alain Hambuckers, Franck Trolliet, Marie Dury, Alexandra-Jane Henrot, Kristof Porteman, et al. "Towards a More Realistic Simulation of Plant Species with a Dynamic Vegetation Model Using Field-Measured Traits: The Atlas Cedar, a Case Study." Forests, 2022, 13 (3), pp.446. 10.3390/f13030446 hal-03664015

HAL Id: hal-03664015

https://hal.umontpellier.fr/hal-03664015

Submitted on 10 May 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à 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.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Article

# Towards a More Realistic Simulation of Plant Species with a Dynamic Vegetation Model Using Field-Measured Traits: The Atlas Cedar, a Case Study

Alain Hambuckers <fnref n="1" />,* D, Franck Trolliet <fnref n="2" />, Marie Dury <fnref n="2" />, Alexandra-Jane Henrot <fnref n="2" />, Kristof Porteman <fnref n="1" />, Yassine El Hasnaoui ³D, Jan Van den Bulcke ⁴D, Tom De Mil ⁵D, Cécile C. Remy ⁶D, Rachid Cheddadi ⁷D and Louis François ²D

<figure><img src="image_4.png" /></figure>

check for
updates

**Citation:** Hambuckers, A.; Trolliet, F.; Dury, M.; Henrot, A.-J.; Porteman, K.; El Hasnaoui, Y.; Van den Bulcke, J.; De Mil, T.; Remy, C.C.; Cheddadi, R.; et al. Towards a More Realistic Simulation of Plant Species with a Dynamic Vegetation Model Using Field-Measured Traits: The Atlas Cedar, a Case Study. *Forests* **2022**, *13*, 446. https://doi.org/10.3390/f13030446

Academic Editor: Nick Strigul

Received: 17 January 2022

Accepted: 9 March 2022

Published: 11 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_5.png" /></figure>

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** Improving the model-based predictions of plant species under a projected climate is essential to better conserve our biodiversity. However, the mechanistic link between climatic variation and plant response at the species level remains relatively poorly understood and not accurately developed in Dynamic Vegetation Models (DVMs). We investigated the acclimation to climate of *Cedrus atlantica* (Atlas cedar), an endemic endangered species from northwestern African mountains, in order to improve the ability of a DVM to simulate tree growth under climatic gradients. Our results showed that the specific leaf area, leaf C:N and sapwood C:N vary across the range of the species in relation to climate. Using the model parameterized with the three traits varying with climate could improve the simulated local net primary productivity (NPP) when compared to the model parameterized with fixed traits. Quantifying the influence of climate on traits and including these variations in DVMs could help to better anticipate the consequences of climate change on species dynamics and distributions. Additionally, the simulation with computed traits showed dramatic drops in NPP over the course of the 21st century. This finding is in line with other studies suggesting the decline in the species in the Rif Mountains, owing to increasing water stress.

**Keywords:** dynamic vegetation modelling; specific leaf area; sapwood nitrogen; leaf nitrogen; acclimation; net primary productivity; *Cedrus atlantica* (Endl.) Manetti ex Carrière

## 1. Introduction

Anthropogenic climate change causes range shifts and might lead to the extinction of numerous plant and animal species [1–3]. The subsequent redistribution of species and the alteration of community composition may impact the functioning of whole ecosystems and, therefore, affect the production of ecosystem services and human well-being [4]. Understanding how species are responding to climate change is paramount to help biodiversity conservationists and to anticipate and adapt to the array of socio-ecological consequences

<footnote n="1">¹ Behavioural Biology Unit, UR-SPHERES, University of Liège, 4020 Liège, Belgium; porteman.kristof@gmail.com</footnote>

<footnote n="2">² Unit for Modelling of Climate and Biogeochemical Cycles, UR-SPHERES, University of Liège, 4000 Liège, Belgium; trollietf@googlemail.com (F.T.); mariedury@hotmail.com (M.D.); henrot_alexandra@yahoo.fr (A.-J.H.); louis.francois@uliege.be (L.F.)</footnote>

<footnote n="3">³ Geophysics and Natural Hazards Laboratory, GEOPAC Research Center, Mohammed V University, Agdal, Rabat BP 703, Morocco; elhasnaoui.yassin@gmail.com</footnote>

<footnote n="4">⁴ UGCT—UGent-Woodlab, Department of Environment, Ghent University, 9000 Gent, Belgium; jan.vandenbulcke@ugent.be</footnote>

<footnote n="5">5 TERRA Teaching and Research Center, Gembloux Agro-Bio Tech, University of Liège, 5030 Gembloux, Belgium; tom.demil@uliege.be</footnote>

<footnote n="6">⁶ Institute of Geography, Augsburg University, 86135 Augsburg, Germany; cecile.remy2@gmail.com</footnote>

<footnote n="7">7 CNRS-UM-IRD, Institut des Sciences de l'Evolution, Université Montpellier, 34095 Montpellier, France; rachid.cheddadi@umontpellier.fr</footnote>

<footnote>* Correspondence: alain.hambuckers@uliege.be</footnote>
of climate change. However, the mechanistic link between climatic variation and plants'
response at the species level remains relatively poorly understood [5].

Global Dynamic Vegetation Models (DVMs) are powerful tools to integrate water and carbon cycles, soil and species ecophysiology, and hence, provide reliable estimates of biomass and primary productivity changes. DVMs have increasingly been used to study the response of vegetation to past climates and projected climate change, e.g., [6–8]. To increase the accuracy of the DVMs and our ability to apply them in the field, we need to use them at the species level and scale, which is limited by the knowledge of the species features [9,10]. DVMs require a set of traits that characterize morphology and physiology of the plant species. These traits describe the exposure thresholds to environmental factors such as water stress or minimum temperature to break dormancy [11] or morphophysiological characteristics such as the stomatal conductance or the root depth [12]. Threshold traits could be extracted from environmental factor maps using species distribution samples [13]. Morpho-physiological traits could be directly quantified by sampling in nature but could also need experiments [14]. An alternative method to obtain traits is based on a Bayesian inversion method from the validation dataset of the model outputs [15,16].

However, traits vary not only between species but may also vary within species in
response to environmental factor exposure [17]. These variations might determine, among
others, different aspects of acclimation of species to local conditions and, hence, have to
be taken into account to accurately infer the effect of ongoing global climate change on
plant species [18,19]. Thus, trait variations are key components for the parametrization of
DVMs [20,21] for evaluating the species' ability to persist within a plant community [22,23].

The western Mediterranean region, especially in northern Africa, is highly sensitive to climate change, as exemplified by several drought events and increasing temperature in recent decades [24–26], and projected future decrease in winter rain fall combined with increase in temperature and evaporation [24]. While the Mediterranean region is viewed among biodiversity hotspots that deserve particular attention for conservation efforts [27], it is only recently that the effects of climate change on some important tree species of the region have been assessed [25, 26, 28–30]. The Atlas cedar (Cedrus atlantica (Endl.) Manetti ex Carrière) is a long-lived majestic tree species, classified as endangered on the IUCN red list [31]. Its extant distribution is limited to mountainous regions of Morocco and Algeria where it is threatened by anthropogenic activities. Semi-nomadic pastoral tradition, widespread across the Rif and Atlas Mountains, continues to generate an intensive grazing pressure that severely prevents natural regeneration [26,32]. Recent increased aridity and temperature led to a decrease in growth rate and an increase in mortality [25, 26, 28, 29]. Notably, the range of the Atlas cedar has shrunk by 75% in the Rif Mountains where populations currently persist in isolated microrefugia [30]. UNESCO has recently designated the Atlas Cedar Biosphere Reserve to ensure the conservation of this ecologically and culturally rich region.

The climatic conditions found across the Moroccan cedar's range show strong varia-
tions. Precipitation and temperature regimes are influenced by air masses from the Atlantic
Ocean, the Mediterranean and the Sahara Desert. Such climatic regimes are ideal to investi-
gate trees' ecophysiological response and acclimation to climate and to test a DVM's ability
to simulate the growth of a mountain tree species under different climatic gradients. Such
a study would be of particular interest because it would allow us to evaluate our capacity
to make realistic prognosis with a DVM, also for the future of threatened species such as
the Atlas cedar and possibly to highlight which actions or research would be needed to
improve the model outputs. Otherwise, it is difficult to obtain validation data for a DVM,
such as net primary productivity (NPP) and traits values, at the same time and for a sample
of locations. Here, we collect biomass, NPP data and three important traits for DVM
parametrization: specific leaf area (SLA), leaf C:N and sapwood C:N of the Atlas cedar.
Those traits are easy to estimate and potentially vary across the range of the species, as
well as the biomass and the NPP. We then test whether the use of trait values estimated
at sampling sites improve the ability of a DVM to simulate local NPP and biomass, and
whether the spatial variations in the three traits would improve the simulated range of the Atlas cedar.

## 2. Materials and Methods

### 2.1. Study Area

We collected the three plant traits of the Atlas cedar (Cedrus atlantica (Endl.) Manetti ex Carrière), in Morocco at three sites of the Rif mountains in March 2017, eight sites of the Middle Atlas Mountains in November 2017 and five sites in November 2018 in both ranges (Figure 1, Table 1). Sites were located at altitudes between 1488 and 2205 m above sea level and encompassed a range of contrasted annual and seasonal regimes of temperatures and precipitations.

<figure><img src="image_6.png" /><figcaption>**Figure 1.** Northern Morocco with the border of Algeria (on the right), extant distribution of the Atlas cedar (light green areas) and sampling locations (blue dots).</figcaption></figure>

<table><caption>**Table 1.** Location, altitude, sampling period, numbers of covered plots and of sampled trees.</caption><thead><tr><th>Site</th><th>Region</th><th>Longitude</th><th>Latitude</th><th>Altitude (m asl)</th><th>Period</th><th>Plot Number</th><th>Tree Number</th></tr></thead><tbody><tr><td>1</td><td>Rif</td><td>-4.9481</td><td>35.0295</td><td>1524</td><td>Spring 2017</td><td>2</td><td>44</td></tr><tr><td>2</td><td>Rif</td><td>-4.7064</td><td>34.9597</td><td>1551</td><td>Spring 2017</td><td>2</td><td>31</td></tr><tr><td>3</td><td>Rif</td><td>-4.3356</td><td>34.8703</td><td>1768</td><td>Spring 2017</td><td>1</td><td>15</td></tr></tbody></table>
<table><caption>**Table 1.** *Cont.*</caption><thead><tr><th>Site</th><th>Region</th><th>Longitude</th><th>Latitude</th><th>Altitude (m asl)</th><th>Period</th><th>Plot Number</th><th>Tree Number</th></tr></thead><tbody><tr><td>4</td><td>Middle Atlas</td><td>-5.2564</td><td>32.9622</td><td>2205</td><td>Autumn 2017</td><td>2</td><td>37</td></tr><tr><td>5</td><td>Middle Atlas</td><td>-5.4316</td><td>32.9430</td><td>1673</td><td>Autumn 2017</td><td>2</td><td>35</td></tr><tr><td>6</td><td>Middle Atlas</td><td>-5.2750</td><td>32.9442</td><td>1960</td><td>Autumn 2017</td><td>2</td><td>56</td></tr><tr><td>7</td><td>Middle Atlas</td><td>-5.4290</td><td>32.9722</td><td>1581</td><td>Autumn 2017</td><td>3</td><td>24</td></tr><tr><td>8</td><td>Middle Atlas</td><td>-5.4767</td><td>32.9735</td><td>1590</td><td>Autumn 2017</td><td>2</td><td>34</td></tr><tr><td>9</td><td>Middle Atlas</td><td>-5.4539</td><td>32.8101</td><td>1488</td><td>Autumn 2017</td><td>8</td><td>15</td></tr><tr><td>10</td><td>Middle Atlas</td><td>-5.3812</td><td>32.9539</td><td>1697</td><td>Autumn 2017</td><td>2</td><td>16</td></tr><tr><td>11</td><td>Middle Atlas</td><td>-5.2317</td><td>32.8888</td><td>1980</td><td>Autumn 2017</td><td>1</td><td>31</td></tr><tr><td>12</td><td>Rif</td><td>-4.4243</td><td>34.8233</td><td>1794</td><td>Autumn 2018</td><td>2</td><td>48</td></tr><tr><td>13</td><td>Rif</td><td>-4.5416</td><td>34.8519</td><td>1820</td><td>Autumn 2018</td><td>1</td><td>37</td></tr><tr><td>14</td><td>Middle Atlas</td><td>-5.0287</td><td>33.5529</td><td>1749</td><td>Autumn 2018</td><td>2</td><td>38</td></tr><tr><td>15</td><td>Middle Atlas</td><td>-5.1120</td><td>33.4972</td><td>1717</td><td>Autumn 2018</td><td>2</td><td>44</td></tr><tr><td>16</td><td>Middle Atlas</td><td>-5.1995</td><td>33.4122</td><td>1755</td><td>Autumn 2018</td><td>2</td><td>41</td></tr></tbody></table>

## 2.2. Morpho-Physiological Trait Estimates

We followed standardized protocols [33]. For SLA and leaf C:N, we collected as far as possible 2 twigs per tree from 10 trees in each plot, selecting twigs exposed to direct sunlight. The initial objective was to sample at least two plots with 20 trees per plot but the number of plots per site varied according to tree density and opportunity (Table 1). To determine sapwood C:N, we took a single branch with a diameter larger than 5 cm (it was not possible to take stem slices as suggested [33]), in each plot. The twigs and the branch were collected in the lowest part of the canopy with a 4 m telescopic handle if necessary. To estimate SLA, we isolated 5 sub-samples of fresh needles from a composite sample from the field samples and took pictures using a flat-bed scanner. We put each sub-sample in a paper envelop for drying in an oven at 70 °C to constant weight and weighed the dry mass. Each scan was analyzed with the ImageJ software (https://imagej.net, accessed on 31 March 2021 [34]). For leaf C:N, we prepared 5 composite sub-samples of needles as for SLA. For sapwood C:N, we prepared sub-samples by cutting chips with a chisel. N concentrations were determined following a standard Kjeldahl method. We considered C concentration as 52.8% in leaves [35] and as 50.1% in sapwood [36].

## 2.3. Climate Data

We used CRU (Climate Research Unit, (http://www.cru.uea.ac.uk/, accessed on 12 May 2021) monthly climatic anomalies (0.5° dataset [37]) combined with WorldClim version 2 monthly climate data for 1971–2000 (30 arc sec dataset [38]) to obtain high-resolution (~1 km) transient climatic inputs over the period 1901–2016. Potential evapotranspiration was computed by the DVM with the Penman formula.

## 2.4. *Traits–Climate Factor Relationships*

Using the coordinates of the sampling sites, we linked the trait values with the climate annual values (temperature, daily temperature range, precipitations, relative humidity, potential evapotranspiration, wind speed, sunshine duration, using the 'raster' R package [39]). We first tested linear mixed models with traits as dependent variables, climate factors as fixed effects and season as a random effect. We tested the significance of the random effect with a likelihood ratio test. Since the random effect was never significant, we used linear models with fixed effects only. To select the best model in terms of the Akaike information criterion (AIC), we conducted an exhaustive search of climate factor effects but limiting the maximal number to be included to 5, using 'glmulti' R package [40]. In a preliminary attempt, we evaluated the linear and quadratic effects and pairwise interactions, and different time frames, the overall means over the 30 last years, only the annual values of 2016 or the means over the 3 previous years. The models with the linear effects and
pairwise interactions produced better fits than models with quadratic effects. Additionally, the climate factor means over the last 3 years, i.e., corresponding to the development of the sampled leaves, gave better fit in terms of AIC than the values of a single year or the means over the last 30 years.

## 2.5. Field Net Primary Productivity Estimates

To obtain growth between 2012 and 2016, we estimated aboveground growth. We fixed belowground growth as 0.3 of the aboveground growth. For aboveground growth, we evaluated (1) wood increments, by forest inventories and (2) the amount of foliage per unit of surface, by dividing the leaf area index (LAI) by SLA, and making the assumptions (1) that its amount was constant over the course of the years and (2) that the renewal time was 3 years. Forest inventories were conducted by establishing, as far as possible, several circular plots with radius varying between 10 and 15.80 m to approximately include 20 cedar trees with DBH ≥ 10 cm. Most of the sites were largely dominated by cedar, except one site (#9) where cedars were mostly scattered among other tree species. For site #9, we kept the surface to 315 m² (10 m radius) and increased the number of plots to approach the target number of specimens. In each plot, we measured tree height ($h_{tot}$) with a clinometer (SILVA 70498, Bromma, Sweden) and the diameter at breast height (DBH) of all the trees with DBH ≥ 10 cm with a forester tape; we also recorded the slope of the plot for area correction. By applying a specific allometric equation for total volume (VolTot, Equation (1)) from [41]:

$$VolTot = \frac{h_{tot} \times c_{130}^2}{4\pi \times \left(1 - \frac{1.3}{h_{tot}}\right)^2} \times \left(a + b \times \frac{\sqrt{c_{130}}}{h_{tot}} + c \times \frac{h_{tot}}{c_{130}}\right) \quad (1)$$

with c130 the circumference computed from DBH, a the constant (0.340), b the robustness (1.756) and c the taper (0.002), and taking into account the proportion of basal area of the Atlas cedar among plot total basal area, we obtained the aboveground wood volumes per unit of surface of each tree in each plot. The estimation of tree growth was based on lateral and height growth. For lateral growth, we cored all trees at breast height with a 5 mm-diameter increment borer (HAGLÖF, Langsle, Sweden) and determined the ring widths using a standard measuring table equipped with a binocular magnifier and recording software. These values allowed us to estimate DBH of the previous years by subtractions of the width of the annual ring. For height growth, we computed a linear relationship between observed heights and the log of the circumferences with the data of each site. This allowed us to compute the annual height increments, using the circumferences obtained thanks to the radial annual increments (Table S1, Figure S1). We estimated the variation of wood density in the rings of a sample of 44 cores for rings up to 50 years old by X-ray-computed tomography. All cores were mounted in a custom-made sample holder [42] and scanned with the Nanowood CT scanner from the UGent Centre for X-ray Tomography (Ghent University, Belgium) [43]: 2500 projections were acquired over 360°, each with an exposure time of 1 s. Density profiles and values per ring were then extracted from the 3D images of the cores using the Densitometry Toolbox [44].

Density was related to the linear and the quadratic effects of log of ring width and tree DBH in a linear mixed model, with tree identity as a random effect (Table S2, Figure S2) using the 'lme4' R package [45]. We then used the model to predict the density of each ring omitting the random effect. The differences in volume between the successive years were multiplied by wood density to obtain the aboveground net primary productivity of wood.

To estimate LAI, we took 13 hemispheric canopy pictures per plot with a SAMYANG (Changwon, Korea) 8 mm f/3.5 CS II Fisheye lens mounted on a NIKON (Groot-Bij-gaarden, Belgium) D3100 camera. The pictures were analyzed with the software Hemisfer [46,47], following the settings of [48]. The amount of foliage per unit of surface was obtained by dividing LAI by SLA.
## 2.6. Vegetation Model Simulations

We used the CARAIB (CARbon Assimilation in the Biosphere) dynamic vegetation model [49] to simulate the growth of the Atlas cedar. The functioning of CARAIB has been widely described [7,11,30]. Our simulations consider only one tree species (the Atlas cedar) in the overstory and C3 herbs in the understory. We used the climate data obtained as described above. Atmospheric CO₂ concentrations are prescribed from time series of global values [50,51]. The soil texture was obtained from the global HWSD v1.21 (30 arc sec) database (FAO/IIASA/ISRIC/ISS-CAS/JRC, Harmonized World Soil Database (version 1.2), 2012). Exposure thresholds were obtained by superimposing the above CRU/WorldClim 2 gridded climatology interpolated at a 30 arc sec resolution with a sample of distribution in natural range of the Atlas cedar, and then selecting the prescribed quantile in the distribution of the climatic variables. The range of the study extended between -6.438° and -3.004° of longitude and between 32.254° N and 35.896° N of latitude.

We conducted spatial simulations over the natural range of the Atlas cedar. Four simulations were carried out, varying the set of traits. The first simulation (*PFT traits*) used conifer PFT values (SLA = 100 sq. cm/g; leaf C:N = 55 g/g; sapwood C:N = 400 g/g [52]). The second one (*Field mean traits*) used the means of the field values determined in this study (SLA = 53.1 sq. m/g; leaf C:N = 46.2 g/g, sapwood C:N = 423 g/g). The third one (*Computed traits*) used spatial extrapolations of the traits obtained by applying the traits-climate factor relationships over the study area but limiting the predicted values to the ranges of the observed values. We also conducted simulations at site scale to retrieve NPP and biomass estimates for each of the sites where field sampling had been carried out using the local estimates (*Site traits*).

The potential distribution of the Atlas cedar simulated with CARAIB is discussed and compared with the modern geographical distribution of the species. As a threshold of presence, we selected the site with the lowest biomass field values, recorded in the most southern site of our sampling, in the Middle Atlas, where only small trees were scattered in other native species forming a low forest.

# 3. Results

We set up statistical relationships between the three morpho-physiological trait values and climate factors. We analyze the differences between the DVM outputs parametrized with standard parametrization, i.e., with the plant functional type values (*PFT traits* simulation), with mean values of field estimates (*Field mean traits*) or with local estimates of the traits derived from the relationships with climate factors (*Computed traits*). We then compare NPP and biomass DVM outputs with the field estimates.

## 3.1. *The Traits and Their Relationships with Climate Factors*

The three traits estimated in the field considerably varied between the sites, particularly those of sapwood C:N (Figure 2, Table S3). SLA values were in the lower range of the values reported in [53] for Gymnosperms but in [54], a nitrogen value in *Cedrus deodora* is given, which combined with mean C content gives leaf C:N of 52.8 g/g, rather close to the mean value of this study. The observed sapwood C:N values were very high compared to the value for temperate Gymnosperms given in [36], but they were in the range of the values compiled for Gymnosperms in [55,56]. The variations are within the ranges found in other studies. SLA inside species variation among 1300 individuals belonging to 383 species accounts for 30% of the total [57]. SLA varies between 10 and 32 sq. cm/g among 10 species [58]. Leaf N variation inside species contributes to 50 % of total variation among 442 samples of 223 species [59] and leaf C:N ranges between 20.2 and 53.4 g/g among 10 species [60]. Additionally, sapwood C:N values in [61] comprise between 84.7 and 1360.8 g/g among 59 species.
<figure><img src="image_7.png" /><figcaption>**Figure 2.** Boxplots of the morpho-physiological traits estimated in field.</figcaption></figure>

Additionally, we found significant relationships between the traits and the climate factors, with R-squared between 0.6551 and 0.9086 (Figure 3). All the tested climate factors appeared at least once in the relationships and most of the effects were pairwise interactions (Table S4).

<figure><img src="image_8.png" /><figcaption>**Figure 3.** Relationships between observed traits (see also Table S3) and predicted values with models shown in Table S4.</figcaption></figure>

The relationships between traits and climate factors allow us to produce successions of maps of trait distribution over the study area for each year (examples in Figure 4). The traits smoothly vary over the course of the years, probably owing to the fact that the climate factors were averaged over three years (Figure 5).
<figure><img src="image_9.png" /><figcaption>**Figure 4.** Examples of predicted distribution of the 3 traits (for 2016 growth year) over the study area with sampling locations (red dots). (A): specific leaf area, (B): leaf C:N, (C): sapwood C:N).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 5.** Evolution in the course of the years of the predicted specific leaf area (SLA), leaf C:N and sapwood C:N, over the pixel of two sampling locations (black: Middle Atlas #9, brown: Rif #12).</figcaption></figure>

## 3.2. Comparisons of the DVM Outputs According to Parametrization

The comparisons of the NPP simulations from 1901 to 2016 produced with *PFT traits*, *field mean traits computed traits* or *site traits* showed contrasted levels and variabilities in time and space (Figure 6 and Videos S1–S3) but, nevertheless, that NPP is predicted higher in the Rif than in the Middle Atlas. *PFT trait* simulations showed the lowest NPP level and variability. *Field mean trait* simulation was intermediate while *computed trait* simulation gave the highest values and variability. These observations were confirmed over the sampling sites, particularly with *site trait* parametrization where the differences between sites were larger (Figure 7). Otherwise, the simulated ranges were much larger than the extant range (Figure 6). Finally, the biomass distribution produced by the model under the four parametrizations reflects the NPPs, but simulated values were much smaller than the estimated ones (Figure 8 and Table S5). The simulations also sometimes showed sudden local drops in biomass (Videos S4–S6).
<figure><img src="image_11.png" /><figcaption>Figure 6. *Cont.*</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 6.** Examples of simulated net primary productivity annual distributions for selected years. Comparisons of simulation outputs using *PFT traits*, *field mean traits* or *computed traits*, with sampling locations (blue squares), and contours of the extant distribution (red lines).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 7.** Simulated net primary productivity (NPP) between 1900 and 2016 over the sampling locations (Middle Atlas in blue colors, Rif in reds) using *PFT traits*, *field mean traits*, *computed traits* and *site traits*.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 8.** Predicted biomass in 2016. Comparisons of simulation outputs using *PFT traits*, *field mean traits* or *computed traits* with sampling locations (blue squares), and contours of the extant distribution (red lines). The threshold of 11.7 kg C/sq. m corresponding to the lowest recorded field value (site #9) was only reached in the Rif with the field and the *computed trait* simulations.</figcaption></figure>
3.3. Comparisons of DVM Outputs with Field Estimates

The field estimations of NPP considerably varied between the sites (Table S5). The comparison of the simulated NPP values with the field estimates revealed the poor response of the model parametrized with *PFT traits*. The simulations seemed improved with *field mean traits* and more with *computed traits* or *site traits* (Figure 9). However, this was only partially confirmed by the performance indicators. For instance, the highest R-squared value was obtained with the *computed trait* simulation for 2016 but the root mean square and the mean bias were also higher (Figure 10). Otherwise, the relationships between simulated and field values of biomass were insignificant.

<figure><img src="image_15.png" /><figcaption>Figure 9. *Cont.*</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 9.** Comparison of net primary productivity (g C/sq. m/yr) model outputs for years 2012 to 2016 with sites estimates, when using *PFT traits*, *field mean traits*, *computed traits* or *site traits*.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 10.** *Cont.*</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 10.** Performance indicators of net primary productivity model outputs for years 2012 to 2016 (comparisons with field estimates), when using *PFT traits*, *field mean traits*, *computed traits* or *site traits*: R-squared, root mean squared errors (RMSE, g C/sq. m/yr) and mean biases (g C/sq. m/yr).</figcaption></figure>

**4. Discussion**

In accordance with our hypothesis, the estimations of the three selected traits showed large variations between the sampling sites (Figure 2), notably for sapwood C:N. We identify significant but complex relationships between the traits and the climate factors over the 3 previous years (Figure 3, Table S4). These relationships were better than those obtained with 30-year climate means or with current year values. They allowed us to produce smooth annual trait maps, which were used in the *computed trait* simulation (Figures 4 and 5). NPP estimated at the field sampling sites also considerably varied (by a factor of 170, Table S5). The four simulations performed with CARAIB (*PFT traits*, *field mean traits* or *computed traits*) produced comparable ranges but were dramatically larger than the extant range. The use of the *computed traits* increased the range of the temporal and the spatial responses of the model (Figures 6–8). Comparisons of the four simulations with the field-estimated NPP suggest an improved response of the model parametrized with *computed traits* or *site traits*, but this was not confirmed by the statistics (Figures 9 and 10).

The differences in NPP estimations between the sampling sites have been explored. Differences in site fertility may often play a role because the Atlas cedar occurs in pure or in quasi-monospecific stands. At margins when competition with other species increases, site fertility probably plays a minor role [62]. The Atlas cedar growth responds to temperature, precipitation and drought stress [29] and the capacity for regeneration relies on winter temperature, which affects the germination rate [63]. We cannot exclude some estimation deviations related to insufficient sampling, uneven radial growth, distribution of wood density, and estimation of the height. Nevertheless, DBH responds well to climate variation [64]. Also, aboveground biomass estimates obtained with allometric equations Combined with tree ring data give similar results to successive estimations in permanent plots [65]. In addition, we varied wood density with ring width and DBH, since this factor is recognized as more important than the selection of allometric equation to estimate carbon immobilization [64]. Using the allometric equations of [35] rather than those provided by [41] has a minor impact on the results. In addition, most of our estimations were in the range of the extrapolated standing wood biomass produced over one century by [66] in Algerian cedar forests. Using estimations based on [35] to obtain branches, leaves and belowground NPP, we found for the Algerian Atlas cedar stands an NPP between 18.5 and 605 g C/sq. m/yr. In our study, the most productive sites have high biomass and trees with
extremely large rings, while those sites with the lowest productivity have very small and scattered trees, smaller than the average with finer rings. The latter may not be considered as cedar forests.

Our study shows that the three traits varied spatially and that these variations may be related to recent potential acclimation of the Atlas cedar to the ongoing climate change in the Mediterranean. SLA increases with temperature and water availability and decreases with light intensity and duration [53]. In evergreen conifers, N is stored in the youngest needles as photosynthetic proteins [67] and leaf N decreases with increasing temperature, drought and irradiance [19,68–70]. Leaf N is positively correlated with wood N [71] but N is not stored in the wood of conifers [67]. Sapwood N does not seem related to drought or temperature [72] and its relationship to soil N is not established [71,73], probably because of its low concentration. For instance, relationship of N sapwood with soil N only appears in the outer stem after N fertilization [74]. Combined with the proportion of sapwood to heartwood, which positively scales radial growth to satisfy the conductive needs and adjusts to tree height and leaf mass [75,76], sapwood N is nevertheless an important component to compute respiratory cost.

We found relationships between traits and climate factors (Figure 3). A single climate factor or aggregated variables have little predictive power, which may explain the data dispersion, e.g., [77]. We assume that the selected models have some predictive power in the range of the observed climate. However, the studied sites do not cover the whole range of conditions over the simulated area and years, which probably limit their relevance.

SLA is considered as a major factor for plant growth rate, although its relationship to maximized photosynthetic function is not well understood [33,57,78]. As SLA increases, the number of photosynthetic cells per unit of plant weight increases too, which boosts photosynthetic rate and carbon uptake. Additionally, high growth rate is reached at the expense of structural tissues and stress protection and notably water stress. SLA triggers N leaf concentration, when N is available [68]. As SLA, leaf N (also leaf C:N) is a major morpho-physiological trait favoring photosynthesis assimilation because it allows the synthesis of enzymes and particularly of Rubisco, but it is also a driver of dark respiration [79]. Sapwood N (also sapwood C:N) is another important trait in trees because it reflects the importance of wood parenchyma [80] and, hence, of wood respiratory costs, which might be significant in the carbon budget, even during the dormant season [81]. The main factors controlling tree performance are the availability of light, water and nutrients; however, performances also rely on the capacity of trees to collect these resources and their efficiency in their use. These features are mediated by intraspecific traits. It is concluded that individual traits are poor predictors of individual performance variation [57,82] because many traits affect more than a single aspect of the physiology and that acclimation, like adaptation, would finally be a combination of opposite requirements. In the DVM CARAIB, SLA allows the conversion of the assimilated carbon into leaf area index [16], which is a fundamental factor of the canopy sub-model computing light assimilation [49]. The leaf C:N ratio controls the dark respiration rate and the maximum carboxylation velocity in the Farquhar et al. photosynthesis model [83], and the sapwood C:N controls sapwood respiration [49]. A major feature of the DVMs is their ability to integrate ecophysiological variables that have opposite effects.

Simulated NPP was strongly influenced by the three traits values as previously observed [9]. The annual and spatial variability notably increased when we used the *field-mean*, the *computed* or the *site traits*. However, the gain in performances of the model was modest. Process-based models lack the ability to reproduce the full range of productivity variability of natural or artificial plant systems. For instance, the simulation of wheat yield for several years across a European gradient indicated the low performances of a set of process-based models [84]. Comparisons of gross primary productivity outputs of several DVMs at a global scale showed low ability to reproduce interannual variations [85]. CARAIB simulating vegetation over Western Europe responded to the strong climate gradients but tended to underestimate the highest NPP values and overestimate the lowest ones [11].
When several PFTs are simulated together over an area, the traits differences between plant types allow geographic turnover in the simulation, mimicking trait acclimation. This is the reason why the DVMs show reasonable agreement with observational data. On the contrary, when a single species is simulated, there is only one trait set and no acclimation mimicry is possible. Adjusting traits (height, SLA, leaf C:N, root depth) to species may markedly modulate the simulated NPP [9]. Additionally, adjusting traits of plant types (SLA, leaf nitrogen and phosphorus, and wood density) to local conditions provided more accurate model predictions in a tropical altitudinal gradient only thanks to solar radiation [86,87]. Beyond SLA and N concentrations, the DVMs may also benefit from improved parametrization for other traits. Notably, stomatal conductance related with water stress is an important trait because it controls the internal CO₂ concentration and photosynthesis dark reaction. Plants also acclimatize to water availability by adjustment of the root to shoot ratio [88,89], of stomatal regulation sensitivity [90] and of leaf features, such as trichomes, epicuticular wax, cuticle, hypoderm, and stomata density [91–97]. Thus, in our *computed* and *site trait* simulations, we assume that the increased variations of NPP may be balanced by variations of other traits.

Concerning the range prediction in this study, it is worth noting that four simulations produced comparable results in terms of NPP and biomass (Figures 6 and 8), the differences being mainly in the levels and the variations, but not in the distributions. This is because of the constraints exerted by the bioclimatic thresholds. With the minimal NPP value (20 g C/sq. m/yr observed at site #9), the simulated potential ranges are much larger than the modern range of the Atlas cedar, which confirms that CARAIB predicts the fundamental climate niche because it is constrained by climate extremes [12]. To compute the realized niche and to delineate the actual range, the model needs to take into account the intensity of the main biotic interactions. DVMs are able to partially simulate interspecific competition, i.e., competition for light and water, but this requires correct parametrization for each species, which leads us back to the main issue of this study, while natural grazing, parasitism and the pressure generated by anthropic activities require the coupling of the DVMs with other sources of information.

**5. Conclusions**

In this study, we show that traits vary across the range of the Atlas cedar in relation-
ship with climate, that their use improves the DVM simulations of local NPP, and that
adjusting these traits with climate could lead to better performances. However, the benefit
of predicting traits with multiple regressions driven by climate factors was modest and
a more sophisticated approach such as a more mechanistic approach for each trait with
climate factors would be desirable. Notwithstanding the limits of our results, it should
also be underlined that the levels obtained for simulated NPP are in the correct range of
the observed values despite the absence of any calibration phase. Thus, quantifying the
influence of climate on traits and including these variations in DVMs may be of paramount
importance for evaluating the consequences of global climate change on species dynamics
and plant community persistence. Besides SLA and leaf and sapwood C:N, additional
traits concerning physiological processes such as stomatal regulation, or elements of plant
structure such as the root to shoot ratio, should be adjusted or modelled in the DVMs to
allow for more reliable simulations.

So far, our best simulation shows the Atlas cedar NPP collapses over several years of
the 21st century in the Rif mountains (Figure 7, *computed traits*), which has to be viewed in
parallel with a probable reduction in available soil water and field observed episodes of
decline due to water stress [30]. Even if these dramatic modelling results are perfectible, they
are nevertheless congruent. However, most of all, they also concern the most productive
stands that we encountered and what is considered to be modern microrefugia, essential
for the future conservation of the species. This is disturbing because it suggests that even
these scattered populations would not be sheltered from climate change.
**Supplementary Materials:** The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/f13030446/s1, Figure S1: Relationships between *Cedrus atlantica* circumference at breast height and height per site; Figure S2: Relationship between ring width and wood density of 44 cores of *Cedrus atlantica*; Table S1: Linear models of *Cedrus atlantica* height as a function of log of circumference for each site; Table S2: Linear mixed model of *Cedrus atlantica* wood density as a function of linear and quadratic effects of the log of ring width and tree diameter at breast height as fixed effects, and tree identity as random effect; Table S3: Traits estimated for each site; Table S4: Linear models of three traits (specific leaf area, leaf C:N and sapwood C:N) as a function of climate factors (means over 2014, 2015 and 2016 of precipitation: prc, of daily temperature range: dte, of temperature: tem, of sunshine hours: shr, wind speed: wnd, of potential evapotranspiration: pet, of relative humidity: rhu) effects (linear effects and pairwise interactions). Best models (with max-imum 5 terms for sample sizes = 16) after exhaustive searches and selections on the basis of AIC and F tests; Table S5: Estimations of mean biomass for 2013 to 2017 (kg C/sq. m) and of net primary productivity (NPP) for years 2012 to 2016 (g C/sq. m/yr); Videos S1–S3: Simulated net primary productivity over the study area between 1900 and 2016 with *PFT traits*, *Field mean traits* or *Computed traits*; Videos S4–S6: Simulated biomass over the study area between 1900 and 2016 with *PFT traits*, *Field mean traits* or *Computed traits*.

**Author Contributions:** Conceptualization, A.H. and L.F.; methodology, A.H. and L.F.; data acquisition, A.H., L.F., F.T., R.C. and Y.E.H.; data analyses, A.H., K.P., F.T., J.V.d.B. and T.D.M.; DVM modelling, L.F., M.D. and A.-J.H.; original draft preparation, A.H. and F.T.; visualization, AH. and R.C.; manuscript revision, A.H., R.C., C.C.R., J.V.d.B., and T.D.M.; funding acquisition, R.C., L.F. and A.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the F.RS.-FNRS under research grant FRS-FNRS X.3041.17—VULPES-ULg. It is a contribution to the Belmont Forum project VULPES ANR-15-MASC-0003. BOF Special Research Fund support the UGCT Center of Expertise (BOF.EXP.2017.0007) and J.V.d.B. (BOF Starting Grant, BOFSTG2018000701).

**Acknowledgments:** We are grateful to the Department of Water and Forests of Morocco which allowed the sampling of the trees and provided assistance for field operations. The authors thank an anonymous referee for helpful comments that improved the paper.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

## References

1. Walther, G.R.; Post, E.; Convey, P.; Menzel, A.; Parmesan, C.; Beebee, T.J.C.; Fromentin, J.M.; Hoegh-Guldberg, O.; Bairlein, F. Ecological responses to recent climate change. *Nature* **2002**, *416*, 389–395. [CrossRef] [PubMed]

2. Thomas, C.D.; Cameron, A.; Green, R.E.; Bakkenes, M.; Beaumont, L.J.; Collingham, Y.C.; Erasmus, B.F.N.; Ferreira De Siqueira, M.; Grainger, A.; Hannah, L.; et al. Extinction risk from climate change. *Nature* **2004**, *427*, 145–148. [CrossRef] [PubMed]

3. Steinbauer, M.J.; Grytnes, J.-A.; Jurasinski, G.; Kulonen, A.; Lenoir, J.; Pauli, H.; Rixen, C.; Winkler, M.; Bardy-Durchhalter, M.; Barni, E.; et al. Accelerated increase in plant species richness on mountain summits is linked to warming. *Nature* **2018**, *556*, 231–234. [CrossRef]

4. Pecl, G.T.; Araújo, M.B.; Bell, J.D.; Blanchard, J.; Bonebrake, T.C.; Chen, I.C.; Clark, T.D.; Colwell, R.K.; Danielsen, F.; Evengård, B.; et al. Biodiversity redistribution under climate change: Impacts on ecosystems and human well-being. *Science* **2017**, *355*, eaai9214. [CrossRef]

5. Meng, H.; Wei, X.; Franklin, S.B.; Wu, H.; Jiang, M. Geographical variation and the role of climate in leaf traits of a relict tree species across its distribution in China. *Plant Biol.* **2017**, *19*, 552–561. [CrossRef] [PubMed]

6. Dai, Z.; Johnson, K.D.; Birdsey, R.A.; Hernandez-Stefanoni, J.L.; Dupuy, J.M. Assessing the effect of climate change on carbon sequestration in a Mexican dry forest in the Yucatan Peninsula. *Ecol. Complex.* **2015**, *24*, 46–56. [CrossRef]

7. Henrot, A.J.; Utescher, T.; Erdei, B.; Dury, M.; Hamon, N.; Ramstein, G.; Krapp, M.; Herold, N.; Goldner, A.; Favre, E.; et al. Middle Miocene climate and vegetation models and their validation with proxy data. *Palaeogeogr. Palaeocl.* **2017**, *467*, 95–119. [CrossRef]

8. Wu, Z.; Ahlström, A.; Smith, B.; Ardö, J.; Eklundh, L.; Fensholt, R.; Lehsten, V. Climate data induced uncertainty in model-based estimations of terrestrial primary productivity. *Environ. Res. Lett.* **2017**, *12*, 064013. [CrossRef]

9. Dury, M.; Mertens, L.; Fayolle, A.; Verbeeck, H.; Hambuckers, A.; François, L. Refining species traits in a dynamic vegetation model to project the impacts of climate change on tropical trees in Central Africa. *Forests* **2018**, *9*, 722. [CrossRef]
10. Maréchaux, I.; Langerwisch, F.; Huth, A.; Bugmann, H.; Morin, X.; Reyer, C.P.O.; Seidl, R.; Collalti, A.; Dantas de Paula, M.; Fischer, R.; et al. Tackling unresolved questions in forest ecology: The past and future role of simulation models. *Ecol. Evol.* **2021**, *11*, 3746–3770. [CrossRef]

11. Dury, M.; Hambuckers, A.; Warnant, P.; Henrot, A.; Favre, E.; Ouberdous, M.; François, L. Responses of European forest ecosystems to 21 st century climate: Assessing changes in interannual variability and fire intensity. *Iforest* **2011**, *4*, 82–99. [CrossRef]

12. Raghunathan, N.; François, L.; Dury, M.; Hambuckers, A. Contrasting climate risks predicted by dynamic vegetation and ecological niche-based models applied to tree species in the Brazilian Atlantic Forest. *Reg. Environ. Change* **2019**, *19*, 219–232. [CrossRef]

13. Schneider, F.D.; Morsdorf, F.; Schmid, B.; Petchey, O.L.; Hueni, A.; Schimel, D.S.; Schaepman, M.E. Mapping functional diversity from remotely sensed morphological and physiological forest traits. *Nat. Commun.* **2017**, *8*, 1441. [CrossRef] [PubMed]

14. Carmona, P.; Rota, C.; Azcárate, F.M.; Peco, B. More for less: Sampling strategies of plant functional traits across local environmental gradients. * Funct. Ecol.* **2015**, *29*, 579–588. [CrossRef]

15. Hartig, F.; Dyke, J.; Hickler, T.; Higgins, S.I.; O'Hara, R.B.; Scheiter, S.; Huth, A. Connecting dynamic vegetation models to data—An inverse perspective. *J. Biogeogr.* **2012**, *39*, 2240–2252. [CrossRef]

16. Minet, J.; Laloy, E.; Tychon, B.; François, L. Bayesian inversions of a dynamic vegetation model at four European grassland sites. *Biogeosciences* **2015**, *12*, 2809–2829. [CrossRef]

17. Shipley, B.; De Bello, F.; Cornelissen, J.H.C.; Laliberté, E.; Laughlin, D.C.; Reich, P.B. Reinforcing loose foundation stones in trait-based plant ecology. *Oecologia* **2016**, *180*, 923–931. [CrossRef] [PubMed]

18. Heilmeier, H. Functional traits explaining plant responses to past and future climate changes. *Flora* **2019**, *254*, 1–11. [CrossRef]

19. Dalle Fratte, M.; Brusa, G.; Pierce, S.; Zanzottera, M.; Cerabolini, B.E.L. Plant trait variation along environmental indicators to infer global change impacts. *Flora* **2019**, *254*, 113–121. [CrossRef]

20. Moran, E.V.; Hartig, F.; Bell, D.M. Intraspecific trait variation across scales: Implications for understanding global change responses. *Glob. Change Biol.* **2016**, *22*, 137–150. [CrossRef]

21. Berzaghi, F.; Wright, I.J.; Kramer, K.; Oddou-Muratorio, S.; Bohn, F.J.; Reyer, C.P.O.; Sabaté, S.; Sanders, T.G.M.; Hartig, F. Towards a new generation of trait-flexible vegetation models. *Trends. Ecol. Evol.* **2020**, *35*, 191–205. [CrossRef] [PubMed]

22. Bloomfield, K.J.; Cernusak, L.A.; Eamus, D.; Ellsworth, D.S.; Colin Prentice, I.; Wright, I.J.; Boer, M.M.; Bradford, M.G.; Cale, P.; Cleverly, J.; et al. A continental-scale assessment of variability in leaf traits: Within species, across sites and between seasons. *Funct. Ecol.* **2018**, *32*, 1492–1506. [CrossRef]

23. Gillison, A.N. Plant functional indicators of vegetation response to climate change, past present and future: I. Trends, emerging hypotheses and plant functional modality. *Flora* **2019**, *254*, 12–30. [CrossRef]

24. Touchan, R.; Anchukaitis, K.J.; Meko, D.M.; Kerchouche, D.; Slimani, S.; Ilmen, R.; Hasnaoui, F.; Guibal, F.; Julio Camarero, J.; Sánchez-Salguero, R.; et al. Climate controls on tree growth in the Western Mediterranean. *Holocene* **2017**, *27*, 1429–1442. [CrossRef]

25. Kherchouche, D.; Kalla, M.; Gutierrez, E.; Briki, A.; Hamchi, A. Drought and dieback of atlas cedar (*Cedrus atlantica* Manetti) in the Belezma mountains (Algeria). *Sécheresse* **2013**, *24*, 129–137. [CrossRef]

26. Slimani, S.; Derridj, A.; Gutierrez, E. Ecological response of *Cedrus atlantica* to climate variability in the massif of Guetiane (Algeria). *For. Syst.* **2014**, *23*, 448–460. [CrossRef]

27. Myers, N.; Mittermeler, R.A.; Mittermeler, C.G.; Da Fonseca, G.A.B.; Kent, J. Biodiversity hotspots for conservation priorities. *Nature* **2000**, *403*, 853–858. [CrossRef]

28. Linares, J.C.; Taïqui, L.; Camarero, J.J. Increasing drought sensitivity and decline of atlas cedar (*Cedrus atlantica*) in the Moroccan Middle Atlas forests. *Forests* **2011**, *2*, 777–796. [CrossRef]

29. Linares, J.C.; Taïqui, L.; Sangüesa-Barreda, G.; Seco, J.I.; Camarero, J.J. Age-related drought sensitivity of Atlas cedar (*Cedrus atlantica*) in the Moroccan Middle Atlas forests. *Dendrochronologia* **2013**, *31*, 88–96. [CrossRef]

30. Cheddadi, R.; Henrot, A.; Francois, L.; Boyer, F.; Bush, M.; Carre, M.; Coissac, E.; De Oliveira, P.; Ficetola, F.; Hambuckers, A.; et al. The role of modern microrefugia in preserving mountain tree species: *Cedrus atlantica* in Morocco. *Front. Ecol. Environ.* **2017**, *5*, 114. [CrossRef]

31. Thomas, P. *Cedrus atlantica*. The IUCN Red List of Threatened Species 2013: e.T42303A2970716. Available online: https://www.iucnredlist.org/species/42303/2970716 (accessed on 10 December 2021).

32. Coudel, M.; Aubert, P.-M.; Aderghal, M.; Hély, C. Pastoral and woodcutting activities drive *Cedrus atlantica* Mediterranean forest structure in the Moroccan Middle Atlas. *Ecol. Appl.* **2016**, *26*, 574–586. [CrossRef] [PubMed]

33. Cornelissen, J.H.C.; Lavorel, S.; Garnier, E.; Diaz, S.; Buchmann, N.; Gurvich, D.E.; Reich, P.B.; ter Steege, H.; Morgan, H.D.; Van der Heijden, M.G.A.; et al. A handbook of protocols for standardised and easy measurement of plant functional traits worldwide. *Austral. J. Bot.* **2003**, *51*, 335–380. [CrossRef]

34. Schneider, C.A.; Rasband, W.S.; Eliceiri, K.W. NIH Image to ImageJ: 25 years of image analysis. *Nat. Methods* **2012**, *9*, 671–675. [CrossRef] [PubMed]

35. Durkaya, B.; Durkaya, A.; Makineci, E.; Ülküdür, M. Estimation of above-ground biomass and sequestered carbon of Taurus Cedar (*Cedrus libani* L.) in Antalya, Turkey. *Iforest* **2013**, *6*, 278–284. [CrossRef]
36. Martin, A.R.; Gezahegn, S.; Thomas, S.C. Variation in carbon and nitrogen concentration among major woody tissue types in temperate trees. *Can. J. For. Res.* **2015**, *45*, 744–757. [CrossRef]

37. Harris, I.; Jones, P.D.; Osborn, T.J.; Lister, D.H. Updated high-resolution grids of monthly climatic observations—The CRU TS3.10 Dataset. *Int. J. Climatol.* **2014**, *34*, 623–642. [CrossRef]

38. Fick, S.E.; Hijmans, R.J. WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. *Int. J. Climatol.* **2017**, *37*, 4302–4315. [CrossRef]

39. Hijmans, R.J.; van Etten, J. Raster: Geographic data Analysis and modeling (R Package). Available online: https://CRAN.R-project.org/package=raster (accessed on 1 April 2020).

40. Calcagno, V.; de Mazancourt, C. glmulti: An R package for easy automated model selection with (generalized) linear models. *J. Stat. Softw.* **2010**, *34*, 1–29. [CrossRef]

41. Deleuze, C.; Morneau, F.; Renaud, J.-P.; Vivien, Y.; Rivoire, M.; Santenoise, P.; Longuetaud, F.; Mothe, F.; Hervé, J.C.; Vallet, P. Estimer le volume total d'un arbre, quelles que soient l'essence, la taille, la sylviculture, la station. *Rendez-Vous Tech.* **2014**, *44*, 22–32. Available online: https://www.onf.fr/outils/ressources/d805d212-4584-4a8d-ae3a-c3a896a386a4/++versions++/1/++paras++/2/++ass++/1/++i18n++data:fr?_=_1543929504.543255&download=1 (accessed on 22 February 2019).

42. Van den Bulcke, J.; Wernersson, E.L.G.; Dierick, M.; Van Loo, D.; Masschaele, B.; Brabant, L.; Boone, M.N.; Van Hoorebeke, L.; Haneca, K.; Brun, A.; et al. 3D tree-ring analysis using helical X-ray tomography. *Dendrochronologia* **2014**, *32*, 39–46. [CrossRef]

43. Dierick, M.; Van Loo, D.; Masschaele, B.; Van Den Bulcke, J.; Van Acker, J.; Cnudde, V.; Van Hoorebeke, L. Recent micro-CT scanner developments at UGCT. *Nucl. Instrum. Methods Phys. Res. B* **2014**, *324*, 35–40. [CrossRef]

44. De Mil, T.; Vannoppen, A.; Beeckman, H.; Van Acker, J.; Van den Bulcke, J. A field-to-desktop toolchain for X-ray CT densitometry enables tree ring analysis. *Ann. Bot.* **2016**, *117*, 1187–1196. [CrossRef] [PubMed]

45. Bates, D.; Mächler, M.; Bolker, B.; Walker, S. Fitting Linear Mixed-Effects Models Using lme4. *J. Stat. Softw.* **2015**, *67*, 48. [CrossRef]

46. Schleppi, P.; Conedera, M.; Sedivy, I.; Thimonier, A. Correcting non-linearity and slope effects in the estimation of the leaf area index of forests from hemispherical photographs. *Agric. Forest Meteorol.* **2007**, *144*, 236–242. [CrossRef]

47. Thimonier, A.; Sedivy, I.; Schleppi, P. Estimating leaf area index in different types of mature forest stands in Switzerland: A comparison of methods. *Eur. J. For. Res.* **2010**, *129*, 543–562. [CrossRef]

48. Bequet, R.; Kint, V.; Campioli, M.; Vansteenkiste, D.; Muys, B.; Ceulemans, R. Influence of stand, site, and meteorological variables on the maximum leaf area index of beech, oak and Scots pine. *Eur. J. For. Res.* **2012**, *131*, 283–295. [CrossRef]

49. Warnant, P.; François, L.; Strivay, D.; Gérard, J.-C. CARAIB: A global model of terrestrial biological productivity. *Global Biogeochem. Cycles* **1994**, *8*, 255–270. [CrossRef]

50. Meinshausen, M.; Smith, S.J.; Calvin, K.; Daniel, J.S.; Kainuma, M.L.T.; Lamarque, J.; Matsumoto, K.; Montzka, S.A.; Raper, S.C.B.; Riahi, K.; et al. The RCP greenhouse gas concentrations and their extensions from 1765 to 2300. *Clim. Change* **2011**, *109*, 213–241. [CrossRef]

51. Dlugokencky, E.J.; Lang, P.M.; Masarie, K.A.; Crotwell, A.M.; Crotwell, M.J. Atmospheric Carbon Dioxide Dry Air Mole Fractions from the NOAA. In ESRL Carbon Cycle Cooperative Global Air Sampling Network. Available online: ftp://aftp.cmdl.noaa.gov/data/trace_gases/co2/flask/surface/ (accessed on 1 April 2020).

52. Popova, S.; Utescher, T.; Gromyko, D.V.; Mosbrugger, V.; Herzog, E.; François, L. Vegetation change in Siberia and the northeast of Russia during the cenozoic cooling: A study based on diversity of plant functional types. *Palaios* **2013**, *28*, 418–432. [CrossRef]

53. Poorter, H.; Niinemets, U.; Poorter, L.; Wright, I.J.; Villar, R. Causes and consequences of variation in leaf mass per area (LMA): A meta-analysis. *New. Phytol.* **2009**, *182*, 565–588. [CrossRef]

54. Perry, E.; Hickman, G.W. A survey to determine the leaf nitrogen concentrations of 25 landscape tree species. *J. Arboric.* **2001**, *27*, 152–159. [CrossRef]

55. Meerts, P. Mineral nutrient concentrations in sapwood and heartwood: A literature review. *Ann. For. Sci.* **2002**, *59*, 713–722. [CrossRef]

56. Reich, P.B.; Tjoelker, M.G.; Pregitzer, K.S.; Wright, I.J.; Oleksyn, J.; Machado, J.L. Scaling of respiration to nitrogen in leaves, stems and roots of higher land plants. *Ecol. Lett.* **2008**, *11*, 793–801. [CrossRef] [PubMed]

57. Poorter, L.; Castilho, C.V.; Schietti, J.; Oliveira, R.S.; Costa, F.R.C. Can traits predict individual growth performance? A test in a hyperdiverse tropical forest. *New. Phytol.* **2018**, *219*, 109–121. [CrossRef]

58. Verbeeck, H.; Betehndoh, E.; Maes, W.H.; Hubau, W.; Kearsley, E.; Buggenhout, L.; Hufkens, K.; Huygens, D.; Van Acker, J.; Beeckman, H.; et al. Functional leaf trait diversity of 10 tree species in congolese secondary tropical forest. *J. Trop. For. Sci.* **2014**, *26*, 409–419.

59. Dong, N.; Colin Prentice, I.; Evans, B.J.; Caddy-Retalic, S.; Lowe, A.J.; Wright, I.J. Leaf nitrogen from first principles: Field evidence for adaptive variation with climate. *Biogeosciences* **2017**, *14*, 481–495. [CrossRef]

60. Luo, T.; Luo, J.; Pan, Y. Leaf traits and associated ecosystem characteristics across subtropical and timberline forests in the Gongga Mountains, Eastern Tibetan Plateau. *Oecologia* **2005**, *142*, 261–273. [CrossRef] [PubMed]

61. Martin, A.R.; Erickson, D.L.; Kress, W.J.; Thomas, S.C. Wood nitrogen concentrations in tropical trees: Phylogenetic patterns and ecological correlates. *New. Phytol.* **2014**, *204*, 484–495. [CrossRef]

62. Guedira, A.; Satrani, B.; Douira, A. Effet de la densité de *Cytisus grandiflorus* (Brot.) DC. sur la fertilité du sol en azote et son influence sur la croissance d'une plantation de cèdre (*Cedrus atlantica* Manetti). *Acta Bot. Gall.* **2008**, *155*, 375–387. [CrossRef]
63. Takos, I.; Merou, T. Effect of storage conditions and seed treatment on germination of *Cedrus deodara* Loud. and C. libani A. Rich. *Silvae Genet.* **2001**, *50*, 205–208. [CrossRef]

64. Bouriaud, O.; Bréda, N.; Dupouey, J.L.; Granier, A. Is ring width a reliable proxy for stem-biomass increment? A case study in European beech. *Can. J. For. Res.* **2005**, *35*, 2920–2933. [CrossRef]

65. Dye, A.; Plotkin, A.B.; Bishop, D.; Pederson, N.; Poulter, B.; Hessl, A. Comparing tree-ring and permanent plot estimates of aboveground net primary production in three eastern U.S. forests. *Ecosphere* **2016**, *7*, e01454. [CrossRef]

66. Bentouati, A.; Oujdehih, B. Première étude de la croissance et de la productivité du Cèdre de l'Atlas (*Cedrus atlantica* Manetti) dans le massif de Bélezma (Aurés—Algérie). *For. Méditerr.* **1999**, *20*, 115–119. Available online: http://www.foret-mediterraneenne.org/upload/biblio/FORET_MED_1999_2_115-119.pdf (accessed on 23 September 2016).

67. Millard, P.; Grelet, G.A. Nitrogen storage and remobilization by trees: Ecophysiological relevance in a changing world. *Tree Physiol.* **2010**, *30*, 1083–1095. [CrossRef] [PubMed]

68. Meziane, D.; Shipley, B. Direct and indirect relationships between specific leaf area, leaf nitrogen and leaf gas exchange. Effects of irradiance and nutrient supply. *Ann. Bot.* **2001**, *88*, 915–927. [CrossRef]

69. Reich, P.B.; Oleksyn, J. Global patterns of plant leaf N and P in relation to temperature and latitude. *Proc. Natl. Acad. Sci. USA* **2004**, *101*, 11001–11006. [CrossRef] [PubMed]

70. Chen, Y.; Han, W.; Tang, L.; Tang, Z.; Fang, J. Leaf nitrogen and phosphorus concentrations of woody plants differ in responses to climate, soil and plant growth form. *Ecography* **2013**, *36*, 178–184. [CrossRef]

71. Heineman, K.D.; Turner, B.L.; Dalling, J.W. Variation in wood nutrients along a tropical soil fertility gradient. *New. Phytol.* **2016**, *211*, 440–454. [CrossRef]

72. Cui, Z.Y.; Yang, Z.J.; Xu, D.P.; Xi, R.C.; Zhang, N.N.; Liu, X.J.; Hong, Z. Stem respiration and chemical composition in *Dalbergia odorifera* plantations differing in soil moisture content. *Aus. J. For. Sci.* **2017**, *134*, 347–365.

73. Kotowska, M.M.; Wright, I.J.; Westoby, M. Parenchyma abundance in wood of evergreen trees varies independently of nutrients. *Front. Plant Sci.* **2020**, *11*, 86. [CrossRef]

74. Alcubilla, M.; Aufsess, H.V.; Heibl, R.; Rehfuess, K.E. Stickstoffdüngungsversuche in einer Fichtenwuchsstockung (*Picea abies* Karst.) auf devastierter Kalkmergel-Rendzina: III. Einflüsse von Ernährungszustand und Inokulation auf die Inhaltsstoffe von Bast und Holz. *Forstwiss. Centralbl.* **1988**, *107*, 205–218. [CrossRef]

75. Koch, G.W.; Fredeen, A.L. Transport challenge in tall trees. In *Vascular Transport in Plants*; Holbrook, N.M., Zwieniecki, M.A., Eds.; Academic Press: Cambridge, MA, USA, 2005. [CrossRef]

76. Morais, M.C.; Pereira, H. Heartwood and sapwood variation in *Eucalyptus globulus* Labill. trees at the end of rotation for pulpwood production. *Ann. For. Sci.* **2007**, *64*, 665–671. [CrossRef]

77. Steane, D.A.; McLean, E.H.; Potts, B.M.; Prober, S.M.; Stock, W.D.; Stylianou, V.M.; Vaillancourt, R.E.; Byrne, M. Evidence for adaptation and acclimation in a widespread eucalypt of semi-arid Australia. *Biol. J. Linn. Soc.* **2017**, *121*, 484–500. [CrossRef]

78. Dijkstra, P. Cause and effect of differences in specific leaf area. In *Causes and Consequences of Variation in Growth Rate and Productivity of Higher Plants*; Lambers, H., Cambridge, M.L., Konings, H., Eds.; SPB Academic Publishing: The Hague, The Netherlands, 1989; pp. 125–140.

79. Reich, P.B.; Ellsworth, D.S.; Walters, M.B. Leaf structure (specific leaf area) modulates photosynthesis-nitrogen relations: Evidence from within and across species and functional groups. * Funct. Ecol.* **1998**, *12*, 948–958. [CrossRef]

80. Pruyn, M.L.; Gartner, B.L.; Harmon, M.E. Storage versus substrate limitation to bole respiratory potential in two coniferous tree species of contrasting sapwood width. *J. Exp. Bot.* **2005**, *56*, 2637–2649. [CrossRef] [PubMed]

81. Maier, C.A.; Zarnoch, S.J.; Dougherty, P.M. Effects of temperature and tissue nitrogen on dormant season stem and branch maintenance respiration in a young loblolly pine (*Pinus taeda*) plantation. *Tree Physiol.* **1998**, *18*, 11–20. [CrossRef] [PubMed]

82. Guimarães, Z.T.M.; dos Santos, V.A.H.F.; Nogueira, W.L.P.; Martins, N.O.D.A.; Ferreira, M.J. Leaf traits explaining the growth of tree species planted in a Central Amazonian disturbed area. *For. Ecol. Manag.* **2018**, *430*, 618–628. [CrossRef]

83. Farquhar, G.D.; von Caemmerer, S.; Berry, J.A. A biogeochemical model of photosynthetic CO₂ assimilation in leaves of C3 species. *Planta* **1980**, *149*, 78–90. [CrossRef]

84. Pirttioja, N.; Carter, T.R.; Fronzek, S.; Bindi, M.; Hoffmann, H.; Palosuo, T.; Ruiz-Ramos, M.; Tao, F.; Trnka, M.; Acutis, M.; et al. Temperature and precipitation effects on wheat yield across a European transect: A crop model ensemble analysis using impact response surfaces. *Clim. Res.* **2015**, *84*, 87–105. [CrossRef]

85. Chen, M.; Rafique, R.; Asrar, G.R.; Bond-Lamberty, B.; Ciais, P.; Zhao, F.; Reyer, C.P.O.; Ostberg, S.; Chang, J.; Ito, A.; et al. Regional contribution to variability and trends of global gross primary productivity. *Environ. Res. Lett.* **2017**, *12*, 105005. [CrossRef]

86. Fyllas, N.M.; Bentley, L.P.; Shenkin, A.; Asner, G.P.; Atkin, O.K.; Díaz, S.; Enquist, B.J.; Farfan-Rios, W.; Gloor, E.; Guerrieri, R.; et al. Solar radiation and functional traits explain the decline of forest primary productivity along a tropical elevation gradient. *Ecol. Lett.* **2017**, *20*, 730–740. [CrossRef] [PubMed]

87. Fyllas, N.M.; Christopoulou, A.; Galanidis, A.; Michelaki, C.Z.; Giannakopoulos, C.; Dimitrakopoulos, P.G.; Arianoutsou, M.; Gloor, M. Predicting species dominance shifts across elevation gradients in mountain forests in Greece under a warmer and drier climate. *Reg. Environ. Change* **2017**, *17*, 1165–1177. [CrossRef]

88. Schwinning, S.; Ehleringer, J.R. Water use trade-offs and optimal adaptations to pulse-driven arid ecosystems. *J. Ecol.* **2001**, *89*, 464–480. [CrossRef]
89. Li, C.; Wang, K. Differences in drought responses of three contrasting *Eucalyptus microtheca* F. Muell. populations. *For. Ecol. Manag.* **2003**, *179*, 377–385. [CrossRef]

90. Zhou, S.X.; Medlyn, B.E.; Prentice, I.C. Long-term water stress leads to acclimation of drought sensitivity of photosynthetic capacity in xeric but not riparian *Eucalyptus* species. *Ann. Bot.* **2016**, *117*, 133–144. [CrossRef] [PubMed]

91. Latimer, J.G.; Severson, R.F. Effect of mechanical and moisture-stress conditioning on growth and cuticle composition of broccoli transplants. *J. Am. Soc. Hort. Sci.* **1997**, *122*, 788–791. [CrossRef]

92. Grill, D.; Tausz, M.; Pöllinger, U.; Jiménez, M.S.; Morales, D. Effects of drought on needle anatomy of *Pinus canariensis*. *Flora* **2004**, *199*, 85–89. [CrossRef]

93. Bacelar, E.A.; Santos, D.L.; Moutinho-Pereira, J.M.; Gonçalves, B.C.; Ferreira, H.F.; Correia, C.M. Immediate responses and adaptive strategies of three olive cultivars under contrasting water availability regimes: Changes on structure and chemical composition of foliage and oxidative damage. *Plant Sci.* **2006**, *170*, 596–605. [CrossRef]

94. Pospíšilová, J.; Synková, H.; Haisel, D.; Semorádová, Š. Acclimation of plantlets to ex vitro conditions: Effects of air humidity, irradiance, CO₂ concentration and abscisic acid (a Review). *Acta Hortic.* **2007**, *748*, 29–38. [CrossRef]

95. Navarro, M.; Ayax, C.; Martinez, Y.; Laur, J.; El Kayal, W.; Marque, C.; Teulières, C. Two EguCBF1 genes overexpressed in *Eucalyptus* display a different impact on stress tolerance and plant development. *Plant. Biotechnol. J.* **2011**, *9*, 50–63. [CrossRef]

96. Tapia, G.; Acuña, H.; Inostroza, L. Plant drought tolerance: Some genetics and agronomics relevant aspects for breeding in forage species. In *Droughts: New Research*; Neves, D.N., Sanz, J.D., Eds.; Nova Science Publishers: Hauppauge, NY, USA, 2013; pp. 157–187.

97. Tapia, G.; Morales-Quintana, L.; Parra, C.; Berbel, A.; Alcorta, M. Study of nsLTPs in *Lotus japonicus* genome reveal a specific epidermal cell member (LjLTP10) regulated by drought stress in aerial organs with a putative role in cutin formation. *Plant Mol. Biol.* **2013**, *82*, 485–501. [CrossRef] [PubMed]