Increasing terrestrial ecosystem carbon release in response to autumn cooling and warming

1College of Global Change and Earth System Science, Beijing Normal University, Beijing, China. 2Department of Physical Geography and Ecosystem Science, Lund University, Lund, Sweden. 3Regional Climate Group, Department of Earth Sciences, University of Gothenburg, Gothenburg, Sweden. 4State Key Laboratory of Desert and Oasis Ecology, Xinjiang Institute of Ecology and Geography, Chinese Academy of Sciences, Urumqi, China. 5College of Water Sciences, Beijing Normal University, Beijing, China. 6School of Atmospheric Sciences, Sun Yat‐Sen University, Guangzhou, China. 7School of Geography and Tourism, Qufu Normal University, Rizhao, China. 8School of Geography, Beijing Normal University, Beijing, China. 9Key Laboratory of Water Cycle and Related Land Surface Processes, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing, China. 10Ministry of Ecology and Environment Center for Satellite application on Ecology and Environment, Beijing, China. 11Advanced Institute of Natural Sciences, Beijing Normal University at Zhuhai, Zhuhai, China. ✉e-mail: hebin@bnu.edu.cn; chenyn@ms.xjb.ac.cn Terrestrial ecosystems in the middle and high latitudes play a critical role in regulating the global carbon cycle and atmospheric CO2 concentration1,2. In recent decades, the overall climate warming has strongly enhanced ecosystem productivity by extending the growing season and increasing photosynthesis, although this warming effect on ecosystem carbon uptake has weakened since the latter half of the 1990s3–5. The Northern Hemisphere (NH) ecosystems act as a carbon sink in spring but a net carbon source in autumn due to the asymmetric response of ecosystem productivity and respiration in response to warming2. In spring, climate warming more strongly promotes photosynthesis than respiration, resulting in net carbon uptake, but an opposite condition occurs in autumn. A recent study reported widespread autumn cooling in the NH since 20046, and the ecosystem respiration is anticipated to decrease with cooling. Therefore, we hypothesize that the autumn carbon release would be suppressed by the emerging cooling. However, for those regions that experienced continued warming, an enhanced autumn carbon loss would be expected. To test this hypothesis, we use independent simulations and observation-derived estimates of net ecosystem exchange (NEE) to comprehensively examine the interannual variation of ecosystem carbon fluxes in response to different autumn temperature changes. The NEE datasets include nine simulations from the TRENDY ensemble of global vegetation models (Supplementary Table 1)7, one simulation from empirical models based on eddy-covariance observational data using machine learning algorithms (FLUXCOM)8–10 and three estimates from atmospheric CO2 inversion models: CarboScope11, CAMS12 and NISMON13. In addition, in situ atmospheric CO2 concentration records and eddy-covariance CO2 flux observations (Supplementary Table 2) were used as observational evidence to indicate net land–atmosphere carbon exchange14. Satellite-based normalized difference vegetation index (NDVI) and contiguous solar-induced chlorophyll fluorescence (CSIF), as well as a near-infrared reflectance(NIRv-) based gross primary productivity (GPP) product (defined as NIRv_GPP), were also used to indicate the ecosystem productivity change.

T errestrial ecosystems in the middle and high latitudes play a critical role in regulating the global carbon cycle and atmospheric CO 2 concentration 1,2 . In recent decades, the overall climate warming has strongly enhanced ecosystem productivity by extending the growing season and increasing photosynthesis, although this warming effect on ecosystem carbon uptake has weakened since the latter half of the 1990s [3][4][5] . The Northern Hemisphere (NH) ecosystems act as a carbon sink in spring but a net carbon source in autumn due to the asymmetric response of ecosystem productivity and respiration in response to warming 2 . In spring, climate warming more strongly promotes photosynthesis than respiration, resulting in net carbon uptake, but an opposite condition occurs in autumn. A recent study reported widespread autumn cooling in the NH since 2004 6 , and the ecosystem respiration is anticipated to decrease with cooling. Therefore, we hypothesize that the autumn carbon release would be suppressed by the emerging cooling. However, for those regions that experienced continued warming, an enhanced autumn carbon loss would be expected.
To test this hypothesis, we use independent simulations and observation-derived estimates of net ecosystem exchange (NEE) to comprehensively examine the interannual variation of ecosystem carbon fluxes in response to different autumn temperature changes. The NEE datasets include nine simulations from the TRENDY ensemble of global vegetation models (Supplementary Table 1) 7 , one simulation from empirical models based on eddy-covariance observational data using machine learning algorithms (FLUXCOM) [8][9][10] and three estimates from atmospheric CO 2 inversion models: CarboScope 11 , CAMS 12 and NISMON 13 . In addition, in situ atmospheric CO 2 concentration records and eddy-covariance CO 2 flux observations (Supplementary Table 2) were used as observational evidence to indicate net land-atmosphere carbon exchange 14 . Satellite-based normalized difference vegetation index (NDVI) and contiguous solar-induced chlorophyll fluorescence (CSIF), as well as a near-infrared reflectance-(NIRv-) based gross primary productivity (GPP) product (defined as NIRv_GPP), were also used to indicate the ecosystem productivity change.

Widespread cooling has halted overall autumn warming
Air temperature trends derived from the Climate Research Union (CRU) 15 for the period from 1981 to 2018 indicate an apparent shift from warming to weak cooling in autumn around 2004 north of 25° N (Fig. 1a,b). The turning point of temperature that occurred in 2004 (P = 0.001) was confirmed by both piecewise linear regression and a moving t test (Methods). Before 2004, the average autumn temperature increased at a rate of 0.45 ± 0.14 °C decade -1 (Mann-Kendall, P = 0.004), which we define as the warming period (WP), whereas post-2004, there was a weak decreasing trend at a rate of −0.07 ± 0.13 °C decade -1 (P = 0.593), defined as the cooling period (CP). Cooling trends were widely observed in the high latitudes of North America and central Eurasia, accounting for about 51% of the total study area. This cooling pattern is consistent with the study of ref. 6 and was attributed to the strengthening of the Pacific Decadal Oscillation and Siberian High.
For the convenience of description, we designated the regions where the autumn temperature shows a decreasing trend in the period 2004-2018 relative to the period 1982-2003 as the cooling Increasing terrestrial ecosystem carbon release in response to autumn cooling and warming areas (CAs). The warming areas (WAs) were defined as regions outside the CAs ( Supplementary Fig. 1). To further verify whether the spatial partition is reasonable, we used the empirical orthogonal function (EOF) method to extract the dominant temporal and spatial characteristics of temperature change ( Supplementary  Fig. 2). We selected the top three EOFs to represent the dominant spatiotemporal distribution features. The first mode represents an overall warming pattern. The second mode reflects mainly the spatial distribution of temperature change affected by interdecadal climate fluctuations such as the Pacific Decadal Oscillation 6,16,17 . The third mode is roughly consistent with our partition of CAs and WAs, probably representing an interannual change in temperature itself after removing long-term global warming and interdecadal climate fluctuations.
Before the turning point in 2004, the CAs underwent warming at a rate of 0.40 ± 0.19 °C decade -1 (Mann-Kendall, P = 0.052), but in recent years it shifted to a strong cooling trend at a rate of -0.69 ± 0.22 °C decade -1 (P = 0.008) (Fig. 1c). The WAs account for 49% of the total land area and maintained a warming trend over the whole period under consideration but exhibited a stronger warming rate after the turning point. In addition to the average temperature, we also examined the maximum and minimum temperature trends, and observed a similar cooling pattern, with a stronger cooling rate for the maximum temperature (−0.76 ± 0.23 °C decade -1 , P = 0.005) than for the minimum temperature (−0.61 ± 0.22 °C decade -1 , P = 0.016) ( Supplementary Fig. 3). Thus, the decrease of temperature during the day in the last decade is faster than that during the night in the CAs.

Increasing land carbon exchange in recent decades
Multiple independent simulations of ecosystem carbon fluxes all indicate an increasing autumn NEE during the CP (2004-2018) relative to the WP , suggesting an increasing loss of ecosystem carbon in autumn (Fig. 2a). Significant correlations were observed among different sources of NEEs ( Supplementary Fig. 4), suggesting different models agree with each other in the simulation of the interannual variation of NEE. During the WP, NEE increased with a relatively weak positive trend, while a strong upward trend was observed during the CP. The NEE derived from process-based models (TRENDY S2 simulation) and an empirical model (FLUXCOM) both suggest a significant NEE increase from the WP to the CP, with the mean increasing rate from 0.03 ± 0.03 PgC decade -1 (Mann-Kendall, P = 0.326) to 0.11 ± 0.04 PgC decade -1 (P = 0.020) and from 0.01 ± 0.01 PgC decade -1 (P = 0.603) to 0.06 ± 0.02 PgC decade -1 (P = 0.018), respectively (Fig. 2b). A significant difference (P < 0.01) in NEE trends between the WP and CP was suggested by the Student's t test. The NEEs from atmospheric CO 2 inversion models (CarboScope, CAMS and NISMON) also indicate increased NEE rates during the CP relative to the WP, with significant trends found for CAMS (0.35 ± 0.09 PgC decade -1 , P = 0.001) and NISMON (0.17 ± 0.08 PgC decade -1 , P = 0.059 for the CP). To highlight the trend changes of NEE, the ratios between the slopes of NEE in the later and earlier periods were calculated. The largest relative increase in NEE trend was found with FLUXCOM, followed by CarboScope and TRENDY (Fig. 2b). In addition, the NEE averaged over all data sources suggest that the mean autumn CO 2 release offset about 24.8% of spring and summer carbon uptake in the WP and 25.2% in the CP ( Supplementary Fig. 5), implying an increasing role of autumn carbon loss in the annual carbon exchange.
To test the robustness of our findings, we further divided the study area into four sub-regions according to the significance of temperature change after 2004 and tested the NEE trend for each region ( Supplementary Fig. 6). We observed that significant cooling resulted in an increase of NEE release compared with that before 2004, and significant warming maintained the increased rate of NEE. In addition, we also identified the temperature turning point at each grid point and examined the NEE change under two different conditions of temperature shift, the continued warming and the shift from warming to cooling after the turning point (Methods). Both of these two types of temperature shifts led to an increase of NEE after the temperature turning point ( Supplementary Fig. 7), suggesting that our findings are not subject to the selection of the turning point of temperature change.
The ecosystem exchanges in the CAs and WAs were studied separately to discern their contributions to the overall NEE trend in the NH (Supplementary Figs. 8 and 9). After the temperature turning point, the NEE trends in both the CAs and WAs increased significantly, with stronger NEE trends in the WAs (from 0.01 ± 0.03 to 0.07 ± 0.04 PgC decade -1 , significantly different at P < 0.01 in TRENDY according to the Student's t test) than in the CAs (from 0.02 ± 0.01 to 0.03 ± 0.01 PgC decade -1 , significantly different at P < 0.01 according to the Student's t test). The relative contributions of the CAs and WAs to the overall NEE trend over the whole study region were 27% and 63%, respectively, implying that the observed enhanced NEE trends were driven mainly by carbon exchange in the WAs.
Evidence from atmospheric CO 2 observations further confirms the increasing ecosystem carbon release. We used observations of CO 2 from Point Barrow to calculate the upward zero-crossing date of CO 2 in autumn (AZC; Methods), which has generally been recognized as an integrative indicator of autumn atmosphere-land carbon exchange over the northern middle and high latitudes [18][19][20] , with earlier AZC indicating larger NEE. When correlating the AZC from Point Barrow with various sources of NEE data averaged northwards of 30° N, strong negative correlations were observed (Supplementary Fig. 10), suggesting that the different models performed well in capturing the variation of carbon exchange over the studied region. In response to warming, the AZC has occurred earlier due to a larger increase of total ecosystem respiration (TER) compared with GPP 2,21 . The AZC shows a persistent advancing trend from 1982 to 2018 ( Supplementary Fig. 10b), suggesting an increasing loss of CO 2 from northern middle-and high-latitude ecosystems in autumn. In addition, the carbon release in autumn, defined as the difference between the CO 2 concentration on the day of the AZC and the CO 2 concentration in the first week of September, also shows a greater upward trend during the CP than the WP, confirming the increase in autumn carbon loss over the second period (2004-2018). We furthermore analysed eddy-covariance measurements of CO 2 fluxes from 48 flux towers with 10 yr continuous data from the Integrated Carbon Observation System 2018 and the FLUXNET Network 2015 (Supplementary Table 2). A majority of the sites have observations around 2000 corresponding to the CP. An increasing autumn NEE was observed in 9 out of 17 sites in the CAs and in 19 out of 31 sites in the WAs ( Supplementary Fig. 11a). Because CO 2 emissions from fires and land use can also affect NEE signals, we used the fourth-generation global fire emissions database to determine the possible impact of fire emissions on the NEE trend. Above 25° N, the average autumn CO 2 emissions from fires (September to November) was 0.50 PgC, accounting for only 4.6% of the mean NEE calculated from the TRENDY simulations ( Supplementary Fig. 12). In addition, during the period 2004-2015, CO 2 emissions from fires did not change significantly. Therefore, we suggest that the influence of fire emissions on NEE in autumn was negligible. We also analysed the NEE changes under the TRENDY S3 mode, which considers emissions from land-use change ( Supplementary Fig. 13). The S2 and S3 simulations show a comparable change in trends during both the CP and the WP, implying a limited contribution of land-use change to the increased NEE.

Terrestrial processes caused increased autumn carbon release
To explore possible reasons for the increased autumn CO 2 release, we examined the NEE changes in terms of GPP and TER derived from the TRENDY simulations. The GPP and TER trends for the entire area north of 25° N, the CAs and the WAs were determined separately for the WP and CP ( Supplementary Fig. 6). We attributed the trend in the net carbon release during the autumn CP in the NH to the additive effect of a decrease in GPP growth (from 0.30 ± 0.08 to 0.26 ± 0.07 PgC decade -1 ) and an increase in TER growth (from 0.34 ± 0.07 to 0.37 ± 0.08 PgC decade -1 ). Although the carbon release increased in both the CAs and WAs, the NEE increases were driven by different contributions from changes in GPP and TER. In the CAs, the growth rates of both GPP and TER decreased with decreasing temperature, but the decrease in GPP (0.13 ± 0.05 to 0.09 ± 0.06 PgC decade -1 ) was greater than that in TER (0.15 ± 0.05 to 0.12 ± 0.02 PgC decade -1 ), leading to an increase of CO 2 release. In the WAs, the growth rate of GPP (0.17 ± 0.05 PgC decade -1 ) did not change between the WP and CP, but the TER trend increased from 0.18 ± 0.03 PgC decade -1 for the WP to 0.25 ± 0.06 PgC decade -1 for the CP, resulting in enhanced carbon release.
Despite a similar increase in NEE during the CP, the changes of GPP and TER revealed by carbon fluxes from FLUXCOM are somewhat different from those shown by the TRENDY simulations ( Supplementary Fig. 9). FLUXCOM indicates a shift in the GPP trend from an increase to a decrease in the NH and the CAs and an unchanged growth rate in the WAs after 2004. The increased NEE during the CP shown by FLUXCOM is due mainly to the larger decrease in the GPP trend (from 0.02 ± 0.02 to −0.04 ± 0.03 PgC decade -1 ) than in the TER trend (from 0.04 ± 0.03 to 0.02 ± 0.03 PgC decade -1 ). The difference in carbon fluxes, especially the difference in GPP trends between FLUXCOM and TRENDY models, is probably because the FLUXCOM models do not consider the fertilization effect of carbon dioxide 10 , which has been recognized as an important process in the ecosystem carbon cycle. Despite the discrepancy in GPP trends between TRENDY and FLUXCOM, they both reveal that the contrasting response of GPP and TER to temperature led to a systematic increase in NEE in the NH.

Response of land carbon flux to autumn temperature variation
To explore how the temperature change regulated the carbon exchange processes, we performed partial correlations between temperature and carbon fluxes (GPP, TER and NEE) from TRENDY simulations while controlling for the effects of soil moisture (SM), vapour-pressure deficit (VPD) and cloudiness ( Fig. 3 and Supplementary Table 3). Cloudiness was used here to indicate changes in solar radiation. Following previous investigations 5, 19 , we detrended all the variables before the partial correlation analysis using a linear regression method. The temperature was significantly positively correlated with GPP and TER in both the WP and the CP, especially in the higher latitudes, implying a dominant influence of temperature on GPP and TER variations and thereby on NEE (Fig. 3). Despite the strong correlations between temperature and both GPP and TER, the NEE-temperature correlations are relatively weak, as has also been reported by ref. 22 . This discrepancy can be explained by the offset effect existing between GPP and TER. For example, warming would promote both GPP and TER, but its impact on NEE depends on the balance between GPP and TER. We calculated the autumn NEE-temperature correlation relying on eddy-covariance measurements from 48 flux towers. Both positive and negative relationships were observed, and significant correlations were found in only a few stations (Supplementary Fig. 11b).
During both periods, the regional mean GPP and TER maintained a high correlation with temperature in the CAs. In the WAs, however, the significant positive correlation between GPP and temperature in the WP changed to a weak negative correlation in the CP (Supplementary Table 3). This shift was associated with a persistent significant positive GPP-SM correlation during both periods, probably because a large proportion of the WAs features water-limited drylands, such as Inner Mongolia, the north of Africa and the southwest of America. We speculate that the disappeared GPP-temperature correlation during the CP could be attributed to the warming-induced heat stress on vegetation growth as reported by ref. 4 or high VPD constraint on plant photosynthesis associated with high temperatures, as indicated by satellite observations (Supplementary Table 5). We observed a strong dependence of the GPP-temperature correlation on mean autumn temperature ( Supplementary Fig. 14a). With the increase in mean autumn temperature, the GPP-temperature correlations decreased from significant positive correlations to weak negative ones, implying a warming constraint on the GPP-temperature relationship. Simultaneously, we also observed an increasing GPP-SM correlation with the increase of mean autumn temperature ( Supplementary Fig. 14b), probably implying an enhanced water demand with warming conditions 4,23,24 . The change of the TER-temperature relationship is similar to the GPP-temperature relationship in the WAs during the CP, demonstrating a reduced temperature control of TER, which was partially replaced by a strong SM constraint, as indicated by the significant TER-SM correlation (P < 0.05). Combined with the observed increased GPP-SM correlation, we speculate that the emerging strong TER-SM correlation during the CP was caused by increasing soil moisture, which maintains the growth of ecosystem productivity and in turn increases the plant-derived carbon inputs to soils for respiration [25][26][27] . The larger increase in TER than GPP in autumn during the CP could be partly explained by the increased carbon input from spring and summer, as observed increased growth rate of GPP in spring and summer relative to the WP in the WAs (Supplementary Fig. 15). Thus, the warming and wetting conditions stimulated a larger increase of TER than GPP during CP relative to WP (Supplementary Fig. 8).
The preceding analysis suggests that the cooling suppressed both the GPP and TER growths in the CAs while the persistent wetting maintained the GPP growth and stimulated even more rapid growth of TER in the WAs (Supplementary Fig. 15). A reduced temperature impact and enhanced water regulation on the GPP and TER in the WAs during the CP was also found in an analysis based on FLUXCOM data (Supplementary Fig. 16 and Supplementary Table 4), confirming the change in the climatic control of the carbon exchange processes.

ecosystem productivity change in autumn from satellite data
The preceding analysis suggests a slowdown of autumn GPP growth, primarily in the CAs during the CP. Three sets of indicators, including the satellite-based NDVI derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) (2004-2018), CSIF and NIRv_GPP were used to further examine the ecosystem productivity change. During the WP, NIRv_GPP indicates that the NH experienced widespread GPP increases in autumn (Fig. 4a), which can be explained by the delayed autumn phenology caused by warming as reported by previous studies 28,29 . During the CP, we found a strong decline in productivity in central Asia and southern Europe in the CAs, which are roughly consistent with the GPP trends simulated by FLUXCOM and TRENDY. Although browning was also revealed by MODIS NDVI in some CAs during the CP, greening trends were predominant over the NH. Similar to MODIS NDVI, the CSIF data reveal widespread increases in plant photosynthesis over the NH during the CP. Differences in trends among satellite observations are expected because a range of factors could trigger uncertainties in inversions of optical remote sensing. We also cannot simply attribute the difference between satellite observations and model simulations to either side, highlighting the urgent need for improving both satellite observations and model simulations to reach robust conclusions.
The partial correlations between satellite indices and climatic variables were also examined by controlling for the influence of SM, cloudiness and VPD ( Supplementary Fig. 17 and Supplementary Table 5). The results partly support our findings, such as (1) the widespread positive correlations in the CAs between satellite indices and temperature, especially for MODIS NDVI, imply a constraint of temperature change on ecosystem productivity; (2) the absence of a positive impact of temperature on ecosystem productivity in the WAs during the CP (Supplementary Table 5). We failed to observe a strong impact of SM on ecosystem productivity; however, a relatively strong negative impact from VPD was found in the WAs during the CP. Further analysis suggests that both temperature and SM were positively coupled with VPD ( Supplementary Fig. 18), implying that the warming and wetting conditions during the CP in the WAs promoted VPD, and thus the high VPD dominates the change of the ecosystem productivity. Our previous study also revealed a strong impact of VPD on global ecosystem productivity and carbon exchange 30 . However, a more in-depth study is needed to isolate the independent influence of VPD on ecosystem productivity from the impacts of temperature and SM.
Previous studies have suggested that the northern ecosystems, stimulated by climate warming, act as a net carbon sink in spring but as a source in autumn 2,21 . To our knowledge, little attention has been paid in the literature to the impact of autumn cooling on ecosystem carbon exchange. In this study, we show that multiple independent lines of evidence point to an increasing carbon loss in autumn in response to cooling over the northern middle-and high-latitude ecosystems. Cooling can reduce ecosystem respiration but also inhibits productivity. Ecosystems that experienced warming also exhibited increased carbon release due to increased ecosystem respiration and limited productivity increase because the warming and wetting conditions were more favourable to TER growth than to GPP increase. If the temperature in the cooling areas continues to drop and the temperature in the warming areas continues to rise in the coming decade, the carbon release in the autumn will probably play a more dominant role in offsetting the carbon absorption in the spring and summer, which could lead to a weakened carbon absorption capacity of the northern ecosystems. methods Climate data. Monthly climate data (air temperature at 2 m and cloudiness) with a spatial resolution of 0.5° were obtained from the CRU Time Series 4.0. 15 We extracted data from 1982 to 2018 to match the time series of satellite vegetation observations. The VPD was calculated as the difference between saturated water-vapour pressure and actual water-vapour pressure 31 . Temperature and vapour-pressure data used for the VPD calculation were obtained from CRU.
Soil moisture data. The daily root-zone soil moisture with a spatial resolution of 0.25° for the period 1980-2018 was obtained from the Global Land Evaporation Amsterdam Model (GLEAM v.3.3a) 32 . The dataset is based on radiation and air temperature from a reanalysis, a combination of gauge-based, reanalysis-based and satellite-based precipitation and satellite-based vegetation optical depth.
Fire emission data. Monthly carbon emissions from biomass burning were obtained from the fourth-generation Global Fire Emission Database 33 . This dataset has a spatial resolution of 0.25° and provides global data on the burning area and emissions on three-hourly, daily and monthly timescales and estimates the contributions of different fire types. Emissions data can be obtained for different substances, such as carbon (C), dry matter (DM), carbon dioxide (CO 2 ), carbon monoxide (CO) and methane (CH 4 ).
Satellite vegetation greenness data. The satellite-based NDVI archived from the MODIS NDVI dataset with a spatial resolution of 0.5° and a temporal resolution of 16 days was used here to detect vegetation greenness changes. In addition, the solar-induced chlorophyll fluorescence product was used as a proxy of vegetation photosynthesis. We furthermore used the four-day clear-sky CSIF time series (2000-2019) with a spatial resolution of 0.05° × 0.05° from ref. 34 (https://osf.io/8xqy6/).

GPP based on NIRv.
The NIRv is a newly developed satellite vegetation index combining NDVI and near-infrared band reflectivity of vegetation and is recognized as a proxy of GPP 35,36 . We obtained the 0.05° NIRv_GPP from 1982 to 2018 from ref. 37 . This product was produced by upscaling the relationships between NIRv and observed GPP to the global scale and was judged to perform well in capturing interannual trends of GPP 37 .
Atmospheric CO 2 data. In situ observations of daily CO 2 concentration at Point Barrow were obtained from the National Oceanic and Atmospheric Administration/Earth System Research Laboratory network. According to analyses of atmospheric transport and mixing processes, the CO 2 signals detected at Barrow are suggested to be an integrated measure of carbon fluxes over both the high latitudes and the middle latitudes 20 .

Ecosystem carbon fluxes.
Simulations of ecosystem carbon fluxes (GPP, TER and NEE) derived from process-based model simulations (TRENDY), empirical models based on flux tower observations (FLUXCOM) and atmospheric CO 2 inversion models were jointly used for the investigation of net ecosystem carbon exchange over the northern middle and high latitudes.
The TRENDY dataset is an ensemble of dynamic global vegetation model (DGVM) simulations that are forced by CRU-National Centers for Environmental Prediction historical climate and CO 2 inputs 38 . The DGVMs use a bottom-up approach to simulate terrestrial CO 2 fluxes (for example, GPP, TER and NEE), and were extensively used to explore the mechanisms driving changes in carbon uptake and fluxes. The simulated GPP, TER and NEE from nine models of TRENDYv.8 (Supplementary Table 1) were used in this study. The S2 experiment, which considered the effect of both observed changes of CO 2 and climate on ecosystem carbon fluxes, was selected for studying the changes of ecosystem carbon fluxes before and after the temperature shift.
The FLUXCOM dataset is an upscaling product using empirical models forced by eddy-covariance data from 224 flux towers, remote sensing data and climate data [8][9][10] . It provides estimates of global energy and carbon fluxes (http:// www.fluxcom.org/). The empirical models were trained by three machine learning algorithms, including Random Forests, Artificial Neural Networks and Multivariate Adaptive Regression Spline, and thus provide a series of estimates of global carbon fluxes. We used the FLUXCOM carbon fluxes data driven by the European Centre for Medium-Range Weather Forecasts Reanalysis v.5 (ERA5) climate reanalysis from 1979 to 2018.
The atmospheric CO 2 inversion datasets provide estimates of NEE over land from long-term atmospheric CO 2 measurements using atmospheric transport models. Three atmospheric CO 2 inversion products were used here: monthly net biome production with a spatial resolution of 3.75° × 2.5° from the JENA CarboScope (version s76_vo2020) for the period 1976-2019, long-term global CO 2 fluxes estimated by the NICAM-based Inverse Simulation for Monitoring CO 2 (NISMON-CO 2 ) between 1990 and 2019 and the Copernicus Atmosphere Monitoring Service 12 (CAMS v.19r1) dataset between 1979 and 2019.
Eddy-covariance CO 2 observation data. The eddy-covariance measurements of carbon fluxes from tower sites were obtained from the Integrated Carbon Observation System 2018 and the FLUXNET Network 2015. We selected 48 eddy-covariance CO 2 observation sites with 10 yr continuous data (Supplementary Table 2) located north of 25° N and extracted temperature and NEE data from September to November to explore the change of ecosystem carbon exchange in autumn.

NEE estimation.
The monthly NEE was estimated as the difference between TER and GPP. The autumn (September to November) GPP and TER derived from TRENDY and FLUXCOM over the study region were obtained by aggregating GPP and TER from each grid cell weighted by the grid-cell area. The NEE derived from atmospheric CO 2 inversions was directly used and compared against those from TRENDY and FLUXCOM. To compare the NEE before and after the temperature turning point, we divided the NEE time series into two periods: 1982-2003 and 2004-2018. Calculation of the AZC. We used observations of CO 2 from Point Barrow to characterize the trends in the zero-crossing date of CO 2 (downward in spring and upward in autumn). These trends roughly correspond to the beginning of net carbon uptake in spring and the beginning of net carbon release in autumn. According to the method of ref. 39 , we obtained the detrended seasonal CO 2 curve by separating the seasonal cycle from the long-term trend and short-term variations, fitting a function consisting of a quadratic polynomial for the long-term trend and four harmonics for the annual cycle to the daily data. The residuals from this function fit are then obtained. A 1.5-month and a 390-day full-width half-maximum-value averaging filter were used for the digital filtering of residuals to remove the short-term variations and the long-term trend, respectively. Then we got the zero-crossing dates when the detrended seasonal CO 2 curve crosses the zero line from positive to negative and negative to positive, respectively.
The autumn carbon release is calculated as the amount of CO 2 released between the autumn zero-crossing date and the first week of September following ref. 21 .
Identification of turning point of temperature. We used the piecewise linear regression method to determine the turning point of the mean autumn (September to November) temperature during 1982-2018 over the area north of 25° N. In addition, a moving t-test method was used to verify the turning-point identification. Then, the temporal trends of the mean autumn temperature before and after the turning point were calculated using the Mann-Kendall non-parametric trend test method, and the confidence intervals were determined using Sen's slope statistics. According to the temperature trends before and after the turning point, we further identified the CAs as where the autumn temperature shows a decreasing trend after the turning point (2004) relative to that before the turning point, and WAs as regions outside the CAs. To maintain spatial integrity and continuity, we ignored the significance of the temperature trend when dividing the CAs and WAs.
To verify that our analysis is not affected by the division of the time period and regions, we also identified the temperature turning point at each grid point using the piecewise linear regression method and then extracted those grid points with significant temperature change and significant NEE change (P < 0.1) before and after the turning point. In addition, to ensure an adequate length of the time series for the trend analysis, only grid points with at least 15 yr of time series of temperature records and NEE data before and after the turning point were selected. For NEE, only TRENDY and CAMS had long enough time series to satisfy this criterion. Then, the selected grid points were classified into two groups, the warming group (shift from significant warming to a larger significant warming rate after the turning point) and the cooling group (shift from significant warming to significant cooling after the turning point), and the NEE trends before and after temperature turning point were compared.

Statistical analysis.
Trend statistics and slopes of other climatic variables, carbon fluxes (GPP, TER, and NEE) and satellite indices during autumn were also determined using Mann-Kendall non-parametric trend test and Sen's slope statistics. The Student's t test was used to determine the significance of the trend difference in carbon fluxes during different periods.
The EOF decomposition method was used to explore the temporal and spatial variation of autumn temperature change. This method can determine the main components of variables on a time-space scale and provide the variance contribution of different distribution characteristics 40 . Here, we decomposed autumn temperatures from 1982 to 2018 using EOF and selected the top three distribution modes according to the variance contribution degree, which represents the main spatial pattern and temporal changes of temperature from 1982 to 2018.
The actual value of NEE differs greatly among different sources of simulations. Consistent with previous studies, different sources of NEEs were normalized using the Z-score method to highlight the interannual variability and trend of NEE change.
Pearson's correlation analysis and partial correlation analysis were performed to examine the relationship between climatic variables and carbon fluxes (GPP, TER and NEE) and satellite indices. We detrended all the variables before the partial correlation analysis using a linear regression method. The significance of correlations was assessed at P < 0.1.