# Diagnosing destabilization risk in global land carbon sinks

Marcos Fernández-Martínez, Josep Peñuelas, Frederic Chevallier, Philippe Ciais, Michael Obersteiner, Christian Rödenbeck, Jordi Sardans, Sara Vicca, Hui Yang, Stephen Sitch, et al.

► To cite this version:

Marcos Fernández-Martínez, Josep Peñuelas, Frederic Chevallier, Philippe Ciais, Michael Obersteiner, et al. Diagnosing destabilization risk in global land carbon sinks. Nature, 2023, 615, pp.848-853. 10.1038/s41586-023-05725-1. hal-04041037

HAL Id: hal-04041037

https://hal.science/hal-04041037v1

Submitted on 2 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
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 © IN COPYRIGHT license

<figure><img src="image_2.png" /></figure>
1 **Diagnosing destabilisation risk in global land carbon sinks**

2 **Authors:** Marcos Fernández-Martínez¹,²,³, Josep Peñuelas²,⁴, Frederic Chevallier⁵, Philippe

3 Ciais⁵, Michael Obersteiner⁶,⁷, Christian Rödenbeck⁸, Jordi Sardans²,⁴, Sara Vicca¹, Hui Yang⁵,

4 Stephen Sitch⁹, Pierre Friedlingstein¹⁰, Vivek K. Arora¹¹, Daniel S. Goll⁵, Atul K. Jain¹², Danica

5 L. Lombardozzi¹³, Patrick C. McGuire¹⁴, Ivan A. Janssens¹.

6 **Affiliations:**

7 ¹PLECO (Plants and Ecosystems), Department of Biology, University of Antwerp, 2610 Wilrijk,

8 Belgium.

9 ²CREAF, E08193 Bellaterra (Cerdanyola del Vallès), Catalonia, Spain

10 ³BEECA-UB, Department of Evolutionary Biology, Ecology and Environmental Sciences,

11 University of Barcelona, E08028 Barcelona, Catalonia, Spain

12 ⁴CSIC, Global Ecology Unit, CREAF-CSIC-UAB, Bellaterra, 08193 Barcelona, Catalonia,

13 Spain

14 ⁵Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ,

15 Université Paris-Saclay, Gif-sur-Yvette, France

16 ⁶International Institute for Applied Systems Analysis (IIASA), Laxenburg, Austria

17 ⁷School of Geography and the Environment, University of Oxford, Oxford, UK

18 ⁸Department of Biogeochemical Systems, Max Planck Institute for Biogeochemistry, Jena,

19 Germany

20 ⁹College of Life and Environmental Sciences, University of Exeter, Exeter, UK

21 ¹⁰College of Engineering, Mathematics, and Physical Sciences, University of Exeter, Exeter,

22 United Kingdom

23 ¹¹Canadian Centre for Climate Modelling and Analysis, Climate Research Division,

24 Environment and Climate Change Canada, Victoria, BC, Canada

25 ¹²Department of Atmospheric Sciences, University of Illinois, Urbana, IL 61801, USA

26 ¹³Climate and Global Dynamics Laboratory, National Center for Atmospheric Research,

27 Boulder, CO 80302, United States of America
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**

Global net land carbon uptake, or net biome production (NBP), has increased during recent decades¹. Whether its temporal variability and autocorrelation have changed during this period, however, remains elusive, even though an increase in both could indicate an increased potential for a destabilised carbon sink²,³. Here we investigate the trends and controls of net terrestrial carbon uptake and its temporal variability and autocorrelation from 1981 to 2018 using two atmospheric inversion models, the amplitude of the seasonal cycle of atmospheric CO₂ concentration derived from nine monitoring stations distributed across the Pacific Ocean and dynamic global vegetation models. We find that annual NBP and its interdecadal variability increased globally while temporal autocorrelation decreased. We observe a separation of regions characterized by increasingly variable NBP, associated with warm regions and increasingly variable temperatures, lower and weaker trends in NBP, and regions where NBP became stronger and less variable. Plant species richness presented a concave-down parabolic spatial relationship with NBP and its variability at the global scale while nitrogen deposition generally increased NBP. Increasing temperature and its increasing variability appears as the most important drivers of declining and increasingly variable NBP. Our results show increasing variability of NBP regionally that can be mostly attributed to climate change, 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

Positive carbon-climate feedbacks have the potential to accelerate climate change and might compromise the attainability of ambitious climate targets such as those set by the Paris agreement⁴. Terrestrial ecosystems are key to the functioning of the global carbon (C) cycle and have increased their productivity and net C uptake during recent decades primarily due to CO₂ fertilisation and forest regrowth¹,⁵⁻⁷. However, land-use change, nutrient limitations and increasing droughts and fires are constraining this potential to sequester C⁸⁻¹¹. Identifying processes that might destabilise net land C uptake (or net biome production, NBP) is of paramount importance for understanding and managing the global C cycle.

Destabilisation is the process of losing stability: that is, when a system loses its ability to return to equilibrium following a disturbance. A symptom of this may be increased variability, as the system spirals away from its current equilibrium point towards a new one. Destabilisation of a dynamical system is usually accompanied by concomitant increases in temporal variability and autocorrelation (i.e., the correlation between consecutive time steps [AR1] also related to reduced resilience²) when AR1 is already positive, because anomalous states of the system create ripples that can get amplified through time rather than compensated as when AR1 is negative. Consequently, increasing temporal variability and autocorrelation have been shown to be potential early-warning signals for abrupt shifts in ecosystems²,³,¹²,¹³. To date, changes in NBP temporal variability and autocorrelation have not yet been investigated even though increased stress and changes in the frequency and intensity of extreme weather events are likely to alter the temporal patterns of NBP¹⁴⁻¹⁶ due to i) a cumulative negative effect of extreme events on ecosystem functioning, ii) increasing climate variability and, iii) decreasing ecosystem resilience due to increased stress. We, thus, hypothesised that regions experiencing increasing trends in climate variability will also experience increasing variability in their NBP, and that increase in variability may indicate that a destabilisation of the net land C uptake is occurring. This finding could serve as an early warning signal for abrupt shifts in ecosystem's functions² that might lead to regime shifts in the Earth's biosphere¹⁷,¹⁸. Even if these changes do not occur at the global scale, increasing temporal variability and autocorrelation in several regions of the globe (e.g. Amazon basin, boreal ecosystems with permafrost) could have a profound impact on the global C balance and a knock-on effect on other ecosystem functions.

Climate is the primary control on NBP in terrestrial ecosystems worldwide in space and time¹⁹,²⁰, together with soil nutrient availability¹⁰,²¹, atmospheric nitrogen (N) deposition, land use and management, and increasing atmospheric CO₂¹. Generally, gross C fluxes (such as photosynthesis and respiration) are larger in the tropics where high temperature coincides with
sufficient precipitation to enable long growing seasons. However, net C uptake tends to be higher in temperate regions due to higher nutrient availability¹⁹,²¹⁻²³. Nutrient availability and N deposition have indeed been shown to increase net land C uptake²¹,²². Another factor believed to be important in determining ecosystem functioning is plant biodiversity²⁴. A large body of evidence indicates its role in promoting ecosystem productivity and stability²⁵,²⁶. Plant biodiversity, however, has been included far less in studies of ecosystem C cycling, presumably due to the difficulty of acquiring good data on species diversity. The few studies that included biodiversity, however, showed relatively modest correlations with C fluxes²⁷⁻²⁹. The role of biodiversity in the global terrestrial NBP has not yet been explored even though the Earth's biosphere is losing biodiversity at an unprecedented speed³⁰ and those changes are expected to alter ecosystem functioning.

The aim of this study was to quantify the trends in global NBP and its intradecadal temporal variability (quantified by the proportional variability index, PV³¹,³²) and autocorrelation (AR1) to see whether changes in these variables suggest that NBP at global or regional scales is being destabilised. We further tested whether regions with increasing variability and AR1 in NBP showed differential trends in annual NBP. We additionally investigated the spatial correlation between NBP metrics and plant biodiversity (derived from a global map of plant species richness³³), atmospheric total N deposition³⁴, climate³⁵, and land use (land-use harmonisation² maps). To realise these objectives, we used estimates of NBP derived from the two longest CO₂ atmospheric inversions (CAMS, CarboScope) and the amplitude of the seasonal cycle of atmospheric CO₂ derived from nine monitoring stations distributed from south to north of the Pacific Ocean for 1981 to 2018. We additionally compared these results with the output from an ensemble of 12 dynamic global vegetation models (TRENDY) to explore how well these state-of-the-art models predict the spatial and temporal patterns in NBP simulated by atmospheric inversion models.

## 115 Results

### *Temporal patterns of net land C uptake*

Global NBP derived from atmospheric inversions increased from 5.6 ± 2.0 gC m⁻² y⁻¹ during 1981 – 1990 (mean ± standard error) to 13.8 ± 1.4 gC m⁻² y⁻¹ during 2009 – 2018 over the global land area excluding Antarctica (Figure 1). This represents an overall increase of 145% and an annual linear increase of 0.24 ± 0.08 gC m⁻² y⁻² (*P*<0.001). Both atmospheric inversions used (CAMS and CarboScope) identified annual increases in NBP (Supplementary Figure 1 and Supplementary Figure 2, 0.32 ± 0.09, and 0.18 ± 0.09 gC m⁻² y⁻² respectively, *P*<0.001, N=38). The positive trend in NBP shown by inversions is similar in magnitude to the trend identified by the TRENDY ensemble (0.10 ± 0.07 gC m⁻² y⁻², *P*<0.001; 147% when comparing
the periods 1981 – 1990 and 2009 – 2018). Both observations and models, however, showed a flattening trend during the last decade (Figure 1c). In parallel with the increased NBP throughout the entire study period and in agreement with previous literature³⁶,³⁷, the amplitude of the seasonal atmospheric CO₂ concentration increased by 0.027 ± 0.004 ppm y⁻¹ (*P*<0.001, 16.6% during the study period) (Supplementary Figure 3a).

Interannual NBP variability (NBP<sub>PV</sub>) derived from the combination of atmospheric inversions increased globally by 7.2% over the entire period (*P*<0.001) (Figure 1f) (CAMS 14.9%, *P*<0.001; CarboScope: -1.3%, *P*<0.001). The origin of the discrepancy between the inversions remains unclear. Though the inversions differ in whether they use a growing number of measurement stations (CAMS) or a constant station set (CarboScope), the discrepancy is unlikely to be related to this because a test inversion with a growing station set in the CarboScope inversion did not yield an increase in NBP<sub>PV</sub>. Temporal variability in the amplitude of the seasonal atmospheric CO₂ concentration also increased (by 35%, *P*<0.001) during the study period (Supplementary Figure 3b). TRENDY NBP did not simulate this increase in variability, showing a negative trend of -1.7% in NBP<sub>PV</sub> similar to that from the CarboScope inversion.

Global NBP<sub>AR1</sub> significantly decreased over time for NBP derived from atmospheric inversions (Figure 1i) and for the average monthly CO₂ concentrations across the Pacific Ocean (Supplementary Figure 3c). TRENDY NBP<sub>AR1</sub>, also showed a significant negative trend representing a reduction of 2.9% over the study period. Trends of contrasting sign, however, were significant amongst atmospheric measurement stations and atmospheric inversions (CAMS: -10.0%, *P*<0.001; CarboScope: 4.7%, *P*<0.001). Several regions showed increases in the temporal autocorrelation between consecutive months (NBP<sub>AR1</sub>) derived from atmospheric inversions that are of a similar magnitude (ca. 0.2 over three decades) to those previously suggested to precede abrupt shifts in climate datasets and simulations²,³. Additional information on how changes in NBP<sub>PV</sub> and NBP<sub>AR1</sub> may affect annual NBP and their limitations can be found in Supplementary text, Section 1.

Our analyses identified several regions of potential concern given their concomitant increase in NBP<sub>PV</sub> and NBP<sub>AR1</sub> (Figure 2a), such as eastern Africa, the Mediterranean region, the west coasts of North and Central America, India, and southeast Asia. Regions with increasing NBP<sub>PV</sub> and NBP<sub>AR1</sub> had statistically lower NBP and experienced a much less pronounced increase in NBP over time (0.15 ± 0.06 gC m⁻² y⁻¹, *P*<0.001) compared to regions where NBP<sub>PV</sub> and NBP<sub>AR1</sub> decreased (0.73 ± 0.12 gC m⁻² y⁻¹, *P*<0.001) (Figure 2b-d, Supplementary Figure 4). Stronger increases in NBP over regions with decreasing NBP<sub>PV</sub> and NBP<sub>AR1</sub>,
159

compared to those where they increased, were also evident from both atmospheric inversions when analysed separately (Supplementary Figure 5).

160

161

Global temperature and precipitation increased during the study period (Supplementary Figure 6a and d) even though there was considerable spatial variability in those trends (Supplementary Figure 7a and d). Temperature and precipitation interannual variability, however, decreased significantly at the global scale, despite the reported increase in extreme weather¹⁴. Temporal autocorrelation of monthly temperature and precipitation, instead, increased slightly, albeit only significantly for temperature. A concomitant increase in temporal variability and autocorrelation of temperature was evident in several regions (e.g., Eurasia, Australia, Central America) but those were not so obvious for precipitation (Supplementary Figure 8).

162

163

164

165

166

167

168

169

170

## Controls of NBP, variability and autocorrelation

171

172

173

174

175

176

NBP, derived from atmospheric inversions, had a concave-down parabolic relationship with biodiversity (Figure 3), increasing from low to intermediate values of biodiversity and decreasing at high biodiversity (Supplementary Figure 9a). Our analyses also identified a significant positive interaction between biodiversity and N deposition accounting for NBP. The concave-down relationship between biodiversity and NBP included mainly positive NBP values across regions with high atmospheric N deposition (850 mg N m⁻² y⁻¹) and mainly negative NBP values at low atmospheric N deposition (100 mg N m⁻² y⁻¹) (Supplementary Figure 9a). Similarly, the positive correlation between N deposition and NBP was stronger in regions with higher plant biodiversity (Supplementary Figure 9b). Opposite relationships emerged when analysing TRENDY NBP.

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

NBP<sub>PV</sub> derived from inversions also showed a concave-down parabolic relationship with plant biodiversity, peaking at intermediate to high values of biodiversity (Supplementary Figure 9d). We also found an interaction between plant biodiversity and N deposition in their relationship with NBP<sub>PV</sub>: the relationship between NBP<sub>PV</sub> and biodiversity differed between areas receiving low N deposition and areas receiving high N deposition (Supplementary Figure 9d). No relationship between N deposition and NBP<sub>PV</sub> was found in regions with high biodiversity, whereas a positive relationship occurred in regions with low biodiversity (Supplementary Figure 9e). Interannual variability in temperature and precipitation were positively correlated with NBP<sub>PV</sub> (Figure 3). In this case, the patterns for the TRENDY ensemble and the inversion models matched very well. Again, NBP<sub>AR1</sub> derived from inversions had a concave-down parabolic relationship with biodiversity (Supplementary Figure 9g), but here no interaction occurred between the effects of biodiversity and N deposition, which was negatively correlated with NBP<sub>AR1</sub>. NBP<sub>AR1</sub> was also positively correlated with temperature
194

AR1, the only result that match those from TRENDY NBP<sub>AR1</sub> (**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, NBP<sub>PV</sub> and NBP<sub>AR1</sub> derived from atmospheric inversions

198

were all correlated with N deposition, climate, and land use (**Figure 3** and **Supplementary Figure 11**).

199

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 NBP<sub>PV</sub> and NBP<sub>AR1</sub> (**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 NBP<sub>PV</sub> were spatially associated with low increases in annual

206

temperature (**Supplementary Figure 11e**). Increases in NBP<sub>PV</sub> were also found to be more

207

likely in regions showing increases in temperature temporal variability (**Supplementary Figure 11f**).

208

Increasing NBP<sub>AR1</sub> 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 NBP<sub>AR1</sub>, but not for NBP<sub>PV</sub>. 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 NBP<sub>PV</sub> and NBP<sub>AR1</sub>.

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

NBP<sub>PV</sub> derived from atmospheric inversions, however, was very similar to the one emerging
229
230
231
232

from TRENDY, which suggests that factors other than biodiversity may be driving this relationship. The positive relationship between biodiversity and NBP<sub>AR1</sub> found here has never been reported before and further research is needed to understand the mechanisms behind this relationship.

233
234
235
236
237
238
239
240
241
242
243
244
245
246

The concave-down parabolic relationships of biodiversity with NBP and NBP<sub>PV</sub> differ from the majority of biodiversity-productivity and stability relationships reported in the literature: positive asymptotic for productivity and negative for variability<sup>38,39</sup>. This difference in the biodiversity-NBP relationship may result from two opposing ecosystem processes, photosynthesis and respiration, because both are likely to be enhanced by biodiversity<sup>27,40,41</sup>. The concave-down relationship would then suggest that the positive effect of biodiversity on respiration overshadows the positive effect of photosynthesis in regions with high biodiversity. However, our biodiversity data, like in similar studies<sup>39</sup>, was restricted to species richness and did not include information on species abundance, their individual contribution to NBP, or traits to allow the calculation of actual species diversity or functional diversity, often better indicators of ecosystem functioning than species richness<sup>42</sup>. Unfortunately, this information is not available at the global scale. Including actual diversity in future analyses could lead to different results to those reported here. Additionally, future efforts are needed to understand how biodiversity loss, not included here, will impact global carbon balance.

247
248
249
250
251
252
253
254
255
256
257

Interestingly, our analyses indicated that the effect of biodiversity on NBP depended on atmospheric N deposition, and vice versa (Figure 3). The effect of atmospheric N deposition on NBP was mainly positive but was stronger in regions with higher biodiversity, further supporting the premise that biodiversity promotes ecosystem functions such as N uptake<sup>43</sup>. Regions with higher N deposition also had larger increases in NBP over time (Supplementary Figure 11a), supporting previous findings suggesting a stronger CO<sub>2</sub> fertilisation effect in regions with higher N deposition<sup>1</sup>. On the other hand, N deposition was negatively related to trends in NBP<sub>PV</sub>, NBP<sub>AR1</sub> and the aggregated trend of NBP<sub>PV-AR1</sub> (Supplementary Figure 12), which suggests that N addition may ameliorate nutrient imbalances derived from increasing atmospheric CO<sub>2</sub> concentrations<sup>8</sup> and prevent ecosystem functioning from becoming more variable.

258

## Climate and changing land C uptake

259
260
261
262
263

Our results clearly indicate that NBP decreased in regions with warm climates and where warming has been most pronounced (Supplementary Figure 11). These findings support results from previous studies¹ and further suggests that increasing droughts and heat waves limit C sequestration by terrestrial ecosystems⁴⁴. Our analyses revealed that higher NBP<sub>PV</sub>, NBP<sub>AR1</sub> and their trends were associated, respectively, with climates showing higher and
264

increasing temporal variability and autocorrelation (**Supplementary Figure 9**). We also found that concomitantly increasing trends in $NBP_{PV-AR1}$ were positively related to increasing temperature variability and autocorrelation (**Supplementary Figure 12**). These results support our initial hypothesis stating that climate change may be the main contributing factor for changing the temporal behaviour of land C sinks.

265

266

267

268

269

## Implications of altered land C sinks

270

271

272

273

274

275

276

277

278

279

280

281

Compelling evidence from atmospheric inversions (**Figure 1**) and the annual amplitude of the seasonal CO₂ cycle (**Supplementary Figure 3**) suggests that global net land C uptake is now larger and, most likely, more variable than three decades ago, while its temporal autocorrelation ($NBP_{AR1}$), instead, has significantly decreased. The increase in global NBP and the reduction in $NBP_{AR1}$ were well identified by TRENDY NBP, including the evidence of a flattening trend in NBP during the last two decades. Our results for the increase in global NBP and a potential recent saturation are consistent with previous findings of increases in global productivity and NBP and a recent decline of the CO₂ fertilisation effect¹,⁵,⁸,²⁰. Even though the main emerging spatial controls for CAMS and CarboScope atmospheric inversions in NBP were in agreement (**Figures S8 and S10**), trends at the global scale for $NBP_{PV}$ and $NBP_{AR1}$ were significantly different (**Supplementary Figure 13**), which calls for caution in the interpretation of these results due to their uncertainty.

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

We observed a bifurcation of regions characterised by concomitant increases and decreases in $NBP_{PV}$ and $NBP_{AR1}$ (**Figure 2**). This bifurcation relates to differences in their mean annual temperature, their increase in annual temperature, their average N deposition loads and their percentage of forests and crops (**Supplementary Figure 12**). The observed increase in the variability and autocorrelation of net C uptake in several regions of the planet (**Figure 2a**) is concerning because of the implications it can have on the stability of their ecosystems. First, increasing $NBP_{PV}$ and $NBP_{AR1}$ is indicative of increasing variability and reducing resilience³,⁴⁵ potentially in many ecosystem processes: from photosynthesis and respiration²⁷ to cascading effects on animals and decomposers⁴⁶. Second, an increase in variability in NBP implies that ecosystem C balance is less predictable over time, which is troublesome for projections of future climate change. Third, regions showing a combined increase in $NBP_{PV}$ and $NBP_{AR1}$ had a consistently lower NBP and lower increases in NBP than regions where $NBP_{PV}$ and $NBP_{AR1}$ decreased (**Figure 2b-d**), suggesting that these increases $NBP_{PV}$ and $NBP_{AR1}$ reflect dynamical instability in the carbon system rather than some other cause (such as changes in the external forcing of the system). This destabilization may jeopardise future C sequestration. Our analyses, however, could not determine the mechanisms driving the observed changes in $NBP_{PV}$ and $NBP_{AR1}$. More importantly, we could not determine whether the observed increases
299 NBP<sub>PV</sub> and NBP<sub>AR1</sub> 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 CO<sub>2</sub>
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 functioning<sup>2,3</sup>. The increase of NBP<sub>PV</sub> and NBP<sub>AR1</sub> 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 (2009).
313

314 3. Dakos, V. et al. Slowing down as an early warning signal for abrupt climate change. *Proc. Natl. Acad. Sci. U. S. A.* **105**, 14308–14312 (2008).
315

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 (2016).
319

320 6. Bastos, A. et al. Contrasting effects of CO2 fertilization, land-use change and warming on seasonal amplitude of Northern Hemisphere CO2 exchange. *Atmos. Chem. Phys.* **19**, 12361–12375 (2019).
321
322

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).
23. Ciais, P. et al. Five decades of northern land carbon uptake revealed by the interhemispheric CO2 gradient. *Nature* **568**, 221–225 (2019).

24. Tilman, D., Lehman, C. L. & Thomson, K. T. Plant diversity and ecosystem productivity: Theoretical considerations. *Proc. Natl. Acad. Sci.* **94**, 1857–1861 (1997).

25. de Mazancourt, C. et al. Predicting ecosystem stability from community composition and biodiversity. *Ecol. Lett.* **16**, 617–625 (2013).

26. Sakschewski, B. et al. Resilience of Amazon forests emerges from plant trait diversity. *Nat. Clim. Chang.* **6**, 1032–1036 (2016).

27. Fernández-Martínez, M. et al. The role of climate, foliar stoichiometry and plant diversity on ecosystem carbon balance. *Glob. Chang. Biol.* **26**, 7067–7078 (2020).

28. Musavi, T. et al. Stand age and species richness dampen interannual variation of ecosystem-level photosynthetic capacity. *Nat. Ecol. Evol.* **1**, 0048 (2017).

29. Anderegg, W. R. L. et al. Hydraulic diversity of forests regulates ecosystem resilience during drought. *Nature* **561**, 538–541 (2018).

30. Díaz, S. et al. IPBES, 2019. *Summary for policymakers of the global assessment report on biodiversity and ecosystem services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services*. (2019) doi:10.1111/padr.12283.

31. Heath, J. P. Quantifying temporal variability in population abundances. *Oikos* **115**, 573–581 (2006).

32. Fernández-Martínez, M., Vicca, S., Janssens, I. A., Martín-Vide, J. & Peñuelas, J. The consecutive disparity index, D, as measure of temporal variability in ecological studies. *Ecosphere* **9**, e02527 (2018).

33. Kreft, H. & Jetz, W. Global patterns and determinants of vascular plant diversity. *Proc Natl Acad Sci U S A* **104**, 5925–5930 (2007).

34. Ackerman, D. E., Chen, X. & Millet, D. B. Global nitrogen deposition (2°×2.5° grid resolution) simulated with GEOS-Chem for 1984-1986, 1994-1996, 2004-2006, and 2014-2016. (2018) doi:10.13020/D6KX2R.

35. Harris, I., Jones, P. D. D., Osborn, T. J. J. & Lister, D. H. H. Updated high-resolution grids of monthly climatic observations - the CRU TS3.10 Dataset. *Int. J. Climatol.* **34**, 623 – 642 (2013).
390
391

36. Graven, H. D. et al. Enhanced seasonal exchange of CO2 by northern ecosystems since 1960. *Science* **341**, 1085–9 (2013).

392
393

37. Wang, K. et al. Causes of slowing-down seasonal CO2 amplitude at Mauna Loa. *Glob. Chang. Biol.* **26**, 4462–4477 (2020).

394
395

38. Tilman, D., Reich, P. B. & Knops, J. M. H. Biodiversity and ecosystem stability in a decade-long grassland experiment. *Nature* **441**, 629–32 (2006).

396
397

39. Liang, J. et al. Positive biodiversity-productivity relationship predominant in global forests. *Science* (80-.). **354**, aaf8957-aaf8957 (2016).

398
399

40. Gessner, M. O. et al. Diversity meets decomposition. *Trends Ecol. Evol.* **25**, 372–380 (2010).

400
401
402

41. Peguero, G. et al. Fast attrition of springtail communities by experimental drought and richness–decomposition relationships across Europe. *Glob. Chang. Biol.* **25**, 2727–2738 (2019).

403
404

42. Díaz, S. & Cabido, M. Vive la différence: Plant functional diversity matters to ecosystem processes. *Trends Ecol. Evol.* **16**, 646–655 (2001).

405
406

43. Cardinale, B. J. Biodiversity improves water quality through niche partitioning. *Nature* **472**, 86–91 (2011).

407
408

44. Ciais, P. et al. Europe-wide reduction in primary productivity caused by the heat and drought in 2003. *Nature* **437**, 529–533 (2005).

409
410

45. Scheffer, M. *Critical Transitions in Nature and Society*. (Princeton University Press, 2009).

411
412

46. Ostfeld, R. & Keesing, F. Pulsed resources and community dynamics of consumers in terrestrial ecosystems. *Trends Ecol. Evol.* **15**, 232–237 (2000).

413
**Figure 1: Global distribution of net biome production (NBP), its interannual variability (NBP<sub>PV</sub>) and temporal autocorrelation (NBP<sub>AR1</sub>), and their trends from 1981 to 2018.** Panel a shows the average global distribution of NBP (positive values indicate C sinks) derived from CAMS and CarboScope atmospheric inversions. Panel b shows the temporal trends in NBP (ΔNBP), and panel c shows the temporal change in NBP derived from atmospheric inversions and 12 DGVMs (TRENDY). Panels d, e, and f, and g, h, and i mimic panels a, b, and c, respectively, for interannual variability of NBP (NBP<sub>PV</sub>) and temporal autocorrelation (NBP<sub>AR1</sub>) calculated over an 11-year moving window. Thick lines in panels c, f, and i represent smooth trends estimated with a local regression (shaded areas show the standard error of the trend), and the coefficients at the top of the panels indicate robust Theil-Sen's slopes (± standard error of the slope) with their corresponding p-values, estimated using the non-parametric Wilcoxon Mann-Whitney rank sum test. Values in panels c, f and i, represent spatially-weighted global means. **, *P*<0.01; ***, *P*<0.001. Data to reproduce this figure can be found in Source Data 1.

(See **Figure 1** on the following page)

429
<figure><img src="image_3.png" /></figure>
**Figure 2: Regions with concomitantly increasing NBP<sub>PV</sub> and NBP<sub>AR1</sub> present lower NBP and a lower increase of NBP over time.** Panel a shows regions of the world where concomitant increases (red) and decreases (blue) in NBP<sub>PV</sub> and NBP<sub>AR1</sub> have occurred during the study period. Regions with contrasting trends in NBP<sub>PV</sub> and NBP<sub>AR1</sub> appear in grey. Panel b shows the aggregated annual mean NBP of regions with increasing (red) and decreasing (blue) NBP<sub>PV</sub> and NBP<sub>AR1</sub> (see **Supplementary Figure 4** for information on their global contribution). Smooth trends were estimated with a local regression (shaded indicate the standard error of the trend), and the coefficients at the top of the panel indicate robust Theil-Sen's slopes (± standard error of the slope) and their associated significance (***, *P*<0.001), estimated using the non-parametric Wilcoxon 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 NBP<sub>PV</sub> and NBP<sub>AR1</sub>. Coefficients are standardised and include the 95% credible intervals (CI). See **Methods** for details on how the aggregated trend of NBP<sub>PV</sub> and NBP<sub>AR1</sub> was calculated. Data to reproduce this figure can be found in Source Data 2.

<figure><img src="image_4.png" /></figure>
**Figure 3:** Contribution of biodiversity, N deposition, climate, and land use to the global terrestrial C balance (NBP), its interannual variability (NBP<sub>PV</sub>), its temporal autocorrelation (NBP<sub>AR1</sub>), and their trends. The colour scale indicates the strength of the relationship (standardised β coefficients) between each predictor (bottom) and the response variable (left). Black dots indicate that 95% of the posterior distributions differed from 0. Biodiversity was fitted as a second-order polynomial to account for nonlinearities (biodiversity + biodiversity<sup>2</sup>; the second term indicates the change in the slope in biodiversity as biodiversity increases). Hashed areas indicate relationships not included in the regression models. See **Methods** for further information on model fitting. INV and TRD indicate the averages for all atmospheric inversions (2) and DGVMs (12), respectively. Predictors MAT<sub>x̄</sub>, PV, AR1 and MAP<sub>x̄</sub>, PV, AR1 indicate mean, PV, or AR1 metrics of MAT and MAP matching the metric of the response variable (e.g. NBP<sub>INV</sub> ~ MAT<sub>x̄</sub>; NBP<sub>PV-INV</sub> ~ MAT<sub>PV</sub>). The same applies for ΔMAT<sub>PV</sub>, ΔMAT<sub>AR1</sub>, ΔMAP<sub>PV</sub>, ΔMAP<sub>AR1</sub>. MAT, mean annual temperature; MAP, mean annual precipitation; PV, proportional variability; AR1, temporal autocorrelation with previous month; Δ, trend in a given variable.

<figure><img src="image_5.png" /></figure>
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 CO₂ 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 inversions⁵¹. 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 CO₂, 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 (NBP<sub>PV</sub>)<sup>31,32</sup>, and its monthly temporal autocorrelation at lag 1
484 (NBP<sub>AR1</sub>) for all of the abovementioned datasets. All three indices (NBP, NBP<sub>PV</sub> and NBP<sub>AR1</sub>)
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,
487 following: $PV = \frac{2 \sum m}{n(n-1)}$, where “n” indicates the length of the variable, where “m” is
488 calculated as: $m = 1 - \frac{\min(z_i, z_j)}{\max(z_i, z_j)}$, 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 $z_i$ vs observation $z_j$). 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.
Unlike other metrics of temporal variability, the PV index provides estimates of temporal variability that are independent of the mean of the time series and that have been proven to be robust even when comparing non-normally distributed data sets, thereby overcoming the mathematical drawbacks of similar variability indices such as the standard deviation or the coefficient of variation⁵². NBP<sub>PV</sub> per pixel was estimated as the average of the interannual variabilities of all months (e.g. NBP<sub>PV</sub> = (NBP<sub>PV January</sub> + ... + NBP<sub>PV December</sub>)/12). By using this method, NBP<sub>PV</sub> reflects not only the interannual variability of the annual NBP, but also the variation occurring amongst equal months across the years. Using this method, a year with average NBP yet with anomalously low NBP in one period offset by anomalously high NBP in another, would thus result in high interannual variability rather than low interannual variability. AR1 was estimated using the residuals of generalised additive models (GAM) used to remove the trend and seasonal cycle of the data. In these models response variable was monthly NBP and the predictors were the month of the year (i.e. a factor of 12 levels) and the year, the latter included as a spline smoothed term to account for nonlinear trends over time (**Supplementary Figure 14**). We used the *mgcv* R package to fit the GAM models⁵³. An 11-year moving average from NBP<sub>PV</sub> and NBP<sub>AR1</sub> per pixel was then calculated for all data sets to investigate trends in interannual variability and temporal autocorrelation. Our sensitivity analysis that showed that our results were consistent despite the selection of different window lengths (7, 11 and 15 years, **Supplementary Figure 13**, **Supplementary Table 1**). Robust Theil-Sen's trends⁵⁴ were then calculated for NBP, NBP<sub>PV</sub> and NBP<sub>AR1</sub> (ΔNBP, ΔNBP<sub>PV</sub>, and ΔNBP<sub>AR1</sub>) per pixel using the *mblm* function in R⁵⁵. Wilcoxon rank sum tests were performed to test the significance of the Theil-Sen's trends. Even though we present results for both atmospheric inversions (CAMS, CarboScope), we combined their results into one fusion dataset in order to further highlight those results for which both inversion models agreed (**Supplementary Figure 10 and Supplementary Figure 12**). The combination was performed by calculating the average value of each of the abovementioned variables per pixel. We followed the same approach to provide a TRENDY dataset combining all simulations following similar studies¹.

Trends in average global NBP, NBP<sub>PV</sub>, and NBP<sub>AR1</sub> (from 1981 to 2018 for NBP, from 1991 to 2018 for NBP<sub>PV</sub>, and NBP<sub>AR1</sub>) were also estimated using the Theil-Sen's approach. We also combined both atmospheric inversions and the TRENDY simulations to provide the average results derived from atmospheric inversions and process-based models. In this case, the combination of products was performed by calculating the average between products per year and then estimating the temporal trends. We additionally calculated an aggregated metric of trends in NBP<sub>PV</sub> and NBP<sub>AR1</sub>. To do so, we first calculated the ratio between each value and the maximum absolute value of each of the variables (ΔNBP<sub>PV</sub> and
528

$\Delta$NBP<sub>AR1</sub>) and then we combined them by summing their relative values per pixel. Next, we selected those pixels in which both, $\Delta$NBP<sub>PV</sub> and $\Delta$NBP<sub>AR1</sub>, showed a positive trend and those that showed a negative trend, and investigated their respective trends in NBP in the same way that we proceeded for global NBP.

529

530

531

532

Results obtained from the atmospheric inversions were compared with those from nine stations monitoring atmospheric CO₂ concentrations distributed from north to south of the Pacific Ocean that comprised the full period of this study. We selected this subset of monitoring stations to minimise the influence of anthropogenic emissions in the signal of atmospheric CO₂ concentration. These data were downloaded from the Scripps CO₂ programme: https://scrippsco2.ucsd.edu/data/atmospheric_co2/<sup>56</sup>. We calculated the annual amplitude in CO₂ concentration (maximum minus minimum) using monthly data following<sup>20,36</sup> as a proxy of the net global C uptake capacity. An 11-year moving average was then used to determine whether the annual amplitude of CO₂ concentrations, its interannual variability (PV), and the monthly temporal autocorrelation (AR1) of CO₂ concentrations changed between 1981 and 2018. Theil-Sen slopes were also used to calculate trends in CO₂ amplitude, its interannual variability, and their monthly temporal autocorrelation.

533

534

535

536

537

538

539

540

541

542

543

544

545

## Drivers of NBP

546

547

548

549

550

551

552

553

554

555

556

557

558

559

560

561

562

We used temperature and precipitation data from the Climatic Research Unit TS4.03 data set<sup>35</sup> to calculate mean annual temperature (MAT), mean annual precipitation (MAP), their temporal variabilities (MAT<sub>PV</sub> and MAP<sub>PV</sub>, average of interannual monthly temporal variability), their monthly temporal autocorrelations (MAT<sub>AR1</sub> and MAP<sub>AR1</sub>), as well as the temporal trends of all these metrics ($\Delta$MAT, $\Delta$MAP, $\Delta$MAT<sub>PV</sub>, $\Delta$MAP<sub>PV</sub>, $\Delta$MAT<sub>AR1</sub>, and $\Delta$MAP<sub>AR1</sub>) for each pixel following the same procedure established for NBP (see above). In order to investigate the controls of the spatial variability in aggregated trends in NBP<sub>PV-AR1</sub>, we also calculated the aggregated indices for $\Delta$MAT<sub>PV-AR1</sub> and $\Delta$MAP<sub>PV-AR1</sub> following the same methodology used for NBP (see above). Land-use changes were extracted from land-use harmonisation<sup>2</sup> maps (LUH2, http://luh.umd.edu/data.shtml). We calculated the percent coverages of forests, croplands, and urban areas per pixel, and the change in these percentages between 1981 and 2015. We calculated mean total atmospheric N deposition per pixel derived from<sup>34</sup>, covering the study period. Biodiversity data were extracted from an interpolated gridded global map of vascular plant biodiversity<sup>33</sup>, providing one datum per pixel including information on the current number of vascular plant species. We could not include the effects of biodiversity loss in our analyses because global gridded time series of plant biodiversity are not available.
We determined how NBP, NBP<sub>PV</sub>, NBP<sub>AR1</sub>, and their trends (ΔNBP, ΔNBP<sub>PV</sub>, ΔNBP<sub>AR1</sub> and the combined ΔNBP<sub>PV-AR1</sub>) were spatially correlated with their drivers using spatial generalised linear mixed models with a Leroux conditional autoregressive prior<sup>57</sup> and the S.CARleroux function in the CARBayes R package<sup>58</sup>. All response variables had one datum per pixel, representing aggregated information over the entire study period (e.g., NBP indicates mean annual NBP from 1981 to 2018). Models predicting spatial variability in NBP, NBP<sub>PV</sub>, and NBP<sub>AR1</sub> included biodiversity as a second-order polynomial function to account for the nonlinearities usually found between biodiversity and productivity<sup>59</sup> (biodiversity + biodiversity<sup>2</sup>), mean atmospheric N deposition, the interaction between biodiversity and N deposition, and the percentages of the area covered by forests and agricultural and urban areas as predictors. MAT and MAP were also included for NBP models, MAT<sub>PV</sub>, MAP<sub>PV</sub> for NBP<sub>PV</sub> models, and MAT<sub>AR1</sub>, MAP<sub>AR1</sub> for NBP<sub>AR1</sub> models. We built an additional model including the interaction between biodiversity and mean and variability in climate (biodiversity:MAT + biodiversity:MAP + biodiversity:MAT<sub>PV</sub> + biodiversity:MAP<sub>PV</sub>) and the interactions MAT:MAP and MAT<sub>PV</sub>:MAP<sub>PV</sub> to further test that the relationship between NBP and biodiversity was not spuriously emerging due to its relationship with climate (Supplementary information – Model summaries, section 8). As indicated by our results, our analyses were able to successfully discern the effect of biodiversity from that of climate (Figure 3). The fact that the relationship between mean NBP and biodiversity emerging from inversions and the TRENDY ensemble (with no parameterisation for biodiversity) were substantially different increases the likelihood of this relationship not being spurious. Models predicting ΔNBP also included the trends in MAT, MAP and the trends in land-use change as predictors. Similarly, models predicting ΔNBP<sub>PV</sub> and ΔNBP<sub>AR1</sub> included trends in temperature and precipitation (ΔMAT and ΔMAP) and their temporal variability (ΔMAT<sub>PV</sub>, ΔMAP<sub>PV</sub>) or autocorrelation (ΔMAT<sub>AR1</sub>, ΔMAP<sub>AR1</sub>) and trends in land-use changes, but not average MAT or MAP or their temporal variability or autocorrelation. Models predicting the combined trends of ΔNBP<sub>PV-AR1</sub> included a second-order biodiversity as a second-order polynomial, N deposition, average and trends in climate, the aggregated trends of MAT<sub>PV-AR1</sub> and MAP<sub>PV-AR1</sub>, average land-use and its temporal trends. We additionally tested whether NBP and ΔNBP were spatially related to increasing ΔNBP<sub>PV-AR1</sub> by fitting two models in which NBP and ΔNBP were the response variables and the aggregated metric of trends in ΔNBP<sub>PV-AR1</sub> (see description above) was their predictor. Models were fitted using normalised variables (mean=0, sd=1), but effect plots were rescaled to their original units to facilitate the interpretation of the results. All analyses were performed using R statistical software v.3.6.3<sup>60</sup>.
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. CO₂ 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. 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-gas-inversion

642 CarboScope: http://www.bgc-jena.mpg.de/CarboScope/

643 Atmospheric CO₂ 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

CO₂ program for providing the records of atmospheric CO₂.

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