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

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.

1. Introduction

Anthropogenic climate change causes range shifts and might lead to the extinction of numerous plant and animal species [1,2,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 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,7,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,25,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,29,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 variations. 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 investigate 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. Location, altitude, sampling period, numbers of covered plots and of sampled trees.

Site	Region	Longitude	Latitude	Altitude (m asl)	Period	Plot Number	Tree Number
1	Rif	−4.9481	35.0295	1524	Spring 2017	2	44
2	Rif	−4.7064	34.9597	1551	Spring 2017	2	31
3	Rif	−4.3356	34.8703	1768	Spring 2017	1	15
4	Middle Atlas	−5.2564	32.9622	2205	Autumn 2017	2	37
5	Middle Atlas	−5.4316	32.9430	1673	Autumn 2017	2	35
6	Middle Atlas	−5.2750	32.9442	1960	Autumn 2017	2	56
7	Middle Atlas	−5.4290	32.9722	1581	Autumn 2017	3	24
8	Middle Atlas	−5.4767	32.9735	1590	Autumn 2017	2	34
9	Middle Atlas	−5.4539	32.8101	1488	Autumn 2017	8	15
10	Middle Atlas	−5.3812	32.9539	1697	Autumn 2017	2	16
11	Middle Atlas	−5.2317	32.8888	1980	Autumn 2017	1	31
12	Rif	−4.4243	34.8233	1794	Autumn 2018	2	48
13	Rif	−4.5416	34.8519	1820	Autumn 2018	1	37
14	Middle Atlas	−5.0287	33.5529	1749	Autumn 2018	2	38
15	Middle Atlas	−5.1120	33.4972	1717	Autumn 2018	2	44
16	Middle Atlas	−5.1995	33.4122	1755	Autumn 2018	2	41

). 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.

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)$$

(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.

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).

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).

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).

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.

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 (Figure 4 and Figure 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 (Figure 6, Figure 7 and Figure 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 (Figure 9 and Figure 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,69,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,92,93,94,95,96,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 (Figure 6 and Figure 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 relationship 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.