Artiﬁcial Light at Night Advances Spring Phenology in the United States

: Plant phenology is closely related to light availability as diurnal and seasonal cycles are essential environmental cues for organizing bio-ecological processes. The natural cycles of light, however, have been dramatically disrupted by artiﬁcial light at night (ALAN) due to recent urbanization. The inﬂuence on plant phenology of ALAN and its spatial variation remain largely unknown. By analyzing satellite data on ALAN intensity across the United States, here, we showed that ALAN tended to advance the start date of the growing season (SOS), although the overall response of SOS to ALAN was relatively weak compared with other potential factors (e.g., preseason temperature). The phenological impact of ALAN showed a spatially divergent pattern, whereby ALAN mainly advanced SOS at climatically moderate regions within the United States (e.g., Virginia), while its effect was insigniﬁcant or even reversed at very cold (e.g., Minnesota) and hot regions (e.g., Florida). Such a divergent pattern was mainly attributable to its high sensitivity to chilling insufﬁciency, where the advancing effect on SOS was only triggered on the premise that chilling days exceeded a certain threshold. Other mechanisms may also play a part, such as the interplay among chilling, forcing and photoperiod and the difference in species life strategies. Besides, urban areas and natural ecosystems were found to suffer from similar magnitudes of inﬂuence from ALAN, albeit with a much higher baseline ALAN intensity in urban areas. Our ﬁndings shed new light on the phenological impact of ALAN and its relation to space and other environmental cues, which is beneﬁcial to a better understanding and projection of phenology changes under a warming and urbanizing future.


Introduction
Changes in phenology, the timing of seasonally recurring biological events (e.g., leaf unfolding and senescence), have been observed in response to changes in environmental cues, including temperature [1][2][3], precipitation [4], and daylength [5]. In the context of global climate and land use/land cover change, an understanding of the drivers and mechanisms of phenology is critical for anticipating and mitigating the potential knock-on impacts on biodiversity, ecosystem functioning, and terrestrial carbon cycling [6,7].
Plant phenology is closely related to light as the diurnal and seasonal cycle of light availability is one of the crucial environmental drivers organizing ecosystem functioning [8,9]. It has been well documented that changes in natural cycle of light, measured by either daylength or photoperiod, can exert influences on plant phenology directly (e.g., advancing leaf-out dates with a longer daylength) and indirectly (e.g., regulating the sensitivity of other phenology cues) [10][11][12].
Due to recent global urbanization and advance in lighting technology, the natural cycle of light, however, has been radically altered by artificial light at night (ALAN) emitted from human settlements and infrastructure, such as street lamps, advertising billboards and buildings [13,14]. The influence of ALAN is extremely extensive. Not only does ALAN occur in the vicinity of urban areas, by means of direct irradiance and indirect light scattering (e.g., skyglow), it has also encroached deep into natural ecosystems and even remote protected areas up to hundreds of kilometers away from the suburbs, especially under cloudy conditions [15,16]. Over 23% of the Earth's land surface, particularly Europe (88%) and the United States (50%), respectively, has been reported to be under the influence of ALAN [17].
The encroachment of artificial light has profoundly disrupted natural regime of the cycles of light and darkness across space, time and wavelengths [18,19]. It interfered with information flow in flora and fauna, introducing misleading environmental signals. A recent surge in research highlights that the widespread ALAN has caused a broad array of bio-ecological consequences from molecular to ecosystem level, in such as plant phenology, organizational physiology (e.g., hormone level), life history traits (e.g., predation), activity pattern (e.g., bird migration), and community characters (e.g., diversity) [20][21][22][23].
In contrast to substantial investigations on animals, less attention has been paid to how ALAN influences plants, especially plant phenology. Thus far, the phenological impact of ALAN and its mechanisms remain largely unclear. Existing studies are limited to small sets of vulnerable species using field observations [24,25] or environment-controlled experiments (e.g., growing chamber studies) [26]. These studies appeared to be strongly influenced by research methodologies and target species, thus often leading to localized and even controversial conclusions. Moreover, most studies are constrained to population or species level, while knowledge at the ecosystem level is scarce. It hinders the understanding of the spatial pattern of ALAN's effect, as well as its underlying mechanism associated with regional climate characteristics. Given that global urban areas are continue to expand rapidly and double the current extent by the middle of the 21th century [27,28], the subsequent impacts of ALAN are projected to further increase both in magnitude and extent. An in-depth insight into the phenological impact of ALAN is thus imperative and valuable.
Here, we provide spatially explicit investigations into the influence of ALAN on spring phenology in the conterminous United States (CONUS). The overarching objectives of this study are to answer the followings questions: (1) How is spring phenology affected by ALAN and in relation to other phenology cues? (2) Does the impact on spring phenology vary spatially and, if so, what are the underlying mechanisms explaining the spatial pattern? (3) How does the ALAN's impact in urban areas differ from that in natural ecosystems?

Materials and Methods
We focused our study to the conterminous United States (CONUS) for the period of 2001-2013 ( Figure 1). The main datasets we used comprised satellite-based land surface phenology data, nighttime light images, land cover maps and climate variables. Further details are specified as follows.    MODIS Land Cover Dynamics Product (MCD12Q2, Collection 6) offers annually land surface phenology metrics with a 500-m resolution. Because of its promising quality, MCD12Q2 has been widely used to explore the phenology response to environment changes [29]. We extracted the greenup date, an indicator of the start of the growing season (SOS), of 2001-2013 from the MCD12Q2 product. We excluded the pixel with a SOS earlier than the 50th day or later than the 180th day of the corresponding year to remove the impact from outliers [30]. Given that MCD12Q2 product contains large uncertainties in densely populated urban regions, to test the robustness of our results, we extracted satellite-based SOS with two additional methods, i.e., HANTS-Maximum method [31] and Spline-Midpoint method [32] from Collection 6 MOD13A1 16-day NDVI time series. We aggregated all the resulting SOS into a 1 × 1 km grid to match the spatial resolution of other datasets used in this study (e.g., ALAN and climate data). Unless stated otherwise, datasets used in this study were acquired and exported from Google Earth Engine (GEE).

Nighttime Light Images
Nighttime light (NTL) images from the Defense Meteorological Satellite Program-Operational Linescan System (DMSP-OLS) was utilized to present ALAN intensity of CONUS during 2001-2013. DMSP-OLS data provides an annually-averaged NTL intensity with 1-km spatial resolution and a low light detection limit down to 5 × 10 −10 nWatts/ cm 2 /sr [33]. We employed methods from Elvidge et al. [34] and Liu et al. [35] to interannually calibrate images from different DMSP-OLS sensors in order to reduce the data inconsistency incurred by the lack of an on-board calibration system.
Due to the inherent sensor defects, DMSP-OLS data is affected background noise and pseudo signals by the blooming effect (i.e., where there is actually no light) [36,37]. We removed pixels with low ALAN intensity with a conservative threshold (Digital number, DN ≤ 7), therefore ensuring the remaining ALAN signals were true NTL signals rather than noises or pseudo lights. Another major drawback of DMSP-OLS data is that it only records the relative NTL intensity in DN values other than the actual NTL radiance [34]. Meanwhile, the data inconsistency cannot be completely addressed even after calibration and correction. Two additional NTL datasets were employed and applied with the same analysis in Section 2.2 to examine whether our results were biased by the drawback of DMSP-OLS data: (1) Visible Infrared Imaging Radiometer Suite (VIIRS, 2012-2018) nighttime light images from Suomi-NPP satellite. VIIRS, a new generation of NTL sensor (2012) with an on-board calibration system, provides 500-m monthly NTL images in consistent and non-saturated radiance values [38,39]. We resampled the VIIRS images into 1-km resolution and aggregated them into annually averaged NTL images from 2012 to 2018. (2) Long-term NTL time series from harmonized DMSP-OLS and VIIRS data (DMSP-VIIRS, 2001-2018) [40]. The DMSP-VIIRS data provides annually averaged NTL with 1-km spatial resolution from 2001-2018. This harmonized NTL data reduces the crosssensor difference between DMSP-OLS and VIIRS, expanding the temporal range of NTL time series. Note that the NTL intensity of the harmonized NTL data is still the relative NTL intensity.
At last, to examine whether the satellite-based NTL intensity can be used to present night sky brightness, we downloaded the artificial light sky brightness map from (https:// dataservices.gfz-potsdam.de/contact/showshort.php?id=escidoc:1541893&contactform, last access: 1 January 2021), which quantifies the artificial sky brightness based on light pollution propagation software and satellite and ground luminosity measurements [17].

Climate Data
The Daymet data is a daily climate variable product in 1-km grid. The following driving factors of phenology were extracted from Daymet: minimum air temperature, maximum air temperature, precipitation, and incident shortwave radiation. We obtained the climate zone map developed by the United States Department of Energy (DOE-developed climate zone, https://www.osti.gov/dataexplorer/biblio/dataset/1617637). DOE-developed climate zone map classifies the US into 8 climate zones (CZ1 to CZ8) at county level according to regional temperature regime, including very hot (CZ1), hot (CZ2), warm (CZ3), mixed (CZ4), cool (CZ5), cold (CZ6), very cold (CZ6), and subarctic (CZ8). CZ1 and CZ8 were not included in the analysis because they were not overlapped with our study area.

Land Cover Maps
Copernicus Global Land Cover Layers (CGLS-LC100, Collection 2), a 100-m resolution land cover map, was utilized to calculate the fraction of impervious surface, forest and cropland within each 1-km pixel. We confined our study area to forest-dominated pixels with a forest fraction larger 60% based on two considerations. First, urban areas often have a very limited vegetation fraction so setting the threshold to 60% ensured a sufficient extent of urban forest was included in the analysis. Second, most of the artificial light is received by canopy layer and forest has a more evident phenology variation [41,42]. Testing with different forest fraction thresholds, from 40% to 80%, returned similar results ( Figure S1). Cropland (fraction > 50%) was excluded because its phenology cycle is mainly manipulated by human management [43]. The 500-m land cover map from MODIS MCD12Q1 product was also obtained to examine the robustness of our analysis [44].

Analysis
We applied partial correlation analysis to unveil the impact of ALAN on SOS (R ALAN ), and its relation to preseason temperature (R T ), which is a foremost cue of SOS [45]. A negative R ALAN (R ALAN − ) indicates that increasing ALAN advances SOS, while a positive R ALAN (R ALAN + ) means increasing ALAN intensity has a delay effect on SOS. With the help of partial correlation analysis, we were able to eliminate the influence of other covariant environment cues of SOS, including preseason temperature [11], precipitation and shortwave radiation [4]. Yearly mean precipitation and shortwave radiation data were derived from Daymet product. Detailed calculation of preseason temperature was described as follows: We calculated the partial correlation coefficient (R T ) between mean air temperature (from T min and T max ) and SOS in a 1-month interval preceding SOS, e.g., 1, 2 . . . , k months (k = 1, 2, 3, . . . ) [43,45]. We limited the maximum length of preseason from November to the multi-year mean SOS value. The length (k months) that maximized the absolute value of R T was assigned as the preseason length and was used to calculate preseason temperature.
To investigate whether the divergent spatial pattern of R ALAN was related to chilling insufficiency, we calculated the chilling requirement of each pixel. Chilling requirement is usually determined as the length of chilling days (CD) when air temperature is within a certain effective range [46]. Temperature slightly higher than the freezing point, i.e., between 0 to 5 • C, is a well-acknowledged range that effectively fulfils the chilling requirement [47]. We adopted this scheme to obtain the accumulated number of CD between 1 November and the date of SOS. For a certain pixel, CD of each year was derived using the following equation (Equation (1)): where t is the t-th day from t 0 (1 November) to t SOS (date of SOS); T t is the mean temperature of day t; t 0 and t SOS indicates 1 November and the date of SOS, respectively. To quantitatively measure the relationship between CD and R ALAN , we grouped the R ALAN into 36 CD levels at a 1-day interval from 10-45 days. We applied a linear regression Remote Sens. 2021, 13, 399 5 of 16 model to capture the response of logarithmic transformed CD and R ALAN value of each CD level. Furthermore, to compare the phenological impact of ALAN in urban areas and natural ecosystems, we adopted the commonly used threshold (impervious surface fraction > 20%) to differentiate urban areas and natural ecosystems (non-urban areas), using both CGLS and MCD12Q1 land cover map products. In addition to the urban-non-urban comparison, we also examined the response of SOS to ALAN at different baseline ALAN intensity levels with 5-DN intervals.

Phenological Impact of ALAN
After excluding other covariant factors, our results overall suggest a statistically negative partial correlation (Wilcoxon signed-rank test, p < 0.001; Figure 2) between SOS and ALAN with an R ALAN value of −0.19 (median, similarly hereinafter). It indicates that the occurrence and increase of ALAN would advance satellite-based spring phenology (i.e., SOS). A high variation in R ALAN was observed (standard deviation = 0.28), ranging from −1 to 0.63 (excluding outliers). For the majority of pixels (71.7%), higher ALAN was associated with an advanced SOS with an R ALAN − of −0.30. By comparison, only 28.3% of pixels showed a delay effect on SOS, whilst the magnitude of the delay effect (R ALAN + = 0.17) was significantly weaker than that of advancing effect (p < 0.001). Preseason air temperature is well-acknowledged for its foremost control on phenology events, including SOS [48]. As expected, the overall phenological respo ALAN was weaker than to preseason temperature (RT = −0.60). The advancing effe warming temperature was observed in the bulk of pixels (91.3%) and had a much association with SOS (RT − = −0.61) than RALAN, whereas their counterparts, i.e., RT + a LAN + did not differ greatly in correlation magnitude (RT + = 0.14). This finding suggest albeit with an overall advancing effect, ALAN is a minor cue of SOS and its capabi influencing spring phenology is relatively weaker than temperature.
Our findings are largely consistent with previous studies based on field observ [24,25,49] and manipulated experiments [18,50,51], which highlight the associatio tween ALAN and SOS but also a great diversity among species and sites. For examp leveraging field phenology observations and ALAN intensity measured by DMSP data, it was found a 7.5 days advance in budburst date of tree species in response unit of ALAN increase [49]. By comparison, a large variation, from four days earlie Preseason air temperature is well-acknowledged for its foremost control on spring phenology events, including SOS [48]. As expected, the overall phenological response to ALAN was weaker than to preseason temperature (R T = −0.60). The advancing effect of a warming temperature was observed in the bulk of pixels (91.3%) and had a much closer association with SOS (R T − = −0.61) than R ALAN , whereas their counterparts, i.e., R T + and R ALAN + did not differ greatly in correlation magnitude (R T + = 0.14). This finding suggests that, albeit with an overall advancing effect, ALAN is a minor cue of SOS and its capability in influencing spring phenology is relatively weaker than temperature. Our findings are largely consistent with previous studies based on field observations [24,25,49] and manipulated experiments [18,50,51], which highlight the association between ALAN and SOS but also a great diversity among species and sites. For example, by leveraging field phenology observations and ALAN intensity measured by DMSP-OLS data, it was found a 7.5 days advance in budburst date of tree species in response to per unit of ALAN increase [49]. By comparison, a large variation, from four days earlier to 12 days later, was observed in the flowering date of grass species exposed to LED lighting in a semi-natural environment [18]. The differences between our findings and others can be resulted from the differences in research methodologies, such as the data used (e.g., satellitederived phenology vs. field observation-based phenology), research objects (e.g., forest vs. grass species), and research settings (e.g., real-world setting vs. environment-controlled condition).
Potential bias can be incurred by DMSP-OLS data for its representativeness of the real night sky brightness. What night light satellite detects is mainly the zenith-directed light emissions at cloud-free nights [52], such that the light scattered into the night sky (e.g., skyglow) can be underestimated. To examine this issue, we compared DMSP-OLS data with a night sky luminance map calibrated by light pollution propagation model and satellite and ground measurements [17]. A moderate correlation between DMSP-OLS and night sky luminance was found (R 2 = 0.63, p < 0.001; Figure 3), suggesting the appropriateness of using DMSP-OLS to present ALAN intensity.  Besides, systematic errors of DMSP-OLS data and the uncertainties in phenology data extracted from MCQ12Q2 and the intercalibration process of DMSP-OLS data may also influence the validity our results. We thus applied the same analysis to two additional ALAN maps from DMSP-VIIRS data (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) and VIIRS data (2012-2018) and two additional phenology estimates using HANTS-maximum and spline-midpoint methods. As demonstrated in Figure 4, results derived from four comparative groups were generally consistent to DMSP-OLS based result, underpinning the robustness of our results. Among these, using the phenology data estimated by HANTS-Maximum method yielded the highest RALAN value (−0.22), while replacing DMSP-OLS with VIIRS data for ALAN mapping produced the highest variation (standard deviation = 0.37).
It should be noted that only DMSP-OLS was used for the subsequently analysis, instead of DMSP-VIIRS data or VIIRS data, because of three concerns: (1) Data availability. Subjected to the availability of NTL data and the corresponding phenology and climate data, VIIRS data were only available for 6.12 ± 1.17 yr. for each pixel, while DMSP-OLS data is available for 12.2 ± 3.28 yr. The partial correlation analysis was vulnerable to the impact of short-term observations [53]. (2) Overpass time. DMSP-OLS (7:30 pm) has a more ideal overpass time than VIIRS (1:30 Besides, systematic errors of DMSP-OLS data and the uncertainties in phenology data extracted from MCQ12Q2 and the intercalibration process of DMSP-OLS data may also influence the validity our results. We thus applied the same analysis to two additional ALAN maps from DMSP-VIIRS data (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) and VIIRS data (2012-2018) and two additional phenology estimates using HANTS-maximum and spline-midpoint methods. As demonstrated in Figure 4, results derived from four comparative groups were generally consistent to DMSP-OLS based result, underpinning the robustness of our results. Among these, using the phenology data estimated by HANTS-Maximum method yielded the highest R ALAN value (−0.22), while replacing DMSP-OLS with VIIRS data for ALAN mapping produced the highest variation (standard deviation = 0.37). Remote Sens. 2021, 13, x FOR PEER REVIEW 8 of 17

Divergent Spatial Pattern of the Phenological Impact of ALAN
Satellite data products allow us to look into the spatial pattern of the phenological impact of ALAN. Figure 5a demonstrates the spatial distribution of RALAN over the CO-NUS. Since most of the forest and human settlements are located in the east part of the CONUS, our results are thus mainly distributed in the entire Northeast, as well as the east part of the Midwest and South. Comparing RALAN + and RALAN − pixels, we found a noticeably clustered spatial pattern. This pattern was corroborated by global Moran's I test (Moran's I index = 0.43, p < 0.001), a spatial autocorrelation measurement, indicating the distribution of RALAN + and RALAN − pixels was more spatially clustered other than randomly scattered ( Figure 5b). As illustrated by latitudinal transects, most of RALAN − pixels were clustered between 30°N and 45°N, while RALAN + pixels mainly spread over low (<30°N) and high (>45°N) latitude regions (Figure 6b). Similar results were found in the four comparative analysis with different input datasets (Figure 5c-f). It should be noted that only DMSP-OLS was used for the subsequently analysis, instead of DMSP-VIIRS data or VIIRS data, because of three concerns: (1) Data availability. Subjected to the availability of NTL data and the corresponding phenology and climate data, VIIRS data were only available for 6.12 ± 1.17 yr. for each pixel, while DMSP-OLS data is available for 12.2 ± 3.28 yr. The partial correlation analysis was vulnerable to the impact of short-term observations [53]. (2) Overpass time. DMSP-OLS (7:30 p.m.) has a more ideal overpass time than VI-IRS (1:30 a.m.) because the peak lighting hour is approximately from sunset to 22:00 p.m. [54]. Using DMSP-OLS data enables us to evaluate the phenological response to ALAN when most artificial lights are on. (3) Despite a longer temporal coverage of DMSP-VIIRS data, the uncertainties of crosssensor calibration algorithm to generate DMSP-VIIRS data may propagate to and confound our analysis [55].

Divergent Spatial Pattern of the Phenological Impact of ALAN
Satellite data products allow us to look into the spatial pattern of the phenological impact of ALAN. Figure 5a demonstrates the spatial distribution of R ALAN over the CONUS. Since most of the forest and human settlements are located in the east part of the CONUS, our results are thus mainly distributed in the entire Northeast, as well as the east part of the Midwest and South. Comparing R ALAN + and R ALAN − pixels, we found a noticeably clustered spatial pattern. This pattern was corroborated by global Moran's I test (Moran's I index = 0.43, p < 0.001), a spatial autocorrelation measurement, indicating the distribution of R ALAN + and R ALAN − pixels was more spatially clustered other than randomly scattered ( Figure 5b). As illustrated by latitudinal transects, most of R ALAN − pixels were clustered between 30 • N and 45 • N, while R ALAN + pixels mainly spread over low (<30 • N) and high (>45 • N) latitude regions (Figure 6b). Similar results were found in the four comparative analysis with different input datasets (Figure 5c-f).
Furthermore, we grouped R ALAN into climate zones developed by the United States Department of Energy (DOE-developed climate zones), which categorize the United States into 8 climate regimes, from the hottest climate zone 1 (CZ1) to the coldest climate zone 8 (CZ8), according to regional temperature characters [56]. It further demonstrated the spatial pattern of ALAN's impact on SOS and revealed its linkage to regional temperature gradient (Figure 6a)   Furthermore, we grouped RALAN into climate zones developed by the United States Department of Energy (DOE-developed climate zones), which categorize the United  Furthermore, we grouped RALAN into climate zones developed by the United States Department of Energy (DOE-developed climate zones), which categorize the United States into 8 climate regimes, from the hottest climate zone 1 (CZ1) to the coldest climate

Urban. vs. Non-Urban. Area
In addition to the regional spatial pattern, we probed into how ALAN's effect varied over urban areas and non-urban natural ecosystems. We classified our study area into urban areas and nature ecosystem using MODIS and CGLS land cover map products, separately, and compared the corresponding R ALAN values. Considering that ALAN is more intensive in urban areas, we expected SOS had a stronger response to ALAN in urban areas than natural ecosystems. Surprisingly, only slight R ALAN differences were observed between urban areas (MODIS: −0.180; CGLS: −0.199) and nature ecosystems (MODIS: −0.174; CGLS: −0.173) (Figure 7a). We also found that R ALAN values remained at a similar level regardless of the baseline ALAN magnitude (Figure 7c). It suggests that, despite a relatively low magnitude of the ALAN encroaching into natural ecosystems, it is enough to trigger phenology changes similar to the effect in brightly-lit urban areas. This phenomenon may be because light is mainly served as an information source, rather than an energy resource, for photoperiodism related processes [57]. Organisms use the changes in the natural light cycles as indicators of time and place and to detect changes in environment [58]. Hence, even though ALAN is at low intensities and/or of short duration in natural ecosystems, it is still able to incur marked phenological effect [42,57]. mately equal amount of RALAN + (51.0%) and RALAN − (49.0%) pixels. For very cold regions (CZ7), like New Hampshire, Vermont, Minnesota, and eastern Wisconsin, the impact of ALAN was even inversed to cause a delayed SOS (p < 0.001), whilst the number of RALAN + pixels (63.3%) surpassed RALAN − pixels (36.7%). These above-mentioned phenomena underlined a spatially divergent response of SOS to ALAN between climatically moderate regions (CZ3 to CZ6) and hot (CZ2)/very cold regions (CZ7).

Urban vs. Non-Urban Area
In addition to the regional spatial pattern, we probed into how ALAN's effect varied over urban areas and non-urban natural ecosystems. We classified our study area into urban areas and nature ecosystem using MODIS and CGLS land cover map products, separately, and compared the corresponding RALAN values. Considering that ALAN is more intensive in urban areas, we expected SOS had a stronger response to ALAN in urban areas than natural ecosystems. Surprisingly, only slight RALAN differences were observed between urban areas (MODIS: −0.180; CGLS: −0.199) and nature ecosystems (MODIS: −0.174; CGLS: −0.173) (Figure 7a). We also found that RALAN values remained at a similar level regardless of the baseline ALAN magnitude (Figure 7c). It suggests that, despite a relatively low magnitude of the ALAN encroaching into natural ecosystems, it is enough to trigger phenology changes similar to the effect in brightly-lit urban areas. This phenomenon may be because light is mainly served as an information source, rather than an energy resource, for photoperiodism related processes [57]. Organisms use the changes in the natural light cycles as indicators of time and place and to detect changes in environment [58]. Hence, even though ALAN is at low intensities and/or of short duration in natural ecosystems, it is still able to incur marked phenological effect [42,57].

Mechanisms of the Divergent Spatial Pattern
Our study found a spatially divergent pattern in the phenological response of ALAN between climatically moderate and hot/very cold regions, respectively. To our knowledge, this study is for the first time to reveal such a spatial pattern, whilst there is barely any wellestablished theory to explain the possible causes of it. Here, we proposed three potential and non-exclusive explanations for the phenomenon.

The Interplay among Chilling, Forcing and Photoperiod
We first attribute the advancing effect of ALAN in climatically moderate regions to the interplay among the intensity and length of winter temperature (chilling), warm spring temperature (forcing), and photoperiod (daylength). Photoperiod exerts an impact on phenology by influencing dormancy release through its sophisticated interaction with other environmental cues [46]. Studies utilizing in situ observations and environmentcontrolled experiments have found that insufficient winter chilling (or spring forcing) can be compensated, at least partly, by a longer photoperiod and additional forcing (or chilling) [10,11]. As a result, this interplay mechanism between photoperiod and the other two environmental cues, i.e., chilling and forcing, leads to an earlier SOS when a longer photoperiod occurs [59]. The ALAN lights up the night sky after sunset, and thus it could be regarded as a factor prolonging the original photoperiod (daylength) [18]. Note that in this study what we measured was how SOS responded to the change in ALAN intensity rather than the change of the duration that plants perceived ALAN. Due to data unavailability, it is currently impossible to measure the duration that the plants are exposed to ALAN during the night at a large scale [41]. We assume that if a pixel has a higher ALAN intensity, it is likely to indicate a longer overall photoperiod of the plants therein. The rationale behind this assumption is that within a 1 × 1 km pixel, artificial light from a brighter source can spread and penetrate to a larger extent [60]. Therefore, a higher proportion of plants within the pixel will be exposed to ALAN and experience a prolonged photoperiod. Based on our assumption and the aforementioned interplay mechanism, we ascribe the advancing effect to the prolonged photoperiod induced by ALAN increase.

Chilling Insufficient
Second, explanation for ALAN's delay effect in hot and very cold regions could be rather complicated. As mentioned above, a longer photoperiod and additional forcing compensate for insufficient chilling, resulting in an earlier SOS. Nevertheless, the interplay mechanism has been reported to be non-linear [53]-that is, excessive chilling insufficiency may outweigh the compensating capacity of a longer photoperiod and additional forcing, consequently constraining the advancing effect on SOS and even turning it into the opposite [9,30,61]. To test whether this "chilling insufficiency" theory held true for the relationship between ALAN and SOS, we calculated the number of effective chilling days (CD, 0 < T < 5 • C) of each pixel and grouped them into DOE-developed climate zones. Our result demonstrated that hot region (CZ2, 20.5 ± 4.3 days) and very cold region (CZ7, 21.7 ± 3.1 days) had fewer number of chilling days than climatically moderate regions (CZ3 to CZ6, 30.7 ± 5.1 days, p < 0.001; Figure 8, red bar). It implies that hot and very cold regions are inclined to suffer more from chilling insufficiency than climatically moderate regions. The insufficient chilling can be caused by the following reasons: in lower temperate and subtropical regions like CZ2, global warming has resulted in warm winters in recent years, making it even hard to meet the chilling requirement [53,62]. For very cold regions (CZ7), the winter temperature is far below the effective range for fulfilling the chilling requirement (0 < T < 5 • C) [47] and the low temperature may also pose extra stress to the ALAN-SOS relationship (Figure 8, blue bar) [63].
Furthermore, we calculated the R ALAN of each CD group at a 1-CD interval and explored the relationship between CD and R ALAN by linear regression analysis (Figure 9). The logarithmic transformed CD (ln(CD)) and R ALAN was found negatively related with a slope of −0.35 (R 2 = 0.82, p < 0.001). For most pixels whose CD value exceeds 19.3 days, the sign of R ALAN remained negative, while R ALAN changed into a positive sign when CD was below this threshold. This phenomenon suggests that the advancing effect of ALAN only functions on the premise that chilling requirement is met, at least to a certain degree. When CD is below a certain threshold, the chilling insufficiency will offset the ALAN-induced SOS advancement and lead to an insignificant or reversed ALAN-SOS relationship. Comparing the sign of R T (R T + and R T − ) for pixels showing a delay effect of ALAN (R ALAN + ), it turned out that a higher temperature could still result in an earlier SOS for more than 98% of R ALAN + pixels. Again, it shows that ALAN is not such a primary phenology driver as preseason temperature and its impact on SOS is sensitive to chilling requirement, while the impact of temperature is more tolerant against chilling insufficiency.
These findings supported our hypothesis that the delay effect of ALAN could be incurred by insufficiency of winter chilling.
Remote Sens. 2021, 13, x FOR PEER REVIEW 12 of 17 (CZ7), the winter temperature is far below the effective range for fulfilling the chilling requirement (0 < T < 5 °C) [47] and the low temperature may also pose extra stress to the ALAN-SOS relationship (Figure 8, blue bar) [63]. Furthermore, we calculated the RALAN of each CD group at a 1-CD interval and explored the relationship between CD and RALAN by linear regression analysis (Figure 9). The logarithmic transformed CD (ln(CD)) and RALAN was found negatively related with a slope of −0.35 (R 2 = 0.82, p < 0.001). For most pixels whose CD value exceeds 19.3 days, the sign of RALAN remained negative, while RALAN changed into a positive sign when CD was below this threshold. This phenomenon suggests that the advancing effect of ALAN only functions on the premise that chilling requirement is met, at least to a certain degree. When CD is below a certain threshold, the chilling insufficiency will offset the ALANinduced SOS advancement and lead to an insignificant or reversed ALAN-SOS relationship. Comparing the sign of RT (RT + and RT − ) for pixels showing a delay effect of ALAN (RALAN + ), it turned out that a higher temperature could still result in an earlier SOS for more than 98% of RALAN + pixels. Again, it shows that ALAN is not such a primary phenology driver as preseason temperature and its impact on SOS is sensitive to chilling requirement, while the impact of temperature is more tolerant against chilling insufficiency. These findings supported our hypothesis that the delay effect of ALAN could be incurred by insufficiency of winter chilling.  Furthermore, we calculated the RALAN of each CD group at a 1-CD interval and explored the relationship between CD and RALAN by linear regression analysis (Figure 9). The logarithmic transformed CD (ln(CD)) and RALAN was found negatively related with a slope of −0.35 (R 2 = 0.82, p < 0.001). For most pixels whose CD value exceeds 19.3 days, the sign of RALAN remained negative, while RALAN changed into a positive sign when CD was below this threshold. This phenomenon suggests that the advancing effect of ALAN only functions on the premise that chilling requirement is met, at least to a certain degree. When CD is below a certain threshold, the chilling insufficiency will offset the ALANinduced SOS advancement and lead to an insignificant or reversed ALAN-SOS relationship. Comparing the sign of RT (RT + and RT − ) for pixels showing a delay effect of ALAN (RALAN + ), it turned out that a higher temperature could still result in an earlier SOS for more than 98% of RALAN + pixels. Again, it shows that ALAN is not such a primary phenology driver as preseason temperature and its impact on SOS is sensitive to chilling requirement, while the impact of temperature is more tolerant against chilling insufficiency. These findings supported our hypothesis that the delay effect of ALAN could be incurred by insufficiency of winter chilling.

Differences in Species' Life Strategies
The difference in species life strategies can also play a role in inducing the divergent response to ALAN. Studies have revealed that early successional species, such as hazel, poplars and birch, adopt a riskier life strategy to make the best use of growing resources and to maximize the growing season [7,64]. The start of the growing season of these species only relies on one or two key phenology drivers, i.e., forcing and/or chilling, and is less sensitive to other minor cues like ALAN and photoperiod. However, late successional species have evolved to be more conservative against environment changes to hinder early growth in unfavorable conditions (e.g., warm winter with high frost risk). These species have a strict requirement for environment cues and therefore can be both temperature-dependent and ALAN-dependent (photoperiod-dependent) [10,59]. In addition, tree cutting experiments reported that the leaf-out date of species from high latitudes is independent of photoperiod [5]. Spring phenology events are mainly driven by chilling and forcing rather than photoperiod in regions with long winter (like regions of CZ7 in this study). On the contrary, species originated from climatically moderate regions (like CZ3 to CZ6) may rely on additional information from photoperiodism [65].

Implications
Our study provides empirical evidence about how spring phenology is noticeably influenced by ALAN encroachment. We showed that ALAN tended to advance SOS, and this effect was statistically significant but relatively weaker than the effect caused by preseason temperature. We found a spatially divergent pattern of the effect in climatically moderate regions and hot/very cold regions and attributed it to: (1) the interplay among forcing, chilling and photoperiod; (2) chilling insufficiency, as we found the phenological impact of ALAN was sensitive to the accumulated chilling days; (3) the differences in life strategies of species. The results have an implication to understand the interaction mechanism between phenology shifts and environmental stimuli. It helps to better forecast how phenology will be influenced by and feedback to the warming climate and the environment being modified by human activities. Revealing the spatial variation of the impact of ALAN and its linkage to other environmental cues of phenology is a complement to in situ observation and experiment-based studies. Our results imply that previous studies about daylength and photoperiod may underestimate its actual potential in shifting phenology, especially for places where ALAN is prevailing. Hence, the long-lasting controversy on how ALAN (or photoperiod) affects phenology can be a result of not taking ALAN into consideration (for studies related to photoperiod) and/or the divergent spatial pattern of the ALAN-SOS relationship (for studies related to photoperiod and ALAN). We foresee a promising potential in improving the existing temperature-photoperiod based phenology models by incorporating the impact of ALAN [66].
Furthermore, our result showed that urban areas and natural ecosystems were subject to similar magnitude of phenological impact from ALAN, although the ALAN intensity in natural ecosystems was much lower than that in urban areas. This finding calls for further awareness of the phenological and ecological impact of ALAN in natural ecosystems because most of the areas under the influence of ALAN are natural ecosystems rather than human-dominated regions. In our case, about 87.7% of areas affected by ALAN were natural land, whereas urban areas only accounted for a small percentage (12.3%) of ALAN-affected areas (Figure 7b). On a global scale, human settlement only occupies 3% of the Earth's land surface, but the resulting influences caused by ALAN will be rather far-reaching beyond urban footprint, where about 20% of the Earth's surface under ALAN's influence is natural land. Such widespread influence of ALAN is projected to continue to exacerbate in terms of both magnitude (1.8% yr −1 ) and extent (2.2% yr −1 ) [52]. Hence, it is imperative to propel lighting regulation and a better lighting design to reduce unnecessary light emission in urban areas and to minimize light pollution encroachment in natural ecosystems.

Limitations
The following issues should be noted and require further scrutiny in future studies. First, our results are based on coarse resolution satellite images, for which the mix pixel issue is inevitable. The phenology data extracted from long-term field observations can be further used. It can also help to explore the potential impact of species and taxa on the ALAN-phenology interaction. However, field phenology observation data were not considered for this study due to their limited temporal availability. For example, the species-site combination of the USA National Phenology Network data (USA-NPN; 2009), on average, only has 2.01 ± 1.87 yr. of valid observations. Less than 5% of the combinations have more than 6-year observations. Using short-term observations will bias our results by making the partial correlation analysis sensitive to the start and end values [53]. In future studies, it would be beneficial to incorporate long-term satellite and field phenology observations, together with other possible phenology cues (e.g., urban heat island effects), to have a more elaborated investigation into the phenology changes induced by human activities.
Second, it is challenging thus far to measure the long-term sky brightness changes at a large spatial scale. Even though a good correlation is found between ALAN intensity derived from DMSP-OLS and night sky brightness, satellite-based ALAN intensity inevitably underestimates what is actually perceived by plants. Likewise, at this moment, there are barely any data that can provide temporally exhaustive information about the duration of ALAN throughout the night. Due to this limitation, we are only able to analyze the association between ALAN and SOS based on the assumption that a higher ALAN intensity will lead to a longer ALAN duration. These issues could be addressed when more data sets become available, e.g., night light satellite with a short revisit time or a geostationary orbit [67,68].
Third, in our analysis chilling days were calculated as the number of days with temperature between 0 and 5 • C, which is a commonly-adopted approach [47,69]. A recent study using in-situ phenology records of 30 perennials [70], however, found that using the temperature between 0 and 5 • C was an invalid method to calculate chilling days, because it could not reflect the negative relationship between chilling days and heat requirement. Instead, temperature below 5 • C is a more suitable way. It is an issue worthwhile further investigation and validation since it is closely linked with the explanation of causes of the spatial divergent response of SOS to ALAN.

Conclusions
The natural cycle of light, as one of the key components organizing the ecosystem, has been dramatically disrupted by the sharp increase and expansion of artificial light at night due to global urbanization. This study leveraged nighttime light remote sensing images and satellite-based phenology data to explore the phenological responses of ALAN, reveal its relation to space and other environmental cues, and explain the underlying mechanisms. We shed spatially explicit insights on how the start date of the growing season (SOS) of plant was affected by ALAN, and found that ALAN generally tended to advance SOS (R ALAN = −0.19). We revealed that the effect had a spatially divergent pattern between climatically moderate regions and very hot/cold regions. Further analysis showed the divergent pattern was mainly attributable to the high sensitivity to chilling insufficiency, while the interplay among chilling, forcing and photoperiod and the difference in species life strategies also played a role. Our findings not only evoke further awareness of the phenological and ecological impacts of artificial light at night, but help us to better understand how phenology will be influenced by and feedback to the continuously changing climate and the environment being modified by human activities.