Global variability of belowground autotrophic respiration in terrestrial ecosystems

Abstract. Belowground autotrophic respiration (RA) is one of the largest, but highly uncertain carbon flux components in terrestrial ecosystems. It has not been explored globally before and still acted as a “black box” in global carbon cycling. Such progress and uncertainty motivate a development of global RA dataset and understand its spatial and temporal pattern, causes and responses to future climate change. This study used Random Forest to study RA's spatial and temporal pattern at the global scale by linking the updated field observations from Global Soil Respiration Database (v4) with global grid temperature, precipitation and other environmental variables. Globally, mean RA was 43.8 ± 0.4 Pg C a−1 with a temporally increasing trend of 0.025 ± 0.006 Pg C a−1 over 1980–2012. Such increment trend was widely spread with 58 % global land areas. For each 1 °C increase in annual mean temperature, global RA increased by 0.85 ± 0.13 Pg C a−1, and it was 0.17 ± 0.03 Pg C a−1 for 10 mm increase in annual mean precipitation, indicating a positive feedback of RA to future climate change. At a global scale, precipitation was the main dominant climatic drivers of the spatial pattern of RA, accounting for 56 % of global land areas with widely spread globally, particularly in dry or semi-arid areas, followed by shortwave radiation (25 %) and temperature (19 %). Different temporal patterns for varying climate zones and biomes indicated uneven response of RA to future climate change, challenging the perspective that the parameters of global carbon stimulation independent on climate zones and biomes. The developed RA database, the missing carbon flux component that is not constrained and validated in terrestrial ecosystem models and earth system models, will provide insights into understanding mechanisms underlying the spatial and temporal variability of belowground carbon dynamics. RA database also has great potentials to serve as a benchmark for future data-model comparisons. The RA product is freely available at https://doi.org/10.6084/m9.figshare.7636193.



1 Introduction
Belowground autotrophic respiration (RA) mainly originated from plant roots, mycorrhizae, and other micro-organisms in the rhizosphere directly relying on labile carbon component leaked from roots (Hanson et al., 2000;Tang et al., 2016;Wang and Yang, 2007). Thus, RA reflects the photosynthesis derived carbon respired back to the atmosphere by roots and regulates the net photosynthetic production 45 allocation to belowground tissues (Högberg et al., 2002). RA is one main component of soil respiration (Hanson et al., 2000), and soil respiration represents the second largest source of carbon fluxes from soil to the atmosphere (after gross primary production, GPP) in the global carbon cycle (Raich and Schlesinger, 1992). Globally, RA could amount roughly up to 54 Pg C yr -1 (1 Pg = 10 15 g, calculating RA as an approximate ratio of 0.5 of soil respiration, more details in Hanson et al., 2000) according to 50 different estimates of global soil respiration (Bond-Lamberty, 2018), which is almost 5 times of the carbon release from human activities (Le Quéré et al., 2018). However, the contribution of RA to soil respiration varied greatly from 10% to 90% across biomes, climate zones and among years (Hanson et al., 2000), leading to the strong spatial and temporal variability in RA. Thus, whether RA varies with ecosystem types or climate zones remains an open question at the global scale .

55
Consequently, an accurate estimate of RA and its spatial-temporal dynamics are critical to understand the response of terrestrial ecosystems to global carbon cycling and climate change.
Due to the difficulties in separation and measurement of RA at varying spatial scales and its diurnal, seasonal and annual variabilities, RA becomes one of the largest but highly uncertain carbon flux components in terrestrial ecosystems. Although individual site measurements of RA have been conducted 60 across ecosystem types and biomes, the globally spatial and temporal pattern of RA has not been explored and still acts as a "black box" in global carbon cycling . This "black box" is not well constrained and validated, because most terrestrial ecosystem models and earth system models were commonly calibrated and validated against eddy covariance measurements of net ecosystem carbon exchange (Yang et al., 2013). Such progress and uncertainty motivate a development of global RA dataset 65 from observations and understand its spatial and temporal pattern, causes and responses to future climate change. Despite of the general agreement that global soil respiration increased during last several decades Bond-Lamberty and Thomson, 2010b;Zhao et al., 2017), how global RA responding to climate change is far from certain because of different temperature sensitivities of RA 4 across terrestrial ecosystems (Liu et al., 2016;Wang et al., 2014). Therefore, reducing RA uncertainty 70 and clarifying its response to climate change, particularly to temperature and precipitation, is essential for global carbon allocation and future projection of the impact of climate change on global terrestrial carbon cycle.
Although several studies have globally estimated soil respiration and its response to climate variables (Bond-Lamberty and Thomson, 2010b;Hursh et al., 2017;Zhao et al., 2017), such efforts have not been 75 conducted for global RA directly. Although Hashimoto et al. (2015) indirectly derived RA via the difference between total soil respiration and heterotrophic respiration, it probably led to uncertainty due to include the temperature and precipitation as the only model drivers and a low model efficiency (32%).
Besides temperature and precipitation, other variables, e.g. soil water, carbon and nitrogen content, are additionally critical factors regulating RA, and those factors generally varied with biomes and climate 80 zones. Consequently, Hashimoto et al. (2015) may not reflect the key processes affecting RA, such as soil nutrient limitation.
On the other hand, the climate-derived models usually explain < 50% variability of soil respiration (Bond-Lamberty and Thomson, 2010b;Hashimoto et al., 2015;Hursh et al., 2017), which might be another uncertainty source. Recent studies have included more variables and field observations to 85 improve the prediction ability of linear and non-linear models (Jian et al., 2018b;Zhao et al., 2017), however, it may propagate error because of the overfitting and autocorrelation among these variables (Long and Scott, 2006). Random Forest (RF, Breiman, 2001), a machine learning approach, could overcome these issues based on the hierarchical structure, and be insensitive to outliers and noise compared to single classifiers (Breiman, 2001;Tian et al., 2017). RF uses a large number of ensemble 90 regression trees but a random selection of predictive variables (Breiman, 2001). RF only requires two free parameter settings: the number of variables sampled as candidates for each split and the number of trees. The performance of the RF model usually is not sensitive to the number of trees and number of variables. Moreover, RF regression can deal with a large number of features, which could help feature selection based on the variable importance and can avoid overfitting (Jian et al., 2018b). RF has been 95 widely used for carbon fluxes modelling in recent years (Bodesheim et al., 2018;Jung et al., 2017). Therefore, we applied RF algorithm to retrieve global RA based on the updated RA field observations from the most updated global soil respiration dataset (SRDB v4, doi: 10.5194/bg-7-1915-2010, Bond-5 Lamberty and Thomson, 2010a with the linkage of other global variables (see "materials and methods" part) in this study for the first time, aiming to: (1) develop a global RA product using field observations 100 across the globe (named data-derived RA); (2) estimate RA's spatial and temporal pattern; (3) identify the main driving factors of RA's spatial and temporal variabilities; (4) compare with the previous RA estimates from Hashimoto et al. (2015). The outcome of this study will advance our understanding of global RA and its spatial and temporal variabilities. The proposed RA product is expected to serve as a benchmark for global vegetation models and its role in global carbon cycling. It will also advance our 105 knowledge of the co-variation of RA with climate, soil and vegetation factors, further link the empirical observations temporally and spatially to bridge the knowledge gap among local, regional and global scales. To control the data quality, annual RA observations were filtered that: (1) annual RA was directly reported 115 in publications indicated by "years of data" of SRDB; (2) the start and end years were recorded in literatures or expanded from "years of data" of SRDB; (3) soil respiration measurements with Alkali absorption and soda lime were not included due to the potential underestimate of respiration rate with the increasing pressure inside chamber (Pumpanen et al., 2004); (4) observations with treatments of nitrogen addition, air/soil warming, and rain/litter exclusion were not included, except cropland; (5) 120 potential problems observations (labelled by "Q10", "Q11", "Q12", "Q13" and "Q14") were excluded.
Finally, this study included a total of 449 field observations ( Fig. 1

Vegetation, climate and soil data
Global land cover with a half degree resolution was obtained from MODIS land cover (https://glcf.umd.edu/data/lc/). Monthly grid data of temperature, precipitation, diurnal temperature 130 range and potential evapotranspiration at 0.5 o resolution were obtained from CRU TS V4.01 from 1901 to 2016 (https://crudata.uea.ac.uk) (Harris et al., 2014). Monthly shortwave radiation (SWR), Palmer Drought Severity Index (PDSI) and soil water content at 0.5 o resolution were from NOAA/ESRL Physical Sciences Division (https://www.esrl.noaa.gov) (Kalnay et al., 1996). Soil organic carbon content with a resolution of 250 m was downloaded from soil grid data (https://soilgrids.org) (Hengl et al., 2017), 135 and soil nitrogen content was from Spatial Data Access Tool (https://webmap.ornl.gov/ogc/index.jsp), while monthly nitrogen deposition data with a resolution of 0.5 o were downloaded from the Earth System Models of GISS-E2-R, CCSM-CAM3.5 and GFDL-AM3, covering since 1850s (https://www.isimip.org). The monthly global variables were first aggregated to year scale and then resampled to a 0.5 o resolution using bilinear interpolation for those variables without a 0.5 o resolution.

140
These variables could represent different aspects controlling RA variability. For instance, temperature, precipitation and soil water content are most important variables controlling plant photosynthesis, which is the primary carbon source of RA (Högberg et al., 2002;Högberg et al., 2001). Finally, global variables of each given site extracted by coordinates corresponding with annual RA estimates from the SRDB.

145
In this study, a RF model was trained with the 11 variables listed in Table S1 by caret by linking RandomForest package in R 3.4.4 (Kabacoff, 2015), then the trained model was implemented to estimate grid RA at 0.5° × 0.5° resolution over 1980-2012. The performance of RF was assessed by a 10-fold cross-validation (CV). A 10-fold CV suggested that the whole dataset was subdivided into 10 parts with approximately an equal number of samples. The target values for each of these 10 parts were predicted 150 on the training using the remaining nine parts. Two statistics were employed in model assessment: modelling efficiency (R 2 ) and root mean square error (RMSE) (Yao et al., 2018). The 10-fold CV result showed that RF performed well and could capture the spatial-and temporal-pattern of RA ( Fig. S1 in supplementary materials).

155
This study applied Theil-Sen linear regression to estimate temporal trend analysis of RA and its driving variables for each grid cell. The Theil-Sen estimator is a median-based non-parametric slope estimator, which has been widely used for spatial analysis of time series carbon flux analysis (Forkel et al., 2016;Zhang et al., 2017). Mann-kendall non-parametric test was applied for the significant change trend in RA and its driving factors for each grid cell (p < 0.05).

Relationships between RA and climate variables
In this study, mean annual temperature, mean annual precipitation and mean annual shortwave radiation were considered as the most important proxies driving RA. The relationships between RA and temperature, precipitation and shortwave radiation were analyzed by partial correlation for each grid cell.
The absolute value of the correlation coefficient of these three variables was used in RGB combination 165 to indicate the dominant factors of RA.

The comparison map profile method
In order to compare with the solely global RA product generated by Hashimoto et al. (2015), which was estimated by a climate-driven model using temperature and precipitation only and obtained from the public available dataset (http://cse.ffpri.affrc.go.jp/shojih/data/index.html), this study applied the 170 comparison map profile (CMP) method. CMP was developed based on absolute distance (D) and crosscorrelation coefficient (CC) through multiple scales (Gaucherel et al., 2008). D and CC reflect the similarity of data values and spatial structure of two images with the same size, respectively (Gaucherel et al., 2008). Low D and higher CC reflects goodness between the compared images, and vice versa. The 8 D among moving windows of two compared images was calculated by equation (1) (Gaucherel et al., 175 2008): ̅̅̅and ̅ are averages calculated over two moving windows. This study used 3×3 to 41×41 pixels.
Finally, the mean D was averaged for different scales.
The CC was calculated by equation (2) (Gaucherel et al., 2008): (2) (3) Where xij and yij are the pixel values at row i and column j of two moving windows of the two compared images, respectively. N represents the number of pixels for each moving window, while σx and σy are the standard deviation calculated from the two moving windows. Finally, like D calculations, CC was 185 calculated as the mean of different scales. The data-derived RA in this study presented a great globally spatial variability during 1980-2012 ( Fig.   2a and Fig. 3). Largest RA fluxes commenced from tropical regions, particularly in Amazon tropical and Southeast areas, where generally have a high RA > 700 g C m -2 yr -1 . Following the tropical areas, 195 subtropics, e.g. South China, East America, and humid temperate areas, e.g. North America, West and Middle Europe, had typical moderate RA fluxes of 400-600 g C m -2 yr -1 . By contrast, the relative low RA fluxes occurred in the areas with sparse vegetation cover, cold and dry climate, e.g. boreal and tundra, which had low temperature and short growing season. Besides, dry or semi-arid areas, e.g. Northwest

Spatial patterns of RA
China and Middle East, also had typical low RA fluxes below 200 g C m -2 yr -1 , where were often limited 200 by water availability. The most significant RA inter-annual variability (expressed by standard deviation, Fig. 2c) was found 205 in topical or subtropical regions with above 80 g C m -2 yr -1 , while most areas remained less variable with less than 40 g C m -2 yr -1 . Latitudinally, zonal mean RA increased from cold and dry biomes (Tundra and semi-arid) to warm and humid biomes (temperate and tropical forests, Fig. 3), reflecting from more to less environmental limitations. RA varied from 112±21 g C m -2 yr -1 at about 70 o N to 552±101 g C m -2 yr -10 occurred in tropical areas and middle Australia. However, it is worth noted that some clear differences between data-derived and Hashimoto RA existed (Fig. 4): specifically, there was a remarkable difference of above 300 g C m -2 yr -1 for South Amazon and 200 g C m -2 yr -1 for subtropical China. Although most areas between data-derived RA and Hashimoto RA expressed high and positive correlations, some areas, such as Middle East, West Russia and East America and North Japan, showed negative correlations.  Hashimoto RA during 1980 The trend of data-derived RA showed heterogeneous patterns in spatial (Fig. 5). A total of 58% of global areas experienced an increasing trend during 1980-2012 (calculating from cell areas), and 33% of these areas showed a significant change (p < 0.05). Generally, the change trend for the majority areas 230 was from -4 to 4 g C m -2 yr -2 , while the most striking increasing change occurred in East Russia, tropical, and Eastern regions in Africa with an increasing trend of above 5 g C m -2 yr -2 . Similarly, 77% of global areas of Hashimoto RA had an increasing trend, 46% of which were statistically significant (p < 0.05).

| Total RA and its temporal trend
11 235 Figure 6 Annual variability of belowground autotrophic respiration (RA) for this study (a) and Hashimoto RA (b) from 1980 to 2012. The grey area represents 95% confidence interval.
Mean global RA was 43.8±0.4 Pg C yr -1 during 1980-2012, varying from 42.9 Pg C yr -1 in 1992 to 44.9 Pg C yr -1 in 2010, with a significant trend of 0.025±0.006 Pg C yr -2 despite of high annual variabilities (0.06% yr -1 , p < 0.001, Fig. 6a). Similarly, a rising trend was also observed for Hashimoto 240 RA, however, its annual increasing trend (0.073±0.009 Pg C yr -2 , p < 0.001, Fig. 6b) was higher than that of data-derived RA. Annual mean of Hashimoto RA was 40.5±0.9 Pg C yr -1 . RA and its trend was also evaluated for three climate zones (boreal, temporal and tropical areas based on Köppen-Geiger climate classification) and eight major biomes (boreal forest, cropland, grassland, 250 savannas, shrubland, temperate forest, tropical forest and wetland, Fig. 7). Tropics had highest RA of 15.6±0.2 Pg C yr -1 , followed by temperate regions with 9.3±0.1 Pg C yr -1 , and boreal areas represented the lowest RA of 6.7±0.1 Pg C yr -1 . These three climate zones were main contributors of global RA, accounting for 72%. Temporally, considerable RA inter-annual variability of these three climate zones existed (Fig. S2). Specifically, RA in tropical and boreal zones showed a significantly increasing trend 255 from 1980 to 2012, with an increasing rate of 0.013±0.003 and 0.008±0.002 Pg C yr -2 , respectively.
However, RA in temperate zones presented a slightly decreasing trend of -0.003±0.001 Pg C yr -2 (p = 0.048) although strong variability was observed.
In terms of biomes, tropical forest had the highest RA, followed by the widely distributed cropland and savannas (Fig. 7), while wetland had the lowest RA due to its limited land cover. RA showed a 260 significantly increasing trend during 1980-2012 (ps < 0.01) in majority biomes, except temperate forest, savannas and wetland. RA in tropical forests, boreal forests and cropland increased by 0.0076±0.0015, 0.0047±0.0016, 0.0036±0.0014 Pg C yr -2 , respectively. Compared to data-derived RA, Hashimoto RA for the three climate zones and eight biomes generally produced similar change patterns, although the magnitude difference existed (Figs. 7, S2 and S3). However, there were significant increasing trends of 265 RA in temperate zones, temperate forest, savannas and wetland of Hashimoto RA, which were not observed in data-derived RA.
RA was significantly correlated with temperature anomaly (R 2 = 0.59, p < 0.001) and precipitation anomaly (R 2 = 0.50, p < 0.001, Fig. 8). On average, RA increased by 0.85±0.13 Pg C yr -2 for 1 o C increment in mean annual temperature, and 0.17±0.03 Pg C yr -2 for 10 mm increase in mean annual 270 precipitation. However, different biomes and climate zones showed uneven responses to the temperature and precipitation change ( Fig. S4 and S5). For example, no significant correlations were found between RA in temperate zone/savannas/wetland and temperature anomaly, while other climate zones and biomes were significantly correlated with temperature/precipitation anomaly. 13 275 Figure 8 The relationships between belowground autotrophic respiration (RA) and the anomaly of temperature and precipitation for this study (a, b) and Hashimoto (c, d) respectively. The anomaly was calculated as the difference between temperature or precipitation of corresponding year to the mean of 1980-2012. *** means significant level at 0.001. The dominant environmental factor was examined with partial regression coefficients when regressing RA against annual mean temperature, annual mean precipitation and shortwave radiation. Latitudinally, 285 higher mean annual temperature, precipitation and shortwave radiation were associated with higher RA in the major latitudinal gradients (positive partial correlations, Fig. S6). Spatially, the dominant environmental factor varied greatly globally (Fig. 9). Precipitation was the most important dominant factor for the spatial pattern of RA among the three environmental controls, covering about 56% of global 14 land areas (Fig. S7), which was widely distributed globally, particularly in dry or semi-arid areas, such 290 as Northwest China, Southern Africa, Middle Australia and America. Temperature dominated about 19% of global land areas, which mainly occurred in tropical Africa, Southern Amazon rainforests, Siberia and partly tundra. The rest land area (25%) was dominated by shortwave radiation, primarily covering boreal areas above 50 o N, Eastern America and middle and Eastern Russian. Similarly, precipitation was also the most important dominant factor for Hashimoto RA, dominating about 77% land area, while 295 temperature and shortwave radiation dominated 13% and 10% land area. However, their spatial patterns varied greatly compared to data-derived RA. For example, temperature was the main dominant factor for most area of Australia for Hashimoto RA, while data-derived RA indicated that precipitation and shortwave radiation dominated such areas (Fig. 9).

310
Globally, mean annual RA amounted to 43.8±0.4 Pg C yr -1 from 1980 to 2012 (Fig. 6). It was slightly higher than Hashimoto RA (40.5±0.9 Pg C yr -1 ), and there was great divergence of spatial and temporal patterns (see discussion part in "Comparison with Hashimoto RA"). Due to no direct estimate on global RA, this study compared other RA estimates using total soil respiration multiplied by the proportion of RA or heterotrophic respiration. Bond-Lamberty et al. (2018) proposed that the global average proportion 315 of heterotrophic respiration ranged from 0.54 to 0.63 over 1990-2014 and global soil respiration was 67 to 108 Pg C yr -1 using different approaches and datasets Bond-Lamberty (2018); (Bond-Lamberty and Thomson, 2010b;Hashimoto et al., 2015;Hursh et al., 2017;Jian et al., 2018b) , thus global RA varied from 25 to 51 Pg C yr -1 . RA estimate in this study fell in this range. Similarly, during 1980 increased by 0.025±0.006 Pg C yr -2 , which may be related to the increasing photosynthesis due to global 320 warming and CO2 fertilization effects, which could increase carbon availability in plant-derived substrate inputs into the soil (e.g. root exudates and biomass) for both root metabolism (Piñeiro et al., 2017;Zhou et al., 2016). Such annual increase accounted for about 25% of global soil respiration increase (0.09 and 0.1 Pg C yr -2 ) (Bond-Lamberty and Thomson, 2010b;Hashimoto et al., 2015), suggesting that about one quarter of the total soil respiration increment due to climate change came from RA. With 1 o C increase 325 in global mean temperature, RA will increase by 0.85±0.13 Pg C yr -2 and 0.17±0.03 Pg C yr -2 for 10 mm increase in precipitation, which indicated that carbon fluxes from RA might positively feedback to future global climate change, which was typically characterized by increasing temperature and changes in precipitation (IPCC, 2013). However, RA increment varied with climate zones and ecosystem types (Figs. S2 and S3), which was similar to previous findings Jian et al., 2018a), who 330 found that total soil respiration or RA varied with climate zones or ecosystem types. These differences may be related to regional heterogeneity and plant functional trait. For example, regional temperature significantly differed from global averages (Huang et al., 2012), with much faster change in high-latitude regions (Hartmann et al., 2014), and semi-arid dominated the trend and variability of global land CO2 sink (Ahlström et al., 2015). Therefore, the regionally uneven responses of RA to climatic variables 335 highlights the urgent need to account for regional heterogeneity when studying the effects of climate change on ecosystem carbon dynamics in future. RA estimate in this study also has important indications of carbon allocation from photosynthesis. The immediate carbon substrates for RA were primarily derived from recent photosynthesis (Högberg et al., 2001;Subke et al., 2011). Strong correlation between photosynthesis and RA demonstrated the evidence for their close coupling relationships (Chen et al., 340 2014;Kuzyakov and Gavrichkova, 2010). Globally, GPP was about 125 Pg C yr -1 during last few decades (Bodesheim et al., 2018;Zhang et al., 2017). Thus, roots respired more than one third carbon from GPP, suggesting that except the carbon used for constructing belowground tissues, a large proportion of carbon will be returned back to atmosphere respired by roots. However, it should be noted that through root respiration, soil nutrients for vegetation growth will be required, which may affect the RA flux.

Dominant factors
Spatially, the driving factors for RA varied greatly. Temperature and shortwave radiation were the main driving factors for high latitudinal areas above 50 o N (Figure 9a). This result was not surprising because RA was positively correlated with temperature or photosynthesis (indirectly reflecting the solar radiation) (Chen et al., 2014;Tang et al., 2016), and high latitudinal regions was always limited by 350 temperature or energy, leading to low RA as well (Fig. 3a).
Globally, precipitation was the most important factor, covering about 56% of land area (Fig. 9a and   S7). Precipitation was always considered as a proxy for soil water content (Hursh et al., 2017;Yao et al., 2018), and such wide dominance of precipitation on RA was related to the mechanisms of the soil water availability driving RA. First, soil water exists in form of ice when temperature is below zero, that plant 355 and soil microbes could not directly use for growth or respiration. This could be observed in some boreal areas where precipitation was the dominant factor of RA (Fig. 9a). Second, too high or too low soil water content (e.g. flooding and drought) could limit the mobility of substrates and carbon input to belowground, which could affect RA. Yan et al. (2014) found that soil respiration decreased once soil water content was below a lower (14.8 %) or above an upper (26.2%) threshold in a poplar 360 plantation. Similarly, Gomez-Casanovas et al. (2012) also found that RA decreased when soil water content was above 30%. These results seemed to support the finding in this study. Third, the relationship between soil water content and RA or total soil respiration is more complex than the relationship between temperature and soil respiration. Numerous formula, such as linear (Tang et al., 2016), polynomial (Moyano et al., 2012), logarithmic (Schaefer et al., 2009), quadratic (Hursh et al., 2017) models have 365 been widely applied to describe the relationship between soil water content and soil respiration. The multifarious relationships between soil water content and RA may occur because soil water content affect RA in multiple ways. Meanwhile, seasonal variability of precipitation and soil water content is often correlated with temperature (Feng and Liu, 2015), making the relationship between soil water content and RA more complex.

370
Similarly, the dominance of precipitation in Hashimoto was also widely observed (Fig. 8), dominating 77% of land area (Fig. S7). Although this percentage was 17% higher than data-derived RA, both results demonstrated that global RA of the majority land cover was dominated by precipitation. However, it is noticeable that dominant environmental factor controlling spatial carbon fluxes gradient may differ among different years (Reichstein et al., 2007), e.g. climate extreme and disturbance. 375 17

Comparison with Hashimoto RA
Globally, total data-derived RA was slightly higher than Hashimoto RA, however, great divergence was observed both spatially and temporally (Fig. 6), particularly in tropical regions, where data-derived RA was much lower than Hashimoto RA (Fig. 3). These differences could be attributed to several reasons.
First, two RA products had different land cover areas, especially in desert areas in North Africa, where 380 existed very sparse or no vegetation. If data-derived RA was masked by Hashimoto RA, global RA was 39.6±0.4 Pg C yr -1 , which was pretty close to Hashimoto RA (Fig. S8). Second, different predictors and algorithms were applied for data-derived RA and Hashimoto RA prediction. Besides temperature and precipitation, RA was also affected by soil nutrient, carbon substrate supply, belowground carbon allocation, site disturbance and other variables (Chen et al., 2014;Hashimoto et al., 2015;Tang et al., 385 2016;Zhou et al., 2016). Hashimoto RA was calculated from the difference between total soil respiration and heterotrophic respiration, which were predicted by a simple climate-driven model using temperature and precipitation only (Hashimoto et al., 2015). Thus Hashimoto RA could not reflect its soil nutrient and other environmental constrains. To overcome such limitations, beside temperature and precipitation, this study included soil water content, soil nitrogen and soil organic carbon as proxies for environmental 390 and nutrient constraints of RA and considered the interactions among these variables using RF. This study indeed improved the model efficiency to 0.52 for RA prediction (Fig. S1), which was 0.32 for Hashimoto soil respiration (Hashimoto et al., 2015). The simple climate-model for Hashimoto soil respiration could be its advantages and limitations (Hashimoto et al., 2015). Third, the empirical model (the relationship between total soil respiration and heterotrophic respiration) deriving Hashimoto RA originated from 395 forest ecosystems (Bond-Lamberty et al., 2004;Hashimoto et al., 2015), which may bring uncertainties to other ecosystems. For example, the difference between data-derived RA and Hashimoto RA varied up to 350 g C yr -1 in South, North Amazon areas and Madagascar, where the savannas widely distributed (Fig. 4), thus Hashimoto RA might not capture the spatial and temporal pattern of RA for non-forest ecosystems. Including more environmental variables and improving algorithm could be a good option to 400 reduce the uncertainty in modelling RA.

Advantages, limitations and uncertainties
Generally, this study had four main advantages to estimate global RA: first, this study, to our knowledge, was the first attempt to model RA using a large number of empirical field observations, and 18 estimate spatial and temporal patterns globally. While most previous studies mainly focused on global 405 soil respiration, which was not partitioned into RA and heterotrophic respiration globally (Hursh et al., 2017;Jian et al., 2018b;Zhao et al., 2017). Second, this study used an up-to-date field observational database developed from SRDB up to the end 2018. This new updated database included a total of 449 field observations (Fig. 1). These observations had a wide coverage range of global terrestrial ecosystems and represented all major biomes and climate zones. Third, the global terrestrial ecosystems were 410 separated into eight biomes, including boreal forest, cropland, grassland, savannas, shrubland, temperate forest, tropical forest and wetland. The total RA and its inter-annual variability were evaluated for each of the eight biomes ( Fig. S3 and S4). Besides, total RA and its inter-annual variability was also assessed for three climate zones -boreal, temperate and tropical zones (Figs. S2 and S5), according to the Köppen-Geiger climate classification system (Peel et al., 2007). These were important climate zones, contributing 415 72% of global RA. Different temporal change trends across biomes and climate zones also further indicated an uneven response of RA to climate change. Fourth, this study used a RF algorithm to model and map global RA with the linkage of climate and soil predictors. The results showed that RF could accurately estimate the relationships between annual RA and predictors (Fig. S1). Compared to linear regressions for soil respiration prediction (because no global RA prediction before this study) with a 420 model efficiency less than 35% (Bond-Lamberty and Thomson, 2010b;Hashimoto et al., 2015;Hursh et al., 2017), RF algorithm achieved a much higher model efficiency to 52%, which indeed improved the RA modelling and reduced the uncertainties.
Although data-derived global RA could serve as a benchmark for global carbon cycling modelling, and this study had filled the data-gaps of global RA, limitations and uncertainties still remained in few 425 aspects. First, although we conducted a data quality control in this study, a lack of reliable approach to separate RA and heterotrophic respiration may lead to an uncertainty of RA values. There are several approaches, e.g. trenching, stable or radioactive isotope, gridding (Bond-Lamberty et al., 2004;Högberg et al., 2001;Hanson et al., 2000), however, each of these approaches has its own limitations. For example, trenching has been widely applied to partition RA and heterotrophic respiration due to easy operation 430 and low cost, on the other hand, heterotrophic respiration may be increased due to the termination of water uptake by roots and the decomposition of remaining dead roots in trenching plots (Hanson et al., 2000;Tang et al., 2016). Commonly, RA was calculated from the difference between total soil respiration and heterotrophic respiration, thus the trenching approach might lead to an underestimation of RA. In our dataset, a total of 254 RA observations were estimated by trenching approach, while the rest RA 435 observations were estimated by other separation approaches, e.g. isotope, radiocarbon, mass balance.
Thus, inconsistent separation approaches could be another source of uncertainty of RA values.
Second, due to the limited observations of RA at a daily or monthly scale, this study only predicted RA at an annual scale. Although there was no direct study to compare the difference of RA upscaling from daily or monthly and annual scale, substantial difference of soil respiration upscaling from daily or 440 monthly and annual scales (Jian et al., 2018b) indirectly illustrated the potential difference of RA upscaling from different timescales.
Third, the effects of rising atmospheric CO2 on root growth was not explicitly represented in this study, although CO2 fertilization effects could partly be represented in the increase temperature. While the magnitude of CO2 fertilization effects on photosynthesis is still uncertain (Gray et al., 2016), RF or other 445 machine learning approaches are encouraged to quantify the uncertainties due to CO2 fertilization.
Fourth, this study did not consider the effects of human activities and historical changes in biomes on RA. However, important changes may occur in tropical forest, grassland and cropland during last several decades due to human activities (Hansen et al., 2013;Klein Goldewijk et al., 2011). Thus, changes in biomes should be included in future global RA and carbon cycling modelling. However, the lack of such 450 data is the main constrain of detecting the effects of biome change on RA.
Finally, uneven coverage of observations in the updated database would be another source of uncertainties. Although our dataset had a wide range of land cover, the observational sites mainly distributed in China, Europe and North America and were dominated by forests. There was a great lack of observations in areas, such as Africa, Austria and Russia, and biomes, such as tropical forest, shrubland, 455 wetland and cropland. Consequently, RA observations caused bias of RF model toward the regions with more observations. Therefore, including more observations in these areas and biomes should largely increase our capability to assess the spatial and temporal patterns of global RA and contribute to improve the global carbon cycling modelling to future climate change.

Conclusions
Although data-derived global RA may serve as a benchmark for ecosystem models, no such study has assessed global variability in RA with a large number of empirical observations that can help bridge the 465 knowledge gap between local, regional, and global scales. This study has filled this knowledge gap by linkage of field observations and global variables using RF algorithm, providing an annual global RA product at a of 0.5 o × 0.5 o resolution from 1980 to 2012. Currently, robust findings include: (1) Annual mean RA was 43.8±0.4 Pg C yr -1 with a temporally increasing trend of 0.025±0.006 Pg C yr -2 over 1980-2012, indicating an increasing carbon return from the roots to the atmosphere; (2) uneven temporal and 470 spatial variabilities in varying climate zones and biomes indicated their uneven temperature sensitivities to future climate change, challenging the perspective that the parameters of global carbon stimulation independent on climate zones and biomes; (3) precipitation dominated the spatial variabilities of RA.
However, further improvements in the approach should overcome shortcomings from reduced data availability and the mismatch in spatial resolution between covariates and in situ RA.

475
Author contributions. XT, SF, WZ and SG designed the research and collected the data, XT, WZ and SG contributed to the data processing and analysis. XT, WZ, SG, GC and LS wrote the manuscript, and all authors contributed to review the manuscript.
Competing interests. The authors declare that they have no conflict of interest.