This is the accepted version of the article: Fernández Martínez, Marcos, et al. Global trends in carbon sinks and their relationships with CO₂ temperatura in Nature climate change, vol. 9 (Jan. 2019), p. 73-79. Available at: https://dx.doi.org/10.1038/s41558-018-0367-7 Published under a “All rights reserved” license. 1 Global trends in carbon sinks and their relationships with CO2 and 2 temperature 3 Authors: M. Fernández-Martínez*1. J. Sardans2,3, F. Chevallier4, P. Ciais4, M. Obersteiner5, S. 4 Vicca1, J. G. Canadell6, A. Bastos4, P. Friedlingstein7, S. Sitch7, S.L. Piao8,9, I.A. Janssens1, J. 5 Peñuelas2,3. 6 Affiliations: 7 1 Centre of Excellence PLECO (Plant and Vegetation Ecology), Department of Biology, 8 University of Antwerp, 2610 Wilrijk, Belgium. 9 2 CSIC, Global Ecology Unit, CREAF-CSIC-UAB, Cerdanyola del Vallès, 08193 Barcelona, 10 Catalonia, Spain 11 3 CREAF, Cerdanyola del Vallès, 08193 Barcelona, Catalonia, Spain 12 4 Laboratoire des Sciences du Climat et de l’Environnement, CEA CNRS UVSQ, 91191 Gif-sur- 13 Yvette, France 14 5 International Institute for Applied Systems Analysis, Schlossplatz 1, 2361 Laxenburg, Austria 15 6 Global Carbon Project, CSIRO Oceans and Atmosphere, Canberra, ACT 2601, Australia 16 7 College of Engineering, Computing and Mathematics, University of Exeter, Exeter EX4 4QF, 17 UK 18 8Sino-French Institute of Earth System Sciences, College of Urban and Environmental 19 Sciences, Peking University, Beijing 100871, China 20 9 Institute of Tibetan Plateau Research, Chinese Academy of Sciences, Beijing 100085, China 21 *Correspondence to: M. Fernández-Martínez, marcos.fernandez-martinez@uantwerpen.be 22 23 24 Elevated CO2 increases photosynthesis and, potentially, net ecosystem production 25 (NEP) meaning greater CO2 uptake. Climate, nutrients, and ecosystem structure, 26 however, influence the effect of increasing CO2. Here, we analysed global NEP from 27 MACC-II and Jena CarboScope atmospheric-inversions and 10 dynamic global 28 vegetation models (TRENDY), using statistical models to attribute the trends in NEP to 29 its potential drivers: CO2, climatic variables and land-use change. We find that 30 increasing CO2 was consistently associated with increased NEP (1995-2014). 31 Conversely, increasing temperatures were negatively associated with NEP. Using the 32 two atmospheric inversions and TRENDY, the estimated global sensitivities for CO2 33 were 6.0 ± 0.1, 8.1 ± 0.3 and 3.1 ± 0.1 Pg C per 100 ppm (~1 ºC increase), and -0.5 ± 34 0.2, -0.9 ± 0.4 and -1.1 ± 0.1 Pg C °C-1 for temperature. These results indicate a 35 positive CO2 effect on terrestrial C sinks that is constrained by climate warming. 36 37 In recent decades, terrestrial ecosystems have been absorbing 15–30% of all 38 anthropogenic CO2 emissions1,2. Direct and indirect anthropogenic impacts on the 39 biosphere, however, can alter terrestrial sinks in the short and long terms3–6. Identifying 40 the factors that affect the capacity of the biosphere to absorb carbon (C) and 41 quantifying the magnitude of the sensitivity of this C sink to its driving factors helps to 42 increase confidence in future projections of the coupled C cycle/climate system. 43 Increasing plant growth is a robust response to increasing CO2 concentrations under 44 experimental conditions (CO2 fertilization effect)7,8. The extent to which increases in 45 CO2 can enhance large-scale photosynthesis and ultimately net ecosystem production 46 (NEP) remains uncertain5,7. Detecting this effect in the real world is much more difficult 47 than under controlled experiments. However, recent efforts using eddy-covariance- 48 based data and statistical models have been successful in detecting positive effects of 49 CO2 on water-use efficiency (WUE)9, photosynthesis, and NEP5. 50 The potential positive effect of elevated CO2 on productivity could be influenced by 51 global warming6 and altered precipitation patterns10 since both water availability and 52 temperature are strong drivers of photosynthesis and respiration worldwide11–13. Land- 53 use change also alters the capacity of the biosphere to sequester C because land use 54 causes a drastic change in C turnover and productivity. Atmospheric deposition of 55 nitrogen (N) and sulphur (S) from the use of fossil fuels and fertilisers may also alter 56 ecosystem biodiversity, function, productivity and NEP5,14–17. N deposition is usually 57 positively correlated with ecosystem productivity and NEP17–19. Conversely, S 58 deposition may reduce ecosystem carbon sinks, this has rarely been investigated in 59 field studies20,21 and absent from global models. Soil acidification, caused by acid 60 deposition, of N and S, often decreases the availability of soil nutrients22 and potentially 61 reduces NEP23. 62 The observations underlying the driver analysis of NEP described above were largely 63 limited to temperate and boreal study sites, making it difficult to assess global 64 scalability. Additionally, until recently, the only way to assess terrestrial C sink was from 65 ensembles of dynamic global vegetation models (DGVMs) or as a residual sink, by 66 subtracting atmospheric and ocean sinks to the estimates of CO2 emissions. Currently, 67 inversion models, as well as long-term remotely sensed data24, can be used to test the 68 generality of the patterns derived from ground-based measurements. Inversion models 69 provide continuous gridded estimates for the net flux of land-atmosphere CO2 70 exchange (i.e. NEP) with global coverage25,26. The gridded NEP results from 71 inversions, combined with CO2-concentration records, gridded fields for climate, land- 72 use change, and atmospheric deposition, are arguably the best observation-based data 73 to attempt a first empirical study of the combined effects of CO2, changes in climate 74 and land use, and atmospheric N and S deposition on terrestrial NEP patterns at the 75 global scale. Given that previous site level studies revealed that increasing CO2 is a 76 dominant driver of trends in NEP, we expect that it will also be the dominant driver at 77 larger spatial scales and across the globe. 78 Here we investigate if the trends of NEP from the two most widely used multi-decadal 79 inversion models (MACC-II and Jena CarboScope) and DGVMs (TRENDY) from 1995 80 to 2014 are related to increasing atmospheric CO2 and changing climate (temperature, 81 precipitation, and drought). Additionally, the effect of land-use on NEP at the global 82 scale was investigated using statistical models to assess the sensitivity of NEP to the 83 abovementioned predictors. We also analysed the effect of changing rates of 84 atmospheric deposition of oxidised and reduced N and S on NEP, combined with 85 increasing CO2 and changing climate and land use, over Europe and the USA. 86 CO2 and climate effects on global NEP 87 Global land (excluding Antarctica) mean annual NEP was 2.3 ± 0.9, 2.3 ± 1.5 and 1.6 ± 88 0.5 Pg C y-1 (mean ± 1σ), respectively, for MACC-II, Jena CarboScope and the 89 TRENDY ensemble during the period 1995–2014, similar in magnitude to the recent 90 global carbon budget2. Both inversions and the TRENDY ensemble showed an overall 91 positive trend in NEP from 1995 to 2014. The estimated NEP increased by (mean ± 92 1SE) 116.9 ± 6.1 Tg C y-1 for the MACC-II dataset, by 178.0 ± 8.1 Tg C y-1 for the Jena 93 CarboScope dataset, and by 22.5 ± 3.1 Tg C y-1 for the TRENDY ensemble (Figure 1). 94 This supports the increases in the global carbon budget2, with a lower increase of the 95 DGVMs than those shown by the inversion models. The large differences between 96 inversion models and DGVMs may arise because of the lack of information on river 97 fluxes, inadequate parameterisations concerning land management and degradation in 98 the process models or because of potential biases in inversion models. Both inversion 99 model datasets produced similar trends for many parts of the world, an increasing NEP 100 for Siberia, Asia, Oceania, and South America, and a decreasing NEP for the southern 101 latitudes of Africa. Differences between inversions emerged for Europe and North 102 America, possibly because Jena CarboScope inversion uses a larger spatial error 103 correlation of prior fluxes than MACC-II or because of other inversion settings2. 104 However, their different flux priors did not drive differences in the trends between both 105 datasets, given that priors did not change over the studied period. Jena CarboScope 106 showed largely positive trends for Europe and largely negative trends for North 107 America; MACC II showed more variation in the trends for both continents. The trends 108 identified by the TRENDY ensemble agreed with atmospheric inversions for the 109 northernmost latitudes, indicating an increase in C-sink capacity, but differed from 110 those in many other regions. 111 Our analyses on temporal contributions, using the temporal anomalies of our 112 predictors, attributed the increases in global NEP to increasing CO2 but found a 113 consistent negative impact of temperature on NEP, which limited the positive effect of 114 increasing CO2 (Figure 1). These results were consistent for both datasets and most of 115 the DGVMs of the TRENDY ensemble. The predictors used in this study explained a 116 modest proportion of the variance in NEP, in contrast to the variance explained by 117 spatial variability (i.e., the pixel), which was rather high (Supplementary Information 118 (SI), Section 2). Unknown contributions to trends in NEP, the difference between all 119 contributions and the observed trend, were very close to zero for the analyses on 120 inverse models and the TRENDY ensemble (Figure 1). This result suggests that trends 121 were very well captured by our analyses, indicating that the methodology was able to 122 disentangle spatial from temporal variability. The sensitivity of NEP to increasing CO2 123 averaged 0.45 ± 0.01, 0.61 ± 0.03 and 0.23 ± 0.01 g C m-2 ppm-1 for MACC-II, Jena 124 CarboScope and TRENDY, respectively (Table 1), representing sensitivities over the 125 entire terrestrial surface of 60.4 ± 1.2, 81.4 ± 3.4 and 30.7 ± 1.2 Tg C ppm-1, 126 respectively. Despite lower temporal attributions for temperature than CO2, the 127 sensitivity of NEP to temperature was high, at -3.8 ± 1.1, -6.4 ± 2.9 and -8.1 ± 0.9 g C 128 m-2 y-1 °C-1 for the MACC-II, Jena CarboScope and TRENDY models, respectively, 129 equivalent to global sensitivities of -515.7 ± 152.4, -859.2 ± 386.3 and -1088.0 ± 118.1 130 Tg C °C-1, respectively. Trends in NEP and the effect of CO2 and temperature on NEP 131 significantly differed in magnitude amongst the datasets used, however, they all point 132 towards the same conclusion: global NEP has increased during the study period and 133 increasing CO2 has been the most likely driving factor despite increasing temperatures 134 constraining this positive effect. The exact magnitude of the effect of increasing CO2 135 and temperatures on global carbon cycle remains to be established 136 Spatial variability on CO2 and climate change effects on NEP 137 Our statistical models for the MACC-II and Jena CarboScope datasets indicated that 138 the positive effect of CO2 on NEP was higher in regions with higher annual precipitation 139 and that this positive effect increased with increasing temperatures (Figure 2, SI 140 Section 1.1). In contrast, our analyses using the TRENDY ensemble did not show a 141 significant interaction between CO2 and precipitation or with temperature, highlighting 142 the different behaviour in the DGVMs compared to inversion models. We also found a 143 significant positive interaction between mean annual temperature and CO2 for Jena 144 CarboScope and TRENDY. However, the same interaction was negative for MACC-II. 145 On the other hand, increasing temperatures reduced NEP in warm regions but 146 increased NEP in cold regions (Figure 2). 147 The analyses on temporal contributions performed for inversion and TRENDY NEP 148 averaged over latitudinal bands (boreal, >55°; temperate, 35-55°; subtropical, 15-35°; 149 and tropical, 15°N-15°S), further supported previous results obtained at the global 150 scale (Table 2, SI Sections 2.2–2.7). Increasing CO2 was the main factor accounting 151 for increasing trends in NEP, with a consistent positive temporal contribution for almost 152 all latitudinal bands considered and for all three datasets. However, contributions 153 estimated from the TRENDY ensemble were generally lower than those of the 154 inversion models. Proportionally, increasing CO2 accounted for more than 90% of the 155 trends in NEP in MACC-II and Jena CarboScope datasets. For the TRENDY ensemble, 156 the estimated contribution of CO2 to the trends in global NEP was more than 2.7 times 157 higher than the estimated trends. Increasing temperatures had a negative effect for all 158 latitudinal bands for the inversion models, but most effects were not statistically 159 significant and need to be interpreted as such. Instead, our analyses for the TRENDY 160 ensemble indicated a significant negative effect for all latitudinal bands, except for the 161 temperate southern hemisphere. Similarly, the proportional contribution of temperature 162 to the trends in NEP was less than 10% for the inversion models, but accounted for 163 almost 95% of the trends estimated using the TRENDY ensemble. These results 164 suggest that the parameterisation of temperature in the DGVMs does not accurately 165 reproduce the estimation of the inverse models. 166 Despite all regions presented, on average, positive trends, the tropical regions clearly 167 had the highest contribution, across all three datasets, to global NEP trends accounting 168 for almost half of the increase (Table 2). Similarly, the tropical regions had the highest 169 sensitivity to CO2 increase, accounting for more than half of the total global sensitivity 170 (Table 1). A similar pattern was found for temperature, although the sign of the 171 contribution was positive for MACC-II but negative for Jena CarboScope and TRENDY. 172 The contribution of the southern hemisphere to the global trends was very modest 173 compared to the contribution of the northern hemisphere using all datasets. Our results 174 using the MACC-II dataset showed that subtropical, temperate and boreal regions of 175 the northern hemisphere accounted for 44.2% of the global trends in NEP, while only 176 9.5% was attributed to subtropical and temperate regions of the southern hemisphere. 177 Using the Jena CarboScope dataset these regions accounted for 63.3% and 6.1%, 178 respectively. Differences on the regional attributions between inversion models may 179 emerge from the different interhemispheric transport models or other inversion 180 settings2. Results from the TRENDY ensemble were more extreme, because they 181 indicated a negative contribution of the subtropical and temperate regions to the global 182 trends in NEP. Differences between the global estimates (trends and contributions of 183 CO2 and temperature) and the sum of every region were low for all datasets. 184 Contribution of other variables to the trends in NEP (precipitation, drought, land-use 185 change, and unknown variables) were on average also low for most of the latitudinal 186 bands, despite the variability amongst datasets (Table 2). 187 Atmospheric deposition 188 The MACC-II and Jena CarboScope datasets showed that NEP increased over Europe 189 and the USA by 0.45 ± 0.13 and 0.68 ± 0.16 g C m-2 y-1, respectively (Figure S1). Our 190 temporal contribution analyses suggested that increasing atmospheric CO2 in both 191 datasets contributed significantly to increasing NEP. NEP sensitivity to CO2 was more 192 than two-fold higher in the Jena CarboScope than the MACC-II dataset (Table S1), 193 similar to the temporal contributions, at 0.22 ± 0.06 and 0.46 ± 0.07 g C m-2 y-1 ppm-1 194 for the MACC-II and Jena CarboScope models, respectively. The temporal contribution 195 of decreasing NOX deposition to NEP differed between the two datasets; the 196 contribution was positive for MACC-II and negative for Jena CarboScope. Our analyses 197 consequently estimated a negative sensitivity of NEP to NOX for the MACC-II dataset 198 but a positive sensitivity for the Jena CarboScope dataset. Additionally, neither MACC- 199 II, nor Jena CarboScope indicated a strong impact of land use change. 200 These statistical models indicated that, in both datasets, the positive effect of CO2 on 201 NEP was higher in regions with higher NRED deposition but lower in regions with high S 202 deposition (means for MACC-II and annual anomalies for Jena CarboScope; see SI 203 section 2.8). The results for NOX deposition, however, differed between the models. 204 The positive effect of CO2 on NEP for the MACC-II dataset was constrained by the 205 annual anomalies of NOX but was higher for the Jena CarboScope dataset. We also 206 estimated an overall negative but not significant sensitivity of NEP to S deposition for 207 both inversion models. 208 CO2 fertilisation and global NEP 209 The positive effect of atmospheric CO2 on NEP must originate from a stronger positive 210 effect on photosynthesis than on the sum of all respiratory processes. Increasing 211 atmospheric CO2 concentrations have been widely reported to increase ecosystem 212 photosynthesis, mainly by two mechanisms: i) increasing carboxylation rates and 213 decreasing photorespiration27, and ii) decreasing stomatal conductance and therefore 214 increasing WUE9,28, which would theoretically increase photosynthesis under water 215 limitation. An increase in GPP by either mechanism may thus account for the higher 216 NEP due to increasing atmospheric CO2. A recent global analysis suggested that most 217 of the GPP gains from CO2 fertilization are associated with ecosystem WUE29. The 218 positive interaction between CO2 and annual precipitation that we found may not 219 support this hypothesis (Figure 2), given that plants living under wet conditions are 220 usually less efficient in water use. However, plants having higher water availability may 221 benefit from increasing CO2 more than those suffering drought because photosynthesis 222 would not be water-limited. 223 Our estimates of global NEP sensitivity to CO2 were 0.45 ± 0.01, 0.61 ± 0.03 and 0.23 224 ± 0.01 g C m-2 ppm-1 (globally 60.4 ± 1.2, 81.4 ± 3.4 and 30.7 ± 3.4 Tg C ppm-1) for the 225 MACC-II, Jena CarboScope and TRENDY datasets, respectively, but these estimates 226 varied amongst the latitudinal bands and were inconsistent between datasets (Table 227 1). These estimates were similar to those reported in CO2-enrichment FACE 228 experiments30, despite the fact that FACE values were calculated for a much higher 229 CO2 range for which the effect of CO2 may saturate31. However, they were much lower 230 than the 4.81 ± 0.52 g C m-2 ppm-1 reported in a study using eddy-covariance flux 231 towers for a similar period5. The much larger areas analysed by the inverse models 232 than the footprints covered by the eddy-covariance flux towers, and FACE 233 experiments, may explain these differences between the estimates. Flux towers are 234 usually located in relatively homogenous, undisturbed ecosystems, while each pixel in 235 the inverse model aggregates information from several ecosystems (and even biomes), 236 often including non-productive land such as bare soil or cities. 237 Our results indicated that the variability of the estimates of NEP sensitivity to CO2 238 amongst the latitudinal bands might be associated with differences in climate and 239 atmospheric N and S deposition. The two atmospheric inversion models indicated that 240 the effect of CO2 fertilisation was stronger in wet climates (high annual precipitation) 241 (Figure 2), supporting the estimates provided by the latitudinal bands, with the highest 242 sensitivity estimates for the tropical band (Table 1). However, analyses based on the 243 TRENDY ensemble did not show the same results. The positive effect of CO2 tended to 244 increase with temperature anomalies in both inversion models, but, again, the DGVMs 245 did not show the same behaviour. These differences between inversion models and 246 process-based models suggest that DGVMs still fail to capture some of the interactions 247 occurring in nature. The MACC-II and Jena CarboScope datasets further agreed on a 248 stronger positive effect of increasing CO2 in regions with higher NRED deposition, which 249 confirms previous studies suggesting that the effect of CO2 fertilisation is stronger in 250 nitrogen-rich sites32–34. 251 Climate, land-use and C sinks 252 Climatic warming clearly had a secondary effect on the trends in NEP from 1995 to 253 2014. The MACC-II, Jena CarboScope and TRENDY datasets estimated that NEP 254 decreased globally by around -0.5 ± 0.2, -0.9 ± 0.4 and -1.1 ± 0.1 Pg C for every 255 degree of increase in the Earth’s temperature. Assuming that a CO2 increase of 100 256 ppm is equivalent to an increase of global temperature of 1 ºC, the effect of the 257 increasing CO2 concentrations largely outweighs the negative effect of increasing 258 temperature on NEP (global estimates: 6.0 ± 0.1, 8.1 ± 0.3 and 3.1 ± 0.1 Pg C for a 259 100 ppm of CO2 increase according to MACC-II, Jena CarboScope and TRENDY). The 260 difference, though, is much lower for TRENDY than for the inversion models, having a 261 higher negative impact of temperature and a lower positive effect of CO 2. This 262 difference in the effects of temperature and CO2 may explain the lower trends observed 263 in TRENDY datasets compared to MACC-II and Jena CarboScope. It also suggests 264 that a different parameterisation of temperature, CO2 and their interaction may be 265 needed on DGVMs to capture the observed trends in the inversion models. 266 The quasi monotonically increasing atmospheric CO2 concentrations have been more 267 important than temperature in driving NEP trends. Increasing temperature, however, 268 did not have the same effect on NEP around the world. The analyses of both inverse 269 models indicated that increasing temperatures had a positive effect on NEP only in cold 270 regions (when MAT ≤ 1.5, 9 and -5.9 ºC for MACC-II and Jena CarboScope and 271 TRENDY respectively, when CO2 = 400 ppm, see SI section 2.1, and Figure 2). 272 These findings support previous literature reporting a positive effect between 273 temperature increase and NEP in temperate and boreal forests35. Instead, the general 274 negative effect of temperature on NEP could be due to a greater stimulation of Re than 275 photosynthesis by higher temperatures36,37. The potential benefit to C sequestration of 276 increased photosynthesis would then be negated by a greater increase in Re. 277 Increasing temperatures can also be linked to heat waves and drier conditions, which 278 may decrease GPP more than Re38. 279 The effects of land-use change on NEP trends differed greatly amongst the datasets, 280 both at the global scale and when using latitudinal bands. Our statistical models 281 identified several significant relationships between NEP and land-use change, but the 282 large differences in effects (direction and magnitude) amongst the datasets preclude 283 drawing firm conclusions. The coarse resolution of analysis likely blurred the effects of 284 land-use change on the NEP trends. 285 Our study highlights the dominant role of rising atmospheric CO2 concentrations 286 triggering an increase in land C sinks over the entire planet from 1995 to 2014, with the 287 tropics accounting for around half of this increase in NEP despite being only around 288 22% of the global land (excluding Antarctica, Table 2). Therefore, preserving tropical 289 ecosystems should be a global priority in order to mitigate anthropogenic CO2 290 emissions. Temperature has diminished the capacity of terrestrial ecosystems to 291 sequester C, which jeopardises future C sink capacity in light of global warming. So far, 292 our results suggest that the benefit of increasing atmospheric concentrations of CO 2 293 are still compensating the negative ones of temperature rise, in terms of C 294 sequestration. However, if it has not started to change already6, this pattern may 295 eventually reverse with saturation of land C sinks5,31 or because warm ecosystems tend 296 to decrease NEP as temperature rises (Figure 2). Additionally, the comparison 297 between model results indicated that the DGVMs were unable to reproduce several 298 features of the global land C sinks observed in inversion models. Process-based earth 299 system models will need to improve their parameterisation to capture these features in 300 order to better predict the future of land C sinks. 301 302 References: 303 1. Canadell, J. G. et al. Contributions to accelerating atmospheric CO2 growth from 304 economic activity, carbon intensity, and efficiency of natural sinks. Proc. Natl. 305 Acad. Sci. U. S. A. 104, 18866–70 (2007). 306 2. Le Quéré, C. et al. Global Carbon Budget 2017. Earth Syst. Sci. Data 10, 405– 307 448 (2018). 308 3. Ciais, P. et al. Europe-wide reduction in primary productivity caused by the heat 309 and drought in 2003. Nature 437, 529–533 (2005). 310 4. Crowther, T. W. et al. Quantifying global soil carbon losses in response to 311 warming. Nature 540, 104–108 (2016). 312 5. Fernández-Martínez, M. et al. Atmospheric deposition, CO2, and change in the 313 land carbon sink. Sci. Rep. 7:9632, 1–13 (2017). 314 6. Peñuelas, J. et al. Shifting from a fertilization-dominated to a warming dominated 315 period. Nat. Ecol. Evol. 1, 1438–1445 (2017). 316 7. Ainsworth, E. A. & Long, S. P. What have we learned from 15 years of free-air 317 CO2 enrichment (FACE)? A meta-analytic review of the responses of 318 photosynthesis, canopy properties and plant production to rising CO2. New 319 Phytol. 165, 351–71 (2005). 320 8. Medlyn, B. E. et al. Using ecosystem experiments to improve vegetation models. 321 Nat. Clim. Chang. 5, 528–534 (2015). 322 9. Keenan, T. F. et al. Increase in forest water-use efficiency as atmospheric 323 carbon dioxide concentrations rise. Nature 499, 324–327 (2013). 324 10. Alexander, L. et al. Climate Change 2013: The Physical Science Basis - 325 Summary for Policymakers. Fifth Assessment Report (Intergovernmental Panel 326 on Climate Change, 2013). 327 11. Fernández-Martínez, M. et al. Spatial variability and controls over biomass 328 stocks, carbon fluxes and resource-use efficiencies in forest ecosystems. Trees, 329 Struct. Funct. 28, 597–611 (2014). 330 12. Beer, C. et al. Terrestrial gross carbon dioxide uptake: global distribution and 331 covariation with climate. Science (80-. ). 329, 834–8 (2010). 332 13. Luyssaert, S. et al. CO 2 balance of boreal, temperate, and tropical forests 333 derived from a global database. Glob. Chang. Biol. 13, 2509–2537 (2007). 334 14. de Vries, W. & Posch, M. Modelling the impact of nitrogen deposition, climate 335 change and nutrient limitations on tree carbon sequestration in Europe for the 336 period 1900–2050. Environ. Pollut. 159, 2289–2299 (2011). 337 15. Wamelink, G. W. W. et al. Modelling impacts of changes in carbon dioxide 338 concentration, climate and nitrogen deposition on carbon sequestration by 339 European forests and forest soils. For. Ecol. Manage. 258, 1794–1805 (2009). 340 16. Wamelink, G. W. W. et al. Effect of nitrogen deposition reduction on biodiversity 341 and carbon sequestration. For. Ecol. Manage. 258, 1774–1779 (2009). 342 17. de Vries, W., Du, E. & Butterbach-Bahl, K. Short and long-term impacts of 343 nitrogen deposition on carbon sequestration by forest ecosystems. Curr. Opin. 344 Environ. Sustain. 9–10, 90–104 (2014). 345 18. Luyssaert, S. et al. The European carbon balance. Part 3: forests. Glob. Chang. 346 Biol. 16, 1429–1450 (2010). 347 19. Magnani, F. et al. The human footprint in the carbon cycle of temperate and 348 boreal forests. Nature 447, 848–50 (2007). 349 20. Thomas, R. B., Spal, S. E., Smith, K. R. & Nippert, J. B. Evidence of recovery of 350 Juniperus virginiana trees from sulfur pollution after the Clean Air Act. Proc. Natl. 351 Acad. Sci. U. S. A. 110, 15319–24 (2013). 352 21. Oulehle, F. et al. Major changes in forest carbon and nitrogen cycling caused by 353 declining sulphur deposition. Glob. Chang. Biol. 17, 3115–3129 (2011). 354 22. Truog, E. Soil Reaction Influence on Availability of Plant Nutrients1. Soil Sci. 355 Soc. Am. J. 11, 305 (1946). 356 23. Fernández-Martínez, M. et al. Nutrient availability as the key regulator of global 357 forest carbon balance. Nat. Clim. Chang. 4, 471–476 (2014). 358 24. Zhu, Z. et al. Greening of the Earth and its drivers. Nat. Clim. Chang. 6, 791–795 359 (2016). 360 25. Chevallier, F. et al. CO 2 surface fluxes at grid point scale estimated from a 361 global 21 year reanalysis of atmospheric measurements. J. Geophys. Res. 115, 362 D21307 (2010). 363 26. Rödenbeck, C., Houweling, S., Gloor, M. & Heimann, M. CO2 flux history 1982– 364 2001 inferred from atmospheric data using a global inversion of atmospheric 365 transport. Atmos. Chem. Phys. 3, 1919–1964 (2003). 366 27. Aber, J. et al. Forest Processes and Global Environmental Change: Predicting 367 the Effects of Individual and Multiple Stressors. Bioscience 51, 735 (2001). 368 28. Prentice, I. C., Heimann, M. & Sitch, S. The carbon balance of the terrestrial 369 biosphere: Ecosystem models and Atmospheric observations. Ecol. Appl. 10, 370 1553–1573 (2000). 371 29. Cheng, L. et al. Recent increases in terrestrial carbon uptake at little cost to the 372 water cycle. Nat. Commun. 8, 110 (2017). 373 30. Norby, R. J., Warren, J. M., Iversen, C. M., Medlyn, B. E. & McMurtrie, R. E. 374 CO2 enhancement of forest productivity constrained by limited nitrogen 375 availability. Proc. Natl. Acad. Sci. U. S. A. 107, 19368–73 (2010). 376 31. Norby, R. J. et al. Forest response to elevated CO2 is conserved across a broad 377 range of productivity. Proc. Natl. Acad. Sci. U. S. A. 102, 18052–18056 (2005). 378 32. Van Groenigen, K. J. et al. The Impact of Elevated Atmospheric CO2 on Soil C 379 and N Dynamics. Ecol. Stud. 187, 374–391 (2006). 380 33. Terrer, C. et al. Mycorrhizal association as a primary control of the CO₂ 381 fertilization effect. Science 353, 72–4 (2016). 382 34. McCarthy, H. R. et al. Re-assessment of plant carbon dynamics at the Duke 383 free-air CO2enrichment site: interactions of atmospheric [CO2] with nitrogen and 384 water availability over stand development. New Phytol. 185, 514–528 (2010). 385 35. Hyvönen, R. et al. The likely impact of elevated [CO2], nitrogen deposition, 386 increased temperature and management on carbon sequestration in temperate 387 and boreal forest ecosystems: a literature review. New Phytol. 173, 463–80 388 (2007). 389 36. Ryan, M. G. Effects of climate change on plant respiration. Ecol. Appl. 1, 157– 390 167 (1991). 391 37. Amthor, J. S. Scaling CO2 Photosynthesis Relationships from the Leaf to the 392 Canopy. Photosynth. Res. 39, 321–350 (1994). 393 38. Wu, Z., Dijkstra, P., Koch, G. W., Peñuelas, J. & Hungate, B. a. Responses of 394 terrestrial ecosystems to temperature and precipitation change: a meta-analysis 395 of experimental manipulation. Glob. Chang. Biol. 17, 927–942 (2011). 396 39. Chevallier, F. et al. Toward robust and consistent regional CO2 flux estimates 397 from in situ and spaceborne measurements of atmospheric CO2. Geophys. Res. 398 Lett. 41, 1065–1070 (2014). 399 40. Olivier, J. G. J. & Berdowski, J. J. M. in The Climate System (eds. Berdowski, J., 400 Guicherit, R. & Heij, B. J.) 33–78 (2001). 401 41. Sitch, S. et al. Recent trends and drivers of regional sources and sinks of carbon 402 dioxide. Biogeosciences 12, 653–679 (2015). 403 42. Harris, I., Jones, P. D. D., Osborn, T. J. J. & Lister, D. H. H. Updated high- 404 resolution grids of monthly climatic observations - the CRU TS3.10 Dataset. Int. 405 J. Climatol. 34, online, update (2013). 406 43. Vicente-serrano, S. M., Beguería, S. & López-Moreno, J. I. A Multiscalar Drought 407 Index Sensitive to Global Warming: The Standardized Precipitation 408 Evapotranspiration Index. J. Clim. 23, 1696–1718 (2010). 409 44. Zuur, A., Ieno, E., Walker, N., Saveliev, A. & Smith, G. Mixed effects models and 410 extensions in ecology with R. (Springer science, 2009). 411 45. Mathias, J. M. & Thomas, R. B. Disentangling the effects of acidic air pollution, 412 atmospheric CO2, and climate change on recent growth of red spruce trees in 413 the Central Appalachian Mountains. Glob. Chang. Biol. 24, 3938–3953 (2018). 414 46. R Core Team. R: A Lenguage and Environment for Stasitical Computing. (2016). 415 47. Barton, K. MuMIn: Multi-model inference. R package version 1.17.1. 416 http://CRAN.R-project.org/package=MuMIn. (2015). 417 48. Nakagawa, S. & Schielzeth, H. A general and simple method for obtaining R 2 418 from generalized linear mixed-effects models. Methods Ecol. Evol. 4, 133–142 419 (2013). 420 49. Breheny, P. & Burchett, W. Visualization of Regression Models Using visreg, R 421 package version 2.2-0. (2015). 422 Correspondence and requests for materials should be addressed to: 423 Marcos Fernández-Martínez: marcos.fernandez-martinez@uantwerpen.be 424 Acknowledgements 425 This research was supported by the Spanish Government project CGL2016-79835-P 426 (FERTWARM), the European Research Council Synergy grant ERC-2013-726 SyG- 427 610028 IMBALANCE-P, and the Catalan Government project SGR 2017-1005. M.F-M 428 and S.V. are a postdoctoral fellows of the Research Foundation – Flanders (FWO). 429 J.G.C. thanks the support of the National Environmental Science Programme ESCC 430 Hub. We thank Christian Röedenbeck for his advice and for distributing Jena 431 CarboScope and all the modellers that contributed to the TRENDY project. 432 433 Author Contributions 434 M.F-M., J.S., I.A.J., and J.P. conceived, analyzed and wrote the paper. F.C., P.F., and 435 S.S., provided data. All authors contributed substantially to the writing and discussion 436 of the paper. 437 438 Figure captions 439 Figure 1: Global trends in NEP and their contributing factors. Global temporal 440 contributions of CO2, climate and land-use change to the trends in NEP (annual 441 change) are shown on the right side of each panel. The difference between the 442 modelled temporal contributions and the trends (shaded) has been treated as an 443 unknown contribution to the temporal variation in NEP. Statistically significant (P < 444 0.01) temporal variations of the predictors are shown in square brackets. Error bars 445 indicate 95% confidence intervals. The boxplots in panel c indicate the estimated 446 contributions of the 10 DVGMs used in the TRENDY ensemble. Units are ppm y-1 for 447 CO2, °C y-1 for temperature, mm y-2 for precipitation, standard deviation for SPEI, and 448 percentage of land-use cover per pixel for forests, crops, and urban areas. See the 449 Materials and Methods section for information about the methodology used to calculate 450 the contributions. Significance levels: *, P < 0.01; **, P < 0.005; ***, P < 0.001. 451 Figure 2: Plots showing the estimated effects of the interactions of the statistical 452 models. The graphs show interactions between CO2 and climate (mean annual 453 precipitation [MAP] and temperature [MAT], and annual anomalies in temperature 454 [MAT.an]) on NEP for the MACC-II and Jena CarboScope inversion models and the 455 TRENDY ensemble. Shaded bands indicate the 95% confidence intervals of the 456 slopes. Non-significant interactions are indicated by “n.s.”. 457 Table 1: Global and latitudinal analyses of sensitivity of NEP to changes in 458 atmospheric CO2 concentrations and mean annual temperature. The “%” columns 459 indicate the contribution of the latitudinal band to the global estimate. Differences are 460 calculated as the difference between the sum of all latitudinal bands and the global 461 estimate. Bold coefficients differ significantly from 0 at the 0.01 level. Empty cells 462 indicate that anomalies in temperature were not a significant predictor in the models 463 predicting NEP. Units are Tg C y-1 ppm-1 for CO2 and Tg C y-1 C-1 for temperature. 464 Table 2: Global and latitudinal trends and temporal contributions of changes in 465 atmospheric CO2 concentrations and mean annual temperature to NEP trends. 466 The “%” columns indicate the percentage of contribution of each latitudinal band to the 467 global estimate. Columns “Cont.” show the percentage of contribution of CO2 and 468 temperature to the trends in NEP. Column “Other” shows the difference between the 469 NEP trend and the sum of contributions of CO2 and temperature. If different from zero, 470 it indicates that other factors are contributing to the trends in NEP. The “differences” 471 rows are calculated as the difference between the sum of all latitudinal bands and the 472 global estimate. NH and SH indicate Northern and Southern Hemispheres, 473 respectively. Bold coefficients differ significantly from 0 at the 0.01 level. Empty cells 474 indicate that anomalies in temperature were not a significant predictor in the models 475 predicting NEP. Units are Tg C y-1 for trends, Tg C y-1 ppm-1 for CO2 and Tg C y-1 C-1 476 for temperature. Errors were calculated using the error propagation method. See the 477 Materials and Methods section for information about the methods used to calculate the 478 contributions. 479 480 481 482 Figure 1 483 484 485 Figure 2 486 487 488 Table 1 CO2 % Temperature % MACC NH >55° 8.5 ± 0.4 14.1 -35.3 ± 24.1 6.8 NH 35-55° 14.7 ± 1.3 24.3 -132.0 ± 259.9 25.6 NH 15-35° -5.0 ± 1.4 -8.3 NH 15-SH 15° 31.9 ± 0.7 52.9 101.9 ± 216.6 -19.8 SH 15-35° 2.2 ± 0.9 3.7 -150.2 ± 131.3 29.1 SH 35-55° 0.6 ± 0.3 1.0 -13.4 ± 49.3 2.6 Global 60.4 ± 1.2 -515.7 ± 152.4 Difference -7.4 ± 2.6 -12.3 286.6 ± 397.4 -55.6 JENA NH >55° -0.3 ± 1.0 -0.3 -49.8 ± 48.2 5.8 NH 35-55° 11.1 ± 3.9 13.6 -213.6 ± 558.1 24.9 NH 15-35° 26.3 ± 2.7 32.3 -268.7 ± 400.0 31.3 NH 15-SH 15° 54.2 ± 3.6 66.6 -697.6 ± 1136.5 81.2 SH 15-35° 5.4 ± 0.9 6.6 -167.0 ± 133.9 19.4 SH 35-55° 0.2 ± 0.0 0.3 Global 81.4 ± 3.4 -859.2 ± 386.3 Difference 15.4 ± 6.9 19.0 -537.4 ± 1390.2 62.5 TRENDY NH >55° 2.8 ± 0.1 9.0 17.3 ± 7.3 -1.6 NH 35-55° 5.8 ± 0.5 19.0 -251.1 ± 79.3 23.1 NH 15-35° 5.9 ± 0.6 19.4 -368.8 ± 51.9 33.9 NH 15-SH 15° 16.6 ± 1.1 54.2 -1612.2 ± 213.4 148.2 SH 15-35° 4.6 ± 1.2 14.9 -379.2 ± 141.1 34.9 SH 35-55° 0.3 ± 0.2 1.0 -36.8 ± 18.1 3.4 Global 30.7 ± 1.2 -1088.0 ± 118.1 Difference 5.4 ± 2.1 17.5 -1542.7 ± 298.0 141.8 489 490 Table 2 Trends % CO2 % Cont. Temp % Cont. Other MACC NH >55° 20.1 ± 1.2 17.2 17.0 ± 0.8 14.1 84.4 -1.2 ± 0.8 11.5 -5.9 4.3 ± 1.7 NH 35-55° 17.5 ± 5.0 15.0 29.2 ± 2.7 24.3 166.6 -1.7 ± 3.2 16.1 -9.4 -10.0 ± 6.5 NH 15-35° 14.0 ± 3.1 12.0 -9.9 ± 2.8 -8.3 -71.0 0.0 23.9 ± 4.1 NH 15- 55.4 ± 2.7 47.4 63.5 ± 1.5 52.9 114.6 0.9 ± 1.9 -8.9 1.6 -9.0 ± 3.6 SH 15° SH 15-35° 7.6 ± 1.4 6.5 4.4 ± 1.9 3.7 57.6 -2.3 ± 2.0 22.2 -29.8 5.5 ± 3.1 SH 35-55° 2.3 ± 0.6 2.0 1.2 ± 0.7 1.0 49.9 -0.3 ± 1.0 2.5 -11.2 1.4 ± 1.3 Global 116.9 ± 6.1 120.1 ± 2.3 102.7 -10.3 ± 3.0 -8.8 7.1 ± 7.2 Difference 0.0 ± 9.1 0.0 -14.8 ± 5.2 -12.3 5.8 ± 5.4 -56.6 JENA NH >55° 13.8 ± 2.2 7.7 -0.5 ± 2.1 -0.3 -3.8 -1.7 ± 1.7 9.9 -12.4 16.0 ± 3.5 NH 35-55° 49.8 ± 5.9 28.0 22.0 ± 7.7 13.6 44.1 -2.7 ± 6.9 15.4 -5.3 30.5 ± 11.9 NH 15-35° 49.2 ± 4.0 27.6 52.3 ± 5.3 32.3 106.2 -5.0 ± 7.4 29.0 -10.2 1.9 ± 10.0 NH 15- 80.4 ± 5.1 45.2 107.7 ± 7.1 66.6 133.9 -5.7 ± 9.2 32.9 -7.0 -21.6 ± 12.7 SH 15° SH 15-35° 10.4 ± 1.3 5.8 10.7 ± 1.7 6.6 103.1 -2.8 ± 2.2 16.2 -26.9 2.5 ± 3.1 SH 35-55° 0.5 ± 0.1 0.3 0.4 ± 0.1 0.3 87.2 0.1 ± 0.1 Global 178.0 ± 8.1 161.8 ± 6.8 90.9 -17.2 ± 7.7 -9.7 33.4 ± 13.1 Difference 26.1 ± 12.2 14.7 30.7 ± 13.8 19.0 -0.6 ± 16.0 3.4 TRENDY NH >55° 9.3 ± 0.6 41.4 5.5 ± 0.3 9.0 59.0 0.6 ± 0.2 -2.7 6.1 3.3 ± 0.7 NH 35-55° 9.4 ± 1.3 41.5 11.6 ± 0.9 19.0 124.0 -3.0 ± 0.9 13.9 -31.6 0.7 ± 1.8 NH 15-35° 3.3 ± 1.3 14.9 11.8 ± 1.1 19.4 352.9 -7.9 ± 1.0 36.9 -235.0 -0.6 ± 2.0 NH 15- 10.1 ± 2.3 45.0 33.0 ± 2.1 54.2 326.2 -17.2 ± 1.8 80.8 -170.2 -5.7 ± 3.6 SH 15° SH 15-35° -13.7 ± 1.8 -60.9 0.5 ± 0.1 0.9 -3.8 -0.3 ± 0.1 1.6 2.5 -13.9 ± 1.8 SH 35-55° -1.0 ± 0.4 -4.7 0.6 ± 0.5 1.0 -55.4 -0.7 ± 0.4 3.5 70.4 -0.9 ± 0.7 Global 22.5 ± 3.1 61.0 ± 2.5 270.7 -21.3 ± 2.2 -94.7 -17.1 ± 4.5 Difference -5.2 ± 4.7 -22.9 2.1 ± 3.6 3.4 -7.3 ± 3.2 34.0 491 492 493 Methods 494 Datasets 495 NEP data 496 We used gridded global monthly NEP data for 1995–2014 from two inversion models: i) 497 the MACC (Monitoring Atmospheric Composition and Climate) CO2 (http://www.gmes- 25,39 498 atmosphere.eu/catalogue/) database, version v14r2 and ii) the Jena CarboScope 499 database version s93_v3.7 using a constant network of towers (http://www.bgc- 26 500 jena.mpg.de/CarboScope/) . The MACC CO2 atmospheric inversion system relies on 501 the variational formulation of Bayes’ theorem to analyse direct measurements of CO 2 502 concentrations from 130 sites around the globe for 1979-2014. Optimised fluxes were 503 calculated at a global horizontal resolution of 3.75 × 1.875° (longitude, latitude) and a 504 temporal resolution of eight days, separately for daytime and night-time. The underlying 505 transport model was run with interannually varying meteorological data from the 506 ECMWF ERA-Interim reanalysis. The Jena inversion model estimates the interannual 507 variability of CO2 fluxes based on raw CO2 concentration data from 50 sites. The model 508 uses a variational approach with the TM3 transport model (4 × 5°, using interannually 509 varying winds). Prior terrestrial fluxes were obtained from a modelled mean biospheric 510 pattern and fossil-fuel emissions from the EDGAR emission database40. We also used 511 NEP data from an ensemble of 10 dynamic global vegetation models (DGVMs) 512 compiled by the TRENDY project (version 4, models CLM4.5, ISAM, JSBACH, JULES, 513 LPJG, LPX, OCN, ORCHIDEE, VEGAS, and VISIT) to see if results obtained from 514 atmospheric inversions data match those obtained with DGVMs simulations41. We used 515 the output from simulation experiment S3, which was run with varying atmospheric CO 2 516 and changing land use and climate41. 517 Meteorological, land-use change and atmospheric CO2 data 518 We extracted gridded temperature and precipitation time series from the Climatic 42 519 Research Unit TS3.23 dataset . We also used the SPEI (Standardised Precipitation- 520 Evapotranspiration Index) drought index43 from the global SPEI database 521 (http://SPEI.csic.es/database.html) as a measure of drought intensity (positive values 522 indicate wetter than average meteorological conditions, negative values indicate drier 523 than average conditions). We used annual SPEI1 (monthly SPEI averaged over a 524 year). Mean annual temperature (MAT) and precipitation (MAP) and SPEI were 525 calculated for each year and pixel. We used land-use change maps from land-use 526 harmonisation2 (LUH2, http://luh.umd.edu/data.shtml) and calculated the percent 527 coverages of forests, croplands, and urban areas per pixel, so we could further 528 estimate whether they increased or decreased from 1995 to 2014. We used the data 529 for atmospheric CO2 concentration from Mauna Loa Observatory provided by the 530 Scripps Institution of Oceanography (Scripps CO2 programme). 531 Data for N and S deposition 532 Annual data for N (oxidised N [NOX] from NO3- and reduced N [NRED] from NH4+) and S 533 (SO4-) wet deposition were extracted from: i) the European Monitoring and Evaluation 534 Programme (EMEP) with a spatial resolution of 0.15 × 0.15° for longitude and latitude, 535 ii) the MSC-W chemical-transport model developed to estimate regional atmospheric 536 dispersion and deposition of acidifying and eutrophying N and S compounds over 537 Europe, and iii) the National Atmospheric Deposition Program (NADP) covering the 538 USA with a spatial resolution of 0.027 × 0.027° for longitude and latitude. We used only 539 data for wet deposition because the NADP database only contained records for dry 540 deposition for 2000. Analyses focused on atmospheric deposition and were restricted 541 to Europe and the USA because temporal gridded maps of atmospheric deposition 542 were not available for other regions. Maps of atmospheric deposition for the regional 543 analyses were adjusted to the resolution of the C-flux maps (3.75 × 1.875° for the 544 MACC-II model and 4 × 5° for the Jena CarboScope model for longitude and latitude). 545 Statistical analyses 546 Gridded, global and regional trend detection on NEP 547 To determine how NEP has changed from 1995 to 2014, we first calculated the trends 548 for each pixel in both inversion models and an average dataset of the TRENDY 549 ensemble using linear regressions with an autoregressive and moving-average 550 (ARMA) (autoregressive structure at lag p=1, and no moving average q=0) correlation 551 structure to account for temporal autocorrelation. Trends over larger areas (e.g. the 552 entire world, latitudinal bands), either for NEP or the predictor variables, were 553 calculated using generalised linear mixed models (GLMMs) with random slopes, 554 including also random intercepts44 (e.g. NEP ~ year). We used pixel as the random 555 factor (affecting the intercepts and slopes of the year), and an ARMA (p=1, q=0) 556 correlation structure. All average trends shown were calculated using this methodology. 557 Calculation of temporal contributions on trends of NEP 558 The temporal contributions of increasing CO2, climate (MAT, MAP, and SPEI), and 559 land-use change (forests, croplands, and urban areas) to the observed trends in NEP 560 were assessed for the MACC-II, Jena CarboScope, and TRENDY datasets for the 561 entire world. We repeated the analysis for five latitudinal bands to determine if the 562 contributions of CO2, climate, and land-use change were globally consistent using 563 MACC-II, Jena CarboScope, and the mean ensemble of the TRENDY datasets. For the 564 MACC-II and Jena CarboScope datasets, we also determined the temporal contribution 565 of atmospheric deposition of N (NOX and NRED) and S to the trends in NEP in a 566 combined analysis that also included CO2, climatic, and land-use trends. This latter 567 analysis was restricted to Europe and the USA due to the lack of atmospheric- 568 deposition time series for the rest of the world. 569 The temporal contributions of the predictor variables were calculated following the 570 methodology established in references5,45, as follows: 571 i) using a GLMM with an autocorrelation structure for lag 1 (AR1) and using the pixel as 572 the random factor affecting only the intercept, we fitted full models for NEP as a 573 function of CO2, mean MAT per pixel, annual anomaly of MAT, mean MAP per pixel, 574 annual anomaly of MAP, the annual SPEI, and mean percentage of forested, cropped, 575 and urban areas per pixel and their annual anomalies. We included the first-order 576 interaction terms between CO2 and all predictors and between the mean values and 577 the anomalies for all predictors (except SPEI, which interacted with mean MAT and 578 MAP). When the interaction term between the means and the anomalies (e.g. MAT 579 mean × MAT anomaly) was included, the model estimated the effect of the anomaly as 580 a function of the average value. This implies a change in the effect of increasing or 581 decreasing the anomalies, depending on the mean for the site (e.g. increasing 582 temperature may have a positive effect in cold climates but a negative effect in warmer 583 climates). For models including atmospheric deposition, we also included the 584 interaction between climatic variables and CO2 and the interactions between the means 585 and the annual anomalies of atmospheric deposition (NOX, NRED, and S). The models 586 were fitted using maximum likelihood to allow the comparison of models with different 587 fixed factors. 588 ii) We used the stepwise backwards-forwards model selection (stepAIC function in R46) 589 from the full models, using the lowest Bayesian information criterion (BIC), to obtain the 590 best model. The amount of the variance explained by the models was assessed using 47 591 the r.squaredGLMM function in R (MuMIn package: ) following the method of 592 Nakagawa and Schielzeth (2013). Model residuals met the assumptions required in all 593 analyses (normality and homoscedasticity of residuals). 594 iii) We then used the selected models to predict the changes of the response variables 595 during the study period (1995–2014). We first extracted the observed trend (mean ± 596 SEM, standard error of the mean) in NEP using raw data with GLMMs with an AR1 597 autocorrelation structure. We then calculated the trend of NEP predicted by the final 598 model and the trends of NEP predicted by the same model while maintaining the 599 temporally varying predictors (i.e., anomalies) constant one at a time (e.g. MAT 600 anomalies were held constant using the median per pixel, while all other predictors 601 changed based on the observations). The difference between the predictions for the 602 final model and when one predictor was controlled was assumed to be the contribution 603 of that predictor variable to the change in NEP. The differences between all individual 604 contributions and the observed trend in NEP were treated as unknown contributions. 605 Calculation of sensitivities of NEP to temporal predictors 606 Finally, we calculated the average sensitivities of NEP to the predictor changes by 607 dividing the temporal contributions of each predictor of delta NEP by their temporal 608 trends. Spatial variability on the effects of temporal predictors to NEP were assessed 609 using the GLMMs fitted to estimate the temporal contributions of the predictors. To 610 visualise the interactions we used the R package visreg49. All errors were calculated 611 using the error-propagation method using the following two equations, for additions and 612 subtractions: 𝜺𝑪 = √(𝜺𝑨)𝟐 + (𝜺𝑩)𝟐 ; and for multiplications and divisions: 𝜺𝑪 = 𝜺𝑨 𝟐 𝜺𝑩 𝟐 613 𝑪√( 𝑨 ) + ( 𝑩 ) ; where ɛ indicates the error associated to each value (A, B or C). To 614 calculate global and regional estimates we multiplied the model outputs, in units of gC 615 m-2, times land area. We considered the land Earth surface area to be 134375000 km2 616 excluding the Antarctic region. Land area for the different latitudinal bands used were: 617 >55º N, 23818000 km2; 35 to 55º N, 31765000 km2; 15 to 35º N, 29213000 km2; 15º S 618 to 15º N, 29926000 km2; 15 to 35º S, 17308000 km2; and 35 to 55º S, 2345600 km2. 619 620 Data availability 621 The authors declare that the data supporting the findings of this study are publically 622 available in the webpages provided in the article. The TRENDY simulations are 623 available from the corresponding author upon request. 624