This is the accepted version of the article: Wang, Songhan; Ju, Weimin; Peñuelas, Josep; [et al.]. «Urban−rural gradients reveal joint control of elevated CO2 and temperature on extended photosynthetic seasons». Nature ecology and evolution, Vol. 3 (March 2019), p. 1076–1085. DOI 10.1038/s41559-019-0931-1 This version is avaible at https://ddd.uab.cat/record/218190 under the terms of the license 1 Urban−rural gradients reveal joint control of elevated CO2 and temperature on 2 extended photosynthetic seasons 3 Songhan Wang1, 2, 3, Weimin Ju1, 2, Josep Peñuelas4, 5, Alessandro Cescatti6, Yuyu Zhou7, 4 Yongshuo Fu8, Alfredo Huete9, Min Liu10, 11, Yongguang Zhang1, 2, 3* 1 5 International Institute for Earth System Sciences, Nanjing University, Nanjing, China. 2 6 Jiangsu Center for Collaborative Innovation in Geographical Information Resource Development and 7 Application, Nanjing, China. 3 8 Collaborative Innovation Center of Novel Software Technology and Industrialization, Nanjing, 210023, 9 China. 4 10 CSIC, Global ecology Unit CREAF-CSIC-UAB, Bellaterra 08193, Catalonia, Spain. 5 11 CREAF, Bellaterra 08193, Catalonia, Spain. 6 12 European Commission, Joint Research Centre, Directorate for Sustainable Resources, Ispra, Italy. 7 13 Department of Geological and Atmospheric Sciences, Iowa State University, Ames, IA 50011, USA. 8 14 College of water sciences, Beijing Normal University, Beijing, China. 9 15 School of Life Sciences, University of Technology Sydney, NSW 2007, Australia. 10 16 Institute of Eco-Chongming, Shanghai 200062, China. 11 17 Shanghai Key Laboratory for Urban Ecological Processes and Eco-Restoration, School of Ecological 18 and Environmental Sciences, East China Normal University, Shanghai 200241, China. 19 *Corresponding author: 20 Yongguang Zhang 21 International Institute for Earth System Sciences 22 Nanjing University 23 163 Xianlin Avenue, Qixia District, Nanjing, Jiangsu Province, China 24 Email: yongguang_zhang@nju.edu.cn 25 Tel. +86-(0)25-89681569 1 26 Abstract 27 Photosynthetic phenology has large effects on the land-atmosphere carbon exchange. 28 Due to limited experimental assessments, a comprehensive understanding of the 29 variations of photosynthetic phenology under future climate and its associated controlling 30 factors is still missing, despite its high sensitivities to climate. Here we develop an 31 approach that uses cities as natural laboratories, since plants in urban areas are often 32 exposed to higher temperature and carbon dioxide (CO2) concentration, which reflect 33 expected future environmental conditions. Using more than 880 urban–rural gradients 34 across the Northern Hemisphere (≥30°), combined with concurrent satellite retrievals of 35 sun-induced chlorophyll fluorescence (SIF) and atmospheric CO2, we investigated the 36 combined impacts of elevated CO2 and temperature on photosynthetic phenology at large 37 scale. Results show that in urban conditions of elevated CO2 and temperature, vegetation 38 photosynthetic activity began earlier (−5.6 ± 0.7 days), peaked earlier (−4.9 ± 0.9 days) 39 and ended later (4.6 ± 0.8 days) than in neighboring rural areas, with a striking two- to 40 four-fold higher climate sensitivity than greenness phenology. The earlier start and peak 41 of season were sensitive to both the enhancements of CO2 and temperature, whereas the 42 delayed end of season was mainly attributed to CO2 enrichments. We used these 43 sensitivities to project phenology shifts under four “Representative Concentration 44 Pathway” (RCP) climate scenarios, predicting that vegetation will have prolonged 45 photosynthetic seasons in the coming two decades. This observation-driven study 46 indicates that realistic urban environments, together with SIF observations, provide a 2 47 promising way for vegetation physiological studies under future climate change. 48 3 49 Terrestrial ecosystems absorb approximately one quarter of the anthropogenic carbon 50 dioxide (CO2) released into the atmosphere, which effectively mitigates climate warming1. 51 The timing and duration of vegetation photosynthetic activity (i.e., photosynthetic 52 phenology), play an important role in the carbon cycle of land ecosystems2,3 and are 53 strongly affected by climate change4-7. Understanding the response of photosynthetic 54 phenology to increasing atmospheric CO2 concentrations and global warming is, 55 therefore, essential to better understand the future dynamics of the terrestrial carbon 56 cycle8,9 and formulate land-based strategies to mitigate climate change10. Current 57 experimental assessments of the effects of climate warming and CO2 fertilisation on plant 58 phenology have generally been restricted to manipulated experiments11,12 with very 59 limited spatial extent and species coverage, which may ultimately led to ambiguous and 60 controversial findings13. In particular, the role of CO2 fertilisation on tree phenology 61 remains largely unknown in real natural settings. Moreover, many field experiments have 62 been recently terminated because of their high cost14, thus limiting consistent long-term 63 observations and reinforcing the need for effective, low-cost alternative monitoring 64 approaches. 65 In this analysis, we make use of contrasting conditions in urban–rural environments, 66 serving as “natural laboratories” to conduct a novel, and first time investigation on the 67 combined effects of future global warming and CO2 fertilisation on vegetation phenology 68 at a very large scale14. Because of the urban heat island effect15 and the larger CO2 69 emissions from fossil-fuel combustion16, vegetation in urban areas is typically exposed to 4 70 higher temperatures and CO2 concentrations than in rural areas and are therefore 71 experiencing growing conditions projected for the future17. Assuming that phenology 72 responds similarly to altered temperature and atmospheric CO2 in urban and rural areas, 73 the urban–rural gradients can be considered as “natural experiments” with similar 74 photoperiod and weather, and at lower costs14 than manipulation experiments. Using the 75 space-for-time substitution concept, we therefore used the physiological dynamics of 76 vegetation across urban–rural gradients to investigate the effects of future global warming 77 and CO2 fertilisation on photosynthetic phenology. Data were collected from 880 urban 78 clusters and their adjacent rural areas at mid- to high-latitudes in the Northern 79 Hemisphere (NH, ≥30°) (Supplementary Fig. 1). Photosynthetic phenology was extracted 80 from the Orbiting Carbon Observatory-218 (OCO-2) sun-induced chlorophyll 81 fluorescence (SIF, a proxy for photosynthesis) data, which are highly correlated with 82 photosynthetic rates and can monitor actual photosynthetic dynamics19. We used the 83 traditional, satellite greenness phenology for comparison, which was extracted from 84 Moderate Resolution Imaging Spectroradiometer (MODIS) Enhanced Vegetation Index 85 (EVI)20. The urban–rural surface CO2 gradients were obtained through a conversion 86 approach with OCO-2 column-averaged CO2 mixing ratios (XCO2) data (see Methods). 87 The primary objective of our analysis was to determine the differences in photosynthetic 88 phenology across urban–rural gradients and explore the combined effect of atmospheric 89 CO2 enrichment and warming on plant photosynthetic phenology at a very large scale. 90 Phenological differences and environmental gradients were then used to quantify the 5 91 sensitivities of the processes and ultimately to project shifts in phenophases under 92 contrasting climate scenarios. 93 Results 94 We found that vegetation showed a larger than expected anticipation and 95 prolongation of the photosynthetic season within urban environments. On average, the 96 SIF-based photosynthetic activity began 5.6 days earlier, peaked 4.9 days earlier, ended 97 4.6 days later and lasted 10.2 days longer in urban than in rural areas (Fig. 1). These 98 variations were two- to four-times larger than those in greenness phenology derived from 99 MODIS EVI, due to the significantly higher urban–rural ∆SIF than ∆EVI during spring 100 and autumn (Supplementary Fig. 2). The timing of each phenological indicator we 101 examined (SOS, start of growing season; POS, peak of growing season; EOS, end of 102 growing season; GSL, length of growing season) was generally highly correlated between 103 urban areas and neighbouring rural counterparts, but with significant systematic 104 deviations, as indicated by the slopes of the regression lines (Supplementary Fig. 3). The 105 SOS and POS occurred earlier in urban areas as suggested by the slopes of larger than 106 one (i.e., 1.01–1.03 and 1.00–1.02; p < 0.001) for the regression lines, whereas the EOS 107 occurred later in urban areas as indicated by the slopes of lower than one (i.e., 0.98–0.99; 108 p < 0.001). On average within cities, SOS and POS were therefore advanced, while EOS 109 was delayed and GSL was extended for the phenology of SIF. Similar patterns were 110 found for the phenology of EVI, but the magnitudes of these changes were substantially 111 smaller than that of SIF (Supplementary Fig. 4) (about four-fold for ∆SOS and ∆POS, 6 112 two-fold for ∆EOS and three-fold for ∆GSL; Fig. 1). The larger differences of 113 photosynthetic phenology in urban–rural gradients than the differences of greenness 114 phenology were not affected by the mismatching in the spatial resolutions of SIF and EVI 115 data, as also tested after resampling to the same resolution (Supplementary Fig. 5). These 116 results suggest a remarkable shift in the response of photosynthetic phenology to urban 117 environmental conditions in the mid- to high-latitudes of the NH. These variations cannot 118 be completely explained by the variations in greenness phenology based on VIs or leaf 119 area index (LAI) and probably involve physiological responses to environmental drivers. 120 Vegetation in urban areas generally had earlier photosynthetic SOS and POS, later 121 EOS and longer GSL in the NH across various latitudes and climatic backgrounds (Fig. 122 2). On average, more than 70% of the urban clusters had longer photosynthetic seasons 123 than rural buffers. ∆SOS and ∆POS were lower than zero days across the whole NH, 124 while most of the ∆EOS were greater than zero days, suggesting that the trends of earlier 125 start/peak and delayed end of photosynthesis in urban environments are common across 126 northern latitudes. Nevertheless, some spatial patterns of phenological shifts along 127 urban–rural gradients were also evident; ∆SOS tended to decrease slightly as latitude 128 increased, ∆EOS and ∆GSL had the highest values at mid-latitudes, whereas ∆POS had 129 no clear trend (second column in Fig. 2). The spatial patterns of SIF ∆SOS were not fully 130 consistent with previous studies conducted with EVI, which showed an increasing trend 131 of EVI ∆SOS with the latitude21 due to the larger urban heat island effects in 132 high-latitudes15,22. This divergence possibly suggests that the spatial distributions of SIF 7 133 urban–rural phenological gradients cannot be explained by temperature alone. 134 Importantly, urban–rural phenological differences were more significant in warm and dry 135 regions (i.e. areas with high mean annual air temperature and low annual precipitation) 136 than in other regions (third column in Fig. 2). Thus, the phenological differences in arid 137 climate zones were significantly higher than that in temperate and boreal areas 138 (Supplementary Fig. 6). These results suggest a general phenomenon that urban 139 environmental conditions extend the photosynthetic seasons and that the impacts in warm 140 and dry regions are greater than those in other climate regions. 141 Urban–rural phenological difference is likely driven by several interacting factors 142 including the enhancements of temperature and CO2, the geolocation of a city, the urban 143 size and climatic background. Partial correlation analysis showed that SIF ∆SOS and 144 ∆POS were mainly controlled by the daytime urban heat island (∆LSTday) effects and 145 urban–rural CO2 gradients (∆CO2), of which ∆LSTday was the most significant (Fig. 3). 146 SIF ∆EOS was mainly attributed to the enhanced CO2 only, resulting in similarly 147 significant effects of temperature and CO2 on ∆GSL (Fig. 3). In order to attribute the 148 effects of environmental drivers on photosynthetic phenology, we further analysed the 149 associations between elevated daytime temperature or atmospheric CO2 and urban–rural 150 phenological differences (Fig. 4). The onset and peak time of photosynthesis were both 151 advanced significantly as atmospheric CO2 concentration increased (p = 0.013 and 0.021, 152 respectively) and as temperature increased (p = 0.019 and 0.006, respectively). The end 153 of photosynthetic activity was significantly delayed as CO2 concentrations increased (p = 8 154 0.028), whereas its association with temperature was not significant (p = 0.259). Our 155 analyses suggest that not only daytime temperature but also atmospheric CO2 have large 156 impacts on photosynthetic phenology of urban plants and, therefore, may potentially 157 inform on the future phenological shifts of natural vegetation under scenarios of warming 158 and atmospheric CO2 enrichment. 159 Finally, using observations along the urban–rural gradients as natural laboratories, 160 we projected the shifts in photosynthetic phenology for the next two decades using the 161 space-for-time substitution. Based on forward stepwise regression models (see Methods), 162 the sensitivities of ∆SOS, ∆POS and ∆GSL to daytime temperature were –1.55 ± 0.60, 163 –2.11 ± 0.78 and 3.31 ± 1.10 days·°C–1, respectively (Fig. 5a). The sensitivities of ∆SOS, 164 ∆POS, ∆EOS and ∆GSL to surface CO2 concentration were –0.15 ± 0.07, –0.18 ± 0.09, 165 0.17 ± 0.07 and 0.32 ± 0.11 days·ppm–1, respectively (Fig. 5b). Given the observed 166 sensitivities and four different “Representative Concentration Pathway” (RCP) scenarios 167 of surface temperature and CO2 concentration10 (i.e., RCP 2.6, 4.5, 6.0 and 8.5, 168 Supplementary Fig. 7), we projected the possible shifts of photosynthetic SOS, POS, 169 EOS and GSL from 2015 to 2035 (Fig. 5c-f). Results show that the phenophase shifts 170 were lowest in the RCP 2.6 scenario and highest in the RCP 8.5 scenario. Assuming 171 constant sensitivities for the next two decades, under the range of climate pathways 172 foreseen by RCP scenarios during 2015–2035, the photosynthetic SOS and POS would 173 advance at a rate of 2.8–5.7 and 3.3–6.9 days·decade–1, respectively; and the EOS would 174 delay at a rate of 3.1–5.7 days·decade–1, prolonging the growing season by a rate of 9 175 5.9–11.4 days·decade–1 (Fig. 5g). It should be noted that the responses of vegetation 176 phenology to temperature are not stationary23 and are likely to decline in the near 177 decades9,24. Besides, the sensitivity of vegetation phenology to CO2 may not be static 178 either in the coming decades, since the projected future concentrations of CO2 could 179 exceed the range currently experienced by the urban vegetation (Fig. 4 & Supplementary 180 Fig. 7b). Therefore, our analysis of the projected photosynthetic phenology based on the 181 static sensitivity may overestimate the future shifts of vegetation phenophases. This 182 caveat, however, can be solved in the future by the same procedure introduced in this 183 study and continued satellite SIF observations. 184 Discussion 185 We found systematic differences between urban–rural gradients of photosynthetic 186 and greenness phenology (Fig. 1). The urbanization effects on vegetation photosynthetic 187 phenology were more than twofold higher than those on greenness. This phenomenon 188 could be explained by the significant differences between season cycle of EVI and SIF, 189 namely EVI shows an earlier spring increase and a later autumn falling than SIF 190 (Supplementary Fig. 2a). This difference has also been demonstrated by a recent study25. 191 The timing of carbon assimilation in spring lags behind the leaf burst, which varies with 192 foliar structure and longevity26. The photosynthetic activity shuts down before the 193 reduction of leaf chlorophyll and leaf abscission because of the limitation of light 194 availability27-29. Another possible reason is that urban phenology indicators from VIs are 10 195 suffering from mixed-pixel effects30, which can bias the gradients of urban–rural 196 greenness phenology. Besides, VIs may have noisy signals from soil background, such as 197 soil color, artificial green building and other non-photosynthetic active materials, 198 especially in urban areas31-33. On the other hand, originating from the vegetation 199 photosynthetic pigments, SIF is less sensitive to the soil background34-36. 200 Our analyses suggest that the urban–rural photosynthetic phenological differences 201 are mainly controlled by elevated atmospheric CO2 concentration and daytime 202 temperature, while the differences of greenness phenology are mostly depending on the 203 urban size (Fig. 3). Previous studies using dozens of cities at regional scale suggest a 204 positive correlation between earlier greenup and surface temperature enhancement of 205 urban areas, such as in eastern North America37,38 and China21,39,40. However, our analysis 206 show that the urban–rural greenness phenological differences are mainly correlated to the 207 urban size in the Northern Hemisphere, which is in line with a recent study conducted 208 with thousands of cities in the conterminous United States41. The positive relationship 209 between urban–rural greenness phenological gradients and urban size (an indicator of 210 urbanization), may result from the shifts of vegetation percentages and regional climate 211 conditions induced by urban expansion41. On the contrary, the earlier start and peak of 212 photosynthesis activity inferred from SIF observations have a significant correlation with 213 elevated temperature in urban areas (Fig. 4a&b). This discrepancy might partly relate to 214 the different sensitivity of VIs and SIF to temperature increases, either from satellite 215 data25 or previous site-level study42. More importantly, SIF has a positive correlation with 11 216 GPP both in leaf-level and site-level19,43,44, which is a more direct indicator of vegetation 217 physiological status than VIs45. Seasonality of SIF from ground-based measurements46, 218 tower-based observations47 and satellite-based data25 generally has high correlations with 219 that of GPP. Therefore SIF-based photosynthetic phenology in studies of “urban 220 laboratories” could provide a better proxy of the actual vegetation photosynthetic changes 221 in response to the climate change than VIs, which provide proxies of vegetation 222 greenness. 223 Using 880 cities, our results confirm that vegetation photosynthesis activities in 224 spring are triggered by daytime temperature, not by nighttime temperature8 (Fig. 3). 225 However, the sensitivity of SOS to daytime temperature in this study (1.55±0.60 days 226 per °C) is lower than previous studies for boreal and temperate forests based on longtime 227 satellite VIs or ground observations8,9,48-52. This may reflect the declining warming effects 228 on spring leaf unfolding9 and weakening temperature control on the interannual 229 variations of vegetation spring carbon uptake53. The reduced sensitivity of spring leaf 230 unfolding to temperature may partly result from the reduced chilling accumulation in 231 urban areas9, in which the winter temperatures may become insufficiently low to meet the 232 requirement of chilling due to the urban heat island effects. Photoperiod may also 233 influence the warming effects on spring phenology, although its impact remains unclear 234 and is under debate5,23,54,55. We also found that the photosynthetic POS have an earlier 235 trend (about 3.3–6.9 days per decade) in the near decades, which is consistent with the 236 POS trend of δ13C data (about 4.3 ± 2.9 days per decade)56, and are larger than the POS 12 237 trends from NDVI data or ground phenology measurements48,56,57. This phenomenon 238 highlights the advantages of SIF for tracking vegetation seasonal photosynthetic activities 239 and carbon uptake. Our analysis also reveals a delayed end of vegetation photosynthesis 240 in urban areas (Fig. 1). Other than the earlier start of SOS and its response to global 241 warming, our understanding of the drivers for delayed EOS is still limited. 242 Temperature12,58, photoperiod59, elevated atmospheric CO260 and spring phenology61 were 243 claimed to be associated to this delayed vegetation senescence. Nevertheless, our analysis 244 suggests a pronounced delay in the end of vegetation photosynthesis in urban 245 environments and these delays are highly correlated to the elevated atmospheric CO2 246 concentration. 247 Current studies generally focus on the responses of phenophases to warming9,13, 248 ignoring the importance of additional environmental drivers. From this point of view our 249 analysis highlights that photosynthetic phenological shifts are not only controlled by 250 daytime temperature but also by atmospheric CO2. CO2 fertilisation seems to have a large 251 effect on the advancement of photosynthetic activity, particularly under warm and dry 252 conditions; and the delayed end of photosynthesis is mainly attributed to the elevated 253 atmospheric CO2 rather than temperature, which is consistent with field experiments62-64. 254 The delayed EOS may be because photosynthesis could be sustained for longer period in 255 the absence of sink limitations under elevated atmospheric CO2 concentrations64. Direct 256 evaluations of the future projected phenophase shifts are difficult, however, given a 257 global increase in CO2 concentration of 16.5 ppm·decade−1 and temperature of 13 258 0.16 °C·decade−1 in the past three decades65 (1982–2011), our estimates of the prolonged 259 photosynthetic growing seasons (5.5 ± 2.1 days·decade−1) are significantly higher than 260 current estimates (2.2 to 3.9 days·decade−1) based on long-term satellite vegetation 261 indices66,67 and the in situ experimental results based on greenness observations from a 262 recent study5. Our results suggest that the sensitivity of photosynthetic phenology to 263 climate is strikingly higher than that of traditional greenness phenology, and therefore, 264 that the current methods based on satellite VIs and in situ greenness observations are 265 likely to underestimate the extension of carbon sink period under climate change. 266 Crucially, our results show that, with increasing atmospheric CO2, the growing seasons 267 are likely to expand in the next two decades, which in turn might contribute to an 268 increasing trend of the terrestrial carbon sink and generate an important negative 269 feedback in the climate system. 270 Yet it should be noted that some limitations still remain in our analysis. Landscape 271 configuration, species composition68, atmospheric deposition (e.g., nitrogen and 272 phosphorus), air pollutants (e.g., ozone), management practices and hydrological regimes 273 may also affect urban–rural phenological differences. These factors may have complex 274 and to some extent counterbalanced impacts on vegetation photosynthesis in urban 275 environments. For example, urban areas generally receive higher rates of atmospheric 276 nitrogen and phosphorus deposition69, which may positively contribute to urban–rural 277 phenological differences (Supplementary Fig. 8b). On the other hand, urban areas 278 typically show higher daytime ozone concentrations that can be detrimental to 14 279 photosynthesis70, and these urban–rural ozone differences may negatively contribute to 280 the urban–rural phenological differences (Supplementary Fig. 8c). Besides, vegetation 281 species and human management practices may vary between urban and rural areas; and 282 apart from the temperature and CO2, water stress indices such as precipitation, vapour 283 pressure deficit (VPD) and soil moisture, may also have influences on the vegetation 284 phenology71-73. Urban–rural phenological gradients in areas with low precipitation are 285 generally more pronounced (Supplementary Fig. 8a), because vegetation in urban areas is 286 likely to be affected by irrigation and thus not limited by water supply. To test whether 287 our results were affected by the irrigation of urban vegetation, we repeated the analysis 288 after excluding the cities in the arid climate zones, assuming that the irrigation effect on 289 urban–rural phenology gradients could be more significant in these areas. Results also 290 showed that the extended photosynthetic seasons in urban areas were both controlled by 291 the elevated temperature and CO2 (Supplementary Fig. 9a), but with slightly lower 292 sensitivities. Therefore, our predictions of the GSL shifts based on all of the 880 cities 293 may be slightly overestimated with a rate of about 0.8–1.2 days·decade–1 (Supplementary 294 Fig. 9b). Moreover, to verify the robustness of our analysis, we further calculated the 295 phenological differences between the 10-km and 30-km rural buffers and related them 296 with temperature and CO2 gradients. Results from the gradients between 10-km and 297 30-km rural buffers showed higher correlations with temperature and CO2 gradients than 298 those from the original urban–rural gradients (Supplementary Fig. 10). These higher 299 correlations not only confirm the robustness of our analysis, but also suggest that the 15 300 confounding factors (e.g., vegetation species, irrigation and nutrient supply) do not affect 301 the controlling role of temperature and CO2 in vegetation photosynthetic phenology, 302 given that the vegetation phenology in 10-km rural buffers are less influenced by these 303 factors than urban areas. However, these confounding factors can hardly be included in 304 our analysis due to the lack of high-resolution data at global scale. 305 The novel insights of this study are derived at large scale from the contrasting 306 behaviour of plants along the urban–rural gradients. The sharp differences in 307 environmental conditions along these gradients can represent a unique “natural laboratory” 308 with different atmospheric CO2 concentration and temperature that may ultimately inform 309 on plant behaviours under future climates and atmospheric composition. The 310 macro-environmental conditions of open urban laboratories are widely reachable in the 311 globe and can therefore lead to more general and robust assessments of vegetation 312 responses to key environmental drivers across different environmental conditions in a 313 global change perspective. Ultimately, our findings shed new light on the behaviour of 314 plants under global change, help to reduce the uncertainties of terrestrial ecosystems 315 models74, and increase our understanding of climate–vegetation interactions. 316 16 317 Methods 318 Urban–rural gradients 319 We focused our study on urban clusters and their adjacent rural buffers at mid- to 320 high-latitudes in the Northern Hemisphere (NH, ≥30°) from 2015 to 2017 321 (Supplementary Fig. 1), where vegetation is highly seasonal and sensitive to climate61,66. 322 The extents of the urban areas were mapped using a cluster-based method that estimates 323 optimal thresholds using nighttime stable-light data from the Defense Meteorological 324 Satellite Program/Operational Linescan System (DMSP/OLS)75. Sensitivity analyses and 325 comparisons with other global urban-area products indicated that this product is reliable 326 and highly accurate76. We used this data product (for 2013) instead of extracting urban 327 areas from land-cover data sets, because these cluster-based urban extents clearly define 328 the boundaries between urban and rural areas, regardless of the internal spatial 329 heterogeneity, which is more similar to the real environment that the vegetation 330 experiences. Although the static urban clusters (2013) and our study period (2015–2017) 331 differed by approximately three years, the overall results should not be affected due to 332 almost stable urban clusters in North America and Europe41 and the slightly increasing 333 urban areas in Asia77. 334 We created a series of buffers extending 10, 20, 30 and 40 km outward from each 335 urban cluster perimeter to select an appropriate threshold for SIF data. Because of the 336 limited spatial coverage of the OCO-2 SIF data, we used an interval of 10 km to ensure 337 that each buffer had a sufficient number of observations. An example of the urban cluster 17 338 and its corresponding rural buffers was shown in Supplementary Fig. 11. Pixels that were 339 water body or crops were excluded from this analysis (Supplementary Fig. 11a). 340 Croplands were excluded because their photosynthetic dynamics are strongly affected by 341 human management. We also excluded the pixels in the rural areas which had elevations 342 greater than ±50 m of the average elevation of urban pixels (for both SIF and EVI data, 343 Supplementary Fig. 11b)21. Given that the number of OCO-2 SIF observations per year 344 was low, we opted for the use of SIF data in each urban cluster and its corresponding 345 rural buffers from 2015 to 2017 as one completed phenological cycle. OCO-2 Footprints 346 with viewing zenith angles (VZA) less than 20° in glint and target modes are used78. We 347 selected urban clusters with at least eight SIF observations in both the urban cluster and 348 each corresponding buffer to ensure a successful calculation of phenology. This resulted 349 in a total of 880 urban clusters (Supplementary Fig. 1). 350 Data 351 We used the OCO-2 SIF Lite product (B8100) at 757 nm for 2015–2017. OCO-2 is a 352 sun-synchronous polar-orbit satellite launched in July 2014, with an altitude of 705 km 353 and a descending node at approximately 13:30 local time18. OCO-2 leads the “A-Train” 354 satellite constellation with a repeating cycle of approximately 16 days. The instrument 355 collects high-resolution radiance spectra in the O2-A band (757–775 nm), which can be 356 used to retrieve fluorescence signals from vegetation. Retrievals based on the in-filling of 357 solar Fraunhofer lines in narrow spectral windows around 757 and 771 nm were 358 conducted using the singular vector decomposition method. The accuracy of single 18 359 measurements at 757 nm after a series of bias-correction and quality-control steps is 360 approximately 15-25% of the typical peak values of SIF18. However, because there are 361 dozens of single measurements in each urban cluster and the corresponding rural buffer, 362 most of the noise in the retrieved SIF has been reduced following the 1/√n law, which 363 means the precision errors in the average of n samples will get reduced by a factor of 364 1/√n79. Moreover, the relative uncertainties of single measurement in the urban clusters 365 and rural buffers are almost the identical (Supplementary Fig. 12); therefore, the retrieval 366 errors of SIF should have not large impacts on the urban–rural phenological patterns. We 367 chose the SIF data at 757 nm because the signals are approximately 1.5-fold greater at 368 757 than those at 771 nm. 369 The data for Aqua MODIS EVI collection-6 (MYD13A2; 16-day composite) with a 370 spatial resolution of 1 km for 2015–2017 were used to extract greenness phenology 371 information. The EVI data partially eliminate the effects of canopy background and 372 remain sensitive to small changes in vegetation activity compared to NDVI (normalized 373 vegetation index) data20. Thus, the EVI data are considered to be more appropriate for use 374 in areas with sparse vegetation coverage, such as urban clusters21,38. We used the Aqua 375 EVI data to maintain consistency with the OCO-2 overpass time (both occur at 376 approximately 13:30 local time). To make a comparison, the MODIS EVI data are 377 selected according to the same OCO-2 overpasses. To test whether the observed 378 urban–rural phenological patterns is an artefact of the difference in spatial resolutions 379 between OCO-2 SIF and MODIS EVI, we firstly resampled both SIF and EVI data to a 19 380 spatial resolution of 5 km and reproduced the urban–rural phenological differences. The 381 results also demonstrated larger differences of photosynthetic phenology in urban–rural 382 gradients than the differences of greenness phenology (Supplementary Fig. 5) and thus 383 demonstrated the robustness of our analysis. 384 We obtained the urban–rural XCO2 gradients using reprocessed Lite files from the 385 OCO-2 product (version 8r) for 2015–2017. The data were bias-corrected using 386 observations from the Total Carbon Column Observing Network (TCCON). OCO-2 also 387 collects high-resolution radiance spectra at the 1.61 and 2.06 μm bands, except for the 388 O2-A band, which can be used to retrieve XCO280 based on the Atmospheric CO2 389 Observations from Space algorithm. Comparisons with the measurements at TCCON 390 sites demonstrated that the precision of the Lite products was approximately 1.5 ppm81. 391 We used the Aqua MODIS collection-6 land-surface temperature (LST) data set 392 (MYD11A2; 8-day composite), with a spatial resolution of 1 km for 2015–2017 to 393 calculate the urban–rural LST gradient. We chose LST data from the Aqua satellite 394 because its data were acquired at 13:30 and 01:30 local time each day, approximately 395 representing the diurnal and nocturnal temperatures, respectively (i.e., maximum and 396 minimum temperatures). We also used the air surface temperature (AST) data from 397 TerraClimate dataset82, with a spatial resolution of 4 km to calculate the urban–rural AST 398 gradient. 399 We also extracted the mean annual air temperature (2 m above the surface) for each 400 urban cluster from the analytical data set based on both the Global Historical Climatology 20 401 Network (version 2) and the Climate Anomaly Monitoring System (GHCN-CAMS), with 402 a spatial resolution of 0.5°. The annual precipitation of each urban area was calculated 403 based on the monthly Integrated Multi-satellite Retrievals (IMERG) data set from global 404 precipitation measurements (GPMs), with a spatial resolution of 0.1°. Digital elevation 405 models (DEMs) at a spatial resolution of 30 arc seconds from the GTOPO30 data set and 406 the land-cover data set for 2015 from the European Space Agency (ESA) Climate Change 407 Initiative were also used in this study. We also used the global gridded nitrogen dioxide 408 (NO2) and ozone (O3) observations from Ozone Monitoring Instrument (OMI). We 409 calculated the mean LST changes for 2015–2035 at mid- to high-latitudes in the NH 410 (≥30°) from four Earth system models (CSIRO-Mk3.6.083, GFDL-CM384, GISS-E2-H85 411 and NorESM1-M86) within CMIP5. Projected CO2 concentration changes from 2015 412 were obtained from the RCP Database. Detailed information on the data sets is presented 413 in Supplementary Table 1. 414 Extraction of phenological data 415 We used a double-sigmoidal function to fit the annual cycles of the SIF and EVI 416 observations87. Choosing an appropriate method to fit noisy SIF data is important because 417 of the wide variation in phenology when estimated with different methods. We chose the 418 double-sigmoidal method for three reasons. First, the speeds of spring leaf out and 419 autumn senescence are not always the same: leaf out is generally faster than senescence88. 420 The double-sigmoidal method can identify the actual biophysical dynamics of vegetation 421 by setting two amplitudes for spring and autumn87. Second, functional fitting methods 21 422 have advantages for estimating phenology with noisy data89. The OCO-2 data are 423 acquired at different time intervals and are relatively noisy, and the number of 424 observations is small. Thus some fitting approaches, such as polynomial fitting, logistic 425 fitting and harmonic analysis, cannot be used or cannot converge to the global optimum. 426 Third, the double-sigmoidal method has been successfully applied to SIF data in a 427 previous study33. The fitting equation is56,87: a2 a3 428 y (t ) = a1 + − 1 + exp( − d1 (t − b1 )) 1 + exp( − d 2 (t − b2 )) (1) 429 where y(t) represents the observed SIF or EVI for a given day of year (DOY; t); ɑ1 430 represents the value in the winter dormant period; ɑ2 represents the value at the spring 431 and early summer plateau; ɑ3 represents the value at the late summer and autumn plateau; 432 b1 and b2 are the DOY mid-points of transitions for spring leaf out and autumn 433 senescence, respectively; and d1 and d2 are the corresponding slope coefficients of these 434 transitions. 435 The fitting equation was adapted to our SIF and EVI time series, and the coefficients 436 were optimized using genetic algorithms. A simple weighting scheme was applied to the 437 noisy data. We assigned a weight of two if a central point was within ±50% of the median 438 for a moving window of three points33. A maximum of 20 fitting attempts with different 439 initial values were made until the fitting procedure converged to the global optimum. 440 Each attempt had a maximum of 2000 iterations. The R2 values of SIF fitting were lower 441 than those of EVI, but most were larger than 0.8 (Supplementary Fig. 13). The R2 values 442 were generally high for the EVI data. The phenological indicators were then determined 22 443 using56: 4.562 SOS = b1 − 2d1  1  POS = 1.317  − + b1   a1 444  4.562 EOS = b2 + (2) 2d 2 GSL = EOS − SOS 445 where SOS, POS, EOS and GSL represent the start, peak, end and length of the growing 446 season, respectively. An example of the fitting method and the estimation of phenological 447 indicators is shown in Supplementary Fig. S14a. We excluded some extreme values using 448 a loose threshold to reduce uncertainties and slightly adjusted the thresholds based on 449 previous studies 21,38,41 due to the differences between the data and methods. SOSs earlier 450 than the 30th or later than the 180th DOY, and EOSs earlier than the 240th or later than the 451 350th DOY, were thus excluded from our analysis. To test whether our results were 452 affected by the weighting scheme, we also calculated the phenological indicators after 453 changing the moving-window and without using the weighting scheme. Results showed 454 no significant differences compared to the original method (Supplementary Fig. S14 b-c). 455 The urban–rural phenological gradients were also evident using these different weighting 456 schemes (Supplementary Fig. S14 d-e). 457 Conversion of urban–rural XCO2 gradients to surface CO2 gradients 458 The XCO2 data are column-averaged CO2 mixed ratios and differ in magnitude from 459 CO2 concentrations near the surface. Therefore, in order to investigate the effect of CO2 23 460 fertilization on vegetation phenology, we need to convert the urban–rural XCO2 gradients 461 to surface CO2 gradients. Although atmospheric transport models would be the ideal 462 approach to fulfill this purpose, there are still many precluding issues. Firstly, high-spatial 463 resolution prior fluxes of CO2, including the anthropogenic emissions and terrestrial 464 carbon uptake, as well as the high-resolution meteorological data (wind speed and wind 465 direction), still have large uncertainties, especially in urban areas with high-variable 466 emission sources (such as traffic pollution) and high degrees of surface heterogeneity90,91. 467 Secondly, as one of the main factors that affect the variations of surface CO2 468 concentration in finer scales, the heterogeneity of vertical profile in urban areas induced 469 by buildings can hardly be accurately estimated and illustrated92. Thirdly, continuous 470 high-spatiotemporal ground and satellite CO2 measurements are still insufficient to 471 provide accurate constraints and evaluations for atmospheric transport models to simulate 472 the high-spatial resolution surface CO2 concentrations at large scales93. Therefore, global 473 near-surface CO2 concentrations from atmospheric transport models or flux-inversion 474 systems generally have relatively low spatial resolutions94, which cannot meet our needs 475 for deriving urban–rural gradients at large scales. Statistical methods are therefore more 476 appropriate for converting XCO2 gradients to surface CO2 gradients. 477 We used two approaches to determine the conversion factors. First, we used a 478 near-surface CO2 data set with a comprehensive spatial coverage. The data were obtained 479 in Shanghai during April and May of 2014 and covered the entire urban area, with a total 480 of 172 sample points95 (Supplementary Fig. 15). The slopes of the near-surface CO2 24 481 measurements and OCO-2 XCO2 observations in the same city were calculated 482 (Supplementary Fig. 16a-b). The conversion factor could then be obtained (K = 21.2). 483 Second, we collected a number of urban–rural surface CO2 gradients from previous 484 studies (Supplementary Table 2) and then calculated the XCO2 gradients of these urban 485 clusters and compared them via linear regression (Supplementary Fig. 16c), which 486 indicated that the conversion factor was approximately 28.5. These data from previous 487 literatures have an average of seven years gap with the satellite observations, and thus the 488 observations of urban–rural surface CO2 gradients may be lower than those of present, 489 which possibly leads to a higher estimation of the K factor. Therefore, we used a K factor 490 of 25, which is the average value of these two approaches. Then we converted the 491 urban–rural XCO2 gradients to surface CO2 gradients through this conversion factor. This 492 conversion factor was derived via statistical methods that may have some uncertainties. 493 For example, the K factor derived from the first approach was 21.2 ± 8.3; while that from 494 the second approach was 28.5 ± 9.5 (mean ± standard error). To test the robustness of our 495 results according to different K values, we repeated our analysis using a wide range of K 496 factors (from 19 to 30), which was the intersection of ranges of the K factors estimated 497 from these two approaches. Results showed that the partial correlations between 498 urban–rural phenological gradients and CO2 gradients were not affected by K values, 499 given that this conversion factor was linear (Supplement Fig. 17a). Although the 500 predicted GSL had a declining trend when the K was increasing (Supplement Fig. 17b), 501 our predicted GSL shifts (5.9–11.4 days·decade–1) based on a K factor of 25 were 25 502 comparable to the average value of these shifts. These results suggested that the joint 503 control of elevated CO2 and temperature on urban–rural photosynthetic phenological 504 gradients was robust across different K values. 505 Analysis 506 The phenological differences, LST gradients and CO2 gradients between urban and 507 rural areas for each city were first calculated as: 508 ΔP =Purban − Pr ural (3) 509 where ∆P represents the urban–rural gradients for the phenological indicators (SOS, POS, 510 EOS and GSL), LST and CO2 concentration; and Purban and Prural represent these values in 511 urban clusters and their corresponding rural buffers, respectively. The LST differences 512 included the mean spring LST (January to May, corresponding to SOS), mean summer 513 LST (June to August, corresponding to POS), mean autumn and winter LST (September 514 to December, corresponding to EOS) and mean annual LST (corresponding to GSL)38. 515 Due to the limited OCO-2 observations, we used the whole year XCO2 observations and 516 the conversion factor (K) to derive the urban–rural CO2 gradients. The differences in 517 photosynthetic (SIF) and greenness (EVI) phenology along the urban–rural gradients for 518 different rural buffers (i.e., 10, 20, 30 and 40 km) were calculated (Supplementary Fig. 4). 519 Previous studies have shown that the mean footprint of urban areas on phenology was 520 about 20–25 km, which should be able to reflect the background vegetation 521 phenology21,38,96. In this analysis, we found that the largest urban–rural phenological 522 differences occurred between urban clusters and 30-km rural buffers. Therefore, we 26 523 therefore used the 30-km rural buffers as representatives in the subsequent analyses. 524 The spatial distribution of photosynthetic phenology differences based on latitude, 525 mean annual air temperature and precipitation were then analysed. We also analysed the 526 phenological differences in different climate zones based on the Köppen–Geiger climate 527 zones classification (http://people.eng.unimelb.edu.au/mpeel/koppen.html). To determine 528 the main controlling drivers of urban–rural phenological gradients, we determined the 529 partial correlations of ∆P with latitude, longitude, urban size, mean annual air 530 temperature, annual precipitation, CO2 gradient (∆CO2) and LST gradient (∆LST). The 531 partial correlation (two-tailed) of each factor was determined while controlling the other 532 factors. In order to reduce the stochastic error, the LST and CO2 gradients are binned 533 every 0.15°C and 1.5 ppm, respectively. Then we determined the main factors controlling 534 each phenological indicator by analysing the associations between urban–rural 535 phenological gradients with temperature and CO2 enhancements. Furthermore, we 536 derived their sensitivities with a forward stepwise regression model (Supplementary 537 Table 3). This method automatically repeats the procedures of forward selection and 538 backward elimination. In each step, an explanatory variable is considered to either add or 539 subtract from the set of variables based on a sequence of F tests. Finally, we estimated the 540 near-term projected shifts in photosynthetic phenology by multiplying the sensitivities 541 with LST and CO2 concentration predictions during 2015–2035 based on four RCP 542 scenarios (RCP 2.6, 4.5, 6.0 and 8.5). 543 Finally, we analysed the uncertainty of projected phenophase shifts as a consequence 27 544 of using LST other than air surface temperature (AST) to characterise the urban–rural 545 temperature gradients. Vegetation phenology theoretically should respond to AST more 546 than to LST. As shown in previous studies on the impact of land cover on temperature97, 547 the spatial gradients of AST are likely smaller than that of LST (Supplementary Fig. 18a), 548 which may translate in a higher sensitivity of GSL to AST than to LST (Supplementary 549 Fig. 18b). Therefore, projected shifts of GSL based on urban–rural AST gradients might 550 be larger than that based on LST gradients (Supplementary Fig. 18c; approximately 551 0.5–0.8 days·decade–1). However, current high-resolution AST datasets are generally 552 based on observations from weather stations, which cannot fully capture the small scale 553 urban–rural heterogeneity because of the low density of stations. Due to the limitation in 554 the current AST datasets, analysis based on this latter variable may actually have larger 555 uncertainties than those based on LST and will need more investigations in the future. 556 Code availability 557 The codes used to estimate the phenological indicators in this study are available at 558 https://drive.google.com/drive/folders/1yzcoRAjjubiLDqlg6zbLUCfE_m1mHAsL?usp=s 559 haring. 560 Data availability 561 OCO-2 SIF and XCO2 data are available at https://disc.gsfc.nasa.gov/. MODIS EVI 562 and LST data are available at https://ladsweb.modaps.eosdis.nasa.gov/. Precipitation data 563 can be obtained from https://pmm.nasa.gov. GHCN-CAMS air temperature data are 564 available at https://www.esrl.noaa.gov. CCI Land cover data are available at 28 565 http://maps.elie.ucl.ac.be/CCI/viewer/download.php. GTOPO30 DEM data are available 566 at https://earthexplorer.usgs.gov/. NO2 and O3 data are available at 567 http://www.temis.nl/index.php. Projected CO2 concentrations can be obtained from the 568 RCP Database (http://www.iiasa.ac.at/web-apps/tnt/RcpDb). Model results and the urban 569 clusters are available at 570 https://drive.google.com/drive/folders/1yzcoRAjjubiLDqlg6zbLUCfE_m1mHAsL?usp=s 571 haring. 572 References 573 1 Canadell, J. G. et al. Contributions to accelerating atmospheric CO2 growth from 574 economic activity, carbon intensity, and efficiency of natural sinks. Proceedings of 575 the national academy of sciences 104, 18866-18870 (2007). 576 2 Keenan, T. F. et al. Net carbon uptake has increased through warming-induced 577 changes in temperate forest phenology. Nature Climate Change 4, 598 (2014). 578 3 Peñuelas, J. & Filella, I. Phenology feedbacks on climate change. Science 324, 579 887-888 (2009). 580 4 Richardson, A. D. et al. Climate change, phenology, and phenological control of 581 vegetation feedbacks to the climate system. Agricultural and Forest Meteorology 582 169, 156-173 (2013). 583 5 Richardson, A. D. et al. Ecosystem warming extends vegetation activity but 584 heightens vulnerability to cold temperatures. Nature, 585 doi:10.1038/s41586-018-0399-1 (2018). 29 586 6 Peñuelas, J. & Filella, I. Responses to a warming world. Science 294, 793-795 587 (2001). 588 7 Piao, S. et al. Net carbon dioxide losses of northern ecosystems in response to 589 autumn warming. Nature 451, 49 (2008). 590 8 Piao, S. et al. Leaf onset in the northern hemisphere triggered by daytime 591 temperature. Nature communications 6, 6911 (2015). 592 9 Fu, Y. H. et al. Declining global warming effects on the phenology of spring leaf 593 unfolding. Nature 526, 104-107 (2015). 594 10 IPCC. Climate change 2014: synthesis report. Contribution of Working Groups I, 595 II and III to the fifth assessment report of the Intergovernmental Panel on Climate 596 Change. (IPCC, 2014). 597 11 Calfapietra, C. et al. Challenges in elevated CO2 experiments on forests. Trends in 598 plant science 15, 5-10 (2010). 599 12 Fu, Y. H. et al. Larger temperature response of autumn leaf senescence than spring 600 leaf‐out phenology. Global change biology (2017). 601 13 Wolkovich, E. M. et al. Warming experiments underpredict plant phenological 602 responses to climate change. Nature 485, 494 (2012). 603 14 Calfapietra, C., Peñuelas, J. & Niinemets, Ü. Urban plant physiology: 604 adaptation-mitigation strategies under permanent stress. Trends in plant science 605 20, 72-75 (2015). 606 15 Peng, S. et al. Surface urban heat island across 419 global big cities. 30 607 Environmental science & technology 46, 696-703 (2011). 608 16 Schwandner, F. M. et al. Spaceborne detection of localized carbon dioxide 609 sources. Science 358, eaam5782 (2017). 610 17 Zhao, S., Liu, S. & Zhou, D. Prevalent vegetation growth enhancement in urban 611 environment. Proceedings of the National Academy of Sciences 113, 6313-6318, 612 doi:10.1073/pnas.1602312113 (2016). 613 18 Sun, Y. et al. Overview of Solar-Induced chlorophyll Fluorescence (SIF) from the 614 Orbiting Carbon Observatory-2: Retrieval, cross-mission comparison, and global 615 monitoring for GPP. Remote Sensing of Environment 209, 808-823 (2018). 616 19 Guanter, L. et al. Global and time-resolved monitoring of crop photosynthesis 617 with chlorophyll fluorescence. Proceedings of the National Academy of Sciences 618 111, E1327-E1333 (2014). 619 20 Huete, A. et al. Overview of the radiometric and biophysical performance of the 620 MODIS vegetation indices. Remote sensing of environment 83, 195-213 (2002). 621 21 Zhou, D., Zhao, S., Zhang, L. & Liu, S. Remotely sensed assessment of 622 urbanization effects on vegetation phenology in China's 32 major cities. Remote 623 Sensing of Environment 176, 272-281 (2016). 624 22 Zhou, D., Zhao, S., Liu, S., Zhang, L. & Zhu, C. Surface urban heat island in 625 China's 32 major cities: Spatial patterns and drivers. Remote Sensing of 626 Environment 152, 51-61 (2014). 627 23 Chuine, I., Morin, X. & Bugmann, H. Warming, photoperiods, and tree phenology. 31 628 Science 329, 277-278 (2010). 629 24 Fu, Y. H., Campioli, M., Deckmyn, G. & Janssens, I. A. Sensitivity of leaf 630 unfolding to experimental warming in three temperate tree species. Agricultural 631 and Forest Meteorology 181, 125-132 (2013). 632 25 Jeong, S.-J. et al. Application of satellite solar-induced chlorophyll fluorescence 633 to understanding large-scale variations in vegetation phenology and function over 634 northern high latitude forests. Remote sensing of environment 190, 178-187 635 (2017). 636 26 Kikuzawa, K. Phenological and morphological adaptations to the light 637 environment in two woody and two herbaceous plant species. Functional Ecology 638 17, 29-38 (2003). 639 27 Daumard, F. et al. A field platform for continuous measurement of canopy 640 fluorescence. IEEE Transactions on geoscience and Remote Sensing 48, 641 3358-3368 (2010). 642 28 Suni, T. et al. Interannual variability and timing of growing‐season CO2 643 exchange in a boreal forest. Journal of Geophysical Research: Atmospheres 108 644 (2003). 645 29 Medvigy, D., Jeong, S. J., Clark, K. L., Skowronski, N. S. & Schäfer, K. V. 646 Effects of seasonal variation of photosynthetic capacity on the carbon fluxes of a 647 temperate deciduous forest. Journal of Geophysical Research: Biogeosciences 648 118, 1703-1714 (2013). 32 649 30 Chen, X., Wang, D., Chen, J., Wang, C. & Shen, M. The mixed pixel effect in 650 land surface phenology: A simulation study. Remote Sensing of Environment 211, 651 338-344 (2018). 652 31 Filella, I., Penuelas, J., Llorens, L. & Estiarte, M. Reflectance assessment of 653 seasonal and annual changes in biomass and CO2 uptake of a Mediterranean 654 shrubland submitted to experimental warming and drought. Remote Sensing of 655 Environment 90, 308-318 (2004). 656 32 Hilker, T. et al. Remote sensing of photosynthetic light-use efficiency across two 657 forested biomes: Spatial scaling. Remote Sensing of Environment 114, 2863-2874 658 (2010). 659 33 Walther, S. et al. Satellite chlorophyll fluorescence measurements reveal large‐ 660 scale decoupling of photosynthesis and greenness dynamics in boreal evergreen 661 forests. Global change biology 22, 2979-2996 (2016). 662 34 Porcar-Castell, A. et al. Linking chlorophyll a fluorescence to photosynthesis for 663 remote sensing applications: mechanisms and challenges. Journal of experimental 664 botany 65, 4065-4095 (2014). 665 35 Baker, N. R. Chlorophyll fluorescence: a probe of photosynthesis in vivo. Annu. 666 Rev. Plant Biol. 59, 89-113 (2008). 667 36 Norton Alexander, J., Rayner Peter, J. & Koffi Ernest, N. Assimilating 668 solar-induced chlorophyll fluorescence into the terrestrial biosphere model 669 BETHY-SCOPE: Model description and information content. 1-26 (2017). 33 670 37 White, M. A., Nemani, R. R., Thornton, P. E. & Running, S. W. Satellite evidence 671 of phenological differences between urbanized and rural areas of the eastern 672 United States deciduous broadleaf forest. Ecosystems 5, 260-273 (2002). 673 38 Zhang, X., Friedl, M. A., Schaaf, C. B., Strahler, A. H. & Schneider, A. The 674 footprint of urban climates on vegetation phenology. Geophysical Research 675 Letters 31 (2004). 676 39 Han, G. & Xu, J. Land surface phenology and land surface temperature changes 677 along an urban–rural gradient in Yangtze River Delta, China. Environmental 678 management 52, 234-249 (2013). 679 40 Cong, N. et al. Spring vegetation green-up date in China inferred from SPOT 680 NDVI data: A multiple model analysis. Agricultural and Forest Meteorology 165, 681 104-113 (2012). 682 41 Li, X. et al. Response of vegetation phenology to urbanization in the 683 conterminous United States. Global change biology 23, 2818-2830 (2017). 684 42 Niu, S. et al. Seasonal hysteresis of net ecosystem exchange in response to 685 temperature change: patterns and causes. Global Change Biology 17, 3102-3114 686 (2011). 687 43 Van der Tol, C., Verhoef, W. & Rosema, A. A model for chlorophyll fluorescence 688 and photosynthesis at leaf scale. Agricultural and forest meteorology 149, 96-105 689 (2009). 690 44 Sun, Y. et al. OCO-2 advances photosynthesis observation from space via 34 691 solar-induced chlorophyll fluorescence. Science 358, eaam5747 (2017). 692 45 Zarco-Tejada, P., Morales, A., Testi, L. & Villalobos, F. Spatio-temporal patterns 693 of chlorophyll fluorescence and physiological and structural indices acquired 694 from hyperspectral imagery as compared with carbon fluxes measured with eddy 695 covariance. Remote Sensing of Environment 133, 102-115 (2013). 696 46 Yang, X. et al. Solar‐induced chlorophyll fluorescence that correlates with 697 canopy photosynthesis on diurnal and seasonal scales in a temperate deciduous 698 forest. Geophysical Research Letters 42, 2977-2987 (2015). 699 47 Joiner, J. et al. The seasonal cycle of satellite chlorophyll fluorescence 700 observations and its relationship to vegetation phenology and ecosystem 701 atmosphere carbon exchange. Remote Sensing of Environment 152, 375-391 702 (2014). 703 48 CaraDonna, P. J., Iler, A. M. & Inouye, D. W. Shifts in flowering phenology 704 reshape a subalpine plant community. Proceedings of the National Academy of 705 Sciences 111, 4916-4921 (2014). 706 49 Menzel, A. et al. European phenological response to climate change matches the 707 warming pattern. Global change biology 12, 1969-1976 (2006). 708 50 Thompson, R. & Clark, R. Is spring starting earlier? The Holocene 18, 95-104 709 (2008). 710 51 Miller-Rushing, A. J. & Primack, R. B. Global warming and flowering times in 711 Thoreau's Concord: a community perspective. Ecology 89, 332-341 (2008). 35 712 52 Vitasse, Y. et al. Leaf phenology sensitivity to temperature in European trees: Do 713 within-species populations exhibit similar responses? Agricultural and forest 714 meteorology 149, 735-744 (2009). 715 53 Piao, S. et al. Weakening temperature control on the interannual variations of 716 spring carbon uptake across northern lands. Nature Climate Change 7, 359 717 (2017). 718 54 Körner, C. & Basler, D. Phenology under global warming. Science 327, 719 1461-1462 (2010). 720 55 Way, D. A. & Montgomery, R. A. Photoperiod constraints on tree phenology, 721 performance and migration in a warming world. Plant, Cell & Environment 38, 722 1725-1736 (2015). 723 56 Gonsamo, A., Chen, J. M. & Ooi, Y. W. Peak season plant activity shift towards 724 spring is reflected by increasing carbon uptake by extratropical ecosystems. 725 Global change biology (2017). 726 57 Xu, C., Liu, H., Williams, A. P., Yin, Y. & Wu, X. Trends toward an earlier peak 727 of the growing season in Northern Hemisphere mid‐latitudes. Global change 728 biology 22, 2852-2860 (2016). 729 58 Wu, C. et al. Contrasting responses of autumn-leaf senescence to daytime and 730 night-time warming. Nature Climate Change, 1 (2018). 731 59 Marchin, R. M., Salk, C. F., Hoffmann, W. A. & Dunn, R. R. Temperature alone 732 does not explain phenological variation of diverse temperate plants under 36 733 experimental warming. Global change biology 21, 3138-3151 (2015). 734 60 Sigurdsson, B. D. Elevated CO2 and nutrient status modified leaf phenology and 735 growth rhythm of young Populus trichocarpa trees in a 3-year field study. Trees 736 15, 403-413 (2001). 737 61 Liu, Q. et al. Delayed autumn phenology in the Northern Hemisphere is related to 738 change in both climate and spring phenology. Global change biology 22, 739 3702-3711 (2016). 740 62 Cleland, E. E., Chiariello, N. R., Loarie, S. R., Mooney, H. A. & Field, C. B. 741 Diverse responses of phenology to global changes in a grassland ecosystem. 742 Proceedings of the National Academy of Sciences 103, 13740-13744 (2006). 743 63 Jach, M. E. & Ceulemans, R. Effects of elevated atmospheric CO2 on phenology, 744 growth and crown structure of Scots pine (Pinus sylvestris) seedlings after two 745 years of exposure in the field. Tree physiology 19, 289-300 (1999). 746 64 Taylor, G. et al. Future atmospheric CO2 leads to delayed autumnal senescence. 747 Global Change Biology 14, 264-275 (2008). 748 65 NOAA National Centers for Environmental Information, State of the Climate: 749 Global Climate Report for May 2018, published online June 2018, retrieved on 750 June 24, 2018 from https://www.ncdc.noaa.gov/sotc/global/201805. 751 66 Jeong, S. J., HO, C. H., GIM, H. J. & Brown, M. E. Phenology shifts at start vs. 752 end of growing season in temperate vegetation over the Northern Hemisphere for 753 the period 1982–2008. Global change biology 17, 2385-2399 (2011). 37 754 67 Barichivich, J. et al. Large‐scale variations in the vegetation growing season and 755 annual cycle of atmospheric CO2 at high northern latitudes from 1950 to 2011. 756 Global change biology 19, 3167-3183 (2013). 757 68 Buyantuyev, A. & Wu, J. Urbanization diversifies land surface phenology in arid 758 environments: interactions among vegetation, climatic variation, and land use 759 pattern in the Phoenix metropolitan region, USA. Landscape and Urban Planning 760 105, 149-159 (2012). 761 69 Decina, S. M., Templer, P. H. & Hutyra, L. R. Atmospheric Inputs of Nitrogen, 762 Carbon, and Phosphorus across an Urban Area: Unaccounted Fluxes and Canopy 763 Influences. Earth's Future 6, 134-148 (2018). 764 70 Gregg, J. W., Jones, C. G. & Dawson, T. E. Urbanization effects on tree growth in 765 the vicinity of New York City. Nature 424, 183 (2003). 766 71 Zhang, X., Friedl, M. A., Schaaf, C. B., Strahler, A. H. & Liu, Z. Monitoring the 767 response of vegetation phenology to precipitation in Africa by coupling MODIS 768 and TRMM instruments. Journal of Geophysical Research: Atmospheres 110 769 (2005). 770 72 Do, F. C. et al. Environmental influence on canopy phenology in the dry tropics. 771 Forest Ecology and Management 215, 319-328 (2005). 772 73 Peñuelas, J. et al. Complex spatiotemporal phenological shifts as a response to 773 rainfall changes. New Phytologist 161, 837-846 (2004). 774 74 Fisher, J. B., Huntzinger, D. N., Schwalm, C. R. & Sitch, S. Modeling the 38 775 terrestrial biosphere. Annual Review of Environment and Resources 39, 91-123 776 (2014). 777 75 Zhou, Y. et al. A cluster-based method to map urban area from DMSP/OLS 778 nightlights. Remote Sensing of Environment 147, 173-185 (2014). 779 76 Zhou, Y. et al. A global map of urban extent from nightlights. Environmental 780 Research Letters 10, 054011 (2015). 781 77 Huang, X., Schneider, A. & Friedl, M. A. Mapping sub-pixel urban expansion in 782 China using MODIS and DMSP/OLS nighttime lights. Remote Sensing of 783 Environment 175, 92-108 (2016). 784 78 Zhang, Y. et al. On the relationship between sub-daily instantaneous and daily 785 total gross primary production: Implications for interpreting satellite-based SIF 786 retrievals. Remote Sensing of Environment 205, 276-289, 787 doi:10.1016/j.rse.2017.12.009 (2018). 788 79 Frankenberg, C. et al. Prospects for chlorophyll fluorescence remote sensing from 789 the Orbiting Carbon Observatory-2. Remote Sensing of Environment 147, 1-12 790 (2014). 791 80 Crisp, D. et al. The on-orbit performance of the Orbiting Carbon Observatory-2 792 (OCO-2) instrument and its radiometrically calibrated products. Atmospheric 793 Measurement Techniques 10, 59 (2017). 794 81 Wunch, D. et al. Comparisons of the Orbiting Carbon Observatory-2 (OCO-2) 795 XCO2 measurements with TCCON. Atmospheric Measurement Techniques 10, 39 796 2209 (2017). 797 82 Abatzoglou, J. T., Dobrowski, S. Z., Parks, S. A. & Hegewisch, K. C. 798 TerraClimate, a high-resolution global dataset of monthly climate and climatic 799 water balance from 1958–2015. Scientific data 5, 170191 (2018). 800 83 Collier, M. et al. in International Congress on Modelling and 801 Simulation–MODSIM. 802 84 Griffies, S. M. et al. The GFDL CM3 coupled climate model: characteristics of 803 the ocean and sea ice simulations. Journal of Climate 24, 3520-3544 (2011). 804 85 Schmidt, G. A. et al. Configuration and assessment of the GISS ModelE2 805 contributions to the CMIP5 archive. Journal of Advances in Modeling Earth 806 Systems 6, 141-184 (2014). 807 86 Bentsen, M. et al. The Norwegian earth system model, NorESM1-M—Part 1: 808 Description and basic evaluation of the physical climate. Geosci. Model Dev 6, 809 687-720 (2013). 810 87 Gonsamo, A., Chen, J. M. & D’Odorico, P. Deriving land surface phenology 811 indicators from CO2 eddy covariance measurements. Ecological indicators 29, 812 203-207 (2013). 813 88 Elmore, A. J., Guinn, S. M., Minsley, B. J. & Richardson, A. D. Landscape 814 controls on the timing of spring, autumn, and growing season length in mid‐ 815 Atlantic forests. Global Change Biology 18, 656-674 (2012). 816 89 Hird, J. N. & McDermid, G. J. Noise reduction of NDVI time series: An empirical 40 817 comparison of selected techniques. Remote Sensing of Environment 113, 248-258 818 (2009). 819 90 Duren, R. M. & Miller, C. E. Measuring the carbon emissions of megacities. 820 Nature Climate Change 2, 560 (2012). 821 91 Bréon, F. et al. An attempt at estimating Paris area CO 2 emissions from 822 atmospheric concentration measurements. Atmospheric Chemistry and Physics 15, 823 1707-1724 (2015). 824 92 Chen, F. et al. The integrated WRF/urban modelling system: development, 825 evaluation, and applications to urban environmental problems. International 826 Journal of Climatology 31, 273-288 (2011). 827 93 Chevallier, F. et al. Toward robust and consistent regional CO2 flux estimates 828 from in situ and spaceborne measurements of atmospheric CO2. Geophysical 829 Research Letters 41, 1065-1070 (2014). 830 94 Peters, W. et al. An atmospheric perspective on North American carbon dioxide 831 exchange: CarbonTracker. Proceedings of the National Academy of Sciences 104, 832 18925-18930 (2007). 833 95 Liu, M. et al. Spatial variation of near-surface CO2 concentration during spring in 834 Shanghai. Atmospheric Pollution Research 7, 31-39 (2016). 835 96 Zhou, D., Zhao, S., Zhang, L., Sun, G. & Liu, Y. The footprint of urban heat 836 island effect in China. Scientific reports 5, 11160 (2015). 837 97 Alkama, R. & Cescatti, A. Biophysical climate impacts of recent changes in 41 838 global forest cover. Science 351, 600-604, doi:10.1126/science.aac8083 (2016). 839 840 42 841 Acknowledgements 842 This research was financially supported by the National Key R&D Program of China 843 (2016YFA0600202), the Strategic Priority Research Program of the Chinese Academy of 844 Sciences under Grants XDA19040500, Jiangsu Provincial Natural Science Fund for 845 Distinguished Young Scholars of China (BK20170018), International Cooperation and 846 Exchange Programs between NSFC and DFG (41761134082) and General Program of 847 National Science Foundation of China (41671421). JP would like to acknowledge the 848 financial support from the European Research Council Synergy grant 849 ERC-SyG-2013-610028 IMBALANCE-P. AH acknowledges financial support from 850 Australian Research Council Discovery Program grant DP170101630. S.H.W was 851 supported by the Postgraduate Research & Practice Innovation Program of Jiangsu 852 Province (KYCX18_0037) and the Key Research Program of the Chinese Academy of 853 Sciences (Grant NO., KFZD-SW-310). 854 Author contributions 855 Y.G.Z. designed the research; S.H.W. performed the analysis; S.H.W., Y.G.Z. and W.M.J. 856 drafted the paper; J.P., and A.C. contributed to the interpretation of the results and to the 857 writing of the paper; A.H., Y.Y.Z. and Y.S.F. contributed to the writing of the paper; Y.Y.Z. 858 and M. L. provided the data. 859 Competing interests 860 The authors declare no competing interests. 861 43 862 Fig. 1. Mean urban–rural phenological differences based on SIF and EVI. SOS, start of 863 growing season; POS, peak of growing season; EOS, end of growing season; GSL, length 864 of growing season. The operator ∆ represents the differences in phenological indicators 865 between urban areas and rural buffers (e.g., ∆SOS = SOSurban − SOSrural). The error bars 866 represent standard errors of the mean (n = 880). 867 Fig. 2. Spatial distributions of ∆SOS (a), ∆POS (b), ∆EOS (c) and ∆GSL (d) for SIF (n = 868 880). The first column represents the urban effects on phenological metrics for various 869 urban clusters. The second column represents phenological differences based on 870 latitudinal profiles (the shaded areas represent the standard errors of the mean); whereas 871 the third column shows the average metrics for various climatic backgrounds (mean 872 annual air temperature and annual precipitation). 873 Fig. 3. Partial correlations of SIF (a) and EVI (b) urban–rural phenological differences 874 (∆SOS, ∆POS, ∆EOS and ∆GSL) with explanatory variables. Lat: latitude; Lon: longitude; 875 Size: log 10 of the urban size; Air T: mean annual air temperature; Pre: mean annual 876 precipitation; ∆CO2: urban–rural CO2 difference; ∆Tday: daytime urban–rural land surface 877 temperature difference; ∆Tnight: night-time urban–rural land surface temperature difference; 878 SOS: start of growing season; POS: peak of growing season; EOS: end of growing season; 879 GSL: length of growing season. *: p < 0.05; **: p < 0.01 (n = 880). 880 Fig. 4. Associations between elevated temperature (a–d), CO2 (e–h) and urban–rural 881 photosynthetic phenological gradients. ∆CO2, urban–rural CO2 gradient; ∆LSTday, 882 urban–rural gradient of daytime LST. The LST and CO2 gradients are binned every 44 883 0.15°C and 1.5 ppm, respectively; each dot represents a bin. The solid lines indicate 884 significant trends and the shaded areas represent the 95% confidence intervals, and p 885 values represent the significance (Student’s t-tests, n = 880). 886 Fig. 5. Projected shifts in photosynthetic phenology during the next two decades. 887 Sensitivities of photosynthetic phenology to temperature (a) and CO2 (b). The error bars 888 represent standard errors of the mean. c–f, Shifts of photosynthetic SOS, POS, EOS and 889 GSL from 2015 to 2035, as calculated from the sensitivities shown in a–b and LST and 890 CO2 concentration predictions during 2015–2035 from four “Representative 891 Concentration Pathway” (RCP) scenarios. The shaded areas represent the standard errors. 892 g, Projected shifts of photosynthetic phenology during the next two decades in per decade 893 (days·decade–1). The negative values represent the advances of SOS and POS, and 894 positive values represent a delayed EOS and a prolonged GSL. The error bars represent 895 standard errors of the mean. 45