Shifts in Plant Phenology and Its Responses to Climate Warming in Three Temperate Cities of China during 1963–2020

Abstract

The advance of spring phenology and the delay of autumn phenology caused by global warming have been documented by many studies. However, most research has focused on natural areas, with limited studies conducted on phenological observations in urban environments. Here, we selected the first flowering date (FFD), first leaf date (FLD), and leaf coloring date (LCD) at three sites (Beijing, Harbin, and Mudanjiang) from the China Phenological Observation Network. We analyzed the phenological changes of 84 species between 1963–1991 and 1992–2020 to examine their response to urban warming. We then quantified the correlations and regressions between phenological events and preseason temperature. The results show the following: (1) Among the three sites, the mean FFD and FLD were earliest in Beijing, while the mean LCD occurred earliest in Harbin and latest in Beijing. (2) FFD and FLD exhibited a significant trend towards earlier occurrences at all three sites, while LCD showed a significant delay trend except for the Mudanjiang site. Specifically, at the Beijing, Harbin, and Mudanjiang sites, the mean FFD advanced by 8.32 days, 6.11 days, and 2.60 days in the latter period (p < 0.05), whereas the mean FLD advanced by 11.30 days, 7.21 days, and 5.02 days (p < 0.05), respectively. (3) In Beijing, Harbin, and Mudanjiang, both FFD and FLD were significantly negatively correlated with preseason temperature. However, no consistent relationship was observed between LCD and preseason temperature. These results enhance our understanding of the response of plant phenology to urban warming.

1. Introduction

Vegetation phenology is a sensitive indicator of climate change [1,2]. Phenological changes can alter ecosystem structure and function by leading to phenological mismatches of interacting species [3], helping the success of invasive plant species [4], and increasing the effective duration of plant assimilation [5]. Meanwhile, phenological changes also provide feedback to the climate system by influencing the water, energy, and carbon cycle [2,6]. Understanding of the phenological changes under historical climate variations could provide insights into how ecosystems might respond to future climate change. Various environmental factors have been identified as influencing vegetation phenology, including temperature [7], photoperiod [8], moisture [9], and CO₂ [10]. Among all these factors, temperature reveals a primary control. Before initiating spring phenological events, plants first experience a period of endodormancy that needs to be broken by cold temperature and then ecodormancy that is accelerated by warm temperature [11]. Elevated autumn temperature was also found to slow down the senescence rate because plants receive more time for photosynthesis and nutrient reallocation under a suitable environment [12,13].

Phenological events are critical timings in plant developmental stages and are sensitive to climate change [6]. In the context of global warming, advanced spring phenological events (such as the leaf-out and flowering date) and delayed autumn phenological events (such as leaf coloring and leaf fall) have been generally reported by previous studies using ground-based phenology observation data across continents [7,14,15,16]. For spring phenological events, a comprehensive study based on more than 125,000 data series from phenological networks in 21 European countries during 1971–2000 demonstrated that the average advance of spring (summer) phenology was 2.5 d/decade [7]. In East Asia, the start of the growing season for Ginkgo biloba L. was found to occur 4 days earlier during 1953–2000 according to the phenological observation data [14]. For autumn phenological events, the leaf senescence dates were delayed by 2.6 d/decade in 1982–2011 on the basis of a meta-analysis in China [15]. Similarly, in the USA, leaf coloring was also found to be delayed by an average of 3.6 d/decade [17]. However, these observations and results were mainly derived from natural ecosystems. Considering that urban warming tends to be stronger than natural vegetation, phenological changes in cities may be more significant than those in natural ecosystems. For example, the length of the growing season in urban core areas was found to be increased by 9.9 ± 0.77 days compared to rural counterparts according to a satellite-based study in 343 Chinese cities [18]. Although there have been studies using remote sensing data to investigate the trends of phenological metrics in major cities of China [19,20], in situ observational studies in cities are still lacking.

In our study, we selected three phenological observation sites located within the urban area of Beijing, Harbin, and Mudanjiang as the study area to investigate the temporal changes of plant phenology and its response to urban warming. These sites record a long and continuous time series of first leaf dates (FLD), first flowering dates (FFD), and leaf coloring dates (LCD) from 1963 to 2020. We first divided the 1963–2020 period into two 29-year periods (1963–1991 and 1992–2020), and calculated the average FLD, FFD, and LCD for each period, as well as the differences between the two 29-year periods. Then, we conducted the correlation and regression analysis between phenophases (FLD, FFD, and LCD) and preseason temperature. In summary, our aims were (1) to identify the distinctions between the phenophases (FLD, FFD, and LCD) in different time periods (1963–1991 and 1992–2020); (2) to quantify the response of phenophases (FLD, FFD, and LCD) to preseason temperature.

2. Data and Methods

2.1. Study Site

We selected Beijing (40°2′ N, 116°33′ E), Harbin (45°70′ N, 126°63′ E), and Mudanjiang (44°43′ N, 129°67′ E) as the study area (Figure 1a) according to the principle of longer observation and better data continuity. The three stations have a temperate monsoon climate with a hot, rainy summer and a cold, dry winter. The coldest month is January, and the hottest month is July. The average temperature from 1963 to 2020 is 12.5 °C, 4.5 °C, and 4.4 °C for Beijing, Harbin, and Mudanjiang, respectively (Figure 1b). The average annual precipitation is 546 mm, 517 mm, and 533 mm for Beijing, Harbin, and Mudanjiang, respectively. The main vegetation type is deciduous broad-leaved forest.

2.2. Data Sources

All phenological data for this study are from the China Phenological Observation Network (CPON), which was established in 1963. We investigated three phenophases, including the first flowering date (FFD), first leaf date (FLD), and leaf coloring date (LCD). Observations of each phenophase conform to uniform observation criteria and guidelines [21]. FFD was defined as the date when a fixed individual has at least 3 flowers fully opens. FLD was defined as the date when a fixed individual with at least 3 leaves fully opens. LCD was defined as the date when an individual shows yellow leaves on 95%–100% of its crowns. FFD and FLD are considered representative indicators of spring phenology, while LCD represents autumn phenology [22,23]. Additionally, FLD and LCD can be interpreted as the start and end dates of the photosynthesis period, respectively [23,24]. The phenological data span from 1963 to 2020, and the effective observation years ranged between 36 and 52 years among sites because there are a few missing data in certain years due to discontinuous observations (

Table A1. Summary of phenological data from Beijing stations investigated in this study. N_FFD1, N_FLD1, N_LCD1, number of observation years during the period 1963–1991; N_FFD2, N_FLD2, N_LCD2, number of observation years during the period 1992–2020; MFLD, mean timing of first leaf dates; MLCD, mean timing of leaf coloring date; MFFD, mean timing of first flowering date (all from 1963 to 2020); LF, life form; T, tree species; S, shrub species; V, vine species. The slash denotes that the data are less than 10 years old and have been excluded.

	Species	Family	LF	NFFD1	NFFD2	MFFD	NFLD1	NFLD2	MFLD	NLCD1	NLCD2	MLCD
1	Fraxinus chinensis	Fraxinus	T	18	24	4/19	24	26	4/14	22	25	10/19
2	Castanea mollissima	Castanea	T	26	25	5/30	26	26	4/19	25	25	11/2
3	Syringa reticulata subsp. amurensis	Syringa	T	22	19	5/18	24	19	4/2	21	18	11/5
4	Prunus persica ‘Duplex’	Prunus	T	18	21	4/8	17	20	4/11	12	16	11/3
5	Platycladus orientalis	Platycladus	T	22	23	3/28	-	-	-	-	-	-
6	Ailanthus altissima	Ailanthus	T	24	20	5/24	25	22	4/18	17	21	10/22
7	Robinia pseudoacacia	Robinia	T	25	26	5/2	25	26	4/14	23	23	10/30
8	Pyrus betulifolia	Pyrus	T	17	21	4/14	17	19	4/8	16	20	10/28
9	Juniperus chinensis	Juniperus	T	23	21	4/1	-	-	-	-	-	-
10	Salix matsudana	Salix	T	26	19	4/2	24	19	3/29	20	18	11/12
11	Albizia julibrissin	Albizia	T	26	18	6/9	26	18	5/1	22	14	11/1
12	Styphnolobium japonicum	Styphnolobium	T	23	20	7/14	23	20	4/25	18	13	11/14
13	Rosa xanthina	Rosa	S	14	20	4/23	14	21	4/2	10	18	11/10
14	Cotinus coggygria var. cinereus	Cotinus	T	17	19	4/21	16	18	4/12	17	19	10/25
15	Rhodotypos scandens	Rhodotypos	T	17	18	4/16	17	17	4/7	16	11	11/13
16	Populus × canadensis	Populus	T	23	26	3/26	25	24	4/7	23	21	10/31
17	Diospyros lotus	Diospyros	T	22	25	5/17	23	25	4/17	16	20	10/28
18	Forsythia suspensa	Forsythia	S	18	21	3/28	18	19	4/8	15	19	11/4
19	Styphnolobium japonicum ‘Pendula’	Styphnolobium	T	16	19	7/11	18	21	4/13	-	-	-
20	Koelreuteria paniculata	Koelreuteria	T	16	19	5/31	19	21	4/11	17	20	10/27
21	Populus tomentosa	Populus	T	22	26	3/17	21	26	4/9	14	18	11/6
22	Paeonia × suffruticosa	Paeonia	S	26	25	4/20	25	26	3/31	19	23	10/29
23	Paulownia fortunei	Paulownia	T	16	16	4/20	15	12	4/23	-	-	-
24	Catalpa bungei	Catalpa	T	18	18	4/27	19	19	4/12	-	-	-
25	Prunus × yedoensis	Prunus	T	26	20	4/14	26	19	4/11	25	19	10/31
26	Morus alba	Morus	S	24	23	4/22	24	23	4/20	21	21	10/30
27	Prunus davidiana	Prunus	T	27	26	3/22	27	26	4/1	18	19	11/6
28	Punica granatum	Punica	T	16	20	5/27	17	21	4/23	13	19	10/30
29	Diospyros kaki	Diospyros	T	23	23	5/14	23	23	4/15	21	22	10/29
30	Metasequoia glyptostroboides	Metasequoia	T	-	-	-	21	20	4/11	19	20	11/3
31	Viburnum farreri	Viburnum	T	16	11	3/30	14	10	3/29	-	-	-
32	Salix matsudana ‘Pendula’	Salix	T	24	20	4/6	26	21	3/30	24	20	11/13
33	Firmiana simplex	Firmiana	T	23	23	6/20	24	25	4/27	21	24	10/29
34	Malus × micromalus	Malus	T	26	21	4/13	26	21	3/31	20	17	11/11
35	Populus simonii	Populus	T	13	15	4/3	14	15	4/6	11	14	11/4
36	Prunus armeniaca	Prunus	T	26	26	4/2	24	25	4/13	22	23	11/6
37	Ginkgo biloba	Ginkgo	T	-	-	-	16	20	4/10	16	20	10/30
38	Jasminum nudiflorum	Jasminum	T	10	19	3/17	-	-	-	-	-	-
39	Pinus tabuliformis	Pinus	T	25	21	4/26	25	22	4/26	-	-	-
40	Ulmus pumila	Ulmus	T	25	18	3/16	26	18	4/8	19	18	11/11
41	Prunus triloba	Prunus	T	21	26	4/6	20	25	4/8	20	24	10/30
42	Yulania denudata	Yulania	T	26	21	3/30	26	21	4/12	25	19	11/1
43	Sorbaria sorbifolia	Sorbaria	S	15	18	6/8	16	20	3/30	14	14	11/8
44	Syringa oblata	Syringa	S	24	26	4/10	24	26	3/31	20	21	11/3
45	Cercis chinensis	Cercis	S	24	26	4/14	23	26	4/14	22	25	10/28
46	Wisteria sinensis	Wisteria	V	24	26	4/21	24	26	4/21	22	23	11/8
47	Lagerstroemia indica	Lagerstroemia	S	22	25	7/5	24	25	4/26	21	24	10/28
48	Yulania liliiflora	Yulania	S	26	20	4/4	25	20	4/11	24	18	10/29
49	Prunus × subhirtella	Prunus	T	16	16	4/10	15	15	4/11	15	15	10/30

,

Table A2. Summary of phenological data from Harbin stations investigated in this study. N_FFD1, N_FLD1, N_LCD1, number of observation years during the period 1963–1991; N_FFD2, N_FLD2, N_LCD2, number of observation years during the period 1992–2020; MFLD, mean timing of first leaf dates; MLCD, mean timing of leaf coloring date; MFFD, mean timing of first flowering date (all from 1963 to 2020); LF, life form; T, tree species; S, shrub species; V, vine species. The slash denotes that the data are less than 10 years old and have been excluded.

	Species	Family	LF	NFFD1	NFFD2	MFFD	NFLD1	NFLD2	MFLD	NLCD1	NLCD2	MLCD
1	Syringa reticulata subsp. amurensis	Syringa	T	12	18	6/7	12	18	4/23	12	10	9/10
2	Forsythia mandschurica	Forsythia	S	11	18	4/19	11	18	5/2	11	12	9/28
3	Philadelphus schrenkii	Philadelphus	S	14	18	5/28	14	18	4/29	12	10	10/5
4	Pinus koraiensis	Pinus	T	-	-	-	10	17	6/1	-	-	-
5	Juglans mandshurica	Juglans	T	17	17	5/18	17	17	5/10	-	-	-
6	Viburnum opulus subsp. calvescens	Viburnum	S	13	18	5/25	13	18	4/28	13	17	10/8
7	Lonicera maackii	Lonicera	S	16	18	5/24	16	18	4/30	-	-	-
8	Viburnum burejaeticum	Viburnum	S	13	18	5/17	13	18	4/26	13	13	10/9
9	Morus alba	Morus	S	11	18	5/20	12	18	5/18	12	15	9/27
10	Prunus sibirica	Prunus	T	15	18	5/1	15	18	5/10	-	-	-
11	Caragana arborescens	Caragana	S	12	18	5/19	12	18	5/8	-	-	-
12	Fraxinus mandshurica	Fraxinus	T	11	18	5/3	14	18	5/13	14	15	9/29
13	Xanthoceras sorbifolium	Xanthoceras	T	11	18	5/19	11	18	5/17	-	-	-
14	Flueggea suffruticosa	Flueggea	S	10	17	6/19	10	18	5/18	10	16	10/2
15	Ulmus pumila	Ulmus	T	17	18	4/16	17	18	5/3	16	12	10/8
16	Prunus triloba	Prunus	T	16	18	5/2	16	18	5/7	14	14	10/2
17	Pinus sylvestris var. mongolica	Pinus	T	-	-	-	10	16	5/28	-	-	-
18	Sorbaria sorbifolia	Sorbaria	S	15	16	7/9	17	18	4/18	-	-	-
19	Catalpa ovata	Catalpa	T	12	16	6/24	13	18	5/21	-	-	-
20	Syringa oblata	Syringa	S	18	18	5/8	18	18	5/2	14	11	10/13

and

Table A3. Summary of phenological data from Mudanjiang stations investigated in this study. N_FFD1, N_FLD1, N_LCD1, number of observation years during the period 1963–1991; N_FFD2, N_FLD2, N_LCD2, number of observation years during the period 1992–2020; MFLD, mean timing of first leaf dates; MLCD, mean timing of leaf coloring date; MFFD, mean timing of first flowering date (all from 1963 to 2020); LF, life form; T, tree species; S, shrub species; V, vine species. The slash denotes that the data are less than 10 years old and have been excluded.

	Species	Family	LF	NFFD1	NFFD2	MFFD	NFLD1	NFLD2	MFLD	NLCD1	NLCD2	MLCD
1	Acer tataricum subsp. ginnala	Acer	T	15	28	5/24	14	27	5/2	14	29	9/24
2	Prunus padus	Prunus	T	15	28	5/7	14	28	4/22	11	28	9/28
3	Salix babylonica	Salix	T	13	29	4/30	13	27	4/23	11	26	10/21
4	Picea koraiensis	Picea	T	-	-	-	11	28	5/8	-	-	-
5	Lespedeza bicolor	Lespedeza	T	14	28	7/6	-	-	-	10	23	9/20
6	Styphnolobium japonicum	Styphnolobium	T	10	10	7/7	11	17	5/17	10	11	9/24
7	Phellodendron amurense	Phellodendron	T	11	29	6/7	-	-	-	-	-	-
8	Prunus salicina	Prunus	T	16	27	5/6	12	25	5/4	11	18	10/3
9	Spiraea salicifolia	Spiraea	S	13	22	6/18	12	18	4/21	-	-	-
10	Rubus komarovii	Rubus	S	11	26	6/4	11	27	4/30	-	-	-
11	Malus baccata	Malus	T	15	28	5/12	14	26	4/24	11	22	10/3
12	Crataegus pinnatifida var. major	Crataegus	T	16	28	5/27	15	28	4/29	11	26	9/28
13	Rhamnus davurica	Rhamnus	S	11	27	5/26	10	25	5/1	-	-	-
14	Caragana arborescens	Caragana	S	11	28	5/16	-	-	-	-	-	-
15	Acer saccharum	Acer	T	12	29	4/26	11	28	4/27	-	-	-
16	Euonymus alatus	Euonymus	S	14	29	5/19	13	28	4/29	11	19	9/29
17	Populus pseudosimonii	Populus	T	10	24	4/27	10	25	5/2	-	-	-
18	Populus simonii	Populus	T	14	20	4/30	14	26	5/5	-	-	-
19	Rhododendron dauricum	Rhododendron	S	10	29	4/29	-	-	-	-	-	-
20	Larix gmelinii	Larix	T	11	29	5/11	11	29	4/26	10	29	10/10
21	Prunus pseudocerasus	Prunus	T	15	20	5/4	13	18	5/1	11	16	10/12
22	Ulmus pumila	Ulmus	T	12	25	4/22	11	21	5/8	-	-	-
23	Prunus triloba	Prunus	T	13	12	5/3	12	11	4/30	11	12	10/5
24	Quercus dentata	Quercus	T	14	23	5/13	14	23	5/3	12	19	9/30
25	Sorbaria sorbifolia	Sorbaria	S	12	18	7/16	12	19	4/13	-	-	-
26	Corylus heterophylla	Corylus	S	-	-	-	10	28	5/1	-	-	-
27	Catalpa ovata	Catalpa	T	11	27	6/30	11	27	5/16	-	-	-
28	Tilia amurensis	Tilia	T	11	28	6/28	-	-	-	-	-	-

). A total of 84 species (belonging to 58 genera), including 67 species of trees, 27 species of shrubs, and 1 vine, were analyzed in three cities.

The climate data, including daily mean temperatures and precipitation, for the period from 1963 to 2020 were sourced from the China Meteorological Data Sharing Service System (https://data.cma.cn/, accessed on 27 September 2024).

2.3. Methods

2.3.1. Statistical Method of Phenological Changes

In order to analyze the impact of climate change on phenology at three sites over the past 60 years, the study period was divided into two periods (1963–1991 and 1992–2020) with the same duration. A comparison of temperature data from the Beijing, Harbin, and Mudanjiang sites between the two periods revealed that the mean monthly average temperatures at all three locations during 1992–2020 were higher than those during 1963–1991 by 1.53 °C, 1.43 °C, and 1.23 °C, respectively. The number of months with significant differences (p < 0.05) was eleven in Beijing and ten for both Harbin and Mudanjiang. Harbin and Mudanjiang experienced higher precipitation during 1992–2020 compared to that in the earlier period by an increase of 2.64 mm and 11.85 mm, respectively. However, Beijing station witnessed a decrease in precipitation by about 5.04 mm during the later period as opposed to its earlier period. The number of months with significant differences (p < 0.05) in monthly precipitation was three for Beijing station, five for Harbin station, and six for Mudanjiang (Figure 2). The significant difference in the monthly mean temperatures between the two periods provides a basis for analyzing the response of phenology to global warming.

The phenophase is represented by the day of the year (DOY). We only analyzed the species whose cumulative observation years are greater than or equal to 10 years. The average FFD, FLD, and LCD of the two periods and the whole period (1963–2020) were counted, the deviation between the two periods (later period minus earlier period) was analyzed, and the frequency distribution diagram and box plot were drawn. An unequal variance t-test was performed to test whether the deviations differed significantly from zero.

2.3.2. Analysis of the Relationship between Phenology and Preseason Temperature

Plant phenology is influenced by fluctuations in mean temperature during the period preceding phenological events [25]. Numerous studies have demonstrated a significant correlation between plant phenology and the mean temperature of the month in which the plant’s multi-year mean phenological onset occurs, as well as the two preceding months [23,26]. However, the optimum period (OP) for temperature influence cannot be limited to a single day. It is generally accepted that plant phenology is affected by temperatures from at least two weeks prior. Therefore, in this study, the OP was calculated using temperatures from 15 to 120 days before the phenological event. The OP was utilized in this study to examine the correlation between phenology and temperature, as defined by [14,23]

OP = [BP,EP],

(1)

where EP is the end date of OP (in DOY), defined as the average date of the phenophase from 1963 to 2020. BP is the beginning date of the period (in DOY). Then, the Pearson’s correlation coefficient (R) between the phenophase (FFD, FLD, or LCD) and average temperature during [BP, EP] was calculated by moving each BP from EP-15 to EP-120 days by a step length of 1 day. The period [BP, EP] with the highest absolute value of correlation coefficient was designated as the OP. The BP of LCD was restricted to occur after the summer solstice [13]. The regression slope for each phenophase (FFD, FLD, LCD) was calculated in relation to the corresponding average temperature within the OP, allowing for the quantification of the temperature sensitivity of each phenophase. The Adjusted R² and RMSE for each phenophase (FFD, FLD, LCD) was calculated in relation to the corresponding average temperature within the OP using the Sigmoidal model (Equation (2), Figure A1) and Linear model (Equation (3), Figure A2). The significance of the correlation coefficients, regression slopes, Adjusted R², and RMSE values were assessed using a two-sample t-test.

$$f(x)={\displaystyle \frac{c}{1+e^{(a+bx)}}}+d$$

(2)

where c and d represent the horizontal asymptotes, while b is the growth rate parameter. The inflection point of the curve is located at x = −a/b.

$$\mathrm{y}=i\mathrm{x}+j$$

(3)

where i is the slope, and j is the y-intercept.

3. Results

3.1. Phenological Variation across Sites

There are significant differences in the mean FFD, FLD, and LCD of each site from 1963 to 2020 (Figure 3). The FFD range across three sites spanned 75 to 198 days, and more than 50% of the species (52.75%) were distributed in DOY 105–150. The FFD of the Beijing site was 24.77 days and 26.58 days earlier than Harbin and Mudanjiang. The FLD range was 75–197 days, and over 80% of species (80.90%) were distributed in DOY 89–125. The mean FLD was earliest in Beijing, followed by Mudanjiang, and the latest in Harbin. The mean LCD of Harbin was slightly earlier than that of Mudanjiang, and the mean LCD of Beijing was 31.14 days later than that of Harbin.

3.2. Phenological Change between 1963–1991 and 1992–2020

Based on the phenological differences observed across all species between the two periods (1963–1991 and 1992–2020), both FFD and FLD advanced at all three sites, while LCD was delayed at all sites except for Mudanjiang (Figure 4). Specifically, the median FFD at Beijing, Harbin, and Mudanjiang in 1992–2020 was earlier by 8.10 days, 2.99 days, and 2.18 days compared to 1963–1991, respectively. The median FLD in 1992–2020 was earlier by 11.15 days, 7.54 days, and 1.15 days, respectively, at these sites compared to the period 1963–1991. The median LCD at Beijing and Harbin in 1992–2020 was later by 5.10 and 7.05 days, respectively, compared to 1963–1991, while it was earlier by 0.17 days at Mudanjiang.

Changes in the frequency distribution of the three phenological phases from 1992–2020 and 1963–1991 reveal an overall advancement in FFD (Figure 5) and FLD (Figure 6), while a delay was observed in LCD at Beijing, Harbin, and Mudanjiang (Figure 7). Specifically, Beijing (85.11%) recorded the highest proportion of species with significantly advanced FFD, followed by Mudanjiang (57.69%), and the lowest proportion was found in Harbin (22.22%). The magnitude of advancement followed the order of Beijing (8.32 days), Harbin (6.11 days), and Mudanjiang (2.60 days). Similarly, for FLD, Beijing (2.60 days) had the highest proportion of significantly advanced species, followed by Harbin (75.00%), and Mudanjiang (8.70%) had the lowest. The magnitude of advancement also followed the same order of Beijing (11.30 days), Harbin (7.21 days), and Mudanjiang (5.02 days). In contrast, Harbin (81.82%) had the highest proportion of species with significantly delayed LCD, followed by Beijing (46.34%), while Mudanjiang showed no clear pattern of change. The magnitude of delay in LCD was greatest in Harbin (11.66 days), followed by Beijing (8.43 days) and Mudanjiang (0.92 days). These findings further underscore the notable advancing trend in spring and summer phenophases and the delaying trend in autumn phenophases.

3.3. Relationship between Phenophases and Temperature

Figure 8 shows the OP corresponding to vegetation phenology at different sites, revealing that FLD and FFD are generally significantly correlated with OP temperature (p < 0.05). In contrast, the correlation between LCD and OP temperature is not significant for most species (Figure 6). The OP duration of FFD ranged from 17 days to 120 days across the three sites, with the mean OP duration at Harbin (50.81 days), Beijing (51.55 days), and Mudanjiang (59.10 days) in ascending order. For FLD, the OP duration varied from 16 days to 120 days across all sites, with a mean OP of 42.43 days, 42.98 days, and 51.20 days at Mudanjiang, Beijing, and Harbin, respectively. For LCD, the OP duration ranged from 16 days to 111 days, with a mean OP duration of 38.33 days, 46.00 days, and 46.80 days at Harbin, Mudanjiang, and Beijing, respectively.

A comparison of the Adjusted R² values between the Sigmoidal and Linear model revealed that, with the exception of the significant difference (p < 0.05) in Adjusted R² for FFD and OP temperatures at the Mudanjiang site, the differences in Adjusted R² between the two models were not statistically significant in other cases (Figure A3). Additionally, no significant differences were observed between the RMSEs derived from the Linear and Sigmoidal models. Notably, the RMSEs for all phenological periods (FFD, FLD, and LCD) in relation to OP temperatures were consistently smaller for the Linear model across all three sites (Figure A4). Consequently, the results of the Linear model were employed in subsequent analyses.

An analysis of the correlation between vegetation phenology and OP temperature revealed that at Beijing, Harbin, and Mudanjiang, FFD and FLD exhibited negative correlations with OP temperature, whereas the correlation between LCD and OP temperature varied significantly among sites (Figure 7). Specifically, the mean correlation coefficients between FFD and OP temperature were −0.76 ± 0.23 (mean ± SD), −0.69 ± 0.22, and −0.53 ± 0.27, respectively, with 89.36%, 88.89%, and 38.46% of species showing significance at the 5% level. For FLD, the mean correlation coefficients were −0.79 ± 0.13, −0.71 ± 0.12, and −0.81 ± 0.10, with 97.83%, 100.00%, and 100.00% of species significant at a 0.05 level. In contrast, the mean correlation coefficients between LCD and OP temperature were 0.21 ± 0.18, 0.49 ± 0.22, and −0.43 ± 0.19, respectively, with only 12.20%, 54.55%, and 15.38% of species displaying a significant correlation coefficient (p < 0.05).

The temperature sensitivity (regression slope) of vegetation phenology was significantly different among sites and species (Figure 9). Significant advancements predominantly characterized the responses of FFD and FLD to OP temperatures. At the Beijing, Harbin, and Mudanjiang sites, the temperature sensitivity of FFD ranged from −0.52 to −6.99 days °C⁻¹, −1.94 to −7.17 days °C⁻¹, and −2.33 to −18.02 days °C⁻¹, respectively. For FLD, the temperature sensitivity ranged from −1.44 to −6.99 days °C⁻¹, −2.42 to −4.47 days °C⁻¹, and −2.44 to −3.95 days °C⁻¹, respectively. In contrast, the response of LCD to OP temperatures did not exhibit a consistent pattern among sites. At Beijing, Harbin, and Mudanjiang, the proportions of species with positive slopes (indicating advancement with increasing OP temperatures) were 71.43%, 91.67%, and 28.57%, respectively. The magnitudes of these responses ranged from 0.0037 to 4.16, 0.73 to 6.36, and 0.57 to 1.35 days °C⁻¹, respectively. For the remaining species, LCD advanced by 0.35–2.31 days, 2.79–2.79 days, and 0.85–7.50 days, respectively, with a 1 °C increase in OP temperatures. Overall, the mean temperature sensitivity across all species at the three sites was −5.90, −3.28, and 0.35 days °C⁻¹ for FFD, FLD, and LCD, respectively.

4. Discussion

In terms of spatial pattern, our results support the bioclimatic law proposed by Hopkins, which states that spring phenology would be delayed but autumn phenology would be advanced by latitude north [27]. Concretely, the median FLD and FFD in Beijing were earlier than in Mudanjiang, followed by Harbin in both observation periods. This latitudinal gradient of spring phenology might be attributed to the response of plants to varying environmental conditions. However, the gradient of FLD and FFD across latitudes between Beijing and Harbin (4.36 days for FFD and 4.53 days for FLD per 1° N) was higher than Hopkins’ hypothesis (4 days per 1° N). Compared to the other cities of China (2.46 to 3.48 days per 1° N), North America (2.62 days per 1° N), or the whole Northern Hemisphere (0.7 to greater than 1.0 days per 1° N), the spring phenology has been more dramatically affected by latitude gradient in this study [28,29,30]. In addition, the LCD in Beijing was significantly earlier than in the other two cities, but the disparity between the other two cities was not conspicuous. These inconsistent quantifications in geographic gradient might first be caused by the difference across observed species themselves [31]. Another reason is that spring phenology could be dramatically advanced and autumn phenology could be delayed by high urbanization intensity [19].

Consistent with previous studies, widespread advancements in the FLD and FFD and inapparent delay in LCD were found in this study [15,32,33,34]. However, the magnitude of the shifts in spring and autumn phenology varied among stations. For example, our results showed that the advancing trends of spring phenology in Beijing across the two periods (deviation in phenophases/29 years) were the most significant (0.28 days year⁻¹ for FFD and 0.38 days year⁻¹ for FLD), which were close to a previous study in China (0.22 to 0.57 days year⁻¹) but slightly higher than those observed in Europe (0.24 days year⁻¹) [15,34]. Nevertheless, the greatest changes in LCD were observed in Harbin (0.40 days year⁻¹ delay), which was much higher than both the previous study of China (0.19–0.25 delay per year) and Europe (0.036 delay per year) [15,34]. This is because the phenological changes in urban areas would be more dramatic than in natural environments, with 0.50 days year⁻¹ much earlier in spring and 0.78 days year⁻¹ much later in autumn [35]. It is worth noting that the phenological shifts in Mudanjiang were more slight than in other cities (0.08 days year⁻¹ for FLD, 0.04 days year⁻¹ for FFD, and 0.01 days year⁻¹ for LCD). This may be attributed to the weaker urban heat island effect in Mudanjiang compared to the other two cities.

Based on the linear regression between the preseason temperature and phenophases, we calculated the temperature sensitivity of spring and autumn phenology. Both temperature sensitivities of FLD and FFD in all three sites were negative overall, which is consistent with the results in Europe and North America [15,36]. In a previous study, spring temperature sensitivity in most of the regions of the mid- and high-latitude Northern Hemisphere showed positive clines across latitudes [37]. However, in this study, there is no evidence to suggest the temperature sensitivity of spring phenology changes along the latitudinal gradient. It indicates that spring temperature sensitivity for a single phenology site might be affected by multiple factors such as biodiversity, radiation, precipitation, soil nitrogen, etc. [38]. In our study area, the Beijing site is the Summer Palace, and the Harbin site is the Heilongjiang Forest Botanical Garden. Precipitation and drought events were not considered, as both sites are artificially irrigated. Furthermore, Figure 2 indicates that there are few significant differences in mean precipitation between 1963–1991 and 1992–2020. Consequently, it is not feasible to analyze the effect of precipitation on phenology in our study. Additionally, some studies suggest that the spring phenology of most vegetation is not sensitive to photoperiod [39,40], and the start of growing season for vegetation in the Northern Hemisphere is not significantly affected by radiation intensity [41]. Therefore, this study did not incorporate the effects of photoperiod or radiation into the analysis of phenological patterns.

In addition, different from the previous study, i.e., plants that have earlier phenophases are more sensitive to temperature than later plants, our results did not show any significant relationship between temperature sensitivities of spring phenology and the timing of phenophases [42]. Unlike spring phenology, the direction of LCD response to preseason temperature was not consistent across sites, where the sensitivity of LCD in Beijing and Harbin was positive, but in Mudanjiang, it was negative. Due to the different urbanization processes, autumn phenology in North China cities exhibited a greater delay (3.1 days to less than 1 day) than in Northeast China [18]. The other reason could be the inconsistent standard in LCD recording caused by subjective observations, which introduced potential errors in the linear regression model [43]. Further, autumn phenological changes are affected by complex factors, and the leaf coloring response may not be linear [44].

5. Conclusions

By utilizing FFD, FLD, and LCD data from the China Phenological Observation Network and daily temperature data, we analyzed the phenological changes in three cities (Beijing, Harbin, and Mudanjiang) between 1963–1991 and 1992–2020 to examine their response to urban warming. The results show the following: (1) FFD and FLD occurred earliest in Beijing, while LCD occurred earliest at Mudanjiang and latest in Beijing. (2) In Beijing, Harbin, and Mudanjiang, the median FFD (FLD) advanced by 8.10 (11.15), 2.99 (7.54), and 2.18 (1.15) days, respectively, during the period 1992–2020 compared to 1963–1991. Meanwhile, the median LCD was delayed by 5.10 days in Beijing and 7.05 days in Harbin but advanced by 0.17 days in Mudanjiang. (3) In Beijing, Harbin, and Mudanjiang, both FFD and FLD were significantly negatively correlated with preseason temperature. However, no consistent relationship was observed between LCD and preseason temperature. Our study suggests that vegetation phenology in cities will experience more pronounced changes than in natural environments due to the urban heat island effect.