# Respiratory loss during late-growing season determines the net carbon dioxide sink in northern permafrost regions

Zhihua Liu, John S. Kimball, Ashley P. Ballantyne, Nicholas C. Parazoo, Wen J. Wang, Ana Bastos, Nima Madani, Susan M. Natali, Jennifer D. Watts, Brendan M. Rogers, et al.

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

## To cite this version:

Zhihua Liu, John S. Kimball, Ashley P. Ballantyne, Nicholas C. Parazoo, Wen J. Wang, et al., Respiratory loss during late-growing season determines the net carbon dioxide sink in northern permafrost regions. Nature Communications, 2022, 13, 10.1038/s41467-022-33293-x. insu-03824353

HAL Id: insu-03824353

https://insu.hal.science/insu-03824353

Submitted on 21 Oct 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Respiratory loss during late-growing season determines the net carbon dioxide sink in northern permafrost regions

Received: 26 January 2022

Accepted: 12 September 2022

Published online: 26 September 2022

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

Check for updates

Zhihua Liu¹,²,¹²✉, John S. Kimball¹,³✉, Ashley P. Ballantyne³,⁴,¹²✉,
Nicholas C. Parazoo⁵, Wen J. Wang⁶✉, Ana Bastos⁷, Nima Madani⁸✉,
Susan M. Natali⁸, Jennifer D. Watts⁸, Brendan M. Rogers⁹, Philippe Ciais⁹,⁴,
Kailiang Yu⁹,⁴, Anna-Maria Virkkala⁹,⁸, Frederic Chevallier⁹,⁴,
Wouter Peters⁹,¹¹, Prabir K. Patra¹⁰ & Naveen Chandra¹⁰

Warming of northern high latitude regions (NHL, > 50 °N) has increased both photosynthesis and respiration which results in considerable uncertainty regarding the net carbon dioxide (CO₂) balance of NHL ecosystems. Using estimates constrained from atmospheric observations from 1980 to 2017, we find that the increasing trends of net CO₂ uptake in the early-growing season are of similar magnitude across the tree cover gradient in the NHL. However, the trend of respiratory CO₂ loss during late-growing season increases significantly with increasing tree cover, offsetting a larger fraction of photosynthetic CO₂ uptake, and thus resulting in a slower rate of increasing annual net CO₂ uptake in areas with higher tree cover, especially in central and southern boreal forest regions. The magnitude of this seasonal compensation effect explains the difference in net CO₂ uptake trends along the NHL vegetation- permafrost gradient. Such seasonal compensation dynamics are not captured by dynamic global vegetation models, which simulate weaker respiration control on carbon exchange during the late-growing season, and thus calls into question projections of increasing net CO₂ uptake as high latitude ecosystems respond to warming climate conditions.

The northern high latitudes (NHL, > 50°N) are experiencing dramatic changes in carbon cycling, evidenced by an increase in the annual terrestrial net CO₂ uptake and in the amplitude of the seasonal cycles of atmospheric CO₂ over the past five decades¹⁻³, but the mechanisms underlying these changes remain highly uncertain. Net CO₂ uptake results from the imbalance between the much larger gross fluxes of plant photosynthesis and ecosystem respiration, which have asynchronous responses to seasonal climatic and environmental change⁴⁻⁶. For example, increased plant photosynthetic CO₂ fixation during the growing season⁷,⁸ may be offset by enhanced respiratory CO₂ release in the fall and/or winter⁹,¹⁰. Such offsets in net CO₂ uptake

among seasons (i.e., seasonal compensation) complicates the detection of the climate-carbon feedbacks at longer time scales over NHL ecosystems⁴,¹¹. Further, the seasonal compensation in CO₂ uptake may vary among biomes given the different sensitivity of above- and belowground carbon cycle processes to climate and environmental controls over changing NHL permafrost conditions¹²⁻¹⁴. Therefore, understanding the magnitude, trends, and spatial patterns of seasonal net CO₂ uptake and its underlying mechanisms is important to address a fundamental question of whether net CO₂ exchange has changed over the NHL, especially in rapidly changing permafrost regions¹⁵⁻¹⁷.
To gain insight into the trends and mechanisms of seasonal CO₂ exchange in the NHL, we addressed the following three major questions. First, what are the trends in net CO₂ uptake and how do they relate to climate, vegetation, and environmental gradients across the NHL? To answer this question, we analyzed an ensemble mean of long-term atmospheric CO₂ inversions (ACIs, *n* = 6, 1980–2017) and a network of Eddy Covariance (EC) observations, each with at least three years of continuous measurements (*n* = 48 sites, 426 site-years, 1990–2017, Fig. S1, Table S1). Second, what are the mechanisms underlying different net CO₂ uptake trends in the NHL? We hypothesized that the different net CO₂ uptake trends were driven by: (H1) different temperature sensitivities of vegetation primary productivity or (H2) the degree to which seasonal net CO₂ uptake is compensated by seasonal respiration losses. We then used structural equation modeling to explore climatic, environmental and vegetation controls on seasonal CO₂ dynamics. Lastly, how well do the latest land surface models replicate the NHL seasonal CO₂ dynamics? To answer this question, we compared the observationally constrained estimates with an ensemble (*n* = 10, 1980–2017) of Dynamic Global Vegetation Models (DGVMs) from the TRENDY intercomparison project (see Methods).

Here, we show that positive trends of annual net CO₂ uptake are smaller with increasing tree cover. We attribute this pattern to the increased magnitude of seasonal compensation due to larger respiratory CO₂ losses during late-growing season at greater levels of tree cover, rather than different temperature sensitivities of vegetation productivity across the NHL. Our synthesis investigation of extensive and diverse empirical datasets shows no evidence that NHL permafrost regions have become net sources of CO₂ to the atmosphere; rather, we find that NHL permafrost regions appear to be gaining C through complex interactions between above and belowground processes that respond differently to seasonal climate. However, current DGVMs are unable to reproduce the different trends of net CO₂ uptake along permafrost-vegetation gradients, potentially calling into question current projections of CO₂ dynamics in the NHL.

## Results

### Trend of net CO₂ uptake is strongly associated with tree cover in the NHL

The ensemble mean of NHL ACIs (inversions) is highly correlated (*r* = 0.78, *p* < 0.001) with independent estimates of net CO₂ exchange at the global scale (Fig. S2)¹⁸, thus giving us confidence to detect regional trends of net CO₂ uptake emerging from the ‘noise’ associated with individual model errors and uncertainties. Estimates from ACIs showed that the NHL contributed roughly 28.7% of the mean global land CO₂ sink (0.67 ± 0.28 PgC yr⁻¹, where positive values denote net C uptake), 17% of the respective global trend (10.4 + 2.02 TgC yr⁻²), and 12.5% of the interannual variability in the global land sink, consistent with another recent model-based assessment¹⁹. From 1980 to 2017, about half (50.3%) of the NHL showed significant (*p* < 0.05) increases in annual net CO₂ uptake, mostly in the permafrost region dominated by tundra shrub and graminoids (Fig. 1a–c, Fig. S1). In contrast, only 4.6% of the NHL showed significant (*p* < 0.05) decreases in annual net CO₂ uptake, mostly in non-permafrost regions dominated by forest (Fig. 1a–c, Fig. S1). The trend of net CO₂ uptake most strongly correlated with tree cover (TC, R² = 0.90, *p* < 0.001), followed by mean annual temperature (R² = 0.77, *p* < 0.001) and permafrost extent (P, R² = 0.27, *p* = 0.033). This suggests that the distribution of forest vegetation, broadly shaped by climate and environmental conditions, is a primary control on net CO₂ uptake trends in the NHL. Despite the fact that net CO₂ uptake in forested regions increased at a slower rate, those regions are still stronger CO₂ sinks than non-forest tundra in colder permafrost regions (Fig. S3). CO₂ source regions are mainly located in the tundra of Alaska and Northeast Siberia, similar to patterns shown by recent analyses¹⁰,²⁰, but the net CO₂ source area is shrinking due to a strong positive trend in net CO₂ uptake (Fig. S3).

To reduce the pixel-level uncertainty in estimating net CO₂ uptake using ACIs, we aggregated net CO₂ uptake to larger regions based on tree cover and permafrost extent, which are negatively correlated in the NHL (R = -0.79, *p* < 0.001, Fig. S4). Regional aggregation of the trends showed that net CO₂ uptake is increasing at a significantly faster rate in tundra (TC < 50%) permafrost (*P* > 10%) regions than in tree-dominated (TC > 50%) non-permafrost regions (*P* < 10%). However, the trends of net CO₂ uptake are not statistically different between low tree cover (TC < 30%) in continuous (*P* > 90%) permafrost regions and intermediate tree cover (TC = 30-50%) in discontinuous (10% < *P* < 90%) permafrost regions (Fig. S4). Therefore, we focused on contrasting trends of net CO₂ uptake between short-vegetated permafrost (TC < 50%) and forest dominated non-permafrost (TC > 50%) regions. The ACI ensemble showed that net CO₂ uptake in permafrost regions has increased at a significantly faster rate (0.58 ± 0.086 gC m⁻² yr⁻²) than in non-permafrost regions (0.13 ± 0.11 gC m⁻² yr⁻²) from 1980 to 2017 (Fig. S4e). This has caused permafrost regions to switch from being CO₂ neutral from 1980 to 2000 (3.16 ± 6.51 gC m⁻² yr⁻¹) to a CO₂ sink after 2000 (15.01 ± 6.14 gC m⁻² yr⁻¹).

We rule out several factors that potentially confound the observed increase in net CO₂ uptake in short-vegetated permafrost regions by: (1) considering the uncertainties of ACI estimates resulting from the variance across individual ACIs, partitioning of fluxes between regions, and time-dependent differences in ACI spread (supplementary text); (2) assessing individual ACIs, where all of the inversions showed increasing net CO₂ uptake in the permafrost region (Fig. S5); (3) conducting a sensitivity test on the time period considered (Fig. S6); and (4) showing that both large-scale patterns and site-level EC measurements had an increasing net CO₂ uptake over permafrost regions, despite occurring at different rates due to scale mismatch between ACIs and EC footprints (Fig. S7); (5) assessing spatial and seasonal consistency of trends from individual ACIs (Fig. S8); and (6) verifying that the addition of more ACIs after 2000 did not significantly alter the trends (Fig. S9 and S6). Therefore, despite generally large uncertainties among ACIs, our results are robust against outliers and in agreement with independent observations from EC data. Therefore, based on the most current atmospheric inversion estimates, the carbon sink strength of shrub and graminoid-dominated permafrost regions has been increasing significantly faster than tree-dominated non-permafrost regions in the NHL.

### Mechanisms underlying different net CO₂ uptake trends along vegetation gradients

Plant photosynthesis is generally thermally-limited and warming-induced relaxation of cold temperature constraints on ecosystem productivity is a key driver for terrestrial CO₂ sink dynamics²¹ and enhanced CO₂ uptake in the NHL²². We therefore tested the first hypothesis (H1) that different trends in net CO₂ uptake along the tree cover gradient were due to differences in vegetation productivity. Physiologically, the larger temperature constraints and amplified rate of warming further north²³ may result in a greater temperature response of vegetation productivity and therefore explain the higher trend of net CO₂ uptake in tundra permafrost regions. Long-term satellite observations of NDVI and GPP derived from light use efficiency models confirm a generally positive trend in annual vegetation productivity along the tree cover gradient (Fig. 2a, c, x-axis), consistent with documented NHL greening trends⁸. We also find a widespread positive correlation between productivity and net CO₂ uptake (Fig. 2b, d), suggesting that increasing net CO₂ uptake is driven by warming-induced increases in productivity. However, the trends of productivity and net CO₂ uptake showed a non-significant (*p* > 0.05) correlation along the tree cover gradient (Fig. 2a, c). In fact, trends in annual productivity increased with tree cover (Fig. S10), and were therefore inconsistent with the trends of net CO₂ uptake shown
<figure><img src="image_5.png" /></figure>

**Fig. 1** | Net CO₂ uptake is increasing faster in regions of lower tree cover and greater permafrost extent from 1980 to 2017 across the NHL (> 50°N) according to an ensemble of ACIs. The maps show regional patterns of a percent tree cover (%TC) based on 35-year (1982-2016) TC layers observed by the Advanced Very High-Resolution Radiometer (AVHRR)¹²; b percent of permafrost extent (%P) based on the average over 1997–2017 from European Space Agency (ESA) Climate

Change Initiative (CCI) Permafrost project data. c The annual trend of net CO₂ uptake (gC m⁻² yr⁻²) from 1980 to 2017 from the ensemble mean of six ACIs in the NHL. Points indicate trends are significant at 0.05 level. **d–f** Correlation between trends of net CO₂ uptake with tree cover, mean annual temperature (MAT), and permafrost extent, respectively.

by ACIs. There are also no apparent trends in annual productivity increase and net CO₂ uptake along the NHL permafrost and temperature gradients (results not shown). Therefore, warming-induced relaxation of temperature constraints appears to be an important driver of increased productivity and net CO₂ uptake across NHL ecosystems, but cannot explain the different trends of net CO₂ uptake seen by ACIs along regional tree cover, temperature, and permafrost gradients.

We then tested the second hypothesis (H2) that seasonal compensation of net CO₂ uptake is the main mechanism explaining the different trends of annual net CO₂ uptake along the tree cover gradients. Across the NHL, accelerating net CO₂ uptake in the early growing season (EGS; May–August) was partially compensated by an accelerating net CO₂ release in the late-growing season (LGS; Sep–Oct), thus resulting in weaker trends in annual net CO₂ uptake (Fig. 3, Fig S11). While the trends of net CO₂ uptake increase in the EGS are similar regardless of tree cover, the trends of net CO₂ release in the LGS increased significantly with greater tree cover further south (Fig. 3d), and thus resulted in a much slower rate of net annual CO₂ uptake and a larger magnitude of seasonal compensation of net CO₂ uptake with greater tree cover (Fig. 1d).

Aggregating over larger regions allows us to more confidently quantify the seasonal compensation of net CO₂ uptake. While the trends of net CO₂ uptake increase were not significantly different

between NHL permafrost and non-permafrost regions, increased emissions in the LGS are significantly lower in permafrost (−1.41 ± 0.15 gC m⁻² yr⁻²) than in non-permafrost (−2.37 ± 0.19 gC m⁻² yr⁻²) regions (Fig. S4e), therefore resulting in a much higher annual rate of net CO₂ uptake increase in permafrost. Consequently, short-vegetated permafrost regions showed a positive sensitivity of annual net CO₂ uptake to mean annual temperature (101 ± 32 TgC yr⁻¹ K⁻¹, *p* < 0.001) whereas tree-dominated non-permafrost regions had a non-significant negative sensitivity (−14 ± 18 TgC yr⁻¹ K⁻¹, *p* = 0.46) in the NHL (Fig. S12). These results indicate that net CO₂ uptake continued to increase in response to warming in permafrost regions, but the positive temperature response has stalled in non-permafrost regions across the NHL due to enhancement in late growing season respiration.

DGVMs have been widely used to explore climate-vegetation-carbon feedbacks and provide recommendations for policy makers¹⁸. Although net biome production (NBP) simulated by DGVMs showed that annual net CO₂ uptake is increasing, there are no apparent model trends of net CO₂ uptake along the tree cover gradients across the NHL (Fig. 4). This is mainly because DGVMs failed to simulate increased CO₂ release in LGS, thus underestimating the seasonal compensation of net CO₂ uptake compared with the observationally-constrained data (Fig. 4), which hinders their ability to capture the changing seasonal atmospheric CO₂ amplitude¹ and its attribution². We try to identify the reasons for this DGVM underestimation below.
<figure><img src="image_6.png" /></figure>

**Fig. 2 |** Weak negative correlations between trends of productivity and net CO₂ uptake suggesting productivity alone cannot explain the different trends of net CO₂ uptake along the tree cover gradient. **a,** **c** show correlation between trends of GIMMS NDVI or LUE GPP and net CO₂ uptake trends along the 5% tree cover interval gradient (blue to green), respectively. Insets show the correlation at

pixel level. **b** and **d** are spatial patterns of correlation (r) between GIMMS NDVI or LUE GPP and net CO₂ uptake, respectively. Points indicate trends are significant at 0.05 level. Insets show the distribution of r between regions with TC > = 50% and TC < 50%.

Climate and environmental controls on seasonal net CO₂ uptake compensation

Our analysis identified the seasonal compensation of net CO₂ uptake as key to understanding the climate-carbon feedback in the NHL. We then used structural equation models (SEMs) to gain insights into the underlying carbon cycle processes and environmental controls contributing to NHL seasonal carbon cycle dynamics. SEM is a multivariate, hypothesis-driven method that is based on a structural model representing a hypothesis about the causal relations among interacting variables. For the EGS (Fig. 5a–c), we considered the relative influences of air temperature (Air T), photosynthetically active radiation (PAR), soil moisture content (SM), percent of tree cover (%TC), and percent of permafrost extent (%P) on gross CO₂ fluxes (i.e., productivity and respiration). All data (EC, ACIs, and TRENDY simulations) showed that net CO₂ uptake was more strongly regulated by ecosystem productivity, which was itself primarily controlled by spring Air T and %TC. This was consistent with previous analyses indicating that temperature-controlled photosynthetic activity and increase in woody vegetation cover were among the major drivers of productivity and net CO₂ uptake in the EGS⁴,²⁴.

In the LGS (Fig. 5d, f), we tested whether net CO₂ release was controlled by moisture stress-induced productivity declines¹¹ or increased respiration²⁵, either from warming-enhanced respiration rate or increased labile carbon carried over from enhanced

productivity early in the growing season (PreGPP)¹⁰,²⁶. Observationally-constrained estimates of net CO₂ fluxes (EC and ACIs) showed that respiration had a stronger control on net CO₂ uptake in the LGS. In contrast, TRENDY simulations showed that productivity continued to have stronger control on net CO₂ uptake in the LGS. All data showed that EGS productivity had the strongest influence on LGS respiration, followed by temperature. %P had a negative influence on LGS respiration, suggesting warmer conditions result in higher LGS respiration. Based on the observationally-constrained data, net CO₂ uptake in the LGS was primarily regulated by respiration, which was mainly controlled by increased labile carbon from enhanced early season productivity and temperature. Therefore, larger increases in productivity as tree cover increases (Fig. S10), together with warming, explain the differences in late growing season respiration along the tree cover gradient (Fig. S13). However, given the strong influence of temperature on EGS productivity (Fig. 5a–c), this indirect effect of warming (via increased input of labile carbon) may also have substantial controls on the increased respiratory CO₂ loss during LGS. Additional analysis also confirmed that warming has both stronger direct and legacy effects on LGS net CO₂ release in the NHL non-permafrost forested regions (Fig. S15), and therefore contributes to a larger seasonal amplitude of net CO₂ fluxes, smaller annual net CO₂ uptake sensitivity to spring temperature, and slower rate of increase in net CO₂ uptake compared to short-vegetated permafrost regions.
<figure><img src="image_7.png" /><figcaption>**Fig. 3 |** Seasonal compensation explains the different trends of net CO₂ uptake along tree cover and permafrost gradients in the NHL. **a–c** show increasing net CO₂ uptake in the early-growing season (a; EGS: May-Aug) and increasing net CO₂ release in the late-growing season (b; LGS: Sep-Oct) and winter (c; Win: Nov-Apr).</figcaption></figure>

d shows the variation in net CO₂ uptake seasonal trends along the 5% tree cover interval gradient in the NHL. Pixels with grey color represent environmental conditions that do not exist.

Discussion

Recent observational²⁷, experimental²⁸,²⁹, and model simulation studies³⁰ have raised concerns that warming-induced acceleration of soil decomposition in permafrost regions may result in a regime shift of arctic ecosystems from a net sink to a net source of CO₂ to the atmosphere³¹,³². In contrast, our analysis suggests that net CO₂ uptake has insofar increased in NHL permafrost regions at a faster rate than in tree dominated non-permafrost regions and that most of these permafrost tundra regions have actually shifted from CO₂ sources to sinks, thus becoming important contributors to the increased CO₂ sink in the northern hemisphere over the past few decades³,³³. Such increases in net CO₂ uptake are primarily driven by warming-enhanced early growing season productivity and increases in shrub/tree vegetation cover³⁴⁻³⁷. Due to the coarse resolution of ACI data, we cannot identify particular ecosystem types contributing to the enhanced net CO₂ uptake, but our results suggest that ecosystems underlaid by intermediate levels of permafrost extent may be a particularly strong contributor to enhanced net C uptake (Fig. 1f) partially due to higher rate of woody cover expansion (Fig. S14). In addition, our analysis indicates northeast Eurasia as a particularly strong contributor to the net CO₂ uptake in the NHL (Fig. 1e), possibly because of higher rates of climate warming³⁸, less water limitation on ecosystem productivity in continuous permafrost areas, regional wetting from enhanced atmospheric moisture transport³⁹,⁴⁰, higher light use efficiency of deciduous larch forests⁴¹, regional expansion of tree and shrub cover⁴²,⁴³, sustained nitrogen deposition⁴⁴, and less severe fire disturbance compared to arctic-boreal North America⁴⁵,⁴⁶. However, as tree migration will likely lag climate change, warming-induced species reassembly and permafrost degradation will significantly alter the above-

vs belowground carbon cycle and its interaction⁴⁷, especially in climate-sensitive permafrost regions with intermediate tree cover⁴⁸.

On the other hand, the net CO₂ uptake in NHL forest dominated regions, such as mixed and evergreen needleleaf boreal forest (Fig. S1), has been increasing at a slower rate, primarily due to enhanced late-season CO₂ release as a result of the warming-induced increase in respiration (Figs. 3 and S13) and emerging moisture limitations on productivity¹¹. Under future warmer and drier conditions, enhanced late-season respiratory CO₂ losses, together with an intensifying fire regime, and moisture stress-induced decreases in productivity, may fully offset or even exceed the net CO₂ gain from photosynthesis, switching the NHL forested region to a predominant net CO₂ source. Our analysis showed that both productivity and respiration are important drivers of the seasonal amplitude of net ecosystem CO₂ uptake in the NHL, but with varied magnitudes and mechanisms across gradients of tree cover, climate, and permafrost. A warming-induced increase of early growing season productivity is the major driver for enhanced net CO₂ uptake in the NHL (Fig. 3), and this mechanism has been supported by model-based analysis² and advancement of the spring zero-crossing date⁴⁹. On the other hand, increasing respiratory CO₂ emission in late growing season of forested ecosystems compensate the CO₂ uptake in early growing season (Fig. 3), and this mechanism is supported by delay of the fall zero-crossing date²⁵,⁴⁹ and airborne and field observations of fall CO₂ fluxes⁹,¹⁰. Moreover, increasing CO₂ respiration during the late growing season with greater tree cover produces a larger seasonal compensation of net CO₂ uptake, indicating boreal forest regions may be a stronger driver for the enhanced atmospheric CO₂ amplitude than arctic tundra. Therefore, our findings help to identify key regions, processes, and seasons to
<figure><img src="image_8.png" /><figcaption>**Fig. 4** | DGVMs cannot capture the decreasing trend of net CO₂ uptake along the NHL tree cover gradient, mainly because the DGVMs failed to simulate increased CO₂ emission in the LGS. a Correlation between trends of net CO₂ uptake with tree cover using the TRENDY ensemble. b–e show net CO₂ uptake in</figcaption></figure>

the early-growing season (EGS: May-Aug), late-growing season (LGS: Sep-Oct), winter (Win: Nov-Apr), and annually along the tree cover and permafrost gradient in the NHL, respectively. Pixels with grey color represent environmental conditions that do not exist.

<figure><img src="image_9.png" /><figcaption>**Fig. 5** | Climate and environmental controls on net CO₂ uptake vary among seasons, and DGVMs need to better capture the respiration processes during the late-growing season. Structural equation models (SEMs) of the relative influence of component fluxes (i.e., gross primary productivity (GPP) and total ecosystem respiration (TER)) on net CO₂ uptake (i.e., net ecosystem exchange: NEE) and climatic and environmental controls including air temperature (AirT), soil moisture (SM), photosynthetically active radiation (PAR), preseason GPP (PreGPP), percent tree cover (%TC), and percent of permafrost extent (%P) on CO₂ fluxes in the NHL. Only significant paths are shown (*p* < 0.05) with the standardized</figcaption></figure>

regression coefficient, and thicker (thinner) lines indicating stronger (weaker) relationships. Numbers in brackets after CO₂ fluxes (i.e., GPP, TER, NEE) indicate the total endogenous (dependent) variation explained by all exogenous (independent) variables in the regression model. Positive (green arrows) and negative (yellow arrows) coefficients indicate respective positive and negative relationships. In the figures, **a**, **b**, **c** are EC, ACI, and TRENDY results, respectively, for the early growing season CO₂ cycle; whereas **d**, **e**, **f** are EC, ACI, and TRENDY results, respectively, for the late growing season CO₂ cycle.
explain the changes in atmospheric CO₂ seasonal amplitude and reconcile previous studies on whether the increasing NHL atmospheric CO₂ seasonal amplitude is due to enhanced productivity² or warming-induced early winter respiration¹⁰.

Our results underscore the importance of the late-growing season CO₂ cycle in regulating seasonal and annual CO₂ dynamics in the NHL. The late-growing season net CO₂ balance has shown contrasting responses to warming²⁵,⁴⁹⁻⁵¹, due to complex interactions among temperature, moisture availability, and radiation on ecosystem productivity¹¹,⁵²,⁵³ and respiration¹⁰,²⁵,²⁶. Our results indicate that the late growing season net CO₂ balance was strongly influenced by respiration, which is sensitive to change of early season productivity and vegetation change induced by warming (Fig. 5) and may increasingly dominate the future carbon balance of NHL ecosystems. However, the current DGVMs underestimated the late-growing season respiration control on net CO₂ loss and therefore failed to capture the full magnitude of seasonal CO₂ dynamics in the high latitudes (Figs. 4 and 5). This helps explain why DGVM-based assessments tend to underestimate the seasonal atmosphere CO₂ amplitude in the NHL¹,⁵⁴,⁵⁵ and attribute long-term changes largely to boreal productivity². A satellite-driven carbon-flux model (SMAP L4C) that accounted for both atmosphere and soil moisture related controls on net CO₂ exchange showed similar behavior as the DGVMs (Fig. S15). To better simulate late-growing season respiration and its impact on the net CO₂ balance, these models may need to include additional processes, such as direct and lagged effects among temperature, moisture availability, and productivity¹¹, and the resultant influences on decomposition and turnover of soil organic carbon⁵⁶,⁵⁷.

In summary, there is currently no evidence that permafrost-dominated regions are, as a whole, a net source of CO₂ to the atmosphere. On the contrary, net CO₂ uptake over tundra dominated permafrost regions has been increasing at a faster rate than tree dominated non-permafrost regions in the NHL since the 1980s. As a result, the NHL permafrost region is becoming a strong contributor to the northern terrestrial carbon sink. The faster rate of CO₂ sink increase in NHL regions is mainly due to enhanced CO₂ uptake in the early growing season, which is outpacing carbon emissions in the late growing season. On the other hand, the larger seasonal compensation in net CO₂ uptake as tree cover increases is an important driver of the increasing seasonal atmospheric CO₂ amplitude in the NHL. Current DGVMs need to better simulate respiration processes in response to changing landscape conditions in order to capture the full magnitude of seasonal CO₂ dynamics and generate more accurate predictions of climate-carbon feedbacks in northern high latitudes. Although we do not see any evidence yet of positive carbon-climate feedbacks in high latitude ecosystems, thresholds in CO₂ dynamics may occur as these ecosystems transition in response to the changing state of permafrost.

## Methods

We focused on the Northern High Latitudes (NHL, latitude > 50°N, excluding Greenland) due to their importance for carbon (CO₂-C, the same hereafter)-climate feedbacks in the Earth system. To minimize the potential human influence on the CO₂ cycle, we excluded areas under agricultural management (croplands, cropland/natural vegetation mosaic, and urban types), and considered only pixels of natural vegetation defined from the MODIS MCD12Q1 (v006) based IGBP land cover classification. Our main focus was the NHL permafrost region because permafrost plays a critical role in the ecology, environment, and society in the NHL. Permafrost, or permanently frozen ground, is defined as ground (soil, sediment, or rock) that remains at or below 0 °C for at least two consecutive years. The occurrence of permafrost is primarily controlled by temperature and has a strong effect on hydrology, soils, and vegetation composition and structure. Based on the categorical permafrost map from the International Permafrost

Association⁵⁸, the permafrost region (excluding permanent snow/ice and barren land), including sporadic (10–50%), discontinuous (50–90%), and continuous (>90%) permafrost, encompasses about 15.7 × 10⁶ km², accounts for 57% of the NHL study dominion, and is dominated by tundra (shrubland and grass) and deciduous needleleaf (i.e., larch) forest that is regionally abundant in Siberia. The NHL non-permafrost region covers about 11.9 × 10⁶ km² and is dominated by mixed and evergreen needleleaf boreal forests (Fig. S1).

## Atmospheric CO₂ inversions (ACIs)

ACIs provide regionally-integrated estimates of surface-to-atmosphere net ecosystem CO₂ exchange (NEE<sub>ACI</sub>) fluxes by utilizing atmospheric CO₂ concentration measurements and atmospheric transport models⁵⁹. ACIs differ from each other mainly in their underlying atmospheric observations, transport models, spatial and temporal flux resolutions, land surface models used to predict prior fluxes, observation uncertainty and prior error assignment, and inversion methods. We used an ensemble mean of six different ACI products, each providing monthly gridded NEE<sub>ACI</sub> at 1-degree spatial resolution, including Carbon-Tracker 2019B (2000-2019, CT2019)⁶⁰, Carbon-Tracker Europe 2020 (2000-2019, CTE2020)⁶¹, Copernicus Atmosphere Monitoring Service (1979-2019, CAMS)⁶², Jena CarboScope (versions s76_v4.2 1976-2017, and s85_v4.2 1985-2017)⁶³,⁶⁴, and JAMSTEC (1996-2017)⁶⁵. The monthly gridded ensemble mean NEE<sub>ACI</sub> at 1-degree spatial resolution was calculated using the available ACIs from 1980-2017. Monthly ACI ensemble mean NEE<sub>ACI</sub> data were summed to seasonal and annual values, and used to calculate the spatial and temporal trends of net CO₂ uptake, and to investigate its relationship to climate and environmental controls.

## Productivity dataset

Direct observations of vegetation productivity do not exist at a circumpolar scale. We therefore used two long-term gridded satellite-based estimates of vegetation productivity, including gross primary production (GPP) derived using a light use efficiency (LUE) approach (LUE GPP, 1982–1985)²¹,⁶⁶ and satellite observations of Normalized Difference Vegetation Index (NDVI) from the Global Inventory Modeling and Mapping Studies (GIMMS NDVI, 1982–1985)⁶⁷. LUE GPP (monthly, 0.5° spatial resolution, 1982–2015) is calculated from satellite observations of NDVI from the Advanced Very High-Resolution Radiometer (AVHRR; 1982 to 2015) combined with meteorological data, using the MOD17 LUE approach. LUE GPP has been extensively validated with a global array of eddy-flux tower sites⁶⁸⁻⁷⁰ and tends to provide better estimates in ecosystems with greater seasonal variability at high latitudes. Following⁶⁶,⁷¹, we used the ensemble mean of GPP estimates from three of the most commonly used meteorological data sets: National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) reanalysis; NASA Global Modeling and Assimilation Office (GMAO) Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2); and European Center for Medium-Range Weather Forecasting (ECMWF). GIMMS NDVI (bimonthly, 1/12 spatial resolution, 1982–2015) provides the longest satellite observations of vegetation “greenness”, and is widely used in studies of phenology, productivity, biomass, and disturbance monitoring as it has proven to be an effective surrogate of vegetation photosynthetic activity⁷².

The gridded GPP data were resampled to 1-degree resolution at monthly time scales, to be consistent with NEE<sub>ACI</sub>, and used to test (H1) whether greater temperature sensitivity of vegetation productivity explains the different trends in net CO₂ uptake across the NHL. LUE GPP was also used to calculate monthly total ecosystem respiration (TER) as the difference between GPP and NEE<sub>ACI</sub> (i.e., TER<sub>residual</sub> = GPP-NEE<sub>ACI</sub>) from 1982-2015, as global observations of respiration do not exist. The NEE<sub>ACI</sub>, GPP and TER<sub>residual</sub> were used as observation-constrained top-down CO₂ fluxes to investigate mechanisms
underlying the seasonal CO₂ dynamics in the structural equation modeling and additional decision tree-based analysis.

## Eddy Covariance (EC) measurements of bottom-up CO₂ fluxes

A total of 48 sites with at least three years of data representing the major NHL ecosystems were obtained from the FLUXNET2015 database (Table S1 and Fig. S1). EC measurements provide direct observations of net ecosystem CO₂ exchange (NEE) and estimate the GPP and TER flux components of NEE using other climate variables. Daily GPP and TER were estimated as the mean value from both the nighttime partitioning method<sup>73</sup> and the light response curve method<sup>74</sup>. More details on the flux partitioning and gap-filling methods used are provided by<sup>75</sup>. Daily fluxes were summed into seasonal and annual values and used to compare with trends from ACIs (Fig. S7), to estimate the climate and environmental controls on the CO₂ cycle in the pathway analysis (Fig. 5), and to calculate the net CO₂ uptake sensitivity to spring temperature (Fig. S14).

## Ensemble of dynamic global vegetation models (TRENDY simulations)

The TRENDY intercomparison project compiles simulations from state-of-the-art dynamic global vegetation models (DGVMs) to evaluate terrestrial energy, water, and net CO₂ exchanges<sup>76</sup>. The DGVMs provide a bottom-up approach to evaluate terrestrial CO₂ fluxes (e.g., net biome production [NBP]) and allow deeper insight into the mechanisms driving changes in carbon stocks and fluxes. We used monthly NBP, GPP, and TER (autotrophic + heterotrophic respiration; Ra + Rh) from ten TRENDY v7 DGVMs<sup>76</sup>, including CABLE-POP, CLM5.0, OCN, ORCHIDEE, ORCHIDEE-CNP, VISIT, DLEM, LPJ, LPJ-GUESS, and LPX. We analyzed the “S3” simulations that include time-varying atmospheric CO₂ concentrations, climate, and land use. All simulations were based on climate forcing from the CRU-NCEPv4 climate variables at 6-hour resolution. CO₂ flux outputs were summarized monthly at 1-degree spatial resolution from 1980 to 2017. Monthly ensemble mean NBP, GPP, and TER were summed to seasonal and annual values, and then used to compare with observation-constrained ACI top-down CO₂ fluxes (Figs. 4 and 5).

## Satellite data-driven carbon flux estimates (SMAP L4C)

We also used a much finer spatio-temporal simulation of carbon fluxes from the NASA Soil Moisture Active Passive (SMAP) mission Level 4 Carbon product (L4C) to quantify the temperature and moisture sensitivity of NHL CO₂ exchange<sup>77</sup>. The SMAP L4C provides global operational daily estimates of NEE and component CO₂ fluxes for GPP and TER at 9 km resolution since 2015; whereas, an offline version of the L4C model provides a similar Nature Run (NR) carbon flux record over a longer period (2000-present), but without the influence of SMAP observational inputs. The L4C model has been calibrated against FLUXNET tower CO₂ flux measurements and shows favorable performance and accuracy in high latitude regions<sup>4,77</sup>. In this analysis, daily gridded CO₂ fluxes at 9-km resolution from the L4C NR record were summed to seasonal and annual values, and used to calculate the sensitivity of net C uptake in response to spring temperature (Fig. S14).

CO₂ fluxes in this analysis are defined with respect to the biosphere so that a positive value indicates the biosphere is a net sink of CO₂ absorbed from the atmosphere. The different data products described above use different terminology (e.g., NEE, NBP) with slightly different meanings; however, they all provide estimates of net land-atmosphere CO₂ exchange<sup>78</sup>.

## Climate, tree cover, permafrost, and soil moisture data

Monthly gridded air temperatures at 0.5-degree spatial resolution from 1980 to 2017 were obtained from the Climate Research Unit (CRU

TS v4.02) at the University of East Anglia<sup>79</sup>. Air temperature was summarized at seasonal and annual scales to calculate temperature sensitivities of net CO₂ uptake and to investigate the mechanism underlying the seasonal CO₂ dynamics.

Percent tree cover (%TC) at 0.05-degree spatial resolution was averaged over a 35-year (1982-2016) period using annual %TC layers derived from the Advanced Very High-Resolution Radiometer (AVHRR) (Fig. 1a)<sup>42</sup>. %TC was binned using 5% TC intervals to assess its relation to net CO₂ uptake, or aggregated at a regional scale (e.g., TC > 50% or TC < 50%) to contrast variation of net CO₂ uptake or used to explore the mechanism underlying the seasonal CO₂ dynamics.

We used two permafrost maps, including a continuous permafrost extent map from European Space Agency's (ESA) Climate Change Initiative (CCI) Permafrost project (Permafrost_CCI, 1-km spatial resolution)<sup>80</sup> and a categorical permafrost zone map from IPA (International Permafrost Association, Permafrost_IPA, vector data)<sup>58</sup>. The Permafrost_CCI product was derived from a thermal model driven and constrained by satellite observations of land surface temperature. The percent of permafrost extent (%P) for each year was calculated as the yearly fraction of permafrost-underlain (ground temperatures at 2 m depth < 0) and permafrost-free area within a pixel. The %P was averaged over 1997-2017 and aggregated to 0.05-degree spatial resolution (Fig. 1b). The %P metric corresponds well with Permafrost_IPA permafrost zones, which distinguish isolated (0-10%), sporadic (10-50%), discontinuous (50-90%) and continuous permafrost (90-100%). In this analysis, the %P map was binned at 5% intervals to show its relation to net CO₂ uptake or used to explore the mechanism underlying the seasonal CO₂ dynamics. The permafrost map was also aggregated into continuous permafrost (ConP, P > 90%), discontinuous permafrost (DisconP, 10% < P < 90%), and non-permafrost (NoP, P < 10%) zones for regional analysis.

We used the ESA CCI soil moisture product (SM, v4.5) produced from combined satellite active and passive microwave remote sensing observations<sup>81</sup> at daily, 0.25-degree spatial resolution from 1980 to 2017. The ESA CCI SM product was developed using data derived from C-band scatterometers suitable for SM retrieval, such as European Remote Sensing Satellites (ERS-1/2) and METOP, as well as the use of data from multi-frequency microwave radiometers such as the Scanning Multichannel Microwave Radiometer (SMMR), Special Sensor Microwave Imager (SSM/I), Microwave Imager (TMI), Advanced Microwave Scanning Radiometer (AMSR-E), and Windsat<sup>81</sup>.

The ESA CCI SM product characterizes surface (0-5 cm depth) SM conditions that are highly variable at daily time scales. We first aggregated the daily gridded SM estimates into monthly mean gridded SM to reduce the temporal variability and spatial gaps. Monthly gridded SM estimates were then aggregated into seasonal and annual averages and used to explore their relationships with the land-atmosphere CO₂ fluxes.

## Net CO₂ trends along vegetation, climate, and permafrost gradients

We binned %TC and %P into 5% intervals, and annual mean air temperature into 1-degree intervals. The net CO₂ uptake for the early-growing season (EGS: May-August), late-growing season (LGS: September-October), winter (November-April), and annual periods from 1980 to 2017 was summarized using the ensemble mean of the ACIs (NEE<sub>ACI</sub>) at each binned interval. The seasonal and annual mean net CO₂ uptake for each interval of %TC, %P, and temperature was then regressed against years using linear regression. The slope of the regression was interpreted as the net CO₂ uptake trend (gC m⁻² yr⁻²). Finally, trends of net CO₂ uptake at seasonal and annual scales were plotted against %TC, %P, and air temperature to understand the trend and seasonality of net CO₂ uptake along the vegetation, climate, and permafrost gradient (Figs. 1 and 3).
## Regional analysis

To reduce the pixel-scale uncertainties of net CO₂ uptake using ACIs data, we also calculated the trends of net CO₂ uptake at regional scales, which were classified by %TC and %P (Fig. S4). Using %TC, the NHL was divided into low (<= 30%), intermediate (30–50%), and high (>50%) tree cover regions. Using %P, the NHL was divided into continuous (ConP, P > 90%), discontinuous (DisconP, 10% < P < 90%), and non-permafrost (NoP, P < 10%) regions. %TC and %P are highly correlated, such that short-vegetated regions (TC < 50%) are primarily overlying permafrost (ConP and DisconP), whereas the tree-dominated regions (TC > 50%) are primarily in non-permafrost (Fig. S4). Spatial and temporal patterns for the net CO₂ uptake trend were calculated seasonally and annually from 1980 to 2017 using the ensemble mean of the ACIs (NEE<sub>ACI</sub>) over different NHL regions. Seasonal and annual mean net CO₂ uptake for each region was regressed against years using linear regression. The slope of the regression was interpreted as the net CO₂ uptake trend (gC m⁻² yr⁻², Fig. S4).

## Robustness analysis

Since trends of net CO₂ uptake are not statistically different between low tree cover (< 30% tree cover in ConP regions) and intermediate tree cover (30-50% tree cover in DisconP regions) (Fig. S4), we aggregated these two regions into a short-vegetated (TC < 50%) permafrost region. We contrasted the net CO₂ uptake between short-vegetated (TC < 50%) permafrost and tree-dominated (TC > 50%) non-permafrost regions. To confirm the net CO₂ uptake trend between permafrost and non-permafrost regions obtained by the ACIs, we performed a series of additional analyses as follows.

**Individual ACI analysis.** We repeated the trend analysis on permafrost and non-permafrost regions for each individual ACI. All ACIs showed that net CO₂ uptake is increasing faster in the permafrost region than in the non-permafrost region, although only three ACIs (CT2019B, CTE2020, CAMS) showed significantly faster trends over the available data periods (Fig. S5).

**Random years and length analysis.** We randomly selected starting years and length of years (> = 10 years) for the ensemble mean ACIs, and then repeated the trend analysis. Results consistently showed that net CO₂ uptake is increasing at a much faster rate in permafrost regions than in non-permafrost regions (Fig. S6).

**Site-level analysis.** We conducted a site-level comparison between EC observations and the ACI ensemble at NHL tower site locations since the 1990s. Both the EC measurements and ACI ensembles confirmed that the net CO₂ uptake is increasing with decreasing tree cover over and increasing permafrost extent, such that the short-vegetated permafrost region (%TC < 50%) had a higher rate of net CO₂ uptake increase than the tree-dominated non-permafrost region (%TC > 50%) during the 1990s – 2010s (Fig. S7). This result is also consistent with an independent EC-based analysis showing that forest sties had a slower increasing rate of net CO₂ uptake than other ecosystem types (e.g., tundra, wetland) in northern permafrost regions¹³. Expectedly, the net CO₂ uptake trends calculated by the EC observations are much higher than the ACIs because (1) the ACIs represent a much larger spatial footprint than the ecosystem-level EC measurements, and therefore average out the local-scale variability, and (2) some episodic ecosystem processes, such as fire disturbance, were not accounted for by the EC observations.

**ACI uncertainty analysis.** The uncertainties among different ACIs may also affect trend estimates. Uncertainty in ACIs estimates may be due to (i) spread across inversions, (ii) differences among inversions in partitioning of fluxes between permafrost and non-permafrost regions, and (iii) time-dependent differences in inversion spread.

Therefore, we used a Generalized Linear Mixed effects Model (GLMM) to estimate trends, considering (i)-(iii) as random effects. The GLMM showed that even after accounting for the uncertainty due to inversion spread, the rate of net CO₂ uptake in permafrost regions is still significantly faster than non-permafrost regions (see supplementary text for full details).

**Spatial trend agreement analysis.** to assess the spatial consistency of trends of net CO₂ uptake derived from individual ACIs, we calculated the number of ACIs that showed similar trends (either decreasing or increasing). If the majority of ACIs (i.e., >= 5 of 6 ACIs) showed similar trends, we considered the trend in these areas to be true. We reported the trends within these areas here. Generally, we found more areas in the EGS showing similar trends than the other seasons, reflecting larger uncertainties in understanding the carbon cycle in LGS and winter season. In the EGS, 76% of the non-permafrost region and 77% of the permafrost region showed increasing net CO₂ uptake. Annually, only 46% of the non-permafrost region and 51% of the permafrost region showed increasing net CO₂ uptake because of high uncertainties in the LGS and winter (Fig. S8).

**Time-series analysis.** to see if the trend analysis was affected by including more ACIs in the record after year 2000, we compared the trends calculated from the two longest ACIs (i.e., CAMS, and Jena CarboScope (s76)) with trends calculated from the 3 shortest ACIs (i.e., CT2020B, CTE2020, and JAMSTEC). Results showed that permafrost regions still have a significantly higher rate of net CO₂ uptake using the two longest continuous ACIs since 1980 (CAMS and Jena CarboScope (s76)). The trend of net CO₂ uptake is significantly higher after 2000, which may elevate trends in the permafrost region since the 1980s (Fig. S9). Therefore, including more ACIs after 2000 may change the magnitude of trends in permafrost regions, but does not alter our finding that net CO₂ uptake increased faster in permafrost than non-permafrost regions.

## Comparison with Global Carbon Budget 2020 (GCB2020)

The global land CO₂ sink and trend between the ensemble of ACIs and the Global Carbon Budget 2020 (GCB2020), which was estimated from the multi-model mean of 16 DGVMs¹⁸, were also compared from 1980 to 2017 (Fig. S2). We note that this comparison is not a fully independent representation of the land sink, as the ACIs and GCB2020 land-sink estimates used similar fossil fuel emissions, ocean fluxes, and land-use change emissions, and some TRENDY model calculations served as a-priori land-sink for individual ACIs.

## Comparison with TRENDY NBP simulations

To see if current state-of-art land surface simulation models reproduce the trend and seasonality of net C uptake relative to the ACIs, the same analysis of net CO₂ uptake spatial and temporal patterns and trends was applied to the NBP simulations from TRENDY (Fig. 4).

## Trend and seasonality of productivity

To test H1 (whether greater temperature sensitivity of vegetation productivity in the NHL explains the different trends in net CO₂ uptake along the gradient of tree cover), we first calculated the trend and seasonality of productivity using the two productivity products described above, following the same trend analysis procedure used for the ACIs (Fig. S10). The analysis showed that the trend and seasonality are similar using different productivity proxies despite different algorithms, observations, and driving variables among products. Both annual NDVI and GPP data showed positive trends in the NHL (Fig. S10), consistent with greening trends reported from other studies. Furthermore, trends in productivity generally increased with tree cover (Fig. S10), in direct contrast to the decreasing trends of net CO₂ uptake with increase in tree cover.
To further test HI, we calculated the correlations between the trends in productivity and net CO₂ uptake along the tree cover gradient at both pixel level and 5% intervals (Fig. 2a, c). We also calculated the correlation between pixel-level time-series net CO₂ uptake and productivity from 1982 to 2015 (Fig. 2b, d).

## Path analysis to explore the mechanism underlying the seasonal CO₂ cycles

Structural equation modeling (SEM) is a powerful multivariate technique to evaluate direct and indirect effects between pre-assumed causal relationships within multivariate data⁸². SEM combines two statistical methods, confirmatory factor analysis and path analysis, and aims to find the causal relationship among variables by creating a path diagram. SEM is well suited to test the complex interaction among CO₂ cycles and controls.

Based on a priori expectations, we constructed one structural equation model (SEM) for each season (EGS and LGS) to test the relative influence of component CO₂ fluxes (i.e., productivity [GPP] and respiration [TER]) on the net CO₂ flux (NEE), and the climate and environmental controls on the NEE component CO₂ fluxes. Our goal was not to precisely predict the spatial and temporal CO₂ cycle variability, but rather to illuminate and quantify the relative influence of the major climatic and environmental controls on CO₂ fluxes. We constructed different models to explain EGS and LGS CO₂ cycles because of the varied climate and environmental controls across seasons over the NHL⁴. For EGS, we tested the relative influences of air temperature (Air T), photosynthetic active radiation (PAR), soil moisture (SM), percent tree cover (%TC), and percent permafrost extent (%P) on vegetation photosynthetic activity and respiration, and the resultant influence on net CO₂ uptake⁴,²⁴. For LGS, we tested the relative influences of Air T, PAR, SM, %TC, %P, and labile carbon carried over from the early growing season (PreGPP) on vegetation photosynthetic activity and respiration, and the resultant influence on the net CO₂ uptake results⁴,¹⁰,²⁶,⁵² (Fig. 5).

SEMs were fit using the sem function of the lavaan package in R⁸³. The performances of the SEMs were evaluated using a combination of the chi-square statistic (where $\chi^2 \le 2$ and $p > 0.05$ indicate a good fitting model), Bentler's comparative fit index (CFI, where $CFI \approx 1$ indicates a good fitting model), and the root mean square error of approximation (RMSEA; where $RMSEA \le 0.05$ and $p > 0.1$ indicate a good fitting model). The standardized regression coefficient can be interpreted as the relative influences of exogenous (independent) variables. The $R^2$ indicates the total variation in an endogenous (dependent) variable explained by all exogenous (independent) variables.

## Direct and legacy effects of temperature on seasonal net CO₂ uptake

Because landscape thawing and snow conditions regulate the onset of vegetation growth and influence the seasonal and annual CO₂ cycles in the NHL²⁴,⁸⁴, we also analyzed the legacy effects of spring (May–Jun) temperature on seasonal net CO₂ uptake. We regressed seasonal and annual net CO₂ uptake from the site-level EC observations, regional-level ACI ensemble, and the TRENDY NBP ensemble against spring (May–June) air temperature. For EC observations, net CO₂ uptake (i.e., NEE) and air temperature were summarized from site-level measurements. For the ACIs and TRENDY ensemble, net CO₂ uptake (i.e., NEE<sub>ACI</sub> and NBP) was summarized as regional means from the ACIs and TRENDY ensemble outputs, and air temperature was summarized as regional means from CRU temperature. The slope of the regression line was interpreted as the spring temperature sensitivity of the CO₂ cycle. Simple linear regression was used here mainly due to the strong influence of spring temperature on the seasonal and annual CO₂ cycle in NHL ecosystems³⁰. Temperature sensitivity ($\gamma$: g C m⁻² day⁻¹ K⁻¹) is the change in net CO₂ flux (g C m⁻² day⁻¹) in response to a 1-degree

temperature change. The sensitivity of net CO₂ uptake to warm spring anomalies was calculated for different seasons (EGS, LGS, and annual) and regions (i.e., permafrost and non-permafrost), and the T-test was used to test for the difference in $\gamma$ among different regions, seasons, and datasets. Similarly, direct effects of temperature on net CO₂ uptake were calculated using the same season data (Fig. S14).

Observationally-constrained estimates (EC and ACIs) showed that the sensitivity of net CO₂ uptake in the EGS to spring temperature is positive ($\gamma > 0$) and not statistically different ($p > 0.05$) between permafrost and non-permafrost regions ($\gamma_{ACI}^{np} = 0.125 \pm 0.020 \text{ gC m}^{-2} \text{ day}^{-1} \text{ K}^{-1}$; $\gamma_{EC}^{np} = 0.052 \pm 0.013 \text{ gC m}^{-2} \text{ day}^{-1} \text{ K}^{-1}$). In contrast, the sensitivity of net CO₂ uptake in LGS to spring temperature is negative ($\gamma < 0$) and significantly smaller in magnitude in permafrost regions ($\gamma_{ACI}^{p} = -0.054 \pm 0.0064 \text{ gC m}^{-2} \text{ day}^{-1} \text{ K}^{-1$; $\gamma_{EC}^{p} = -0.0157 \pm 0.0084 \text{ gC m}^{-2} \text{ day}^{-1} \text{ K}^{-1}$) than non-permafrost regions ($\gamma_{ACI}^{np} = -0.093 \pm 0.0070 \text{ gC m}^{-2} \text{ day}^{-1} \text{ K}^{-1}$; $\gamma_{EC}^{np} = -0.061 \pm 0.013 \text{ gC m}^{-2} \text{ day}^{-1} \text{ K}^{-1}$). Net CO₂ uptake in the LGS is also negatively correlated to LGS temperature in the non-permafrost region ($\gamma_{ACI}^{np} = -0.056 \pm 0.0062 \text{ gC m}^{-2} \text{ day}^{-1} \text{ K}^{-1}$; $\gamma_{EC}^{np} = -0.061 \pm 0.011 \text{ gC m}^{-2} \text{ day}^{-1} \text{ K}$), but showed a varied response in permafrost regions ($\gamma_{ACI}^{p} = -0.027 \pm 0.0046 \text{ gC m}^{-2} \text{ day}^{-1} \text{ K}^{-1$; $\gamma_{EC}^{p} = 0.024 \pm 0.011 \text{ gC m}^{-2} \text{ day}^{-1} \text{ K}^{-1}$, Fig. S14). Overall, spring warming has stronger legacy effects on LGS net CO₂ release in the NHL non-permafrost regions, and therefore contributes to a larger seasonal amplitude of net CO₂ fluxes, smaller annual net CO₂ uptake sensitivity to spring temperature, and slower rate of increase in net CO₂ uptake compared to permafrost regions.

Consistent with regional estimates from ACIs, NBP simulated by DGVMs showed a positive sensitivity of EGS net CO₂ uptake to spring temperature ($\gamma_{TRENDY}^{np} = 0.044 \pm 0.025 \text{ gC m}^{-2} \text{ day}^{-1} \text{ K}^{-1$; $\gamma_{TRENDY}^{p} = 0.111 \pm 0.0075 \text{ gC m}^{-2} \text{ day}^{-1} \text{ K}^{-1}$), contributing to an annual net CO₂ uptake in the NHL. However, simulated LGS NBP showed non-significant ($\gamma_{TRENDY}^{p} = 0.0034 \pm 0.0036 \text{ gC m}^{-2} \text{ day}^{-1} \text{ K}^{-1}$) or positive ($\gamma_{TRENDY}^{np} = 0.0201 \pm 0.0061 \text{ gC m}^{-2} \text{ day}^{-1} \text{ K}^{-1}$) sensitivities to spring warming, failed to capture the seasonal compensation in net CO₂ uptake in both permafrost and non-permafrost regions, which therefore partially explained the lack of trends in NBP along the NHL tree cover gradient (Fig. 4).

## Data availability

All data analyzed in this study are publicly available. Carbon tracker is obtained from NOAA Earth System Research Laboratory (https://www.esrl.noaa.gov/gmd/ccgg/carbontracker/), Carbon Tracker Europe from Wageningen University (http://www.carbontracker.eu/), Jena CarboScope is from MPG (http://www.bgc-jena.mpg.de/CarboScope/), and CAMS from ECMWF (http://apps.ecmwf.int/datasets/data/cams-ghg-inversions/). SMAP L4C data are from the National Snow and Ice Data Center (https://nsidc.org/). TRENDY simulation is obtained from http://dgvm.ceh.ac.uk/index.html. Climate data are from Climate Research Unit (https:// rudata.uea.ac.uk/cru/data/hrg/). EC data is from FLUXNET2015 Dataset (https://fluxnet.org/data/fluxnet2015-dataset/).

## Code availability

All computer codes used in this study are available via GitHub at https://github.com/liuzh811/PermafrostCarbonCycle, and also archived at Zenodo (https://zenodo.org/badge/latesttdoi/510523596).

## References

1. Graven, H. D. et al. Enhanced seasonal exchange of CO₂ by northern ecosystems since 1960. *Science* **341**, 1085–1089, https://doi.org/10.1126/science.1239207 (2013).

2. Forkel, M. et al. Enhanced seasonal CO₂ exchange caused by amplified plant productivity in northern ecosystems. *Science* **351**, 696–699, https://doi.org/10.1126/science.aac4971 (2016).
3. Ciais, P. et al. Five decades of northern land carbon uptake revealed by the interhemispheric CO₂ gradient. *Nature* **568**, 221–225, https://doi.org/10.1038/s41586-019-1078-6 (2019).

4. Liu, Z. et al. Increased high-latitude photosynthetic carbon gain offset by respiration carbon loss during an anomalous warm winter to spring transition. *Glob. Chang Biol.* **26**, 682–696, https://doi.org/10.1111/gcb.14863 (2020).

5. Xia, J. Y. et al. Terrestrial carbon cycle affected by non-uniform climate warming. *Nat. Geosci.* **7**, 173–180, https://doi.org/10.1038/Ngeo2093 (2014).

6. Helbig, M. et al. Warming response of peatland CO₂ sink is sensitive to seasonality in warming trends. *Nat. Clim. Change* **12**, 743–749, https://doi.org/10.1038/s41558-022-01428-z (2022).

7. Myers-Smith, I. H. et al. Complexity revealed in the greening of the Arctic. *Nat. Clim. Change* **10**, 106–117, https://doi.org/10.1038/s41558-019-0688-1 (2020).

8. Piao, S. et al. Characteristics, drivers and feedbacks of global greening. *Nat. Rev. Earth Environ.* **1**, 14–27, https://doi.org/10.1038/s43017-019-0001-x (2019).

9. Natali, S. M. et al. Large loss of CO₂ in winter observed across the northern permafrost region. *Nat. Clim. Change* **9**, 852–857, https://doi.org/10.1038/s41558-019-0592-8 (2019).

10. Commane, R. et al. Carbon dioxide sources from Alaska driven by increasing early winter respiration from Arctic tundra. *Proc. Natl Acad. Sci. USA* **114**, 5361–5366, https://doi.org/10.1073/pnas.1618567114 (2017).

11. Buermann, W. et al. Widespread seasonal compensation effects of spring warming on northern plant productivity. *Nature* **562**, 110–114, https://doi.org/10.1038/s41586-018-0555-7 (2018).

12. Mekonnen, Z. A. et al. Arctic tundra shrubification: a review of mechanisms and impacts on ecosystem carbon balance. *Environ. Res. Lett.* **16**, 053001 (2021).

13. Li, Z.-L. et al. Changes in net ecosystem exchange of CO₂ in Arctic and their relationships with climate change during 2002–2017. *Adv. Clim. Change Res.* **12**, 475–481 (2021).

14. Welp, L. R. et al. Increasing summer net CO₂ uptake in high northern ecosystems inferred from atmospheric inversions and comparisons to remote-sensing NDVI. *Atmos. Chem. Phys.* **16**, 9047–9066, https://doi.org/10.5194/acp-16-9047-2016 (2016).

15. Schuur, E. A. et al. Climate change and the permafrost carbon feedback. *Nature* **520**, 171–179, https://doi.org/10.1038/nature14338 (2015).

16. McGuire, A. D. et al. Dependence of the evolution of carbon dynamics in the northern permafrost region on the trajectory of climate change. *Proc. Natl Acad. Sci. USA* **115**, 3882–3887, https://doi.org/10.1073/pnas.1719903115 (2018).

17. McGuire, A. D. et al. An assessment of the carbon balance of Arctic tundra: comparisons among observations, process models, and atmospheric inversions. *Biogeosciences* **9**, 3185–3204, https://doi.org/10.5194/bg-9-3185-2012 (2012).

18. Friedlingstein, P. et al. Global Carbon Budget 2020. *Earth Syst. Sci. Data* **12**, 3269–3340, https://doi.org/10.5194/essd-12-3269-2020 (2020).

19. Tagesson, T. et al. Recent divergence in the contributions of tropical and boreal forests to the terrestrial carbon sink. *Nat. Ecol. Evol.* **4**, 202–209, https://doi.org/10.1038/s41559-019-1090-0 (2020).

20. Virkkala, A.-M. et al. Statistical upscaling of ecosystem CO₂ fluxes across the terrestrial tundra and boreal domain: regional patterns and uncertainties. *Glob. Change Biol.* **27**, 4040–4059, https://doi.org/10.1111/gcb.15659 (2021).

21. Running, S. W. et al. A continuous satellite-derived measure of global terrestrial primary production. *Bioscience* **54**, 547–560 (2004).

22. Liu, J., Wennberg, P. O., Parazoo, N. C., Yin, Y. & Frankenberg, C. Observational constraints on the response of high-latitude northern forests to warming. *AGU Adv.* **1**, e2020AV000228 (2020).

23. Serreze, M. C. & Barry, R. G. Processes and impacts of Arctic amplification: a research synthesis. *Glob. Planet. Change* **77**, 85–96 (2011).

24. Parazoo, N. C. et al. Spring photosynthetic onset and net CO₂ uptake in Alaska triggered by landscape thawing. *Glob. Chang Biol.* **24**, 3416–3435, https://doi.org/10.1111/gcb.14283 (2018).

25. Piao, S. et al. Net carbon dioxide losses of northern ecosystems in response to autumn warming. *Nature* **451**, 49–52, https://doi.org/10.1038/nature06444 (2008).

26. Jeong, S. J. et al. Accelerating rates of Arctic carbon cycling revealed by long-term atmospheric CO₂ measurements. *Sci. Adv.* **4**, eaao1167 (2018).

27. Belshe, E. F., Schuur, E. A. & Bolker, B. M. Tundra ecosystems observed to be CO₂ sources due to differential amplification of the carbon cycle. *Ecol. Lett.* **16**, 1307–1315, https://doi.org/10.1111/ele.12164 (2013).

28. Plaza, C. et al. Direct observation of permafrost degradation and rapid soil carbon loss in tundra. *Nat. Geosci.* **12**, 627 (2019).

29. Mauritz, M. et al. Nonlinear CO₂ flux response to 7 years of experimentally induced permafrost thaw. *Glob. Chang Biol.* **23**, 3646–3666, https://doi.org/10.1111/gcb.13661 (2017).

30. Turetsky, M. R. et al. Carbon release through abrupt permafrost thaw. *Nat. Geosci.* **13**, 138–143, https://doi.org/10.1038/s41561-019-0526-0 (2020).

31. Schadel, C. et al. Potential carbon emissions dominated by carbon dioxide from thawed permafrost soils. *Nat. Clim. Change* **6**, 950–953, https://doi.org/10.1038/Nclimate3054 (2016).

32. Hugelius, G. et al. Large stocks of peatland carbon and nitrogen are vulnerable to permafrost thaw. *Proc. Natl Acad. Sci. USA* **117**, 20438–20446, https://doi.org/10.1073/pnas.1916387117 (2020).

33. Ballantyne, A. P., Alden, C. B., Miller, J. B., Tans, P. P. & White, J. W. Increase in observed net carbon dioxide uptake by land and oceans during the past 50 years. *Nature* **488**, 70–72, https://doi.org/10.1038/nature11299 (2012).

34. Zhu, Z. C. et al. Greening of the Earth and its drivers. *Nat. Clim. Change* **6**, 791, https://doi.org/10.1038/Nclimate3004 (2016).

35. Bjorkman, A. D. et al. Status and trends in Arctic vegetation: evidence from experimental warming and long-term monitoring. *Ambio* **49**, 678–692, https://doi.org/10.1007/s13280-019-01161-6 (2020).

36. Myers-Smith, I. H. et al. Eighteen years of ecological monitoring reveals multiple lines of evidence for tundra vegetation change. *Ecol. Monogr.* **89**, e01351, https://doi.org/10.1002/ecm.1351 (2019).

37. Bjorkman, A. D. et al. Plant functional trait change across a warming tundra biome. *Nature* **562**, 57–62, https://doi.org/10.1038/s41586-018-0563-7 (2018).

38. Jansen, E. et al. Past perspectives on the present era of abrupt Arctic climate change. *Nat. Clim. Change* **10**, 714–721, https://doi.org/10.1038/s41558-020-0860-7 (2020).

39. Pulliainen, J. et al. Patterns and trends of Northern Hemisphere snow mass from 1980 to 2018. *Nature* **581**, 294–298, https://doi.org/10.1038/s41586-020-2258-0 (2020).

40. Lian, X. et al. Summer soil drying exacerbated by earlier spring greening of northern vegetation. *Sci. Adv.* **6**, eaax0255 (2020).

41. Welp, L. R., Randerson, J. T. & Liu, H. P. The sensitivity of carbon fluxes to spring warming and summer drought depends on plant functional type in boreal forest ecosystems. *Agric. For. Meteorol.* **147**, 172–185 (2007).

42. Song, X. P. et al. Global land change from 1982 to 2016. *Nature* **560**, 639–643, https://doi.org/10.1038/s41586-018-0411-9 (2018).
43. Frost, G. V. & Epstein, H. E. Tall shrub and tree expansion in Siberian tundra ecotones since the 1960s. *Glob. Chang Biol.* **20**, 1264–1277, https://doi.org/10.1111/gcb.12406 (2014).

44. Liu, X. et al. Enhanced nitrogen deposition over China. *Nature* **494**, 459–462, https://doi.org/10.1038/nature11917 (2013).

45. Rogers, B. M., Soja, A. J., Goulden, M. L. & Randerson, J. T. Influence of tree species on continental differences in boreal fires and climate feedbacks. *Nat. Geosci.* **8**, 228–234, https://doi.org/10.1038/ Ngeo2352 (2015).

46. Liu, Z., Ballantyne, A. P. & Cooper, L. A. Biophysical feedback of global forest fires on surface temperature. *Nat. Commun.* **10**, 214, https://doi.org/10.1038/s41467-018-08237-z (2019).

47. Rees, W. G. et al. Is subarctic forest advance able to keep pace with climate change? *Glob. Change Biol.* **26**, 3965–3977, https://doi. org/10.1111/gcb.15113 (2020).

48. Scheffer, M., Hirota, M., Holmgren, M., Van Nes, E. H. & Chapin, F. S. Thresholds for boreal biome transitions. *PNAS* **109**, 21384–21389, https://doi.org/10.1073/pnas.1219844110 (2012).

49. Barichivich, J. et al. Large-scale variations in the vegetation growing season and annual cycle of atmospheric CO2 at high northern latitudes from 1950 to 2011. *Glob. Change Biol.* **19**, 3167–3183 (2013).

50. Keenan, T. F. et al. Net carbon uptake has increased through warming-induced changes in temperate forest phenology. *Nat. Clim. Change* **4**, 598–604, https://doi.org/10.1038/ Nclimate2253 (2014).

51. Ueyama, M., Iwata, H. & Harazono, Y. Autumn warming reduces the CO2 sink of a black spruce forest in interior Alaska based on a nine-year eddy covariance measurement. *Glob. Chang Biol.* **20**, 1161–1173, https://doi.org/10.1111/gcb.12434 (2014).

52. Zhang, Y., Parazoo, N. C., Williams, A. P., Zhou, S. & Gentine, P. Large and projected strengthening moisture limitation on end-of-season photosynthesis. *Proc. Natl Acad. Sci. USA* **117**, 9216–9222, https:// doi.org/10.1073/pnas.1914436117 (2020).

53. Zhang, Y., Commane, R., Zhou, S., Williams, A. P. & Gentine, P. Light limitation regulates the response of autumn terrestrial carbon uptake to warming. *Nat. Clim. Change* **10**, 739–743, https://doi.org/ 10.1038/s41558-020-0806-0 (2020).

54. Thomas, R. T. et al. Increased light-use efficiency in northern terrestrial ecosystems indicated by CO2 and greening observations. *Geophys. Res. Lett.* **43**, 11339–11349, https://doi.org/10.1002/ 2016gl070710 (2016).

55. Peng, S. S. et al. Benchmarking the seasonal cycle of CO2 fluxes simulated by terrestrial ecosystem models. *Glob. Biogeochemical Cycles* **29**, 46–64, https://doi.org/10.1002/2014gb004931 (2015).

56. Crowther, T. W. et al. The global soil community and its influence on biogeochemistry. *Science* **365**, eaav0550 (2019).

57. Endsley, K. A., Kimball, J. S. & Reichle, R. H. Soil Respiration Phenology Improves Modeled Phase of Terrestrial net Ecosystem Exchange in Northern Hemisphere. *J. Adv. Modeling Earth Syst.* **14**, e2021MS002804 (2022).

58. Brown, J., Ferrians Jr, O., Heginbottom, J. & Melnikov, E. *Circum-Arctic map of permafrost and ground-ice conditions*. (US Geological Survey Reston, VA, 1997).

59. Peylin, P. et al. Global atmospheric carbon budget: results from an ensemble of atmospheric CO2 inversions. *Biogeosciences* **10**, 6699–6720, https://doi.org/10.5194/bg-10-6699-2013 (2013).

60. Peters, W. et al. An atmospheric perspective on North American carbon dioxide exchange: CarbonTracker. *Proc. Natl Acad. Sci. USA* **104**, 18925–18930, https://doi.org/10.1073/pnas. 0708986104 (2007).

61. Peters, W. et al. Seven years of recent European net terrestrial carbon dioxide exchange constrained by atmospheric observations. *Glob. Change Biol.* **16**, 1317–1337, https://doi.org/10.1111/j. 1365-2486.2009.02078.x (2010).

62. Chevallier, F. et al. CO2 surface fluxes at grid point scale estimated from a global 21 year reanalysis of atmospheric measurements. *Journal of Geophysical Research-Atmospheres* **115**, https://doi.org/ 10.1029/2010jd013887 (2010).

63. Rödenbeck, C., Conway, T. & Langenfelds, R. The effect of systematic measurement errors on atmospheric CO 2 inversions: a quantitative assessment. *Atmos. Chem. Phys.* **6**, 149–161 (2006).

64. Rödenbeck, C., Houweling, S., Gloor, M. & Heimann, M. CO 2 flux history 1982–2001 inferred from atmospheric data using a global inversion of atmospheric transport. *Atmos. Chem. Phys.* **3**, 1919–1964 (2003).

65. Saeki, T. & Patra, P. K. Implications of overestimated anthropogenic CO 2 emissions on East Asian and global land CO 2 flux inversion. *Geosci. Lett.* **4**, 9 (2017).

66. Kolby Smith, W. et al. Large divergence of satellite and Earth system model estimates of global terrestrial CO2 fertilization. *Nat. Clim. Change* **6**, 306–310, https://doi.org/10.1038/nclimate2879 (2016).

67. Pinzon, J. E. & Tucker, C. J. A non-stationary 1981–2012 AVHRR NDVI3g time series. *Remote Sens.* **6**, 6929–6960 (2014).

68. Heinsch, F. A. et al. Evaluation of remote sensing based terrestrial productivity from MODIS using regional tower eddy flux network observations. *Ieee Trans. Geosci. Remote Sens.* **44**, 1908–1925, https://doi.org/10.1109/Tgrs.2005.853936 (2006).

69. Turner, D. P. et al. Evaluation of MODIS NPP and GPP products across multiple biomes. *Remote Sens. Environ.* **102**, 282–292, https://doi.org/10.1016/j.rse.2006.02.017 (2006).

70. Zhao, M. S., Heinsch, F. A., Nemani, R. R. & Running, S. W. Improvements of the MODIS terrestrial gross and net primary production global data set. *Remote Sens. Environ.* **95**, 164–176, https:// doi.org/10.1016/j.rse.2004.12.011 (2005).

71. Ballantyne, A. et al. Accelerating net terrestrial carbon uptake during the warming hiatus due to reduced respiration. *Nat. Clim. Change* **7**, 148, https://doi.org/10.1038/Nclimate3204 (2017).

72. Pettorelli, N. et al. Using the satellite-derived NDVI to assess eco-logical responses to environmental change. *Trends Ecol. Evol.* **20**, 503–510, https://doi.org/10.1016/j.tree.2005.05.011 (2005).

73. Reichstein, M. et al. On the separation of net ecosystem exchange into assimilation and ecosystem respiration: review and improved algorithm. *Glob. Change Biol.* **11**, 1424–1439, https://doi.org/10.1111/ j.1365-2486.2005.001002.x (2005).

74. Lasslop, G. et al. Separation of net ecosystem exchange into assimilation and respiration using a light response curve approach: critical issues and global evaluation. *Glob. Change Biol.* **16**, 187–208, https://doi.org/10.1111/j.1365-2486.2009.02041.x (2010).

75. Barr, A. G. et al. Use of change-point detection for friction-velocity threshold evaluation in eddy-covariance studies. *Agric. For. Meteorol.* **171–172**, 31–45, https://doi.org/10.1016/j.agrformet.2012. 11.023 (2013).

76. Sitch, S. et al. Recent trends and drivers of regional sources and sinks of carbon dioxide. *Biogeosciences* **12**, 653–679, https://doi. org/10.5194/bg-12-653-2015 (2015).

77. Jones, L. A. et al. The SMAP level 4 carbon product for monitoring ecosystem land-atmosphere CO2 exchange. *IEEE Trans. Geosci. Remote Sens.* **55**, 6517–6532 (2017).

78. Chapin, F. S. et al. Reconciling carbon-cycle concepts, terminol-ogy, and methods. *Ecosystems* **9**, 1041–1050, https://doi.org/10. 1007/s10021-005-0105-7 (2006).

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

80. Obu, J. et al. Northern Hemisphere permafrost map based on TTOP modelling for 2000–2016 at 1 km² scale. *Earth-Sci. Rev.* **193**, 299–316 (2019).

81. Gruber, A., Scanlon, T., van der Schalie, R., Wagner, W. & Dorigo, W. Evolution of the ESA CCI Soil Moisture climate data records and
their underlying merging methodology. *Earth Syst. Sci. Data* **11**, 717–739, https://doi.org/10.5194/essd-11-717-2019 (2019).

82. Fan, Y. et al. Applications of structural equation modeling (SEM) in ecological studies: an updated review. *Ecol. Process.* **5**, 19 (2016).

83. Rosseel, Y. Lavaan: An R package for structural equation modeling and more. Version 0.5–12 (BETA). *J. Stat. Soft.* **48**, 1–36 (2012).

84. Piao, S. L. et al. Weakening temperature control on the interannual variations of spring carbon uptake across northern lands. *Nat. Clim. Change* **7**, 359–363, https://doi.org/10.1038/Nclimate3277 (2017).

## Acknowledgements

Z.L. acknowledges support from the NSFC (41922006), CAS Project for Young Scientists in Basic Research (YSBR-037), and Major Program of Institute of Applied Ecology, Chinese Academy of Sciences (IAEMP202201); B.M.R. acknowledges support from the National Aeronautics and Space Administration (NASA) Arctic-Boreal Vulnerability Experiment (ABOVE) and Carbon Cycle Science programs (NNX17AE13G); J.S.K. acknowledges support from NASA Terrestrial Ecology and MEaSUREs programs (80NSSC19M0114, 80NSSC18K0980); B.M.R. and S.M.N. acknowledge NASA ABoVE support (NNX15AT81A); J.D.W. and A.-M.V. acknowledge support from the Gordon and Betty Moore Foundation. J.D.W. acknowledges support from the NASA New Investigator Program (NNH17ZDA001N). A portion of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. N.C.P., N.M., and Z.L. acknowledge support from the NASA Earth Science Division IDS program; P.C. and A.B. acknowledge supports from the ESA CCI RECCAP2 project (ESRIN/4000123002/18/I-NB); A.P.B. supported by MOPGA ANR 18-MPGA-0007. Copyright 2020. All rights reserved. N.C. is supported by the Arctic Challenge for Sustainability II grant (JPMXD1420318865) of the Ministry of Education, Science, Culture and Sports (MEXT), Japan. CarbonTracker CT2019B results provided by NOAA ESRL, Boulder, Colorado, USA from the website at http://carbontracker.noaa.gov.

## Author contributions

Z.L. conceived the study and carried out the analysis with input from A.P.B. and W.W. A.B. conducted the mixed effects model analysis of net CO₂ trends over different regions. Z.L. wrote the paper with input from

J.K., A.P.B., N.P., W.W., A.B., N.M., S.N., J.W., B.R., P.C., K.Y., A.V., F.C., W.P., P.P., and N.C.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-022-33293-x.

**Correspondence** and requests for materials should be addressed to Zhihua Liu, John S. Kimball, Ashley P. Ballantyne or Wen J. Wang.

**Peer review information** *Nature Communications* thanks the anonymous reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022

<footnote n="1">¹Numerical Terradynamic Simulation Group, WA Franke College of Forestry and Conservation, University of Montana, Missoula, MT, USA. ²CAS Key Laboratory of Forest Ecology and Management, Institute of Applied Ecology, Chinese Academy of Sciences, Shenyang, Liaoning, China. ³Department of Ecosystem and Conservation Sciences, University of Montana, Missoula, MT, USA. ⁴Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, Gif-sur-Yvette, France. ⁵Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA. ⁶Northeast Institute of Geography and Agroecology, Chinese Academy of Sciences, Changchun, Changchun, Jilin, China. ⁷Max Planck Institute for Biogeochemistry, Department of Biogeochemical Integration, Jena, Germany. ⁸Woodwell Climate Research Center, Falmouth, MA, USA. ⁹Meteorology and Air Quality Group, Wageningen University and Research, Wageningen, the Netherlands. ¹⁰Research Institute for Global Change, Japan Agency for Marine-Earth Science and Technology (JAMSTEC), Yokohama, Japan. ¹¹University, Centre for Isotope Research, Groningen, the Netherlands. ¹²These authors contributed equally: Zhihua Liu, Ashley P. Ballantyne. e-mail: liuzh811@126.com; John.Kimball@umontana.edu; ashley.ballantyne@umontana.edu; wangwenj@iga.ac.cn</footnote>