Complex Spatiotemporal Responses of Global Terrestrial Primary Production to Climate Change and Increasing Atmospheric CO2 in the 21st Century

Quantitative information on the response of global terrestrial net primary production (NPP) to climate change and increasing atmospheric CO2 is essential for climate change adaptation and mitigation in the 21st century. Using a process-based ecosystem model (the Dynamic Land Ecosystem Model, DLEM), we quantified the magnitude and spatiotemporal variations of contemporary (2000s) global NPP, and projected its potential responses to climate and CO2 changes in the 21st century under the Special Report on Emission Scenarios (SRES) A2 and B1 of Intergovernmental Panel on Climate Change (IPCC). We estimated a global terrestrial NPP of 54.6 (52.8–56.4) PgC yr−1 as a result of multiple factors during 2000–2009. Climate change would either reduce global NPP (4.6%) under the A2 scenario or slightly enhance NPP (2.2%) under the B1 scenario during 2010–2099. In response to climate change, global NPP would first increase until surface air temperature increases by 1.5°C (until the 2030s) and then level-off or decline after it increases by more than 1.5°C (after the 2030s). This result supports the Copenhagen Accord Acknowledgement, which states that staying below 2°C may not be sufficient and the need to potentially aim for staying below 1.5°C. The CO2 fertilization effect would result in a 12%–13.9% increase in global NPP during the 21st century. The relative CO2 fertilization effect, i.e. change in NPP on per CO2 (ppm) bases, is projected to first increase quickly then level off in the 2070s and even decline by the end of the 2080s, possibly due to CO2 saturation and nutrient limitation. Terrestrial NPP responses to climate change and elevated atmospheric CO2 largely varied among biomes, with the largest increases in the tundra and boreal needleleaf deciduous forest. Compared to the low emission scenario (B1), the high emission scenario (A2) would lead to larger spatiotemporal variations in NPP, and more dramatic and counteracting impacts from climate and increasing atmospheric CO2.


Introduction
Net Primary Productivity (NPP), a balance between photosynthetic carbon (C) uptake (Gross Primary Productivity; GPP) and losses due to plant respiration, represents the net C retained by terrestrial vegetation. It is of particular importance to humans since the largest portion of the food supply comes from terrestrial NPP [1]. NPP is also an important indicator of ecosystem health and services [2,3], and is an essential component of the global C cycle [4]. Terrestrial NPP is sensitive to multiple environmental changes including climate and atmospheric changes [5]. The IPCC Fourth Assessment (AR4) assessment indicated that global average temperature has increased by 0.74uC since the preindustrial times and that this trend is expected to continue through the 21 st century [6]. In addition, atmospheric CO 2 concentration have increased from the pre-industrial level of 280 ppm to the 2005 level of 379 ppm [6]. Comparing the 2090s with the 2000s, under the high emission scenario (A2), global mean temperature would increase by 4.6uC, while global annual precipitation would increase by 16.8%. However, the Representative Concentration Pathways (RCP's) scenarios used in the IPCC Fifth Assessment Report (AR5) IPCC report [7] were created in a different way and span a wider range of 21 st century projections. There are notable differences among the IPCC SRES scenarios and the IPCC RCP scenarios. The B1 scenario is very close to the RCP 4.5 by 2100, but there is lower emissions at the middle of 21 st century [8]. Similarly, the A2 scenario is between the RCP 6.0 and RCP 8.5 scenarios. Projected changes in atmospheric CO 2 concentration showed large increases under the A2 scenario, from 379.6 ppm in the 2000s to 809 ppm in the 2090s, which is equivalent to an overall increase of 113.8%. Such dramatic changes in climate and atmospheric composition would profoundly affect the NPP of terrestrial ecosystems. It is of critical importance to quantitatively analyze the contemporary pattern of global NPP and project to what extent climate change and increasing atmospheric CO 2 in the 21 st century would alter the magnitude and spatiotemporal patterns of NPP across the terrestrial ecosystems [9].
Previous studies reported that climate and increasing atmospheric CO 2 are the primary drivers for changes in global terrestrial NPP [4,[10][11][12]. Enhanced terrestrial NPP in response to increase in temperature [13] and atmospheric CO 2 concentration [14,15] have been suggested across a range of terrestrial ecosystems. Over the recent three decades, climate change has been the major driver of terrestrial NPP [10,16,17], with further benefits from CO 2 fertilization [4,18]. Nemani et al. [10] reported that climate change resulted in a 6% increase in global terrestrial NPP from 1982 to1999, with the largest increase in low-latitude ecosystems. Zhao and Running et al [16] reported that in the recent decade (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009), high temperature has increased water stresses and autotrophic respiration in the Southern Hemisphere resulting in a decline in global NPP by 0.55 PgC. Potter et al. [17], however, reported that rapid climate warming alleviated temperature limitations in high-latitude ecosystems which led to an increase in global terrestrial NPP by 0.14 PgC during 2000-2009. Higher temperature affect plant phenology, promoting an early growth and increasing the C assimilation in temperature-limited regions [13,19,20] due to the acceleration of enzymatic processes. Also, increasing CO 2 has been found to reduce stomatal conductance [14] and increase water use efficiency [21,22]. However, the stimulated effects of temperature on NPP may also be mitigated by increasing soil water stress and respiration rates induced by temperature rise [18,23,24]. In addition, increased water stress reduces nutrient uptake [25] which could potentially lead to a decline in productivity [26].While there is little doubt that climate change and increasing atmospheric CO 2 are the primary drivers of terrestrial NPP for the recent decade, the relative contribution of different drivers in the future is still unclear.
Process-based ecosystem models are effective tools for future projection of terrestrial NPP in response to global change [27]. Various process-based ecosystem models have been developed to estimate NPP response to changes in climate and increasing atmospheric CO 2 concentration at several scales from continental to global for both contemporary and future climatic conditions [4,28]. Previous modeling studies found that climate change resulted in an overall decline in global NPP, but doubling atmospheric CO 2 concentration resulted in an increase in global NPP by 16-25% [4,18]. In a process-based model comparison study, Cramer et al. [29] found differences in global NPP among 17 models (ranging from 44.4 to 66.3 PgC yr 21 ) due largely to how the water balance was represented in models. In a similar 17model comparison study, Friedlingstein et al. [30] found large uncertainties associated with belowground processes that resulted in different responses of NPP to global change factors across models. Thus, realistic historical assessments and future projections of global terrestrial NPP in a rapidly changing climate require more comprehensive models that include ecological, physiological and biogeochemical processes such as changes in phenology, length of growing seasons, nutrient dynamics, and ecohydrological processes.
The purpose of this study is to understand complex responses of terrestrial NPP at latitudinal, biome and global levels to projected climate change and increasing atmospheric CO 2 in the 21 st century. To accomplish this task, we first established the baseline estimate of global terrestrial NPP for the first decade of the 21 st century by using a well-evaluated process-based ecosystem model (the Dynamic Land Ecosystem Model, DLEM [22]) driven by multiple environmental factors. Then we used the DLEM model to examine responses of terrestrial NPP to projected climate change and increasing atmospheric CO 2 during the rest of the 21 st century under the IPCC Special Report on Emissions Scenarios (A2 and B1).The major objectives of this study are: (1) to estimate the contemporary global terrestrial NPP, (2) to project its changing trend in the 21 st century, (3) to attribute the relative contribution of climate, elevated CO 2 , and their interaction; and (4) to investigate the spatiotemporal pattern of global NPP as well as the response of different biomes to climate and CO 2 changes.

Model description
The DLEM is a highly integrated, process-based terrestrial ecosystem model that aims at simulating the structural and functional dynamics of land ecosystems affected by multiple factors including climate, atmospheric compositions (CO 2 , nitrogen deposition, and tropospheric ozone), land use and land cover change, and land management practices (harvest, rotation, fertilization etc). The DLEM has five core components ( Figure 1): 1) biophysics, 2) plant physiology, 3) soil biogeochemistry, 4) dynamic vegetation, and 5) land use and management [22]. This model has been extensively calibrated against various field data covering forest, grassland, and cropland from the Chinese Ecological Research Network, the US Long Term Ecological Research (LTER) sites, the AmeriFlux network and other field sites. Detailed information on how DLEM simulates these processes is available in our published papers [31][32][33][34]. Recently, we updated the model to DLEM 2.0 version, which is characterized by cohort structure, multiple soil layer processes, coupled C, water and nitrogen cycles, multiple greenhouse (GHG) emissions simulation, enhanced land surface processes, and dynamic linkages between terrestrial and riverine ecosystems [34,35]. Below, we briefly describe the simulation of GPP and NPP, calculation of relative CO 2 fertilization effect, input datasets used to drive the DLEM, and global-level evaluation of simulated NPP against satellite data.

Modeling gross primary productivity (GPP) in the DLEM
The DLEM uses a modified Farquhar's model to simulate GPP [36][37][38][39]. The canopy is divided into sunlit and shaded layers. GPP (gC m 22 day 21 ) is calculated by scaling leaf assimilation rates (mmol CO 2 m 22 s 21 ) up to the whole canopy: Where GPP sun and GPP shade are the GPP of sunlit and shaded canopy, respectively; A sun and A shade are assimilation rates of sunlit and shaded canopy; plai sun and plai shade are sunlit and shaded leaf area indices; dayl is daytime length (second) in a day. 12.01610 26 is a constant to change the unit from mmol CO 2 to gram C.
The plai sun and plai shade are estimated as: Where, proj LAI is the projected leaf area index. Using methods similar to Collatz et al. [37], DLEM determines the C assimilation rate as the minimum of three limiting rates, w c , w j , w e , which are functions that represents the assimilation rates as limited by the efficiency of photosynthetic enzymes system (Rubisco-limited), the amount of Photosynthetically Active Radiation (PAR) captured by leaf chlorophyll (light-limited), and the capacity of leaves to export or utilize photosynthesis products (export-limited) for C 3 species, respectively. For C 4 species, w e refer to the Phosphoenolpyruvate (PEP) carboxylase limited rate of carboxylation. The sunlit and the shaded canopy C assimilation rate can be estimated as: where c i is the internal leaf CO 2 concentration (Pa); O i is the O 2 concentration (Pa); C Ã is the CO 2 compensation point (Pa); K c and K o are the Michaelis-Menten constants for CO 2 and O 2 , respectively; a is the quantum efficiency; ø is the absorbed photosynthetically active radiation (W?M 22 ); V max is the maximum rate of carboxylation varies with temperature, foliage nitrogen concentration, and soil moisture [38]: where V max25 is the value at 25uC and a vmax is the temperature sensitivity parameter; f(Tday) is a function of temperature-related metabolic processes [36,37].  Where N leaf is the nitrogen concentration and N leaf_opt is the optimal leaf nitrogen concentration for photosynthesis.
b t is a function, ranging from one to zero that represents the soil moisture and the lower temperature effects on stomatal resistance and photosynthesis. where T min is the daily minimum temperature; w i is the soil water stress of soil layer i; ps i is the soil water potential of soil layer i, which is estimated from the soil volume water content based on equations from Saxton and Rawls [40]; r i is the root fractions distributed in soil layer i; psi_close and psi_open are the plant functional specific tolerance of the soil water potential for stomata overall close and open. The water stress in plants is a function of soil water potential which ranges from 0 to 1. Under no water limitations, the soil water stress of soil layer i (w i ) is equal to 1 where the soil water potential is at its maximum i.e., soil water potential when the stomata is opened (psi_open). Under frequent water stress, however, wi is calculated based on wilting point potential of specific plant functional types and depends on the balance between psi_open and psi_close.
Leaf stomatal resistance and leaf photosynthesis are coupled together through the following [37][38][39]. where r s is the leaf stomatal resistance, m is an empirical parameter, A is the leaf photosynthesis, c s is the leaf surface CO 2 concentration, e s is the leaf surface vapor pressure, e i is the saturated vapor pressure inside leaf, b is the minimum stomata conductance with A = 0, and P atm is the atmospheric pressure. Together in the following equations: where c a is the atmospheric CO 2 concentration, rb is the boundary resistance, e a is the vapor pressure of air, and stomatal resistance is the larger of the two roots of this quadratic [38].

Modeling net primary productivity (NPP) in the DLEM
NPP is the net C gain by vegetation and equals the difference between GPP and plant respiration, which is calculated as: The DLEM estimates maintenance respiration (Mr, unit: gC m 22 day 21 ) and growth respiration (Gr, unit: gC m 22 day 21 ) as a function of assimilated C, surface air temperature and biomass nitrogen content. Gr is calculated by assuming that the fixed part of assimilated C will be used to construct new tissue (for turnover or plant growth). During these processes, 25% of assimilated C is supposed to be used as growth respiration [41]. Maintenance respiration is related to surface temperature and biomass nitrogen content. The following is used to calculate the maintenance respiration of leaves, sapwood, fine roots, and coarse roots: Where i denotes the C pool of different plant parts (leaf, sapwood, fine root, or coarse root); Mr i (gC m 22 day 21 ) is the maintenance respiration of different pools; rf is a parameter indicating growing phase, which is set at 0.5 for the non-growing season and 1.0 for the growing season; R coeff is a plant functional type-specific respiration coefficient; N i (gN m 22 ) is the nitrogen content of pool i; f(T) is the temperature factor and is calculated as follows: Where T is the daily average air temperature for modeling aboveground C pools such as leaves, sapwood, and heartwood or soil temperature for modeling belowground pools such as coarse roots and fine roots.

The effect of CO 2 fertilization
In this study, we further quantified the effects of direct CO 2 fertilization on terrestrial NPP by calculating the 'beta' (b) effect. b effect measures the strength of changes in terrestrial NPP in response to increasing CO 2 concentration as follows: Where, NPP CO2 is the relative contribution of direct CO 2 fertilization on terrestrial NPP under the A2 and B1 scenarios, NPP clm+CO2 is the terrestrial NPP under the climate plus CO 2 experiment and NPP clm is the terrestrial NPP under the climate only experiment. CO 2 concentration (ppm) is the concentration of atmospheric CO 2 under the A2 and B1 scenarios.

Input datasets
The spatially-explicit data sets for driving the DLEM model include time series of daily climate, CO 2 concentration, annual land cover and land use (LCLU), nitrogen deposition, tropospheric ozone, and land management practices (irrigation and nitrogen fertilizer use). Other ancillary data include river network, cropping system, soil property, and topography maps. Contemporary vegetation map include 18 plant functional types ( Figure 2). Cropland and urban distribution datasets were developed by aggregating the 5-arc minute resolution HYDE v3.1 global cropland distribution data [42]. The vegetation map was developed based on global land-cover data derived from Landsat imageries [43], the National Land Cover Dataset 2000 (www.usgs. gov), and the global database of lakes, reservoirs, and wetlands [44]. The vegetation is transient and does not include any disturbance during the course of simulation. Half degree daily climate data (including average, maximum, minimum air temperature, precipitation, relative humidity, and shortwave radiation) were derived from newly available CRU-NCEP climate forcing data (1900-2009, 6-hour, half degree spatial resolution) [45]. The annual nitrogen deposition dataset for the historical period were based on Dentener [46]. Ozone AOT 40 data sets were based on the global AOT 40 index developed by Felzer et al. [47]. The gridded monthly CO 2 concentration data were derived from Multi-scale Synthesis and Terrestrial Model Intercomparison Project (MSTIMP) (http://nacp.ornl.gov/MsTMIP.shtml). Consumption of nitrogen fertilizers from 1961 to 2008 were derived from country level Food and Agriculture Organization of the United Nations (FAO) statistic database (http://faostat.fao.org). We then calculated the annual fertilization rate (gN m 22 ) as the ratio of national fertilizer application amount to total cropland area in each country [48]. The contemporary irrigation map was developed based on LCLU data and global irrigatied fraction map [49,50]. Long-term average climate datasets from 1900 to 1930 were used to represent the initial climate state in 1900.
For future projections, we used two IPCC emission secenarios (A2 and B1) datasets containing atmospheric CO 2 concentration and climate (precipitation and temperature) from the Community Climate System Model (CCSM3) (Figure 3-4). The A2 scenario (high emission scenario) is characterized by rapid population growth and low per capita income, with regionally oriented economic development. The B1 scenario (low emission scenario) describes the same global population as the A2 storyline, but it is less materially intensive in its service and information, economic structure, with emerging clean and resource-efficient technology [6]. The climate datasets were downloaded from the World Climate Research Programme's Coupled Model Intercomparison Project phase 3 (CMIP3) multi-model database (Meehl et al. [51]; www.engr.scu.edu/,emaurer/global_data). These datasets were downscaled as described by Maurer et al. [52] using the biascorrection/spatial downscaling method [53] to a 0.5 degree resolution, based on the 1950-1999 gridded observations of Adam and Lettenmaier [54]. For the future projections (2010-2099), we assumed that nitrogen deposition, ozone pollution, and LCLU remains unchanged from 2009. Temperature and precipitation have been projected to increase substantially during 2010 to 2099 ( Figure

Model parameterization, calibration and evaluation
The DLEM has been parameterized and applied across several regional and continental studies including Asia [31,33,[56][57][58][59][60], the United States [22,34,61] and North America [62,63] using longterm observational data for all defined plant functional types. The calibrated parameter values have been used to drive the model for specific plant functional types ( Figure 5). In this study, we compared DLEM-simulated global estimates of terrestrial NPP with MODIS NPP to evaluate model performance during 2000-2009. We first evaluated DLEM performance at the global level in simulating spatial pattern of terrestrial NPP by comparing DLEMsimulated NPP with MODIS NPP ( Figure 5). We then carried out a grid-to-grid comparison of DLEM-simulated NPP with MODIS product by randomly selecting 6000 grids from MODIS product (10% of the total sampling units). For bare land such as part of Africa, China and Mongolia (Taklamakan desert in West China, Gobi desert in Mongolia), there is no NPP data available from MODIS and we excluded those areas from analysis. At the global scale, the spatial pattern of DLEM-simulated NPP is consistent with that of MODIS NPP ( Figure 5A-B). In addition, we found a good agreement between MODIS and DLEM-simulated NPP for randomly selected 6000 grid points ( Figure 5C).The fitted line between DLEM-simulated and MODIS derived NPP had a slope of 0.73 and a high correlation coefficient (R 2 = 0.68). Our model evaluation for the effect of CO 2 fertilization is available in Lu et al. [60].

Experimental design and model implementation
We designed 262 factorial simulation experiments (2 simula-tions62 scenarios). The two simulation experiments include: (1) climate change only: only climate changes with time and other environmental factors are held constant during the study period (2010-2099); and (2) climate plus CO 2 : both climate and atmospheric CO 2 concentration change during the study period while other environmental factors including LCLU, nitrogen deposition, and tropospheric ozone are held constant at 2009 leve1s. The A2 and B1 climate scenarios were used to drive these two simulations.
The model simulation follows three important stages: an equilibrium simulation stage, a 3000-year spin-up stage, and a transient simulation stage. The model simulation begins with an equilibrium run with long-term average climate data for the period 1901-1930, with the 1900 levels of atmospheric CO 2 concentration, nitrogen deposition, and LCLU map to develop simulation baselines for C, nitrogen, and water pools. However, the tropospheric ozone data is kept at the 1935 level during the equilibrium. The simulation baseline (equilibrium) is approached when the net C exchange between the atmosphere and terrestrial ecosystems is less than 0.1 gC m 22 , the change in soil water pool is less than 0.1 mm, and the difference in soil mineral nitrogen content and nitrogen uptake is less than 0.1 gN m 22 among consecutive years. After the equilibrium run, a 3000-year spin up is carried out using transient climate data and LCLU distribution in the 1900 to eliminate system fluctuations caused by simulation      [64]. To examine whether the contemporary trend will continue, we further projected climate change effects on ecosystem productivity during the rest of the century (2010-2099).  Figure 7A; Table 1), but there would be a declining trend afterwards. Climate change under the A2 scenario would result in an overall decline in terrestrial NPP by 2.51 PgC (4.6%) in the 2090s compared to the 2000s (Figure 8; left panel). The B1 scenario shows an increasing trend with the highest increase in NPP by 1.57 PgC in the 2050s; however, NPP would level off after  Table 1). The magnitude of terrestrial NPP would be highest in low latitude regions (30.1-42.1 PgC yr 21 ), and lowest in high latitude regions (2.23-3.85 PgC yr 21 ). Compared to the 2000 s, the A2 climate would increase terrestrial NPP in midand high-latitudes by 4.6% and 34.8%, respectively, but would decrease terrestrial NPP in low latitude by 11.3% in the 2090s. The B1 climate scenario shows similar temporal pattern for midand high-latitidue regions with an increase in terrestrial NPP by 1.6% and 13.6%, respectively. The low latitude regions, however, show a declining trend with an increase in terrestrial NPP by 2.0% in the 2060s and 1.7% in the 2090s when compared to the 2000s. Under the A2 climate-only scenario, terrestrial NPP would increase in high latitude region by 0.82 PgC while decrease in    Responses of Global Terrestrial NPP to Climate Change and Rising CO 2 PLOS ONE | www.plosone.org warming temperature and likely drought events. The climate-only experiment under the B1 scenario, however, shows no substantial decrease in terrestrial NPP in tropical regions. Our results further show that in a particular year, drought or rainfall deficits combined with increasing temperature can potentially reduce terrestrial NPP. Compared to the average precipitation of the 21 st century, terrestrial ecosystems experienced the largest reduction in global precipitation (Figure 10) of about 68.97 and 43.07 mm yr 21 in 2019 and 2020 under the A2 and B1 scenario, respectively. This annual precipitation deficit resulted in a reduction in global terrestrial NPP by 375 gC m 22 in the tropics and low latitude regions ( Figure 10).

Latitudinal and spatial responses of terrestrial NPP
to climate change and increasing CO 2 concentration. Climate plus CO 2 experiment under both the A2 and B1 scenarios show substantial difference in the magnitude and temporal pattern of terrestrial NPP along latitudes ( Figure 7B-D). The A2 climate plus CO 2 experiment shows a substantial increase in terrestrial NPP in low latitude regions until the 2060s where NPP would increase by 4.87 PgC which is equivalent to an increase of 13.5%. However, terrestrial NPP would start to decline in the 2090s (Table 1). In mid-and high-latitude regions, terrestrial NPP would continue to increase through the 21 st century, with an increase of about 22.0% and 51.9%, respectively by the end of the century (2090s vs. 2000s). The climate plus CO 2 experiment under the B1 scenario shows an increasing trend in terrestrial NPP across all latitudes where NPP would increase by 12.4%, 9.9%, and 19.9% in low-, mid-, and high-latitude regions, respectively. During the 2090 s, the largest magnitude (4.48 PgC) of increase in terrestrial NPP would occur in low latitude regions under the B1 scenario, while the largest rate (51.9%) of increase would occur in high latitude regions, under the A2 scenario.
Large increase in terrestrial NPP would occur in tropical regions especially in central and southern Africa, and the Amazon basin under the A2 scenario, while modest increase under the B1 scenario ( Figure 9B, 9D). The A2 scenario shows large increase in terrestrial NPP by .200 gC m 22 yr 21 , while the B1 scenario shows an increase of 100-200 gC m 22 yr 21 in tropical regions. Large increase in NPP in the tropical regions due to increasing atmospheric CO 2 concentration is primarily because DLEM uses a Farquhar model that shows a higher NPP enhancement at high temperatures under elevated CO 2 . However, the mid-and highlatitude regions show lower NPP enhancement in response to  Table 2). The A2 climate scenario would result in a NPP decrease for major biomes such as tropical broadleaf deciduous and tropical evergreen forest, evergreen shrubs, grasses, woody wetland and croplands in the 2090s. The largest percent decline would occur in tropical broadleaf evergreen forest where NPP would decrease by 104.0 gC m 22 , equivalent to a decrease of 9.3% compared to the contemporary NPP. The largest percent increase under the A2 scenario would occur in tundra and boreal needleleaf deciduous forest by 41.5 gC m 22 and 176 gC m 22 equivalent to an increase of 59.1% and 54.3%; respectively. The B1 climate scenario, however, shows no substantial change in terrestrial NPP across different biomes compared to the A2 scenario. We found the largest percent decline in NPP in woody wetland by 2.4%, and the largest percent increase in NPP in boreal needleleaf deciduous forest by 16.6% in the climate-only simulation under the B1 scenario. Croplands, in particular, show a decline in NPP due to climate variability under both the A2 and B1 scenarios with the largest decline of 9.9% under the A2 scenario. Simulations with climate change and increasing atmospheric CO 2 would result in an increase in tropical broadleaf evergreen forest NPP by 129.6 gC m 22 (which decreases with climate alone), equivalent to a percent increase of 11.5%, indicating that CO 2 ferilization could have a substantial effect on terrestrial NPP. The B1 scenario, however, shows a relatively small increase in terrestrial NPP (7.5-22.8%) across all biomes. Interestingly, croplands show the largest increase in NPP by 45.9 gC m 22 (9.4%) under B1 scenario when climate change and increasing atmospheric CO 2 concentration are included in an experiment. This increase is higher than that under the A2 scenario when climate change is coupled with increasing atmospheric CO 2 concentration.

The CO 2 fertilization effects on terrestrial NPP
Our results show that terrestrial NPP during the 2090s would increase due to increasing atmospheric CO 2 across the globe, but the magnitude of increase varied substantially across different latitudes ( Figure 11A-B). The largest increase in terrestrial NPP due to increasing atmospheric CO 2   Responses of Global Terrestrial NPP to Climate Change and Rising CO 2 Figure 10. Spatial variation in precipitation and NPP estimated as a difference between dry year and long term (2000-2099) mean: precipitation difference between 2019 and long term mean (A) for A2 scenario, NPP difference between 2019 and long term mean (B) for A2 scenario climate-only simulation, precipitation difference between 2020 and long term mean (C) for B1 scenario and NPP difference between 2020 and long term mean (D) for B1 climate-only simulation. doi:10.1371/journal.pone.0112810.g010 Table 2. Decadal mean of terrestrial net primary production (NPP) in the contemporary period (2000-2009) and change in NPP between the 2090s and the 2000s among major biomes. We further calculated the 'beta' (b) effect that measures the strength of changes in terrestrial NPP in response to increasing atmospheric CO 2 concentration to examine the relative contribution of direct CO 2 fertilization on terrestrial NPP (Figure 12). By the end of 21 st century, large increase in global terrestrial NPP would occur under the A2 scenario, while a small increase would occur under the B1 scenario. At the global scale, the effect of direct CO 2 fertilization was 91.4 mgC m 22 /CO 2 (ppm) under the A2 scenario, implying that about 91.4 mgC m 22 would be fixed by plants as NPP by using 1 ppm of atmospheric CO 2 by the 2090s. The effect of CO 2 fertilization, however, would start to decline after the 2080 s, with a net reduction in stimulative effect by 8%.
Across different latitudes, the effect of CO 2 fertilization is highest in low latitudes (127.87 mgC m 22 /CO 2 (ppm)), and lowest in high latitude (15.52 mgC m 22 /CO 2 (ppm)). Although a higher effect of CO 2 feritlization is found under the A2 scenario in low latitude regions, this effect would decline after the 2070s with a net reduction in strength by 12.63 mgC m 22 during the 2090s compared to the 2070s. Interestingly, the strength of CO 2 fertilization would increase continuously under the B1 scenario in low latitude regions through the 21 st century. The strength of CO 2 fertilization is lower in mid-and high-latitude regions compared to low latitude regions where NPP shows a continuous increase in response to CO 2 fertilization through the 21 st century.

Comparison of DLEM-simulated NPP with previous estimates
The DLEM-simulated global terrestrial NPP is 54.57 PgC yr 21 for the period 2000-2009, which is comparable to MODIS-based estimate of 53.5 PgC yr 21 during the 2000s [16], and also falls in the range of 44-66 PgC yr 21 as estimated by 17 global terrestrial biosphere models [29]. The large variation in global NPP estimates among these models are primarily caused by different model representations of nutrient and water constraints on NPP in various terrestrial ecosystems. For instance, both field-and modelbased studies indicated that carbon and nitrogen are closely interacted in many terrestrial ecosystems [60]. Some of these models do not simulate carbon-nitrogen interaction, which generally result in an overestimate of global NPP. In addition, moisture in root zone is critically important for plants especially in areas under frequent water stress [65,66]. Remote sensing algorithms, uses a nonlinear function of increasing temperature to estimate production responses to water limitation, which is not capable of accurately accounting for environmental stresses such as rooting depths. Such function can enhance heat stress greatly in low latitude regions, introducing large uncertainties in NPP estimates.
Our results indicate that high latitude biomes have the greatest potential to increase NPP in response to future climate and CO 2 changes. Mainly due to climate change effect, the boreal needleleaf deciduous forest and tundra would have the largest NPP increase by 79.5% and 78.4% in the 21 st century, respectively ( Table 2). Hill and Henry [67] found that aboveground biomass of tundra sedge community increased on average by 158% in 2005 compared to the 1980s primarily due to climate warming. Our climate dataset indicates boreal forests and tundra ecosystems would experience 3-6uC temperature increase in the 21 st century ( Figure 4). In contrast, the model predicted a decline in NPP in tropical forests in the climate-only experiment under the A2 scenario that predicted rapid temperature increase in the 21 st century. This prediction is supported by Clark et al. [68]'s long-term study in tropical rain forest, which suggests that temperature increase in tropical regions would suppress forests NPP.

Climate change effects on global NPP in the 21 st century
Our climate-only scenario simulations indicate that temperature rise in the 21 st century would first promote and then reduce global terrestrial NPP. We found that under the A2 scenario the turning point of global NPP would occur when the global mean temperature reaches about 16.5uC. As temperature crosses 16.5uC, the terrestrial biosphere would show a net reduction in NPP compared to the contemporary period ( Figure 8; left panel). Our results also show a levelling-off of terrestrial NPP when the global mean temperature reaches about 15uC under both scenarios. Compared with the contemporary global temperature of ,13.5uC, our finding justifies the goal of maintaining global warming rate under 2uC by the Copenhagen Accord [69]. However, the temperature value of 15uC should not be treated as a threshold to judge climate effects on ecosystem NPP at regional scale, due to the complex responses of ecosystems to climate change in different regions [70]. Our analysis shows that NPP would increase rapidly in high latitude region (e.g., by 34% under the A2 scenario) while decrease in low latitude region in response to climate warming ( Figure 9A). The decline in terrestrial NPP in low latitudes is possibly the result of enhanced autotrophic respiration and increased moisture stress caused by rising temperature [27,29,30]. Warming may cause increasing moisture stress resulting in a decline in net nitrogen mineralization, therefore lead to reduction in NPP [25]. In contrast, temperature increase in mid-and high-latitudes can stimulate NPP by alleviating low temperature constraints to plant growth [10], lengthening growing season [13,71], and enhancing nitrogen mineralization and availability [72].
More extreme precipitation regimes are expected to have a substantial effect on terrestrial NPP during the 21 st century, and therefore ecological implications of greater intra-annual variability and extremes in rainfall events have received much attention from the scientific community [73,74]. Our study suggests that precipitation deficit would result in reduction in terrestrial NPP by 375 gC m 22 across eastern United States, South America, Africa, Europe and Southeast Asia under the A2 scenario. The B1 scenario, however, shows a reduction of up to 175 gC m 22 with no net change in most of the areas across the globe. Such reduction in terrestrial NPP due to decrease in precipitation has been reported worldwide [16,75,76]. Although precipitation is projected to increase during the 21 st century under the A2 and B1 scenarios, their spatio-temporal variability is a major factor that determines the magnitude of terrestrial NPP. For instance, we observed a reduction in global precipitation by 8% in 2019 that led to a decline in terrestrial NPP of about 375 gC m 22 dominated largely in the tropics. Mohamed et al. [77] found larger NPP decline in tropical regions indicating greater impacts of drought and high sensitivity and weaker physiological adjustment to climate variability particularly precipitation. Such variability in precipitation will induce more decline of NPP in the tropics. In addition, frequent extreme drought may counteract the effects of anticipated warming and growing season extension and reduce terrestrial production [76] in mid-and high-latitude reigons. Increase in drought events is primiarly the result of less frequent but more intense precipitation events which would likely decrease NPP in wetter ecosystems but increase NPP in drier ecosystems [73]. Therefore, a deeper understanding of the importance of more extreme precipitation patterns relative to other global change drivers such as increasing atmospheric CO 2 concentration, elevated temperature and atmospheric pollution is critical to project future NPP in a changing global environment.

CO 2 fertilization effect and its interaction with climate change
While climate change may have negative (24.6% under the A2 scenario) to small positive (2.2% under the B1 scenario) effect on global terrestrial NPP, our Climate plus CO 2 simulations project notable increase (12% under the B1 scenario to 13.9% under the A2 scenario) in NPP by the end of 21 st century. Majority of laboratory studies of tree growth under high CO 2 concentration show enhanced growth rates, on average, with a 60% increase in productivity as a result of doubling of atmospheric CO 2 [18]. This explains the prominent CO 2 fertilization effect found in low latitude, especially under the A2 scenario ( Figure 12; Table 1 and 2), which would more than compensate the negative effect (2 76.1-2104 gC m 22 yr 21 ) of climate change and result in an increase in NPP (132.2-233.6 gC m 22 yr 21 ) in tropical forest during the 21 st century. Elevated CO 2 was also found to reduce stomatal conductance and increase water use efficiency in waterlimited environments [78]. Our study indicates that CO 2 fertilization effect might help the Mediteranian evergreen shrubland to overcome severe NPP loss (215.1%) in response to the projected future drought (Figure 4) due to the poleward shift of subtropical dry zones [79], and resulted increase in shrubland NPP (6%) by the end of the 21 st century under the A2 scenario (Table 2).
However, field studies also found large uncertainties related to terrestrial ecosystems' responses to CO 2 enrichment. Plant response to elevated atmospheric CO 2 can be modified by increasing temperature [80], soil nutrient deficiency [81], and tropospheric ozone pollution [82]. Temperature rise in the 21 st century could enhance mineralization rate of soil nutrients, which would support higher ecosystem productivity. In low latitudes, however, nutrient leaching rate will also increase due to intensive and increasing precipitation. Our analyses indicate that at the beginning of 21 st century, relative effect of CO 2 fertilization on terrestrial NPP would increase fast; then as a result of CO 2 saturation and gradually intensified water and nutrient limitations, the effect of CO 2 fertilization would reach its maximum potential in the 2070s, level off and even decline afterwards ( Figure 12).
While climate change would result in large increase in terrestrial NPP in mid-and high-latitude regions due to greater growing season extension compared to low-latitude region [83], the effect of increasing CO 2 concentration on plant growth will be stronger in low-latitude region during the 21 st century. This is primarily because of much stronger photosynthetic response under elevated CO 2 at high temperatures [84]. The DLEM uses a Farquhar model to examine the response of terrestrial primary production to climate change and increasing atmospheric CO 2 concentration. The Farquhar model may cause a much stronger CO 2 enhancement at high temperature compared to low temperatures, resulting in higher NPP response to increasing CO 2 concentration in low latitude compared to mid-and high-latitude regions. Hickler et al. [85] evaluated the Farquhar photosynthesis model and found that the optimum temperature for primary production shifts to higher temperatures under elevated CO 2 . The DLEM simulates a stronger NPP response to elevated CO 2 under drier environments because elevated CO 2 reduces the negative effect of drought on plant growth resulting in higher NPP response in low-latitude region. However, mid-and high-latitude regions are primarily thought to be nitrogen limited. Increase in temperature and precipitation in mid-and high-latitude regions would possibly enhance decomposition of soil organic matter resulting in more nitrogen available for plant uptake. The resultant increase in plant nitrogen uptake with further benefits from elevated CO 2 would result in an increase in NPP in mid-and high-latitude ecosystems such as eastern Europe, eastern Russia, and parts of China ( Figure 9B, 9D). Our regional variation in the response of terrestrial NPP to rising CO 2 concentration are similar to Kirschbaum [86] who used a modified Farquhar model to simulate the response of photosynthesis to elevated CO 2 .

Uncertainty and future research needs
Several limitations and uncertainties are inherent in this study regarding input data, model parameterization and simulations. Our goal is to investigate impacts of future climate change and elevated CO 2 concentration on the global terrestrial NPP. We applied one GCMs (CCSM3) climate output together with atmospheric CO 2 data under two emission scenarios (A2 and B1). Discrepancies existing in different global climate models would lead to different estimation of NPP response as projected climate variables change [87]. Further analysis is helpful to explore the uncertainty ranges by adopting climate projections derived from multiple climate models [29]. In addition, terrestrial NPP response to increasing CO 2 concentration in the model is based on calibration and parameterization at several Free-Air Concentration Enrichment (FACE) sites based on Ainsworth and Long [88]. Because current long-running FACE experiments are all located in the temperate zone [89], we have very little empirical information available on the response of NPP to elevated CO 2 in the tropics and high-latitude ecosystems. While parameters were well calibrated based on existing field observations, some processes such as responses of C assimilation/allocation and stomatal conductance to elevated temperature and CO 2 may change due to plant acclimation [90] and dynamic responses of phenology and growing season length [91], which have not been included in the current DLEM simulations. In addition, we only attempted to quantify the terrestrial ecosystem response to seasonal and interannual climatic variability, and increasing atmospheric CO 2 concentration, but did not consider how climate induced functional change in ecosystem processes [92] may affect terrestrial NPP. We should be aware of how other environmental factors (e.g. nitrogen deposition, tropospheric ozone, and LCLU change and land management, fertilization and irrigation of croplands) may work together with climate change and elevated CO 2 to influence terrestrial NPP. Also, disturbance such as timber harvest and cropland abandonment may have a substantial effect on C dynamics at the regional scale [93]. Therefore, future research must take into account additional environmental and human factors such as nitrogen deposition, ozone pollution and land use/land cover change. Furthermore, model representations of ecosystem processes associated with human activities must be improved.

Conclusion and Implications for Climate Change Policy
The DLEM simulations estimate a mean global terrestrial NPP of 54.57 Pg C yr 21 in the first decade of the 21 st century, resulting from multiple environmental factors including climate, atmospheric CO 2 concentration, nitrogen deposition, tropospheric ozone, and LCLU change. Climate change in the 21 st century would either reduce global NPP by 4.6% under the A2 scenario or slightly enhance NPP (2.2%) under the B1 scenario. In response to climate change, global NPP would first increase and then decline after global warming exceeds 1.5uC, which justifies the goal of keeping the global warming rate under 2uC by the Copenhagen Accord (2009). These results also support the Accord Acknowledgement, which states that staying below 2uC may not be sufficient and a review in 2015 will be conducted to assess the need to potentially aim for staying below 1.5uC. Terrestrial NPP at high latitude (60uN-90uN) in the Northern Hemisphere would benefit from climate change, but decline at low latitude (30uS-30uN) and the central USA. The CO 2 fertilization effect would improve global productivity as simulated by the DLEM forced by climate plus CO 2 , showing notable increase in NPP by the end of 21 st century (12% under the B1 scenario to 13.9% under the A2 scenario). Ecosystem responses to increasing CO 2 , however, are complex and involve large uncertainties, especially in the tropics.
The relative CO 2 fertilization effect, i.e. change in NPP on per CO 2 (ppm) bases, was projected to increase quickly, level off in the 2070s, and then even decline by the end of the 2080s, possibly due to CO 2 saturation and nutrient limitation. Compared to the low emission scenario (B1), earth ecosystems would have a less predictable and thus unfavorable future under the high emission scenario (A2) due to large uncertainties related to the global NPP, which are characterized by stronger spatial variation, more complex temporal dynamics, larger differences in biomes' responses to global changes, and more dramatic and counteracting impacts from climate and elevated CO 2 .
Our model projections indicate that high emission scenario (A2) will not necessarily negatively affect global terrestrial NPP in the future. However, the A2 scenario describes a future with high ecological uncertainty, the rest of 21 st century may experience frequent climate extreme events and increasing ecological risks, and a global ecosystem whose behavior and trend are difficult to predict and therefore can hardly be protected with efficiency. The world under the A2 scenario would likely face more serious challenges in food security and water scarcity given a continually growing world human population. To avoid this scenario, we should shift to the low emission scenario.