Vegetation Phenology Influenced by Rapid Urbanization of The Yangtze Delta Region

Impacts of urbanization and climate change on ecosystems are widely studied, but these drivers of change are often difficult to isolate from each other and interactions are complicated. Ecosystem responses to each of these drivers are perhaps most clearly seen in phenology changes due to global climate change (warming climate) and urbanization (heat island effect). The phenology of vegetation can influence many important ecological processes, including primary production, evapotranspiration, and plant fitness. Therefore, evaluating the interacting effects of urbanization and climate change on vegetation phenology has the potential to provide information about the long-term impact of global change. Using remotely sensed time series of vegetation on the Yangtze River Delta in China, this study evaluated the impacts of rapid urbanization and climate change on vegetation phenology along an urban to rural gradient over time. Phenology markers were extracted annually from an 18-year time series by fitting the asymmetric Gaussian function model. Thermal remote sensing acquired at daytime and nighttime was used to explore the relationship between land surface temperature and vegetation phenology. On average, the spring phenology marker was 9.6 days earlier and the autumn marker was 6.63 days later in urban areas compared with rural areas. The spring phenology of urban areas advanced and the autumn phenology delayed over time. Across space and time, warmer spring daytime and nighttime land surface temperatures were related to earlier spring, while autumn daytime and nighttime land surface temperatures were related to later autumn phenology. These results suggest that urbanization, through surface warming, compounds the effect of climate change on vegetation phenology.


Introduction
Vegetation phenology refers to annual reoccurring cycles of plant activity, such as leaf onset and offset in spring and autumn, respectively. These cycles are driven by climatic cues [1] and are therefore an important indicator of trends in annual weather conditions [2,3]. For example, in Europe it was found that a spring air temperature increase of 1 • C has been associated with an advance in the beginning of the growing season by 7 days [4]. Similar results are found in the eastern United States with spring arriving 4-6 days earlier since the mid-1960s [5][6][7][8][9][10]. In China, spring has advanced on average 2.88 days per decade in response to spring warming [1]. Therefore, trends in vegetation

Data
The data used in this research includes information on the spatial and temporal extent of urban land cover, spring and autumn vegetation phenology, and land surface temperature. The urban areas were identified from the Suomi National Polar-orbiting Partnership (SNPP) Visible Infrared Imaging Radiometer Suite (VIIRS) nighttime light data ( Figure 2a) by applying an optimal thresholds method based on the urban land-use product generated by the Resource and Environment Data Center of Chinese Academy of Sciences [18,29,30]. Strong spatial heterogeneity in the spectral characteristics of urban areas, including many mixed pixels at moderate resolution, increase uncertainty in the location of the urban boundary. Nighttime light data, however, produce reliable delineation of urban clusters, and clearly distinguish urban areas from the surrounding suburban and rural areas that lack nighttime lighting [6]. In this study, the urban areas derived from nightlight data from 2000,2005,2010, and 2015 were used to represent the urban area in 2001-2002, 2003-2007, 2008-2012, and 2013-2018, respectively [25].
The vegetation types data (30 m resolution; Figure 2b were acquired from the work of Gong et al. 2012, and included five main categories: mixed forest (MF), evergreen forest (EF), shrubland (SB), grassland (GS), and cropland (CP) [31,32]. Land cover was resampled to 250 m to match the MODIS enhanced vegetation index (EVI) data. The daytime LST at 10:30 AM and nighttime LST at 22:30 PM local time were acquired from the MODIS LST product (MOD11A2, 8-day composite) with a spatial resolution of 1 × 1 km. It has been confirmed that the accuracy of MODIS LST reaches 1K [33].

Data
The data used in this research includes information on the spatial and temporal extent of urban land cover, spring and autumn vegetation phenology, and land surface temperature. The urban areas were identified from the Suomi National Polar-orbiting Partnership (SNPP) Visible Infrared Imaging Radiometer Suite (VIIRS) nighttime light data ( Figure 2a) by applying an optimal thresholds method based on the urban land-use product generated by the Resource and Environment Data Center of Chinese Academy of Sciences [18,29,30]. Strong spatial heterogeneity in the spectral characteristics of urban areas, including many mixed pixels at moderate resolution, increase uncertainty in the location of the urban boundary. Nighttime light data, however, produce reliable delineation of urban clusters, and clearly distinguish urban areas from the surrounding suburban and rural areas that lack nighttime lighting [6]. In this study, the urban areas derived from nightlight data from 2000, 2005, 2010, and 2015 were used to represent the urban area in 2001-2002, 2003-2007, 2008-2012, and 2013-2018, respectively [25].
The vegetation types data (30 m resolution; Figure 2b were acquired from the work of Gong et al. 2012, and included five main categories: mixed forest (MF), evergreen forest (EF), shrubland (SB), grassland (GS), and cropland (CP) [31,32]. Land cover was resampled to 250 m to match the MODIS enhanced vegetation index (EVI) data. The daytime LST at 10:30 a.m. and nighttime LST at 22:30 p.m. local time were acquired from the MODIS LST product (MOD11A2, 8-day composite) with a spatial resolution of 1 × 1 km. It has been confirmed that the accuracy of MODIS LST reaches 1K [33].
Vegetation phenology information was extracted from Moderate Resolution Imaging Spectrometer (MODIS) enhanced vegetation index (EVI) product (MOD13Q1 EVI, 16-day composite, 250 m spatial resolution) for the period 2001-2018. These data reduce the impact of variation in atmospheric scattering compared with the Normalized Difference Vegetation Index (NDVI) [3,34], and have been demonstrated to be suitable for monitoring phenology information in urban areas [2,14,19,20,35,36] Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 18 in the year 2015.
Vegetation phenology information was extracted from Moderate Resolution Imaging Spectrometer (MODIS) enhanced vegetation index (EVI) product (MOD13Q1 EVI, 16-day composite, 250m spatial resolution) for the period 2001-2018. These data reduce the impact of variation in atmospheric scattering compared with the Normalized Difference Vegetation Index (NDVI) [3,34], and have been demonstrated to be suitable for monitoring phenology information in urban areas [2,14,19,20,35,36]

Methods
The methods used to analyze the impacts of urbanization on the vegetation phenology are given in Figure 3 to give a clear representation. First, vegetation phenology markers were derived from MODIS EVI time series data, and urban areas as well as buffer zones were extracted from the nighttime light data. Then, we analyzed the spatial-temporal variation of vegetation phenology in urban, suburban, and rural areas from 2001 to 2018. Finally, the relationship between the urban-rural phenology difference and LST difference was constructed using linear regression analysis, and the vegetation phenology responses to urbanization across vegetation types were analyzed.

Methods
The methods used to analyze the impacts of urbanization on the vegetation phenology are given in Figure 3 to give a clear representation. First, vegetation phenology markers were derived from MODIS EVI time series data, and urban areas as well as buffer zones were extracted from the nighttime light data. Then, we analyzed the spatial-temporal variation of vegetation phenology in urban, suburban, and rural areas from 2001 to 2018. Finally, the relationship between the urban-rural phenology difference and LST difference was constructed using linear regression analysis, and the vegetation phenology responses to urbanization across vegetation types were analyzed.

Asymmetric Gaussians Smoothing
An asymmetric Gaussian function (A-G) was fitted to the MODIS EVI time series data to retrieve phenology markers for spring and autumn (Equation (1)) [10,37]. To fit EVI data accurately, MODIS pixel quality information was used to assess pixel reliability of EVI data. Pixels with higher reliability (higher quality) were given greater weight in the fitting algorithm (following Table 1), which reduced the impact of data noise [38,39]. The asymmetric Gaussian function applied was:

Asymmetric Gaussians Smoothing
An asymmetric Gaussian function (A-G) was fitted to the MODIS EVI time series data to retrieve phenology markers for spring and autumn (Equation (1)) [10,37]. To fit EVI data accurately, MODIS pixel quality information was used to assess pixel reliability of EVI data. Pixels with higher reliability (higher quality) were given greater weight in the fitting algorithm (following Table 1), which reduced the impact of data noise [38,39]. The asymmetric Gaussian function applied was: where x 1 is the Day of Year (DOY) of the maximum or minimum over time, t. The fitted parameters x 2 and x 3 determine the width and flatness of curve when t > x 1 . Similarly, x 4 and x 5 determine the width and flatness of the curve when t < x 1 .

Extraction of Phenology Metrics
Three phenology markers, start of season (SOS), end of season (EOS), and growth season length (GSL), were used to analyze the influence of urbanization on vegetation phenology from MODIS EVI data. SOS is the day of year (DOY) when the fitted curve rises beyond a threshold, defined relative to the season amplitude. EOS is the DOY when the declining side of the fitted curve falls below a threshold, also defined relative to the seasonal amplitude. Compared with rural areas, the mean EVI of urban areas is low, so a fixed threshold, that is not sensitive to seasonal amplitude, cannot be applied to both urban and surrounding areas. We adopted a threshold equal to 0.2 of the seasonal amplitude, which has been shown to be suitable for extracting urban phenology [10,16,26,40] (Figure 4). Pixels with minimal mean vegetation cover (EVI), such as water and entirely hardened surfaces, were removed. These unrealistic extreme values outside of the range of 50 to 180 for SOS and the range of 240 to 330 for EOS were removed from the analysis [26,41]. The GSL is defined as duration between SOS and EOS in days.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 18 to the season amplitude. EOS is the DOY when the declining side of the fitted curve falls below a threshold, also defined relative to the seasonal amplitude. Compared with rural areas, the mean EVI of urban areas is low, so a fixed threshold, that is not sensitive to seasonal amplitude, cannot be applied to both urban and surrounding areas. We adopted a threshold equal to 0.2 of the seasonal amplitude, which has been shown to be suitable for extracting urban phenology [10,16,26,40] ( Figure  4). Pixels with minimal mean vegetation cover (EVI), such as water and entirely hardened surfaces, were removed. These unrealistic extreme values outside of the range of 50 to 180 for SOS and the range of 240 to 330 for EOS were removed from the analysis [26,41]. The GSL is defined as duration between SOS and EOS in days.

Urban Buffer Zones and Phenology Change
We used a range of distances from urban areas to define zones of urbanization (i.e., urban, suburban, and rural), and then analyzed these areas to understand the difference in phenology

Urban Buffer Zones and Phenology Change
We used a range of distances from urban areas to define zones of urbanization (i.e., urban, suburban, and rural), and then analyzed these areas to understand the difference in phenology response across urbanization levels [23]. It was previously reported [15] that the mean distance from urban land cover affecting vegetation phenology is less than 20 km from the urban boundary. Therefore, we adopted an approach that defined distance from the urban edge within 20 km as the suburban area. Phenology in the rural area (20-25 km from urban edge) was adopted as a reference.
To explore the changes of vegetation phenology along the urbanization gradient, suburban areas were further divided into six buffer areas (0-1, 1-2, 2-5, 5-10, 10-15, 15-20 km) ( Figure 5) extending from the urban edge [42]. Response of phenology to urbanization was defined as the phenology differences between urban zones (i.e., urban, suburban) and rural areas. We assumed that the farther away from the urban edge, the lower the level of urbanization. Phenology differences were therefore defined as follows: where ∆P is the difference between phenology markers extracted from (sub)urban buffer zones (P ub ) and the rural reference (P r ). A negative value of ∆P for SOS represents an earlier SOS in the urban areas compared with rural, while a positive value of ∆P for EOS represents a later EOS in the urban areas. The earlier SOS and later EOS correspond to a longer GSL.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 18 where ΔP is the difference among phenology markers between urban or buffer zones and rural areas, and ΔLST is the LST difference between urban or buffer zones and rural areas. Parameter m is the slope of the model. Where m >0, later SOS or EOS (whichever is used to calculate ∆P) are related to higher LST; if m <0, earlier SOS or EOS is related to higher temperatures. Parameter n represents the intercept, signifying the value of ∆P when there is no temperature difference between urban and rural areas.
To provide greater insight into the nature of phenology change over time and between urban and rural areas, we performed the previously described analyses on individual land cover types and on all land covers combined.

Phenology Changes with Time
Spatial patterns in the phenology markers were complex, but even in map form earlier SOS and later EOS were evident in urban areas relative to rural ( Figure 6). These spatial differences were more noticeable earlier in the timeseries, motivating an analysis of urban and rural phenology over time. The Yangtze River Delta has experienced significant phenology changes during 2001-2018, especially in suburban and rural areas. Urban areas exhibited earlier SOS and later EOS and, as a result, the GSL in urban areas was the longest of any area, in any given year ( Figure 7). Conversely, both the latest SOS and earliest EOS were found in rural areas, resulting in the shortest GSL. The SOS in urban areas advanced (became earlier in the year) from 2001-2018 (slope = -0.43 days/year, p < 0.05; Figure 7a). A similar pattern was observed in suburban and rural areas at the rate of −0.85 days/year (p < 0.05) and -0.72 days/year (p < 0.05), respectively. The EOS delayed significantly in suburban (0.52 days/year, p < 0.01) and rural areas (0.64 days/year, p < 0.01). However, there was no trend in the EOS of urban areas, which were more stable and generally later in the year than observed in suburban and rural areas (Figure 7b). The GSL increased in all areas from 2001 to 2018 (Figure 7c). The GSL for suburban

Urbanization Effect on Phenology
The relationships between LST and phenology were estimated by calculating Spearman's correlations between the phenology differences (∆P) and LST differences (∆LST) across the urban and six buffer zones toward rural zones. For those higher correlation coefficients, the linear regression model between phenology difference and LST difference was constructed according to Equation (3).
Remote Sens. 2020, 12, 1783 where ∆P is the difference among phenology markers between urban or buffer zones and rural areas, and ∆LST is the LST difference between urban or buffer zones and rural areas. Parameter m is the slope of the model. Where m > 0, later SOS or EOS (whichever is used to calculate ∆P) are related to higher LST; if m < 0, earlier SOS or EOS is related to higher temperatures. Parameter n represents the intercept, signifying the value of ∆P when there is no temperature difference between urban and rural areas.
To provide greater insight into the nature of phenology change over time and between urban and rural areas, we performed the previously described analyses on individual land cover types and on all land covers combined.

Phenology Changes with Time
Spatial patterns in the phenology markers were complex, but even in map form earlier SOS and later EOS were evident in urban areas relative to rural ( Figure 6). These spatial differences were more noticeable earlier in the timeseries, motivating an analysis of urban and rural phenology over time. The Yangtze River Delta has experienced significant phenology changes during 2001-2018, especially in suburban and rural areas. Urban areas exhibited earlier SOS and later EOS and, as a result, the GSL in urban areas was the longest of any area, in any given year ( Figure 7). Conversely, both the latest SOS and earliest EOS were found in rural areas, resulting in the shortest GSL. The SOS in urban areas advanced (became earlier in the year) from 2001-2018 (slope = −0.43 days/year, p < 0.05; Figure 7a). A similar pattern was observed in suburban and rural areas at the rate of −0.85 days/year (p < 0.05) and −0.72 days/year (p < 0.05), respectively. The EOS delayed significantly in suburban (0.52 days/year, p < 0.01) and rural areas (0.64 days/year, p < 0.01). However, there was no trend in the EOS of urban areas, which were more stable and generally later in the year than observed in suburban and rural areas (Figure 7b). The GSL increased in all areas from 2001 to 2018 (Figure 7c). The GSL for suburban and rural areas increased at the rate of 1.37 days/year (p < 0.01) and 1.36 days/year (p < 0.01), respectively, which was a faster rate than that observed in urban areas (0.4 days/year, p < 0.05).
Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 18 and rural areas increased at the rate of 1.37 days/year (p < 0.01) and 1.36 days/year (p < 0.01), respectively, which was a faster rate than that observed in urban areas (0.4 days/year, p < 0.05).

Phenology Changes along the Urbanization Gradient
Phenology differences (∆P) exhibited heterogeneities along an urban-rural gradient across years in the Yangtze River Delta (Figure 8). The SOS was earlier in urban compared with rural areas in 14 out of 18 years, ranging from −22.3 days in 2006 to 0.3 days in 2003. Conversely, the EOS was later in urban zones (16 out of 18 years) with the differences ranging from 15 days to -0.9 days. On average, SOS advanced 9.6 days and EOS delayed 6.63 days across the urban-rural gradient. This is consistent with previous reports that the differences for SOS and EOS in urban centers relative to rural zones

Phenology Changes along the Urbanization Gradient
Phenology differences (∆P) exhibited heterogeneities along an urban-rural gradient across years in the Yangtze River Delta (Figure 8). The SOS was earlier in urban compared with rural areas in 14 out of 18 years, ranging from −22.3 days in 2006 to 0.3 days in 2003. Conversely, the EOS was later in urban zones (16 out of 18 years) with the differences ranging from 15 days to -0.9 days. On average, SOS advanced 9.6 days and EOS delayed 6.63 days across the urban-rural gradient. This is consistent with previous reports that the differences for SOS and EOS in urban centers relative to rural zones were 11.9 and 5.4 days, respectively [26]. Thus, SOS is more sensitive to the urbanization gradient, compared to EOS [11]. In addition, we also found that the maximum value of EOS tended to occur in the buffered zones adjacent to urban areas, rather than within urban areas.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 18 were 11.9 and 5.4 days, respectively [26]. Thus, SOS is more sensitive to the urbanization gradient, compared to EOS [11]. In addition, we also found that the maximum value of EOS tended to occur in the buffered zones adjacent to urban areas, rather than within urban areas.

Phenology Response to Urbanization Across Vegetation Types
To explore the variations in phenology response to urbanization by vegetation types, phenology metrics were compared between urban zones (i.e. urban, suburban) and rural areas for each vegetation type over the period 2013 to 2018. Phenology differences of different vegetation types from 2013 to 2018 are given in Figure 9. ∆SOS1, ∆EOS1, and ∆GSL1 refer to the differences in SOS, EOS, and GSL, respectively, between urban and rural areas. ∆SOS2, ∆EOS2, and ∆GSL2 are defined as the differences in SOS, EOS, and GSL, respectively, between suburban and rural areas.
The SOS for different vegetation types was found to have distinct differences between urban, suburban, and rural zones (p < 0.05), and was earlier in urban areas by more than 10 days compared with rural areas (Figure 9). Croplands in suburban areas exhibited the earliest SOS difference relative to rural areas (−6.8 days, p < 0.01), followed by the evergreen forests (−4.3 days, p < 0.01) and shrublands (−4.0 days, p < 0.01). On average, across all land cover types, the EOS in suburban areas was later than the EOS observed in urban areas. This is consistent with the previous finding ( Figure 8). Only croplands showed significant differences in EOS between urban and rural areas (3.6 days, p < 0.05). However, the EOS for all vegetation types in suburban areas was greater than that in rural areas (p < 0.05), ranging from 1.1 to 4.8 days.
Urbanization had the greatest impact on the GSL of croplands, with mean differences of 13.9 days (urban) and 10.4 days (suburban) (p < 0.01; Table 2). This is consistent with a previous report, which found the GSL difference in croplands between urban and rural areas was 15.2 days [6]. Comparison of GSL in urban and rural areas, it was found that grassland had lowest GSL increase (8.8 days), followed by mixed forest, evergreen forest and shrubland types (11.4, 11.6, 12.1 days). However, between suburban and rural areas, the GSL difference of grassland, mixed forest, evergreen forest and shrubland types were 4.1, 6.8, 7.2, 8.8 days, respectively. areas (p < 0.05), ranging from 1.1 to 4.8 days.
Urbanization had the greatest impact on the GSL of croplands, with mean differences of 13.9 days (urban) and 10.4 days (suburban) (p < 0.01; Table 2). This is consistent with a previous report, which found the GSL difference in croplands between urban and rural areas was 15.2 days [6]. Comparison of GSL in urban and rural areas, it was found that grassland had lowest GSL increase (8.8 days), followed by mixed forest, evergreen forest and shrubland types (11.4、11.6、12.1 days). However, between suburban and rural areas, the GSL difference of grassland, mixed forest, evergreen forest and shrubland types were 4.1, 6.8、7.2、8.8 days, respectively.   Figure 9.

Correlations between the Urban-Rural Phenology Difference and LST Difference
Our study reveals the change and trend of phenology in the Yangtze River Delta from the perspective of time and space, which can be attributed to the influence of urbanization. According to Spearman's correlations (Table 3), along urban to rural gradients, phenology markers are correlated with LST ( Figure 10). Positive LST anomalies were observed between urban and suburban buffer areas, and rural areas, and these anomalies were closely related to earlier SOS and later EOS. Via linear regression, we observed that ∆SOS and ∆LST in spring was negatively correlated (p < 0.01) and ∆EOS and ∆LST in autumn were positively correlated (p < 0.05). A 1°C increase in daytime LST led to a 4.3 day advance in SOS, and a 1°C increase in nighttime LST led to a 12.1 day advance in SOS. In autumn, a positive correlation between ∆EOS and ∆LST was observed, which means the increase of LST leads to the delay of EOS. A 1°C increase in daytime LST led to a 3.6 day delay of EOS, and a 1°C increase in nighttime LST led to a 5.7 day delay of SOS.

Temporal Profile of Vegetation Phenology
The impacts of global warming and urbanization on vegetation phenology have been gradually recognized and accepted as a consensus. In our study area, advanced SOS, delayed EOS, and prolonged GSL were observed in urban, suburban, and rural areas, which was expected based on previous research [11,16,24,26,43,44]. However, our work uncovered some specific details about the sensitivity of phenology to climate change, urbanization, and the specific role of LST in each of the land cover types studied. Compared with suburban and rural areas, urban areas exhibited the weakest SOS trends over time, suggesting climate change has the least impact in areas that are already fully developed. The suburban areas have the greatest rate of change (−0.85 days/year), reflecting the combined effects of climate warming and increased urban development. In rural areas, where the effect of urbanization is minimized, we also found a strong effect of time on spring phenology (−0.72 days/year) and growing season length (1.36 days/year). Autumn phenology (EOS) changed most rapidly in suburban (0.52 days/year) and rural (0.64 days/year) areas, and least rapidly in urban areas (−0.03 days/year). Indeed, autumn phenology of urban areas appears to have been relatively stable over this period and the temporal trends in growing season length are driven entirely by trends in SOS in urban areas. Over time, there was a strong convergence in spring and autumn timing across the three land cover types. This could be due to either continued urbanization that makes land cover in suburban and rural landscapes function more like urban land cover over time, or a gradual weakening of the effect of climate change on phenology as summer and winter temperatures pass plant physiological thresholds [45].

Impacts of Urbanization on Vegetation Phenology
While it is difficult to definitively separate the effects of urbanization and climate change on phenology in this rapidly developing region, temporal changes in the spatial pattern of phenology along the urban-rural gradient show that phenology changed the most in suburban areas closest to the urban core where the pressure of development is greatest (Figures 6 and 7). Our conceptual model is that the effect of urbanization will be weakest in urban and rural areas where the rate of development is slowest, and strongest in suburban areas that developed most rapidly. Conversely, the dominate driver of temporal change in urban and rural areas might be climate change, because these areas were either already built out (urban) or remained relatively distant from urban land cover (rural). The pattern of delayed SOS, earlier EOS, and shorter GSL with increasing distance from urban areas was sometimes close to linear, but there was generally some evidence of a curve linear response, the shape of which changed over time (Figure 8). Early in the time series (2001), the SOS and EOS converged to the rural value 15 km from the urban edge. This resulted in a concave downward (SOS) or upward (EOS) pattern (Figure 8, year 2001). However, over time, the curve-linear pattern changed, with SOS becoming concave downward and EOS becoming concave upward. The best explanation of this pattern is continued urbanization in the mid-range distances from the urban edge driving more rapid changes in phenology in these suburban areas. These spatio-temporal patterns are a window into the interacting effects of urbanization and climate change, and highlight the rapid changes experienced by vegetation in suburban areas over this period.
There has been considerable work that uses remote sensing to describe the effect of the urban heat island on phenology, with most studies finding that vegetation phenology varies significantly along the urbanization gradient [11,16,19,24,26,28,44,46,47]. One study quantified the effect of urbanization as occurring within 6 km of the urban edge [24]. In another study, the effect of urbanization was measured from the urban center and detected out to 32 km [28]. However, remote sensing time series are made longer each year and many previous studies did not attempt to examine patterns over time (i.e., were limited to an urban-rural comparison). Leveraging the temporal dimension, we observed that in 2001 the effect of urbanization extended nearly 15 km from the urban edge. However, by 2018 the effect extended to 25 km from the urban edge, and these effects were seen most strongly in the EOS timing (Figures 6 and 7).

Phenology of Different Vegetation Types
We found that vegetation type affects the impact of urban and suburban land cover on phenology, which is consistent with the previous findings. Similar to the effect of land cover on temporal trends (Figure 7), the effect of vegetation type was strongest for EOS. While SOS timing exhibited the expected trend towards earlier springs with increasing urbanization for all vegetation types (Figure 8), EOS was similar for urban and suburban areas in mixed forest, evergreen forest, shrubland, and grassland. Only cropland showed the expected earlier EOS in suburban areas relative to urban. This is consistent with the previous finding of White et al [13]. Croplands occupy a large proportion in many cities, and have been subject to urban expansion over the past several decades [17,48]. From this work, it is clear that remaining croplands in urban and suburban landscapes enjoy a longer growing season. These trends might be expected to lead to multiple cropping cycles, and associated changes in fertilization rates, crop types, and other practices.
Responses of different types of vegetation to climate warming and urbanization have been noted in previous work. A study of the phenology of leaves and male cones of Chinese pine found that the rural-urban gradient drove spatial variation in phenology in the range of 0.36-1.92 (leaves) and 0.41-0.66 days/km (cones) [46]. Similar results were found in plant species in the Seoul Capital Area, specifically, that spring budburst sensitivity is greatest in Prumus mume, while the most dominant autumn leaf coloring sensitivity was observed in Acer palmatum [11]. Limited by the spatial resolution of MODIS data used here, it is difficult to separate the effect of warming and urbanization on different plant species in this study. However, the differences we did detect among vegetation types are consistent with past findings.

Relations between LST and LSP
Urbanization is expected to influence LST through the removal of vegetation and replacement with hardened surfaces with lower evaporative cooling and higher heat capacity. It has been found that the urban heat island intensity contributed greatly to the increase of urbanization effects on SOS [25] and UHI effects exist in many cities [24]. Therefore, we considered LST a critical variable to analyze in the context of phenology, similar to past work [24,49,50]. We focused our analysis on the SOS, EOS, and LST differences between urban, suburban, and rural areas ( Figure 10). In the spring, we found a negative correlation between ∆SOS and ∆LST, and in the autumn, a positive correlation between ∆EOS and ∆LST. However, both ∆SOS and ∆EOS were least sensitive to ∆LST day (compared with ∆LST night ), and the effect of ∆LST day clearly was reduced close to urban areas. Only at a distance of greater than 5 km did ∆SOS become correlated with ∆LST day , largely due to a lack of any substantial gradient in ∆LST day over this distance from urban areas. In contrast, ∆SOS and ∆EOS were more sensitive to ∆LST night , and this sensitivity continued to the urban area. The dominance of the effect of ∆LST night on phenology is therefore pervasive in our results. This suggests that vegetation is more sensitive to nighttime temperature, which would be expected if hardened surfaces were keeping air temperatures warmer through the night. These findings were comparable to past results [11], in which it was found that in spring a 1°C warming at night resulted in an earlier spring budburst of 2.76 days, and in the autumn, a 1°C warming at night was related to a later autumn leaf coloring of 2.55 days.
Prolonged GSL may lead to cooler surfaces on average due to the fact that vegetated surfaces are substantially cooler than non-vegetated surfaces, which would represent a negative feedback process on phenological changes related to climate warming. It was found that the quality and character of vegetation patches in Boston USA affect local heat island intensity [51]. However, at the scale of the MODIS pixel, the extent and temporal variation of the impact of prolonged GSL on the LST is difficult to characterize.

Conclusions
The impacts of urbanization on vegetation phenology were observed over space and time in the rapidly developing region of the Yangtze delta. As urbanization progressed across the region, land surface temperatures increased. This change was greater in more rural areas, and changes in land surface temperatures (primarily nighttime temperatures) drove a change in spring and autumn phenology. The continued urbanization is enlarging the urban heat island, and making once rural areas more similar to urban areas in terms of their vegetation phenology. An earlier SOS, a later EOS, and a longer GSL are associated with the increase in LST along the urban-rural gradient direction. Different vegetation types demonstrated different phenology characteristics, but the differences between vegetation types was not as stark as the differences between urban, suburban, and rural areas. By comparing the urban-rural gradient over time, we attempted to separate the effects of urbanization from climate change. While this challenge proved difficult, we observed interesting patterns representing the spread of the UHI from urban areas. This increasing UHI effect drove greater changes in suburban and rural landscapes than in urban areas. There is still work to be done to understand how these observations might inform the management of urban landscapes, possibly through planting different types of vegetation, or adapting agricultural practices to the longer growing seasons experienced in urban areas. In general, higher spatial resolution data will be necessary for most of these applications.