Biogeographic historical legacies in the net primary productivity of Northern Hemisphere forests

Timo Conradi,* Koenraad Van Meerbeek, Alejandro Ordonez and Jens-Christian Svenning Abstract It has been suggested that biogeographic historical legacies in plant diversity may influence ecosystem functioning. This is expected because of known diversity effects on ecosystem functions, and impacts of historical events such as past climatic changes on plant diversity. However, empirical evidence for a link between biogeographic history and present-day ecosystem functioning is still limited. Here, we explored the relationships between Late-Quaternary climate instability, speciespool size, local species and functional diversity, and the net primary productivity (NPP) of Northern Hemisphere forests using structural equation modelling. Our study confirms that past climate instability has negative effects on plant functional diversity and through that on NPP, after controlling for present-day climate, soil conditions, stand biomass and age. We conclude that global models of terrestrial plant productivity need to consider the biogeographical context to improve predictions of plant productivity and feedbacks with the climate system.


INTRODUCTION
A research priority of global relevance is to develop robust predictions of climate-change impacts on the functioning of ecosystems. Crucial for this goal is to understand the drivers of terrestrial plant productivity, the process of plant carbon (C) uptake and conversion to biomass, which is the basis of the dynamics of the biosphere and its feedbacks with the climate system. Consequently, net primary productivity (NPP) is an important output of Earth System Models (ESMs) and Dynamic Global Vegetation Models (DGVMs), mechanistic models that have become the primary tool for forecasting the future states and functioning of the biosphere. However, current state-of-the-art ESMs do not reproduce well the spatial pattern of NPP dynamics as documented by Earth observation satellites (Smith et al. 2015), implying that ESMs cannot reliably predict future changes. Reasons for this could be that important mechanisms are not well represented in ESMs or that the models are not well parametrised because they ignore biogeographic contingencies in ecosystem functioning Moncrieff et al. 2016;Higgins 2017). If ecosystem functioning is contingent on biogeographic history, then models parametrised for single biogeographic regions cannot extrapolate globally and, vice versa, predictions from globally parametrised models may deviate strongly from local observations even though the mechanisms that drive NPP were correctly represented (Moncrieff et al. 2016). The implication would be that, as Higgins (2017) put it, "ecologists need to broaden from biophysical determinism towards an acknowledgement of biogeographical context" when modelling ecosystem functioning at large spatial scales.
Biogeographic history could influence NPP because it has shaped present-day patterns in species and functional diversity , and there is accumulating evidence that diversity drives NPP. For example, experiments and observational studies in grasslands and forests demonstrate a positive effect of plant species richness on NPP (Vil a et al. 2007;Grace et al. 2016;Liang et al. 2016;Oehri et al. 2017;van der Plas 2019). This effect has been attributed to higher complementarity in resource use (Hooper et al. 2005), a reduced density-dependent abundance of herbivores and pathogens (Civitello et al. 2015), facilitation, and a greater chance that productive species are included in species-rich communities (Hooper et al. 2005). A signal of biogeographic history in the richness-NPP relationship is conceivable because local species richness is influenced by the size of the regional species pool (Lalibert e et al. 2014;Conradi & Kollmann 2016;Ricklefs & He 2016), and pool size is shaped by biogeographic history (Ricklefs 2004;Belmaker & Jetz 2015;Eiserhardt et al. 2015;Jim enez-Alfaro et al. 2018). Quaternary climate oscillations, in particular, have left a significant legacy in present-day species richness patterns (Latham & Ricklefs 1993;Eiserhardt et al. 2015) that could contribute to geographic variation in NPP.
Rather than simply through species richness, biogeographic history could influence NPP via its legacy effect on plant functional trait diversity. Functional trait diversity within a plant community may be a better predictor of NPP than species richness because it more directly reflects niche complementarity (Morin et al. 2011;Paquette & Messier 2011), facilitation potential (Dawson 1993) and differences in herbivore and pathogen defence (Carmona et al. 2011). Plant functional trait diversity can be strongly influenced by biogeographic history. European and North American angiosperm assemblages, for example, exhibit lower functional diversity in areas subject to high Quaternary glacial-interglacial climate change velocities , 2016a, probably reflecting non-random extinctions related to frost tolerance and correlated traits as climate repeatedly cooled and warmed during the Pleistocene (Svenning 2003;Bhagwat & Willis 2008;Eiserhardt et al. 2015). Despite evidence for historical effects on diversity and diversity effects on NPP, it has yet to be assessed if NPP depends on biogeographic history. If true, this would have implications for how researchers represent biogeographic factors in models of NPP.
In this study, we explore if biogeographic history influences the NPP of Northern Hemisphere forests via its effects on species and functional diversity. We focus on temperate and boreal forest stands in Europe, western North America and eastern North America. These three regions have distinct tree floras and were affected to different extents by extinctions during the Neogene and Pleistocene climate cooling (Latham & Ricklefs 1993;Svenning 2003;Eiserhardt et al. 2015;Ordonez & Svenning 2018). In addition, range maps exist for the complete tree floras of all three regions (Montoya et al. 2007), and there is good coverage of functional trait data. This makes these forests suitable study objects to explore effects of biogeographic history on NPP and how species pools and functional trait diversity mediate these effects. Figure 1 shows a new conceptual model of forest NPP that we evaluated in this study, representing the general hypothesis that effects of biogeographic history on regional and local species richness and functional trait diversity will be important determinants of NPP. The model is derived from specific hypotheses about the drivers of diversity patterns based on the biogeographic and diversity-productivity literature: (1) Climate-change velocity from the Last Glacial Maximum (LGM) to present has influenced the size of local species pools. Climate-change velocity [km year À1 ] quantifies how fast a hypothetical plant must spread over the surface of the Earth during two time points to maintain constant climate (Loarie et al. 2009;Corlett & Westcott 2013).
LGM-to-present climate-change velocity represents the magnitude of local climatic changes between the LGM and present (and between glacial and interglacial maxima more generally) and is expected to exert negative influences on species-pool size because extinctions and colonisation lags should be higher in places with high climate-change velocity (Dynesius & Jansson 2000;Jansson 2003).
(2) A region effect exists on species-pool size. The region effect (Ricklefs & He 2016) encompasses differences in speciation and extinction among Europe, western North America and eastern North America that we could not directly quantify in this study, some of which dating back deeper in time than the Pleistocene (Svenning 2003;Eiserhardt et al. 2015). For example, the region effect could reflect differences in the accessibility of climate refugia among the three regions (Latham & Ricklefs 1993) that are not captured by local LGM-to-present climate-change velocity.
(3) Species-pool size positively influences the tree species richness of individual forest stands as shown in recent studies (Ricklefs & He 2016;Jim enez-Alfaro et al. 2018). (4) The functional diversity of forest stands increases with species richness. This is expected unless species in a community are completely functionally redundant. Consequently, most functional diversity metrics increase with species richness, at least up until a richness of c. 20 species (Lalibert e & Legendre 2010; Mouchet et al. 2010). The stands in our dataset have a species richness between 1 and 15, which are typical values for temperate and boreal forest stands (e.g. Paquette & Messier 2011). (5) The functional diversity of forest stands influences their NPP. This is expected due to a number of mechanisms mentioned above.
Using structural equation modelling, we show that the NPP of Northern Hemisphere forest stands is influenced by functional diversity, and via this effect by biogeographic history, as functionally diverse stands are associated with areas that have been climatically stable during the glacial-interglacial oscillations, likely reflecting that their species pools have been less pruned by extinctions and colonisation lags. These results demonstrate that biogeographic historical legacies in plant diversity influence ecosystem functioning.

Data compilation
To test the above hypotheses, we used the global terrestrial plant production dataset of Michaletz et al. (2014). This highquality dataset has been used in previous studies of forest productivity (Michaletz et al. 2014(Michaletz et al. , 2018 and contains information on the NPP, biomass, age and geographic coordinates of forest stands. We used only the North American and European stands of the dataset ( Fig. 2; n = 61 after removing two plots without species lists [see below]) because range maps for the complete tree flora (height > 4 m) were only available for these continents (Montoya et al. 2007). We obtained an updated version of the Montoya et al. (2007) range maps in June 2010, available upon request from these authors. These maps were used to compute the species-pool size of individual forest stands. This was accomplished by projecting the range maps on a 10-km Behrmann equal-area grid and then counting all species whose ranges overlapped in a grid cell with one of our NPP stands. Because North America has two distinct biogeographic regions of temperate tree flora (Latham & Ricklefs 1993), we grouped the data into three regions: western North America, eastern North America and Europe. One forest stand had wrong coordinates in the terrestrial plant production dataset and fell into the Atlantic Ocean (row ID 60 in Michaletz et al. 2014). After checking the original reference (Whittaker & Woodwell 1969), we corrected its latitude to 40.4 decimal degrees.
The LGM-to-present temperature velocities at the location of each forest stand were extracted from velocity maps of Sandel et al. (2011). We see these velocities as representative of the repeated late-Quaternary glacial-interglacial velocities more generally. For the computation of temperature velocity, Sandel et al. (2011) used the average LGM temperature prediction of two climate models, CCSM3 (Collins et al. 2006;Otto-Bliesner et al. 2006) and MIROC3.2 (k-2004model developers 2004, because the LGM temperature predictions were highly correlated. The terrestrial plant production dataset does not contain information on the species richness of the forest stands. As we were interested in how species-pool size affected forest NPP via stand richness and functional trait diversity, we compiled tree species lists for each stand from the original studies contained in the terrestrial plant production dataset (see Table S1). Four original studies only reported the presence of 'Carya spp.', implying that more than one species of this genus were present in these NPP plots (row IDs 1225, 1227 and 1228 in Michaletz et al. 2014). These must be two or more of C. cordiformis, C. glabra, C. ovata and C. tomentosa, as these are the only Carya species whose ranges overlap the respective NPP plot locations. Dale et al. (1990) also list exactly these species in another description of the same sites (Oak-Hickory forests at Walker Branch, TN; USA), and we thus assumed that all four Carya species were present in all three plots.
We collected data on five functional traits (leaf area, specific leaf area [SLA], maximum height, wood density, seed   (Westoby 1998;Chave et al. 2009;Moles et al. 2009;Carmona et al. 2011). We expected tree communities with a diverse combination of these traits to exhibit higher levels of complementarity in light use (as far as height is considered; Sapijanskas et al. 2014) and to have higher and more stable productivity over time. The latter is expected due to the ability of functionally diverse stands to compensate for reduced productivity or high mortality of individual species by increased production and recruitment of other species when environmental conditions fluctuate temporarily (Allan et al. 2011), and because of a reduced likelihood that specific consumers or pathogens reach high abundances when tree species with different seed, wood and leaf traits are present. Growth form (shrub or tree) was included as additional trait to represent a number of architectural differences between shrubs and trees that may be important for productivity, but are not captured by the above-mentioned traits. For example, in comparison with trees, shrubs often have higher foliage density, relative crown width and resprouting vigour, lower crown volume and area, and different diameter-height allometries (Zizka et al. 2014). If species-level trait information was not available, we used the average value for the genus, and if that was also not available, we took the average value for the family. We used the BIEN R package (Maitner et al. 2018) to download data from BIEN and the references for these data are listed in Appendix S1. The trait data were used to calculate Functional Dispersion (FDis; Lalibert e & Legendre 2010). FDis measures the spread of trait values around the centre of the multivariate trait space of a community and has thus been used to quantify niche complementarity (Paquette & Messier 2011;Finegan et al. 2015;Ratcliffe et al. 2016). We computed FDis based on the Gower distances between the species, calculated from the untransformed trait values.
Finally, mean annual temperature and annual precipitation were extracted from the WorldClim global climate layers with 1-km spatial resolution (Hijmans et al. 2005). Soil pH and soil cation exchange capacity (CEC) for six soil horizons representing the upper 100 cm of the soil were extracted from 250-m resolution grids (Hengl et al. 2017) and the median values were computed across these six horizons. These data were included in the statistical modelling to control for the effects of present-day environmental conditions. At the coordinates of one plot, soil pH and CEC were not available, and we used predictive mean matching, implemented in the mice R package (van Buuren & Groothuis-Oudshoorn 2011) to impute these two missing values based on soil pH, CEC, mean annual temperature, annual precipitation, region, NPP, stand biomass and age. Table S2 shows the value of all the variables for all study plots used in this study.

Data analysis
We used Bayesian piecewise structural equation modelling to examine the hypothesised causal relationships between biogeographic history and the NPP of forest stands (Lefcheck 2016;B€ urkner 2018). The model controlled for the effects of present-day environmental conditions, stand age and biomass (see Fig. 3a). Specifically, the model assumed that.
• Stand biomass influences NPP because biomass controls total leaf area, which drives photosynthesis (Michaletz et al. 2014).
• Stand age influences NPP directly through age-related declines in growth (Gower et al. 1996), but also indirectly and positively via its effect on biomass, as old stands have accumulated more biomass (Michaletz et al. 2018).
• Mean annual temperature and precipitation as well as soil conditions influence NPP directly through constrains on physiological rates (Berry & Bjorkman 1980;Evans 1989;Medrano & Flexas 2002), but also indirectly via constraining the biomass of forest stands (Michaletz et al. 2014(Michaletz et al. , 2018. • Mean annual temperature and precipitation have an additional indirect effect on NPP that is mediated by their effects on species-pool size. A correlation between presentday climate and broad-scale variation in species-pool size has been reported in previous studies (Currie & Paquin 1987;O'Brien 1993), albeit the causal mechanisms behind this relationship remain elusive (Currie et al. 2004).
Soil CEC was excluded from the analysis because it was highly collinear with mean annual precipitation (Pearson r = 0.72). Beyond a threshold of |r| >0.7, collinearity begins to severely distort model estimation (Dormann et al. 2013). NPP was log-transformed and continuous variables were standardised to mean = 0 and SD = 1 to aid model convergence (except for biomass, which was scaled to SD = 1 only). Pool size and species richness were not transformed and were modelled with a negative binomial distribution to account for overdispersion, and biomass was modelled with a Gamma distribution, both with a log link function. The model was fit using the R package brms (B€ urkner 2017) with 2 chains, 10,000 iterations and a warm-up of 1000 runs. Using the standard priors, the model converged with b R values close to 1 (Gelman and Rubin's diagnostic) and effective sample sizes > 5000 for all estimated effects. Population-level effects were assumed to be significant if the 95% credible interval of an effect did not include zero. Model validation was performed using approximate leave-oneout cross-validation (LOOIC; Vehtari et al. 2017) using the loo package (Vehtari et al. 2019).
The Pareto shape k parameter was > 0.7 for two observations, indicating outliers with high influence and inaccurate LOO approximations (Vehtari et al. 2017) (see detailed model statistics in Table S3 in Supporting Information). One of these observations had much higher NPP (4370 g m À2 year À1 ) than any other of the forest plots in the dataset (second largest NPP: 2520 g m À2 year À1 ; median NPP: 1304 g m À2 year À1 ), and the NPP value is unlikely to be correct, given the young age of this forest stand (26 years). In fact, even when the global terrestrial plant production dataset of Michaletz et al. (2014) with > 1200 plots is considered, this stand still exceeds the second largest NPP by > 1000 g m À2 year À1 . We removed this observation from all analyses because it belonged to a range of NPP values that is not otherwise represented in our dataset and the global dataset, but had a big influence on the relationship between FDis and NPP and was thus flagged by the LOOIC validation. We additionally removed the second observation with Pareto k > 0.7, but this had little quantitative effects on parameter estimates and did not change their direction or significance (see Table S3). In the main text, we present results from models without these two observations. Trace plots were used to assess mixing and convergence of the two chains (Appendix S2). The model was further validated with posterior predictive checks (Appendix S3), using the bayesplot R package (Gabry & Mahr 2019).
We also tested alternative models (see Table S3). In the first one, we included direct effects of temperature, precipitation, soil pH, stand biomass and age on species richness (Fig. 3b). Trace plots and posterior predictive checks for this model are shown in Appendices S4 and S5. In the second alternative model, we included effects of (log-transformed) community mean trait values on NPP, and effects of temperature velocity and region on community mean trait values. Because of high Figure 3 Structural Equation Models depicting direct and indirect drivers of forest net primary productivity (NPP). In (a), species-pool size is the only driver of richness, whereas in (b) effects of climate, soils and stand structure on richness are also evaluated. The width of the arrows is proportional to their relative effect size (see Table 1). Dashed arrows show non-significant effects. Blue arrows represent positive effects and red arrows represent negative effects. The black arrow is for the categorical variable 'Region'. Age = age of forest stand; pH = soil pH; MAP = mean annual precipitation; MAT = mean annual temperature; Velo Temp = Last Glacial Maximum-to-present temperature velocity; Region = biogeographic region (i.e. Europe, eastern North America or western North America); Richness = tree species richness of forest stand; FDis = Functional Dispersion of forest stand. Note that the estimates for the species-pool effect on species richness and the estimate for the effect of species richness on FDis are not comparable with the other estimates. This is because the units of species-pool size and species richness are species, whereas the other variables were scaled prior to modelling. collinearity (r > 0.7) with other community mean trait values, leaf area and wood density were excluded from this model. Specifically, mean wood density was highly correlated with three mean traits: height (r = À0.82), SLA (r = 0.80) and seed mass (r = 0.78). Mean leaf area was highly correlated with mean SLA (r = 0.78) and FDis (r = 0.77). We retained mean SLA instead of mean leaf area because (1) we would have to remove FDis otherwise, which is a variable we were interested in, and (2) because SLA may be more informative of resource investment in growth than leaf area. Finally, the third alternative model added to the second model the effects of temperature, precipitation and soil pH on community mean trait values. Both models with community mean trait values could not be fit with our dataset. Each time we removed an observation with Pareto k value > 0.7, new observations were flagged as problematic (k > 0.7), indicating that these two models may be too complex for our dataset. The two models with community mean traits were thus discarded.

RESULTS
The results of our models revealed a clear signal of biogeographic history in the NPP of Northern Hemisphere forests stands after controlling for present-day climate, soils, stand biomass and stand age ( Fig. 3; Table 1; Table S3). LGMto-present temperature velocity had a significant negative effect on species-pool size that was similar in magnitude to present-day mean annual temperature. Species pools were also significantly larger in eastern North America than in western North America and Europe. Species-pool size trickled down to influence stand-level tree species richness, and more species-rich stands also had higher FDis. Finally, elevated FDis led to higher NPP of forest stands (Fig. 4). Note that in Table 1 (and Table S3) and Fig. 3, the estimates for the species-pool effect on species richness and the estimate for the effect of species richness on FDis are not comparable with the other estimates. This is because the Model with pool size as single driver of species richness (Fig. 3a) Model with additional drivers of species richness (Fig. 3b units of species-pool size and species richness are species, whereas the other variables were scaled prior to modelling (see methods). Stand biomass was the strongest determinant of NPP, with greater NPP in biomass-rich stands. Stand age exerted a negative direct effect on NPP, but because older stands had a greater biomass, there was also a positive indirect effect of age on NPP (Fig. 3a). We found no significant direct effects of present-day climate on NPP, but climate did have indirect effects. This was because species pools were larger in warmer areas, and species pools had positive indirect effects on NPP that were mediated by stand richness and FDis, as described above. However, species-pool size and thus richness and FDis were more strongly driven by biogeographic factors than by present-day climate (Table 1). Soil pH did not influence NPP and biomass.
The link between biogeographic history and NPP was also confirmed in the alternative SEM in which species richness was not only influenced by species-pool size but also by climate, soil and stand structure variables (Fig. 3b). Despite adding additional drivers of species richness, the estimated effect of species-pool size on species richness was unchanged and remained significantly positive, albeit confidence intervals were slightly larger ( Table 1). The difference in LOOIC values between the alternative model and our main model was < 2, and there is thus no statistical reason to favour one model over the other (cf. Burnham & Anderson 2002). We thus present both models in Fig. 3. Raw bivariate plots of the relationships between the response variables and explanatory variables with significant effects in the SEMs are shown in Appendix S6.

DISCUSSION
Our study demonstrates that the NPP of temperate forest stands is contingent on biogeographic history. Specifically, our results reveal a link between late-Quaternary climate instability and region history with NPP that is mediated by their effects on species pools and the species and functional diversity of forest stands (Fig. 3). The results confirm the hypothesis that present-day ecosystem functioning reflects a long history of environmental dynamics via their cumulative effects on biological diversity , and suggest that global models of terrestrial plant productivity need to consider the biogeographic context.
Region was particularly important in explaining the size of the species pools from which local forest stands could assemble (Table 1). This region effect may reflect differences in speciation and extinction through deep time (Ricklefs & He 2016) that could be explained by possible differences in the areal extent of the regions through time (Belmaker & Jetz 2015), topography effects on diversification or differential connectivity to the subtropics (Latham & Ricklefs 1993;Ricklefs 2004), but certainly reflects the contrasting extents to which the respective floras have been pruned during Neogene and Pleistocene cooling (Latham & Ricklefs 1993;Svenning 2003). However, there may be further differences between regions that were unaccounted for in our analysis, such as climate and soil parameters or biotic interactions, that are not well captured in the available global data. In addition to the broad-scale differences in region history, local LGM-topresent temperature velocity also exerted significant and negative effects on species-pool size. This is in agreement with previous studies showing LGM climate effects on tree species richness patterns within continents (Montoya et al. 2007;Svenning & Skov 2007), what may be attributed to more frequent extinctions and colonisation lags in areas with high climate instability (Dynesius & Jansson 2000;Jansson 2003).
These historical effects on species pools trickled down to influence the species richness of single forest stands. This is in line with recent studies that have identified species pool size as a key determinant of local richness (Lalibert e et al. 2014;Ricklefs & He 2016;Jim enez-Alfaro et al. 2018), suggesting that local richness may often be unsaturated and constrained by broad-scale historical factors rather than local environmental conditions (Mateo et al. 2017). Species-rich stands also had higher functional diversity (Table 1; Fig. 3). Eastern North American stands in particular had consistently high functional diversity (Fig. 4). Some of the genera that occurred in these stands also occurred in the European stands (e.g. Acer, Betula, Quercus), but were often represented with multiple species each (see Table S1), in contrast to usually only one in the European stands. In addition, eastern North American stands contained genera that went extinct in Europe during the Plio-Pleistocene cooling events such as Carya, Liriodendron, Nyssa and Tsuga. The extinctions in Europe were deterministic with respect to phylogeny and plant traits (Svenning 2003;Eiserhardt et al. 2015), which provides an explanation for the reduced functional diversity of the European stands. Some of the western North American stands also had a high functional diversity (Fig. 4), namely, those with a diverse mix Effect of plant functional diversity on forest net primary productivity (NPP). Shown is the partial relationship of Functional Dispersion (FDis) with NPP after controlling for the effects of climate, soils, stand biomass and age. Units are residual deviations from predicted partial scores. The adjusted R 2 is from a linear regression model. Both variables were scaled to zero mean and unit variance prior to modelling. of broadleaf and coniferous tree species (see Table S1), again including genera lost from Europe during the Plio-Pleistocene.
Higher functional diversity was found to significantly enhance NPP when controlling for environmental conditions, stand biomass and age (Fig. 4). The positive effect of functional diversity on NPP is consistent with the idea that higher complementarity in resource use enables increased resource exploitation at the ecosystem level and thus higher NPP (Hooper et al. 2005;Paquette & Messier 2011). For example, a stand with tree species of different height or the presence of shrubs may have enhanced stand-level light interception through denser packing of vertical space (Sapijanskas et al. 2014). In addition, the presence of a range of seed, leaf and wood characteristics within a forest may prevent high abundances of consumers and pathogens (Civitello et al. 2015), and may contribute to stabilise NPP across years with variable environmental conditions (Allan et al. 2011). It is conceivable that we even underestimated the importance of functional diversity. We had no trait data for traits such as crown and root structure, which may have provided even better estimates of complementarity in light interception or exploitation of soil nutrients and water (Hardiman et al. 2011;Jucker et al. 2015).
Although we focus here on biogeographic factors, the most important drivers of NPP were stand biomass and age. Biomass controls the leaf mass available for photosynthesis and is thus an important driver of NPP. By controlling for biomass, we found stand age to negatively affect NPP, which is likely due to age-related declines in tree growth (Gower et al. 1996). However, because stand biomass was larger in older stands, age had a positive indirect effect on NPP. Notably, present-day climate did not influence NPP directly in our forests stands, even though our dataset covered a broad mean annual precipitation (556-2183 mm) and temperature (À0.1 to 13.7°C) gradient. Climate influenced NPP indirectly, however, via its positive effect on species-pool size and functional diversity. The NPP data analysed in this study is a subset of the terrestrial plant production dataset analysed by Michaletz et al. (2014Michaletz et al. ( , 2018. Using different ways to analyse the data, these authors also found no to little direct effects of climate on NPP. They concluded that climate controls NPP only indirectly via its effects on stand structure rather than by influencing the kinetics of photosynthesis and respiration, in line with what others had shown before (e.g. Chapin 2003; K€ orner 2006; but cf. Chu et al. 2016).
Our finding of a biogeographic contingency in terrestrial plant productivity has important implications for modelling NPP at large spatial scales. Global models of NPP need to both represent functional diversity as an important direct driver of NPP and consider that the functional diversity at a site is contingent on that site's biogeographic history and thus cannot be predicted from contemporary environmental conditions alone. Neither is currently the case in most state-of-the-art ESMs and DGVMs. Most models represent the global functional diversity of plants by a predefined number of plant functional types (PFTs). The implicit assumption is that there is no functional variation within PFTs, e.g. that all summer-green trees are functionally identical. However, the results of this and previous studies suggest that within-PFT functional diversity is an important source of variation in NPP (Paquette & Messier 2011), and that there is a strong influence of biogeographic history on the spatial distribution of this within-PFT functional diversity . Trait-based DGVMs (e.g. Pavlick et al. 2013;Scheiter et al. 2013) and vegetation demographic models (Fisher et al. 2018) may provide a way forward, but there is currently no workflow for calibrating these models to account for the uneven distribution of functional diversity. One reason is that we currently do not have sufficient data to adequately characterise species and functional diversity in most parts of the world; this may improve in the future as databases on species distributions and traits continue to grow, e.g. BIEN (http://bien.nceas.ucsb.edu/bien/) or TRY (Kattge et al. 2011). A second reason is that the focus with trait-based DGVMs (Pavlick et al. 2013;Scheiter et al. 2013) has been on exploring what the "optimal" trait combinations are, and not how biogeographic history may have moved the system away from this "optimal" trait spectrum (Higgins 2017).
Echoing Higgins (2017), we propose two solutions to account for biogeographic contingencies in large-scale models of NPP. One way is to parametrise regional subsets (e.g. grid tiles in global models) of such models separately to account for historical idiosyncrasies of the subsets (Higgins 2017), or to calibrate models separately for recognised biogeographic regions (see e.g. Moncrieff et al. 2016). The optimal size of these regions will depend on the study question and geographic extent, but could follow existing delimitations of biogeographic regions (Good 1947;Takhtajan 1986;Cox 2001). A second more ambitious solution may be to build mechanistic models of how traits or individual species assemble from regional species pools to local ecosystems. Mechanistic predictions of ecosystem assembly and functioning could be achieved by linking physiology-based distribution models that simulate plant growth (e.g. Higgins et al. 2012) with demographic range models (e.g. Schurr et al. 2012) and models of trait-based selection in local assemblages (e.g. Scheiter et al. 2013). Massive global databases of species occurrence records and functional trait data such as BIEN and TRY have been compiled in recent years and are still growing. Together with advancements in hierarchical and inverse modelling strategies, this means that process-based range and niche modelling for many thousands of species with different traits is on the horizon (Evans et al. 2016). The implication is that we may soon have the data and modelling tools to jointly represent biogeographic, plant physiological and community ecological theories in a new class of mechanistic models that may lead to enhanced predictions of terrestrial plant productivity.
Another implication of our results is that climate changedriven extinctions in trees may lead to long-lasting reductions in functional diversity and associated ecosystem productivity, highlighting that future extinctions driven by anthropogenic climate change may cast long shadows on the functioning of ecosystems.
JCS also considers this work a contribution to his VILLUM Investigator project "Biodiversity Dynamics in a Changing World" funded by VILLUM FONDEN (grant 16549) and the TREECHANGE project funded by Independent Research Fund Denmark | Natural Sciences (grant 6108-00078B). AO also considers this work a contribution to his AUFF-starting grant (AUFF-F-201 8-7-8).

DATA ACCESSIBILITY STATEMENT
No new data were used.