How did the characteristics of the growing season change during the past 100 years at a steep river basin in Japan?

The effects of climate change on plant phenological events such as flowering, leaf flush, and leaf fall may be greater in steep river basins than at the horizontal scale of countries and continents. This possibility is due to the effect of temperature on plant phenology and the difference between vertical and horizontal gradients in temperature sensitivities. We calculated the dates of the start (SGS) and end of the growing season (EGS) in a steep river basin located in a mountainous region of central Japan over a century timescale by using a degree-day phenological model based on long-term, continuous, in situ observations. We assessed the generality and representativeness of the modelled SGS and EGS dates by using phenological events, live camera images taken at multiple points in the basin, and satellite observations made at a fine spatial resolution. The sensitivity of the modelled SGS and EGS dates to elevation changed from 3.29 days (100 m)−1 (−5.48 days °C−1) and −2.89 days (100 m)−1 (4.81 days °C−1), respectively, in 1900 to 2.85 days (100 m)−1 (−4.75 days °C−1) and −2.84 day (100 m)−1 (4.73 day °C−1) in 2019. The long-term trend of the sensitivity of the modelled SGS date to elevation was −0.0037 day year−1 per 100 m, but the analogous trend in the case of the modelled EGS date was not significant. Despite the need for further studies to improve the generality and representativeness of the model, the development of degree-day phenology models in multiple, steep river basins will deepen our ecological understanding of the sensitivity of plant phenology to climate change.


Introduction
In temperate and subarctic climate zones, plant phenological events such as flowering, leaf flushing, leaf colouring, and leaf fall show spatiotemporal characteristics along latitudinal, longitudinal, and elevational gradients (e.g., [1][2][3][4][5][6][7][8][9][10]). In particular, the timing of the start (SGS) and end of the growing season (EGS) are more sensitive to elevational gradients than to latitudinal gradients [8]. This difference in sensitivity is caused by the difference in the sensitivity of temperature to gradients associated with elevation versus gradients associated with latitude and longitude, because temperature affects the SGS and EGS dates (e.g., [11,12]). For instance, the sensitivities of temperature to elevation and latitude are 154 m˚C −1 and about 122 km C −1 , respectively [13]. The sensitivity to elevation suggests that the effect of climate change on plant phenology may be apparent in a steep river basin within a distance of no more than a few tens of kilometres, much shorter than the latitudinal scales of countries and continents (hundreds to thousands of kilometres). For instance, a 3˚C increase in annual average temperature would shift the climate by an amount corresponding to a decrease in elevation of about 460 m. For this reason, examination of the relationship between SGS and EGS dates versus climate change along an elevational gradient is a useful way to estimate the spatiotemporal changes of SGS and EGS dates in a warmer climate. Previous studies examined the spatiotemporal distribution of the timing of leaf flushing, leaf colouring, and leaf fall along an elevational gradient using data from in situ observations [14][15][16][17], phenological images [7,18,19], herbarium records [10], phenological information published on the Internet [9], data from satellite observations [6,12,[20][21][22][23], and modelling [24,25]. The scales of the targeted regions have been national (Germany [15], Slovakia [17], and China [12,22,24]), continental (Europe [14] and North America [7,10]), and global [6]. In contrast, previous studies that focused on the scale of river basins (10−100 km) involved relatively limited spatiotemporal scales (a mountainous region in central Japan and 12 years [26]; a mountainous region in central Japan and 6 years [27]; a mountainous region in central Japan and 3 years [19]; a mountainous region in central Japan and 1 year [28,29]; the Pyrenees in France and 3 years [20]; the Pyrenees in France and 2 years [11]; and the Alps in Germany [18,30]). Four considerations may account for these relatively small spatiotemporal scales. First, the steep basins that are found in the mountainous region of central Japan (the Japan Alps) and the Alps, in which there are large differences in elevation within a short horizontal distance (1000-2000 m in elevation versus a horizontal distance of 10-50 km), are not widely distributed in the world. Second, long-term, in situ observations at multiple points in a river basin are labour-intensive and expensive. Third, satellite sensors with a spatial resolution of 500-1000 m (MODIS [Moderate Resolution Imaging Spectroradiometer] sensors mounted on Terra and Aqua satellites and the VEGETATION sensor mounted on the SPOT [Satellite pour l'Observation de la Terre] satellite) cannot accurately detect the spatiotemporal variability of the SGS and EGS dates because of the heterogeneity of vegetation [31] and microtopography. Fourth, data obtained from these satellite observations have been limited to only the past 20 years.
Taking account of these issues will require carrying out the following tasks in an integrated manner: (1) the long-term SGS and EGS dates should be calculated using a degree-day phenology model that is based on long-term, continuous biometeorological observations [1,16,[32][33][34][35]; (2) the accuracy of the modelled SGS and EGS dates should be assessed by using longterm, continuous phenological images taken at multiple sites [31,36]; and (3) the spatiotemporal distribution of the SGS and EGS dates should be thoroughly evaluated around validation sites by analysing satellite observations with a fine spatial resolution (e.g., 10 m) [31]. In particular, to evaluate the long-term effects of climate change, we need to examine the spatiotemporal variability of SGS and EGS dates for at least the past 120 years, during which continuous observations of temperature have been recorded with modern-day technology. Making connections between plots (by in situ observations) and across regional scales (by satellite observations) [37] as well as performing century-scale phenological observations at the community scale [38] are especially important tasks that have been carried out to a limited extent in previous studies.
Accordingly, we calculated the SGS and EGS dates from 1900 to 2020 by using a degree-day phenological model that was developed for a cool-temperate, deciduous, broad-leaved forest site located in the upper reaches of the basin [32]. We then assessed the accuracy of the modelled dates by using daily images of vegetation phenology taken in the upper reaches of the basin [39,40], daily live camera images taken in the lower reaches [41], and satellite observations with a spatial resolution of 10 m. The goals of this study were (1) to examine the characteristics of the spatiotemporal variability of the SGS and EGS dates in a steep river basin under the influence of climate change on a century timescale and to identify what caused the variability, and (2) to assess the generality and representativeness of a degree-day model of vegetation phenology.

Daily mean air temperature
Air temperatures are recorded by AMeDAS stations installed by the Japan Meteorological Agency at about 840 sites in Japan. The interval between sites is about 21 km (https://www. jma.go.jp/jma/kishou/know/amedas/kaisetsu.html, accessed 6 July 2021). However, there are only 12 AMeDAS stations at elevations above 1000 m a.s.l. [13]. It is therefore difficult to obtain temperature data at multiple points in a steep river basin such as the Daihachiga River basin. In fact, there was only one AMeDAS station within our target area (i.e., the Takayama- AMeDAS station). However, meteorological observations have been conducted at a weather station about 500 m south of the TKY-tower (TKY-weather station; 36˚08 0 34@N 137˚25 0 20@E, 1342. a.s.l.) since January 1996 ( [32]; https://www.green.gifu-u.ac.jp/takayama/, accessed 6 July 2021).
We assumed that daily mean temperature data were not typically available at high elevations. We estimated daily mean temperature at TKY-tower by using a temperature lapse rate of 0.6˚C (100 m) −1 [15,26,51,52] and daily mean temperatures observed at Takayama-AMe-DAS from 1 January 1990 to 31 May 2020 (https://www.data.jma.go.jp/gmd/risk/obsdl/index. php, accessed 6 July 2021). We assumed that the heat island effect on temperature caused by urbanization was relatively low at Takayama-AMeDAS. Despite merging of municipalities in 2005, the population of the city of Takayama increased from 63,520 in 1920 to only 89,182 in 2015 (https://www.city.takayama.lg.jp/_res/projects/default_project/_page_/001/011/771/ 18gou_siryou2.pdf, accessed 6 July 2021). In the Discussion, we consider the effect on the calculated SGS and EGS dates of the seasonal gap between the estimated daily mean temperature at TKY-tower and the daily mean temperature that was actually observed at TKY-weather station.
To evaluate the characteristics of climate change, we examined the timeseries of 30-year monthly mean temperatures (climatological mean) and their standard deviations (SDs) during March−May and September−November at Takayama-AMeDAS calculated at intervals of 10 years. In accordance with the Japanese Meteorological Agency, the climatological mean was defined as the average during 30 years at intervals of 10 years. We assumed that the year-toyear variability of monthly mean temperatures (i.e., the SDs) had an important effect on the modelled SGS and EGS dates.

Degree-day model of vegetation phenology
We used a degree-day model of vegetation phenology that we developed by analysing the relationship between daily mean temperature and daily phenological images from 2004 to 2011 at TKY-tower [32]. The timing of leaf flush (i.e., SGS date) as well as leaf colouring and leaf fall (i.e., EGS date) in Japan are explained mainly by temperatures [8,9,53,54]. The SGS date was defined as the first day when the cumulative effect temperature (CET) was greater than 255.4˚C. This CET SGS was calculated with Eq (1), where we set the start date to be 1 January and the threshold temperature for the CET (T t,SGS ) to be 2˚C. In contrast, the EGS date was defined as the first day when the CET was less than −375.1˚C. This CET SGS was calculated with Eq (2), where we set the start date to be 1 August and the threshold temperature for CET (T t, EGS ) to be 18˚C [32].
https://doi.org/10.1371/journal.pone.0255078.g001 dormancy and a heating requirement after the release of bud dormancy [34,56,57]. However, we did not incorporate a chilling requirement for the release of bud dormancy, because the basin of the Daihachiga River is sufficiently cold in winter (December to March) to satisfy the requirement for a dormancy period. Climatological means during winter ranged from −1.4˚C in January to +2.9˚C in March (https://www.data.jma.go.jp/obd/stats/etrn/index.php, accessed 6 July 2021). We therefore expected to reduce the degrees of freedom of the degree-day model of vegetation phenology by not including a chilling requirement.
First, we calculated the SGS dates from 1900 to 2020 at TKY-tower and the EGS dates from 1900 to 2019 at Takayama-AMeDAS. Second, we examined the long-term linear trends from 1900 to 2019/2020 and timeseries of 30-year monthly SDs, calculated at intervals of 10 years, of modelled SGS and EGS dates at TKY-tower, Takayama-AMeDAS, and their difference. The climatological mean of temperatures was defined as the 30-year average temperature at intervals of 10 years. Finally, to evaluate the effect of climate change, we examined the correlations between the SDs of 30-year monthly mean temperatures at timesteps of 10 years and the SDs of 30-year modelled SGS and EGS dates at TKY-tower, Takayama-AMeDAS, and their difference.

Camera images
To assess the generality of the modelled SGS and EGS dates, we used daily images of phenology taken at TKY-tower [39,40] and live camera images taken at two points near Takayama [58,59]. At Koito-Pottery and Kaji-bridge, the latest live camera images are publicly available on an hourly or 3-minute timeframe from the Web sites of the Koito Pottery (https://koitoyaki.com/livecam/, which we have not accessed since March 2021) and Takayama Printing Co, Ltd (https://www.takayama-dp.com/live/, accessed 6 July 2021). At Koito-Pottery, the live camera images showed mainly deciduous trees close to the cameras and distant deciduous and evergreen forests. At Kaji-bridge, the live camera images showed mainly deciduous trees close to the cameras and distant deciduous forests. We assumed that there was no difference between the SGS and EGS dates at Takayama-AMeDAS, Koito-Pottery, and Kaji-bridge, because the difference in elevation among the three sites is very small (max. 45 m). We have downloaded live camera images by running shell scripts and the crontab command on a Linux personal computer since December 2017. We obtained permission to download images from the Koito Pottery and Takayama Printing Co, Ltd.

Satellite data
To assess the spatial representativeness of the modelled SGS and EGS dates, we used surface reflectance data in the red, green, and blue (RGB) bands with a 10-m spatial resolution observed by SENTINEL-2A/2B satellites (https://sentinel.esa.int/web/sentinel/missions/ sentinel-2, accessed 6 July 2021). The modelled SGS and EGS dates may have been affected by the spatial heterogeneity and distribution of tree species as well as by microtopography. For this reason, analysis of SENTINEL-2A/2B satellite observations was a useful way to fill spatial gaps between modelled dates in a degree-day model of phenology and in situ observed data at a validation site [33]. We collected data observed within 5 days before and after modelled SGS and EGS dates at TKY-tower and Takayama-AMeDAS because the time interval between SENTINEL-2A and 2B satellite observations is 5 days (by both satellites). Despite the short observation period of the SENTINEL-2A/2B satellites (since 2015 for SENTINEL-2A and 2017 for SENTINEL-2B), the SENTINEL-2A and 2B satellites are the only non-commercial satellites that have provided accurate geographical distributions of plant phenology in a steep river basin with a high spatiotemporal resolution.
We downloaded satellite data from the CREODIAS web site (https://creodias.eu, accessed 6 July 2021). We collected data since 2018 because the SENTINEL-2A and 2B satellites were launched in 2015 and 2017, respectively. We collected only data with less than or equal to 20% cloud coverage. We used RGB composite images and calculated the green-red vegetation index (GRVI; [60]). In deciduous broad-leaved forests, a GRVI of 0 indicates both the start of leaf flushing [60] and the time of maximum leaf colouring and leaf fall [61]. In this study, atmospheric corrections were performed by using the SNAP (Sentinel Application Platform) toolboxes (http://step.esa.int/main/toolboxes/snap/, accessed 6 July 2021). We did not apply further cloud masking to the satellite data because there was no cloud contamination in our target river basin.
The SDs of the modelled SGS dates during 1940−1969 and 1980−2009 at TKY-tower and Takayama-AMeDAS, respectively, were large compared to the corresponding SDs during 1900−1929 and 1910−1939, respectively. The SD of the differences in dates between sites was larger from 1950−1979 to 1990−2019 than from 1900−1929 to 1920−1959 (Fig 3A). In contrast, the SDs of the modelled EGS dates at TKY-tower and Takayama-AMeDAS were larger from 1960−1989 to 1990−2019 than from 1900−1929 to 1940−1969. The interannual variability of the SDs of the difference in EGS dates between sites was smaller than that of the SDs of the modelled EGS dates at the sites (Fig 3B).

Air temperature
Timeseries of 30-year monthly mean temperatures, calculated at intervals of 10 years (Fig 4A  and 4C), showed significantly positive rates of change in March, April, May, September, October, and November (Table 1).
Compared with May, the SD of the monthly mean temperature in March and April was always large (Fig 4B). In contrast, compared with September and October, the SD of the 30-year monthly mean temperature in November was large from 1900−1929 to 1930−1959. In addition, the SDs have gradually increased since 1930 (Fig 4D). Table 2 shows the relationship between the SDs of the 30-year monthly mean temperatures, calculated at intervals of 10 years (Fig 4B and 4D), and the 30-year modelled SGS ( Fig 3A) and EGS (Fig 3B) dates at TKY-tower, Takayama-AMeDAS, and their difference. At TKY-tower, the SDs of the 30-year modelled SGS dates were correlated significantly with those of the monthly mean temperatures in April and May. In contrast, at TKY-tower and Takayama-AMeDAS, the SDs of the 30-year modelled EGS dates were significantly correlated with those of the monthly mean temperatures in September and October.

Assessment via camera images and satellite observations
At TKY-tower, images of phenology taken on modelled SGS dates in each year showed almost the same condition of the canopy surface (Fig 5A). This condition was almost the same as the condition that defined the SGS date based on an analysis of RGB values extracted from images of phenology by Nagai et al. (2013) [32,Fig 2]. In addition, the modelled SGS dates from 2004 to 2011 correlated significantly with the observed SGS dates shown in Table 1 of Nagai et al. (2013) (Spearman's rank correlation coefficient, ρ = 0,94, p < 0.001) [32]. In contrast, images taken on modelled EGS dates in each year showed small variations in the conditions of leaf colouring and leaf fall on the canopy surface ( Fig 5B). However, those images included the same conditions as the images of phenology associated with the EGS date defined by Nagai et al. (2013) [32,Fig 2]. In addition, the modelled EGS dates from 2004 to 2011 correlated significantly with the observed EGS dates shown in Table 1 (Fig 6A). In contrast, despite variations among individual trees, images taken on the modelled EGS dates at Takayama-AMeDAS showed leaf colouring and leaf fall (Fig 6B). Table 3 shows the relationships between the modelled SGS and EGS dates and the closest satellite-observed dates. Except for 2018, we could not obtain satellite data with less than or equal to 20% cloud coverage at about the times of the four modelled dates. The satelliteobserved dates were 2 or 3 days later than the modelled SGS and EGS dates. On the satelliteobserved date closest to the modelled SGS date at Takayama-AMeDAS (3 days later than the modelled SGS date), leaf flush had not occurred widely throughout the river basin (Fig 7A). The GRVI of deciduous, broad-leaved forests was less than or equal to 0.0 (Fig 7B). On the satellite-observed date closest to the modelled SGS at TKY-tower (2 days later than the modelled SGS date), leaf flush had expanded throughout the basin (Fig 7C), and the GRVI in deciduous, broad-leaved forests had increased to 0.2-0.4. In addition, the GRVI around TKY-tower had   increased to 0.1-0.4 (Fig 7D). In contrast, on the satellite-observed date closest to the modelled EGS at TKY-tower (3 days later than the modelled EGS date), leaf colouring and leaf fall had occurred at elevations greater than 1000 m (Fig 8A). The GRVI in deciduous, broad-leaved forests was less than or equal to 0.0 ( Fig 8B). On the satellite-observed date closest to the modelled EGS at Takayama-AMeDAS (2 days later than the modelled EGS date), leaf colouring and leaf fall had advanced throughout the basin (Fig 8C), and the area characterized by GRVIs less than or equal to 0.0 had expanded in the basin (Fig 8D).

Sensitivity of modelled SGS and EGS dates to elevation
The results of our analysis of the sensitivity of modelled SGS and EGS dates to elevation (Fig 2) implied the following relationships. The modelled SGS date at Takayama . However, the yearly difference of the modelled SGS dates at Takayama-AMeDAS and TKY ranged from −15 days (earlier than previous year) to 21 days (later than previous year) and from −16 days to 23 days, respectively (Fig 2A). These differences are equivalent to elevation changes of −526 m (shift to lower elevation) to +737 m (shift to higher elevation) and to elevation changes of −561 m to +807 m, respectively (in the case of 2019). The yearly difference of the modelled EGS dates at Takayama-AMeDAS and TKY ranged from −10 days (earlier than previous year) to 17 days (later than previous year) and from −12 days to 15 days, respectively (Fig 2B). These differences were equivalent to elevation changes of +352 m to −599 m and +423 m to −528 m, respectively (in the case of 2019). These comparisons suggest that deciduous, broad-leaved forests in the Daihachiga River basin might have acclimated to the year-to-year variability of climate during the past 120 years by changing their leaf flush and leaf fall phenology. In other Table 2. Summary of correlations (Spearman's rank, ρ) between SDs of 30-year monthly mean temperature at a time-step of 10 years (Fig 4B and 4D) and that of 30-year modelled SGS (Fig 3A) and EGS dates (Fig 3B)  words, the warmer spring temperature triggered earlier SGS dates, and the warmer autumn temperature triggered later EGS dates. However, this acclimation may not apply to other river basins where annual mean temperatures are much higher than the temperatures of the Daihachiga River basin. Nagai et al. (2020) [63] reported that the correlation between the observed first flowering date of Yoshino cherry and latitude has decreased since 1980 in Japan. They explained the lower correlation by reasoning that higher temperatures in winter had delayed the release from dormancy and the subsequent first-flowering dates in areas where the annual mean temperature was high. To evaluate the significance of the chilling requirement for dormancy release, we will need to examine the relationships between elevation and the dates of SGS and EGS in other steep river basins where annual mean temperatures are much higher than the temperatures in the Daihachiga River basin. Our result was similar to results previously reported for the SGS dates of oak and larch and for the flowering of horse chestnut, alder, cocksfoot grass, goat willow, rye, and small-leaved lime ( Table 4). The sensitivity we observed, however, was large compared with the sensitivity of the flowering of Norway spruce and common wine grapes, the SGS of beech, and other SGS and EGS dates inferred from satellite observations (absolute values) ( Table 4). The following reasoning may account for these similarities and differences. The degree-day phenology model that we used [32] was developed at TKY, where the dominant tree species are oak (Q. crispula) and birch (B. ermanii and B. platyphylla) [37,49,50]. Consequently, our results were similar to previously reported results for oak. In contrast, the timing and patterns of leaf flush and leaf fall in beech, which is sparsely distributed around TKY-tower, differs from those of oak and birch. For instance, the leaf flush of beech is earlier than that of oak and birch [64]. These phenological characteristics may explain the different sensitivities of the SGS and EGS dates to elevation. The characteristics of the sensitivities of flowering, leaf flush, and leaf fall to temperature have been reported in many previous studies [10,11,65,66]. The differences in reported sensitivities suggest that use of a single degree-day model of phenology to estimate SGS and EGS dates may be misleading with respect to those dates in a steep river basin where there is high diversity and an uneven distribution of tree species.
The following lines of reasoning may account for the differences between our results and previously reported satellite-based results of the sensitivity of SGS and EGS dates to elevation. First, satellite sensors observe average plant phenology within a square with sides of 500-1000    [14] in Table 2 Beginning of flowering  [14] in Table 2 Beginning of flowering  [14] in Table 2 Beginning of flowering  [14] in Table 2 Beginning of flowering  [14] in Table 2 First flowers open Germany Goat willow In situ observation 1971 −2000 2.51 ± 0.16 day (100 m) −1 [14] in Table 2 Full flowering Switzerland Black elder In situ observation 1971 −2000 3.12 ± 0.21 day (100 m) −1 [14] in Table 2 First flowers open Austria Rye In situ observation 1971 −2000 3.22 ± 0.39 day (100 m) −1 [14] in Table 2 First flowers open Germany Small-leaved lime In situ observation 1971 −2000 4.03 ± 0.24 day (100 m) −1 [14] in Table 2 Full flowering Switzerland Common grape wine In situ observation 1971 −2000 1.13 ± 0.26 day (100 m) −1 [14] in Table 2 SGS m. Such a square may include many tree species characterized by a wide range of the timing and patterns of leaf flush and leaf fall. Second, the SGS and EGS dates defined by satellite observations may differ from those defined by in situ observations. Third, data are missing from a time-series of satellite observations when there is cloud contamination and atmospheric noise. Finally, the characteristics of microtopography cannot be detected by satellite sensors with a spatial resolution of 500-1000 m. Miura et al. [67] showed the utility of analysing a time-series of an index of vegetation based on observations by the AHI (Advanced Himawari Imager) sensor mounted on the geostationary Himawari-8 satellite for phenological observations at TKY. Despite its coarse spatial resolution (1000 m), the AHI sensor has a very high temporal resolution (every 2.5 min around Japan). Important tasks for the future include integrative evaluation of calculations of SGS and EGS dates with a degree-day model, analysis of SENTINEL-2A/2B-satellite observations with a fine spatial resolution but low observation frequency, and analysis of Himawari-8 satellite observations with a high observation frequency but coarse spatial resolution.

Interannual variability of modelled SGS and EGS dates
The difference between modelled SGS dates at Takayama-AMeDAS and TKY-tower decreased by 0.32 day decade −1 ( Fig 2C). Several different mechanisms may have accounted for this pattern. First, the daily mean temperature between March and May increased ( Fig 4A). As a result, the date when the daily mean temperature exceeded the threshold temperature for CET SGS (T t,SGS = 2˚C) advanced, especially at TKY-tower. Second, the date when the CET reached CET SGS (= 255.4˚C) also advanced. Ziello et al. [14] reported advances of 0.065 ± 0.028 day year −1 per 100 m for the full flowering of the common alder and 0.049 ± 0.020 day year −1 per 100 m for the beginning of the flowering of Norway spruce, and delay of 0.025 ± 0.011 day year −1 per 100 m for the opening of the first flower of black elder in Europe. Our result of advance of 0.0037 day year −1 per 100 m for the SGS was smaller than those of common alder and Norway spruce. In contrast, there was no significant long-term trend of the difference between the modelled EGS dates at Takayama-AMeDAS and TKYtower ( Fig 2D). This lack of a trend may have resulted from the fact that the daily mean temperature from September to November increased ( Fig 4C). The date when the daily mean temperature was below the threshold temperature for CET EGS (T t,EGS = 18˚C) at those sites was therefore delayed. Because of this delay, the date when the CET reached CET EGS (= −375.1˚C) was also delayed. The difference between the Takayama-AMeDAS and TKY-tower results was larger for the modelled SGS dates than for the modelled EGS dates (Fig 2C and 2D). This disparity may have been caused by the relationship between the interannual variability of the SGS and EGS dates and that of the monthly mean temperature. The interannual variability of the modelled SGS date correlated with that of monthly mean temperature at TKY-tower but not at Takayama-AMeDAS (Table 2). In contrast, at both sites, the interannual variability of the modelled EGS dates correlated with that of monthly mean temperature ( Table 2). These relationships suggest that earlier SGS dates in the lower reaches of a river basin may not be synchronized with the earlier SGS dates in the upper reaches of the same basin.

Sensitivity of temperature to elevation
We used a temperature lapse rate of 0.60˚C (100 m) −1 , which has been commonly used in several previous studies [15,26,51,52]. For comparison, Ueno et al. [68] reported temperature lapse rates in the mountainous region of central Japan, which includes the steep river basin targeted in this study, of 0.59˚C (100 m) −1 during March-May, 0.55˚C (100 m) −1 during June-August, 0.55˚C (100 m) −1 during September-November, and 0.58˚C (100 m) −1 during December-February. In contrast, studies in Europe indicated that the temperature lapse rates in the Pyrenees and Slovakia were only 0.44-0.51˚C (100 m) −1 in 2005 and 2006 [11], but they averaged 0.59˚C (100 m) −1 between 1981 and 2010 [23].
To assess the accuracy of the temperature lapse rate of 0.6˚C (100 m) −1 , we examined the differences between the observed daily mean temperature at TKY-weather station (1342 m a.s. l.) and that estimated from the observed daily mean temperature at Takayama-AMeDAS (560 m a.s.l.) and a temperature lapse rate of 0.6˚C (100 m) −1 . The differences displayed a seasonal pattern (Fig 9; Table 5). The differences (estimated − observed) were −0.62˚C during January-March and −1.25˚C during October-December. The temperatures estimated with a lapse rate  Table 5.
https://doi.org/10.1371/journal.pone.0255078.g009 of 0.6˚C (100 m) −1 therefore underestimated the observed temperatures at TKY-weather station. However, the differences in April, May, and September, which are important times to calculate SGS and EGS dates, were comparatively small and ranged from −0.03 to −0.40˚C. Furthermore, the differences between June and August were positive and ranged from 0.41 to 0.55˚C. The seasonality of the temperature lapse rates is primarily a result of the seasonality of humidity [30,[68][69][70]. The small differences between the estimated and observed daily mean temperatures (0.03-0.40˚C) in April, May, and September may have caused errors of a few days in the modelled SGS and EGS dates. Consideration should also be given to the effect of temperature inversions on the SGS date in a river basin [30].
These observations suggest that the temperature lapse rate of 0.6˚C (100 m) −1 may slightly underestimate the temperature in high-elevation areas, except during summer. However, images of phenology taken on modelled SGS and EGS dates in each year at TKY-tower showed almost the same condition of the canopy surface, although there were small variations in leaf colouring and leaf fall (Fig 5). This result suggests that interannual variability of SGS and EGS dates at a site can be evaluated in a relative sense with a constant lapse rate. Important tasks for the future will include meteorological observations in mountainous regions (i.e., greater than 1000 m a.s.l.) where there are few weather stations such as the AMeDAS [13] and examination of the temperature lapse rate in each river basin.

Further acquisition of ground truth data
We used phenology and live camera images taken in the upper and lower reaches of the Daihachiga River basin as validation data (Figs 5 and 6). However, evaluation of the uncertainty of the spatial representation of phenological events by our model has been constrained by the amount of data available to us. The development of degree-day models of phenology to estimate SGS and EGS dates will require many ground truth observations of phenology at multiple points along an elevational gradient. In particular, reducing the uncertainties in the areal coverages of different tree species captured in phenological and live-camera images will require many observations at multiple points. We have recently been able to easily acquire live camera images taken at multiple tourist sites and riversides. Silva et al. [71] reported that phenological information can be extracted from analysis of images uploaded on Twitter. However, those images are poorly suited for detecting plant phenology if the percentage of vegetation in them is low, and the purpose of capturing the images is completely different (e.g., weather monitoring, tourism, disaster prevention, and recreation). Such images may therefore be difficult to use for quantitative analyses (e.g., extraction of RGB values from daily images of phenology; e.g., [31,72,73]). In addition, it is difficult to obtain such images in forests and mountainous regions, where there are few requirements related to traffic and societal needs (e.g., tourism and disaster prevention).
Websites on the Internet that can be easily accessed for uploading images with geotags (e.g., Mapillary at https://www.mapillary.com, accessed 6 July 2021) can be used to upload and download images with properties such as GPS (Global Positioning System) coordinates, date, time, and contributor. This capability suggests that images taken along roads within a river basin and periodically archived could be used to evaluate the spatiotemporal variability of SGS and EGS dates. In fact, Ishihara et al. [29] used daily movie-camera images taken along roads in the Daihachiga River basin to document phenology. One of us (K.N.N.) has already used a drive recorder to capture images of the Daihachiga River basin and has uploaded the images with GPS logs to Mapillary. Important tasks for the future include comparing periodically archived ground-truth images, data recorded by SENTINEL-2A/2B-satellites, and modelled SGS and EGS dates.

Conclusions
We used a degree-day phenology model to evaluate the spatiotemporal variability of the SGS and EGS dates in the Daihachiga River Basin on a century timescale and the relationship of that variability to changes of temperature. We found that (1) the sensitivity of the modelled SGS and EGS dates to elevation changed from 3.29 days (100 m) −1 (−5.48 days˚C −1 ) and −2.89 days (100 m) −1 (4.81 days˚C −1 ), respectively, in 1900 to 2.85 days (100 m) −1 (−4.75 days C −1 ) and −2.84 day (100 m) −1 (4.73 day˚C −1 ), respectively, in 2019, and (2) the long-term trend of the sensitivity of the modelled SGS date to elevation was −0.0037 day year −1 per 100 m, but the analogous trend in the case of the modelled EGS date was not significant. Despite the use of daily phenology and live camera images as well as satellite data with a high spatial resolution to assess the generality and representativeness of the modelled SGS and EGS dates, limitations of our study still remain. The next steps will require (1) examination of the relationship between elevation and SGS and EGS dates in other steep river basins where annual mean temperatures are much higher than the temperatures in the Daihachiga River basin, (2) integrative evaluation of the SGS and EGS dates based on a degree-day phenological model and satellite observations with a fine spatial resolution and high observation frequency, (3) acquisition of an accurate basin-specific temperature lapse rate in each target basin, and (4) further acquisition of ground truth data such as live camera images and images with geotags on the Internet. Plant phenology was more sensitive to elevation than to latitude in a local area (e.g., within 2˚latitude × 3˚longitude) [9], where we could ignore the effect of phenological plasticity associated with climate change. Hence, the development of degree-day phenological models in multiple steep river basins with different climate conditions will deepen our ecological understanding of the sensitivity of spring and autumn phenology to future climate change.