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

HAL
open science

# Air temperature optima of vegetation productivity across global biomes.

Mengtian Huang, Shilong Piao, Philippe Ciais, Josep Peñuelas, Xuhui Wang, Trevor F Keenan, Shushi Peng, Joseph A Berry, Kai Wang, Jiafu Mao, et al.

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

## To cite this version:

Mengtian Huang, Shilong Piao, Philippe Ciais, Josep Peñuelas, Xuhui Wang, et al., Air temperature optima of vegetation productivity across global biomes., Nature Ecology & Evolution, 2019, 3 (5), pp. 772-779. 10.1038/s41559-019-0838-x. hal-02181076

HAL Id: hal-02181076

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

Submitted on 2 Oct 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_4.png" /></figure>

# Air temperature optima of vegetation productivity across global biomes

Mengtian Huang¹, Shilong Piao¹,²,³,*, Philippe Ciais⁴, Josep Peñuelas⁵,⁶, Xuhui Wang¹, Trevor F. Keenan⁷,⁸, Shushi Peng¹, Joseph A. Berry⁹, Kai Wang¹, Jiafu Mao¹⁰, Ramdane Alkama¹¹, Alessandro Cescatti¹¹, Matthias Cuntz¹², Hannes De Deurwaerder¹³, Mengdi Gao¹, Yue He¹, Yongwen Liu¹, Yiqi Luo¹⁴, Ranga B. Myneni¹⁵, Shuli Niu¹⁶, Xiaoying Shi¹⁰, Wenping Yuan¹⁷, Hans Verbeeck¹³, Tao Wang²,³, Jin Wu¹⁸, and Ivan A. Janssens¹⁹

¹Sino-French Institute for Earth System Science, College of Urban and Environmental Sciences, Peking University, Beijing 100871, China.

²Key Laboratory of Alpine Ecology and Biodiversity, Institute of Tibetan Plateau Research, Chinese Academy of Sciences, Beijing 100085, China.

³Center for Excellence in Tibetan Earth Science, Chinese Academy of Sciences, Beijing 100085, China.

⁴Laboratoire des Sciences du Climat et de l'Environnement, CEA CNRS UVSQ, Gif-sur-Yvette 91191, France.

⁵CREAF, Cerdanyola del Vallès, Barcelona 08193, Catalonia, Spain.

⁶CSIC, Global Ecology Unit CREAF-CSIC-UAB, Bellaterra, Barcelona 08193, Catalonia, Spain.

⁷Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA.

⁸Department of Environmental Science, Policy and Management, UC Berkeley, Berkeley, California 94720, USA.

⁹Department of Global Ecology, Carnegie Institution for Science, Stanford, California 94305, USA.

¹⁰Climate Change Science Institute and Environmental Sciences Division, Oak Ridge National Laboratory, Oak Ridge, TN, USA.

¹¹European Commission, Joint Research Centre, Institute for Environment and Sustainability, Ispra, Italy.

¹²INRA, Université de Lorraine, AgroParisTech, UMR Silva, 54000 Nancy, France.

Users may view, print, copy, and download text and data mine the content in such documents, for the purposes of academic research, subject always to the full Conditions of use: http://www.nature.com/authors/editorial_policies/license.html#terms

*Correspondence to: slpiao@pku.edu.cn.

**Author contributions:** S.L.P. designed the research; M.T.H. performed the analysis; S.L.P. drafted the paper; and all authors contributed to the interpretation of the results and to the text.

## Data Availability

All data is available in the main text or the supplementary information. All computer codes used in this study can be provided by the corresponding author upon reasonable requests.

**Competing interests:** Authors declare no competing interests.
¹³CAVElab Computational & Applied Vegetation Ecology, Faculty of Bioscience Engineering, Ghent University, Coupure Links 653, B-9000 Gent, Belgium.

¹⁴Department of Biological Sciences, Northern Arizona University, Arizona, 86011, USA.

¹⁵Department of Earth and Environment, Boston University, Boston, Massachusetts 02215, USA.

¹⁶Key Laboratory of Ecosystem Network Observation and Modeling, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing, China.

¹⁷School of Atmospheric Sciences, Center for Monsoon and Environment Research, Sun Yatsen University, Guangzhou 510275, China.

¹⁸Environmental & Climate Sciences Department, Brookhaven National Laboratory, Upton, NY, 11973, USA.

¹⁹Centre of Excellence PLECO (Plant and Vegetation Ecology), Department of Biology, University of Antwerp, Universiteitsplein 1, B-2610 Wilrijk, Belgium.

Abstract

The global distribution of the optimum air temperature for ecosystem-level gross primary productivity ($T_{opt}^{eco}$) 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. $T_{opt}^{eco}$ is consistently lower than the physiological optimum temperature of leaf-level photosynthetic capacity, which typically exceeds 30 °C. The global average $T_{opt}^{eco}$ is estimated to be 23±6 °C, with warmer regions having higher $T_{opt}^{eco}$ values than colder regions. In tropical forests, particularly, $T_{opt}^{eco}$ 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¹⁻⁴. Most leaf-scale temperature response curves show that photosynthetic capacity increases with temperature up to an optimum temperature ($T_{opt}^{leaf}$), which typically occurs in the 30–40°C temperature range⁵,⁶. Above this optimum temperature, foliar photosynthetic capacity sharply declines as electron-transport and Rubisco enzymatic capacities become impaired⁷. Field et al.⁸ first suggested that ecosystem-scale optimum temperature ($T_{opt}^{eco}$) may differ from $T_{opt}^{leaf}$. 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)¹⁰. 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¹² 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⁵,¹³,¹⁴. This direct scaling of temperature responses from leaves to ecosystems partly determines model projections of Gross Primary Productivity (GPP) and CO₂ uptake by terrestrial ecosystems in climatic scenarios. Verifying the existence of T<sub>opt</sub><sup>eco</sup> in real world ecosystems, defining its spatial distribution across and within biomes and understanding the relationships between T<sub>opt</sub><sup>eco</sup>, prevailing air temperature and T<sub>opt</sub><sup>leaf</sup> 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) T<sub>opt</sub><sup>eco</sup> is higher for biomes where air temperature during growth is warmer, (ii) T<sub>opt</sub><sup>eco</sup> is lower than T<sub>opt</sub><sup>leaf</sup> 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 T<sub>opt</sub><sup>eco</sup> above which canopy photosynthesis may decrease with even moderate air temperature warming<sup>15,16</sup>. Here we defined T<sub>opt</sub><sup>eco</sup> as the daytime air temperature at which GPP is highest over a period of several years, thus T<sub>opt</sub><sup>eco</sup> 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 CO₂ flux measurements at 153 globally distributed eddy covariance sites and found that a robust estimate of T<sub>opt</sub><sup>eco</sup> could be derived at 125 out of 153 sites (see Methods). T<sub>opt</sub><sup>eco</sup> values derived from the FLUXNET data range from 8.2°C to 35.8°C (Fig. 1a, Supplementary Table 1). Tropical sites have higher T<sub>opt</sub><sup>eco</sup> values than temperate and boreal sites (Supplementary Fig. 1), implying a dependence of T<sub>opt</sub><sup>eco</sup> on background climate. The FLUXNET multi-site analysis further indicates that across sites T<sub>opt</sub><sup>eco</sup> values are positively correlated with growing-season mean daily maximum air temperature (T<sub>max gs</sub><sup>air</sup>, 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 T<sub>opt</sub><sup>eco</sup> values occur where higher growth temperatures prevail, in support of findings of Baldocchi et al.<sup>17</sup> and Niu et al.<sup>18</sup>.

Since eddy covariance measurements do not have a continuous spatial coverage, we also used satellite observations known to be highly correlated with photosynthetic activity¹⁹, that is, GPP proxies. The first proxy used is the near-infrared reflectance of terrestrial vegetation (NIR<sub>V</sub>, the product of total-scene NIR reflectance (NIR<sub>T</sub>) by the Normalized Difference Vegetation Index (NDVI). NIR<sub>V</sub> was proven to have a high temporal correlation with GPP at
flux tower sites¹⁹. Satellite observations of NIR<sub>T</sub> and NDVI from the Terra Moderate Resolution Imaging Spectroradiometer (MODIS) were used to calculate NIR<sub>V</sub> for the period of 2001–2013 (see Methods). NIR<sub>V</sub>-derived T<sub>opt</sub><sup>eco</sup> is comparable with that estimated from eddy covariance flux tower measurements (Fig. 1b), which gives support to using the NIR<sub>V</sub> proxy for a global mapping of T<sub>opt</sub><sup>eco</sup>. The average T<sub>opt</sub><sup>eco</sup> 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 T<sub>opt</sub><sup>eco</sup> is robust to the choice of a particular climate forcing dataset, or to the method used to estimate T<sub>opt</sub><sup>eco</sup> (Supplementary Fig. 2, see also Methods). Similar results are also found for other GPP proxies (vegetation greenness (NDVI)²⁰, Enhanced Vegetation Index (EVI)²¹, sun-induced vegetation fluorescence (Sun-induced Chlorophyll Fluorescence, SIF)²²), or when daily mean air temperature (T<sub>mean</sub><sup>air</sup>) is used instead of daily maximum air temperature (T<sub>max</sub><sup>air</sup>) to calculate T<sub>opt</sub><sup>eco</sup> (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 T<sub>opt</sub><sup>eco</sup> at the global scale (see Methods).

In order to test the second hypothesis, we compared satellite-derived $T_{opt}^{eco}$ with $T_{opt}^{leaf}$ from the responses of maximum Rubisco-limited carboxylation rates ($V_{cmax}$) to temperature from leaf-scale measurements for 36 species⁵. Note that the $T_{opt}^{leaf}$ 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 $T_{opt}^{eco}$ is indeed lower than $T_{opt}^{leaf}$ (Supplementary Fig. 7). This difference may originate from that $T_{opt}^{eco}$ is additionally limited by high VPD during hot and dry periods⁶ and by soil-moisture deficits during extensive dry episodes²³, 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²⁴. In contrast, leaf-level photosynthesis measurements that determine the temperature response curve of $V_{cmax}$ 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.³¹ and Zhou et al.³². In addition, plant phenology controls leaf age, vitality (photosynthetic rates) and foliar density (e.g. Leaf Area Index, LAI)³³, and may therefore co-determine ecosystem-level temperature limitations and the optimum temperature for canopy photosynthesis³⁴. It is also important to remark when comparing $T_{opt}^{leaf}$ with $T_{opt}^{eco}$ that leaf-scale measurements are often limited to sunlit leaves, which could lead to a positive bias of existing in-situ $T_{opt}^{leaf}$ measurements. Furthermore, the
tree species database used by Kattge & Knorr⁵ from which $T_{opt}^{leaf}$ data were collected does not include any tropical species. This may explain why global models prescribed with $T_{opt}^{leaf}$ give divergent results for tropical biomes.

The relationship between $T_{opt}^{eco}$ and background climate is shown in Fig. 1d. The sampling of leaf-scale studies does not provide consistent evidence about the dependence of $T_{opt}^{leaf}$ on climate, with positive correlations between $T_{opt}^{leaf}$ and growing season air temperature in a set of studies¹, ⁵, ³⁵⁻³⁷ attributed to evolutionary adaptation³⁸, but no clear relationship between $T_{opt}^{leaf}$ and growth temperature³⁹⁻⁴¹. In contrast, $T_{opt}^{eco}$ inferred from satellite GPP proxies in our study increases with $T_{max gs}^{air}$ across the globe. In temperature-precipitation space, the spatial sensitivity of $T_{opt}^{eco}$ to $T_{max gs}^{air}$ (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 $T_{opt}^{eco}$ are smaller than those of $T_{max gs}^{air}$, possibly because hydraulic and phenological limitations further limit $T_{opt}^{eco}$ across spatial gradients. In fact, the spatial sensitivity of $T_{opt}^{eco}$ to $T_{max gs}^{air}$ generally increases with increasing mean annual precipitation (Fig. 1d), even though $T_{opt}^{eco}$ is not significantly correlated with precipitation after controlling for the effect of $T_{max gs}^{air}$ (Fig. 1d). This thermal adaptation of $T_{opt}^{eco}$, suggested by the positive spatial slope of the $T_{opt}^{eco}$-air temperature relation, is also observed across biomes. As shown in Fig. 2, there is a significant positive correlation between $T_{opt}^{eco}$ and $T_{max gs}^{air}$ with a slope of 0.76 across different biomes. Among biomes, the largest mean $T_{opt}^{eco}$ is found in tropical evergreen broad-leaved forest (29±3 °C), and the smallest mean $T_{opt}^{eco}$ (13±3 °C) in cold grasslands covering the Tibetan Plateau (Fig. 2 and Supplementary Fig. 8).

Both model results and very limited observational studies suggest a decrease in canopy photosynthesis of tropical forests at high temperature¹⁵, ⁴²⁻⁴⁵ which led us to formulate the third hypothesis of tropical forests already operating at $T_{opt}^{eco}$ being close to $T_{max gs}^{air}$ implying that canopy photosynthesis may decrease in case of future warming¹⁵, ¹⁶. This hypothesis is verified from the data shown in Fig. 3 (see also Supplementary Fig. 9). $T_{opt}^{eco}$ is indeed slightly (1.4 °C) lower than $T_{max gs}^{air}$ 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 $T_{opt}^{eco}$ and $T_{max gs}^{air}$ (Fig. 3a and Supplementary Fig. 9). Analysis of the 16-day averaged $T_{max}^{air}$ distribution during the period when $T_{opt}^{eco}$ is observed further shows that the rank of $T_{opt}^{eco}$ in the $T_{max}^{air}$ distribution is already
near the highest quantile of $T_{max}^{air}$ (>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.

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⁴⁶. We compared these $T_{max gs}^{air}$ projections with the present-day distribution of $T_{opt}^{eco}$ with a focus on tropical evergreen forests, where optimum temperature is currently just below the limit of $T_{max gs}^{air}$ (see Methods; Fig. 3b and 3c). The key uncertainty in this discussion is, however, whether or not $T_{opt}^{eco}$ will acclimate and follow the increase in $T_{max gs}^{air}$. We therefore looked at possible acclimation from time series of $T_{opt}^{eco}$ 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 $T_{opt}^{eco}$ 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⁴⁷ (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⁴⁸. Also the annual $T_{opt}^{eco}$ derived from flux sites estimates of GPP did not exhibit a positive trend and was not significantly correlated with annual variations of $T_{max gs}^{air}$, although the flux time series are probably too short to properly evaluate trends of $T_{opt}^{eco}$ 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 $T_{max gs}^{air}$ projections from climate models with the current distribution of $T_{opt}^{eco}$. Under this assumption, the average $T_{max gs}^{air}$ of tropical evergreen forests will exceed the current value of $T_{opt}^{eco}$ 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 $T_{opt}^{eco}$ staying above $T_{max gs}^{air}$ 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⁴⁹ as applied in several studies using observed spatial gradients to hindcast temporal changes⁵⁰,⁵¹. Here, we assume that temporal change of $T_{opt}^{eco}$ will evolve proportionally to $T_{max gs}^{air}$, following the spatial temperature sensitivity of $T_{opt}^{eco}$ to $T_{max gs}^{air}$ 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, $T_{max gs}^{air}$ will still surpass $T_{opt}^{eco}$ 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 $T_{opt}^{eco}$ 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 $T_{opt}^{eco}$ 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 $T_{opt}^{eco}$. Furthermore, the acclimation of plants to increasing atmospheric CO$_2$ 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 $T_{opt}^{eco}$ over time is a priority for future studies. Continuous monitoring and dedicated manipulative experiments could improve our understanding of the features of $T_{opt}^{eco}$ 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 ($T_{max}^{air}$) 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 (NIR$_V$), which can differentiate between the confounding effects of background brightness, leaf area and the distribution of photosynthetic capacity with depth in canopies$^{19}$. NIR$_V$ is calculated as the product of total scene NIR reflectance (NIR$_T$) and Normalized Difference Vegetation Index (NDVI)$^{19}$. As a proxy of photosynthesis, NIR$_V$ 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 NIR$_V$ to calculate and map the optimal air temperature for vegetation productivity at an ecosystem scale ($T_{opt}^{eco}$). Following Badgley et al.$^{19}$, we calculated 16-day NIR$_V$ 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 NIR$_V$ 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⁵⁴. 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²⁰,⁵⁵,⁵⁶. 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²¹, and is found to correlate well with estimated GPP on a site-by-site basis⁵⁷. 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²¹.

## 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⁵⁸. Following the previous studies⁵⁸,⁵⁹, 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)²⁰. 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⁶⁰. 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⁶⁰, 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.⁶¹, 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⁶².

## 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⁶³. Generally, the occurrence time of $T_{max}^{surface}$ (2:00–4:00 P.M.) is relatively close to the Aqua overpass time (1:30 P.M.); thus, we assumed that $T_{max}^{surface}$ from MODIS-Aqua is comparable with the daily maximum leaf surface temperature ($T_{max}^{leaf}$).

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 ($T_{max}^{air}$) of the growing season averaged over 2001–2013 was calculated as the current mean growing-season daily maximum air temperature ($T_{maxgs}^{air}$). Information on the growing season was derived from the study by Zhu et al.⁶⁴, which was first determined from the GIMMS Leaf Area Index data set (GIMMS LAI₃g) 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)⁶⁵.

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⁴⁶ to determine the impact of future warming on vegetation productivity (see model list in Supplementary Table 2).
Analysis

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. $T_{max gs}^{air}$ for the late 21st century was estimated using the same temperature thresholds as for the current $T_{max gs}^{air}$. All GCM projections were resampled to a resolution of 1° using first-order conservative interpolation method<sup>66</sup>.

We estimated local $T_{opt}^{eco}$ by examining the temperature response curve of MODIS NIR<sub>V</sub>. Following Yuan et al.<sup>37</sup> and Niu et al.<sup>18</sup>, NIR<sub>V</sub> 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 NIR<sub>V</sub> data as the response of NIR<sub>V</sub> 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 NIR<sub>V</sub>. The $T_{opt}^{eco}$ was then determined from the response curve at which NIR<sub>V</sub> was maximized (Supplementary Fig. 16). Note that $T_{opt}^{eco}$ may not be detected for some pixels where the maximum NIR<sub>V</sub> was only attained at either end of the response curve, accounting for 3.5% of the vegetated areas. Only vegetated areas with detectable $T_{opt}^{eco}$ were shown when mapping the spatial pattern of $T_{opt}^{eco}$. The derivation of $T_{opt}^{eco}$ 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 NIR<sub>V</sub> to calculate $T_{opt}^{eco}$, we also applied nonlinear regression of the photosynthetic temperature response data (Eq.1) to estimate $T_{opt}^{eco}$, which produced similar results (Supplementary Fig. 2):

$$ \text{NIR}_{V(T)} = \text{NIR}_{V(OPT)} - b(T - T_{opt}^{eco})^2 \quad (\text{Eq.1}) $$

where NIR<sub>V(T)</sub> is the NIR<sub>V</sub> value at a daily maximum temperature T, b is a parameter describing the spread of the parabola<sup>48,67</sup>. $T_{opt}^{eco}$ is the vertex of each fit and NIR<sub>V(OPT)</sub> is the NIR<sub>V</sub> value at $T_{opt}^{eco}$. Finally, we used daily mean air temperature ($T_{mean}^{air}$) instead of ($T_{max}^{air}$) to calculate $T_{opt}^{eco}$. In this test, $T_{opt}^{eco}$ derived from $T_{mean}^{air}$ is smaller than $T_{opt}^{eco}$ estimated from $T_{max}^{air}$, but the two variables were strongly spatially correlated (Supplementary Fig. 6).

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

Air temperature, atmospheric vapor pressure deficit (VPD) and solar radiation usually vary in time and space, so that the empirical observation of spatial patterns of $T_{opt}^{eco}$ in this study cannot be unambiguously attributed to air temperature as a single explaining factor of $T_{opt}^{eco}$. 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²⁴. We show that across climatic gradients $T_{opt}^{eco}$ 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 $T_{max gs}^{air}$ under each VPD bin in the regression model of:

$$T_{opt}^{eco} = k_0 + k_1 \times T_{max gs}^{air} + k_2 \times VPD \quad (\text{Eq.2})$$

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 $T_{opt}^{eco}$ and $T_{max gs}^{air}$, we further calculated the partial (‘intrinsic’) sensitivity of $T_{opt}^{eco}$ to $T_{max gs}^{air}$ in each grid point based on the following bilinear regression:

$$T_{opt}^{eco} = k_0 + k_1 \times T_{max gs}^{air} + k_2 \times VPD + k_3 \times VPD \times T_{max gs}^{air} \quad (\text{Eq.3})$$

where the partial sensitivity of $T_{opt}^{eco}$ to $T_{max gs}^{air}$ is defined as $k_1$ in Eq. 3 under each VPD bin. Then we compared the partial sensitivity with the apparent sensitivity of $T_{opt}^{eco}$ to $T_{max gs}^{air}$ estimated using abovementioned linear regression between $T_{opt}^{eco}$ and $T_{max gs}^{air}$ for each VPD bin. As shown in Supplementary Fig. 19, although the apparent sensitivity of $T_{opt}^{eco}$ to $T_{max gs}^{air}$ is generally lower than the partial (‘intrinsic’) sensitivity of $T_{opt}^{eco}$ to $T_{max gs}^{air}$ the apparent sensitivity to $T_{max gs}^{air}$ 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 $T_{opt}^{eco}$ 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 $T_{opt}^{eco}$ is observed for the 16-day averaged Rad distribution. As shown in Supplementary Fig. 20, the Rad value when $T_{opt}^{eco}$ 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 $T_{opt}^{eco}$ is consistent the time of maximum solar radiation. This is because $T_{opt}^{eco}$ in these regions generally appears in summer, which is also the period when solar radiation displays its maximum during the year.

The NIR<sub>V</sub>-derived $T_{opt}^{eco}$ was compared with $T_{opt}^{eco}$ estimated using GPP data from 153 eddy covariance sites. Flux-derived $T_{opt}^{eco}$ 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 $T_{opt}^{eco}$ using NIR<sub>V</sub> datasets. A robust estimate of $T_{opt}^{eco}$ can be derived for 125 sites (Supplementary Table 1). For each site, we calculated the mean and SD of $T_{opt}^{eco}$ across different years. Then, we extracted and averaged $T_{opt}^{eco}$ values within a 3x3 pixel window around each site from the NIR<sub>V</sub>-derived $T_{opt}^{eco}$ map, and calculated the SD of the nine $T_{opt}^{eco}$ values within the window. The relationship between NIR<sub>V</sub>- and flux-derived $T_{opt}^{eco}$ 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 NIR<sub>V</sub>-derived $T_{opt}^{eco}$ is comparable with that estimated independently from measurements of flux-tower eddy covariance (Fig. 1b).

We compared the spatial distribution of $T_{opt}^{eco}$ derived from NIR<sub>V</sub> with the one obtained from NDVI datasets. Consistent spatial patterns of $T_{opt}^{eco}$ are derived from each of the three NDVI datasets (Supplementary Fig. 21). A global composite map of $T_{opt}^{eco}$ (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 $T_{opt}^{eco}$ to a common grid of 8 km before averaging. $T_{opt}^{eco}$ from NDVI datasets generally show a spatial pattern similar to that from NIR<sub>V</sub>, but with smaller NDVI-derived $T_{opt}^{eco}$ values for central Australia and southern South America (Supplementary Fig. 3). We compared the spatial distribution of $T_{opt}^{eco}$ derived from NIR<sub>V</sub> with that from MODIS EVI data during 2001–2013, and found that the EVI derived $T_{opt}^{eco}$ shows very similar spatial pattern to that of NIR<sub>V</sub> derived $T_{opt}^{eco}$ (Supplementary Fig. 4). The distribution of $T_{opt}^{eco}$ derived from NIR<sub>V</sub> and from GOSAT SIF datasets also have similar spatial patterns, even though the NIR<sub>V</sub>-derived $T_{opt}^{eco}$ 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⁶⁸. Much of the effort to date to understand variability in the leaf-level photosynthesis-temperature response has focused on biochemical processes⁶⁸, with $V_{cmax}$ and $J_{max}$ being two major parameters controlling the maximum rates of photosynthesis limited by CO₂ and light, respectively⁶⁹. Therefore, in this study, we compared $T_{opt}^{eco}$ derived using GPP proxies with leaf-scale optimal temperature of maximum Rubisco-limited carboxylation rates ($V_{cmax}$), 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⁷⁰, the optimum temperature of GPP ($T_{opt}^{eco}$) should, in theory, be lower than the optimal temperature of maximum Rubisco-limited carboxylation rates ($V_{cmax}$). For this comparison to be made, we extracted and averaged $T_{opt}^{eco}$ values within a 3×3 pixel window from the NIRV-derived $T_{opt}^{eco}$ 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 $T_{opt}^{eco}$ values across pixels with both the same growing season mean temperature and the same plant functional type as the corresponding site.

$T_{opt}^{eco}$ is different from $T_{opt}^{leaf}$ not only because of respiratory process, but also because air temperature can differ from leaf temperatures⁷¹, which are regulated by leaf traits affecting the leaf energy balance⁷². Because, to our knowledge, global gridded monthly leaf temperature data are not available, we use daily maximum land surface temperature ($T_{max}^{surface}$) from MODIS to calculate $T_{opt\_LST}^{eco}$ in order to illustrate the potential differences between $T_{opt\_LST}^{eco}$ and $T_{opt}^{eco}$. As shown in Supplementary Figure 22, the $T_{opt\_LST}^{eco}$ is similar to $T_{opt}^{eco}$ over tropical savannas. However, over moist tropical forests $T_{opt\_LST}^{eco}$ is lower than $T_{opt}^{eco}$ which can be explained by the lower daytime surface temperature than air temperature as a result of strong evapotranspiration effects⁷¹,⁷³. This ecosystem-dependent difference between $T_{opt\_LST}^{eco}$ and $T_{opt}^{eco}$ suggests that the leaf thermal regulation mechanism through the physiological and morphological changes⁷² is an important ecosystem process shaping spatial variations of $T_{opt}^{eco}$. 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 $T_{opt}^{eco}$ under future warming, we estimated the acclimated $T_{opt}^{eco}$ for vegetation productivity by the end of the 21st century (2091–2100), using recent
IPCC climate projections⁴⁶. To this end, we applied the space-for-time substitution approach⁴⁹, assuming that that temporally $T_{opt}^{eco}$ will evolve proportionally to $T_{max gs}^{air}$ following the spatial temperature sensitivity of $T_{opt}^{eco}$ to $T_{max gs}^{air}$. Given the relatively large uncertainties of precipitation projections, we considered two future precipitation scenarios. For the first scenario, we estimated acclimated $T_{opt}^{eco}$ pixel by pixel using the temperature sensitivity of $T_{opt}^{eco}$ 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 $T_{opt}^{eco}$ was calculated pixel by pixel using the temperature sensitivity of $T_{opt}^{eco}$ under the projected MAP level for 2091–2100. Acclimated $T_{opt}^{eco}$ was averaged across the GCMs under each scenario. Latitudinal variation of future $T_{opt}^{eco}$ was derived by averaging within 1°-latitude bins from future $T_{opt}^{eco}$ maps and then compared with that in future $T_{max gs}^{air}$ summarized by latitude from future $T_{max gs}^{air}$ maps.

## Supplementary Material

Refer to Web version on PubMed Central for supplementary material.

## Acknowledgments:

This study was supported by the Strategic Priority Research Program (A) of the Chinese Academy of Sciences (grant XDA20050101), the National Natural Science Foundation of China (41530528) and the National Key R&D Program of China (2017YFA0604702). This work used eddy covariance data acquired by the FLUXNET community and in particular by the following networks: AmeriFlux (U.S. Department of Energy, Biological and Environmental Research, Terrestrial Carbon Program (DE-FG02-04ER63917 and DE-FG02-04ER63911)), AfriFlux, AsiaFlux, CarboAfrica, CarboEuropeIP, CarboItaly, CarboMont, ChinaFlux, Fluxnet-Canada (supported by CFCAS, NSERC, BIOCAP, Environment Canada, and NRCan), GreenGrass, KoFlux, LBA, NECC, OzFlux, TCOS-Siberia, USCCC. We acknowledge the financial support to the eddy covariance data harmonization provided by CarboEuropeIP, FAO-GTOS-CO, iLEAPS, Max Planck Institute for Biogeochemistry, National Science Foundation, University of Tuscia, Université Laval and Environment Canada and US Department of Energy and the database development and technical support from Berkeley Water Center, Lawrence Berkeley National Laboratory, Microsoft Research eScience, Oak Ridge National Laboratory, University of California-Berkeley, University of Virginia. P.C., J.P., and I.A.J. would like to acknowledge the financial support from the European Research Council Synergy grant ERC-SyG-2013-610028 IMBALANCE-P. P.C. was also supported from the French Agence Nationale de la Recherche (ANR) Con-vergence Lab Changement climatique et usage des terres (CLAND). I.A.J. acknowledges the Methusalem funding of the Flemish Community through the Research Council of the University of Antwerp. T.F.K. was supported by the NASA Terrestrial Ecology Program IDS Award NNH17AE86I. J. Mao and X. Shi are supported by the Terrestrial Ecosystem Science Scientific Focus Area project funded through the Terrestrial Ecosystem Science Program in the Climate and Environmental Sciences Division of the Biological and Environmental Research Program in the US Department of Energy Office of Science. Oak Ridge National Laboratory is supported by the Office of Science of the US Department of Energy under Contract No. DE-AC05-00OR22725. M. C. was supported by a grant overseen by the French National Research Agency (ANR) as part of the "Investissements d'Avenir" program (ANR-11-LABX-0002-01, Lab of Excellence ARBRE).

## References

1. Berry J & Bjorkman O Photosynthetic response and adaptation to temperature in higher plants. Annual Review of Plant Physiology 31, 491–543 (1980).

2. Hughes L Biological consequences of global warming: is the signal already apparent? Trends in Ecology & Evolution 15, 56–61 (2000). [PubMed: 10652556]

3. Niu S et al. Climatic warming changes plant photosynthesis and its temperature dependence in a temperate steppe of northern China. Environmental and Experimental Botany 63, 91–101 (2008).
4. Way DA, Yamori W Thermal acclimation of photosynthesis: on the importance of adjusting our definitions and accounting for thermal acclimation of respiration. Photosynthesis research 119, 89–100 (2014). [PubMed: 23812760]

5. Kattge J, Knorr W Temperature acclimation in a biochemical model of photosynthesis: a reanalysis of data from 36 species. Plant, cell & environment 30, 1176–1190 (2007).

6. Lloyd J, Farquhar GD Effects of rising temperatures and [CO₂] on the physiology of tropical forest trees. Philosophical Transactions of the Royal Society B: Biological Sciences 363, 1811–1817 (2008).

7. Medlyn B et al. Temperature response of parameters of a biochemically based model of photosynthesis. II. A review of experimental data. Plant, Cell & Environment 25, 1167–1179 (2002).

8. Field CB, Randerson JT, Malmström CM Global net primary production: combining ecology and remote sensing. Remote sensing of Environment 51, 74–88 (1995).

9. Niinemets Ü Global-scale climatic controls of leaf dry mass per area, density, and thickness in trees and shrubs. Ecology 82, 453–469 (2001).

10. Williams AP et al. Temperature as a potent driver of regional forest drought stress and tree mortality. Nature Climate Change 3, 292–297 (2013).

11. Tyree MT & Dixon MA Water stress induced cavitation and embolism in some woody plants. Physiologia Plantarum 66, 397–405 (1986).

12. Yin X, Struik P C₃ and C₄ photosynthesis models: an overview from the perspective of crop modelling. NJAS-Wageningen Journal of Life Sciences 57, 27–38 (2009).

13. Smith NG, Malyshev SL, Shevliakova E, Kattge J, Dukes JS, Foliar temperature acclimation reduces simulated carbon sensitivity to climate. Nature Climate Change 6, 407 (2016).

14. Mercado LM et al. Large sensitivity in land carbon storage due to geographical and temporal variation in the thermal response of photosynthetic capacity. New Phytologist 218, 1462–1477 (2018). [PubMed: 29635689]

15. Corlett RT Impacts of warming on tropical lowland rainforests. Trends in Ecology & Evolution 26, 606–613 (2011). [PubMed: 21803440]

16. Doughty CE, Goulden ML Are tropical forests near a high temperature threshold? Journal of Geophysical Research: Biogeosciences 113, G00B07 (2008).

17. Baldocchi D, Falge E, Gu L et al. FLUXNET: A new tool to study the temporal and spatial variability of ecosystem-scale carbon dioxide, water vapor, and energy flux densities. Bulletin of the American Meteorological Society 82, 2415–2434 (2001).

18. Niu et al. Thermal optimality of net ecosystem exchange of carbon dioxide and underlying mechanisms. New Phytologist 194, 775–783 (2012). [PubMed: 22404566]

19. Badgley G, Field CB, Berry JA Canopy near-infrared reflectance and terrestrial photosynthesis. Science Advances 3, e1602244 (2017). [PubMed: 28345046]

20. Pinzon JE, Tucker CJ A non-stationary 1981–2012 AVHRR NDVI₃g time series. Remote Sensing 6, 6929–6960 (2014).

21. Huete A et al. Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sensing of Environment 83, 195–213 (2002).

22. Joiner J, Yoshida Y, Vasilkov A, Middleton E First observations of global and seasonal terrestrial chlorophyll fluorescence from space. Biogeosciences 8, 637–651 (2011).

23. Ma S, Osuna JL, Verfaillie J et al. Photosynthetic responses to temperature across leaf-canopy-ecosystem scales: a 15-year study in a Californian oak-grass savanna. Photosynthesis research 132, 277–291 (2017). [PubMed: 28425026]

24. Stocker BD et al. Quantifying soil moisture impacts on light use efficiency across biomes. New Phytologist 218, 1430–1449 (2018). [PubMed: 29604221]

25. Bunce JA Acclimation of photosynthesis to temperature in eight cool and warm climate herbaceous C3 species: temperature dependence of parameters of a biochemical photosynthesis model. Photosynthesis Research 63, 59–67 (2000). [PubMed: 16252165]
26. Han Q, Kawasaki T, Nakano T & Chiba Y Spatial and seasonal variability of temperature responses of biochemical photosynthesis parameters and leaf nitrogen content within a Pinus densiflora crown. Tree physiology 24, 737-744 (2004). [PubMed: 15123445]

27. Harley P, Thomas R, Reynolds J & Strain B Modelling photosynthesis of cotton grown in elevated CO₂. Plant, Cell & Environment 15, 271-282 (1992).

28. Onoda Y, Hikosaka K & Hirose T The balance between RuBP carboxylation and RuBP regeneration: a mechanism underlying the interspecific variation in acclimation of photosynthesis to seasonal change in temperature. Functional Plant Biology 32, 903-910 (2005).

29. Walcroft A, Le Roux X, Diaz-Espejo A, Dones N & Sinoquet H Effects of crown development on leaf irradiance, leaf morphology and photosynthetic capacity in a peach tree. Tree Physiology 22, 929-938 (2002). [PubMed: 12204849]

30. Wang KY, Kellomäki S & Laitinen K Acclimation of photosynthetic parameters in Scots pine after three years exposure to elevated temperature and CO₂. Agricultural and Forest Meteorology 82, 195-217 (1996).

31. Vaz M et al. Drought-induced photosynthetic inhibition and autumn recovery in two Mediterranean oak species (Quercus ilex and Quercus suber). Tree Physiology 30, 946-956 (2010). [PubMed: 20571151]

32. Zhou L et al. Responses of photosynthetic parameters to drought in subtropical forest ecosystem of China. Scientific reports 5, 18254 (2015). [PubMed: 26666469]

33. Doughty CE, Goulden ML Seasonal patterns of tropical forest leaf area index and CO₂ exchange. Journal of Geophysical Research: Biogeosciences 113, G00B06, doi: 10.1029/2007JG000590 (2008).

34. Wu J et al. Partitioning controls on Amazon forest photosynthesis between environmental and biotic factors at hourly to interannual timescales. Global Change Biology, 23, 1240-1257 (2017). [PubMed: 27644012]

35. Gunderson CA, O'Hara KH, Campion CM, Walker AV, Edwards NT Thermal plasticity of photosynthesis: the role of acclimation in forest responses to a warming climate. Global Change Biology 16, 2272-2286 (2010).

36. Mooney HA, Björkman O, Collatz GJ Photosynthetic acclimation to temperature in the desert shrub, Larrea divaricata. I. Carbon dioxide exchange characteristics of intact leaves. Plant Physiology 61, 406-410 (1978). [PubMed: 16660303]

37. Yuan W et al. Thermal adaptation of net ecosystem exchange. Biogeosciences 8, 1453-1463 (2011).

38. Sage RF, Kubien DS The temperature response of C₃ and C₄ photosynthesis. Plant, Cell & Environment 30, 1086-1106 (2007).

39. Armond PA, Schreiber U & Björkman O Photosynthetic acclimation to temperature in the desert shrub, Larrea divaricata: II. Light-harvesting efficiency and electron transport. Plant Physiology 61, 411-415 (1978). [PubMed: 16660304]

40. Badger M, Björkman O & Armond P An analysis of photosynthetic response and adaptation to temperature in higher plants: temperature acclimation in the desert evergreen Nerium oleander L. Plant, Cell & Environment 5, 85-99 (1982).

41. Atkin O, Scheurwater I & Pons T High thermal acclimation potential of both photosynthesis and respiration in two lowland Plantago species in contrast to an alpine congeneric. Global Change Biology 12, 500-515 (2006).

42. Doughty CE An in situ leaf and branch warming experiment in the Amazon. Biotropica 43, 658-665 (2011).

43. Koch GW, Amthor JS, Goulden ML Diurnal patterns of leaf photosynthesis, conductance and water potential at the top of a lowland rain forest canopy in Cameroon: measurements from the Radeau des Cimes. Tree Physiology 14, 347-360 (1994). [PubMed: 14967691]

44. Tribuzy ES Variações da temperatura foliar do dossel eo seu efeito na taxa assimilatória de CO₂ na Amazônia Central, Universidade de São Paulo (2005).

45. White A, Cannell MG, Friend AD CO₂ stabilization, climate change and the terrestrial carbon sink. Global Change Biology 6, 817-833 (2000).
46. IPCC, Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (Cambridge University Press, Cambridge, UK, 2013).

47. Xu L et al. Temperature and vegetation seasonality diminishment over northern lands. Nature Climate Change 3, 581 (2013).

48. Sendall KM et al. Acclimation of photosynthetic temperature optima of temperate and boreal tree species in response to experimental forest warming. Global Change Biology 21, 1342–1357 (2015). [PubMed: 25354151]

49. Likens GE Long-term Studies in Ecology (Springer, 1989).

50. Jung M et al. Global patterns of land-atmosphere fluxes of carbon dioxide, latent heat, and sensible heat derived from eddy covariance, satellite, and meteorological observations. Journal of Geophysical Research: Biogeosciences 116, G00J07 (2011).

51. Lester RE, Close PG, Barton JL, Pope AJ & Brown SC Predicting the likely response of data-poor ecosystems to climate change using space-for-time substitution across domains. Global Change Biology 20, 3471–3481 (2014). [PubMed: 24832685]

52. Cavaleri MA, Reed SC, Smith WK, Wood TE Urgent need for warming experiments in tropical forests. Global Change Biology 21, 2111–2121 (2015). [PubMed: 25641092]

53. Papale D et al. Towards a standardized processing of Net Ecosystem Exchange measured with eddy covariance technique: algorithms and uncertainty estimation. Biogeosciences 3, 571–583 (2006).

54. Myneni RB, Ramakrishna R, Nemani R, Running SW Estimation of global leaf area index and absorbed PAR using radiative transfer models. IEEE Transactions on Geoscience and remote sensing 35, 1380–1393 (1997).

55. Maisongrande P, Duchemin B, Dedieu G VEGETATION/SPOT: an operational mission for the Earth monitoring; presentation of new standard products. International Journal of Remote Sensing 25, 9–14 (2004).

56. Vermote EF, El Saleous NZ, Justice CO Atmospheric correction of MODIS data in the visible to middle infrared: first results. Remote Sensing of Environment 83, 97–111 (2002).

57. Rahman AF, Sims DA, Cordova VD, El-Masri BZ Potential of MODIS EVI and surface temperature for directly estimating per-pixel ecosystem C fluxes. Geophysical Research Letters 32, 156–171 (2005).

58. Frankenberg C et al. New global observations of the terrestrial carbon cycle from GOSAT: Patterns of plant fluorescence with gross primary productivity. Geophysical Research Letters 38, 351–365 (2011).

59. Lee JE et al. Forest productivity and water stress in Amazonia: observations from GOSAT chlorophyll fluorescence. Proceedings of the Royal Society B Biological Sciences 280, 176–188 (2013).

60. Kottek M, Grieser J, Beck C, Rudolf B, Rubel F World map of the Köppen-Geiger climate classification updated. Meteorologische Zeitschrift 15, 259–263 (2006).

61. Qin J, Yang K, Liang S & Guo X The altitudinal dependence of recent rapid warming over the Tibetan Plateau. Climatic Change 97, 321 (2009).

62. Liu D et al. Contrasting responses of grassland water and carbon exchanges to climate change between Tibetan Plateau and Inner Mongolia. Agricultural and Forest Meteorology 249, 163–175 (2018).

63. Wan Z & Dozier J A generalized split-window algorithm for retrieving land-surface temperature from space. IEEE Transactions on geoscience and remote sensing 34, 892–905 (1996).

64. Zhu Z et al. Greening of the Earth and its drivers. Nature Climate Change, 6, 701–795 (2016).

65. Weedon GP et al. The WFDEI meteorological forcing data set: WATCH Forcing Data methodology applied to ERA-Interim reanalysis data. Water Resources Research 50, 7505–7514 (2014).

66. Jones PW, First-and second-order conservative remapping schemes for grids in spherical coordinates. Monthly Weather Review 127, 2204–2210 (1999).

67. Battaglia M, Beadle C, Loughhead S Photosynthetic temperature responses of Eucalyptus globulus and Eucalyptus nitens. Tree Physiology 16, 81–89 (1996). [PubMed: 14871750]
68. Lin YS, Medlyn BE & Ellsworth DS Temperature responses of leaf net photosynthesis: the role of component processes. Tree Physiology 32, 219–231 (2012). [PubMed: 22278379]

69. Farquhar G. v., von Caemmerer S. v. & Berry J A biochemical model of photosynthetic CO₂ assimilation in leaves of C₃ species. Planta 149, 78–90 (1980). [PubMed: 24306196]

70. Baker NR Photosynthesis and the Environment (Kluwer Academic Publishers, 2006).

71. Helliker BR & Richter SL Subtropical to boreal convergence of tree-leaf temperatures. Nature 454, 511 (2008). [PubMed: 18548005]

72. Michaletz ST et al. The energetic and carbon economic origins of leaf thermoregulation. Nature Plants 2, 16129 (2016). [PubMed: 27548589]

73. Lian X et al. Spatiotemporal variations in the difference between satellite-observed daily maximum land surface temperature and station-based daily maximum near-surface air temperature. Journal of Geophysical Research: Atmospheres 122, 2254–2268 (2017).
<figure><img src="image_5.png" /><figcaption>Fig. 1 | Distribution of ecosystem-scale optimal temperature ($T_{opt}^{eco}$) for vegetation productivity derived from flux tower sites and satellite-based data for near-infrared reflectance of vegetation (NIRV).</figcaption></figure>

a, Relationship between mean annual daily maximum air temperature during the growing season ($T_{max gs}^{air}$) and $T_{opt}^{eco}$ derived from daily measurements of photosynthesis across eddy-covariance sites. Flux-derived $T_{max gs}^{air}$ and $T_{opt}^{eco}$ 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 $T_{opt}^{eco}$ derived from flux data and $T_{opt}^{eco}$ derived from NIRV data. For each site, we extracted and averaged $T_{opt}^{eco}$ values within a 3×3 pixel window around the site from NIRV-derived $T_{opt}^{eco}$ map, and calculated the SD of the nine $T_{opt}^{eco}$ 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 $T_{opt}^{eco}$ for vegetation productivity (left panel), and $T_{opt}^{eco}$ averaged by latitude (right panel). $T_{opt}^{eco}$ 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 $T_{opt}^{eco}$ 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 $T_{opt}^{eco}$ 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 $T_{opt}^{eco}$ summarized by latitude. **d**, $T_{opt}^{eco}$ in the climate space (left panel) and the temperature sensitivity of $T_{opt}^{eco}$ along the precipitation gradient (right panel). Each climate bin is defined by 1-°C intervals of $T_{max gs}^{air}$ 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 $T_{opt}^{eco}$ along the precipitation gradient, calculated as the slope of the linear regression between $T_{opt}^{eco}$ and $T_{max gs}^{air}$ for a given precipitation level. The shaded area indicates the SD of temperature sensitivity of $T_{opt}^{eco}$ estimated by bootstrapping.
<figure><img src="image_6.png" /><figcaption>**Fig. 3 |.** Change with latitude in ecosystem-scale optimal temperature for vegetation productivity ($T_{opt}^{eco}$) and daily maximum air temperature averaged over the growing season ($T_{max gs}^{air}$).</figcaption></figure>

a, Current $T_{opt}^{eco}$ versus current $T_{max gs}^{air}$; b, Future $T_{opt}^{eco}$ versus future $T_{max gs}^{air}$. Current $T_{opt}^{eco}$ and $T_{max gs}^{air}$ are calculated using current temperature for 2001–2013, whereas acclimated $T_{opt}^{eco}$ and future $T_{max gs}^{air}$ 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 $T_{opt}^{eco}$ is determined based on the projected temperature and temperature sensitivity of $T_{opt}^{eco}$ 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 $T_{opt}^{eco}$ or $T_{max gs}^{air}$ summarized by latitude. c, Future $T_{opt}^{eco}$ versus future $T_{max gs}^{air}$ for tropical evergreen forests. ** indicates that $T_{opt}^{eco}$ is significantly lower than $T_{max gs}^{air}$ at P<0.01 in a paired t-test. Error bars indicate ±SD.