Air temperature optima of vegetation productivity across global biomes.

The global distribution of the optimum air temperature for ecosystem-level gross primary productivity ([Formula: see text]) is poorly understood, despite its importance for ecosystem carbon uptake under future warming. We provide empirical evidence for the existence of such an optimum, using measurements of in situ eddy covariance and satellite-derived proxies, and report its global distribution. [Formula: see text] is consistently lower than the physiological optimum temperature of leaf-level photosynthetic capacity, which typically exceeds 30 °C. The global average [Formula: see text] is estimated to be 23 ± 6 °C, with warmer regions having higher [Formula: see text] values than colder regions. In tropical forests in particular, [Formula: see text] is close to growing-season air temperature and is projected to fall below it under all scenarios of future climate, suggesting a limited safe operating space for these ecosystems under future warming.

Understanding how photosynthesis responds to warming has been a focus in plant research in recent decades, and most of the existing knowledge comes from leaf-scale measurements[1-4]. Most leaf-scale temperature response curves show that photosynthetic capacity increases with temperature up to an optimum temperature , which typically occurs in the 30–40°C temperature range[5,6]. Above this optimum temperature, foliar photosynthetic capacity sharply declines as electron-transport and Rubisco enzymatic capacities become impaired[7]. Field et al.[8] first suggested that ecosystem-scale optimum temperature may differ from . At ecosystem scale, elevated air temperatures do limit canopy photosynthesis by other processes than leaf carboxylation rates. For instance, elevated air temperatures may accelerate leaf ageing and increase leaf thickness (phenology; e.g. ref 9), and control stomatal closure, because a higher temperature usually comes with a higher vapor pressure deficit (VPD)[10]. In a more extreme case, warming-induced water stress could suppress canopy photosynthesis though partial hydraulic failure (hydraulics) by cavitation (e.g. ref 11).

Empirical leaf-scale photosynthesis-temperature relationships[12] have been directly incorporated into global ecosystem models, with variants to account for acclimation, i.e. a temporal adjustment of optimum photosynthetic temperature to air temperature during growth[5,13,14]. This direct scaling of temperature responses from leaves to ecosystems partly determines model projections of Gross Primary Productivity (GPP) and CO2 uptake by terrestrial ecosystems in climatic scenarios. Verifying the existence of in real world ecosystems, defining its spatial distribution across and within biomes and understanding the relationships between , prevailing air temperature and are important for evaluating models and understanding the impacts of various targets of climatic warming targets on ecosystem productivity.

In this study, we formulate and test the following hypotheses: (i) is higher for biomes where air temperature during growth is warmer, (ii) is lower than for any given ecosystem, because the above mentioned limitations of stomatal conductance and phenology emerge before temperature begins to impair foliar photosynthetic capacity, and (iii) tropical forests already operate near a high above which canopy photosynthesis may decrease with even moderate air temperature warming[15,16]. Here we defined as the daytime air temperature at which GPP is highest over a period of several years, thus can be empirically determined from productivity observations and proxies (see Methods).

Results and discussion

We first applied this approach on time series of daily GPP derived from CO2 flux measurements at 153 globally distributed eddy covariance sites and found that a robust estimate of could be derived at 125 out of 153 sites (see Methods). values derived from the FLUXNET data range from 8.2°C to 35.8°C (Fig. 1a, Supplementary Table 1). Tropical sites have higher values than temperate and boreal sites (Supplementary Fig. 1), implying a dependence of on background climate. The FLUXNET multi-site analysis further indicates that across sites values are positively correlated with growing-season mean daily maximum air temperature (, see calculation in Methods) (R=0.46, P<0.01, t test), with a spatial linear regression slope of 0.61 °C per °C across sites (Fig. 1a). Overall, these results confirm our first hypothesis stating that higher values occur where higher growth temperatures prevail, in support of findings of Baldocchi et al.[17] and Niu et al.[18].

Fig. 1 |

Distribution of ecosystem-scale optimal temperature for vegetation productivity derived from flux tower sites and satellite-based data for near-infrared reflectance of vegetation (NIRV).

a, Relationship between mean annual daily maximum air temperature during the growing season and derived from daily measurements of photosynthesis across eddy-covariance sites. Flux-derived and were both obtained using observations from flux towers. Error bars indicate ±SD. The dotted gray line represents y=x and the dot line in red is y=0.61x+10.65, which is derived by linear regression with the statistical significance of the slope, or its p-value, given by Student’s t test. b, Relationship between derived from flux data and derived from NIRV data. For each site, we extracted and averaged values within a 3×3 pixel window around the site from NIRV-derived map, and calculated the SD of the nine values within the window. Error bars indicate ±SD. The dotted gray line represents y=x and the dot line in red is y=0.74x+7.10, which is derived by linear regression with the statistical significance of the slope, or its p-value, given by Student’s t test. c, Spatial distribution of for vegetation productivity (left panel), and averaged by latitude (right panel). is determined using NIRV data calculated based on satellite observations from Moderate Resolution Imaging Spectroradiometer (MODIS). Note that only gridded pixels with annual mean NDVI value larger than 0.1 and detectable are shown here. Areas of tropical forests based on current vegetation distribution are indicated by hatching. The circles on the map are colored according to the local value of retrieved from GPP at the location of each flux site. The solid line and shaded area in the right panel indicate the mean and SD, respectively, of summarized by latitude. d, in the climate space (left panel) and the temperature sensitivity of along the precipitation gradient (right panel). Each climate bin is defined by 1-°C intervals of and 100-mm intervals of mean annual precipitation, based on current climate conditions averaged over 2001–2013. The solid line in the right panel represents the temperature sensitivity of along the precipitation gradient, calculated as the slope of the linear regression between and for a given precipitation level. The shaded area indicates the SD of temperature sensitivity of estimated by bootstrapping.

Since eddy covariance measurements do not have a continuous spatial coverage, we also used satellite observations known to be highly correlated with photosynthetic activity[19], that is, GPP proxies. The first proxy used is the near-infrared reflectance of terrestrial vegetation (NIRV, the product of total-scene NIR reflectance (NIRT) by the Normalized Difference Vegetation Index (NDVI). NIRV was proven to have a high temporal correlation with GPP at flux tower sites[19]. Satellite observations of NIRT and NDVI from the Terra Moderate Resolution Imaging Spectroradiometer (MODIS) were used to calculate NIRV for the period of 2001–2013 (see Methods). NIRV-derived is comparable with that estimated from eddy covariance flux tower measurements (Fig. 1b), which gives support to using the NIRV proxy for a global mapping of . The average over the global vegetated areas is estimated to be 23±6 °C (mean ± 1 SD) with large spatial gradients in latitude. As shown in Fig. 1c, maximum values close to 30°C mainly appear over tropical forests, savannas and drylands and minimum values near 10°C prevail at high-latitudes and in mountainous regions (Fig. 1c). This spatial pattern of is robust to the choice of a particular climate forcing dataset, or to the method used to estimate (Supplementary Fig. 2, see also Methods). Similar results are also found for other GPP proxies (vegetation greenness (NDVI)[20], Enhanced Vegetation Index (EVI)[21], sun-induced vegetation fluorescence (Sun-induced Chlorophyll Fluorescence, SIF)[22]), or when daily mean air temperature is used instead of daily maximum air temperature to calculate (Supplementary Figs. 3–6, see also Methods). Note that although the covariance between air temperature, atmospheric VPD and solar radiation may confuse the direct effect of air temperature on vegetation productivity, we verified that neither VPD nor radiation is the dominant factor determining the pattern of at the global scale (see Methods).

In order to test the second hypothesis, we compared satellite-derived with from the responses of maximum Rubisco-limited carboxylation rates (Vcmax) to temperature from leaf-scale measurements for 36 species[5]. Note that the here refers to the temperature optima for leaf-scale (gross) photosynthetic capacity rather than for leaf net photosynthesis, which equals gross photosynthesis minus photorespiration and minus dark respiration (see more details in Methods). We found that is indeed lower than (Supplementary Fig. 7). This difference may originate from that is additionally limited by high VPD during hot and dry periods[6] and by soil-moisture deficits during extensive dry episodes[23], under real world conditions. Under conditions of high temperature, atmospheric VPD increases while soil moisture decreases. Stomatal conductance and hence carbon assimilation rates (GPP at ecosystem-scale) decrease to prevent exceedingly low leaf water potentials and resulting plant tissue damage from cavitation[24]. In contrast, leaf-level photosynthesis measurements that determine the temperature response curve of Vcmax are usually performed in absence of water stress through maintaining relatively low VPD conditions (e.g. ref 25–30) unless the research objective is to investigate drought effect on leaf photosynthetic parameters as the studies by Vaz et al.[31] and Zhou et al.[32]. In addition, plant phenology controls leaf age, vitality (photosynthetic rates) and foliar density (e.g. Leaf Area Index, LAI)[33], and may therefore co-determine ecosystem-level temperature limitations and the optimum temperature for canopy photosynthesis[34]. It is also important to remark when comparing with that leaf-scale measurements are often limited to sunlit leaves, which could lead to a positive bias of existing in-situ measurements. Furthermore, the tree species database used by Kattge & Knorr[5] from which data were collected does not include any tropical species. This may explain why global models prescribed with give divergent results for tropical biomes.

The relationship between and background climate is shown in Fig. 1d. The sampling of leaf-scale studies does not provide consistent evidence about the dependence of on climate, with positive correlations between and growing season air temperature in a set of studies[1, 5, 35–37] attributed to evolutionary adaptation[38], but no clear relationship between and growth temperature[39-41]. In contrast, inferred from satellite GPP proxies in our study increases with across the globe. In temperature-precipitation space, the spatial sensitivity of to (the slope of the linear regression between these two variables) is lower than 1 for any precipitation bin (Fig. 1d), suggesting that spatial gradients of are smaller than those of , possibly because hydraulic and phenological limitations further limit across spatial gradients. In fact, the spatial sensitivity of to generally increases with increasing mean annual precipitation (Fig. 1d), even though is not significantly correlated with precipitation after controlling for the effect of (Fig. 1d). This thermal adaptation of , suggested by the positive spatial slope of the -air temperature relation, is also observed across biomes. As shown in Fig. 2, there is a significant positive correlation between and with a slope of 0.76 across different biomes. Among biomes, the largest mean is found in tropical evergreen broad-leaved forest (29±3 °C), and the smallest mean (13±3 °C) in cold grasslands covering the Tibetan Plateau (Fig. 2 and Supplementary Fig. 8).

Fig. 2 |

Relationship between mean annual daily maximum air temperature during the growing season and ecosystem-scale optimum temperature for vegetation productivity across vegetation types.

The error bars indicate the SDs of for each vegetation type: ENF, evergreen needle-leaved forest; EBF, evergreen broad-leaved forest; DNF, deciduous needle-leaved forest; DBF, deciduous broad-leaved forest; MF, mixed forest; Shrub, closed and open shrublands. The light-gray dotted line represents y=x. The dark-gray dotted line is y=0.76x+6.48 derived by linear regression with the slope value (estimated using Student’s t test) shown in the bottom right. The red dotted line is the flux tower derived slope (0.61) from Fig. 1a. The size of each symbol corresponds to the three categories (< 3%, 3%−10% and > 10%) of occupied vegetated area on land. Error bars indicate ±SD.

Both model results and very limited observational studies suggest a decrease in canopy photosynthesis of tropical forests at high temperature[15, 42–45] which led us to formulate the third hypothesis of tropical forests already operating at being close to implying that canopy photosynthesis may decrease in case of future warming[15, 16]. This hypothesis is verified from the data shown in Fig. 3 (see also Supplementary Fig. 9). is indeed slightly (1.4 °C) lower than over tropical evergreen forests, suggesting a small ‘safety margin’ for canopy photosynthesis under future warming. Note that the “safety margin” could become larger than that suggested by the air temperature data if leaf thermal regulation would acclimate to the warming air temperature (see Methods). In contrast, arctic (north of 65°N) and boreal (50°N-65°N) ecosystems exhibit substantially larger safety margins, i.e. a larger positive difference between and (Fig. 3a and Supplementary Fig. 9). Analysis of the 16-day averaged distribution during the period when is observed further shows that the rank of in the distribution is already near the highest quantile of (>70%) for tropical evergreen forests (Supplementary Fig. 10). Based on this result, one may expect that rising air temperature in the future, irrespective of the indirect effect of increasing VPD, may limit or decrease vegetation productivity in tropical forests, but not in temperate or boreal ecosystems.

Fig. 3 |

Change with latitude in ecosystem-scale optimal temperature for vegetation productivity and daily maximum air temperature averaged over the growing season .

a, Current versus current ; b, Future versus future . Current and are calculated using current temperature for 2001–2013, whereas acclimated and future are first calculated pixel by pixel using temperature for 2091–2100 projected by General Circulation Models (GCMs) under the RCP4.5 scenario and then averaged by latitude. Acclimated is determined based on the projected temperature and temperature sensitivity of using the annual precipitation level predicted for 2091–2100. The solid line and shaded area in each panel indicate the mean and SD, respectively, of or summarized by latitude. c, Future versus future for tropical evergreen forests. ** indicates that is significantly lower than at P<0.01 in a paired t-test. Error bars indicate ±SD.

Global terrestrial daily maximum air temperature is projected to rise by 1.9 °C under the RCP2.6 low warming climate scenario and by 5.6 °C under the RCP8.5 scenario by 2100[46]. We compared these projections with the present-day distribution of with a focus on tropical evergreen forests, where optimum temperature is currently just below the limit of (see Methods; Fig. 3b and 3c). The key uncertainty in this discussion is, however, whether or not will acclimate and follow the increase in . We therefore looked at possible acclimation from time series of retrieved from the Advanced Very High Resolution Radiometer (AVHRR) NDVI, which spans the last 30 years and comprises almost a 1°C temperature range. NDVI-derived did not have a significant trend over the last three decades except for the northern lands (north of 60°N) where warming is more pronounced[47] (Supplementary Fig. 11). This suggests that the recent 1°C warming is not large enough to elicit an acclimation response from some ecosystems, given decadal variability[48]. Also the annual derived from flux sites estimates of GPP did not exhibit a positive trend and was not significantly correlated with annual variations of , although the flux time series are probably too short to properly evaluate trends of related to possible acclimation processes (Supplementary Fig. 12). Because we detected no indication for its existence, we first assumed no acclimation in the comparison of future projections from climate models with the current distribution of . Under this assumption, the average of tropical evergreen forests will exceed the current value of for RCP2.6 by 2.6°C, and by 5.7 °C for RCP8.5 (Fig. 3c). On the other hand, boreal and arctic biomes will still remain within the safety margin, with staying above except under the RCP8.5 high warming scenario (Fig. 3b and Supplementary Fig. 13).

Despite the lack of in situ observational evidence for GPP acclimation to the ongoing warming trend, we tested a simple future acclimation scenario based on the space-for-time substitution approach[49] as applied in several studies using observed spatial gradients to hindcast temporal changes[50,51]. Here, we assume that temporal change of will evolve proportionally to , following the spatial temperature sensitivity of to in Fig. 1d and the indirect effects of temperature increase (e.g. by increasing VPD) are excluded. We took the differences in precipitation levels into account, so that areas that become wetter also exhibit faster acclimation. Even with this assumed acclimation law, will still surpass by 1.7 °C under RCP2.6 and by 2.5 °C under RCP8.5 for tropical evergreen forests (Fig. 3c). Not accounting for precipitation levels in the acclimation rates produced similar results (Supplementary Figs. 14 and 15).

Our global-scale analysis of derived from globally distributed point measurements of eddy covariance and space-borne observations of proxies of vegetation productivity is a first attempt to diagnose the global distribution of ecosystem-scale temperature optima of photosynthesis. It should be noted, however, that hypotheses about that thermal acclimation of are still highly uncertain, because ecosystem adjustments can lag substantially behind the rate of future warming, particularly for forests. More studies using data sets with longer time spans are needed in the future in order to more accurately detect eventual thermal acclimation of . Furthermore, the acclimation of plants to increasing atmospheric CO2 concentration and to changes in other environmental factors (e.g. VPD) was also not considered in the current analyses. Constraining the spatially observed temperature sensitivity of over time is a priority for future studies. Continuous monitoring and dedicated manipulative experiments could improve our understanding of the features of and thermal acclimation in earth system models[52].

Methods

FLUXNET data

The half-hourly eddy-covariance Gross Primary Productivity (GPP) data were obtained from FLUXNET datasets, and were quality-controlled, filtered against low turbulence, and gap-filled using consistent methods, as described by Papale et al.[53]. Only freely available FLUXNET data were used in this study. All the half-hourly GPP data were aggregated into daily accumulated GPP for further estimates of the optimal temperature for vegetation productivity. Daily maximum air temperature was determined as the maximum air temperature value from all the half-hourly air temperature observations. We included only site-years with more than 80% of half-hourly data available. A total of 153 individual FLUXNET sites with 663 site-years of GPP data were used in this study.

Near-infrared reflectance of terrestrial vegetation

Badgley et al.[19] have recently proposed a new approach for estimating vegetation photosynthetic capacity by remote sensing, i.e. the near-infrared reflectance (NIR) of terrestrial vegetation (NIRV), which can differentiate between the confounding effects of background brightness, leaf area and the distribution of photosynthetic capacity with depth in canopies[19]. NIRV is calculated as the product of total scene NIR reflectance (NIRT) and Normalized Difference Vegetation Index (NDVI)[19]. As a proxy of photosynthesis, NIRV is suggested to be strongly correlated with Solar-induced Chlorophyll Fluorescence (SIF), a direct index of photons intercepted by chlorophyll, and shows higher correlation with observed GPP than NDVI[19]. We used satellite-derived NIRV to calculate and map the optimal air temperature for vegetation productivity at an ecosystem scale . Following Badgley et al.[19], we calculated 16-day NIRV for 2001–2013 as the product of Moderate Resolution Imaging Spectroradiometer (MODIS) 16-day NIR reflectance and MODIS 16-day NDVI, both of which were derived from the MOD13A2 Vegetation Index Product with a spatial revolution of 1 km. Only positive NIRV values were used in the analysis.

Normalized Difference Vegetation Index

The Normalized Difference Vegetation Index (NDVI) is a vegetation index defined as the ratio of the difference between NIR and red visible reflectance to their sum, and is widely used to represent vegetation greenness[54]. To account for uncertainties from different satellite datasets, three independent NDVI datasets were utilized, including biweekly NDVI data from Global Inventory Modeling and Mapping Studies (GIMMS) Advanced Very High Resolution Radiometer (AVHRR), 16-day NDVI data from Terra Moderate Resolution Imaging Spectroradiometer (MODIS) and 10-day NDVI data from Satellite Pour l’Observation de la Terre Vegetation (SPOT Vegetation). The three NDVI datasets spanned the last three decades: 1982–2009 for AVHRR NDVI datasets, 2000–2009 for MODIS NDVI datasets, and 1999–2009 for SPOT NDVI datasets, with the spatial resolutions of 8 km, 1 km, and 1 km, respectively. All NDVI datasets have been corrected to reduce the effects of volcanic aerosols, solar angle, and sensor errors[20,55,56]. Pixels with a mean annual NDVI>0.1 were defined as the vegetated area for each dataset.

Enhanced Vegetation Index

The Enhanced Vegetation Index (EVI) is another vegetation index designed to enhance the vegetation signal by minimizing canopy-soil variations and to improve sensitivity over dense vegetation conditions[21], and is found to correlate well with estimated GPP on a site-by-site basis[57]. We used a 16-day EVI dataset for the period of 2000–2009 with a spatial resolution of 1 km from the MOD12A1 Vegetation Index Product. Effects from aerosols, solar angle and sensor error have all been corrected[21].

Sun-induced Chlorophyll Fluorescence

Chlorophylls in plant absorbs shortwave radiation and dissipates excess energy as light or heat. The longwave radiation re-emitted by chlorophylls is referred as chlorophyll fluorescence. Recent studies have reported that remotely sensed Sun-induced Chlorophyll Fluorescence (SIF) could serve as an indicator of photosynthesis rate and it is well correlated with model-simulated GPP[58]. Following the previous studies[58,59], we retrieved SIF from two different retrieval windows, 757 nm and 771 nm, as well as the two polarization states, S and P using observation of Fourier Transform Spectrometer (FTS) on the Japanese Greenhouse gases Observing SATellite (GOSAT)[20]. These diverse SIF samples of were then aggregated into monthly gridded data at a spatial resolution of 2° from June 2009 to June 2012.

Vegetation distribution

We used MODIS land cover with the classification scheme of the International Geosphere-Biosphere Programme (IGBP). The MODIS IGBP land cover data were derived from the MOD12Q1 Land Cover Science Data Product at a spatial resolution of 1 km and an updated digital Köppen-Geiger world map of climatic classification[60]. Within the vegetated area defined by NDVI thresholds, the 17 land cover types were reclassified into 9 vegetation types, namely evergreen needle-leaved forest (ENF), evergreen broad-leaved forest (EBF), deciduous needle-leaved forest (DNF), deciduous broad-leaved forest (DBF), mixed forest (MF), savannas, cropland, grassland and shrubland. Based on the main climates in the world Map of the Köppen-Geiger climatic classification[60], grassland was further subdivided into temperate grasslands, boreal and arctic tundra, and shrubland was further subdivided into temperate and boreal shrubland. The grassland over the Tibetan Plateau was considered separately considering the fact that Tibetan Plateau has an average altitude higher than 4000 m a.s.l.[61], and thus an unique alpine climate. In contrast to temperate grasslands / shrubland where water is major limiting factor for vegetation productivity, alpine ecosystems on the Tibetan Plateau are mainly limited by thermal conditions[62].

Climate dataset

The gridded air temperature and precipitation data for the period of 1982–2013 were obtained from the CRU/NCEP 6-hourly dataset with a spatial resolution of 0.5°. Note that the purpose of this study is to investigate the optimal air temperature for photosynthesis. Optimal leaf temperature is also interesting but not addressed in this study because accurate canopy-integrated measurements of leaf temperatures are not available at the eddy covariance sites and at global scale as gridded datasets. For a discussion about calculation of temperature optimum from air temperature and from surface temperature, we used the remotely sensed land surface temperature (LST), which is inversed from infra-red emissivity measured by MODIS (MYD11A2 version 6). This dataset had an original spatial resolution of 1 km, spanning from July 2002 to December 2014. The error of MODIS LST product, which primarily stems from cloud contamination and emissivity uncertainties, was reported to be less than 3°C[63]. Generally, the occurrence time of (2:00–4:00 P.M.) is relatively close to the Aqua overpass time (1:30 P.M.); thus, we assumed that from MODIS-Aqua is comparable with the daily maximum leaf surface temperature . Corresponding to the temporal resolutions of MODIS, AVHRR and SPOT datasets, respectively, the 6-hourly climate data were aggregated into 16-day, biweekly, and 10-day values before further analyses. Given the different spatial resolutions of satellite observations and climate data, we simply extracted time series of daily maximum air temperature and precipitation from the aggregated CRU/NCEP data for each pixel of the sets of remotely sensed data. The daily maximum air temperature of the growing season averaged over 2001–2013 was calculated as the current mean growing-season daily maximum air temperature . Information on the growing season was derived from the study by Zhu et al.[64], which was first determined from the GIMMS Leaf Area Index data set (GIMMS LAI3g) using a Savitzky-Golay filter and then refined by excluding the ground-freeze period identified by the Freeze/Thaw Earth System Data Record (see details in ref 64). We also documented the temperature thresholds at which the growing season begins and ends for each year. Temperature thresholds were averaged over 2001–2013 for the onset and end of the growing season, respectively. We also used WATCH Forcing Data Methodology to ERA-Interim data with a temporal resolution of three hours (WFDEI)[65].

We also used climate projections for the end of the 21st century (2091–2100) using 20 models that participated in the phase five of Coupled Model Intercomparison Project (CMIP5) under the RCP2.6, RCP4.5 and RCP8.5 scenarios[46] to determine the impact of future warming on vegetation productivity (see model list in Supplementary Table 2). Considering the mismatch between CRU/NCEP datasets and outputs from General Circulation Models (GCMs) for current climate conditions, we generated future temperature and precipitation maps by adding the relative changes in GCM-derived climate projections to the current climate for each pixel. for the late 21st century was estimated using the same temperature thresholds as for the current . All GCM projections were resampled to a resolution of 1° using first-order conservative interpolation method[66].

Analysis

We estimated local by examining the temperature response curve of MODIS NIRV. Following Yuan et al.[37] and Niu et al.[18], NIRV time series throughout the entire monitoring period and the corresponding temperature data were grouped into 1°C-temperature bins for each pixel within vegetated areas, which were defined as regions with a mean annual NDVI value larger than 0.1. We used the 90% quantile of the NIRV data as the response of NIRV within each temperature bin due to the potential influences of other environmental constraints like clouds and droughts. We next calculated the running means of every three temperature bins to develop the temperature response curve of NIRV. The was then determined from the response curve at which NIRV was maximized (Supplementary Fig. 16). Note that may not be detected for some pixels where the maximum NIRV was only attained at either end of the response curve, accounting for 3.5% of the vegetated areas. Only vegetated areas with detectable were shown when mapping the spatial pattern of . The derivation of is robust to the choice of a particular climate-forcing dataset (Supplementary Fig. 2). Instead of using the temperature corresponding to the maximum 90th quantile NIRV to calculate , we also applied nonlinear regression of the photosynthetic temperature response data (Eq.1) to estimate , which produced similar results (Supplementary Fig. 2): where NIRV(T) is the NIRV value at a daily maximum temperature T, b is a parameter describing the spread of the parabola[48,67]. is the vertex of each fit and NIRV(OPT) is the NIRV value at . Finally, we used daily mean air temperature instead of to calculate In this test, derived from is smaller than estimated from , but the two variables were strongly spatially correlated (Supplementary Fig. 6).

We investigated the relationship between and climate variables by averaging in the climate space with 1-°C intervals of mean annual averaged over the growing season and 100-mm intervals of mean annual precipitation (MAP) (Fig. 1d). For each MAP interval, we calculated the ‘apparent’ spatial sensitivity of in response to changes in using bootstrapping method. We performed the linear regression analysis 1000 times by randomly selecting a subset of 80% of the samples from pairs of and within each MAP interval. The mean and SD of the 1000 temperature sensitivity of were subsequently estimated along the MAP gradient.

Air temperature, atmospheric vapor pressure deficit (VPD) and solar radiation usually co-vary in time and space, so that the empirical observation of spatial patterns of in this study cannot be unambiguously attributed to air temperature as a single explaining factor of . Under conditions of high temperature, atmospheric VPD increases, soil moisture decreases with a lag, stomatal conductance and hence carbon assimilation rates (GPP at the ecosystem-scale) decrease to prevent exceedingly low leaf water potentials and resulting plant tissue damage from cavitation[24]. We show that across climatic gradients is systematically higher at high maximum air temperatures but not systematically lower at high VPD conditions (Supplementary Fig. 17). Then we calculated the variance inflation factor (VIF) between VPD and under each VPD bin in the regression model of: As shown in Supplementary Figure 18, we observed that VIF value ranged only between 1.001 and 1.438, suggesting relatively low multicollinearity between VPD and temperature. Even so, to examine whether VPD can substantially affect the relationship between and , we further calculated the partial (‘intrinsic’) sensitivity of to in each grid point based on the following bilinear regression: where the partial sensitivity of to is defined as in Eq. 3 under each VPD bin. Then we compared the partial sensitivity with the apparent sensitivity of to estimated using abovementioned linear regression between and for each VPD bin. As shown in Supplementary Fig. 19, although the apparent sensitivity of to is generally lower than the partial (‘intrinsic’) sensitivity of to the apparent sensitivity to remains positive even when VPD is taken into account excepted under very high VPD bins (higher than ~4.5 kPa) representing less than 1% of the study area. These results indicate that the patterns of are not dominated by high VPD reducing canopy photosynthesis, as an indirect effect of higher air temperature increasing VPD. Moreover, we also calculated the percentiles of downward shortwave solar radiation (Rad) at the time of year when is observed for the 16-day averaged Rad distribution. As shown in Supplementary Fig. 20, the Rad value when was retrieved from global observations were below the 95th percentile in the 16-day Rad distribution for ~80% of the study area, which is mainly in mid and low latitudes such as Africa, India, Australia, eastern Brazil and southern and southwestern of North America. By comparison, for most part of boreal regions, part of south China, southeastern US, as well as part of South America, the timing of is consistent the time of maximum solar radiation. This is because in these regions generally appears in summer, which is also the period when solar radiation displays its maximum during the year.

The NIRV-derived was compared with estimated using GPP data from 153 eddy covariance sites. Flux-derived was determined for each site-year with daily-accumulated GPP and corresponding temperature data from flux tower observations, applying the same method to estimate local using NIRV datasets. A robust estimate of can be derived for 125 sites (Supplementary Table 1). For each site, we calculated the mean and SD of across different years. Then, we extracted and averaged values within a 3×3 pixel window around each site from the NIRV-derived map, and calculated the SD of the nine values within the window. The relationship between NIRV- and flux-derived was reported using a least square linear regression, and the statistical significance of the slope, or its p-value, given by Student’s t test. The results show that NIRV-derived is comparable with that estimated independently from measurements of flux-tower eddy covariance (Fig. 1b).

We compared the spatial distribution of derived from NIRV with the one obtained from NDVI datasets. Consistent spatial patterns of are derived from each of the three NDVI datasets (Supplementary Fig. 21). A global composite map of (Supplementary Fig. 3) was then generated by averaging over estimates derived from the three NDVI datasets. Given the inconsistent spatial resolutions of the different products, we resampled to a common grid of 8 km before averaging. from NDVI datasets generally show a spatial pattern similar to that from NIRV, but with smaller NDVI-derived values for central Australia and southern South America (Supplementary Fig. 3). We compared the spatial distribution of derived from NIRV with that from MODIS EVI data during 2001–2013, and found that the EVI derived shows very similar spatial pattern to that of NIRV derived (Supplementary Fig. 4). The distribution of derived from NIRV and from GOSAT SIF datasets also have similar spatial patterns, even though the NIRV-derived is higher in tropical regions, particularly in cultivated areas of southeastern Brazil (Supplementary Fig. 5).

At leaf scale, the photosynthesis-temperature response is suggested to be primarily controlled by three sets of processes, namely biochemical, respiratory and stomatal processes[68]. Much of the effort to date to understand variability in the leaf-level photosynthesis-temperature response has focused on biochemical processes[68], with Vcmax and Jmax being two major parameters controlling the maximum rates of photosynthesis limited by CO2 and light, respectively[69]. Therefore, in this study, we compared derived using GPP proxies with leaf-scale optimal temperature of maximum Rubisco-limited carboxylation rates (Vcmax), although GPP is in theory more comparable to net photosynthesis, that is, leaf gross photosynthesis minus photorespiration and minus dark respiration. Since photorespiration increases exponentially with temperature[70], the optimum temperature of GPP should, in theory, be lower than the optimal temperature of maximum Rubisco-limited carboxylation rates (Vcmax). For this comparison to be made, we extracted and averaged values within a 3×3 pixel window from the NIRV-derived map around the reported site location (longitude and latitude) of leaf-scale measurements. For leaf-scale measurements without the information of site location, we calculated the average NIRV-derived values across pixels with both the same growing season mean temperature and the same plant functional type as the corresponding site.

is different from not only because of respiratory process, but also because air temperature can differ from leaf temperatures[71], which are regulated by leaf traits affecting the leaf energy balance[72]. Because, to our knowledge, global gridded monthly leaf temperature data are not available, we use daily maximum land surface temperature from MODIS to calculate in order to illustrate the potential differences between and . As shown in Supplementary Figure 22, the is similar to over tropical savannas. However, over moist tropical forests is lower than , which can be explained by the lower daytime surface temperature than air temperature as a result of strong evapotranspiration effects[71,73]. This ecosystem-dependent difference between and suggests that the leaf thermal regulation mechanism through the physiological and morphological changes[72] is an important ecosystem process shaping spatial variations of . In addition, if the difference between leaf temperature and air temperature would become larger in response to warmer air temperatures (i.e. if leaf thermal regulation acclimates to warmer temperature), the “safety margin” of tropical ecosystems would become larger than the air temperature data are currently suggesting. However, the long-term in-situ leaf temperature data required to test this hypothesis independently are currently not available.

To account for potential changes in under future warming, we estimated the acclimated for vegetation productivity by the end of the 21st century (2091–2100), using recent IPCC climate projections[46]. To this end, we applied the space-for-time substitution approach[49], assuming that that temporally will evolve proportionally to following the spatial temperature sensitivity of to . Given the relatively large uncertainties of precipitation projections, we considered two future precipitation scenarios. For the first scenario, we estimated acclimated pixel by pixel using the temperature sensitivity of under the present MAP level, assuming that MAP does not change between the end of the 21st century. For the second scenario, we accounted for in MAP, and the acclimated was calculated pixel by pixel using the temperature sensitivity of under the projected MAP level for 2091–2100. Acclimated was averaged across the GCMs under each scenario. Latitudinal variation of future was derived by averaging within 1°-latitude bins from future maps and then compared with that in future summarized by latitude from future maps.