This is the accepted version of the article: Huang, M.; Piao, S.; Wang, X.; [et al.]. «Air temperature optima of vegetation productivity across global biomes». Nature ecology and evolution, Vol. 3, issue 5 (May 2019), p. 772–779. DOI 10.1038/s41559-019-0838-x This version is avaible at https://ddd.uab.cat/record/216962 under the terms of the license 1 Air temperature optima of vegetation productivity across global biomes 2 Mengtian Huang1, Shilong Piao1,2,3*, Philippe Ciais4, Josep Peñuelas5,6, Xuhui Wang1, Trevor F. 3 Keenan7,8, Shushi Peng1, Joseph A. Berry9, Kai Wang1, Jiafu Mao10, Ramdane Alkama11, 4 Alessandro Cescatti11, Matthias Cuntz12, Hannes De Deurwaerder13, Mengdi Gao1, Yue He1, 5 Yongwen Liu1, Yiqi Luo14, Ranga B. Myneni15, Shuli Niu16, Xiaoying Shi10, Wenping Yuan17, 6 Hans Verbeeck13, Tao Wang2,3, Jin Wu18, Ivan A. Janssens19 7 1 8 Sino-French Institute for Earth System Science, College of Urban and Environmental Sciences, 9 Peking University, Beijing 100871, China. 2 10 Key Laboratory of Alpine Ecology and Biodiversity, Institute of Tibetan Plateau Research, 11 Chinese Academy of Sciences, Beijing 100085, China. 3 12 Center for Excellence in Tibetan Earth Science, Chinese Academy of Sciences, Beijing 100085, 13 China. 14 4 Laboratoire des Sciences du Climat et de l’Environnement, CEA CNRS UVSQ, Gif-sur-Yvette 15 91191, France. 5 16 CREAF, Cerdanyola del Vallès, Barcelona 08193, Catalonia, Spain. 6 17 CSIC, Global Ecology Unit CREAF-CSIC-UAB, Bellaterra, Barcelona 08193, Catalonia, Spain. 7 18 Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA. 8 19 Department of Environmental Science, Policy and Management, UC Berkeley, Berkeley, 20 California 94720, USA. 9 21 Department of Global Ecology, Carnegie Institution for Science, Stanford, California 94305, 22 USA. 10 23 Climate Change Science Institute and Environmental Sciences Division, Oak Ridge National 24 Laboratory, Oak Ridge, TN, USA. 11 25 European Commission, Joint Research Centre, Institute for Environment and Sustainability, 26 Ispra, Italy. 12 27 INRA, Université de Lorraine, AgroParisTech, UMR Silva, 54000 Nancy, France. 13 28 CAVElab Computational & Applied Vegetation Ecology, Facultyof Bioscience Engineering, 29 Ghent University, Coupure Links 653, B-9000 Gent, Belgium. 14 30 Department of Biological Sciences, Northern Arizona University, Arizona, 86011, USA. 15 31 Department of Earth and Environment, Boston University, Boston, Massachusetts 02215, USA. 16 32 Key Laboratory of Ecosystem Network Observation and Modeling, Institute of Geographic 33 Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing, China. 17 34 School of Atmospheric Sciences, Center for Monsoon and Environment Research, Sun Yatsen 35 University, Guangzhou 510275, China. 18 36 Environmental & Climate Sciences Department, Brookhaven National Laboratory, Upton, NY, 37 11973, USA. 19 38 Centre of Excellence PLECO (Plant and Vegetation Ecology), Department of Biology, 39 University of Antwerp, Universiteitsplein 1, B-2610 Wilrijk, Belgium. 40 41 * Correspondence to: slpiao@pku.edu.cn 42 43 44 45 46 The global distribution of the optimum air temperature for ecosystem-level gross primary 47 productivity (𝑻𝒆𝒄𝒐 𝒐𝒑𝒕 ) is poorly understood, despite its importance for ecosystem carbon uptake 48 under future warming. We provide empirical evidence for the existence of such an optimum, 49 using measurements of in situ eddy covariance and satellite-derived proxies, and report its 50 global distribution. 𝑻𝒆𝒄𝒐 𝒐𝒑𝒕 is consistently lower than the physiological optimum temperature of 51 leaf-level photosynthetic capacity, which typically exceeds 30 °C. The global average 𝑻𝒆𝒄𝒐 𝒐𝒑𝒕 is 52 estimated to be 23±6 ºC, with warmer regions having higher 𝑻𝒆𝒄𝒐 𝒐𝒑𝒕 values than colder regions. 53 In tropical forests, particularly, 𝑻𝒆𝒄𝒐 𝒐𝒑𝒕 is close to growing-season air temperature and is 54 projected to fall below it under all scenarios of future climate, suggesting a limited safe 55 operating space for these ecosystems under future warming. 56 57 Understanding how photosynthesis responds to warming has been a focus in plant research in recent 58 decades, and most of the existing knowledge comes from leaf-scale measurements1-4. Most leaf- 59 scale temperature response curves show that photosynthetic capacity increases with temperature up 𝑙𝑒𝑎𝑓 60 to an optimum temperature (𝑇𝑜𝑝𝑡 ), which typically occurs in the 30-40°C temperature range5,6. 61 Above this optimum temperature, foliar photosynthetic capacity sharply declines as electron- 62 transport and Rubisco enzymatic capacities become impaired7. Field et al.8 first suggested that 𝑒𝑐𝑜 𝑙𝑒𝑎𝑓 63 ecosystem-scale optimum temperature (𝑇𝑜𝑝𝑡 ) may differ from 𝑇𝑜𝑝𝑡 . At ecosystem scale, elevated 64 air temperatures do limit canopy photosynthesis by other processes than leaf carboxylation rates. 65 For instance, elevated air temperatures may accelerate leaf ageing and increase leaf thickness 66 (phenology; e.g. ref 9), and control stomatal closure, because a higher temperature usually comes 67 with a higher vapor pressure deficit (VPD)10. In a more extreme case, warming-induced water stress 68 could suppress canopy photosynthesis though partial hydraulic failure (hydraulics) by cavitation 69 (e.g. ref 11). 70 71 Empirical leaf-scale photosynthesis-temperature relationships12 have been directly incorporated 72 into global ecosystem models, with variants to account for acclimation, i.e. a temporal adjustment 73 of optimum photosynthetic temperature to air temperature during growth5,13,14. This direct scaling 74 of temperature responses from leaves to ecosystems partly determines model projections of Gross 75 Primary Productivity (GPP) and CO2 uptake by terrestrial ecosystems in climatic scenarios. 𝑒𝑐𝑜 76 Verifying the existence of 𝑇𝑜𝑝𝑡 in real world ecosystems, defining its spatial distribution across 𝑒𝑐𝑜 77 and within biomes and understanding the relationships between 𝑇𝑜𝑝𝑡 , prevailing air temperature 𝑙𝑒𝑎𝑓 78 and 𝑇𝑜𝑝𝑡 are important for evaluating models and understanding the impacts of various targets of 79 climatic warming targets on ecosystem productivity. 80 𝑒𝑐𝑜 81 In this study, we formulate and test the following hypotheses: (i) 𝑇𝑜𝑝𝑡 is higher for biomes where 𝑒𝑐𝑜 𝑙𝑒𝑎𝑓 82 air temperature during growth is warmer, (ii) 𝑇𝑜𝑝𝑡 is lower than 𝑇𝑜𝑝𝑡 for any given ecosystem, 83 because the above mentioned limitations of stomatal conductance and phenology emerge before 84 temperature begins to impair foliar photosynthetic capacity, and (iii) tropical forests already operate 𝑒𝑐𝑜 85 near a high 𝑇𝑜𝑝𝑡 above which canopy photosynthesis may decrease with even moderate air 𝑒𝑐𝑜 86 temperature warming15,16. Here we defined 𝑇𝑜𝑝𝑡 as the daytime air temperature at which GPP is 𝑒𝑐𝑜 87 highest over a period of several years, thus 𝑇𝑜𝑝𝑡 can be empirically determined from productivity 88 observations and proxies (see Methods). 89 90 Results and discussion 91 We first applied this approach on time series of daily GPP derived from CO2 flux measurements at 𝑒𝑐𝑜 92 153 globally distributed eddy covariance sites and found that a robust estimate of 𝑇𝑜𝑝𝑡 could be 𝑒𝑐𝑜 93 derived at 125 out of 153 sites (see Methods). 𝑇𝑜𝑝𝑡 values derived from the FLUXNET data range 𝑒𝑐𝑜 94 from 8.2°C to 35.8°C (Fig. 1a, Supplementary Table 1). Tropical sites have higher 𝑇𝑜𝑝𝑡 values than 𝑒𝑐𝑜 95 temperate and boreal sites (Supplementary Fig. 1), implying a dependence of 𝑇𝑜𝑝𝑡 on background 𝑒𝑐𝑜 96 climate. The FLUXNET multi-site analysis further indicates that across sites 𝑇𝑜𝑝𝑡 values are 𝑎𝑖𝑟 97 positively correlated with growing-season mean daily maximum air temperature (𝑇𝑚𝑎𝑥 𝑔𝑠 , see 98 calculation in Methods) (R=0.46, P<0.01, t test), with a spatial linear regression slope of 0.61 °C 99 per °C across sites (Fig. 1a). Overall, these results confirm our first hypothesis stating that higher 𝑒𝑐𝑜 100 𝑇𝑜𝑝𝑡 values occur where higher growth temperatures prevail, in support of findings of Baldocchi 101 et al.17 and Niu et al.18. 102 103 Since eddy covariance measurements do not have a continuous spatial coverage, we also used 104 satellite observations known to be highly correlated with photosynthetic activity19, that is, GPP 105 proxies. The first proxy used is the near-infrared reflectance of terrestrial vegetation (NIRV, the 106 product of total-scene NIR reflectance (NIRT) by the Normalized Difference Vegetation Index 107 (NDVI). NIRV was proven to have a high temporal correlation with GPP at flux tower sites 19. 108 Satellite observations of NIRT and NDVI from the Terra Moderate Resolution Imaging 109 Spectroradiometer (MODIS) were used to calculate NIRV for the period of 2001-2013 (see 𝑒𝑐𝑜 110 Methods). NIRV-derived 𝑇𝑜𝑝𝑡 is comparable with that estimated from eddy covariance flux tower 111 measurements (Fig. 1b), which gives support to using the NIRV proxy for a global mapping of 𝑒𝑐𝑜 𝑒𝑐𝑜 112 𝑇𝑜𝑝𝑡 . The average 𝑇𝑜𝑝𝑡 over the global vegetated areas is estimated to be 23±6 °C (mean ± 1 SD) 113 with large spatial gradients in latitude. As shown in Fig. 1c, maximum values close to 30°C mainly 114 appear over tropical forests, savannas and drylands and minimum values near 10°C prevail at high- 𝑒𝑐𝑜 115 latitudes and in mountainous regions (Fig. 1c). This spatial pattern of 𝑇𝑜𝑝𝑡 is robust to the choice 𝑒𝑐𝑜 116 of a particular climate forcing dataset, or to the method used to estimate 𝑇𝑜𝑝𝑡 (Supplementary Fig. 117 2, see also Methods). Similar results are also found for other GPP proxies (vegetation greenness 118 (NDVI)20, Enhanced Vegetation Index (EVI)21, sun-induced vegetation fluorescence (Sun-induced 𝑎𝑖𝑟 119 Chlorophyll Fluorescence, SIF)22), or when daily mean air temperature (𝑇𝑚𝑒𝑎𝑛 ) is used instead of 𝑎𝑖𝑟 𝑒𝑐𝑜 120 daily maximum air temperature ( 𝑇𝑚𝑎𝑥 ) to calculate 𝑇𝑜𝑝𝑡 (Supplementary Figs. 3-6, see also 121 Methods). Note that although the covariance between air temperature, atmospheric VPD and solar 122 radiation may confuse the direct effect of air temperature on vegetation productivity, we verified 𝑒𝑐𝑜 123 that neither VPD nor radiation is the dominant factor determining the pattern of 𝑇𝑜𝑝𝑡 at the global 124 scale (see Methods). 125 𝑒𝑐𝑜 𝑙𝑒𝑎𝑓 126 In order to test the second hypothesis, we compared satellite-derived 𝑇𝑜𝑝𝑡 with 𝑇𝑜𝑝𝑡 from the 127 responses of maximum Rubisco-limited carboxylation rates (Vcmax) to temperature from leaf-scale 𝑙𝑒𝑎𝑓 128 measurements for 36 species5. Note that the 𝑇𝑜𝑝𝑡 here refers to the temperature optima for leaf- 129 scale (gross) photosynthetic capacity rather than for leaf net photosynthesis, which equals gross 130 photosynthesis minus photorespiration and minus dark respiration (see more details in Methods). 𝑒𝑐𝑜 𝑙𝑒𝑎𝑓 131 We found that 𝑇𝑜𝑝𝑡 is indeed lower than 𝑇𝑜𝑝𝑡 (Supplementary Fig. 7). This difference may 𝑒𝑐𝑜 132 originate from that 𝑇𝑜𝑝𝑡 is additionally limited by high VPD during hot and dry periods6 and by 133 soil-moisture deficits during extensive dry episodes23, under real world conditions. Under 134 conditions of high temperature, atmospheric VPD increases while soil moisture decreases. Stomatal 135 conductance and hence carbon assimilation rates (GPP at ecosystem-scale) decrease to prevent 136 exceedingly low leaf water potentials and resulting plant tissue damage from cavitation24. In 137 contrast, leaf-level photosynthesis measurements that determine the temperature response curve of 138 Vcmax are usually performed in absence of water stress through maintaining relatively low VPD 139 conditions (e.g. ref 25-30) unless the research objective is to investigate drought effect on leaf 140 photosynthetic parameters as the studies by Vaz et al.31 and Zhou et al.32. In addition, plant 141 phenology controls leaf age, vitality (photosynthetic rates) and foliar density (e.g. Leaf Area Index, 142 LAI)33, and may therefore co-determine ecosystem-level temperature limitations and the optimum 𝑙𝑒𝑎𝑓 143 temperature for canopy photosynthesis34. It is also important to remark when comparing 𝑇𝑜𝑝𝑡 with 𝑒𝑐𝑜 144 𝑇𝑜𝑝𝑡 that leaf-scale measurements are often limited to sunlit leaves, which could lead to a positive 𝑙𝑒𝑎𝑓 145 bias of existing in-situ 𝑇𝑜𝑝𝑡 measurements. Furthermore, the tree species database used by Kattge 𝑙𝑒𝑎𝑓 146 & Knorr5 from which 𝑇𝑜𝑝𝑡 data were collected does not include any tropical species. This may 𝑙𝑒𝑎𝑓 147 explain why global models prescribed with 𝑇𝑜𝑝𝑡 give divergent results for tropical biomes. 148 𝑒𝑐𝑜 149 The relationship between 𝑇𝑜𝑝𝑡 and background climate is shown in Fig. 1d. The sampling of leaf- 𝑙𝑒𝑎𝑓 150 scale studies does not provide consistent evidence about the dependence of 𝑇𝑜𝑝𝑡 on climate, with 𝑙𝑒𝑎𝑓 151 positive correlations between 𝑇𝑜𝑝𝑡 and growing season air temperature in a set of studies1, 5, 35-37 𝑙𝑒𝑎𝑓 152 attributed to evolutionary adaptation38, but no clear relationship between 𝑇𝑜𝑝𝑡 and growth 𝑒𝑐𝑜 153 temperature39-41. In contrast, 𝑇𝑜𝑝𝑡 inferred from satellite GPP proxies in our study increases with 𝑎𝑖𝑟 𝑒𝑐𝑜 154 𝑇𝑚𝑎𝑥 𝑔𝑠 across the globe. In temperature-precipitation space, the spatial sensitivity of 𝑇𝑜𝑝𝑡 to 𝑎𝑖𝑟 155 𝑇𝑚𝑎𝑥 𝑔𝑠 (the slope of the linear regression between these two variables) is lower than 1 for any 𝑒𝑐𝑜 156 precipitation bin (Fig. 1d), suggesting that spatial gradients of 𝑇𝑜𝑝𝑡 are smaller than those of 𝑎𝑖𝑟 𝑒𝑐𝑜 157 𝑇𝑚𝑎𝑥 𝑔𝑠 , possibly because hydraulic and phenological limitations further limit 𝑇𝑜𝑝𝑡 across spatial 𝑒𝑐𝑜 𝑎𝑖𝑟 158 gradients. In fact, the spatial sensitivity of 𝑇𝑜𝑝𝑡 to 𝑇𝑚𝑎𝑥 𝑔𝑠 generally increases with increasing mean 𝑒𝑐𝑜 159 annual precipitation (Fig. 1d), even though 𝑇𝑜𝑝𝑡 is not significantly correlated with precipitation 𝑎𝑖𝑟 𝑒𝑐𝑜 160 after controlling for the effect of 𝑇𝑚𝑎𝑥 𝑔𝑠 (Fig. 1d). This thermal adaptation of 𝑇𝑜𝑝𝑡 , suggested by 𝑒𝑐𝑜 161 the positive spatial slope of the 𝑇𝑜𝑝𝑡 -air temperature relation, is also observed across biomes. As 𝑒𝑐𝑜 𝑎𝑖𝑟 162 shown in Fig. 2, there is a significant positive correlation between 𝑇𝑜𝑝𝑡 and 𝑇𝑚𝑎𝑥 𝑔𝑠 with a slope 𝑒𝑐𝑜 163 of 0.76 across different biomes. Among biomes, the largest mean 𝑇𝑜𝑝𝑡 is found in tropical 𝑒𝑐𝑜 164 evergreen broad-leaved forest (29±3 °C), and the smallest mean 𝑇𝑜𝑝𝑡 (13±3 °C) in cold grasslands 165 covering the Tibetan Plateau (Fig. 2 and Supplementary Fig. 8). 166 167 Both model results and very limited observational studies suggest a decrease in canopy 168 photosynthesis of tropical forests at high temperature15, 42-45 which led us to formulate the third 𝑒𝑐𝑜 𝑎𝑖𝑟 169 hypothesis of tropical forests already operating at 𝑇𝑜𝑝𝑡 being close to 𝑇𝑚𝑎𝑥 𝑔𝑠 implying that canopy 170 photosynthesis may decrease in case of future warming15, 16. This hypothesis is verified from the 𝑒𝑐𝑜 171 data shown in Fig. 3 (see also Supplementary Fig. 9). 𝑇𝑜𝑝𝑡 is indeed slightly (1.4 °C) lower than 𝑎𝑖𝑟 172 𝑇𝑚𝑎𝑥 𝑔𝑠 over tropical evergreen forests, suggesting a small ‘safety margin’ for canopy 173 photosynthesis under future warming. Note that the “safety margin” could become larger than that 174 suggested by the air temperature data if leaf thermal regulation would acclimate to the warming air 175 temperature (see Methods). In contrast, arctic (north of 65oN) and boreal (50oN-65oN) ecosystems 𝑒𝑐𝑜 176 exhibit substantially larger safety margins, i.e. a larger positive difference between 𝑇𝑜𝑝𝑡 and 𝑎𝑖𝑟 𝑎𝑖𝑟 177 𝑇𝑚𝑎𝑥 𝑔𝑠 (Fig. 3a and Supplementary Fig. 9). Analysis of the 16-day averaged 𝑇𝑚𝑎𝑥 distribution 𝑒𝑐𝑜 𝑒𝑐𝑜 𝑎𝑖𝑟 178 during the period when 𝑇𝑜𝑝𝑡 is observed further shows that the rank of 𝑇𝑜𝑝𝑡 in the 𝑇𝑚𝑎𝑥 𝑎𝑖𝑟 179 distribution is already near the highest quantile of 𝑇𝑚𝑎𝑥 (>70%) for tropical evergreen forests 180 (Supplementary Fig. 10). Based on this result, one may expect that rising air temperature in the 181 future, irrespective of the indirect effect of increasing VPD, may limit or decrease vegetation 182 productivity in tropical forests, but not in temperate or boreal ecosystems. 183 184 Global terrestrial daily maximum air temperature is projected to rise by 1.9 °C under the RCP2.6 185 low warming climate scenario and by 5.6 °C under the RCP8.5 scenario by 210046. We compared 𝑎𝑖𝑟 𝑒𝑐𝑜 186 these 𝑇𝑚𝑎𝑥 𝑔𝑠 projections with the present-day distribution of 𝑇𝑜𝑝𝑡 with a focus on tropical 𝑎𝑖𝑟 187 evergreen forests, where optimum temperature is currently just below the limit of 𝑇𝑚𝑎𝑥 𝑔𝑠 (see 𝑒𝑐𝑜 188 Methods; Fig. 3b and 3c). The key uncertainty in this discussion is, however, whether or not 𝑇𝑜𝑝𝑡 𝑎𝑖𝑟 189 will acclimate and follow the increase in 𝑇𝑚𝑎𝑥 𝑔𝑠 . We therefore looked at possible acclimation from 𝑒𝑐𝑜 190 time series of 𝑇𝑜𝑝𝑡 retrieved from the Advanced Very High Resolution Radiometer (AVHRR) 191 NDVI, which spans the last 30 years and comprises almost a 1°C temperature range. NDVI-derived 𝑒𝑐𝑜 192 𝑇𝑜𝑝𝑡 did not have a significant trend over the last three decades except for the northern lands (north 193 of 60oN) where warming is more pronounced47 (Supplementary Fig. 11). This suggests that the 194 recent 1°C warming is not large enough to elicit an acclimation response from some ecosystems, 𝑒𝑐𝑜 195 given decadal variability48. Also the annual 𝑇𝑜𝑝𝑡 derived from flux sites estimates of GPP did not 𝑎𝑖𝑟 196 exhibit a positive trend and was not significantly correlated with annual variations of 𝑇𝑚𝑎𝑥 𝑔𝑠 , 𝑒𝑐𝑜 197 although the flux time series are probably too short to properly evaluate trends of 𝑇𝑜𝑝𝑡 related to 198 possible acclimation processes (Supplementary Fig. 12). Because we detected no indication for its 𝑎𝑖𝑟 199 existence, we first assumed no acclimation in the comparison of future 𝑇𝑚𝑎𝑥 𝑔𝑠 projections from 𝑒𝑐𝑜 𝑎𝑖𝑟 200 climate models with the current distribution of 𝑇𝑜𝑝𝑡 . Under this assumption, the average 𝑇𝑚𝑎𝑥 𝑔𝑠 𝑒𝑐𝑜 201 of tropical evergreen forests will exceed the current value of 𝑇𝑜𝑝𝑡 for RCP2.6 by 2.6°C, and by 5.7 202 °C for RCP8.5 (Fig. 3c). On the other hand, boreal and arctic biomes will still remain within the 𝑒𝑐𝑜 𝑎𝑖𝑟 203 safety margin, with 𝑇𝑜𝑝𝑡 staying above 𝑇𝑚𝑎𝑥 𝑔𝑠 except under the RCP8.5 high warming scenario 204 (Fig. 3b and Supplementary Fig. 13). 205 206 Despite the lack of in situ observational evidence for GPP acclimation to the ongoing warming 207 trend, we tested a simple future acclimation scenario based on the space-for-time substitution 208 approach49 as applied in several studies using observed spatial gradients to hindcast temporal 𝑒𝑐𝑜 𝑎𝑖𝑟 209 changes50,51. Here, we assume that temporal change of 𝑇𝑜𝑝𝑡 will evolve proportionally to 𝑇𝑚𝑎𝑥 𝑔𝑠 , 𝑒𝑐𝑜 𝑎𝑖𝑟 210 following the spatial temperature sensitivity of 𝑇𝑜𝑝𝑡 to 𝑇𝑚𝑎𝑥 𝑔𝑠 in Fig. 1d and the indirect effects 211 of temperature increase (e.g. by increasing VPD) are excluded. We took the differences in 212 precipitation levels into account, so that areas that become wetter also exhibit faster acclimation. 𝑎𝑖𝑟 𝑒𝑐𝑜 213 Even with this assumed acclimation law, 𝑇𝑚𝑎𝑥 𝑔𝑠 will still surpass 𝑇𝑜𝑝𝑡 by 1.7 °C under RCP2.6 214 and by 2.5 °C under RCP8.5 for tropical evergreen forests (Fig. 3c). Not accounting for 215 precipitation levels in the acclimation rates produced similar results (Supplementary Figs. 14 and 216 15). 217 𝑒𝑐𝑜 218 Our global-scale analysis of 𝑇𝑜𝑝𝑡 derived from globally distributed point measurements of eddy 219 covariance and space-borne observations of proxies of vegetation productivity is a first attempt to 220 diagnose the global distribution of ecosystem-scale temperature optima of photosynthesis. It should 𝑒𝑐𝑜 221 be noted, however, that hypotheses about that thermal acclimation of 𝑇𝑜𝑝𝑡 are still highly uncertain, 222 because ecosystem adjustments can lag substantially behind the rate of future warming, particularly 223 for forests. More studies using data sets with longer time spans are needed in the future in order to 𝑒𝑐𝑜 224 more accurately detect eventual thermal acclimation of 𝑇𝑜𝑝𝑡 . Furthermore, the acclimation of 225 plants to increasing atmospheric CO2 concentration and to changes in other environmental factors 226 (e.g. VPD) was also not considered in the current analyses. Constraining the spatially observed 𝑒𝑐𝑜 227 temperature sensitivity of 𝑇𝑜𝑝𝑡 over time is a priority for future studies. Continuous monitoring and 𝑒𝑐𝑜 228 dedicated manipulative experiments could improve our understanding of the features of 𝑇𝑜𝑝𝑡 and 229 thermal acclimation in earth system models52. 230 231 Methods 232 FLUXNET data The half-hourly eddy-covariance Gross Primary Productivity (GPP) data were 233 obtained from FLUXNET datasets, and were quality-controlled, filtered against low turbulence, 234 and gap-filled using consistent methods, as described by Papale et al.53. Only freely available 235 FLUXNET data were used in this study. All the half-hourly GPP data were aggregated into daily 236 accumulated GPP for further estimates of the optimal temperature for vegetation productivity. 𝑎𝑖𝑟 237 Daily maximum air temperature (𝑇𝑚𝑎𝑥 ) was determined as the maximum air temperature value 238 from all the half-hourly air temperature observations. We included only site-years with more than 239 80% of half-hourly data available. A total of 153 individual FLUXNET sites with 663 site-years of 240 GPP data were used in this study. 241 242 Near-infrared reflectance of terrestrial vegetation Badgley et al.19 have recently proposed a new 243 approach for estimating vegetation photosynthetic capacity by remote sensing, i.e. the near-infrared 244 reflectance (NIR) of terrestrial vegetation (NIRV), which can differentiate between the confounding 245 effects of background brightness, leaf area and the distribution of photosynthetic capacity with 246 depth in canopies19. NIRV is calculated as the product of total scene NIR reflectance (NIRT) and 247 Normalized Difference Vegetation Index (NDVI)19. As a proxy of photosynthesis, NIRV is 248 suggested to be strongly correlated with Solar-induced Chlorophyll Fluorescence (SIF), a direct 249 index of photons intercepted by chlorophyll, and shows higher correlation with observed GPP than 250 NDVI19. We used satellite-derived NIRV to calculate and map the optimal air temperature for 𝑒𝑐𝑜 251 vegetation productivity at an ecosystem scale (𝑇𝑜𝑝𝑡 ). Following Badgley et al.19, we calculated 16- 252 day NIRV for 2001-2013 as the product of Moderate Resolution Imaging Spectroradiometer 253 (MODIS) 16-day NIR reflectance and MODIS 16-day NDVI, both of which were derived from the 254 MOD13A2 Vegetation Index Product with a spatial revolution of 1 km. Only positive NIR V values 255 were used in the analysis. 256 257 Normalized Difference Vegetation Index The Normalized Difference Vegetation Index (NDVI) 258 is a vegetation index defined as the ratio of the difference between NIR and red visible reflectance 259 to their sum, and is widely used to represent vegetation greenness54. To account for uncertainties 260 from different satellite datasets, three independent NDVI datasets were utilized, including biweekly 261 NDVI data from Global Inventory Modeling and Mapping Studies (GIMMS) Advanced Very High 262 Resolution Radiometer (AVHRR), 16-day NDVI data from Terra Moderate Resolution Imaging 263 Spectroradiometer (MODIS) and 10-day NDVI data from Satellite Pour l'Observation de la Terre 264 Vegetation (SPOT Vegetation). The three NDVI datasets spanned the last three decades: 1982- 265 2009 for AVHRR NDVI datasets, 2000-2009 for MODIS NDVI datasets, and 1999-2009 for SPOT 266 NDVI datasets, with the spatial resolutions of 8 km, 1 km, and 1 km, respectively. All NDVI 267 datasets have been corrected to reduce the effects of volcanic aerosols, solar angle, and sensor 268 errors20,55,56. Pixels with a mean annual NDVI>0.1 were defined as the vegetated area for each 269 dataset. 270 271 Enhanced Vegetation Index The Enhanced Vegetation Index (EVI) is another vegetation index 272 designed to enhance the vegetation signal by minimizing canopy-soil variations and to improve 273 sensitivity over dense vegetation conditions21, and is found to correlate well with estimated GPP 274 on a site-by-site basis57. We used a 16-day EVI dataset for the period of 2000-2009 with a spatial 275 resolution of 1 km from the MOD12A1 Vegetation Index Product. Effects from aerosols, solar 276 angle and sensor error have all been corrected21. 277 278 Sun-induced Chlorophyll Fluorescence Chlorophylls in plant absorbs shortwave radiation and 279 dissipates excess energy as light or heat. The longwave radiation re-emitted by chlorophylls is 280 referred as chlorophyll fluorescence. Recent studies have reported that remotely sensed Sun- 281 induced Chlorophyll Fluorescence (SIF) could serve as an indicator of photosynthesis rate and it is 282 well correlated with model-simulated GPP58. Following the previous studies58,59, we retrieved SIF 283 from two different retrieval windows, 757 nm and 771 nm, as well as the two polarization states, S 284 and P using observation of Fourier Transform Spectrometer (FTS) on the Japanese Greenhouse 285 gases Observing SATellite (GOSAT)20. These diverse SIF samples of were then aggregated into 286 monthly gridded data at a spatial resolution of 2o from June 2009 to June 2012. 287 288 Vegetation distribution We used MODIS land cover with the classification scheme of the 289 International Geosphere-Biosphere Programme (IGBP). The MODIS IGBP land cover data were 290 derived from the MOD12Q1 Land Cover Science Data Product at a spatial resolution of 1 km and 291 an updated digital Köppen-Geiger world map of climatic classification60. Within the vegetated area 292 defined by NDVI thresholds, the 17 land cover types were reclassified into 9 vegetation types, 293 namely evergreen needle-leaved forest (ENF), evergreen broad-leaved forest (EBF), deciduous 294 needle-leaved forest (DNF), deciduous broad-leaved forest (DBF), mixed forest (MF), savannas, 295 cropland, grassland and shrubland. Based on the main climates in the world Map of the Köppen- 296 Geiger climatic classification60, grassland was further subdivided into temperate grasslands, boreal 297 and arctic tundra, and shrubland was further subdivided into temperate and boreal shrubland. The 298 grassland over the Tibetan Plateau was considered separately considering the fact that Tibetan 299 Plateau has an average altitude higher than 4000 m a.s.l.61, and thus an unique alpine climate. In 300 contrast to temperate grasslands / shrubland where water is major limiting factor for vegetation 301 productivity, alpine ecosystems on the Tibetan Plateau are mainly limited by thermal conditions62. 302 303 Climate dataset The gridded air temperature and precipitation data for the period of 1982-2013 304 were obtained from the CRU/NCEP 6-hourly dataset with a spatial resolution of 0.5o. Note that the 305 purpose of this study is to investigate the optimal air temperature for photosynthesis. Optimal leaf 306 temperature is also interesting but not addressed in this study because accurate canopy-integrated 307 measurements of leaf temperatures are not available at the eddy covariance sites and at global scale 308 as gridded datasets. For a discussion about calculation of temperature optimum from air 309 temperature and from surface temperature, we used the remotely sensed land surface temperature 310 (LST), which is inversed from infra-red emissivity measured by MODIS (MYD11A2 version 6). 311 This dataset had an original spatial resolution of 1 km, spanning from July 2002 to December 2014. 312 The error of MODIS LST product, which primarily stems from cloud contamination and emissivity 𝑠𝑢𝑟𝑓𝑎𝑐𝑒 313 uncertainties, was reported to be less than 3°C63. Generally, the occurrence time of 𝑇𝑚𝑎𝑥 (2:00- 314 4:00 P.M.) is relatively close to the Aqua overpass time (1:30 P.M.); thus, we assumed that 𝑠𝑢𝑟𝑓𝑎𝑐𝑒 315 𝑇𝑚𝑎𝑥 from MODIS-Aqua is comparable with the daily maximum leaf surface temperature 𝑙𝑒𝑎𝑓 316 (𝑇𝑚𝑎𝑥 ). Corresponding to the temporal resolutions of MODIS, AVHRR and SPOT datasets, 317 respectively, the 6-hourly climate data were aggregated into 16-day, biweekly, and 10-day values 318 before further analyses. Given the different spatial resolutions of satellite observations and climate 319 data, we simply extracted time series of daily maximum air temperature and precipitation from the 320 aggregated CRU/NCEP data for each pixel of the sets of remotely sensed data. The daily maximum 𝑎𝑖𝑟 321 air temperature (𝑇𝑚𝑎𝑥 ) of the growing season averaged over 2001-2013 was calculated as the 𝑎𝑖𝑟 322 current mean growing-season daily maximum air temperature ( 𝑇𝑚𝑎𝑥 𝑔𝑠 ). Information on the 323 growing season was derived from the study by Zhu et al.64, which was first determined from the 324 GIMMS Leaf Area Index data set (GIMMS LAI3g) using a Savitzky-Golay filter and then refined 325 by excluding the ground-freeze period identified by the Freeze/Thaw Earth System Data Record 326 (see details in ref 64). We also documented the temperature thresholds at which the growing season 327 begins and ends for each year. Temperature thresholds were averaged over 2001-2013 for the onset 328 and end of the growing season, respectively. We also used WATCH Forcing Data Methodology to 329 ERA-Interim data with a temporal resolution of three hours (WFDEI)65. 330 331 We also used climate projections for the end of the 21st century (2091-2100) using 20 models that 332 participated in the phase five of Coupled Model Intercomparison Project (CMIP5) under the 333 RCP2.6, RCP4.5 and RCP8.5 scenarios46 to determine the impact of future warming on vegetation 334 productivity (see model list in Supplementary Table 2). Considering the mismatch between 335 CRU/NCEP datasets and outputs from General Circulation Models (GCMs) for current climate 336 conditions, we generated future temperature and precipitation maps by adding the relative changes 𝑎𝑖𝑟 337 in GCM-derived climate projections to the current climate for each pixel. 𝑇𝑚𝑎𝑥 𝑔𝑠 for the late 21 st 𝑎𝑖𝑟 338 century was estimated using the same temperature thresholds as for the current 𝑇𝑚𝑎𝑥 𝑔𝑠 . All GCM 339 projections were resampled to a resolution of 1o using first-order conservative interpolation 340 method66. 341 𝑒𝑐𝑜 342 Analysis We estimated local 𝑇𝑜𝑝𝑡 by examining the temperature response curve of MODIS NIRV. 343 Following Yuan et al.37 and Niu et al.18, NIRV time series throughout the entire monitoring period 344 and the corresponding temperature data were grouped into 1oC-temperature bins for each pixel 345 within vegetated areas, which were defined as regions with a mean annual NDVI value larger than 346 0.1. We used the 90% quantile of the NIRV data as the response of NIRV within each temperature 347 bin due to the potential influences of other environmental constraints like clouds and droughts. We 348 next calculated the running means of every three temperature bins to develop the temperature 𝑒𝑐𝑜 349 response curve of NIRV. The 𝑇𝑜𝑝𝑡 was then determined from the response curve at which NIRV 𝑒𝑐𝑜 350 was maximized (Supplementary Fig. 16). Note that 𝑇𝑜𝑝𝑡 may not be detected for some pixels where 351 the maximum NIRV was only attained at either end of the response curve, accounting for 3.5% of 𝑒𝑐𝑜 352 the vegetated areas. Only vegetated areas with detectable 𝑇𝑜𝑝𝑡 were shown when mapping the 𝑒𝑐𝑜 𝑒𝑐𝑜 353 spatial pattern of 𝑇𝑜𝑝𝑡 . The derivation of 𝑇𝑜𝑝𝑡 is robust to the choice of a particular climate-forcing 354 dataset (Supplementary Fig. 2). Instead of using the temperature corresponding to the maximum 𝑒𝑐𝑜 355 90th quantile NIRV to calculate 𝑇𝑜𝑝𝑡 , we also applied nonlinear regression of the photosynthetic 𝑒𝑐𝑜 356 temperature response data (Eq.1) to estimate 𝑇𝑜𝑝𝑡 , which produced similar results (Supplementary 357 Fig. 2): 𝑒𝑐𝑜 2 358 NIRV(T)= NIRV(OPT) - b(T-𝑇𝑜𝑝𝑡 ) (Eq.1) 359 where NIRV(T) is the NIRV value at a daily maximum temperature T, b is a parameter describing 𝑒𝑐𝑜 360 the spread of the parabola48,67. 𝑇𝑜𝑝𝑡 is the vertex of each fit and NIRV(OPT) is the NIRV value at 𝑒𝑐𝑜 𝑎𝑖𝑟 𝑎𝑖𝑟 𝑒𝑐𝑜 361 𝑇𝑜𝑝𝑡 . Finally, we used daily mean air temperature (𝑇𝑚𝑒𝑎𝑛 ) instead of 𝑇𝑚𝑎𝑥 to calculate 𝑇𝑜𝑝𝑡 . In this 𝑒𝑐𝑜 𝑎𝑖𝑟 𝑒𝑐𝑜 𝑎𝑖𝑟 362 test, 𝑇𝑜𝑝𝑡 derived from 𝑇𝑚𝑒𝑎𝑛 is smaller than 𝑇𝑜𝑝𝑡 estimated from 𝑇𝑚𝑎𝑥 , but the two variables were 363 strongly spatially correlated (Supplementary Fig. 6). 364 𝑒𝑐𝑜 𝑒𝑐𝑜 365 We investigated the relationship between 𝑇𝑜𝑝𝑡 and climate variables by averaging 𝑇𝑜𝑝𝑡 in the 𝑎𝑖𝑟 𝑎𝑖𝑟 366 climate space with 1-°C intervals of mean annual 𝑇𝑚𝑎𝑥 averaged over the growing season (𝑇𝑚𝑎𝑥 𝑔𝑠 ) 367 and 100-mm intervals of mean annual precipitation (MAP) (Fig. 1d). For each MAP interval, we 𝑒𝑐𝑜 𝑎𝑖𝑟 368 calculated the ‘apparent’ spatial sensitivity of 𝑇𝑜𝑝𝑡 in response to changes in 𝑇𝑚𝑎𝑥 𝑔𝑠 using 369 bootstrapping method. We performed the linear regression analysis 1000 times by randomly 𝑒𝑐𝑜 𝑎𝑖𝑟 370 selecting a subset of 80% of the samples from pairs of 𝑇𝑜𝑝𝑡 and 𝑇𝑚𝑎𝑥 𝑔𝑠 within each MAP interval. 𝑒𝑐𝑜 371 The mean and SD of the 1000 temperature sensitivity of 𝑇𝑜𝑝𝑡 were subsequently estimated along 372 the MAP gradient. 373 374 Air temperature, atmospheric vapor pressure deficit (VPD) and solar radiation usually co-vary in 𝑒𝑐𝑜 375 time and space, so that the empirical observation of spatial patterns of 𝑇𝑜𝑝𝑡 in this study cannot be 𝑒𝑐𝑜 376 unambiguously attributed to air temperature as a single explaining factor of 𝑇𝑜𝑝𝑡 . Under conditions 377 of high temperature, atmospheric VPD increases, soil moisture decreases with a lag, stomatal 378 conductance and hence carbon assimilation rates (GPP at the ecosystem-scale) decrease to prevent 379 exceedingly low leaf water potentials and resulting plant tissue damage from cavitation24. We show 𝑒𝑐𝑜 380 that across climatic gradients 𝑇𝑜𝑝𝑡 is systematically higher at high maximum air temperatures but 381 not systematically lower at high VPD conditions (Supplementary Fig. 17). Then we calculated the 𝑎𝑖𝑟 382 variance inflation factor (VIF) between VPD and 𝑇max 𝑔𝑠 under each VPD bin in the regression 383 model of: 𝑒𝑐𝑜 𝑎𝑖𝑟 384 𝑇𝑜𝑝𝑡 = 𝑘0 + 𝑘1 × 𝑇𝑚𝑎𝑥 𝑔𝑠 + 𝑘2 × 𝑉𝑃𝐷 (Eq.2) 385 As shown in Supplementary Figure 18, we observed that VIF value ranged only between 1.001 and 386 1.438, suggesting relatively low multicollinearity between VPD and temperature. Even so, to 𝑒𝑐𝑜 𝑎𝑖𝑟 387 examine whether VPD can substantially affect the relationship between 𝑇𝑜𝑝𝑡 and 𝑇𝑚𝑎𝑥 𝑔𝑠 , we 𝑒𝑐𝑜 𝑎𝑖𝑟 388 further calculated the partial (‘intrinsic’) sensitivity of 𝑇𝑜𝑝𝑡 to 𝑇𝑚𝑎𝑥 𝑔𝑠 in each grid point based on 389 the following bilinear regression: 𝑒𝑐𝑜 𝑎𝑖𝑟 𝑎𝑖𝑟 390 𝑇𝑜𝑝𝑡 = 𝑘0 + 𝑘1 × 𝑇𝑚𝑎𝑥 𝑔𝑠 + 𝑘2 × 𝑉𝑃𝐷 + 𝑘3 × 𝑉𝑃𝐷 × 𝑇𝑚𝑎𝑥 𝑔𝑠 (Eq.3) 𝑒𝑐𝑜 𝑎𝑖𝑟 391 where the partial sensitivity of 𝑇𝑜𝑝𝑡 to 𝑇𝑚𝑎𝑥 𝑔𝑠 is defined as 𝑘1 in Eq. 3 under each VPD bin. Then 𝑒𝑐𝑜 𝑎𝑖𝑟 392 we compared the partial sensitivity with the apparent sensitivity of 𝑇𝑜𝑝𝑡 to 𝑇𝑚𝑎𝑥 𝑔𝑠 estimated using 𝑒𝑐𝑜 𝑎𝑖𝑟 393 abovementioned linear regression between 𝑇𝑜𝑝𝑡 and 𝑇𝑚𝑎𝑥 𝑔𝑠 for each VPD bin. As shown in 𝑒𝑐𝑜 𝑎𝑖𝑟 394 Supplementary Fig. 19, although the apparent sensitivity of 𝑇𝑜𝑝𝑡 to 𝑇𝑚𝑎𝑥 𝑔𝑠 is generally lower than 𝑒𝑐𝑜 𝑎𝑖𝑟 𝑎𝑖𝑟 395 the partial (‘intrinsic’) sensitivity of 𝑇𝑜𝑝𝑡 to 𝑇𝑚𝑎𝑥 𝑔𝑠 the apparent sensitivity to 𝑇𝑚𝑎𝑥 𝑔𝑠 remains 396 positive even when VPD is taken into account excepted under very high VPD bins (higher than 397 ~4.5 kPa) representing less than 1% of the study area. These results indicate that the patterns of 𝑒𝑐𝑜 398 𝑇𝑜𝑝𝑡 are not dominated by high VPD reducing canopy photosynthesis, as an indirect effect of 399 higher air temperature increasing VPD. Moreover, we also calculated the percentiles of downward 𝑒𝑐𝑜 400 shortwave solar radiation (Rad) at the time of year when 𝑇𝑜𝑝𝑡 is observed for the 16-day averaged 𝑒𝑐𝑜 401 Rad distribution. As shown in Supplementary Fig. 20, the Rad value when 𝑇𝑜𝑝𝑡 was retrieved from 402 global observations were below the 95th percentile in the 16-day Rad distribution for ~80% of the 403 study area, which is mainly in mid and low latitudes such as Africa, India, Australia, eastern Brazil 404 and southern and southwestern of North America. By comparison, for most part of boreal regions, 𝑒𝑐𝑜 405 part of south China, southeastern US, as well as part of South America, the timing of 𝑇𝑜𝑝𝑡 is 𝑒𝑐𝑜 406 consistent the time of maximum solar radiation. This is because 𝑇𝑜𝑝𝑡 in these regions generally 407 appears in summer, which is also the period when solar radiation displays its maximum during the 408 year. 409 𝑒𝑐𝑜 𝑒𝑐𝑜 410 The NIRV-derived 𝑇𝑜𝑝𝑡 was compared with 𝑇𝑜𝑝𝑡 estimated using GPP data from 153 eddy 𝑒𝑐𝑜 411 covariance sites. Flux-derived 𝑇𝑜𝑝𝑡 was determined for each site-year with daily-accumulated GPP 412 and corresponding temperature data from flux tower observations, applying the same method to 𝑒𝑐𝑜 𝑒𝑐𝑜 413 estimate local 𝑇𝑜𝑝𝑡 using NIRV datasets. A robust estimate of 𝑇𝑜𝑝𝑡 can be derived for 125 sites 𝑒𝑐𝑜 414 (Supplementary Table 1). For each site, we calculated the mean and SD of 𝑇𝑜𝑝𝑡 across different 𝑒𝑐𝑜 415 years. Then, we extracted and averaged 𝑇𝑜𝑝𝑡 values within a 3×3 pixel window around each site 𝑒𝑐𝑜 𝑒𝑐𝑜 416 from the NIRV-derived 𝑇𝑜𝑝𝑡 map, and calculated the SD of the nine 𝑇𝑜𝑝𝑡 values within the window. 𝑒𝑐𝑜 417 The relationship between NIRV- and flux-derived 𝑇𝑜𝑝𝑡 was reported using a least square linear 418 regression, and the statistical significance of the slope, or its p-value, given by Student’s t test. The 𝑒𝑐𝑜 419 results show that NIRV-derived 𝑇𝑜𝑝𝑡 is comparable with that estimated independently from 420 measurements of flux-tower eddy covariance (Fig. 1b). 421 𝑒𝑐𝑜 422 We compared the spatial distribution of 𝑇𝑜𝑝𝑡 derived from NIRV with the one obtained from NDVI 𝑒𝑐𝑜 423 datasets. Consistent spatial patterns of 𝑇𝑜𝑝𝑡 are derived from each of the three NDVI datasets 𝑒𝑐𝑜 424 (Supplementary Fig. 21). A global composite map of 𝑇𝑜𝑝𝑡 (Supplementary Fig. 3) was then 425 generated by averaging over estimates derived from the three NDVI datasets. Given the 𝑒𝑐𝑜 426 inconsistent spatial resolutions of the different products, we resampled 𝑇𝑜𝑝𝑡 to a common grid of 𝑒𝑐𝑜 427 8 km before averaging. 𝑇𝑜𝑝𝑡 from NDVI datasets generally show a spatial pattern similar to that 𝑒𝑐𝑜 428 from NIRV, but with smaller NDVI-derived 𝑇𝑜𝑝𝑡 values for central Australia and southern South 𝑒𝑐𝑜 429 America (Supplementary Fig. 3). We compared the spatial distribution of 𝑇𝑜𝑝𝑡 derived from NIRV 𝑒𝑐𝑜 430 with that from MODIS EVI data during 2001-2013, and found that the EVI derived 𝑇𝑜𝑝𝑡 shows 𝑒𝑐𝑜 431 very similar spatial pattern to that of NIRV derived 𝑇𝑜𝑝𝑡 (Supplementary Fig. 4). The distribution 𝑒𝑐𝑜 432 of 𝑇𝑜𝑝𝑡 derived from NIRV and from GOSAT SIF datasets also have similar spatial patterns, even 𝑒𝑐𝑜 433 though the NIRV-derived 𝑇𝑜𝑝𝑡 is higher in tropical regions, particularly in cultivated areas of 434 southeastern Brazil (Supplementary Fig. 5). 435 436 At leaf scale, the photosynthesis-temperature response is suggested to be primarily controlled by 437 three sets of processes, namely biochemical, respiratory and stomatal processes68. Much of the 438 effort to date to understand variability in the leaf-level photosynthesis-temperature response has 439 focused on biochemical processes68, with Vcmax and Jmax being two major parameters controlling 440 the maximum rates of photosynthesis limited by CO2 and light, respectively69. Therefore, in this 𝑒𝑐𝑜 441 study, we compared 𝑇𝑜𝑝𝑡 derived using GPP proxies with leaf-scale optimal temperature of 442 maximum Rubisco-limited carboxylation rates (Vcmax), although GPP is in theory more comparable 443 to net photosynthesis, that is, leaf gross photosynthesis minus photorespiration and minus dark 444 respiration. Since photorespiration increases exponentially with temperature70, the optimum 𝑒𝑐𝑜 445 temperature of GPP (𝑇𝑜𝑝𝑡 ) should, in theory, be lower than the optimal temperature of maximum 446 Rubisco-limited carboxylation rates (Vcmax). For this comparison to be made, we extracted and 𝑒𝑐𝑜 𝑒𝑐𝑜 447 averaged 𝑇𝑜𝑝𝑡 values within a 3×3 pixel window from the NIRV-derived 𝑇𝑜𝑝𝑡 map around the 448 reported site location (longitude and latitude) of leaf-scale measurements. For leaf-scale 449 measurements without the information of site location, we calculated the average NIRV-derived 𝑒𝑐𝑜 450 𝑇𝑜𝑝𝑡 values across pixels with both the same growing season mean temperature and the same plant 451 functional type as the corresponding site. 452 𝑒𝑐𝑜 𝑙𝑒𝑎𝑓 453 𝑇𝑜𝑝𝑡 is different from 𝑇𝑜𝑝𝑡 not only because of respiratory process, but also because air 454 temperature can differ from leaf temperatures71, which are regulated by leaf traits affecting the leaf 455 energy balance72. Because, to our knowledge, global gridded monthly leaf temperature data are not 𝑠𝑢𝑟𝑓𝑎𝑐𝑒 456 available, we use daily maximum land surface temperature (𝑇𝑚𝑎𝑥 ) from MODIS to calculate 𝑒𝑐𝑜 𝑒𝑐𝑜 𝑒𝑐𝑜 457 𝑇𝑜𝑝𝑡_𝐿𝑆𝑇 in order to illustrate the potential differences between 𝑇𝑜𝑝𝑡_𝐿𝑆𝑇 and 𝑇𝑜𝑝𝑡 . As shown in 𝑒𝑐𝑜 𝑒𝑐𝑜 458 Supplementary Figure 22, the 𝑇𝑜𝑝𝑡_𝐿𝑆𝑇 is similar to 𝑇𝑜𝑝𝑡 over tropical savannas. However, over 𝑒𝑐𝑜 𝑒𝑐𝑜 459 moist tropical forests 𝑇𝑜𝑝𝑡_𝐿𝑆𝑇 is lower than 𝑇𝑜𝑝𝑡 , which can be explained by the lower daytime 460 surface temperature than air temperature as a result of strong evapotranspiration effects71,73. This 𝑒𝑐𝑜 𝑒𝑐𝑜 461 ecosystem-dependent difference between 𝑇𝑜𝑝𝑡_𝐿𝑆𝑇 and 𝑇𝑜𝑝𝑡 suggests that the leaf thermal 462 regulation mechanism through the physiological and morphological changes72 is an important 𝑒𝑐𝑜 463 ecosystem process shaping spatial variations of 𝑇𝑜𝑝𝑡 . In addition, if the difference between leaf 464 temperature and air temperature would become larger in response to warmer air temperatures (i.e. 465 if leaf thermal regulation acclimates to warmer temperature), the “safety margin” of tropical 466 ecosystems would become larger than the air temperature data are currently suggesting. However, 467 the long-term in-situ leaf temperature data required to test this hypothesis independently are 468 currently not available. 469 𝑒𝑐𝑜 𝑒𝑐𝑜 470 To account for potential changes in 𝑇𝑜𝑝𝑡 under future warming, we estimated the acclimated 𝑇𝑜𝑝𝑡 471 for vegetation productivity by the end of the 21st century (2091-2100), using recent IPCC climate 472 projections46. To this end, we applied the space-for-time substitution approach49, assuming that that 𝑒𝑐𝑜 𝑎𝑖𝑟 473 temporally 𝑇𝑜𝑝𝑡 will evolve proportionally to 𝑇𝑚𝑎𝑥 𝑔𝑠 following the spatial temperature sensitivity 𝑒𝑐𝑜 𝑎𝑖𝑟 474 of 𝑇𝑜𝑝𝑡 to 𝑇𝑚𝑎𝑥 𝑔𝑠 . Given the relatively large uncertainties of precipitation projections, we 𝑒𝑐𝑜 475 considered two future precipitation scenarios. For the first scenario, we estimated acclimated 𝑇𝑜𝑝𝑡 𝑒𝑐𝑜 476 pixel by pixel using the temperature sensitivity of 𝑇𝑜𝑝𝑡 under the present MAP level, assuming that 477 MAP does not change between the end of the 21st century. For the second scenario, we accounted 𝑒𝑐𝑜 478 for in MAP, and the acclimated 𝑇𝑜𝑝𝑡 was calculated pixel by pixel using the temperature sensitivity 𝑒𝑐𝑜 𝑒𝑐𝑜 479 of 𝑇𝑜𝑝𝑡 under the projected MAP level for 2091-2100. Acclimated 𝑇𝑜𝑝𝑡 was averaged across the 𝑒𝑐𝑜 480 GCMs under each scenario. Latitudinal variation of future 𝑇𝑜𝑝𝑡 was derived by averaging within 𝑒𝑐𝑜 𝑎𝑖𝑟 481 1°-latitude bins from future 𝑇𝑜𝑝𝑡 maps and then compared with that in future 𝑇𝑚𝑎𝑥 𝑔𝑠 summarized 𝑎𝑖𝑟 482 by latitude from future 𝑇𝑚𝑎𝑥 𝑔𝑠 maps. 483 484 Data Availability 485 All data is available in the main text or the supplementary information. All computer codes used in 486 this study can be provided by the corresponding author upon reasonable requests. 487 488 References 489 1. Berry, J. & Bjorkman, O. Photosynthetic response and adaptation to temperature in higher 490 plants. Annual Review of Plant Physiology 31, 491-543 (1980). 491 2. Hughes, L. Biological consequences of global warming: is the signal already apparent? 492 Trends in Ecology & Evolution 15, 56-61 (2000). 493 3. Niu, S. et al. Climatic warming changes plant photosynthesis and its temperature 494 dependence in a temperate steppe of northern China. Environmental and Experimental 495 Botany 63, 91-101 (2008). 496 4. Way, D. A., Yamori, W. Thermal acclimation of photosynthesis: on the importance of 497 adjusting our definitions and accounting for thermal acclimation of respiration. 498 Photosynthesis research 119, 89-100 (2014). 499 5. Kattge, J., Knorr, W. Temperature acclimation in a biochemical model of photosynthesis: 500 a reanalysis of data from 36 species. Plant, cell & environment 30, 1176-1190 (2007). 501 6. Lloyd, J., Farquhar, G. D. Effects of rising temperatures and [CO2] on the physiology of 502 tropical forest trees. Philosophical Transactions of the Royal Society B: Biological Sciences 503 363, 1811-1817 (2008). 504 7. Medlyn, B. et al. Temperature response of parameters of a biochemically based model of 505 photosynthesis. II. A review of experimental data. Plant, Cell & Environment 25, 1167- 506 1179 (2002). 507 8. Field, C. B., Randerson, J. T., Malmström, C. M. Global net primary production: combining 508 ecology and remote sensing. Remote sensing of Environment 51, 74-88 (1995). 509 9. Niinemets, Ü. Global-scale climatic controls of leaf dry mass per area, density, and 510 thickness in trees and shrubs. Ecology 82, 453-469 (2001). 511 10. Williams, A. P. et al. Temperature as a potent driver of regional forest drought stress and 512 tree mortality. Nature Climate Change 3, 292-297 (2013). 513 11. Tyree, M. T. & Dixon, M. A. Water stress induced cavitation and embolism in some woody 514 plants. Physiologia Plantarum 66, 397-405 (1986). 515 12. Yin, X., Struik, P. C3 and C4 photosynthesis models: an overview from the perspective of 516 crop modelling. NJAS-Wageningen Journal of Life Sciences 57, 27-38 (2009). 517 13. N. G. Smith, S. L. Malyshev, E. Shevliakova, J. Kattge, J. S. Dukes, Foliar temperature 518 acclimation reduces simulated carbon sensitivity to climate. Nature Climate Change 6, 407 519 (2016). 520 14. Mercado, L. M. et al. Large sensitivity in land carbon storage due to geographical and 521 temporal variation in the thermal response of photosynthetic capacity. New Phytologist 218, 522 1462-1477 (2018). 523 15. Corlett, R. T. Impacts of warming on tropical lowland rainforests. Trends in Ecology & 524 Evolution 26, 606-613 (2011). 525 16. Doughty, C. E., Goulden, M. L. Are tropical forests near a high temperature threshold? 526 Journal of Geophysical Research: Biogeosciences 113, G00B07 (2008). 527 17. Baldocchi, D., Falge, E., Gu, L. et al. FLUXNET: A new tool to study the temporal and 528 spatial variability of ecosystem-scale carbon dioxide, water vapor, and energy flux densities. 529 Bulletin of the American Meteorological Society 82, 2415-2434 (2001). 530 18. Niu et al. Thermal optimality of net ecosystem exchange of carbon dioxide and underlying 531 mechanisms. New Phytologist 194, 775-783 (2012). 532 19. Badgley, G., Field, C. B., Berry, J. A. Canopy near-infrared reflectance and terrestrial 533 photosynthesis. Science Advances 3, e1602244 (2017). 534 20. Pinzon, J. E., Tucker, C. J. A non-stationary 1981-2012 AVHRR NDVI3g time series. 535 Remote Sensing 6, 6929-6960 (2014). 536 21. Huete, A. et al. Overview of the radiometric and biophysical performance of the MODIS 537 vegetation indices. Remote Sensing of Environment 83, 195-213 (2002). 538 22. Joiner, J., Yoshida, Y., Vasilkov, A., Middleton, E. First observations of global and 539 seasonal terrestrial chlorophyll fluorescence from space. Biogeosciences 8, 637-651 (2011). 540 23. Ma, S., Osuna, J. L., Verfaillie, J. et al. Photosynthetic responses to temperature across leaf- 541 canopy-ecosystem scales: a 15-year study in a Californian oak-grass savanna. 542 Photosynthesis research 132, 277-291 (2017). 543 24. Stocker, B. D. et al. Quantifying soil moisture impacts on light use efficiency across biomes. 544 New Phytologist 218, 1430-1449 (2018). 545 25. Bunce, J. A. Acclimation of photosynthesis to temperature in eight cool and warm climate 546 herbaceous C3 species: temperature dependence of parameters of a biochemical 547 photosynthesis model. Photosynthesis Research 63, 59-67 (2000). 548 26. Han, Q., Kawasaki, T., Nakano, T. & Chiba, Y. Spatial and seasonal variability of 549 temperature responses of biochemical photosynthesis parameters and leaf nitrogen content 550 within a Pinus densiflora crown. Tree physiology 24, 737-744 (2004). 551 27. Harley, P., Thomas, R., Reynolds, J. & Strain, B. Modelling photosynthesis of cotton grown 552 in elevated CO2. Plant, Cell & Environment 15, 271-282 (1992). 553 28. Onoda, Y., Hikosaka, K. & Hirose, T. The balance between RuBP carboxylation and RuBP 554 regeneration: a mechanism underlying the interspecific variation in acclimation of 555 photosynthesis to seasonal change in temperature. Functional Plant Biology 32, 903-910 556 (2005). 557 29. Walcroft, A., Le Roux, X., Diaz-Espejo, A., Dones, N. & Sinoquet, H. Effects of crown 558 development on leaf irradiance, leaf morphology and photosynthetic capacity in a peach 559 tree. Tree Physiology 22, 929-938 (2002). 560 30. Wang, K.Y., Kellomäki, S. & Laitinen, K. Acclimation of photosynthetic parameters in 561 Scots pine after three years exposure to elevated temperature and CO2. Agricultural and 562 Forest Meteorology 82, 195-217 (1996). 563 31. Vaz, M. et al. Drought-induced photosynthetic inhibition and autumn recovery in two 564 Mediterranean oak species (Quercus ilex and Quercus suber). Tree Physiology 30, 946-956 565 (2010). 566 32. Zhou, L. et al. Responses of photosynthetic parameters to drought in subtropical forest 567 ecosystem of China. Scientific reports 5, 18254 (2015). 568 33. Doughty, C. E., Goulden, M. L. Seasonal patterns of tropical forest leaf area index and CO2 569 exchange. Journal of Geophysical Research: Biogeosciences 113, G00B06, 570 doi:10.1029/2007JG000590 (2008). 571 34. Wu, J. et al. Partitioning controls on Amazon forest photosynthesis between environmental 572 and biotic factors at hourly to interannual timescales. Global Change Biology, 23, 1240- 573 1257 (2017). 574 35. Gunderson, C. A., O'Hara, K. H., Campion, C. M., Walker, A. V., Edwards, N. T. Thermal 575 plasticity of photosynthesis: the role of acclimation in forest responses to a warming climate. 576 Global Change Biology 16, 2272-2286 (2010). 577 36. Mooney, H. A., Björkman, O., Collatz, G. J. Photosynthetic acclimation to temperature in 578 the desert shrub, Larrea divaricata I. Carbon dioxide exchange characteristics of intact 579 leaves. Plant Physiology 61, 406-410 (1978). 580 37. Yuan, W. et al. Thermal adaptation of net ecosystem exchange. Biogeosciences 8, 1453- 581 1463 (2011). 582 38. Sage, R. F., Kubien, D. S. The temperature response of C3 and C4 photosynthesis. Plant, 583 Cell & Environment 30, 1086-1106 (2007). 584 39. Armond, P. A., Schreiber, U. & Björkman, O. Photosynthetic acclimation to temperature 585 in the desert shrub, Larrea divaricata: II. Light-harvesting efficiency and electron transport. 586 Plant Physiology 61, 411-415 (1978). 587 40. Badger, M., Björkman, O. & Armond, P. An analysis of photosynthetic response and 588 adaptation to temperature in higher plants: temperature acclimation in the desert evergreen 589 Nerium oleander L. Plant, Cell & Environment 5, 85-99 (1982). 590 41. Atkin, O., Scheurwater, I. & Pons, T. High thermal acclimation potential of both 591 photosynthesis and respiration in two lowland Plantago species in contrast to an alpine 592 congeneric. Global Change Biology 12, 500-515 (2006). 593 42. Doughty, C. E. An in situ leaf and branch warming experiment in the Amazon. Biotropica 594 43, 658-665 (2011). 595 43. Koch, G. W., Amthor, J. S., Goulden, M. L. Diurnal patterns of leaf photosynthesis, 596 conductance and water potential at the top of a lowland rain forest canopy in Cameroon: 597 measurements from the Radeau des Cimes. Tree Physiology 14, 347-360 (1994). 598 44. Tribuzy, E. S. Variações da temperatura foliar do dossel eo seu efeito na taxa assimilatória 599 de CO2 na Amazônia Central, Universidade de São Paulo (2005). 600 45. White, A., Cannell, M. G., Friend, A. D. CO2 stabilization, climate change and the terrestrial 601 carbon sink. Global Change Biology 6, 817-833 (2000). 602 46. IPCC, Climate Change 2013: The Physical Science Basis. Contribution of Working Group 603 I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. 604 (Cambridge University Press, Cambridge, UK, 2013). 605 47. Xu, L. et al. Temperature and vegetation seasonality diminishment over northern lands. 606 Nature Climate Change 3, 581 (2013). 607 48. Sendall, K. M. et al. Acclimation of photosynthetic temperature optima of temperate and 608 boreal tree species in response to experimental forest warming. Global Change Biology 21, 609 1342-1357 (2015). 610 49. Likens, G. E. Long-term Studies in Ecology. (Springer, 1989). 611 50. Jung, M. et al. Global patterns of land-atmosphere fluxes of carbon dioxide, latent heat, and 612 sensible heat derived from eddy covariance, satellite, and meteorological observations. 613 Journal of Geophysical Research: Biogeosciences 116, G00J07 (2011). 614 51. Lester, R. E., Close, P. G., Barton, J. L., Pope, A. J. & Brown, S. C. Predicting the likely 615 response of data-poor ecosystems to climate change using space-for-time substitution 616 across domains. Global Change Biology 20, 3471-3481 (2014). 617 52. Cavaleri, M. A., Reed, S. C., Smith, W. K., Wood, T. E. Urgent need for warming 618 experiments in tropical forests. Global Change Biology 21, 2111–2121 (2015). 619 53. Papale, D. et al. Towards a standardized processing of Net Ecosystem Exchange measured 620 with eddy covariance technique: algorithms and uncertainty estimation. Biogeosciences 3, 621 571-583 (2006). 622 54. Myneni, R. B., Ramakrishna, R., Nemani, R., Running, S. W. Estimation of global leaf area 623 index and absorbed PAR using radiative transfer models. IEEE Transactions on Geoscience 624 and remote sensing 35, 1380-1393 (1997). 625 55. Maisongrande, P., Duchemin, B., Dedieu, G. VEGETATION/SPOT: an operational 626 mission for the Earth monitoring; presentation of new standard products. International 627 Journal of Remote Sensing 25, 9-14 (2004). 628 56. Vermote, E. F., El Saleous, N. Z., Justice, C. O. Atmospheric correction of MODIS data in 629 the visible to middle infrared: first results. Remote Sensing of Environment 83, 97-111 630 (2002). 631 57. Rahman, A. F., Sims, D. A., Cordova, V. D., El‐Masri, B. Z. Potential of MODIS EVI 632 and surface temperature for directly estimating per-pixel ecosystem C fluxes. Geophysical 633 Research Letters 32, 156-171 (2005). 634 58. Frankenberg, C. et al. New global observations of the terrestrial carbon cycle from GOSAT: 635 Patterns of plant fluorescence with gross primary productivity. Geophysical Research 636 Letters 38, 351-365 (2011). 637 59. Lee, J. E. et al. Forest productivity and water stress in Amazonia: observations from 638 GOSAT chlorophyll fluorescence. Proceedings of the Royal Society B Biological Sciences 639 280, 176-188 (2013). 640 60. Kottek, M., Grieser, J., Beck, C., Rudolf, B., Rubel, F. World map of the Köppen-Geiger 641 climate classification updated. Meteorologische Zeitschrift 15, 259-263 (2006). 642 61. Qin, J., Yang, K., Liang, S. & Guo, X. The altitudinal dependence of recent rapid warming 643 over the Tibetan Plateau. Climatic Change 97, 321 (2009). 644 62. Liu, D. et al. Contrasting responses of grassland water and carbon exchanges to climate 645 change between Tibetan Plateau and Inner Mongolia. Agricultural and Forest Meteorology 646 249, 163-175 (2018). 647 63. Wan, Z. & Dozier, J. A generalized split-window algorithm for retrieving land-surface 648 temperature from space. IEEE Transactions on geoscience and remote sensing 34, 892-905 649 (1996). 650 64. Zhu, Z. et al. Greening of the Earth and its drivers. Nature Climate Change, 6, 701-795 651 (2016). 652 65. Weedon, G. P. et al. The WFDEI meteorological forcing data set: WATCH Forcing Data 653 methodology applied to ERA-Interim reanalysis data. Water Resources Research 50, 7505- 654 7514 (2014). 655 66. Jones, P. W., First-and second-order conservative remapping schemes for grids in spherical 656 coordinates. Monthly Weather Review 127, 2204-2210 (1999). 657 67. Battaglia, M., Beadle, C., Loughhead, S. Photosynthetic temperature responses of 658 Eucalyptus globulus and Eucalyptus nitens. Tree Physiology 16, 81-89 (1996). 659 68. Lin, Y. S., Medlyn, B. E. & Ellsworth, D. S. Temperature responses of leaf net 660 photosynthesis: the role of component processes. Tree Physiology 32, 219-231 (2012). 661 69. Farquhar, G. v., von Caemmerer, S. v. & Berry, J. A biochemical model of photosynthetic 662 CO2 assimilation in leaves of C3 species. Planta 149, 78-90 (1980). 663 70. Baker, N. R. Photosynthesis and the Environment. (Kluwer Academic Publishers, 2006). 664 71. Helliker, B. R. & Richter, S. L. Subtropical to boreal convergence of tree-leaf temperatures. 665 Nature 454, 511 (2008). 666 72. Michaletz, S. T. et al. The energetic and carbon economic origins of leaf thermoregulation. 667 Nature Plants 2, 16129 (2016). 668 73. Lian, X. et al. Spatiotemporal variations in the difference between satellite-observed daily 669 maximum land surface temperature and station‐based daily maximum near-surface air 670 temperature. Journal of Geophysical Research: Atmospheres 122, 2254-2268 (2017). 671 672 Acknowledgments: This study was supported by the Strategic Priority Research Program (A) of 673 the Chinese Academy of Sciences (grant XDA20050101), the National Natural Science Foundation 674 of China (41530528) and the National Key R&D Program of China (2017YFA0604702). This work 675 used eddy covariance data acquired by the FLUXNET community and in particular by the 676 following networks: AmeriFlux (U.S. Department of Energy, Biological and Environmental 677 Research, Terrestrial Carbon Program (DE-FG02-04ER63917 and DE-FG02-04ER63911)), 678 AfriFlux, AsiaFlux, CarboAfrica, CarboEuropeIP, CarboItaly, CarboMont, ChinaFlux, Fluxnet- 679 Canada (supported by CFCAS, NSERC, BIOCAP, Environment Canada, and NRCan), 680 GreenGrass, KoFlux, LBA, NECC, OzFlux, TCOS-Siberia, USCCC. We acknowledge the 681 financial support to the eddy covariance data harmonization provided by CarboEuropeIP, FAO- 682 GTOS-CO, iLEAPS, Max Planck Institute for Biogeochemistry, National Science Foundation, 683 University of Tuscia, Université Laval and Environment Canada and US Department of Energy 684 and the database development and technical support from Berkeley Water Center, Lawrence 685 Berkeley National Laboratory, Microsoft Research eScience, Oak Ridge National Laboratory, 686 University of California-Berkeley, University of Virginia. P.C., J.P., and I.A.J. would like to 687 acknowledge the financial support from the European Research Council Synergy grant ERC-SyG- 688 2013-610028 IMBALANCE-P. P.C. was also supported from the French Agence Nationale de la 689 Recherche (ANR) Con-vergence Lab Changement climatique et usage des terres (CLAND). I.A.J. 690 acknowledges the Methusalem funding of the Flemish Community through the Research Council 691 of the University of Antwerp. T.F.K. was supported by the NASA Terrestrial Ecology Program 692 IDS Award NNH17AE86I. J. Mao and X. Shi are supported by the Terrestrial Ecosystem Science 693 Scientific Focus Area project funded through the Terrestrial Ecosystem Science Program in the 694 Climate and Environmental Sciences Division of the Biological and Environmental Research 695 Program in the US Department of Energy Office of Science. Oak Ridge National Laboratory is 696 supported by the Office of Science of the US Department of Energy under Contract No. DE-AC05- 697 00OR22725. M. C. was supported by a grant overseen by the French National Research Agency 698 (ANR) as part of the "Investissements d’Avenir" program (ANR-11-LABX-0002-01, Lab of 699 Excellence ARBRE). 700 701 Author contributions: S.L.P. designed the research; M.T.H. performed the analysis; S.L.P. 702 drafted the paper; and all authors contributed to the interpretation of the results and to the text. 703 704 Competing interests: Authors declare no competing interests. 705 706 Figure legends 707 Fig. 1 | Distribution of ecosystem-scale optimal temperature (𝑻𝒆𝒄𝒐 𝒐𝒑𝒕 ) for vegetation productivity 708 derived from flux tower sites and satellite-based data for near-infrared reflectance of 709 vegetation (NIRV). a, Relationship between mean annual daily maximum air temperature during 𝑎𝑖𝑟 𝑒𝑐𝑜 710 the growing season (𝑇𝑚𝑎𝑥 𝑔𝑠 ) and 𝑇𝑜𝑝𝑡 derived from daily measurements of photosynthesis across 𝑎𝑖𝑟 𝑒𝑐𝑜 711 eddy-covariance sites. Flux-derived 𝑇𝑚𝑎𝑥 𝑔𝑠 and 𝑇𝑜𝑝𝑡 were both obtained using observations from 712 flux towers. Error bars indicate ±SD. The dotted gray line represents y=x and the dot line in red is 713 y=0.61x+10.65, which is derived by linear regression with the statistical significance of the slope, 𝑒𝑐𝑜 714 or its p-value, given by Student’s t test. b, Relationship between 𝑇𝑜𝑝𝑡 derived from flux data and 𝑒𝑐𝑜 𝑒𝑐𝑜 715 𝑇𝑜𝑝𝑡 derived from NIRV data. For each site, we extracted and averaged 𝑇𝑜𝑝𝑡 values within a 3×3 𝑒𝑐𝑜 𝑒𝑐𝑜 716 pixel window around the site from NIRV-derived 𝑇𝑜𝑝𝑡 map, and calculated the SD of the nine 𝑇𝑜𝑝𝑡 717 values within the window. Error bars indicate ±SD. The dotted gray line represents y=x and the dot 718 line in red is y=0.74x+7.10, which is derived by linear regression with the statistical significance 𝑒𝑐𝑜 719 of the slope, or its p-value, given by Student’s t test. c, Spatial distribution of 𝑇𝑜𝑝𝑡 for vegetation 𝑒𝑐𝑜 𝑒𝑐𝑜 720 productivity (left panel), and 𝑇𝑜𝑝𝑡 averaged by latitude (right panel). 𝑇𝑜𝑝𝑡 is determined using 721 NIRV data calculated based on satellite observations from Moderate Resolution Imaging 722 Spectroradiometer (MODIS). Note that only gridded pixels with annual mean NDVI value larger 𝑒𝑐𝑜 723 than 0.1 and detectable 𝑇𝑜𝑝𝑡 are shown here. Areas of tropical forests based on current vegetation 724 distribution are indicated by hatching. The circles on the map are colored according to the local 𝑒𝑐𝑜 725 value of 𝑇𝑜𝑝𝑡 retrieved from GPP at the location of each flux site. The solid line and shaded area 𝑒𝑐𝑜 𝑒𝑐𝑜 726 in the right panel indicate the mean and SD, respectively, of 𝑇𝑜𝑝𝑡 summarized by latitude. d, 𝑇𝑜𝑝𝑡 𝑒𝑐𝑜 727 in the climate space (left panel) and the temperature sensitivity of 𝑇𝑜𝑝𝑡 along the precipitation 𝑎𝑖𝑟 728 gradient (right panel). Each climate bin is defined by 1-oC intervals of 𝑇𝑚𝑎𝑥 𝑔𝑠 and 100-mm 729 intervals of mean annual precipitation, based on current climate conditions averaged over 2001- 𝑒𝑐𝑜 730 2013. The solid line in the right panel represents the temperature sensitivity of 𝑇𝑜𝑝𝑡 along the 𝑒𝑐𝑜 𝑎𝑖𝑟 731 precipitation gradient, calculated as the slope of the linear regression between 𝑇𝑜𝑝𝑡 and 𝑇𝑚𝑎𝑥 𝑔𝑠 for 𝑒𝑐𝑜 732 a given precipitation level. The shaded area indicates the SD of temperature sensitivity of 𝑇𝑜𝑝𝑡 733 estimated by bootstrapping. 734 735 Fig. 2 | Relationship between mean annual daily maximum air temperature during the 736 growing season ( 𝑻𝒂𝒊𝒓 𝒎𝒂𝒙 𝒈𝒔 ) and ecosystem-scale optimum temperature for vegetation 737 productivity (𝑻𝒆𝒄𝒐 𝑒𝑐𝑜 𝑎𝑖𝑟 𝒐𝒑𝒕 ) across vegetation types. The error bars indicate the SDs of 𝑇𝑜𝑝𝑡 /𝑇𝑚𝑎𝑥 𝑔𝑠 for 738 each vegetation type: ENF, evergreen needle-leaved forest; EBF, evergreen broad-leaved forest; 739 DNF, deciduous needle-leaved forest; DBF, deciduous broad-leaved forest; MF, mixed forest; 740 Shrub, closed and open shrublands. The light-gray dotted line represents y=x. The dark-gray dotted 741 line is y=0.76x+6.48 derived by linear regression with the slope value (estimated using Student’s t 742 test) shown in the bottom right. The red dotted line is the flux tower derived slope (0.61) from Fig. 743 1a. The size of each symbol corresponds to the three categories (< 3%, 3%-10% and > 10%) of 744 occupied vegetated area on land. Error bars indicate ±SD. 745 746 Fig. 3 | Change with latitude in ecosystem-scale optimal temperature for vegetation productivity (𝑻𝒆𝒄𝒐 𝒐𝒑𝒕 ) and daily 𝑒𝑐𝑜 747 maximum air temperature averaged over the growing season (𝑻𝒂𝒊𝒓 𝑎𝑖𝑟 𝒎𝒂𝒙 𝒈𝒔 ). a, Current 𝑇𝑜𝑝𝑡 versus current 𝑇𝑚𝑎𝑥 𝑔𝑠 ; b, 𝑒𝑐𝑜 𝑎𝑖𝑟 𝑒𝑐𝑜 𝑎𝑖𝑟 748 Future 𝑇𝑜𝑝𝑡 versus future 𝑇𝑚𝑎𝑥 𝑔𝑠 . Current 𝑇𝑜𝑝𝑡 and 𝑇𝑚𝑎𝑥 𝑔𝑠 are calculated using current temperature for 2001-2013, 𝑒𝑐𝑜 𝑎𝑖𝑟 749 whereas acclimated 𝑇𝑜𝑝𝑡 and future 𝑇𝑚𝑎𝑥 𝑔𝑠 are first calculated pixel by pixel using temperature for 2091-2100 750 projected by General Circulation Models (GCMs) under the RCP4.5 scenario and then averaged by latitude. 𝑒𝑐𝑜 𝑒𝑐𝑜 751 Acclimated 𝑇𝑜𝑝𝑡 is determined based on the projected temperature and temperature sensitivity of 𝑇𝑜𝑝𝑡 using the annual 752 precipitation level predicted for 2091-2100. The solid line and shaded area in each panel indicate the mean and SD, 𝑒𝑐𝑜 𝑎𝑖𝑟 𝑒𝑐𝑜 𝑎𝑖𝑟 753 respectively, of 𝑇𝑜𝑝𝑡 or 𝑇𝑚𝑎𝑥 𝑔𝑠 summarized by latitude. c, Future 𝑇𝑜𝑝𝑡 versus future 𝑇𝑚𝑎𝑥 𝑔𝑠 for tropical evergreen 𝑒𝑐𝑜 𝑎𝑖𝑟 754 forests. ** indicates that 𝑇𝑜𝑝𝑡 is significantly lower than 𝑇𝑚𝑎𝑥 𝑔𝑠 at P<0.01 in a paired t-test. Error bars indicate ±SD.