More artificial light delays leaf aging in city trees during autumn

Distinct spatial patterns were observed in both DFS and ALAN intensity across cities. DFS occurred earlier at higher latitudes, while ALAN intensity was greater in developed regions such as North America, Europe, South Korea, and Japan (Supplementary Fig. 1a, c). Temporally, DFS was delayed in 343 cities (75.9%) from 2001 to 2022, with an average trend of 0.77 days per year (Supplementary Fig. 1b). Concurrently, ALAN increased in 430 of the 452 cities (95.1%), with a mean rate of 0.53 nW cm^–2 sr^–1 y^–1 (Supplementary Fig. 1d).

Through an extensive review and synthesis of relevant literature, we proposed two hypotheses to investigate the mechanisms by which increases in ALAN intensity could delay DFS (Supplementary Table): (H1) increased ALAN enhances vegetation photosynthetic carbon fixation, thereby postponing DFS; (H2) ALAN alters the response of DFS to climatic drivers, further affecting its temporal variation. 2

We further used DMT and DMT_ALAN to project future DFS under two climatic scenarios: Shared Socioeconomic Pathways (SSP) 245 and 585 (Methods). Compared to the original DMT, DMT_ALAN projected consistently later DFS dates. During 2081–2100, DMT_ALAN estimated delays of 2.46 ± 9.73 days under SSP245 and 2.3 ± 12.05 days under SSP585 relative to the DMT ensemble mean. Moreover, DMT_ALAN projections indicated significant delaying trends in DFS over 2001–2100, with slopes of 0.147 and 0.224 days per year under SSP245 and SSP585, respectively (P < 0.001; Fig. 4e).

This study integrated in-situ records of DFS from two authoritative networks of phenological observations: the China Phenological Observation Network (CPON)³⁹ and the Pan European Phenology Project (PEP725)⁴⁰. We identified and removed potential outliers using median absolute deviation (MAD) method, which is more robust to skewed distributions and extreme values than standard deviation-based approaches. Specifically, MAD for a site- and species-specific DFS data set (DFS₁, DFS₂, …, DFS_i) was calculated as:1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{MAD}}={\mathrm{median}}\left(|{\mathrm{DFS}}_{i}-{\mathrm{median}}({\mathrm{DFS}})|\right)$$\end{document}MAD=median∣DFSi−median(DFS)∣

Any DFS value that exceeded 2.5 times the MAD was considered an outlier and excluded from further analysis³⁴. We focused on DFS records from 2001 onward. We excluded site-species DFS time series with fewer than 10 years of data to ensure a robust analysis based on sufficient observations. We thereby obtained 5113 DFS records representing 228 species at 10 sites across China from CPON (2001–2018) and 57,881 DFS records representing six species at 1209 sites across Europe from PEP725 (2001–2015). See Supplementary Fig. 2 for details of the distribution of the observation sites.

This study also used the Land Cover Dynamics phenological product (MCD12Q2 Version 6.1) of the Moderate Resolution Imaging Spectroradiometer (MODIS)⁴¹. These phenological metrics are derived from time series of the two-band enhanced vegetation index (EVI2). Spring phenology is defined as the first date when EVI2 exceeds 15% of its amplitude, and autumnal phenology is defined as the last date when EVI2 falls below 15% of its amplitude. These data offer crucial information for analyzing and monitoring vegetation phenological changes. We extracted the Greenup and Dormancy layers from the MCD12Q2 data set for 2001–2022, which represent the leaf-unfolding date (LUD) and DFS, respectively.

Satellite-based DFS estimates rely on seasonal changes in vegetation greenness (i.e. EVI2), which might not perfectly match ground-based measurements because of differences in spatial scale and pixel-mixing impacts. Specifically, MODIS EVI2 mixed pixels often incorporate multiple land cover types (e.g. forest, cropland, bare soil), while ground observations reflect localized homogeneous vegetation conditions. To address these potential biases, we analyzed the satellite and in situ data sets independently rather than through direct integration or comparison⁸.

We used the annual ALAN data derived from the “NPP-VIIRS-like” nighttime light product⁴², which was generated by harmonizing observations from two satellite sensors, the Defense Meteorological Satellite Program Operational Linescan System (DMSP-OLS) and Suomi National Polar-orbiting Partnership Visible Infrared Imaging Radiometer Suite (NPP-VIIRS), using an improved autoencoder-based neural-network model. The NPP-VIIRS-like dataset exhibits strong ability to capture both population density distributions and changes in nighttime illumination across multiple spatial and temporal scales, closely resembling the composited NPP-VIIRS NTL data⁴². In addition, we also employed the original DMSP-OLS⁴³ and NPP-VIIRS nighttime light datasets⁴⁴, as well as another long-term harmonized nighttime light dataset (H-NTL-v2)⁴⁵. Unlike the NPP-VIIRS-like dataset, the H-NTL-v2 dataset was generated by calibrating nighttime light data through the conversion of NPP-VIIRS observations into DMSP-OLS-like data. These complementary datasets were used to assess the consistency of DFS responses to ALAN across different data sources.

To better characterize the seasonal variations of ALAN intensity and its impact on DFS, we determined the monthly ALAN considering the potential duration of ALAN exposure⁴⁶. We first calculated the solar declination-hour angle based on the latitude and day of the year (DOY) for each pixel to calculate the theoretical daily night length under cloud-free conditions. We then calculated the monthly average theoretical night length for each pixel as the duration exposed to ALAN. Finally, the monthly ALAN was calculated as:2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\rm{A}}{\rm{LAN}}_{\rm{month}}=\frac{{\rm{NH}}_{\rm{month}} \times {12}}{{\sum }_{{\rm{month}}={1}}^{12} {\rm{NH}}_{\rm{month}}} \times {\rm{ALAN}}_{\rm{year}}$$\end{document}ALANmonth=NHmonth×12∑month=112NHmonth×ALANyearwhere ALAN_month represents the monthly ALAN value, NH_month represents the average nightly hours per month and ALAN_year represents the yearly ALAN value.

We obtained monthly data for temperature, precipitation, and shortwave radiation for from the TerraClimate data set⁴⁷. For DFS modeling, we used six-hourly temperature data from the CRU JRA v2.4 data set⁴⁸ and converted it into daily temperature data. To predict DFS under future scenarios, we also obtained daily temperature from the CMIP6 NorESM2-MM data set under both SSP245 and 585 scenarios⁴⁹. We used the Global Artificial Impervious Area product to delineate urban boundaries⁵⁰. Cities with built-up areas >100 km² were selected²⁰. Ultimately, 452 major cities in the Northern Hemisphere were identified. Buffer zones covering 100, 200, and 400% of the urban areas were generated to analyze the distance-decay effect of ALAN influence (Supplementary Fig. 4).

We conducted a spatial attribution analysis of the ALAN-DFS relationship by incorporating multiple factors of economic development and vegetation, including the Human Development Index, gross domestic product, per capita GDP, water-use efficiency, light-use efficiency, aboveground biomass, tree density, and canopy height. The data for HDI, GDP, and per capita GDP were obtained from the gridded global data sets for GDP and HDI⁵¹. For future scenarios, changes in per capita GDP relative to current statistical values were set based on country-specific scenarios from the llASA-SSP database⁵². We also used vegetation data, including aboveground biomass⁵³, tree density⁵⁴, and canopy height⁵⁵. The WUE and LUE data represented the average monthly values from 2001 to 2022. Monthly WUE and LUE were calculated as:3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\rm{WUE}}=\,\frac{{\rm{GPP}}}{{\rm{ET}}}$$\end{document}WUE=GPPET4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\rm{LUE}}=\,\frac{{\rm{GPP}}}{{\rm{SW}}\times 0.45\times {\rm{FPAR}}}$$\end{document}LUE=GPPSW×0.45×FPARwhere GPP is gross primary productivit⁵⁶, ET is evapotranspiration⁵⁷, and FPAR is the fraction of photosynthetically active radiation⁵⁸, they are all derived from the MODIS products. SW is shortwave radiation derived from the TerraClimate monthly data set.

For the temporal analysis, we used the time-series data for V_cmax, and SIF as mediators to link ALAN and DFS. V_cmax data from were obtained from the data set of maximum rate of carboxylation (V_cmax25)⁵⁹ and the SIF data were sourced from the global ‘OCO-2’ SIF data set⁶⁰. See Supplementary Table 1 for detailed information on data descriptions, spatial and temporal resolutions, temporal coverage, and data sources for all data sets used in this study.

We applied the Theil-Sen trend analysis to quantify temporal changes in DFS, ALAN, and climate events. As a robust and nonparametric method, the Theil-Sen estimator is resistant to outliers and measurement errors, making it well-suited for long-term time-series analysis. To assess the statistical significance of these trends, we employed the Mann-Kendall test at a 0.05 significance level. This test is widely used because it does not assume normality or linearity and is unaffected by missing data or outliers.

The relationship between DFS and specific influencing factors in this study involved both simple correlation and partial-correlation analyses. We defined the influence of a factor on DFS as the maximum correlation coefficient between the phenological date and the variable values of different preseason lengths during the study period. We determined this influence by first calculating the optimal preseason length for each variable relative to DFS. Specifically, we used an exhaustive approach: starting from the month of the multi-year average phenological date, we extended backward in one-month steps up to six months, calculating the mean of each variable for these preseason periods⁶¹. We then performed a partial-correlation analysis between these means and the phenological time series to identify the preseason length that yielded the highest partial-correlation coefficient. Influencing factors often interact, so a partial-correlation analysis is typically used to isolate the effect of a single factor on DFS by controlling for the influence of other variables. Before performing correlation analyses, linear detrending is usually applied to each variable individually to remove long-term temporal trends, allowing the study to focus on the impact of the interannual variability of influencing factors on the interannual changes in DFS.

We employed a moving average approach to analyze the correlation between ALAN and climate sensitivity. Specifically, during the period from 2001 to 2022, we used an 11-year moving window with a 1-year step to calculate the average values of ALAN and climate sensitivity within each window, resulting in 12 data sets. Based on these moving averages, we then calculated the correlation coefficient to evaluate the relationship between ALAN and climate sensitivity. To avoid potential multicollinearity among the influencing factors, we employed ridge regression. Ridge regression is a linear regression model with L2 regularization, which introduces a penalty term in the loss function to suppress large regression coefficients, thereby effectively alleviating the problem of collinearity among independent variables. This method enhances the stability and generalization ability of the model. In this study, the response variable is DFS, and the predictor variables include preseason ALAN and climatic factors. We used the normalized anomalies of climatic factors, ALAN, and DFS as the inputs for the regression model. The resulting regression coefficients are interpreted as sensitivities of drivers.

To assess the potential causal influence of ALAN on the DFS and to further investigate whether ALAN, the photosynthesis indicators (SIF, V_cmax), and DFS exhibit not only statistical associations but also potential causal relationships, we employed the PCMCI+ algorithm to determine the direction of causality. Compared to traditional methods such as Granger causality that focus on bivariate prediction-based relationships, PCMCI+ can handle high-dimensional multivariate time series and account for complex interdependencies among multiple variables through conditional independence testing. PCMCI consists of two key steps: the PC algorithm and the Momentary Conditional Independence (MCI) test^62,63, which aimed at addressing the common issue of autocorrelation in time series data. In this study, we used the extended PCMCI+ framework, which is capable of identifying both lagged and contemporaneous causal relationships. Prior to applying PCMCI + , we first conducted partial correlation analysis to assess the statistical associations among ALAN, the photosynthesis indicators, and DFS. Subsequently, PCMCI+ was applied to infer whether these associations reflected potential causal structures while accounting for interdependencies among variables. It should be noted that PCMCI+ identifies statistically inferred causal relationships based on conditional independence, which provides evidence for potential causal mechanisms but does not establish definitive causation. All analyses were performed using the Tigramite Python package.

We used interpretable machine learning with Shapley Additive Explanations (SHAP) to identify the key drivers of the spatial distribution of the effects of ALAN, quantified as the sensitivity of DFS to ALAN⁶⁴. Various factors were categorized into three groups: (1) the economic-development factors HDI, GDP, per capita GDP, and ALAN; (2) the climatic factors temperature, precipitation, and shortwave radiation; and (3) the vegetation factors WUE, LUE, aboveground biomass, tree density, and canopy height. Both averages and trends were calculated for multi-year ALAN and the climatic variables. Detailed descriptions of all variables are provided in Supplementary Table 1.

We used these factors as predictor variables to develop an eXtreme Gradient Boosting (XGBoost) model. XGBoost is an optimized gradient boosting decision tree algorithm that is widely used in regression, classification, and ranking tasks due to its high computational efficiency and excellent predictive performance⁶⁵. XGBoost effectively handles the complex nonlinear relationships and interactions commonly found in spatial environmental data and is highly robust to multicollinearity and missing values. In addition, its built-in regularization mechanism helps prevent overfitting, which is particularly important for modeling highly heterogeneous spatial data. Combined with efficient parallel computing capabilities, XGBoost can quickly and reliably process large-scale spatial datasets, enhancing the accuracy and interpretability of attribution analysis. To better interpret the predictions made by XGBoost, we employed the SHAP method to quantify the marginal contribution of each predictor variable to the target variable. SHAP, a unified interpretation framework based on the game-theoretic Shapley value, fairly allocates the contribution of each feature to the model’s predictions. It not only explains individual predictions but also reveals the overall importance of variables, greatly improving the transparency and trustworthiness of complex “black-box” models. These methods were implemented in Python using the “scikit-learn”, “xgboost”, and “shap” packages⁶⁵.

We improved four widely used DFS models, i.e., CDD⁶⁶, DM⁶⁷, SIAM⁶⁸, and DMT^69,70. Their respective versions enhanced with ALAN are termed CDD_ALAN, DM_ALAN, SIAM_ALAN, and DMT_ALAN.

The CDD model estimated DFS based on cumulative cold temperatures, as such conditions are likely key environmental signals triggering leaf coloration or senescence:5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\rm{CDD}}}_{t}=\mathop{\sum }\limits_{t={t}_{0}}^{{t}_{\gamma }}max \left({T}_{{\rm{th}}}-{T}_{(t)},0\right)$$\end{document}CDDt=∑t=t0tγmaxTth−T(t),06\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\rm{DFS}}={t}_{r},{\mathrm{if}},{{\rm{CDD}}}_{t}\ge {{\rm{CDD}}}_{{\rm{th}}}$$\end{document}DFS=tr,if,CDDt≥CDDthwhere CDD_t denotes the cumulative cooling degree days calculated from the start date t₀ to a given day t, while T(t) represents the average daily temperature on day t. The calculation begins on day t₀, defined as the first day when the temperature drops below the critical threshold (T_th), and continues until day t_y, when the total CDD_t reaches a predefined value (CDD_th). The day t_y is then recorded as the predicted DFS, represented by DOY_ty. In our study, the T_th was set within a range of 0 to 50 °C.

Building upon the CDD model, the DM approach integrates both temperature and photoperiod as key factors influencing DFS:7\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${S}_{{\rm{sen}}}(t)={S}_{{\rm{sen}}}(t-1)+{R}_{{\rm{sen}}}(t)$$\end{document}Ssen(t)=Ssen(t−1)+Rsen(t)8\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${R}_{{\rm{sen}}}(t)=\left\{\begin{array}{c}{\left[{T}_{b}-{T}_{(t)}\right]}^{x}\times f{\left[{P}_{(t)}\right]}^{y},{\mathrm{if}}t\ge {{\rm{DOY}}}_{s}\\ 0,{\mathrm{if}}t < {{\rm{DOY}}}_{s}\end{array}\right.$$\end{document}Rsen(t)=Tb−T(t)x×fP(t)y,ift≥DOYs0,ift<DOYs9\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f\left[{P}_{(t)}\right]=\frac{{P}_{(t)}}{{P}_{s}}{\mathrm{or}}\left[{P}_{(t)}\right]=1-\frac{{P}_{(t)}}{{P}_{s}}$$\end{document}fP(t)=P(t)PsorP(t)=1−P(t)Ps10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\rm{DFS}}={t}_{y},{\mathrm{if}}\,{S}_{{\rm{sen}}}(t)\ge {Y}_{{\rm{crit}}}$$\end{document}DFS=ty,ifSsen(t)≥Ycritwhere S_sen(t) represents the state of leaf coloring on day t, while R_sen(t) denotes its rate. T_(t) is the average temperature for day t, with T_b as the base temperature. P_(t) refers to the P_(t) on day t, and P_s is the critical threshold for photoperiod. The parameters x and y range between 0 and 2. Leaf coloration was initiated on DOY_s when the daily temperature dropped to or below T_b and the photoperiod shortened to P_s. The process continued until day t_y, when S_sen(t) accumulated to a predefined threshold (Y_crit). The day t_y was then identified as the DFS.

SIAM and DMT are modified from DM, with Y_crit linearly related to the LUD anomaly in SIAM, and to spring–summer temperature in DMT:11\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Y}_{\mathrm{crit}}=a+b\times {{\rm{LUD}}}_{{\rm{a}}}$$\end{document}Ycrit=a+b×LUDa12\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Y}_{\mathrm{crit}}=a+b\times \mathrm{Tss}$$\end{document}Ycrit=a+b×Tsswhere LUD_a is the LUD anomaly and Tss is the mean spring-summer temperature.

We constructed a model coefficient, ALANin_i, which included the annual ALAN values, and incorporated ALANin_i into the calculation of the forcing rate in the DFS model:13\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\rm{ALANin}}}_{i}\,={\rm{e}}}^{k\times (\frac{{{\rm{ALAN}}}_{i}\,-{{\rm{ALAN}}}_{\max }}{{{\rm{ALAN}}}_{\max }})}$$\end{document}ALANini=ek×(ALANi−ALANmaxALANmax)where ALAN_i represents the ALAN value in year i and ALAN_max represents the maximum ALAN value.

To improve the predictive accuracy of the model, we employed the Particle Swarm Optimization algorithm to optimize key model parameters (such as T_th, T_b, P_s, Y_crit, a, b). This algorithm simulates the movement and collaborative search of particles within the parameter space to find the optimal combination of parameters that minimizes model error in a multi-dimensional space. This algorithm is known for its fast convergence and strong global search capability, making it particularly suitable for nonlinear and multi-parameter coupled problems⁷¹. The resulting optimal parameter set was then used for subsequent model training and validation.

The metrics for evaluating the accuracy between the predicted and observed DFS included the R, the proportion of observations with significant at P < 0.05, RMSE and AIC. AIC is used to balance the model’s goodness of fit and complexity. A smaller AIC value indicates that the model achieves a good fit while maintaining lower complexity, making it more concise and efficient⁷².14\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\rm{AIC}}=n\times {\mathrm{ln}}\left(\frac{\sum _{i=1}^{n}{\left({p}_{i}-{o}_{i}\right)}^{2}}{n}\right)+2k$$\end{document}AIC=n×ln∑i=1npi−oi2n+2kwhere n is the total number of years, k is the number of parameters in the model, p_i and o_i represent the predicted and observed DFS in year i, respectively.

To obtain future ALAN data, we performed ordinary least squares regression (Eq. (15)) on the observed ALAN values using each city’s average ALAN from 2001 to 2022 and annual GDP data. Model accuracy was evaluated using five-fold cross-validation (Supplementary Fig. 18), and the resulting estimates were used for subsequent analyses.15\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\rm{ALAN}}}_{i}=a \times {{\rm{ALAN}}}_{{mean}}+b\times {\mathrm{ln}}\left({{\rm{GDP}}}_{i}\right)+c$$\end{document}ALANi=a×ALANmean+b×lnGDPi+cwhere ALAN_mean is the average ALAN, and GDP_i is the GDP in year i, a, b, and c are the coefficients of the regression equation.

Further information on research design is available in thelinked to this article. Nature Portfolio Reporting Summary