Tree species composition governs urban phenological responses to warming

Domain: earth_science
Published: 2025-04-18
Editorial outcome: major_revision

Abstract
Urban environments are typically warmer than surrounding rural areas, providing a unique setting for studying phenological responses to climate warming. Phenological differences between urban and rural trees are driven by local climate and species composition. Yet, the extent to which species composition influences phenological responses to urbanization remains poorly understood. To address this, we combine manipulative experiments, satellite-derived phenology data, and georeferenced tree occurrence records. Our findings show that, across Northern Hemisphere cities, differences in the temperature sensitivity of spring phenology between urban and rural areas are largely driven by urban-rural variation in species composition, surpassing the effects of preseason temperature. This pattern is particularly pronounced in Asian cities, where urban areas exhibit 0.74 ± 0.24 days/°C higher temperature sensitivity than rural areas. In-depth analyses using experiments and high-resolution satellite imagery from Beijing further demonstrate species-specific phenological responses to urbanization, with urban-dominant species exhibiting higher temperature sensitivity in urban environments compared to rural ones. These findings show that both interspecific variation in temperature sensitivity and species-specific responses to urbanization contribute to the pronounced impact of species composition on urban-rural phenological patterns. Our study underscores the importance of considering species composition when studying phenological responses to climate warming, especially in urban contexts.

Introduction
Human-induced climate change has caused changes in plant phenology1–3, affecting the structure and function of terrestrial ecosystems4–6. In addition, the urban heat island effect has caused many urban areas to exhibit higher temperatures than adjacent rural regions7,8. These urbanized zones offer unique opportunities to study how ecosystems react to human disturbances and rising temperatures9,10. In recent years, urban tree phenology has received increasing attention due to its direct impact on water, energy and carbon fluxes, as well as its role in modifying the albedo effect11–13. Although the magnitude of urban temperature gradients is a key factor driving phenological differences between urban and rural areas, the composition of tree species–each with distinct responses to warming5,11,14–likely plays an equally important role. However, the relative influence of species composition and urban climate on phenological responses to urbanization remains poorly understood.

Given the large impact of temperature on plant physiology, it is key to understand how urbanization affects the temperature sensitivity of vegetation phenology. Recent satellite-derived data suggest that urban warming reduces the sensitivity of spring vegetation phenology to temperature10, meaning that the start of season advances by fewer days per degree of warming in urban areas compared to rural ones. However, studies have so far provided mixed evidence for this15,16, likely due to the substantial interspecific variation in phenological responses to temperature1,11,14. Urban-rural differences in species composition may therefore greatly impact the relationship between local temperatures and tree phenology.

Urban environments often host a wide range of species from diverse climatic and phylogenetic origins. Consequently, urban areas not only differ in species composition from their rural counterparts but also tend to exhibit higher species richness17,18. In natural forests, greater species richness has been shown to reduce phenological sensitivity to temperature, either directly through the presence of temperature-resilient species or indirectly via mechanisms such as deeper root systems and altered soil processes19–21. This raises an important question: do phenological differences between urban and rural areas primarily reflect variations in species identities, or are they driven by differences in species richness? Despite the importance of this question, interspecific variations in phenological sensitivity have not yet been thoroughly quantified in urban phenology studies11,22,23.

In this study, we analyzed ground-sourced and satellite-derived data from 135 largest cities in the Northern Hemisphere to explore how urban-rural differences in species composition, topography, and climate affect spring phenological responsiveness to urbanization. To uncover the mechanisms driving the observed large-scale patterns, we conducted experiments in Beijing, focusing on phenological differences between urban and rural areas across four tree species commonly found in both settings. Of these species, three are predominantly urban, while one is more common in rural areas. Moreover, using high-resolution remote sensing data from Beijing, we compared how urbanization influences spring phenology at both the broader vegetation and individual species levels. We expected that species composition would influence the temperature sensitivity of spring phenology (ST, expressed as change in days per degree of warming, d °C−1) due to interspecific differences in ST or distinct responses of species to urbanization. Therefore, urban-rural differences in species composition should contribute to variation in ST between urban and rural vegetation.

Results
We investigated the temperature sensitivity of tree phenology across urban and rural areas by studying phenological variation in 135 Northern Hemisphere cities24. Using the MODIS Global Vegetation Phenology product25, we obtained start-of-season (SOS) dates from 2003 to 2020 at 500m spatial resolution. For each city, we distinguished between urban and rural areas based on predefined city boundaries and proportionate outward extensions (see “Methods” and Supplementary Fig. 1a)10,26. In 97% of the cities, spring temperatures in urban areas were significantly higher than those in rural counterparts (p<0.1), with an average difference of 0.30 ± 0.02 °C, ranging from 0.1 °C to 2.2 °C (Supplementary Fig. 2). Consistent results were observed for mean annual temperatures (Supplementary Fig. 3). In 66% of cities, the SOS occurred earlier in urban areas than in rural areas, with an average difference of 6.5 ± 0.9 days (Fig. 1a and Supplementary Fig. 4). This finding is consistent with a previous study, which reported a difference of 6.1 ± 0.7 days10.

We calculated the ST of each pixel by analyzing the linear relationship between each year’s SOS and preseason temperature from 2003 to 202015,27, which reflects the change in the SOS date per degree of warming (days °C−1). For 60% of the analyzed pixels, preseason temperature significantly (p<0.1) affected SOS (Supplementary Fig. 5). Subsequent analyses focused on these significant pixels, though results remained consistent when all pixels were included (Supplementary Figs. 5, 6). ST was higher in urban than in rural areas for 51% of cities, with an average difference of 0.69 ± 0.24 d °C−1 (Fig. 1b and Supplementary Fig. 5a, b). Using all pixels, ST was higher in urban than in rural areas for 66% of cities, with an average difference of 0.35 ± 0.11 d °C−1 (Supplementary Figs. 5c, d, 6b). The urban-rural difference in ST was especially pronounced in Asian cities (0.74 ± 0.24 d °C−1) but was minimal in European (0.11 ± 0.27 d °C−1) and American (0.14 ± 0.09 d °C−1) cities (Fig. 1b, Supplementary Figs. 6b, 7), aligning with prior regional studies10,16.

To understand how species composition might influence observed trends, we extracted 3.5 million urban tree observations from the Global Biodiversity Information Facility (GBIF)28, the most comprehensive global urban biodiversity data network currently available. We used the beta diversity index to quantify species composition differences between urban and rural areas of each city29,30. These differences were most pronounced in Asian cities, then American cities, with European cities showing the least variation (Figs. 1c and Supplementary Fig. 8). A higher beta diversity strongly increased urban-rural differences in SOS (ΔSOS, Figs. 1d and Supplementary Fig. 6d) and ST (ΔST, Fig. 1e and Supplementary Figs. 6e, 9a, b). The pronounced ΔST observed in Asia may thus result from substantial differences in species composition between urban and rural areas. This is consistent with prior studies highlighting extensive landscaping efforts in Chinese cities that introduce exotic tree species and enhance urban-rural differences in species composition17,18. As a result, cities with high beta diversity tend to exhibit more pronounced variation in phenological responses between urban and rural areas. Given that greater species diversity is associated with enhanced ecosystem resilience to climate warming21, cities with higher beta diversity, such as those in Asia, may have an advantage in adapting to future climate changes. However, further studies are needed to comprehensively understand how urban beta diversity impacts the resilience of cities to climate change.

To understand the key drivers of phenological variation between urban and rural areas, we quantified the effects of beta diversity, species richness, climate, vegetation characteristics, and population density on ΔSOS and ΔST. Both the differences (Δ) between urban and rural areas and city-level averages were considered. After evaluating multicollinearity among variables using variance inflation factors31,32, we selected nine key variables for our analysis. Using random forest models32, urban-rural differences in preseason temperature (ΔTemperature) was identified as the top driver of ΔSOS (Supplementary Fig. 10a). Beta diversity emerged as the top driver of ΔST, followed by ΔTemperature (Fig. 2a). Urban-rural species richness differences (ΔRichness) had a small effect on ΔSOS and ΔST, suggesting that phenological differences between urban and rural areas are driven primarily by variations in species composition rather than species richness itself. Further examination through recursive partitioning analyses32,33 confirmed beta diversity as the top driver of ΔST, with ΔTemperature and mean annual temperature influencing ΔST under certain conditions (Fig. 2b). Specifically, in cities with high dissimilarity in species composition between urban and rural areas (beta diversity ⩾ 0.73), a warmer mean annual temperature is associated with a larger ΔST. Conversely, in cities with lower beta diversity (<0.73), a larger ΔTemperature leads to an increased ΔST.

To test whether the definition of urban and rural areas affected our findings and mitigate a potential bias due to fewer tree observation records in rural areas, we conducted an additional analysis comparing the inner, highly urbanized areas of cities to their outer, more rural parts within the city limits (inward delineation, Supplementary Fig. 1b). This approach revealed consistent patterns and corroborated the initial analysis, which treated all areas within city boundaries as urban (Supplementary Figs. 9c, d, 10–13). These findings show that differences in species composition between urban and rural areas largely drive variation in spring phenological responsiveness.

To investigate the mechanisms by which beta diversity influences the responsiveness of spring phenology to urbanization, we conducted species-level experiments in Beijing, China. With a population of ~21.9 million, Beijing has undergone extensive urban expansion and warming in recent decades8,34. The Beijing was divided into six concentric zones (Zones 1 to 6) based on its ring roads (Supplementary Fig. 14). The city center (zone 1) of Beijing is 0.4 °C (March to May temperature, Spring T) warmer than its most rural area (zone 6, Fig. 3a). To analyze spring phenology along this urbanization gradient, we used the Landsat-based enhanced vegetation index (EVI) from 2003–2020. The SOS was defined as the date when EVI first crossed 50% of its annual peak (SOS50). To exclude potential elevational effects due to mountains in the west of Beijing, we removed pixels with elevations higher than the mean ± 1.5 standard deviations of all pixels (268 m a.s.l.). In urban areas, the average SOS occurred on May 3, 15 days ahead of the SOS in rural areas on May 18 (Fig. 3b). Consistent with the Northern Hemisphere-wide trends (Fig. 1b), urban areas in Beijing displayed greater ST, increasing from 2.1 to 4.4 d °C−1 compared to rural areas (Fig. 3c).

We assessed the urbanization impact on spring leaf-out at the species level by collecting twigs from three species dominant in urban areas (Ginkgo biloba, Acer truncatum and Platanus occidentalis) and one species common in rural areas (Fraxinus chinensis). Samples were collected from both rural (Jiufeng Mountain, outside zone 6) and urban (Beijing Zoo, in zone 2) locations (see “Methods” and Supplementary Fig. 14). These twigs, collected at two different dormancy stages, were exposed to the full combination of four temperature (5, 10, 15 and 20 °C) and two photoperiod treatments (8 and 16 h) in controlled climate chambers35. Contrary to the vegetation-level observations, no difference in ST was found between urban and rural plants when grouping all species (Fig. 3e), implying that the vegetation-level results are driven by urban-rural variation in species composition. Furthermore, the experiments showed that the effect of urbanization on ST differed significantly among species (Fig. 3e). Except for Platanus, which showed similar ST in urban and rural areas (Urban ST = −6.2 ± 0.4 °C−1; Rural ST = −6.4 ± 0.4 d °C−1), urban-dominant species (Ginkgo and Acer) exhibited higher ST in urban than in rural areas (Ginkgo: Urban ST = −5.7 ± 0.4; Rural ST = −5.1 ± 0.3 d °C−1; Acer: Urban ST = −6.5 ± 0.2; Rural ST = −5.7 ± 0.1 °C−1). In contrast, the rural-dominant species exhibited the opposite response (Fraxinus: Urban ST = −5.5 ± 0.4; Rural ST = −7.6 ± 0.4 d °C−1). ANOVA results confirmed a significant effect of species on ST, and, importantly, also showed a strong interaction between species and urbanization, suggesting that the effect of urbanization on ST depends on the species (Supplementary Table. S1). This supports the fundamental role of species composition in driving variation in ST between urban and rural areas.

Using high-resolution Landsat imagery, we validated our experimental findings by analyzing satellite-derived spring phenology for each of the four tree species in both urban (zone 1 and 2) and rural areas (zones 5 and 6). The satellite-derived results were consistent with the experiment, with Ginkgo and Acer showing significantly higher ST in urban than in rural areas, while Fraxinus showed the opposite trend (Fig. 3e). Additionally, our twig-cutting experiment revealed varied effects of photoperiod and chilling on species’ spring phenology and their degree-day requirements (Supplementary Figs. 15, 16). Specifically, Fraxinus was influenced by day length and unaffected by chilling, whereas Ginkgo and Acer were more strongly affected by chilling, possibly accounting for the observed variation in ST among the species36,37. These results highlight the substantial influence of species-specific phenological differences on the satellite-derived results at the vegetation level.

Based on our species-level observational and experimental data, we propose two mechanisms for how species-specific phenological responsiveness to urbanization may affect the urban-rural variation in vegetation-level ST. First, inherent temperature sensitivities vary among species, meaning that differences in species composition between urban and rural areas directly contribute to differences in temperature sensitivity (Fig. 4a). Second, the effect of urbanization on ST depends on the species. For instance, should urbanization-driven warming reduce ST of species A but increase it in species B (Fig. 4b), vegetation-level ST remains constant across urban and rural areas when both species are equally abundant in urban and rural areas (Fig. 4c, Scenario 1, ST_urban = ST_rural). Conversely, ST will decline in urban areas if species A is predominantly found there (Scenario 2, ST_urban < ST_rural), whereas ST will increase in urban areas if species B is more dominant (Scenario 3, ST_urban > ST_rural).

These findings demonstrate that variation in species-specific temperature sensitivities, combined with species-specific responses to urbanization, drive the observed urban-rural differences in vegetation-level phenological temperature sensitivity. While we aimed to include species representative of the dominant trees in urban and rural settings, we acknowledge that other species may exhibit different patterns, warranting further investigation.

Methods
Analysis of satellite-derived phenology observations at city level

Urban boundary data. The urban area boundaries were determined using the latest version of the Global Urban Boundaries (GUB) dataset (http://data.ess.tsinghua.edu.cn)24. The dataset was obtained from 30m resolution Landsat imagery, employing a hierarchical approach to improve the homogeneity of built-up areas within urban centers while preserving the heterogeneity of built-up areas at the urban fringes24. We extracted data for the 300 large cities across the Northern Hemisphere, spanning from 30° to 90°N latitude. After quality checks and data matching with phenology datasets, a total of 285 cities, with an area between 157 to 5117 km2, were included in further analysis.

Tree species observation data. Our georeferenced species’ occurrences were obtained from the Global Biodiversity Information Facility (GBIF; https://www.gbif.org), an international network and data infrastructure that provides open access to data encompassing all forms of life on Earth28. The GBIF dataset stands as the world’s largest biodiversity data network45 and has been extensively used in biodiversity and phenology studies46, allowing us to conduct this study at a global scale. We obtained tree species occurrences following these steps: i) we removed the records from living specimens and conserved specimens; ii) queried all tracheophyte species and extracted georeferenced data records, excluding any coordinates with an abnormal latitude and longitude range; iii) standardized and harmonized scientific names of tracheophytes using the R package U.Taxonstand47; iv) extracted all tree species occurrence records for the urban and rural areas of the selected cities.

Phenology data. Start-of-season (SOS) dates for the period 2003–2020 were obtained from the MODIS Global Vegetation Phenology product (MCD12Q2 V6)25. This product provides phenological metrics at a spatial resolution of 500m for every vegetated pixel on land, derived from time series of the 2-band Enhanced Vegetation Index (EVI2). SOS was defined as the date when EVI2 first crossed 50% of the segment EVI2 amplitude greenness. We extracted multi-year SOS data for all pixels within the urban-rural extent. To reduce uncertainties, we removed pixels i) in which SOS occurred before March 1 or after May 31, and ii) with phenology information available for fewer than 15 years.

Climate geography data. We obtained daily near-surface air temperature data at 1-km resolution from ref. 48. Additionally, 30-year mean annual temperature and precipitation data were obtained from CHELSA V2.1 (Climatologies at High resolution for the Earth Land Surface Areas, https://chelsa-climate.org/)49, hosted by the Swiss Federal Institute for Forest, Snow, and Landscape Research WSL. Elevation data was derived from EarthEnv (https://www.earthenv.org/topography), which is based on the digital elevation model products of global 250m GMTED2010 and near-global 90m SRTM4.1dev50. Human density data were derived from the Gridded Population of the World product51, using 2020 as the reference year.

Vegetation data. The vegetation classification dataset was derived from the MODIS Land Cover Type (MCD12Q1) Version 6.1 data product (https://lpdaac.usgs.gov/products/mcd12q2v061/). This product provides global land cover types at yearly intervals, with 2011 used as the reference year in this study. Vegetation types, including cultivated croplands and grasslands, were excluded from the analyses10. Enhanced Vegetation Index (EVI) data were derived from Aqua Moderate Resolution Imaging Spectroradiometer (MODIS) Vegetation Indices (MYD13Q1) Version 6.

Statistical analysis. To compare vegetation responses between urban and rural areas, we defined urban areas as the regions within the established city boundaries for each city. Rural areas were defined as a proportionate outward extension of these boundaries, meaning the rural area boundary is symmetrically and evenly expanded in size and shape to cover the same area as the predefined urban portion of the city (outward delineation, Supplementary Fig. 1a)10,26. Alternatively, we used inward delineation, where urban and rural areas were represented by inner and outer polygons of equivalent area within the city boundary (inward delineation, Supplementary Fig. 1b). This dual delineation allowed us to test the robustness of our findings and account for potential biases arising from variations in tree observation records and sampling efforts across urban and rural zones. These paired urban and rural areas were used in subsequent analyses. To avoid potential bias from small sample sizes of tree observations in some cities, we excluded cities with fewer than 500 tree records in either rural or urban areas. Consequently, 135 cities (>1000 tree records) were retained for the outward strategy, and 115 cities for the inward strategy, both of which were used in the final analysis.

To compare the difference in species composition between urban and rural areas, we calculated beta diversity as the Bray-Curtis difference using the number of individuals per species30. The R package vegan 2.6-4 provides functions for calculating beta diversity52.

We defined the preseason as the 90 days preceding the mean SOS date for each pixel53, and calculated the mean temperature during the preseason for each year. We calculated the ST of each pixel by analyzing the linear relationship between each year’s SOS and preseason temperature from 2003 to 202015,27, which reflects the change in the SOS date per degree of warming (days °C-1). We further calculated the mean SOS and ST for both rural and urban areas (SOSrural, SOSurban, ST_rural and ST_urban). The difference in SOS (ΔSOS) and ST (ΔST) was calculated to present the variation in SOS and ST between urban and rural areas. Based on the sample size (number of cities), we grouped the beta diversity index into 5 or 6 subsets at different beta diversity intervals to examine the effects of beta diversity on ΔSOS and ΔST (Supplementary Fig. 9). Variance inflation factor (VIF) analysis was used to examine the multicollinearity among variables32, including beta diversity, and both the urban-rural differences (Δ) and the city-level averages of climate factors (preseason temperature, mean annual temperature, and mean annual precipitation), topography (elevation and area), vegetation (species richness and EVI), and human-related factors (population density) between urban and rural areas. After excluding highly correlated variables, all VIFs were smaller than 5, suggest sufficient independence among the predictors. We selected nine key variables for our analysis and assessed their relative importance using random forest models32. To further explore the interactive effects among the main drivers of ΔST, we conducted recursive partitioning analyses using the ‘partykit’ and ‘ggparty’ packages32,33. Using a decision tree algorithm, we investigated the context dependency for the effects of the variables on absolute ΔST. The minimum node size was set to ~20% (27 data points for the outward strategy and 23 for the inward strategy) of the total city number.

Climate-controlled experiment and high-resolution satellite-based observations at species level

Study area. Beijing is one of the largest cities around the world with a population of 21.9 million (115.7-117.4°E, 39.4-41.6°N, Supplementary Fig. 14), and has experienced extensive expansion and intensive warming in recent decades8,34. Beijing is characterized by a warm, temperate monsoon climate, with an average annual temperature of 11–13 °C and an average annual precipitation of ~600 mm. Precipitation is high in summer, while winters are cold and dry. The dominant vegetation type is broad-leaved deciduous forest. Based on the ring roads of Beijing, we divided the study area into six concentric zones (Zones 1 to 6, Supplementary Fig. 14) to roughly reflect the degrees of urbanization following previous studies34,54.

Vegetation phenology using remote sensing. The SOS was estimated using long-term Landsat images from 2003 to 2020 via the Google Earth Engine (GEE) platform following Li et al. (2017, 2019). We used atmospherically corrected surface reflectance Landsat images (Level 2, Collection 2) from different sensors (Thematic Mapper, Enhanced Thematic Mapper Plus and Operational Land Imager) at 30m spatial and 16-day temporal resolution. Previous studies indicated that the enhanced vegetation index (EVI) improves the estimation of spring phenology relative to the Normalized Difference Vegetation Index (NDVI) in temperate regions55, and we thus used the EVI to quantify spring phenology. We then created a composite of the EVI time series for each pixel56,57, and applied a masking procedure to remove clouds and shadows58.

To capture the seasonal patterns of vegetation dynamics, we selected clear-sky observations of EVI at each pixel (30m × 30m) and ordered them by the day of year (DOY). To minimize the impact of abnormal observations, we applied a moving average method to the continuous observations. A double logistic model was then applied to fit the time series of EVI curves to characterize the seasonal change of vegetation growth as follows:

f(t) = v1 + v2 * (1 / (1 + e^(-m1 * (t - n1))) - 1 / (1 + e^(-m2 * (t - n2))))

where f(t) represents the fitted EVI value at DOY t; v1 and v2 refer to the background and amplitude of EVI over the year; n1 and n2 refer to the dates of maximum increase or decrease in EVI based on sigmoid curves, and m1 and m2 are the slopes that determine the shape of the sigmoid curves. The parameters were estimated using a stepwise statistical approach57.

Consistent with MODIS phenology dataset, the SOS was defined as the date when the EVI first reaches 50% of its annual amplitude (SOS50). We extracted vegetation SOS for Beijing from MODIS and Landsat data, and matched their resolutions to 500m (Supplementary Figs. 17a, b). The SOS derived from Landsat was highly correlated with that from MODIS, indicating consistency between the two datasets (r=0.73, Supplementary Fig. 17c). To exclude potential elevational effects due to mountains in the west of Beijing, we removed pixels with elevations higher than the mean ± 1.5 standard deviations of all pixels (268 m a.s.l.). We used the MODIS Land Cover Type data product to exclude cultivated croplands and grasslands. Daily near-surface air temperatures were extracted from a global gridded dataset with a spatial resolution of 1km covering the period from 2003 to 2020 (https://doi.org/10.25380/iastate.c.6005185)48. The elevation data for Beijing was derived from the 30m resolution NASADEM dataset (https://earthdata.nasa.gov/esds/competitive-programs/measures/nasadem), which was used to exclude mountain pixels.

Species-level phenology using a manipulative experiment. Twigs of four tree species were collected from a rural (Jiufeng Mountain, outside zone 6, 30km away from the city center) and an urban area (Beijing Zoo in the city center in zone 2, Supplementary Fig. 14). The temperature sensors installed at the two collection sites showed that the daily air temperature in the urban area was 4.5 °C higher than in the rural area, with an average annual temperature of 13.6 and 9.1 °C, respectively (Supplementary Fig. 18).

Based on the List of Major Native Tree Species in Beijing and GBIF observations, we selected four common temperate deciduous tree species: three of them dominant in urban areas (Ginkgo biloba, Acer truncatum, and Platanus occidentalis) and one predominant in rural areas (Fraxinus chinensis) (hereafter referred to as Ginkgo, Acer, Fraxinus and Platanus). We selected twelve healthy adult individuals for each of the four species, with six individuals sampled from Jiufeng Mountain and six from Beijing Zoo, yielding a total of 48 tree individuals (6 individuals × 4 species × 2 sites). All the selected trees were mature and had similar diameters at breast height. Eight healthy twigs ~40 cm in length were cut from each individual tree (from the southern side of the canopy) in December 2020 and January 2021, respectively. We disinfected the twigs with commercial hypochlorite solution and then placed them in 395 ml plastic bottles filled with 265 ml cool tap water. After pre-processing, we transferred the twigs to climate chambers, recut them by ~2 cm at the base every two weeks to avoid the occlusion of vessels, and changed the water in the bottles every one to two weeks37.

The experiment consisted of 16 combinations of environmental treatments (2 chilling × 4 temperature × 2 photoperiod) using 768 twigs. Four constant temperature treatments, i.e., 5 °C (T5), 10 °C (T10), 15 °C (T15) and 20 °C (T20) were implemented in four climate chambers. In each chamber, two photoperiod treatments (8 and 16 h) were applied to test whether photoperiod had any effect on budburst time by covering half the space of the chambers with shade cloth every day35. Furthermore, two chilling treatments were applied by collecting twigs earlier and later in winter, causing different exposure of twigs to outdoor winter conditions37,59 (December 31, 2020 [Low chilling, Clow] and January 19, 2021 [High chilling, Chigh]). The environmental treatment of this experiment lasted from January 1, 2020, to June 3, 2021. Temperature sensors (HOBO MX2202, Onset Computer Corporation, Bourne, MA, USA) were used to record the air temperature of each environmental treatment every 30 minutes (Supplementary Fig. 19). Illuminance in all chambers was ~8944 lux (161 μmol m−2 s−1), the average concentration of CO2 was ~442 ppm, the relative air humidity was ~40%.

For each cultivated twig, we monitored the development of buds every three days and recorded the budburst date, which was defined as the date when the bud scale opened and the leaf was partially visible37,60. Spring phenology was then represented using the days to achieve budburst, defined as the number of days from the date of treatment start to the budburst date. The heat requirement (growing degree days, GDD) for each twig to achieve budburst was calculated using the mean daily temperature, with a base temperature set at 0 °C.

Species-level phenology using remote sensing. To validate the results of the twig-cutting experiment, we used 30m resolution Landsat images to extract the SOS for the same species and conducted further species-level analysis. We selected 3-10 patches (dominated by one of the four species used in the twig-cutting experiment) at rural areas (zones 5 and 6; Fenghuang Mountain, 116°05'30“E, 40°06'46“N; Jiufeng Mountain, 116°5'11“E, 40°03'42“N; Xishan Forest Park, 116°12'36“E, 39°58'40“N) and urban areas (zone 2; Beijing Normal University, 116°22'20“E, 39°58'02“N; Beijing Zoo, 116°20'36“E, 39°56'52“N; Yuyuantan Park, 116°19'35“E, 39°55'21“N). Latitude and longitude of the patches were recorded using a hand-held GPS to extract Landsat and climate data for the respective sites. We then extracted the spring phenology dates of each patch for each species using the same methodology used for the vegetation-level phenology analysis.

Statistical analysis. We first analyzed the effects of urbanization on spring temperature (from March to May, Spring T), spring vegetation phenology (SOS) and the sensitivity of spring phenology to temperature (ST) along the urbanization gradients. To compare the vegetation-level phenology results with species-level patterns, we additionally analyzed the effects of urbanization on spring phenology using the observations from the manipulative experiment and species-level vegetation phenology. The temperature sensitivity of spring phenology (ST, defined as the advancement in budburst date per degree warming, days °C-1) was calculated using linear regression analysis35. For the manipulative experiment, the air temperature of each chamber was used to calculate the ST. For the satellite-derived SOS, we defined the preseason as the 90 days preceding the mean SOS date for each pixel53, and calculated the mean temperature during the preseason for the ST calculation. The difference in spring phenology and growing degree days between rural and urban areas was tested using the paired t-test, ANCOVA was used to test for ST differences among treatments. ANOVA was conducted to test the effects of species and urbanization on ST, with species, urbanization, chilling and photoperiod treatments as explanatory variables and ST as response variable. Two models were employed: Model 2 included these variables as predictors, while Model 1 extended to involve the interaction of species and site as predictors (Supplementary Table. 1).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This paper explored how special compositions influence the vegetation phenology in cities. The topic is interesting and\ncrucial since the special differences and compositions have commonly been questioned in urban vegetation phenology\nstudies. Multi-source data, including experiments, satellite-derived phenology data, and georeferenced tree occurrence data,\nwere jointly used to answer these questions. In particular, experiments with twig cutting and climate chambers were made to\nenhance the proposed assumption. The paper is well written with clear logic, and I have outlined comments below, hope\nthey are helpful to improve this paper and would like to review this paper after revision.\n1) The experiment is awesome and highly recommended to answer some essential questions. My comment is the length (or\nduration) of the experiment is adequate to draw the conclusion here. I didn’t find specific information about them, and I may\nwonder if the duration is too short; the results may be biased. Please clarify it in the paper.\n2) Another concern is the MODIS and Landsat phenology data. My question is whether the different sources of datasets may\nresult in varying outputs, as the previous paper highlighted. It is better to do some comparison among them, thus making the\nregional findings and local experiments consistent. Also, the 500m resolution is adequate to reflect the tree species?\n3) Page 105, 1.9(0.9) days between urban and rural areas, seems a little bit tiny. It is better to compare and discuss with\nexisting publications.\n4) Page 106, is the temperature can well characterize the difference in the urban domain and surrounding rural areas? I\nhope to seem some results about the spatial pattern of temperature during the crucial growth stages of vegetation.\n5) Page 165: How to ensure these two species can be used to represent the richness in the urban and surrounding rural\nareas? It is better to provide some numbers to convince us they are dominatent and are good for representation.\n6) Page 138, how to define this value of 0.73. Please clarify.\n7) Page 285, “proportionate”, unclear.\n8) In the Beijing case, are those pixels locating in the mountain (west side of Beijing) excluded before the analysis", "reviewer_2": "The current research investigates how species composition affects phenological responses to urbanization through\nexperiments, satellite-derived phenology data, and georeferenced tree occurrence records. While the topic and experiments\nare intriguing, the findings may not be surprising, as phenological responses often vary significantly by species. Additionally,\nthere may be uncertainties in the evidence due to the different measures used. Here are some major concerns:\n1. The study's primary innovation is its focus on species composition, yet it only includes trees. Other vegetation types, such\nas crops and shrubs/grass, could significantly influence phenological patterns detected by remote sensing.\n2. The paper claims to explore phenological variation along urbanization gradients, but it only compares phenological\nvariation between urban and rural areas. The study does not clarify how urbanization gradients were defined or analyzed.\n3. The definition or calculation of temperature sensitivity of SOS is unclear. It can be assessed either temporally (using long-\nterm SOS data on a per-pixel basis) or spatially (using SOS data from all pixels in a given year), which have different\nimplications.\n4. The experimental design is confusing. The paper states that it selected two urban-dominated and two rural-dominated\nspecies, sampling six individuals from Jiufeng Mountain and six from Beijing Zoo. It is unclear how many individuals of each\nspecies were sampled and whether the species were sampled in both urban and rural areas or only in their respective\nenvironments.\n5. There appears to be a mismatch between the remote sensing and experimental findings. Remote sensing data show a\nrelationship between ΔST and beta diversity across cities, with higher beta diversity associated with larger ΔST. The\nenvironmental experiments in Beijing, however, focus on temperature sensitivity and do not clearly support the remote\nsensing findings or reflect the ΔST-beta diversity relationship.\n6. It is unclear why the study estimated Landsat-based SOS only from 2003 to 2015, rather than the entire study period.\nAdditionally, the rationale for re-estimating SOS with Landsat data in Beijing, as opposed to using MODIS data, is not\nexplained.\n7. Elevation likely influences vegetation phenology, particularly in mountainous cities, but there are no measures to address\npotential uncertainties related to elevation effects.", "reviewer_3": "This manuscript examines the impact of urban gradients and species composition on the temperature sensitivity of spring\nphenology in rural and urban areas. The authors use sophisticated methods to obtain a very robust data set, which showed\nthat the impact of urban gradients is a lesser driver of temperature sensitivity than species composition. While I believe this\ndata set and results may be informative to the scientific community, I was let down the overall framing of the paper and the\nexplanation of how these results help explain phenological trends and/or might be utilized. For example, are these\nparameters used in Earth System Models or could this information be used to by city planners to develop communities that\nwill be resilient to climate warming? More specific comments below.\nINTRODUCTION\nL73-75. I think a clearer description of what the authors mean by \"stronger phenological responses\" and \"decreases in\nphenological responses\" is needed here. Does a decreased response mean that trees will acclimate to increased\ntemperatures and their phenological timing will revert back what was observed 100 years ago?\nOverall, I think more information about the authors' expectations for changes in temperature sensitivity related to species\ncomposition is needed. This is such a broad term, so maybe examples would help. Is the thought that a greater variety of\nspecies will have a lower temperature sensitivity simply because there is a high likelihood that a subset of those species will\nhave an innate resilience to warming? Or maybe that some of those species will undergo increased transpiration rates and\ncause a cooling effect in their community?\nRESULTS AND DISCUSSION\nIn this section you state that, \"The pronounced ΔST in Asia may thus result from substantial differences in species\ncomposition between urban and rural areas. This is consistent with prior studies highlighting extensive landscaping efforts in\nChinese cities that increase urban tree species diversity and enhance urban-rural differences in species composition.\" This\nseems like an interesting finding, but I have trouble understanding how to interpret this result. What does the fact that many\nAsian cities have a high beta diversity, and thus large differences in SOS and ΔST, mean for the resilience of these cities to\nclimate change? How will the response of these cities compare to those with low beta diversity?\nOverall, I found this section to be heavy on the results and light on the discussion or interpretation of those results. Placing\nthese results in the context of real world examples or theoretical scenarios would allow the reader to see their importance."}, "author_response": "[Comment 1] This paper explored how special compositions influence the vegetation\nphenology in cities. The topic is interesting and crucial since the special differences and\ncompositions have commonly been questioned in urban vegetation phenology studies.\nMulti-source data, including experiments, satellite-derived phenology data, and\ngeoreferenced tree occurrence data, were jointly used to answer these questions. In\nparticular, experiments with twig cutting and climate chambers were made to enhance\nthe proposed assumption. The paper is well written with clear logic, and I have outlined\ncomments below, hope they are helpful to improve this paper and would like to review\nthis paper after revision.\n[Response 1] We thank the reviewer for the positive evaluation and constructive\ncomments on our manuscript. Please see below our responses to each comment.\n[Comment 2] The experiment is awesome and highly recommended to answer some\nessential questions. My comment is the length (or duration) of the experiment is\nadequate to draw the conclusion here. I didn’t find specific information about them, and\nI may wonder if the duration is too short; the results may be biased. Please clarify it in\nthe paper.\n[Response 2] Thank you for your positive comment on our experimental design. The\nexperiment started on January 1, 2021, and lasted for more than five months, until June\n3, 2021. This duration and extensive sampling of twig cuttings ensures robust results1–\n3. We clarified the experimental duration in the revised manuscript as follows: \"The\nenvironmental treatment of this experiment lasted from January 1, 2020, to June 3,\n2021.\" (L. 427-428).\n[Comment 3] Another concern is the MODIS and Landsat phenology data. My\nquestion is whether the different sources of datasets may result in varying outputs, as\nthe previous paper highlighted. It is better to do some comparison among them, thus\nmaking the regional findings and local experiments consistent. Also, the 500 m\nresolution is adequate to reflect the tree species?\n[Response 3.1] Following the reviewer's comment, we compared the MODIS and\nLandsat datasets and found that they are highly correlated (r = 0.67, Fig. S18),\nindicating consistency between the two datasets. This comparison has now been\nincorporated in the revised manuscript: “The SOS from Landsat was highly correlated\nwith the SOS from MODIS, indicating consistency between the two datasets (r = 0.67,\nFig. S18).” (L. 389-391)\nFig. S18 | Correlation between the start of season (SOS) from MODIS and Landsat datasets.\n[Response 3.2] We fully agree with the reviewer that the 500 m resolution is insufficient\nto reflect tree species in species-level analyses. This was not done here. Instead, we\nused the higher-resolution Landsat 30 m dataset, which better aligns with field data on\ntree species and provides a more detailed local context. This is now clarified in the\nrevised Methods: “To validate the results of the twig-cutting experiment, we used 30 m\nresolution Landsat images to extract the SOS for the same species and conducted further\nspecies-level analysis.” (L. 442-444)\n[Comment 4] Page 105, 1.9(0.9) days between urban and rural areas, seems a little bit\ntiny. It is better to compare and discuss with existing publications.\n[Response 4] The previously reported smaller SOS difference between urban and rural\nareas was averaged across all cities. In the revised manuscript, we updated our data and\nfound that in 67% of cities, the SOS occurs earlier in urban areas than in rural areas,\nwith an average difference of 6.5 ± 0.9 days. This result aligns closely with a previous\nstudy, which reported a similar difference of 6.1 ± 0.7 days4. We have now included a\ncomparison with this previous study in the revised manuscript: “In 67% of cities, the\nSOS occurred earlier in urban areas than in rural areas, with an average difference of\n6.5 ± 0.9 days (Figs. 1a, S4). This finding is consistent with a previous study, which\nreported a difference of 6.1 ± 0.7 days.” (L. 92-95)\n[Comment 5] Page 106, is the temperature can well characterize the difference in the\nurban domain and surrounding rural areas? I hope to seem some results about the spatial\npattern of temperature during the crucial growth stages of vegetation.\n[Response 5] Previous studies have demonstrated that temperature is a crucial factor in\ndistinguishing between urban and rural areas due to the urban heat island effect4,5. This\nhas also been demonstrated in our study, where the mean annual temperatures in urban\nareas are significantly higher than in their rural counterparts (Fig. S3). Following the\nreviewer’s suggestion, we analyzed the spatial pattern of spring temperature (March to\nMay) and observed consistent findings. This finding has now been highlighted in L. 89-\n92: “In 97% of the cities, spring temperatures in urban areas were significantly higher\nthan those in rural counterparts (p < 0.1), with an average difference of 0.29 ± 0.02°C,\nranging from 0.1°C to 2.2°C (Fig. S2). Consistent results were observed for mean\nannual temperatures (Fig. S3).”\n[figure redacted]\nFig. S2 | Spring temperature for urban and rural areas of the selected cities. Maps showing the\nmulti-year spring temperature for urban areas (a), rural areas (b), and the difference between urban\nand rural areas (c) using the outward delineation method.\n[Comment 6] Page 165: How to ensure these two species can be used to represent the\nrichness in the urban and surrounding rural areas? It is better to provide some numbers\nto convince us they are dominatent and are good for representation.\n[Response 6] Thank you for pointing this out. These four species were selected to\nexplore whether the urban-rural differences in temperature sensitivity of spring\nphenology are due to species-specific responses to urbanization, rather than to represent\nthe species richness in the urban and surrounding rural areas. Based on the List of Major\nNative Tree Species in Beijing and GBIF observations, more than 70% of Ginkgo\nbiloba, Acer truncatum, and Platanus occidentalis were widely recorded in urban\nsettings (urban-dominant species), while nearly 80% of Fraxinus chinensis were\ndistributed in rural landscapes (rural-dominant species). We have now improved the\nprevious classification: “Based on the List of Major Native Tree Species in Beijing and\nGBIF observations, we selected four common temperate deciduous tree species: three\nof them dominant in urban areas (Ginkgo biloba, Acer truncatum, and Platanus\noccidentalis) and one predominant in rural areas (Fraxinus chinensis) (hereafter\nreferred to as Ginkgo, Acer, Fraxinus and Platanus).” (L. 405-409).\nFurthermore, we now further discuss this in the revised manuscript: “While we aimed\nto include species representative of the dominant trees in urban and rural settings, we\nacknowledge that other species may exhibit different patterns, warranting further\ninvestigation.” (L. 214-217)\n[Comment 7] Page 138, how to define this value of 0.73. Please clarify.\n[Response 7] We utilized the beta diversity index to quantify differences in species\ncomposition between urban and rural areas of each city25,26. Beta diversity is calculated\nbased on the abundance of tree species (i.e., the number of individuals per species)26.\nThe beta-diversity index ranges from 0 to 1, with values closer to 1 indicating greater\ndifferences in tree species composition between urban and rural areas, and values closer\nto 0 representing more similar species composition. A beta diversity value of 0.73\nsuggests a relatively high dissimilarity in species composition, meaning that urban and\nrural areas have fewer shared species. We have clarified this in l. 139-141: “Specifically,\nin cities with high dissimilarity in species composition between urban and rural areas\n(beta diversity > 0.73), a warmer mean annual temperature is associated with a larger\nΔS .”\nT\n[Comment 8] Page 285, “proportionate”, unclear.\n[Response 8] The term \"proportionate\" refers to an extension of the rural area boundary\nthat is symmetrically and evenly matched in size and shape to the predefined urban area\n(Fig. S1a). To improve the clarity, we revised the expression as: “To compare vegetation\nresponses between urban and rural areas, we defined urban areas as the regions within\nthe established city boundaries for each city. Rural areas were defined as a proportionate\noutward extension of these boundaries, meaning the rural area boundary is\nsymmetrically and evenly expanded in size and shape to cover the same area as the\npredefined urban portion of the city (outward delineation, Fig. S1a)10,26.” (L. 316-321)\n[figure redacted]\nFig. S1 | Definitions of urban and rural areas based on city boundary delineation methods. a,\nOutward delineation: Urban areas (in red) are defined within the predefined city boundary, while\nrural areas (in light red) represent an outward buffer surrounding the city boundary, covering an area\nequivalent to the urban region. b, Inward delineation: Urban and rural areas are defined as the inner\n(red) and outer (light red) polygons within the city boundary, each covering equivalent areas. This\ninward approach mitigates potential biases due to uneven sampling in rural zones.\n[Comment 9] In the Beijing case, are those pixels locating in the mountain (west side\nof Beijing) excluded before the analysis.\n[Response 9] Thank you for this thoughtful comment. We now excluded mountain\npixels in the west of Beijing and found consistent results (Figs. 3, S15). This has been\nclarified in the revised manuscript (L. 161-162; 166-167; 391-392): “To exclude\npotential elevational effects due to mountains in the west of Beijing, we removed pixels\nlocated above 100 m or 500 m a.s.l.”\n[figure redacted]\nFigure 3. Experimental evidence for the effects of urbanization on spring\nphenology. a-c, Effects of urbanization on spring temperature (Spring T, a), the start-\nof-season (SOS, b) and the temperature sensitivity of SOS (S , c). Pixels 500 m above\nT\nsea level were excluded, respectively. The black lines indicate the rings 2-6 in Beijing,\nwhich divide the study area into six concentric zones (Zone 1 to Zone 6). d, Effects of\nurbanization on spring phenology based on a manipulative experiment at the species\nlevel using twig cuttings and climate chambers. The dots indicate the average time to\nbudburst of each species when growing in a rural or urban area, respectively. The\nsquares and error bars indicate the average budburst dates and its standard deviation\nacross species. e, Effects of urbanization on the sensitivity of spring phenology to\ntemperature using either twig cutting experiments or individual-based satellite images\nat the species level. The temperature sensitivities were estimated from the manipulative\nexperiments (solid bars) or the satellite-based observations (bars with blue lines). The\nerror bars indicate the standard error of the mean. ns, not significant; *, p < 0.1; **,\np < 0.05 and ***, p < 0.01. f, Effects of urbanization on spring phenology using individual-\nbased satellite observations at the species level. To obtain species-level observations,\npatches where one of the four tree species was dominant were selected. Dots refer to\nthe spring phenology of each patch, the semi-violin plots indicate the distribution of\nspring phenology, and the error bar refers to the standard deviation of spring phenology.\nGinkgo, Acer, Fraxinus and Platanus refer to Ginkgo biloba, Acer truncatum, Fraxinus\nchinensis and Platanus occidentalis, respectively.\n[figure redacted]\nFig. S15 | Effects of urbanization gradients on spring temperature (Spring T), the start of the\nvegetation season (SOS), and the temperature sensitivity of SOS (S ). Pixels 100 m above sea level\nT\nwere excluded, respectively. The black lines indicate the rings 2-6 in Beijing, which divide the study\narea into six concentric zones (Zone 1 to Zone 6).\nReviewer #2\n[Comment 1] The current research investigates how species composition affects\nphenological responses to urbanization through experiments, satellite-derived\nphenology data, and georeferenced tree occurrence records. While the topic and\nexperiments are intriguing, the findings may not be surprising, as phenological\nresponses often vary significantly by species. Additionally, there may be uncertainties\nin the evidence due to the different measures used.\n[Response 10] We appreciate Reviewer #2’s acknowledgment of the significance of\nour research topic and the insightful concerns. While phenological responses vary by\nspecies, as noted in our Introduction, prior studies on urban phenology have not\nexamined how urban-rural differences in species composition impact these responses.\nOur study addresses this gap by combining experiments, satellite-derived phenology\ndata and georeferenced tree occurrence records. We believe our findings provide\nvaluable insights into the field.\nFurthermore, we have now revised our introduction to further underscore the\nsignificance of these insights, as detailed in the revised manuscript (L. 50-68). Please\nfind our point-by-point responses to each comment below.\n“Recent satellite-derived data suggest that urban warming reduces the sensitivity of\nspring vegetation phenology to temperature10, meaning that the start of season advances\nby fewer days per degree of warming in urban areas compared to rural ones. However,\nstudies have so far provided mixed evidence for this15,16, likely due to the substantial\ninterspecific variation in phenological responses to temperature1,11,14. Urban-rural\ndifferences in species composition may therefore greatly impact the relationship\nbetween local temperatures and tree phenology.\nUrban environments often host a wide range of species from diverse climatic and\nphylogenetic origins. Consequently, urban areas not only differ in species composition\nfrom their rural counterparts but also tend to exhibit higher species richness17,18. In\nnatural forests, greater species richness has been shown to reduce phenological\nsensitivity to temperature, either directly through the presence of temperature-resilient\nspecies or indirectly via mechanisms such as deeper root systems and altered soil\nprocesses19–21. This raises an important question: do phenological differences between\nurban and rural areas primarily reflect variations in species identities, or are they driven\nby differences in species richness?Despite the importance of this question, interspecific\nvariations in phenological sensitivity have not yet been thoroughly quantified in urban\nphenology studies11,22,23.” (L. 50-67)\n[Comment 2] The study's primary innovation is its focus on species composition,\nyet it only includes trees. Other vegetation types, such as crops and shrubs/grass,\ncould significantly influence phenological patterns detected by remote sensing.\n[Response 11] Our study focused primarily on trees due to their key role in urban\necosystems and the relative ease of identifying tree species diversity. We acknowledge\nthe potential influence of other vegetation types on phenological patterns, which is an\nimportant area for future study. Following the reviewer’s comments, we now excluded\nother vegetation types (including crops and grasslands)4, which showed consistent\nresults (Figs. 1-3, S2-13). Additionally, we also mentioned the need to address other\nvegetation types in future studies in the Discussion: “Additionally, although our study\nfocused on trees due to their ecological importance and the relative ease of assessing\ntheir diversity, future research should consider other vegetation types to gain a more\ncomprehensive understanding of urban phenology.” (L. 260-262).\n[Comment 3] The paper claims to explore phenological variation along urbanization\ngradients, but it only compares phenological variation between urban and rural areas.\nThe study does not clarify how urbanization gradients were defined or analyzed.\n[Response 12] Thank you for highlighting this point. We have adjusted our language\nto accurately reflect our study's approach, which focuses on differences between urban\nand rural areas rather than analyzing urbanization gradients. In urban phenological\nstudies, it is common to define rural comparison zones by establishing a buffer around\nurban areas4,12. To strengthen our findings, we applied two methods for delineating\nurban and rural areas—outward and inward delineation (Fig. S1b). Both approaches\nyielded consistent results, supporting the robustness of our conclusions.\nAdditionally, we clarified our definitions and methods for distinguishing urban and\nrural areas and revised the description of these areas in relation to urbanization gradients\nas follows: “To compare vegetation responses between urban and rural areas, we\ndefined urban areas as the regions within the established city boundaries for each city.\nRural areas were defined as a proportionate outward extension of these boundaries,\nmeaning the rural area boundary is symmetrically and evenly expanded in size and\nshape to cover the same area as the predefined urban portion of the city (outward\ndelineation, Fig. S1a)10,26. Alternatively, we used inward delineation, where urban and\nrural areas were represented by inner and outer polygons of equivalent area within the\ncity boundary (inward delineation, Fig. S1b). This dual delineation allowed us to test\nthe robustness of our findings and account for potential biases arising from variations\nin tree observation records and sampling efforts across urban and rural zones. These\npaired urban and rural areas were used in subsequent analyses.” (L. 316-326)\n[figure redacted]\nFig. S1 | Definitions of urban and rural areas based on city boundary delineation methods. a,\nOutward delineation: Urban areas (in red) are defined within the predefined city boundary, while\nrural areas (in light red) represent an outward buffer surrounding the city boundary, covering an area\nequivalent to the urban region. b, Inward delineation: Urban and rural areas are defined as the inner\n(red) and outer (light red) polygons within the city boundary, each covering equivalent areas. This\ninward approach mitigates potential biases due to uneven sampling in rural zones.\n[Comment 4] The definition or calculation of temperature sensitivity of SOS is unclear.\nIt can be assessed either temporally (using long-term SOS data on a per-pixel basis) or\nspatially (using SOS data from all pixels in a given year), which have different\nimplications.\n[Response 13] The temperature sensitivity of SOS was quantified using long-term\n(2003-2020) SOS data as the regression coefficient, which reflects the change in the\nSOS date per degree of warming (days °C-1). This method has been widely used in\nphenological studies7,13. However, the indirect space-for-time method, which calculates\nthe temperature sensitivity of SOS using data from all pixels in a given year, can\npotentially lead to inaccurate findings because preseason lengths can vary significantly\namong pixels. Following the reviewer's comment, we clarified the definition of\ntemperature sensitivity of SOS in the revised manuscript: “We calculated the S of each\nT\npixel by analyzing the linear relationship between each year’s SOS and preseason\ntemperature from 2003 to 202015,27, which reflects the change in the SOS date per\ndegree of warming (days °C-1).” (L. 97-99; 335-338)\n[Comment 5] The experimental design is confusing. The paper states that it selected\ntwo urban-dominated and two rural-dominated species, sampling six individuals from\nJiufeng Mountain and six from Beijing Zoo. It is unclear how many individuals of each\nspecies were sampled and whether the species were sampled in both urban and rural\nareas or only in their respective environments.\n[Response 14] We apologize for the lack of clarity in our original description. In our\nstudy, we selected four common temperate deciduous tree species, sampling twelve\nindividuals of each species across both urban and rural areas. Specifically, we sampled\nsix individuals from Jiufeng Mountain (rural site) and six from Beijing Zoo (urban site),\nresulting in a total of 48 trees (6 individuals × 4 species × 2 sites). We have revised the\nmanuscript to clarify this experimental design.\nThe updated text reads: “We selected twelve healthy adult individuals for each of the\nfour species, with six individuals sampled from Jiufeng Mountain and six from Beijing\nZoo, yielding a total of 48 tree individuals (6 individuals × 4 species × 2 sites).” (L.\n409-411)\n[Comment 6] There appears to be a mismatch between the remote sensing and\nexperimental findings. Remote sensing data show a relationship between ΔST and beta\ndiversity across cities, with higher beta diversity associated with larger ΔST. The\nenvironmental experiments in Beijing, however, focus on temperature sensitivity and\ndo not clearly support the remote sensing findings or reflect the ΔST-beta diversity\nrelationship.\n[Response 15] Thank you for this comment. Our satellite-based analysis aimed to\nexamine the effect of urban-rural differences in species composition (beta diversity) on\nthe phenological response to temperature (ΔS ). However, our experiments were\nT\nspecifically designed to explore the underlying mechanisms of how higher beta\ndiversity contributes to increased ΔS . The experimental results demonstrated species-\nT\nspecific phenological responses to urbanization: urban-dominant species exhibited\nhigher S in urban areas compared to rural ones, whereas rural-dominant species\nT\nshowed the opposite pattern. Furthermore, we observed a significant interaction\nbetween species type and urbanization, suggesting that the impact of urbanization on\ntemperature sensitivity is species dependent. These findings highlight the role of\nspecies composition in driving variations in temperature sensitivity along urbanization\ngradients, aligning with our remote sensing results.\nTo enhance clarity, we have revised the manuscript to strengthen the connection\nbetween the remote sensing analysis and the experimental findings, stating: “To\ninvestigate the mechanisms by which beta diversity influences the responsiveness of\nspring phenology to urbanization, we conducted species-level experiments in Beijing,\nChina.” (L. 152-154)\n[Comment 7] It is unclear why the study estimated Landsat-based SOS only from 2003\nto 2015, rather than the entire study period. Additionally, the rationale for re-estimating\nSOS with Landsat data in Beijing, as opposed to using MODIS data, is not explained.\n[Response 16] The Landsat-based SOS data for Beijing was previously available only\nfrom 2003 to 2015. In the revised version, we have extended the dataset to cover the\nentire study period (2003–2020), and found consistent results (Fig. 3).\nAs for the rationale for using Landsat instead of MODIS data, the spatial resolution of\nMODIS (500 m) is suitable for estimating city-level temperature sensitivity at a global\nscale, but it lacks the finer resolution needed for species-level analysis. To better\nunderstand the mechanisms driving our large-scale findings, we used the 30 m Landsat\ndataset for a more detailed, species-specific SOS analysis in Beijing. We have clarified\nthis in the manuscript: “To validate the results of the twig-cutting experiment, we used\n30 m resolution Landsat images to extract the SOS for the same species and conducted\nfurther species-level analysis.” (L. 442-444)\n[Comment 8] Elevation likely influences vegetation phenology, particularly in\nmountainous cities, but there are no measures to address potential uncertainties related\nto elevation effects.\n[Response 17] We agree that elevation can significantly impact the SOS, and that it\noften correlates with temperature. To address this, we conducted an additional analysis\nusing a random forest model to assess the relative importance of elevation compared to\nother factors. Our analysis revealed that urban-rural differences in preseason\ntemperature (ΔTemperature) were the primary driver of ΔSOS. In contrast, elevation\nranked as only the fourth or fifth most important variable, depending on the specific\nmethod used (Fig. S10). These findings suggest that while elevation does influence\nphenology, its impact is relatively minor compared to differences in temperature and\nspecies composition.\nFig. S10 | Drivers of urban-rural variation in the start of season (SOS). a, Variable\nimportance metrics from random forest models for the effects of urban-rural differences\nin species composition (beta diversity), species richness (ΔRichness), preseason\ntemperature (ΔTemperature), enhanced vegetation index (ΔEVI), geography\n(ΔElevation), population density (ΔPopulation density), the city-level averages of mean\nannual temperature (MAT) and precipitation (MAP), and city area on ΔS using the\nT\noutward delineation method. b, Same as a, but using the inward delineation method.\nTo further address potential uncertainties related to elevation, we also conducted an\nadditional analysis for Beijing, excluding pixels from the mountainous regions on the\ncity's western side. The results remained consistent, reinforcing the robustness of our\nfindings (see also Response 9). This analysis demonstrates that the observed patterns\nare not significantly influenced by elevation effects.\nReviewer #3\n[Comment 1] This manuscript examines the impact of urban gradients and species\ncomposition on the temperature sensitivity of spring phenology in rural and urban areas.\nThe authors use sophisticated methods to obtain a very robust data set, which showed\nthat the impact of urban gradients is a lesser driver of temperature sensitivity than\nspecies composition. While I believe this data set and results may be informative to the\nscientific community, I was let down the overall framing of the paper and the\nexplanation of how these results help explain phenological trends and/or might be\nutilized. For example, are these parameters used in Earth System Models or could this\ninformation be used to by city planners to develop communities that will be resilient to\nclimate warming? More specific comments below.\n[Response 18] We thank the reviewer for the supportive and constructive comments.\nFollowing the reviewer’s comment, we have extensively revised the discussion. Please\nsee below our responses to each comment.\n[Comment 2] L73-75. I think a clearer description of what the authors mean by\n\"stronger phenological responses\" and \"decreases in phenological responses\" is needed\nhere. Does a decreased response mean that trees will acclimate to increased\ntemperatures and their phenological timing will revert back what was observed 100\nyears ago?\n[Response 19] Thank you for your suggestion. We have clarified the meaning of\n\"phenological response to temperature\" in terms of the advancement of spring\nphenology, expressed in days per degree of warming. For example, during 1980–1994,\nthe phenological response was 4.0 ± 1.8 days/°C, meaning spring phenology advanced\nby 4.0 ± 1.8 days for each degree of warming. This response decreased to 2.3 ± 1.6\ndays/°C in the period 1999–2013, indicating a reduced sensitivity to temperature over\ntime. Similarly, \"stronger phenological responses\" in urban areas imply a greater\nadvancement of spring phenology per degree of warming compared to rural areas. We\nhave revised these sentences for clarity as follows:\n“Recent satellite-derived data suggest that urban warming reduces the sensitivity of\nspring vegetation phenology to temperature10, meaning that the start of season advances\nby fewer days per degree of warming in urban areas compared to rural ones.” (L. 50-\n53)\n[Comment 3] Introduction: Overall, I think more information about the authors'\nexpectations for changes in temperature sensitivity related to species composition is\nneeded. This is such a broad term, so maybe examples would help. Is the thought that\na greater variety of species will have a lower temperature sensitivity simply because\nthere is a high likelihood that a subset of those species will have an innate resilience to\nwarming? Or maybe that some of those species will undergo increased transpiration\nrates and cause a cooling effect in their community?\n[Response 20] We agree that the term \"species composition\" can be broad, and it is\nimportant to clarify its specific implications in our study. While species richness can\ninfluence the temperature sensitivity of the SOS, our findings indicate that urban-rural\ndifferences in beta diversity—reflecting differences in species identities rather than\nspecies counts—are the primary driver.\nIn response to the reviewer’s suggestion, we tested for potential direct effects of species\nrichness differences on the urban-rural variation in the temperature sensitivity of SOS.\nThe results showed that species richness had a minimal effect, while beta diversity\nremained the strongest driver (Fig. 2).\nTo address this point more comprehensively, we have expanded the Introduction as\nfollows: “Recent satellite-derived data suggest that urban warming reduces the\nsensitivity of spring vegetation phenology to temperature10, meaning that the start of\nseason advances by fewer days per degree of warming in urban areas compared to rural\nones. However, studies have so far provided mixed evidence for this15,16, likely due to\nthe substantial interspecific variation in phenological responses to temperature1,11,14.\nUrban-rural differences in species composition may therefore greatly impact the\nrelationship between local temperatures and tree phenology.\nUrban environments often host a wide range of species from diverse climatic and\nphylogenetic origins. Consequently, urban areas not only differ in species composition\nfrom their rural counterparts but also tend to exhibit higher species richness17,18. In\nnatural forests, greater species richness has been shown to reduce phenological\nsensitivity to temperature, either directly through the presence of temperature-resilient\nspecies or indirectly via mechanisms such as deeper root systems and altered soil\nprocesses19–21. This raises an important question: do phenological differences between\nurban and rural areas primarily reflect variations in species identities, or are they driven\nby differences in species richness?Despite the importance of this question, interspecific\nvariations in phenological sensitivity have not yet been thoroughly quantified in urban\nphenology studies11,22,23.” (L. 50-67)\nFigure 2. Drivers of urban-rural variation in temperature sensitivity of the start\nof season (ΔS ). a, Variable importance metrics from random forest models for the\nT\neffects of urban-rural differences in species composition (beta diversity), species\nrichness (ΔRichness), preseason temperature (ΔTemperature), enhanced vegetation\nindex (ΔEVI), geography (ΔElevation), population density (ΔPopulation density), the\ncity-level averages of mean annual temperature (MAT) and precipitation (MAP), and\ncity area on ΔS . b, Recursive partitioning tree for the effects of beta diversity,\nT\nΔTemperature and MAT on ΔS .\nT\n[Comment 4] In this section you state that, \"The pronounced ΔST in Asia may thus\nresult from substantial differences in species composition between urban and rural areas.\nThis is consistent with prior studies highlighting extensive landscaping efforts in\nChinese cities that increase urban tree species diversity and enhance urban-rural\ndifferences in species composition.\" This seems like an interesting finding, but I have\ntrouble understanding how to interpret this result. What does the fact that many\nAsian cities have a high beta diversity, and thus large differences in SOS and ΔST, mean for\nthe resilience of these cities to climate change? How will the response of these cities\ncompare to those with low beta diversity?\n[Response 21] Thank you for this insightful comment. While our study does not\ndirectly assess the resilience of cities to climate change, it highlights the potential\nimplications of urban-rural differences in species composition for phenological\nresponses. In response to the reviewer’s suggestion, we have expanded the\ninterpretation of these findings in the revised manuscript: “The pronounced ΔS\nT\nobserved in Asia may thus result from substantial differences in species composition\nbetween urban and rural areas. This is consistent with prior studies highlighting\nextensive landscaping efforts in Chinese cities that introduce exotic tree species and\nenhance urban-rural differences in species composition17,18. As a result, cities with high\nbeta diversity tend to exhibit more pronounced variation in phenological responses\nbetween urban and rural areas. Given that greater species diversity is associated with\nenhanced ecosystem resilience to climate warming21, cities with higher beta diversity,\nsuch as those in Asia, may have an advantage in adapting to future climate changes.\nHowever, further studies are needed to comprehensively understand how urban beta\ndiversity impacts the resilience of cities to climate change.” (L. 116-126)\n[Comment 5] Overall, I found this section to be heavy on the results and light on the\ndiscussion or interpretation of those results. Placing these results in the context of real\nworld examples or theoretical scenarios would allow the reader to see their importance.\n[Response 22] We appreciate the reviewer’s valuable feedback, which has helped us\nidentify areas where our discussion could be strengthened. In response, we have\nsubstantially expanded the discussion to place our findings in the context of practical\nscenarios and theoretical frameworks, emphasizing their broader implications. The\nrevised section reads:\n“Implication of the findings\nOur study underscores the critical role of species composition in shaping urban-rural\nvariations in the temperature sensitivity of spring phenology, with significant\nimplications for urban ecology and climate change research.\n1. Refining phenological predictions: The strong influence of species composition\nsuggests that satellite-derived phenology along urbanization gradients often fails to\naccurately reflect individual species’ responses to climate change. Integrating species-\nspecific data into phenological studies is essential for improving predictions of\necosystem responses to warming.\n2. Cascading ecological effects: Our findings suggest that community structure strongly\ninfluences plant phenology21,37, potentially leading to interspecific mismatches that can\ndisrupt trophic interactions, such as those between plants and pollinators or predators\nand prey38,39. These mismatches can cascade through ecosystems, affecting population\ndynamics and functions like carbon sequestration and nutrient cycling39,40.\nIncorporating these dynamics into vegetation models is critical for improving\npredictions of future carbon and water cycles.\n3. Urban ecology as a research framework: Urban areas are frequently used as proxies\nfor studying climate change impacts due to their elevated temperatures and CO₂\nconcentrations. However, ignoring species composition may lead to inaccurate\npredictions of ecosystem responses12. By recognizing the importance of community\nstructure, urban ecologists can better leverage cities as “natural laboratories” to explore\necosystem functioning and the ecological consequences of urbanization.\n4. Practical insights for urban planning: Strategic management of urban species\ncomposition offers actionable solutions for mitigating urban heat island effects,\nstabilizing microclimates, and reducing temperature extremes. These nature-based\napproaches enhance urban resilience to climate change and align with broader\nsustainability goals, providing cities with pathways to address ecological challenges\nwhile adapting to warming climates.” (L. 218-245)"}, {"round": 2, "reviewer_comments": {"reviewer_1": "I commend the authors for their efforts in enhancing the quality of the manuscript. Most of my previous concerns have been\nadequately addressed, particularly regarding the effect of mountain pixels and the cross-data comparison of phenology data.\nTo further improve the manuscript, I recommend the following suggestions:\n1. Site Introduction: It would be beneficial to provide a brief introduction to the sites (or pixels) used for the comparison of the\nStart of Season (SOS) between MODIS and Landsat data (Fig. S18). This will help contextualize the analysis for the\nreaders.\n2. Mountain Pixels Removal: Consider using statistical indicators to remove mountain pixels. For example, pixels that are\nhigher than the mean plus 1.5 standard deviations could be excluded. Additionally, clarify if resolutions of 100m or 500m\neffectively characterize these high-altitude pixels. This can be briefly explained or discussed.", "reviewer_2": "I appreciate the author’s efforts in addressing my concerns. The manuscript has been greatly improved, and I would like to\nrecommend the acceptance of the current version.", "reviewer_3": "I have reviewed the revised manuscript and response to previous reviewer comments. The authors were very detailed in\ntheir responses and I believe the updated manuscript is much improved and worthy of publication."}, "author_response": "Response to Reviewer #1\n[Comment 1] I commend the authors for their efforts in enhancing the quality of the\nmanuscript. Most of my previous concerns have been adequately addressed,\nparticularly regarding the effect of mountain pixels and the cross-data comparison of\nphenology data. To further improve the manuscript, I recommend the following\nsuggestions.\n[Response 1] We thank the reviewer for the positive evaluation and constructive\ncomments on our manuscript. Please see below our responses to each comment.\n[Comment 2] Site Introduction: It would be beneficial to provide a brief introduction\nto the sites (or pixels) used for the comparison of the Start of Season (SOS) between\nMODIS and Landsat data (Fig. S18). This will help contextualize the analysis for the\nreaders.\n[Response 2] Following the reviewer’s comment, we have added the vegetation SOS\nmaps based on the MODIS and Landsat datasetsin the appendix figure, and provided a\nbrief introduction to the pixels used for the comparison in the revised manuscript: “We\nextracted vegetation SOS for Beijing from MODIS and Landsat data, and matched their\nresolutions to 500 m (Fig. S17a, b). The SOS derived from Landsat was highly\ncorrelated with that from MODIS, indicating consistency between the two datasets (r =\n0.73, Fig. S17c).”. (L. 390-393)\nFig. S17 | The start of season (SOS) extracted from MODIS (a) and Landsat (b) datasets,\nand their correlation (c).\n[Comment 3] Mountain Pixels Removal: Consider using statistical indicators to\nremove mountain pixels. For example, pixels that are higher than the mean plus 1.5\nstandard deviations could be excluded. Additionally, clarify if resolutions of 100m or\n500m effectively characterize these high-altitude pixels. This can be briefly explained\nor discussed.\n[Response 3] Thank you for your valuable suggestion. Following the reviewer’s\nsuggestion, we have now excluded pixels with elevations higher than the mean ± 1.5\nstandard deviations of all pixels (268 m a.s.l.) and found consistent results (see updated\nFig. 3). This has been updated in the revised manuscript (L. 161–163; 392–394): \"To\nexclude potential elevational effects due to mountains in the west of Beijing, we\nremoved pixels with elevations higher than the mean ± 1.5 standard deviations of all\npixels (268 m a.s.l.)\".\nRegarding spatial resolution, the elevation data for Beijing had a resolution of 30 m,\nwhich can effectively characterize high-altitude pixels. This has been clarified in the\nrevised manuscript: “The elevation data for Beijing was derived from the 30 m\nresolution NASADEM dataset (https://earthdata.nasa.gov/esds/competitive-\nprograms/measures/nasadem), which was used to exclude mountain pixels.” (L. 398–\n400).\nResponse to Reviewer #2\n[Comment 1] I appreciate the author’s efforts in addressing my concerns. The\nmanuscript has been greatly improved, and I would like to recommend the acceptance\nof the current version.\n[Response 4] We thank the reviewer for recommending our manuscript for acceptance\nand for your valuable contributions to our work.\nResponse to Reviewer #3\n[Comment 1] I have reviewed the revised manuscript and response to previous\nreviewer comments. The authors were very detailed in their responses and I believe the\nupdated manuscript is much improved and worthy of publication.\n[Response 5] We thank the reviewer for recommending our manuscript for acceptance\nand for your valuable contributions to our work."}]

Decision letter