Autumn leaf color brightness of Japanese alpine vegetation is projected to decrease under future climate change

TL;DR

Future climate change is expected to reduce the brightness of autumn leaf colors in Japanese alpine vegetation, affecting tourism and ecosystems.

Contribution

This study identifies green-up timing as a key factor influencing autumn leaf color brightness and projects future declines under climate change scenarios.

Findings

Earlier green-up timing correlates with duller autumn leaf colors due to reduced physiological vitality.

Climate models predict a decrease in autumn color brightness by up to 15% by the end of the century.

Projected color reduction could impact tourism and ecosystem productivity.

Abstract

Autumn leaf coloring is an essential cultural ecosystem service, but mechanisms of color brightness have scarcely been focused on, even though the autumn color crisis (leaf fall without coloring) has been reported in Japanese alpine shrubs. To approach the mechanisms, we analyzed the relationship between brightness and environmental conditions and projected possible future changes. Autumn color brightness was measured by repeated camera observations at three Japanese alpine sites. Environmental factors were analyzed using linear mixed model comparisons. Results showed earlier green-up timing was related to a duller autumn leaf color, possibly because a prolonged leaf period creates older leaves and reduces physiological vitality in autumn, affecting anthocyanin biosynthesis. Green-up timing based on snowmelt day and degree-day-based prediction of snowmelt day predicted earlier snowmelt and green-up in the future, with decreases in autumn color brightness. The amount of color reduction varied among climate models and scenarios, from a decline of several percent under MRI-CGCM3 RCP 2.6 to around 15% under MIROC5 RCP8.5 by the end of this century. Projected reductions in autumn color brightness could have an economic impact on tourism, and it could also be linked to modifications in material cycles and the ecosystem’s productivity.

The online version contains supplementary material available at 10.1038/s41598-025-14547-2.

Figures

Response lines of the best model. Different symbols represent areas of interest (AOIs), and filled symbols mean red coloring. Note that the ASH site did not contain yellow coloring AOIs. Regression lines by the linear mixed model 5 in Table [1](#Tab1) are plotted.

Maps of current (1981–2000) and future (2081–2100) projected autumn leaf color brightness and its future retained rate. Differences among GCMs and RCPs are shown in Supplementary Fig. S4.

Density plot of future retention of autumn leaf color brightness among GCMs and RCP scenarios.

The future retention of autumn leaf color brightness (VARI) and the current temperature distribution map. The scatter plot between future retained VARI (Fig. [2](#Fig2)) and current annual mean temperature (1981–2000) for the alpine vegetation cells shows that the future VARI decline is evident in current high temperature cells (i.e., temperature > 3.5 °C). Future VARI change was projected using the MIROC5 RCP8.5 scenario in 2081–2100.

Full text

Introduction

Autumn in temperate to sub-boreal regions in the Northern Hemisphere is marked by red, yellow, orange, and brown leaves. This autumn leaf coloring is an important cultural ecosystem service^1,2^ having commercial importance for tourism industries^3,4^ by creating beautiful landscapes full of color^5,6^. Many websites report observed or projected autumn leaf vistas (e.g., https://smokymountains.com/fall-foliage-map/, accessed on 2025/5/7). The timing of this event is reported to be described by several climatic drivers such as temperature^7–9^ day length, especially at high latitudes^10,11^ and drought^12^. This timing has economic importance when it aligns with holiday seasons or calendar-based festivals^13^. But the phenological timing is changing with climate change^14^ leading to potential changes in economic income^15^.

Aside from a culturally important ecosystem service, autumn leaf coloring is also an important marker of ecosystem functional processes. The phenological productive period (from leaf expansion to senescence) strongly drives above-ground productivity^16,17^. Leaf color pigments also support productivity through physiological protection for photosynthetic systems from several stress factors^18–20^. Such above-ground productivity directly affects below-ground activities through the quantity and quality of litterfall, affecting the rate of decomposition^21,22^ and root exudates^23^. These factors link ecosystem-scale functions such as productivity and the nitrogen and carbon balances^24,25^ suggesting a need to understand detailed mechanisms in autumn leaf coloring.

However, in addition to the limited reports on autumn phenology compared to spring^26^ previous studies on autumn leaf coloring have three limitations. Firstly, they have focused on timing (phenology), ignoring color brightness, despite its primal human appeal^5,6^. Autumn reds and yellows are due to the accumulation of anthocyanins and carotenoids, with a simultaneous reduction of chlorophylls^27^. Temperature and solar radiation determine the physiological response of these pigments and the brightness of autumn leaf color^28–30^. Stress factors such as solar radiation, drought, nutrient deficiency, salinity, and pests and diseases also affect pigment dynamics^19,31,32^. These reports, however, derive mainly from laboratory experiments or individual observations, so the narrow spatiotemporal scale is the second limitation, suggesting the need for a broader ecological-scale analysis. Which climatic factor or stress condition is important for local autumn color brightness is unclear at broad spatial scales, though there are a lot of broad-scale reports on phenological autumn color timings^33,34^. Future climate change could cause both increases and decreases in autumn leaf color brightness, but we found no reference to these in our review. The third limitation is that previous works have been focused mainly on temperate trees, with little attention to alpine vegetation^35^. Human habitation is concentrated in the temperate zone, so autumn leaf coloration there is important, but alpine autumn scenery is also important, especially for tourism. It is reported that autumn leaf coloring in Japanese alpine vegetation is eventually disappearing in some sites and years^35^ so an approach to the mechanisms underlying alpine autumn color brightness is necessary both academically and socially.

The aims of this paper are to clarify the relationship between autumn leaf color brightness of alpine vegetation and environmental conditions, and to project possible future changes in brightness in Japanese alpine vegetation. We used repeated camera observations and climate datasets to construct a linear mixed model to predict autumn color brightness. Using the model and future climate models, we projected future changes in the autumn color brightness on a national scale. We posed two key questions: Which of the analyzed environmental variables is most strongly correlated with the autumn leaf color brightness in Japanese alpine vegetation? Will autumn leaf color brightness increase or decrease under future climate warming scenarios?

Results

Green-up day affected autumn leaf color brightness

After model comparisons between the response (i.e., annual maximum autumn color brightness, which was indexed by the visible atmospherically resistant index (VARI)) and predictor variables (Supplementary Table S1), the model containing green-up day and July mean temperature best explained the annual maximum VARI (model 1, Table 1). The green-up day had a higher standardized parameter value (0.069 ± 0.012, estimated parameter by model 1 and its standard error, n = 160) than July mean temperature (0.056 ± 0.013), and it was selected in all top 10 models, suggesting that it is the most important variable for autumn leaf color brightness in our sites. On the other hand, July temperature was selected in only 3 of the top 10 models. VARI showed vague responses to other explanatory variables (Supplementary Fig. S1), suggesting low explanatory power. Considering relatively weaker effects of the other variables and following wide uncertainty in future projections, we used the green-up day only model (model 5) to predict future changes in autumn leaf color brightness, though the coefficient of determination was the lowest in model 5 among the top 10 models (Table 1). As the green-up day and annual peak VARI had a positive correlation (Fig. 1), early green-up years would have duller autumn leaf colors according to this model.

Table 1. The best 10 model structures and their performances.ModelStructureAIC R ^2^ 1Gup + Temp_7−130.30.512Gup + Prec_8−126.50.513Gup + Temp_7 + GSR_8−124.10.514Gup + Prec_8 + Sun_8−122.80.535Gup−122.20.476Gup + Temp_7 + Sun_8−121.50.517Gup + Prec_8 + Sun_7−120.10.528Gup + GSR_8 + Prec_8−118.90.519Gup + GSR_7 + Prec_8−117.50.5110Gup + Sun_7−116.00.48“_7”, July; “_8”, August. Gup, green-up day; Temp, mean temperature; Tmax, maximum temperature; Prec, precipitation; GSR, global solar radiation; Sun, sunshine duration. AIC, Akaike’s information criteria.

Fig. 1. Response lines of the best model. Different symbols represent areas of interest (AOIs), and filled symbols mean red coloring. Note that the ASH site did not contain yellow coloring AOIs. Regression lines by the linear mixed model 5 in Table 1 are plotted.

Future VARI reduction through earlier snowmelt and green-up

To predict future changes in autumn leaf color brightness, we predicted future green-up days using the correlation between green-up days and snow-melt days in our research site (Supplementary Fig. S2) and a degree-day-based snow projection model (details are in the Materials and methods part). Using the predicted future green-up days, we predicted future autumn color brightness. Under future warming and subsequent advancement of snowmelt and green-up days (Supplementary Figs. S2, S3), alpine VARI was projected to decrease nationally (Fig. 2, Supplementary Fig. S4). The decrease varied widely among GCMs and RCPs, from several percentage points under MRI-CGCM3 RCP2.6 2031–2050 to ~ 15% under MIROC5 RCP8.5 2081–2100 (Fig. 3). In some cells, it reached 25% reduction under MIROC5 RCP8.5 2081–2100 (Figs. 2 and 3). The degree of reduction was consistent with the degree of warming (Supplementary Fig. S3), suggesting greater VARI reduction under severe warming scenario. MRI-CGCM3 RCP2.6 2031–2050 predicted an increase of < 1 °C in annual mean temperature, but MIROC5 RCP8.5 2081–2100 predicted around 5 °C. Warmer conditions would advance snowmelt and green-up days, creating duller autumn leaf colors.

Fig. 2. Maps of current (1981–2000) and future (2081–2100) projected autumn leaf color brightness and its future retained rate. Differences among GCMs and RCPs are shown in Supplementary Fig. S4.

Fig. 3. Density plot of future retention of autumn leaf color brightness among GCMs and RCP scenarios.

Vulnerable areas

Spatial differences in future VARI changes were related to current temperature conditions: sites with greater VARI reduction appeared in warmer areas, especially those whose current mean temperature is > 3.5 °C (Fig. 4). Such warmer sites at lower elevations represent populations of alpine plants that are vulnerable not only in their distribution^36^ but also in autumn leaf color brightness. Vulnerable alpine vegetation was predicted in southern Hokkaido, such as on Mt. Hokkaidō-Komakgatake, Mt. Esan, and Oshima-Ōshima (Figs. 2 and 4). Another typically vulnerable region was predicted in the northern part of the central region, such as in the Echigo-Komagatake mountains.

Fig. 4. The future retention of autumn leaf color brightness (VARI) and the current temperature distribution map. The scatter plot between future retained VARI (Fig. 2) and current annual mean temperature (1981–2000) for the alpine vegetation cells shows that the future VARI decline is evident in current high temperature cells (i.e., temperature > 3.5 °C). Future VARI change was projected using the MIROC5 RCP8.5 scenario in 2081–2100.

Discussion

The green-up day was the most important variable (Table 1; Fig. 1) to predict autumn color brightness in our research sites. This result is inconsistent with previous studies reporting that warmer autumn temperature was correlated with autumn red coloration in tall temperate deciduous trees in North America^28,30^. Other papers reported that solar radiation was correlated with the autumn color brightness, as greater solar radiation in open canopies than in closed canopies increased the amount of carotenoids in European Fagus sylvatica seedlings^31,37^. Although the green-up day was noted as an important factor for autumn phenological timing^33,34^, its effect on color brightness was not previously assessed. The mechanisms behind spring-autumn phenology connection are still under discussion, but earlier spring leaf-out may cause: (1) earlier autumn senescence if leaf life span is fixed^12,38^; (2) water shortage in the early period of the growing season^39^; (3) damage from spring frost^40^ and insect outbreaks^41^; and (4) earlier achievement of maximum carbon content, which advances senescence in autumn^42^. Although further quantitative analyses comparing the effect of green-up day with the other temperature and radiation on autumn color brightness at the other sites are needed, we think the leaf lifespan matters here. Relating to the previous mechanism (1) and (4), the duration of leaf lifespan is prolonged in early green-up years, resulting in older leaves in autumn. Our data show a positive effect of green-up day on peak autumn color day (peak autumn color day ≈ 0.359 × green-up day, using random intercept for each individual detected by the area of interest (AOI)), associating an early green-up day with an early peak autumn color day, but this autumn earliness was weaker than spring earliness (green-up day) (the slope parameter (0.359) was lower than 1), resulting in older leaves in an early green-up year. This early autumn leaf senescence in an early spring year can be listed as a main cause of the correlation between autumn color brightness and green-up day in Japanese alpine regions. Older leaves tend to have decreased physiological performance, possibly reducing the biosynthesis of pigments^43,44^. This prolonged leaf duration effect would be more serious in less productive alpine zones than in warmer regions, due to the limited resources in alpine zones for creating enough defensive and adaptive physiological mechanisms.

Future reduction of alpine autumn color is evident in warmer lower-elevation or latitude sites (Fig. 4), suggesting these areas as vulnerable to a reduction of autumn color brightness. Regions whose current annual mean temperature is > 3.5 °C were suggested to expect a severe VARI decrease (Fig. 4). These hot regions would be in or near cool temperate regions (below alpine zone) according to the species distribution model of Japanese stone pine (Pinus pumila), a typical alpine species, for which the warmth index (Sum of (Tm-5) when Tm is above 5 °C: (Tm: Monthly Mean Temperature)) below 43.2 °C month (i.e., an annual mean temperature of ~ 4 °C according to the mesh climate Japan 2010 dataset) is the first requirement of its distribution^36^. This low-elevation alpine plant distribution^45^ is possibly a result of stress conditions caused by volcanic activity^46^ heavy snowfall^47^ cool microsites such as wind-hole cool spots^48^ or algific talus^49^. Our projection area included these special areas because of their vegetation structure and component species, but caution is needed because some of this low-elevation alpine vegetation would not contain our target species (S. matsumurana, B. ermanii, and A. tschonoskii).

Although the projected color change has wide uncertainty (i.e., VARI decrease of up to 25%) according to the predicted rate of warming among GCMs and RCPs (Fig. 3), the trend of color decrease was consistent. Such a decrease would cause widespread economic loss in the tourism industry. However, this prediction contradicts previous work predicting increases in economic income by Chinese tourists in Japanese temperate regions where a later autumn leaf color season will coincide with Chinese national holidays^15^. Although these authors assessed economic effects by focusing on color timing, it is clear that color brightness is also relevant. In any case, our results demonstrate that phenological perspectives alone are not enough in assessing alpine autumn leaf color tourism, and future changes in autumn color brightness should be considered for this important cultural ecosystem service.

Future reduction of autumn color brightness could also affect ecosystem productivity or the C: N balance, although linkage to a mechanistic process-based vegetation model (especially for alpine vegetation, e.g^50,51^. is needed for its quantitative assessment. The home-field advantage of the decomposing process^52^ should amplify the complexity of ecosystems through spatiotemporal matching of decomposing ability and litter quality and quantity, especially when climate change and related vegetation shifts occur. Our work is a macro-scale, simple projection, but there is potential to connect and broaden research on ecosystem functions.

Our focus on macro-scale national patterns could lead to mismatches in fine-scale predictions. The site cameras (Supplementary Fig. S5) cover limited slope direction and topography. Snow cover tends to differ with slope direction and topography^53,54^ so our projections have uncertainties in relation to these landscape micro-patterns. Local management plans need future fine-scale analyses and projections. Another limitation is local adaptation or phenotypic plasticity on physiological performance^55,56^. The spatial extent of our observation data is limited, then it possibly causes uncertainty, especially for the low-elevation alpine zones. Future modeling should focus on finer spatial scales. This detailed research on individual differences (e.g., phenotype, micro environment, age, etc.) can increase the model accuracy, which is another limitation of our model 5 (R^2^ < 0.5, Table 1). Although we focused on several typical Japanese alpine species, future climate change has the possibility of causing compositional changes^57^. This ecological effect of competitions and distributions should be mentioned further.

Conclusion

Previous phenological studies focused on the timing of autumn leaf coloring. A predicted timing match with holiday seasons under future climate conditions gives a predicted increase in economic income from tourism^15^. However, our predicted decrease in autumn leaf color brightness rather suggests a possibility of decreasing income. More attention to autumn leaf color brightness is needed, for which our broad-scale modeling can be a first step toward revealing the underlying mechanisms. Future work should approach the finer-scale ecophysiological mechanisms to support forest management under future climate conditions. Fine-scale modeling can be used to predict refugia^58^ for bright autumn colors under future conditions, which should be a target area for conservation or tourism use. It might be possible to artificially mimic and manipulate such environments (e.g., by accumulating snow in specific valleys) to retain bright autumn vistas, although careful ethical discussion is needed, as human interventions could affect the integrity of these largely unspoiled alpine regions. Approaches to mechanisms of local adaptation and assisted gene flow from populations in warmer locations^59,60^ should be considered to stabilize autumn color brightness, although ethical discussion is needed for this approach too. For such analyses of ecophysiological mechanisms of autumn leaf color, broader observations on a finer spatial scale are crucial, especially in alpine vegetation. Potential data sources include fine-scale remote-sensing data from satellites (e.g., Dove satellite constellation) and drones^61^ civil science-based observation networks or analysis of data from social network services^62^ and an increase in repeated camera observations in collaboration with mountain huts or park rangers^63^. We hope that this will stimulate such future works.

Materials and methods

Overall flowchart of our analyses was summarized in Supplementary Fig. S6.

Repeated camera observations

We used digital repeat photography data captured at three alpine sites in Japan (details are in Koide, et al.^35^: Asahidake (ASH, 43°39′44″N, 142°49′31″E, 1600 m a.s.l.), in Hokkaido, and two sites in the central alpine area of Honshu, Murodo (MRD, 36°34′40″N, 137°36′04″E, 2450 m a.s.l.), on the Sea of Japan side, and Gokurakudaira (GOK, 35°46′39″N, 137°48′47″E, 2700 m a.s.l.), on the Pacific Ocean side (Supplementary Fig. S5). We downloaded the data for ASH (2016–2020) from the Internet Nature Institute (https://www.sizenken.biodic.go.jp/index_en.php, accessed on 2025/5/7). For GOK (2012–2016) and MRD (2011–2016), we used original high-resolution photographs from the Center for Global Environmental Research (https://db.cger.nies.go.jp/gem/ja/mountain/, accessed on 2025/5/7; note that the resolution is reduced). We used different year range for ASH because the camera and angle settings changed in June 2016 in this site. Camera manufacturer and model were Panasonic BB-SW174WA for ASH, Nikon D7000 for GOK, and Canon EOS 5D Mk2 for MRD. The mean annual temperature at the sites ranges from − 1.8 to + 0.1 °C (ASH: −1.8, GOK: +0.1, MRD: −0.4) and annual precipitation from 1591 to 2714 mm (ASH: 1591.3, GOK: 2714.1, MRD: 2288.7) (1981–2010 average) according to 1-km-mesh climate data^64^. All sites have alpine vegetation in a mosaic of several dominant species, including Pinus pumila, Sorbus matsumurana, Betula ermanii, and Acer tschonoskii. The digital photographs were taken hourly; photographs captured around noon were selected in consideration of sunshine and the slope direction of each site (ASH: 12:00–14:00, GOK: 9:00–16:00, MRD: 10:00–14:00). Ten areas of interest (AOIs) were set (Supplementary Fig. S5), and average values of each of the red, green, and blue channels in each AOI were calculated. We calculated the visible atmospherically resistant index (VARI) as an indicator of the leaf color brightness as VARI = (Red − Green)/(Red + Green − Blue). Note that the numerator was converted from the ordinal VARI index (i.e. Green − Red) to focus on red. We used VARI for its relatively robust performance at explaining autumn leaf color brightness compared to other indices^35^. We used the maximum VARI value of each AOI in each year as the autumn color brightness in the analyses. The images were processed in the free scripting language ZeScript and the free ZeGraph software (http://www.zegraph.com/). To focus on autumn leaf coloring, we selected photos from 1 September to the first snowfall of the season (early to mid-October) for analysis. The number of photographs for each year was around 183, 315, and 250 for ASH, GOK, and MRD, respectively. At ASH, all AOIs were set in red leaf areas. At GOK and MRD, five and six AOIs, respectively, were set in red leaf areas, and the others were set in yellow leaf areas. Sorbus matsumurana is the main red-leafed species at our sites, although Sorbus sambucifolia may also have been present in the AOIs. Betula ermanii and Acer tschonoskii are the main yellow-leafed species.

To detect the effects by spring phenology, we detected the green-up day and the snowmelt day from spring camera observation data (May–July) from each AOI. The snowmelt day was detected by eye, because automated detection misidentifies by observing snow covers behind the vegetation through gaps in bare branches. The green-up day was detected by fitting a sigmoid curve to the seasonal changes in the green ratio (GR) index (GR = Green/(Red + Green + Blue)). The day when the sigmoid-fitted GR growth rate first reached the maximum was selected as the green-up day according to the previous work^63^.

Climate data

Past climatic data (Supplementary Table S1) were gathered from on-site observations or the nearest AMeDAS (Automated Meteorological Data Acquisition System) weather station of the Japan Meteorological Agency (JMA; https://www.data.jma.go.jp/stats/etrn/index.php, accessed on 2025/5/7). Daily temperature, precipitation, sunshine duration, and snow depth data were collected. At GOK (35°46′, 137°48′, 2630 m a.s.l.), we collected hourly temperature data with a HD9817T1R sensor (DeltaOHM; http://ims.shinshu-u.ac.jp/~metims_web/index.php?sokuhou, accessed on 2025/5/7). At MRD, we used temperature data provided by the Tateyama-Murodou mountain villa (http://www.murodou.co.jp/kishou/test.php, accessed on 2025/5/7). For other climatic data, the nearest AMeDAS stations around GOK (the average among Kiso-Fukushima (35°50.4′N, 137°41.3′E, 750 m a.s.l), Ina (35°49.5′N, 137°57.3′E, 633 m a.s.l.), Iijima (35°39.2′N, 137°53.9′E, 728 m a.s.l.) for most variables, and Kaida-Kougen (35°56.3′N, 137°36.1′E, 1130 m a.s.l.) for the maximum snow depth) and MRD (the average among Omachi (36°31.4′N, 137°49.9′E, 784 m a.s.l.) and Kamiiichi (36°40.2′N, 137°25.4′E, 296 m a.s.l.) for most variables, and Omachi for the maximum snow depth) was used. At ASH, we collected temperature and other data from the nearest AMeDAS stations (Shibinai (43°38.6′N, 142°34.9′E, 310 m a.s.l.) for most variables, and Sounkyo (43°45.2′N, 142°55.8′E, 540 m a.s.l.) for the maximum snow depth). Since the elevational range among AOIs at GOK and MRD is wide, we accounted for the elevational temperature difference by using a lapse rate of 0.55 °C/100 m^65^. We collected global solar radiation data from the Agro-Meteorological Grid Square Data (AMGSD) of the National Agriculture and Food Research Organization (https://amu.rd.naro.go.jp/wiki_open/doku.php?id=start2, accessed on 2025/5/7). This dataset is spatially interpolated at a 1-km resolution through statistical modeling. In projecting the spatial distribution of autumn leaf color brightness, we used the AMGSD as current climatic data.

We used climatic projection data of the World Climate Research Program’s Coupled Model Intercomparison Project Phase 5 as future climate data. We used the Social Implementation Program on Climate Change Adaptation Technology (SI-CAT: https://research-er.jp/categories/4160, accessed on 2025/5/7) to statistically downscale the original global climate models (GCMs) to a 1-km scale. We used two GCMs (MIROC5, MRI-CGCM3), two representative concentration pathway scenarios (RCP 2.6, 8.5), and two time-periods (2031–2050, 2081–2100) to test future climatic uncertainties. The two GCMs were made by Japanese research institutes, and MIROC5 shows a higher temperature warming trend (Fig. S3). The baseline condition was set at 1981–2000. The bias in each 1-km cell was corrected first by calculating temperature differences and the ratio of precipitation between current (1981–2000) and future time-periods in the downscaled GCMs, and then by adding the temperature difference to or by multiplying the ratio of precipitation difference by the current AMGSD data.

Prediction of snowmelt day

To estimate future changes in green-up days, we focused on the snowmelt day, which has a strong correlation with the green-up day (Supplementary Fig. S2). We based a snow projection model on the degree-day method^66^ as used to project future distribution changes among Japanese deer^67^. The daily existence of snow cover is predicted from the mean temperature and precipitation in 1-km cells. The snowmelt coefficient was optimized by using SPOT Vegetation satellite data^68^. From the predicted daily snow existence in the present and future periods, we set the snowmelt day as the day following the last snow day. Then we predicted the current and future green-up days from the snowmelt day with the equation relating them (Supplementary Fig. S2).

Analysis

To predict the annual maximum VARI in each AOI, we constructed linear mixed models for every combination among 14 explanatory variables (Supplementary Table S1). The explanatory variables included monthly variables in each July and August (the productive period for alpine vegetation without snow cover). The broader seasonal range (June to September) was also analyzed in the same way, but the results were almost the same (Supplementary Table S2), then we report results in July to August in the main text. We used temperature (mean/maximum/minimum)^28^^[–30^, global solar radiation and total sunshine duration^31,37,69^ total precipitation^28^ maximum snow depth^70^ and the annual green-up day^33,34,71,72^ for their effects on autumn phenology and color amount. We set leaf color type (red/yellow) for each AOI as a fixed parameter, because yellows have lower VARI values than reds. Since our model treats the VARI difference between red and yellow coloring as constant (Fig. 1), we only showed predictions for the red coloring in our maps and graphs. Each AOI was used as a random intercept, as this analysis focused on temporal patterns at each AOI to predict future changes, given the limited climatic and topographic research range of our dataset. All possible combinations among explanatory variables were compared using Akaike’s information criteria (AIC). Models with significant correlations (p < 0.05) among explanatory variables were removed.

Supplementary Information

Below is the link to the electronic supplementary material.

Supplementary Material 1

Supplementary Material 2