This is the accepted version of the journal article: Fernández-Martínez, Marcos; Peñuelas, Josep; Chevallier, Frédéric; [et al.]. «Diagnosing destabilisation risk in global land carbon sinks». Nature, Vol. 615, issue 7954 (March 2023), p. 848–853. DOI 10.1038/s41586-023-05725-1 This version is available at https://ddd.uab.cat/record/285125 under the terms of the license 1 Diagnosing destabilisation risk in global land carbon sinks 2 Authors: Marcos Fernández-Martínez1,2,3, Josep Peñuelas2,4, Frederic Chevallier5, Philippe 3 Ciais5, Michael Obersteiner6,7, Christian Rödenbeck8, Jordi Sardans2,4, Sara Vicca1, Hui Yang5, 4 Stephen Sitch9, Pierre Friedlingstein10, Vivek K. Arora11, Daniel S. Goll5, Atul K. Jain12, Danica 5 L. Lombardozzi13, Patrick C. McGuire14, Ivan A. Janssens1. 6 Affiliations: 1 7 PLECO (Plants and Ecosystems), Department of Biology, University of Antwerp, 2610 Wilrijk, 8 Belgium. 2 9 CREAF, E08193 Bellaterra (Cerdanyola del Vallès), Catalonia, Spain 3 10 BEECA-UB, Department of Evolutionary Biology, Ecology and Environmental Sciences, 11 University of Barcelona, E08028 Barcelona, Catalonia, Spain 4 12 CSIC, Global Ecology Unit, CREAF-CSIC-UAB, Bellaterra, 08193 Barcelona, Catalonia, 13 Spain 14 5 Laboratoire des Sciences du Climat et de l’Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, 15 Université Paris-Saclay, Gif-sur-Yvette, France 6 16 International Institute for Applied Systems Analysis (IIASA), Laxenburg, Austria 7 17 School of Geography and the Environment, University of Oxford, Oxford, UK 8 18 Department of Biogeochmical Systems, Max Planck Institute for Biogeochemistry, Jena, 19 Germany 9 20 College of Life and Environmental Sciences, University of Exeter, Exeter, UK 10 21 College of Engineering, Mathematics, and Physical Sciences, University of Exeter, Exeter, 22 United Kingdom 11 23 Canadian Centre for Climate Modelling and Analysis, Climate Research Division, 24 Environment and Climate Change Canada, Victoria, BC, Canada 12 25 Department of Atmospheric Sciences, University of Illinois, Urbana, IL 61801, USA 13 26 Climate and Global Dynamics Laboratory, National Center for Atmospheric Research, 27 Boulder, CO 80302, United States of America 14 28 Department of Meteorology, Department of Geography & Environmental Science, National 29 Centre for Atmospheric Science, University of Reading, Reading, UK 30 *Correspondence to: M. Fernández-Martínez, m.fernandez@creaf.uab.cat 31 Length of the main text: 3176 words 32 Figures and tables: 3 figures, 0 tables 33 References: 46 (+ 14 in Methods) 34 35 Summary 36 Global net land carbon uptake, or net biome production (NBP), has increased during recent 37 decades1. Whether its temporal variability and autocorrelation have changed during this period, 38 however, remains elusive, even though an increase in both could indicate an increased 39 potential for a destabilised carbon sink2,3. Here we investigate the trends and controls of net 40 terrestrial carbon uptake and its temporal variability and autocorrelation from 1981 to 2018 41 using two atmospheric inversion models, the amplitude of the seasonal cycle of atmospheric 42 CO2 concentration derived from nine monitoring stations distributed across the Pacific Ocean 43 and dynamic global vegetation models. We find that annual NBP and its interdecadal variability 44 increased globally while temporal autocorrelation decreased. We observe a separation of 45 regions characterized by increasingly variable NBP, associated with warm regions and 46 increasingly variable temperatures, lower and weaker trends in NBP, and regions where NBP 47 became stronger and less variable. Plant species richness presented a concave-down 48 parabolic spatial relationship with NBP and its variability at the global scale while nitrogen 49 deposition generally increased NBP. Increasing temperature and its increasing variability 50 appears as the most important drivers of declining and increasingly variable NBP. Our results 51 show increasing variability of NBP regionally that can be mostly attributed to climate change, 52 and that may point to destabilisation of the coupled carbon-climate system. 53 Keywords: carbon cycle, variability, biodiversity, climate, land-use change, early warnings 54 55 Main 56 Positive carbon-climate feedbacks have the potential to accelerate climate change and might 57 compromise the attainability of ambitious climate targets such as those set by the Paris 58 agreement4. Terrestrial ecosystems are key to the functioning of the global carbon (C) cycle 59 and have increased their productivity and net C uptake during recent decades primarily due to 60 CO2 fertilisation and forest regrowth1,5–7. However, land-use change, nutrient limitations and 61 increasing droughts and fires are constraining this potential to sequester C8–11. Identifying 62 processes that might destabilise net land C uptake (or net biome production, NBP) is of 63 paramount importance for understanding and managing the global C cycle. 64 Destabilisation is the process of losing stability: that is, when a system loses its ability to return 65 to equilibrium following a disturbance. A symptom of this may be increased variability, as the 66 system spirals away from its current equilibrium point towards a new one. Destabilisation of a 67 dynamical system is usually accompanied by concomitant increases in temporal variability and 68 autocorrelation (i.e., the correlation between consecutive time steps [AR1] also related to 69 reduced resilience2) when AR1 is already positive, because anomalous states of the system 70 create ripples that can get amplified through time rather than compensated as when AR1 is 71 negative. Consequently, increasing temporal variability and autocorrelation have been shown 72 to be potential early-warning signals for abrupt shifts in ecosystems2,3,12,13. To date, changes 73 in NBP temporal variability and autocorrelation have not yet been investigated even though 74 increased stress and changes in the frequency and intensity of extreme weather events are 75 likely to alter the temporal patterns of NBP14–16 due to i) a cumulative negative effect of extreme 76 events on ecosystem functioning, ii) increasing climate variability and, iii) decreasing 77 ecosystem resilience due to increased stress. We, thus, hypothesised that regions 78 experiencing increasing trends in climate variability will also experience increasing variability 79 in their NBP, and that increase in variability may indicate that a destabilisation of the net land 80 C uptake is occurring. This finding could serve as an early warning signal for abrupt shifts in 81 ecosystem’s functions2 that might lead to regime shifts in the Earth’s biosphere17,18. Even if 82 these changes do not occur at the global scale, increasing temporal variability and 83 autocorrelation in several regions of the globe (e.g. Amazon basin, boreal ecosystems with 84 permafrost) could have a profound impact on the global C balance and a knock-on effect on 85 other ecosystem functions. 86 Climate is the primary control on NBP in terrestrial ecosystems worldwide in space and 87 time19,20, together with soil nutrient availability10,21, atmospheric nitrogen (N) deposition, land 88 use and management, and increasing atmospheric CO21. Generally, gross C fluxes (such as 89 photosynthesis and respiration) are larger in the tropics where high temperature coincides with 90 sufficient precipitation to enable long growing seasons. However, net C uptake tends to be 91 higher in temperate regions due to higher nutrient availability19,21–23. Nutrient availability and N 92 deposition have indeed been shown to increase net land C uptake21,22. Another factor believed 93 to be important in determining ecosystem functioning is plant biodiversity24. A large body of 94 evidence indicates its role in promoting ecosystem productivity and stability25,26. Plant 95 biodiversity, however, has been included far less in studies of ecosystem C cycling, 96 presumably due to the difficulty of acquiring good data on species diversity. The few studies 97 that included biodiversity, however, showed relatively modest correlations with C fluxes27–29. 98 The role of biodiversity in the global terrestrial NBP has not yet been explored even though the 99 Earth’s biosphere is losing biodiversity at an unprecedented speed30 and those changes are 100 expected to alter ecosystem functioning. 101 The aim of this study was to quantify the trends in global NBP and its intradecadal temporal 102 variability (quantified by the proportional variability index, PV31,32) and autocorrelation (AR1) to 103 see whether changes in these variables suggest that NBP at global or regional scales is being 104 destabilised. We further tested whether regions with increasing variability and AR1 in NBP 105 showed differential trends in annual NBP. We additionally investigated the spatial correlation 106 between NBP metrics and plant biodiversity (derived from a global map of plant species 107 richness33), atmospheric total N deposition34, climate35, and land use (land-use harmonisation2 108 maps). To realise these objectives, we used estimates of NBP derived from the two longest 109 CO2 atmospheric inversions (CAMS, CarboScope) and the amplitude of the seasonal cycle of 110 atmospheric CO2 derived from nine monitoring stations distributed from south to north of the 111 Pacific Ocean for 1981 to 2018. We additionally compared these results with the output from 112 an ensemble of 12 dynamic global vegetation models (TRENDY) to explore how well these 113 state-of-the-art models predict the spatial and temporal patterns in NBP simulated by 114 atmospheric inversion models. 115 Results 116 Temporal patterns of net land C uptake 117 Global NBP derived from atmospheric inversions increased from 5.6 ± 2.0 gC m-2 y-1 during 118 1981 – 1990 (mean ± standard error) to 13.8 ± 1.4 gC m-2 y-1 during 2009 – 2018 over the 119 global land area excluding Antarctica (Figure 1). This represents an overall increase of 145% 120 and an annual linear increase of 0.24 ± 0.08 gC m-2 y-2 (P<0.001). Both atmospheric inversions 121 used (CAMS and CarboScope) identified annual increases in NBP (Supplementary Figure 1 122 and Supplementary Figure 2, 0.32 ± 0.09, and 0.18 ± 0.09 gC m-2 y-2 respectively, P<0.001, 123 N=38). The positive trend in NBP shown by inversions is similar in magnitude to the trend 124 identified by the TRENDY ensemble (0.10 ± 0.07 gC m-2 y-2, P<0.001; 147% when comparing 125 the periods 1981 – 1990 and 2009 – 2018). Both observations and models, however, showed 126 a flattening trend during the last decade (Figure 1c). In parallel with the increased NBP 127 throughout the entire study period and in agreement with previous literature36,37, the amplitude 128 of the seasonal atmospheric CO2 concentration increased by 0.027 ± 0.004 ppm y-1 (P<0.001, 129 16.6% during the study period) (Supplementary Figure 3a). 130 Interannual NBP variability (NBPPV) derived from the combination of atmospheric inversions 131 increased globally by 7.2% over the entire period (P<0.001) (Figure 1f) (CAMS 14.9%, 132 P<0.001; CarboScope: -1.3%, P<0.001). The origin of the discrepancy between the inversions 133 remains unclear. Though the inversions differ in whether they use a growing number of 134 measurement stations (CAMS) or a constant station set (CarboScope), the discrepancy is 135 unlikely to be related to this because a test inversion with a growing station set in the 136 CarboScope inversion did not yield an increase in NBPPV. Temporal variability in the amplitude 137 of the seasonal atmospheric CO2 concentration also increased (by 35%, P<0.001) during the 138 study period (Supplementary Figure 3b). TRENDY NBP did not simulate this increase in 139 variability, showing a negative trend of -1.7% in NBPPV similar to that from the CarboScope 140 inversion. 141 Global NBPAR1 significantly decreased over time for NBP derived from atmospheric inversions 142 (Figure 1i) and for the average monthly CO2 concentrations across the Pacific Ocean 143 (Supplementary Figure 3c). TRENDY NBPAR1, also showed a significant negative trend 144 representing a reduction of 2.9% over the study period. Trends of contrasting sign, however, 145 were significant amongst atmospheric measurement stations and atmospheric inversions 146 (CAMS: -10.0%, P<0.001; CarboScope: 4.7%, P<0.001). Several regions showed increases 147 in the temporal autocorrelation between consecutive months (NBPAR1) derived from 148 atmospheric inversions that are of a similar magnitude (ca. 0.2 over three decades) to those 149 previously suggested to precede abrupt shifts in climate datasets and simulations2,3. Additional 150 information on how changes in NBPPV and NBPAR1 may affect annual NBP and their limitations 151 can be found in Supplementary text, Section 1. 152 Our analyses identified several regions of potential concern given their concomitant increase 153 in NBPPV and NBPAR1 (Figure 2a), such as eastern Africa, the Mediterranean region, the west 154 coasts of North and Central America, India, and southeast Asia. Regions with increasing 155 NBPPV and NBPAR1 had statistically lower NBP and experienced a much less pronounced 156 increase in NBP over time (0.15 ± 0.06 gC m-2 y-1, P<0.001) compared to regions where NBPPV 157 and NBPAR1 decreased (0.73 ± 0.12 gC m-2 y-1, P<0.001) (Figure 2b-d, Supplementary 158 Figure 4). Stronger increases in NBP over regions with decreasing NBPPV and NBPAR1, 159 compared to those where they increased, were also evident from both atmospheric inversions 160 when analysed separately (Supplementary Figure 5). 161 Global temperature and precipitation increased during the study period (Supplementary 162 Figure 6a and d) even though there was considerable spatial variability in those trends 163 (Supplementary Figure 7a and d). Temperature and precipitation interannual variability, 164 however, decreased significantly at the global scale, despite the reported increase in extreme 165 weather14. Temporal autocorrelation of monthly temperature and precipitation, instead, 166 increased slightly, albeit only significantly for temperature. A concomitant increase in temporal 167 variability and autocorrelation of temperature was evident in several regions (e.g., Eurasia, 168 Australia, Central America) but those were not so obvious for precipitation (Supplementary 169 Figure 8). 170 Controls of NBP, variability and autocorrelation 171 NBP, derived from atmospheric inversions, had a concave-down parabolic relationship with 172 biodiversity (Figure 3), increasing from low to intermediate values of biodiversity and 173 decreasing at high biodiversity (Supplementary Figure 9a). Our analyses also identified a 174 significant positive interaction between biodiversity and N deposition accounting for NBP. The 175 concave-down relationship between biodiversity and NBP included mainly positive NBP values 176 across regions with high atmospheric N deposition (850 mg N m-2 y-1) and mainly negative 177 NBP values at low atmospheric N deposition (100 mg N m-2 y-1) (Supplementary Figure 9a). 178 Similarly, the positive correlation between N deposition and NBP was stronger in regions with 179 higher plant biodiversity (Supplementary Figure 9b). Opposite relationships emerged when 180 analysing TRENDY NBP. 181 NBPPV derived from inversions also showed a concave-down parabolic relationship with plant 182 biodiversity, peaking at intermediate to high values of biodiversity (Supplementary Figure 183 9d). We also found an interaction between plant biodiversity and N deposition in their 184 relationship with NBPPV: the relationship between NBPPV and biodiversity differed between 185 areas receiving low N deposition and areas receiving high N deposition (Supplementary 186 Figure 9d). No relationship between N deposition and NBPPV was found in regions with high 187 biodiversity, whereas a positive relationship occurred in regions with low biodiversity 188 (Supplementary Figure 9e). Interannual variability in temperature and precipitation were 189 positively correlated with NBPPV (Figure 3). In this case, the patterns for the TRENDY 190 ensemble and the inversion models matched very well. Again, NBPAR1 derived from inversions 191 had a concave-down parabolic relationship with biodiversity (Supplementary Figure 9g), but 192 here no interaction occurred between the effects of biodiversity and N deposition, which was 193 negatively correlated with NBPAR1. NBPAR1 was also positively correlated with temperature 194 AR1, the only result that match those from TRENDY NBPAR1 (Supplementary Figure 9i). 195 Results emerging from individual atmospheric inversions (CAMS and CarboScope) mostly 196 coincided with those reported above (Supplementary Figure 10). 197 Spatial variability in the trends of NBP, NBPPV and NBPAR1 derived from atmospheric inversions 198 were all correlated with N deposition, climate, and land use (Figure 3 and Supplementary 199 Figure 11). The estimated effects of land use and land use change, however, were generally 200 lower than those from N deposition and climate, hence accounting for a smaller proportion of 201 the change at the global scale. N deposition was positively correlated with the trends in NBP 202 and negatively with trends in NBPPV and NBPAR1 (Supplementary Figure 11). Warmer regions 203 were more likely than colder regions to have decreasing trends in NBP, and increasing 204 temperatures contributed to decreasing NBP (Supplementary Figure 11b and c). Regions 205 with the strongest increases in NBPPV were spatially associated with low increases in annual 206 temperature (Supplementary Figure 11e). Increases in NBPPV were also found to be more 207 likely in regions showing increases in temperature temporal variability (Supplementary Figure 208 11f). Increasing NBPAR1 was more likely in regions with increasing temperatures 209 (Supplementary Figure 11h) and mildly with increasing temporal autocorrelation of 210 precipitation (Figure 3). Our analyses using TRENDY reproduced the abovementioned 211 findings well for trends in NBP and NBPAR1, but not for NBPPV. Again, results from individual 212 atmospheric inversions were very similar to those reported here (Supplementary Figure 10). 213 Warm regions presenting a concomitant increase in temperature temporal variability and 214 autocorrelation (symptoms of destabilisation) were more likely to be correlated with similar 215 changes in NBP variability and AR1, and those results were well supported by both 216 atmospheric inversions and their combination (Supplementary Figure 12). Additionally, 217 regions with higher N deposition, larger proportions of forested areas and lower crops were 218 related to concomitantly decreasing NBPPV and NBPAR1. 219 Biodiversity, N deposition and net C uptake 220 NBP derived from atmospheric inversions provided correlational evidence indicating that plant 221 biodiversity may be playing an important role in regulating regional variation in the land C 222 balance and in its temporal variability. TRENDY models, instead, do not include biodiversity in 223 their parameterisation and, hence, any spurious relationship is necessarily driven by factors 224 other than biodiversity. The spatial relationship between biodiversity and NBP clearly differed 225 between atmospheric inversions and TRENDY (Supplementary Figure 9a), hence 226 suggesting that the reported effect of biodiversity when analysing atmospheric inversions may 227 emerge due to a mechanistic effect. The emerging relationship between biodiversity and 228 NBPPV derived from atmospheric inversions, however, was very similar to the one emerging 229 from TRENDY, which suggests that factors other than biodiversity may be driving this 230 relationship. The positive relationship between biodiversity and NBPAR1 found here has never 231 been reported before and further research is needed to understand the mechanisms behind 232 this relationship. 233 The concave-down parabolic relationships of biodiversity with NBP and NBPPV differ from the 234 majority of biodiversity-productivity and stability relationships reported in the literature: positive 235 asymptotic for productivity and negative for variability38,39. This difference in the biodiversity- 236 NBP relationship may result from two opposing ecosystem processes, photosynthesis and 237 respiration, because both are likely to be enhanced by biodiversity27,40,41. The concave-down 238 relationship would then suggest that the positive effect of biodiversity on respiration 239 overshadows the positive effect of photosynthesis in regions with high biodiversity. However, 240 our biodiversity data, like in similar studies39, was restricted to species richness and did not 241 include information on species abundance, their individual contribution to NBP, or traits to allow 242 the calculation of actual species diversity or functional diversity, often better indicators of 243 ecosystem functioning than species richness42. Unfortunately, this information is not available 244 at the global scale. Including actual diversity in future analyses could lead to different results 245 to those reported here. Additionally, future efforts are needed to understand how biodiversity 246 loss, not included here, will impact global carbon balance. 247 Interestingly, our analyses indicated that the effect of biodiversity on NBP depended on 248 atmospheric N deposition, and vice versa (Figure 3). The effect of atmospheric N deposition 249 on NBP was mainly positive but was stronger in regions with higher biodiversity, further 250 supporting the premise that biodiversity promotes ecosystem functions such as N uptake43. 251 Regions with higher N deposition also had larger increases in NBP over time (Supplementary 252 Figure 11a), supporting previous findings suggesting a stronger CO2 fertilisation effect in 253 regions with higher N deposition1. On the other hand, N deposition was negatively related to 254 trends in NBPPV, NBPAR1 and the aggregated trend of NBPPV-AR1 (Supplementary Figure 12), 255 which suggests that N addition may ameliorate nutrient imbalances derived from increasing 256 atmospheric CO2 concentrations8 and prevent ecosystem functioning from becoming more 257 variable. 258 Climate and changing land C uptake 259 Our results clearly indicate that NBP decreased in regions with warm climates and where 260 warming has been most pronounced (Supplementary Figure 11). These findings support 261 results from previous studies1 and further suggests that increasing droughts and heat waves 262 limit C sequestration by terrestrial ecosystems44. Our analyses revealed that higher NBPPV, 263 NBPAR1 and their trends were associated, respectively, with climates showing higher and 264 increasing temporal variability and autocorrelation (Supplementary Figure 9). We also found 265 that concomitantly increasing trends in NBPPV-AR1 were positively related to increasing 266 temperature variability and autocorrelation (Supplementary Figure 12). These results support 267 our initial hypothesis stating that climate change may be the main contributing factor for 268 changing the temporal behaviour of land C sinks. 269 Implications of altered land C sinks 270 Compelling evidence from atmospheric inversions (Figure 1) and the annual amplitude of the 271 seasonal CO2 cycle (Supplementary Figure 3) suggests that global net land C uptake is now 272 larger and, most likely, more variable than three decades ago, while its temporal 273 autocorrelation (NBPAR1), instead, has significantly decreased. The increase in global NBP and 274 the reduction in NBPAR1 were well identified by TRENDY NBP, including the evidence of a 275 flattening trend in NBP during the last two decades. Our results for the increase in global NBP 276 and a potential recent saturation are consistent with previous findings of increases in global 277 productivity and NBP and a recent decline of the CO2 fertilisation effect1,5,8,20. Even though the 278 main emerging spatial controls for CAMS and CarboScope atmospheric inversions in NBP 279 were in agreement (Figures S8 and S10), trends at the global scale for NBPPV and NBPAR1 280 were significantly different (Supplementary Figure 13), which calls for caution in the 281 interpretation of these results due to their uncertainty. 282 We observed a bifurcation of regions characterised by concomitant increases and decreases 283 in NBPPV and NBPAR1 (Figure 2). This bifurcation relates to differences in their mean annual 284 temperature, their increase in annual temperature, their average N deposition loads and their 285 percentage of forests and crops (Supplementary Figure 12). The observed increase in the 286 variability and autocorrelation of net C uptake in several regions of the planet (Figure 2a) is 287 concerning because of the implications it can have on the stability of their ecosystems. First, 288 increasing NBPPV and NBPAR1 is indicative of increasing variability and reducing resilience3,45 289 potentially in many ecosystem processes: from photosynthesis and respiration27 to cascading 290 effects on animals and decomposers46. Second, an increase in variability in NBP implies that 291 ecosystem C balance is less predictable over time, which is troublesome for projections of 292 future climate change. Third, regions showing a combined increase in NBPPV and NBPAR1 had 293 a consistently lower NBP and lower increases in NBP than regions where NBPPV and NBPAR1 294 decreased (Figure 2b-d), suggesting that these increases NBPPV and NBPAR1 reflect 295 dynamical instability in the carbon system rather than some other cause (such as changes in 296 the external forcing of the system). This destabilization may jeopardise future C sequestration. 297 Our analyses, however, could not determine the mechanisms driving the observed changes in 298 NBPPV and NBPAR1. More importantly, we could not determine whether the observed increases 299 NBPPV and NBPAR1 truly reflect dynamical instability in the carbon system, as opposed to some 300 other cause, such as changes in the external forcing of the system (e.g., increasing CO2 301 emissions, N deposition). Hence, regions showing increased variability and autocorrelation 302 should be monitored in detail to properly understand the mechanisms and consequences 303 behind these changes given that increasing variability and autocorrelation have been shown 304 to act as early warning signals preceding abrupt phase transitions in simulations of ecosystem 305 functioning2,3. The increase of NBPPV and NBPAR1 in several regions should serve as an early 306 warning signal of potential future changes that the Earth’s biosphere may be facing. Given the 307 main role of climate change as a driver of these changes in their temporal behaviour, mitigating 308 climate change is needed to prevent further unforeseen changes in land C sinks. 309 References 310 1. Fernández-Martínez, M. et al. Global trends in carbon sinks and their relationships 311 with CO2 and temperature. Nat. Clim. Chang. 9, 73–79 (2019). 312 2. Scheffer, M. et al. Early-warning signals for critical transitions. Nature 461, 53–59 313 (2009). 314 3. Dakos, V. et al. Slowing down as an early warning signal for abrupt climate change. 315 Proc. Natl. Acad. Sci. U. S. A. 105, 14308–14312 (2008). 316 4. Gasser, T. et al. Path-dependent reductions in CO2 emission budgets caused by 317 permafrost carbon release. Nat. Geosci. 11, 830–835 (2018). 318 5. Zhu, Z. et al. Greening of the Earth and its drivers. Nat. Clim. Chang. 6, 791–795 319 (2016). 320 6. Bastos, A. et al. Contrasting effects of CO2 fertilization, land-use change and warming 321 on seasonal amplitude of Northern Hemisphere CO2 exchange. Atmos. Chem. Phys. 322 19, 12361–12375 (2019). 323 7. Pugh, T. A. M. et al. Role of forest regrowth in global carbon sink dynamics. Proc. 324 Natl. Acad. Sci. U. S. A. 116, 4382–4387 (2019). 325 8. Wang, S. et al. Recent global decline of CO2 fertilization effects on vegetation 326 photosynthesis. Science (80-. ). 370, 1295–1300 (2020). 327 9. Peñuelas, J. et al. Assessment of the impacts of climate change on Mediterranean 328 terrestrial ecosystems based on data from field experiments and long-term monitored 329 field gradients in Catalonia. Environ. Exp. Bot. (2017) 330 doi:10.1016/j.envexpbot.2017.05.012. 331 10. Terrer, C. et al. Nitrogen and phosphorus constrain the CO2 fertilization of global plant 332 biomass. Nat. Clim. Chang. (2019) doi:10.1038/s41558-019-0545-2. 333 11. Gatti, L. V. et al. Amazonia as a carbon source linked to deforestation and climate 334 change. Nature 595, 388–393 (2021). 335 12. Carpenter, S. R. & Brock, W. A. Rising variance: A leading indicator of ecological 336 transition. Ecol. Lett. 9, 311–318 (2006). 337 13. Dakos, V., Nes, E. H. & Scheffer, M. Flickering as an early warning signal. Theor. 338 Ecol. 6, 309–317 (2013). 339 14. Sillmann, J., Daloz, A. S., Schaller, N. & Schwingshackl, C. Extreme weather and 340 climate change. Climate Change (Elsevier B.V., 2021). doi:10.1016/b978-0-12- 341 821575-3.00016-5. 342 15. Reichstein, M. et al. Climate extremes and the carbon cycle. Nature 500, 287–95 343 (2013). 344 16. Wang, X. et al. A two-fold increase of carbon cycle sensitivity to tropical temperature 345 variations. Nature 506, 212–5 (2014). 346 17. Barnosky, A. D. et al. Approaching a state shift in Earth’s biosphere. Nature 486, 52–8 347 (2012). 348 18. Buermann, W. et al. Climate-driven shifts in continental net primary production 349 implicated as a driver of a recent abrupt increase in the land carbon sink. 350 Biogeosciences 13, 1597–1607 (2016). 351 19. Luyssaert, S. et al. CO 2 balance of boreal, temperate, and tropical forests derived 352 from a global database. Glob. Chang. Biol. 13, 2509–2537 (2007). 353 20. Peñuelas, J. et al. Shifting from a fertilization-dominated to a warming dominated 354 period. Nat. Ecol. Evol. 1, 1438–1445 (2017). 355 21. Fernández-Martínez, M. et al. Nutrient availability as the key regulator of global forest 356 carbon balance. Nat. Clim. Chang. 4, 471–476 (2014). 357 22. Fernández-Martínez, M. et al. Spatial variability and controls over biomass stocks, 358 carbon fluxes and resource-use efficiencies in forest ecosystems. Trees, Struct. Funct. 359 28, 597–611 (2014). 360 23. Ciais, P. et al. Five decades of northern land carbon uptake revealed by the 361 interhemispheric CO2 gradient. Nature 568, 221–225 (2019). 362 24. Tilman, D., Lehman, C. L. & Thomson, K. T. Plant diversity and ecosystem 363 productivity: Theoretical considerations. Proc. Natl. Acad. Sci. 94, 1857–1861 (1997). 364 25. de Mazancourt, C. et al. Predicting ecosystem stability from community composition 365 and biodiversity. Ecol. Lett. 16, 617–625 (2013). 366 26. Sakschewski, B. et al. Resilience of Amazon forests emerges from plant trait diversity. 367 Nat. Clim. Chang. 6, 1032–1036 (2016). 368 27. Fernández‐Martínez, M. et al. The role of climate, foliar stoichiometry and plant 369 diversity on ecosystem carbon balance. Glob. Chang. Biol. 26, 7067–7078 (2020). 370 28. Musavi, T. et al. Stand age and species richness dampen interannual variation of 371 ecosystem-level photosynthetic capacity. Nat. Ecol. Evol. 1, 0048 (2017). 372 29. Anderegg, W. R. L. et al. Hydraulic diversity of forests regulates ecosystem resilience 373 during drought. Nature 561, 538–541 (2018). 374 30. Díaz, S. et al. IPBES, 2019. Summary for policymakers of the global assessment 375 report on biodiversity and ecosystem services of the Intergovernmental Science‐Policy 376 Platform on Biodiversity and Ecosystem Services. (2019) doi:10.1111/padr.12283. 377 31. Heath, J. P. Quantifying temporal variability in population abundances. Oikos 115, 378 573–581 (2006). 379 32. Fernández-Martínez, M., Vicca, S., Janssens, I. A., Martín-Vide, J. & Peñuelas, J. The 380 consecutive disparity index, D, as measure of temporal variability in ecological studies. 381 Ecosphere 9, e02527 (2018). 382 33. Kreft, H. & Jetz, W. Global patterns and determinants of vascular plant diversity. Proc 383 Natl Acad Sci U S A 104, 5925–5930 (2007). 384 34. Ackerman, D. E., Chen, X. & Millet, D. B. Global nitrogen deposition (2°×2.5° grid 385 resolution) simulated with GEOS-Chem for 1984-1986, 1994-1996, 2004-2006, and 386 2014-2016. (2018) doi:10.13020/D6KX2R. 387 35. Harris, I., Jones, P. D. D., Osborn, T. J. J. & Lister, D. H. H. Updated high-resolution 388 grids of monthly climatic observations - the CRU TS3.10 Dataset. Int. J. Climatol. 34, 389 623 – 642 (2013). 390 36. Graven, H. D. et al. Enhanced seasonal exchange of CO2 by northern ecosystems 391 since 1960. Science 341, 1085–9 (2013). 392 37. Wang, K. et al. Causes of slowing-down seasonal CO2 amplitude at Mauna Loa. Glob. 393 Chang. Biol. 26, 4462–4477 (2020). 394 38. Tilman, D., Reich, P. B. & Knops, J. M. H. Biodiversity and ecosystem stability in a 395 decade-long grassland experiment. Nature 441, 629–32 (2006). 396 39. Liang, J. et al. Positive biodiversity-productivity relationship predominant in global 397 forests. Science (80-. ). 354, aaf8957–aaf8957 (2016). 398 40. Gessner, M. O. et al. Diversity meets decomposition. Trends Ecol. Evol. 25, 372–380 399 (2010). 400 41. Peguero, G. et al. Fast attrition of springtail communities by experimental drought and 401 richness–decomposition relationships across Europe. Glob. Chang. Biol. 25, 2727– 402 2738 (2019). 403 42. Díaz, S. & Cabido, M. Vive la différence: Plant functional diversity matters to 404 ecosystem processes. Trends Ecol. Evol. 16, 646–655 (2001). 405 43. Cardinale, B. J. Biodiversity improves water quality through niche partitioning. Nature 406 472, 86–91 (2011). 407 44. Ciais, P. et al. Europe-wide reduction in primary productivity caused by the heat and 408 drought in 2003. Nature 437, 529–533 (2005). 409 45. Scheffer, M. Critical Transitions in Nature and Society. (Princeton University Press, 410 2009). 411 46. Ostfeld, R. & Keesing, F. Pulsed resources and community dynamics of consumers in 412 terrestrial ecosystems. Trends Ecol. Evol. 15, 232–237 (2000). 413 414 Figure 1: Global distribution of net biome production (NBP), its interannual variability 415 (NBPPV) and temporal autocorrelation (NBPAR1), and their trends from 1981 to 2018. Panel 416 a shows the average global distribution of NBP (positive values indicate C sinks) derived from 417 CAMS and CarboScope atmospheric inversions. Panel b shows the temporal trends in NBP 418 (∆NBP), and panel c shows the temporal change in NBP derived from atmospheric inversions 419 and 12 DGVMs (TRENDY). Panels d, e, and f, and g, h, and i mimic panels a, b, and c, 420 respectively, for interannual variability of NBP (NBPPV) and temporal autocorrelation (NBPAR1) 421 calculated over an 11-year moving window. Thick lines in panels c, f, and i represent smooth 422 trends estimated with a local regression (shaded areas show the standard error of the trend), 423 and the coefficients at the top of the panels indicate robust Theil-Sen’s slopes (± standard error 424 of the slope) with their corresponding p-values, estimated using the non-parametric Wilcoxon 425 Mann-Whitney rank sum test. Values in panels c, f and i, represent spatially-weighted global 426 means. **, P<0.01; ***, P<0.001. Data to reproduce this figure can be found in Source Data 427 1. 428 (See Figure 1 on the following page) 429 430 Figure 1 431 432 Figure 2: Regions with concomitantly increasing NBPPV and NBPAR1 present lower NBP and a lower increase of NBP over time. Panel a 433 shows regions of the world where concomitant increases (red) and decreases (blue) in NBPPV and NBPAR1 have occurred during the study period. 434 Regions with contrasting trends in NBPPV and NBPAR1 appear in grey. Panel b shows the aggregated annual mean NBP of regions with increasing 435 (red) and decreasing (blue) NBPPV and NBPAR1 (see Supplementary Figure 4 for information on their global contribution). Smooth trends were 436 estimated with a local regression (shaded indicate the standard error of the trend), and the coefficients at the top of the panel indicate robust 437 Theil-Sen’s slopes (± standard error of the slope) and their associated significance (***, P<0.001), estimated using the non-parametric Wilcoxon 438 Mann-Whitney rank sum test. Panel c and d show the spatial relationship of NBP and trends in NBP (ΔNBP) with the aggregated trend of NBPPV 439 and NBPAR1. Coefficients are standardised and include the 95% credible intervals (CI). See Methods for details on how the aggregated trend of 440 NBPPV and NBPAR1 was calculated. Data to reproduce this figure can be found in Source Data 2. 441 442 Figure 3: Contribution of biodiversity, N deposition, climate, and land use to the 443 global terrestrial C balance (NBP), its interannual variability (NBPPV), its temporal 444 autocorrelation (NBPAR1), and their trends. The colour scale indicates the strength of the 445 relationship (standardised β coefficients) between each predictor (bottom) and the response 446 variable (left). Black dots indicate that 95% of the posterior distributions differed from 0. 447 Biodiversity was fitted as a second-order polynomial to account for nonlinearities 448 (biodiversity + biodiversity2; the second term indicates the change in the slope in biodiversity 449 as biodiversity increases). Hashed areas indicate relationships not included in the 450 regression models. See Methods for further information on model fitting. INV and TRD 451 indicate the averages for all atmospheric inversions (2) and DGVMs (12), respectively. 452 Predictors MAT𝑥̅ , PV, AR1 and MAP𝑥̅ , PV, AR1 indicate mean, PV, or AR1 metrics of MAT and 453 MAP matching the metric of the response variable (e.g. NBPINV ~ MAT𝑥̅ ,; NBPPV-INV ~ 454 MATPV). The same applies for ΔMATPV, ΔMATAR1, ΔMAPPV, ΔMAPAR1. MAT, mean annual 455 temperature; MAP, mean annual precipitation; PV, proportional variability; AR1, temporal 456 autocorrelation with previous month; Δ, trend in a given variable. 457 458 459 Methods 460 Data sets 461 NBP and atmospheric CO2 data 462 We obtained the global NBP data for 1981–2018 from the two atmospheric-inversion 463 models that provided the longest time series: i) the CAMS Greenhouse Gases Flux 464 Inversions (https://ads.atmosphere.copernicus.eu/cdsapp#!/dataset/cams-global- 465 greenhouse-gas-inversion)47,48 version v18r3, and ii) the Jena CarboScope database 466 version s81oc_v2020 using a constant network of measurement stations (http://www.bgc- 467 jena.mpg.de/CarboScope/)49,50. We also used NBP data from an ensemble of 12 dynamic 468 global vegetation models (DGVMs) run with varying concentrations of atmospheric CO2 and 469 changing land uses and climates. Models compiled by the TRENDY project (version 8, 470 models CABLE-POP, CLASS-CTEM, CLM5.0, ISAM, JSBACH, JULES-ES, LPX-BERN, 471 OCN, ORCHIDEE, ORCHIDEE-CNP, SDGVM and VISIT) were used to test whether the 472 DGVMs also identified the patterns from atmospheric inversions51. Even though previous 473 studies indicate that DGVMs explain trends in NBP adequately, we are the first to test 474 whether their emerging trends in NBP temporal variability and autocorrelation also match 475 those from local observations or atmospheric inversions. We used model results from the 476 simulation experiment S3, which was run with changing atmospheric CO2, land use and 477 climate (see https://blogs.exeter.ac.uk/trendy/protocol/ for more details). We used monthly 478 NBP estimates as the basis for all calculations in this study. We rescaled all atmospheric 479 inversions, the TRENDY model outputs, and the predictors to the same spatial resolution 480 of the coarsest data set for fitting the statistical models (see section below) (i.e. CAMS: 3.75 481 × 1.875°). 482 We calculated average annual NBP per pixel, its temporal variability expressed as the 483 proportional variability index (NBPPV)31,32, and its monthly temporal autocorrelation at lag 1 484 (NBPAR1) for all of the abovementioned datasets. All three indices (NBP, NBPPV and NBPAR1) 485 were calculated as the average for the entire period. The PV index is calculated as the 486 mean proportional variability amongst all possible combinations of values in a time series, 2 ∑𝑚 487 following: 𝑃𝑉 = 𝑛(𝑛−1) , where “n” indicates the length of the variable, where “m” is min(𝑧𝑖 ,𝑧𝑗 ) 488 calculated as: 𝑚 = 1 − max(𝑧𝑖 ,𝑧𝑗 ) , and where “z” represent the individual values from which 489 to calculate all the pairwise comparisons between the observations of the time series (e.g., 490 observation zi vs observation zj). When negative values occurred in a time series, we added 491 a constant to the entire time series equivalent to the minimum absolute value plus one. 492 Unlike other metrics of temporal variability, the PV index provides estimates of temporal 493 variability that are independent of the mean of the time series and that have been proven 494 to be robust even when comparing non-normally distributed data sets, thereby overcoming 495 the mathematical drawbacks of similar variability indices such as the standard deviation or 496 the coefficient of variation52. NBPPV per pixel was estimated as the average of the 𝑁𝐵𝑃𝑃𝑉 𝐽𝑎𝑛𝑢𝑎𝑟𝑦 + [...]+ 𝑁𝐵𝑃𝑃𝑉 𝐷𝑒𝑐𝑒𝑚𝑏𝑒𝑟 497 interannual variabilities of all months (e.g. 𝑁𝐵𝑃𝑃𝑉 = 12 ). By 498 using this method, NBPPV reflects not only the interannual variability of the annual NBP, but 499 also the variation occurring amongst equal months across the years. Using this method, a 500 year with average NBP yet with anomalously low NBP in one period offset by anomalously 501 high NBP in another, would thus result in high interannual variability rather than low 502 interannual variability. AR1 was estimated using the residuals of generalised additive 503 models (GAM) used to remove the trend and seasonal cycle of the data. In these models 504 response variable was monthly NBP and the predictors were the month of the year (i.e. a 505 factor of 12 levels) and the year, the latter included as a spline smoothed term to account 506 for nonlinear trends over time (Supplementary Figure 14). We used the mgcv R package 507 to fit the GAM models53. An 11-year moving average from NBPPV and NBPAR1 per pixel was 508 then calculated for all data sets to investigate trends in interannual variability and temporal 509 autocorrelation. Our sensitivity analysis that showed that our results were consistent despite 510 the selection of different window lengths (7, 11 and 15 years, Supplementary Figure 13, 511 Supplementary Table 1). Robust Theil-Sen’s trends54 were then calculated for NBP, 512 NBPPV and NBPAR1 (ΔNBP, ΔNBPPV, and ΔNBPAR1) per pixel using the mblm function in R55. 513 Wilcoxon rank sum tests were performed to test the significance of the Theil-Sen’s trends. 514 Even though we present results for both atmospheric inversions (CAMS, CarboScope), we 515 combined their results into one fusion dataset in order to further highlight those results for 516 which both inversion models agreed (Supplementary Figure 10 and Supplementary 517 Figure 12). The combination was performed by calculating the average value of each of the 518 abovementioned variables per pixel. We followed the same approach to provide a TRENDY 519 dataset combining all simulations following similar studies1. 520 Trends in average global NBP, NBPPV, and NBPAR1 (from 1981 to 2018 for NBP, from 1991 521 to 2018 for NBPPV, and NBPAR1) were also estimated using the Theil-Sen’s approach. We 522 also combined both atmospheric inversions and the TRENDY simulations to provide the 523 average results derived from atmospheric inversions and process-based models. In this 524 case, the combination of products was performed by calculating the average between 525 products per year and then estimating the temporal trends. We additionally calculated an 526 aggregated metric of trends in NBPPV and NBPAR1. To do so, we first calculated the ratio 527 between each value and the maximum absolute value of each of the variables (ΔNBPPV and 528 ΔNBPAR1) and then we combined them by summing their relative values per pixel. Next, we 529 selected those pixels in which both, ΔNBPPV and ΔNBPAR1, showed a positive trend and 530 those that showed a negative trend, and investigated their respective trends in NBP in the 531 same way that we proceeded for global NBP. 532 Results obtained from the atmospheric inversions were compared with those from nine 533 stations monitoring atmospheric CO2 concentrations distributed from north to south of the 534 Pacific Ocean that comprised the full period of this study. We selected this subset of 535 monitoring stations to minimise the influence of anthropogenic emissions in the signal of 536 atmospheric CO2 concentration. These data were downloaded from the Scripps CO2 537 programme: https://scrippsco2.ucsd.edu/data/atmospheric_co2/56. We calculated the 538 annual amplitude in CO2 concentration (maximum minus minimum) using monthly data 539 following20,36 as a proxy of the net global C uptake capacity. An 11-year moving average 540 was then used to determine whether the annual amplitude of CO2 concentrations, its 541 interannual variability (PV), and the monthly temporal autocorrelation (AR1) of CO2 542 concentrations changed between 1981 and 2018. Theil-Sen slopes were also used to 543 calculate trends in CO2 amplitude, its interannual variability, and their monthly temporal 544 autocorrelation. 545 Drivers of NBP 546 We used temperature and precipitation data from the Climatic Research Unit TS4.03 data 547 set35 to calculate mean annual temperature (MAT), mean annual precipitation (MAP), their 548 temporal variabilities (MATPV and MAPPV, average of interannual monthly temporal 549 variability), their monthly temporal autocorrelations (MATAR1 and MAPAR1), as well as the 550 temporal trends of all these metrics (ΔMAT, ΔMAP, ΔMATPV, ΔMAPPV, ΔMATAR1, and 551 ΔMAPAR1) for each pixel following the same procedure established for NBP (see above). In 552 order to investigate the controls of the spatial variability in aggregated trends in NBPPV-AR1, 553 we also calculated the aggregated indices for ΔMATPV-AR1 and ΔMAPPV-AR1 following the 554 same methodology used for NBP (see above). Land-use changes were extracted from land- 555 use harmonisation2 maps (LUH2, http://luh.umd.edu/data.shtml). We calculated the percent 556 coverages of forests, croplands, and urban areas per pixel, and the change in these 557 percentages between 1981 and 2015. We calculated mean total atmospheric N deposition 558 per pixel derived from34, covering the study period. Biodiversity data were extracted from 559 an interpolated gridded global map of vascular plant biodiversity33, providing one datum per 560 pixel including information on the current number of vascular plant species. We could not 561 include the effects of biodiversity loss in our analyses because global gridded time series 562 of plant biodiversity are not available. 563 Statistical analyses 564 We determined how NBP, NBPPV, NBPAR1, and their trends (ΔNBP, ΔNBPPV, ΔNBPAR1 and 565 the combined ΔNBPPV-AR1) were spatially correlated with their drivers using spatial 566 generalised linear mixed models with a Leroux conditional autoregressive prior57 and the 567 S.CARleroux function in the CARBayes R package58. All response variables had one datum 568 per pixel, representing aggregated information over the entire study period (e.g., NBP 569 indicates mean annual NBP from 1981 to 2018). Models predicting spatial variability in NBP, 570 NBPPV, and NBPAR1 included biodiversity as a second-order polynomial function to account 571 for the nonlinearities usually found between biodiversity and productivity59 (biodiversity + 572 biodiversity2), mean atmospheric N deposition, the interaction between biodiversity and N 573 deposition, and the percentages of the area covered by forests and agricultural and urban 574 areas as predictors. MAT and MAP were also included for NBP models, MATPV, MAPPV for 575 NBPPV models, and MATAR1, MAPAR1 for NBPAR1 models. We built an additional model 576 including the interaction between biodiversity and mean and variability in climate 577 (biodiversity:MAT + biodiversity:MAP + biodiversity:MATPV + biodiversity:MAPPV) and the 578 interactions MAT:MAP and MATPV:MAPPV to further test that the relationship between NBP 579 and biodiversity was not spuriously emerging due to its relationship with climate 580 (Supplementary information – Model summaries, section 8). As indicated by our results, our 581 analyses were able to successfully discern the effect of biodiversity from that of climate 582 (Figure 3). The fact that the relationship between mean NBP and biodiversity emerging 583 from inversions and the TRENDY ensemble (with no parameterisation for biodiversity) were 584 substantially different increases the likelihood of this relationship not being spurious. Models 585 predicting ΔNBP also included the trends in MAT, MAP and the trends in land-use change 586 as predictors. Similarly, models predicting ΔNBPPV and ΔNBPAR1 included trends in 587 temperature and precipitation (ΔMAT and ΔMAP) and their temporal variability (ΔMATPV, 588 ΔMAPPV) or autocorrelation (ΔMATAR1, ΔMAPAR1) and trends in land-use changes, but not 589 average MAT or MAP or their temporal variability or autocorrelation. Models predicting the 590 combined trends of ΔNBPPV-AR1 included a second-order biodiversity as a second-order 591 polynomial, N deposition, average and trends in climate, the aggregated trends of MAT PV- 592 AR1 and MAPPV-AR1, average land-use and its temporal trends. We additionally tested 593 whether NBP and ΔNBP were spatially related to increasing ΔNBPPV-AR1 by fitting two 594 models in which NBP and ΔNBP were the response variables and the aggregated metric of 595 trends in ΔNBPPV-AR1 (see description above) was their predictor. Models were fitted using 596 normalised variables (mean=0, sd=1), but effect plots were rescaled to their original units 597 to facilitate the interpretation of the results. All analyses were performed using R statistical 598 software v.3.6.360. 599 600 Methods references 601 47. Chevallier, F. et al. CO 2 surface fluxes at grid point scale estimated from a global 602 21 year reanalysis of atmospheric measurements. J. Geophys. Res. 115, D21307 603 (2010). 604 48. Chevallier, F. et al. Toward robust and consistent regional CO2 flux estimates from 605 in situ and spaceborne measurements of atmospheric CO2. Geophys. Res. Lett. 606 41, 1065–1070 (2014). 607 49. Rödenbeck, C., Houweling, S., Gloor, M. & Heimann, M. CO2 flux history 1982– 608 2001 inferred from atmospheric data using a global inversion of atmospheric 609 transport. Atmos. Chem. Phys. 3, 1919–1964 (2003). 610 50. Rödenbeck, C., Zaehle, S., Keeling, R. & Heimann, M. How does the terrestrial 611 carbon exchange respond to interannual climatic variations? A quantification based 612 on atmospheric CO2 data. Biogeosciences Discuss. 1–22 (2018) doi:10.5194/bg- 613 2018-34. 614 51. Sitch, S. et al. Recent trends and drivers of regional sources and sinks of carbon 615 dioxide. Biogeosciences 12, 653–679 (2015). 616 52. Fernández‐Martínez, M. & Peñuelas, J. Measuring temporal patterns in ecology: 617 The case of mast seeding. Ecol. Evol. 11, 2990–2996 (2021). 618 53. Wood, S. N. Generalized Additive Models: An introduction with R (2nd edition). 619 (Chapman and Hall/CRC, 2017). doi:10.1201/9781315370279. 620 54. Ohlson, J. A. & Kim, S. Linear valuation without OLS: the Theil-Sen estimation 621 approach. Review of Accounting Studies vol. 20 (2015). 622 55. Lukasz Komsta. mblm: Median-Based Linear Models. (2013). 623 56. Keeling, C. D. et al. Atmospheric CO2 and 13CO2 exchange with the terrestrial 624 biosphere and oceans from 1978 to 2000: observations and carbon cycle 625 implications. in A History of Atmospheric CO2 and its effects on Plants, Animals, 626 and Ecosystems (eds. Ehleringer, J. R., Cerling, T. E. & Dearing, M. D.) 83–113 627 (Springer Verlag, 2005). 628 57. Leroux, B. G., Lei, X. & Breslow, N. Estimation of Disease Rates in Small Areas: A 629 new Mixed Model for Spatial Dependence. in Statistical Models in Epidemiology, 630 the Environment and Clinical Trials (eds. Halloran, M. & Berry, D.) 179–191 631 (Springer-Verlag, 2000). doi:10.1007/978-1-4612-1284-3_4. 632 58. Lee, D. CARBayes: An R package for Bayesian spatial modeling with conditional 633 autoregressive priors. J. Stat. Softw. 55, 1–24 (2013). 634 59. Gonzalez, A. et al. Scaling‐up biodiversity‐ecosystem functioning research. Ecol. 635 Lett. 15, ele.13456 (2020). 636 60. R Core Team. R: A Language and Environment for Statistical Computing. (2020). 637 638 Data availability 639 Data supporting the findings of this study are available in the following open repositories. 640 CAMS: https://ads.atmosphere.copernicus.eu/cdsapp#!/dataset/cams-global-greenhouse- 641 gas-inversion 642 CarboScope: http://www.bgc-jena.mpg.de/CarboScope/ 643 Atmospheric CO2 concentration: https://scrippsco2.ucsd.edu/data/atmospheric_co2/ 644 Data from the TRENDY ensembles can be provided upon request from 645 http://globalcarbonbudget.data.org/index.html. Data to reproduce Figures 1 and 2 are 646 provided as Source Data. Code and data to perform the statistical analyses, calculations 647 and figures is publicly available at Figshare: 648 https://doi.org/10.6084/m9.figshare.17081717.v5. 649 Code availability 650 Code and data to perform the statistical analyses, calculations and figures is publicly 651 available at Figshare: https://doi.org/10.6084/m9.figshare.17081717.v5. 652 Acknowledgements 653 This research was funded by the Spanish Government project PID2019-110521GB-I00, the 654 Fundación Ramón Areces project CIVP20A6621, the Catalan government project 655 SGR2017-1005, and the European Research Council project ERCSyG-2013-610028 656 IMBALANCE-P. M.F-M. was supported by a postdoctoral fellowship of the Research 657 Foundation-Flanders (FWO) and is now supported by a fellowship from ”la Caixa” 658 Foundation (ID 100010434), code: LCF/BQ/PI21/11830010. We acknowledge the Scripps 659 CO2 program for providing the records of atmospheric CO2. 660 Author contribution statement 661 MFM, JP, and IAJ planned and designed the research. FC, CR, SS, PF, VA, DG, AKJ, DLL 662 and PCM provided the data. MFM analysed the data. All authors contributed substantially 663 to the writing of the manuscript. 664 665 Author information statement 666 The authors declare no competing interest. Correspondence and requests for materials 667 should be addressed to: Marcos Fernández-Martínez: m.fernandez@creaf.uab.cat 668