Diverse Responses of Vegetation Dynamics to Snow Cover Phenology over the Boreal Region

Abstract

Snow cover phenology plays an important role in vegetation dynamics over the boreal region, but the observed evidence of this interaction is limited. A comprehensive understanding of the changes in vegetation dynamics and snow cover phenology as well as the interactions between them is urgently needed. To investigate this, we calculated two indicators, the start of the growing season (SOS) and the annual maximum enhanced vegetation index (EVI_max), as proxies of vegetation dynamics using the Moderate Resolution Imaging Spectroradiometer (MODIS) enhanced vegetation index (EVI). Snow cover duration (SCD) and snow cover end date (SCE) were also extracted from MODIS snow cover datasets. Then, we quantified the spatial-temporal changes in vegetation dynamics and snow cover phenology as well as the relationship between them over the boreal region. Our results showed that the EVI_max generally demonstrated an increasing trend, but SOS varied in different regions and vegetation types from 2001 to 2014. The earlier onset of SOS was mainly concentrated in the Siberian boreal region. In the Eurasian boreal region, we observed an advance in the SCE and decrease in the SCD, while in the North American boreal region, the spatial distribution of the trends exhibited substantial heterogeneity. Our results also indicated that the snow cover phenology had significant impacts on the SOS and the EVI_max, but the effects varied in different regions, vegetation types, and climate gradients. Our findings provide strong evidence of the interaction between vegetation dynamics and snow cover phenology, and snow cover should be considered when analyzing future vegetation dynamics in the boreal region.

1. Introduction

Both modelling and observational studies suggest that the mid to high latitudes of the Northern Hemisphere have undergone substantial warming in recent decades [1,2]. There is mounting scientific evidence that the raising temperature has a great impact on different aspects of vegetation, such as phenology, biomass, and productivity [3,4,5,6]. For changes in vegetation phenology, one of the most noteworthy findings was a wide range of warm-induced advances in the start of the growing season (SOS) that occurred in the Northern Hemisphere in the 1980s and 1990s, such as in North America [7], Eurasia [8], Europe [9], and the Arctic [10], although the magnitude of the trend varied among different research. Interestingly, multiple reports from both in situ and satellite observations indicated that since 2000, the advancing trends of the SOS have slowed down, with even a reverse trend occurring in some regions [3,11,12]. Meanwhile, for changes in vegetation greenness, a slowed and even reverse trend was also reported in different regions over the Northern Hemisphere after 2000 [13,14]. Considering that changes in vegetation may also lead to multiple feedbacks to the climate system [15], it is necessary to conduct further research on vegetation dynamics after 2000 to understand global climate change.

Temperature has been widely reported as a key factor that drives the changes in vegetation dynamics at global and regional scales [16,17,18]. However, for areas with seasonal snow cover, in addition to temperature, snow cover also affects vegetation dynamics, but the knowledge of the response of vegetation dynamics to snow cover phenology is still limited [19,20,21]. Over the past several decades, remarkable shifts in snow cover phenology have also been reported in many regions [22,23]. In situ observations indicated that between 1980 and 2006, the snow cover end date (SCE) showed a rapid advance in Eurasia at a rate of 6 ± 5.6 days decade⁻¹, and no trend was found in North America [24]. Satellite observations from 1982 to 2013 showed that the changes in the SCE in the Northern Hemisphere, Eurasia, and North America were −6.12, −7.17, and −3.78 days decade⁻¹, respectively; the snow cover duration (SCD) maintained a shortening trend of 1.04, 1.71, and 1.33 days decade⁻¹ in the Northern Hemisphere, Eurasia, and North America, respectively [22]. The changes in snow cover phenology may have significant impacts on vegetation dynamics for several reasons. First, although there is no direct relationship between snow and vegetation growth [21], changes in snow cover can have cascading effects on the known drivers of vegetation growth, such as the chilling requirement, forcing temperatures, and soil moisture status [13,15,16]. Second, changes in snow cover play a very important role in energy feedback through albedo, surface heating, and water budgeting [23,25]. In seasonal snow-covered regions, snow–vegetation feedback is considered as a key mechanism behind snow-driven growth changes [26]. Therefore, an understanding of the relationship between vegetation dynamics and snow cover phenology is vital for the projection of ecosystem dynamics in the future.

The sharp contrast in albedo between snow cover and vegetation, the potential release of large carbon stocks, and strong feedback on ocean-atmosphere circulation make the boreal region play a very important role in terrestrial ecosystems and regulation of the global climate system [27,28]. Therefore, studying the changes in vegetation dynamics in the boreal region and its relationship with snow phenology is very important to understanding the interaction between snow and vegetation under climate change. To understand the influence of snow in the boreal region on the spring phenology of vegetation, observations from passive microwave data and in situ flux data indicated that early snow melting led to the early arrival of spring recovery and enhanced boreal springtime carbon uptake [29]. Evidence from the community land model, which examines the biogeochemical and biogeophysical processes in terrestrial ecosystems, also showed that an increase in winter precipitation could alleviate the strong advancing trend of spring vegetation growth [27]. These conclusions are important, but additional driving mechanisms of different responses between snow cover phenology and vegetation dynamics need to be further explored.

In this study, based on the Moderate Resolution Imaging Spectroradiometer (MODIS) 0.05° data for the period 2001–2014, we first calculated the SOS and annual maximum enhanced vegetation index (EVI_max), which indicates the maximum level of seasonal photosynthetic capacity [30], as the main indicators of vegetation dynamics over the boreal region. SCD and SCE were also derived to represent the snow cover phenology variability. We then assessed the spatial-temporal trends of vegetation dynamics and snow cover phenology over the boreal region and identified the relationships between them. The objectives of this study are therefore to: (1) Quantify the spatial patterns of vegetation dynamics and snow cover phenology across the boreal region; (2) evaluate the temporal trends of vegetation dynamics and snow cover phenology over the whole boreal region and in different land cover types; (3) analyze the relationships between vegetation dynamics and snow cover phenology over the whole boreal region and in different land cover types; (4) explore the climate controls on the response of vegetation dynamics to snow cover phenology.

2. Materials and Methods

2.1. Study Area

The boreal region with an area of approximately 15.09 × 10⁸ ha covers one of the largest forest biomes on Earth and accounts for approximately 14.5% of the total land area (excluding glaciers) and 30% of the forest area [28,31]. This region forms a broad east-west belt stretching over northern regions of the Eurasian and North American continents, and 60% of this region occurs in Russia, 28% in Canada, and the remaining 12% is distributed in more than 10 countries, such as China, Finland, and Norway. In terms of the latitudinal division, the boreal region is located at latitudes of 45~75° (Figure A1). The mean annual temperature ranges from approximately −5 to +5 °C with annual precipitation ranging from 20 to almost 200 cm [31].

Most of the boreal region is occupied by needleleaf forest, which consists of evergreen needleleaf forests broadly distributed throughout the study area and deciduous needleleaf forests distributed in the East Siberian Taiga. The border to the north is in the transition zone between the boreal forest biome and the subarctic biome. To the south, the needleleaf forests are gradually replaced with various types of vegetation, such as deciduous forest, montane coniferous forest, and grassland, depending on the geographical location (Figure 1).

2.2. Data Source and Data Processing

2.2.1. MODIS EVI Product

To cover such a large area (approximately 15.09 × 10⁸ ha), we chose the MODIS enhanced vegetation index (EVI, MOD13C1) product (Collection 6) from 2001 to 2014, which has a spatial resolution of 0.05° (approximately 5 km) and a temporal resolution of 16 days (

Table A1. Detailed information of the data sources used here.

Dataset	Resolution	Cover Time	Download Website
MOD13C1	0.05°	2001–2014	https://search.earthdata.nasa.gov/search
MOD10C2	0.05°	2000–2014	https://nsidc.org
MCD12C1	0.05°	2001	https://search.earthdata.nasa.gov/search
CRU-TS 4.02	0.5°	2001–2014	https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_4.02/

All MODIS data belong to the Level 3 product, which is expressed in a uniform time-space grid and usually has certain integrity and consistency.

). This product can be acquired from NASA’s Earth Observing System data gateway (https://search.earthdata.nasa.gov/s-earch). The EVI was developed to be more sensitive in high biomass regions and minimize the disturbance of the soil background [32,33]. Furthermore, by utilizing a constrained-view angle-maximum value composite (CV-MVC) method, the EVI dataset was processed to reduce the disturbances from clouds, atmosphere, sensors, and surface bidirectional reflectance [32,34].

2.2.2. MODIS Snow Cover Product

The MODIS snow cover product used in this study was acquired from the National Snow and Ice Data Center (https://nsidc.org). Considering the persistent cloudiness, particularly at high latitudes in the Northern Hemisphere, we selected the 8-day composite data rather than the daily resolution data, which can eliminate cloud obscuration and provide more consistent and cloud-free data. Moreover, to remain consistent with the spatial resolution of the EVI data, we selected MOD10C2 snow cover fraction products to derive the snow cover phenology parameters across the boreal region (Table A1). The MOD10C2 snow cover fraction products are directly created from MOD10A2. In MOD10A2, each 500 m pixel is defined as snow, no snow, cloud, etc. When the data are transformed to 0.05° spatial resolution, the original MOD10A2 pixels and the total numbers of the corresponding pixels are considered to derive the fraction or percentage of snow cover [35]. We selected the MOD10C2 fractional snow cover product from January 2000 to December 2014 to derive the snow cover phenology of the boreal region.

2.2.3. Land Cover and Climate Data

To explore the relationship between vegetation dynamics and snow cover phenology in different vegetation types, we classified the land surface using MODIS land cover products. To remain consistent with the spatial resolution of the snow cover and vegetation data, we selected MCD12C1 products with the same resolution of 0.05° (Table A1). These products were derived using the same algorithm used to produce the Global 500 m Land Cover Type Product (MCD12Q1) and include three classification schemes that describe the land cover characteristics derived from the annual Terra and Aqua MODIS data observations. The product contains the following three classification schemes, which are derived from the supervised decision tree classification: 1. The International Geosphere-Biosphere Programme (IGBP) global vegetation classification scheme, 2. University of Maryland (UMD) Program, and 3. the LAI/FPAR scheme based on MODIS. As the IGBP classification scheme more prominently reflects the importance of natural vegetation in various land cover types [36], we chose the IGBP classification scheme. There are 16 land cover types in the boreal region according to the IGBP classification, excluding evergreen broadleaf forest. After eliminating the nonvegetated area and croplands, we reclassified the land cover types into 7 major land cover types (

Table A2. Harmonized legend used and correspondence with the original product classes.

Classes	Value	Number of Pixels
Deciduous needleleaf forest	3	25,171
Evergreen needleleaf forest	1	37,649
Grasslands	10	9538
Mixed forest	5	63,022
Open shrublands	6, 7	117,121
Savannas	8, 9	80,000
Permanent wetlands	11	8475

The non-vegetated parts and crops were masked out.

) in this study: ENF (evergreen needleleaf forest), DNF (deciduous needleleaf forest), MF (mixed forest), OSH (open shrublands), SA (savannas), GL (grasslands), and WET (permanent wetlands). The distribution of the major vegetation types in the boreal region is depicted in Figure 1.

The monthly mean temperature and precipitation data with a spatial resolution of 0.5° were extracted from the CRU-TS 4.02 climate dataset [37], covering the period from 2001 to 2014 (Table A1). CRU-TS databases have been widely used in recent climate change and vegetation growth dynamic studies [21,38]. In order to match the MODIS products, we used the nearest-neighbor interpolation method to interpolate climate data into the resolution of 0.05°.

2.3. Vegetation Dynamics Metrics

Considering the influence of noise in the original synthesized EVI time series data, many scholars have proposed several effective means to gap-fill and smooth the raw datasets, such as the Savitzky-Golay filtering method [39], the asymmetric Gaussian fitting method [40], the double logistic function [41], and the harmonic analysis of time series method (Hants) [42]. The double logistic method has been proven to be more suitable than other algorithms at high latitudes, such as Savitzky-Golay and asymmetric Gaussian [16,43]. Here, we used the double logistic method to smooth and reconstruct the EVI time series dataset over the boreal region. The basic function of the double logistic algorithm is expressed as:

$$\mathrm{g}\left(t;x_1,\dots ,x_4\right)=\frac{1}{1+e^{\left(\frac{x_1-t}{x_2}\right)}}+\frac{1}{1+e^{\left(\frac{x_3-t}{x_4}\right)}},$$

(1)

where $t$ is the Julian day, $x_1$ and $x_3$ determine the left and right inflection points, respectively, while $x_2$ and $x_4$ indicate the corresponding rates of change.

Declining snow cover may be responsible for the change of EVI associated with the greening trend. To avoid the effects of snow melt, we assigned the negative EVI value to 0 [43,44]. Then, we employed TIMESAT software, which includes a double logistic algorithm and has been widely used to investigate vegetation dynamics [17,44,45]. It allows various vegetation parameters to be derived, including the SOS and EVI_max (the maximum value of the annual EVI time series fitted by the double logistic curve). We used the following parameters: Adaptation strength of 2, number of envelope iterations of 2, seasonal parameter of 1, and amplitude season start and end of 20%, following Karkauskaite et al. [46].

2.4. Snow Cover Phenology Metrics

Based on the seasonal cycle of snow cover over the boreal region, we used the hydrological year instead of the calendar year to derive the snow cover phenology parameters [22,47]. Here, we defined the hydrological year of the boreal region as the period from September 1st of a given year (t − 1) to August 31st of the next year (t). Then, we defined the onset season of snow cover from September 1st of a given year (t − 1) to the last day in February of the following year (t), while the melting season spanned from March to August of the same year (t). Referring to the definition of snow cover parameters in Chen et al. [22], during the onset season, we searched for the interval (i to i + 7) when the snow cover fraction was equal to 0 and subsequently superior to 0. Then, we defined the date (i + 3.5) as the snow cover onset date (SCO). Likewise, during the melting season, we searched for the interval (i to i + 7) when the snow cover fraction was above 0 and subsequently below 0. Then, we defined the date (i + 3.5) as the SCE. The SCD was defined as the interval between SCO and SCE. This method is one of the easiest ways to calculate snow cover phenology, and can avoid the impact of ephemeral snow [22,24,48].

2.5. Mann–Kendall Trend Test

The nonparametric Mann-Kendall (MK) trend test [49,50] was used to analyze the annual trends in both vegetation dynamics and snow cover phenology. Compared with parametric tests, the nonparametric MK trend test is applicable to non-normally distributed data and is not impacted by a few abnormal values, and this method has been widely applied for trend detection [34,51,52]. The Theil-Sen median slope estimator was used to analyze the slope of the annual variations in vegetation dynamics and snow cover phenology. The equation is:

$$\mathsf{\beta }=\mathrm{median}\left(\frac{x_{i-}{x}_j}{i-j}\right),i>j,$$

(2)

where $\mathsf{\beta }>0$ indicates an increasing trend in SCD and EVI_max and an advancing trend in SOS and SCE, and $\mathsf{\beta }<0$ indicates a reverse trend.

2.6. Partial Correlation Analysis

The relationship between multivariate related variables is complex, and there are often different degrees of simple correlation between any two variables [53]. Partial correlation analysis is a statistical method that studies the correlation between two related variables by fixing other related variables, which can eliminate the influence of other variables [38]. This method has been successfully applied to previous studies involving climate change and vegetation phenology [24,38,53]. The partial correlation coefficient was calculated as follows:

$$r_{xy.z}=\frac{r_{xy}-r_{yz}{r}_{xz}}{\sqrt{\left(1-r_{yz}{}^2\right)\left(1-r_{xz}{}^2\right)}},$$

(3)

where $r_{xy.z}$ is the partial correlation coefficient between variable x and y after fixing variable z, and $r_{xy}$, $r_{yz}$, and $r_{xz}$ are the correlation coefficients between two variables. In addition, a t-test was performed to reflect the significance level (p-values) of the relationship.

3. Results

3.1. Spatial Patterns of Vegetation dynamics and Snow Cover Phenology

To analyze the spatial distribution of the vegetation dynamics and snow cover phenology in the boreal region, we calculated the average values of the SOS, EVI_max, SCD, and SCE from 2001 to 2014 (Figure 2). The average SOS ranged from approximately 100 to 160 days (Figure 2a), accounting for 94.6% of the area of the boreal region. The SOS is generally delayed with increasing latitude. In addition, due to the influence of the altitude, climate, and vegetation types, the distribution of the SOS at the same latitude has great spatial variability. Furthermore, the average spring phenology based on vegetation types indicates that the evergreen needle forest has the earliest SOS, while the grasslands and open shrublands have the latest SOS. Compared with the SOS, the EVI_max showed an almost opposite trend in spatial distribution. The average EVI_max mainly ranged from 0.1 to 0.7 (Figure 2b), accounting for 98.6% of the area of the boreal region. The maximum mean EVI_max occurred in the mixed forest and deciduous needle forest, and the minimum mean EVI_max occurred in grasslands (

Table A3. Average and standard deviation (in brackets) of the start of the growing season (SOS, in Julian days), annual maximum EVI (EVI_max), snow cover duration (SCD, days), and snow cover end date (SCE, in Julian days) for different land cover types.

Land Cover Types	SOS	EVImax	SCD	SCE
Deciduous needleleaf forest	130 (8)	0.43 (0.05)	193 (13)	134 (7)
Evergreen needleleaf forest	124 (19)	0.35 (0.08)	174 (19)	126 (14)
Grasslands	139 (21)	0.33 (0.14)	190 (37)	134 (28)
Mixed forest	132 (16)	0.48 (0.08)	165 (21)	120 (14)
Open shrublands	145 (12)	0.36 (0.08)	213 (22)	144 (14)
Savannas	136 (15)	0.38 (0.08)	198 (23)	138 (17)
Permanent wetlands	127 (13)	0.34 (0.10)	190 (18)	134 (14)

).

The average spatial patterns of annual SCD mainly ranged from 130 to 250 days (Figure 2c), accounting for 97.4% of the area of the boreal region. Similar to SOS, SCD increased with latitude, but the trend of the change in SCD with latitude was more uniform. The shortest SCD occurred in mixed forests and evergreen needle forests; the longest SCD mainly occurred in open shrublands. Compared with the SCD, the SCE presented nearly the same trend in spatial distribution, except for some differences in North America. The average SCE mainly ranged from 100 to 160 days (Figure 2d), accounting for 94.4% of the boreal region. The latest SCE was mainly distributed in open shrublands; the earliest SCE was mainly distributed in mixed forests and evergreen needle forests (Table A3).

3.2. Temporal Trends in Vegetation Dynamics and Snow Cover Phenology

Figure 3 shows the spatial patterns of temporal trends in the SOS, EVI_max, SCD, and SCE across the boreal region. Significant trends were determined by the Mann-Kendall test. For the vegetation dynamics, the temporal trends of the SOS and EVI_max showed very different spatial distribution patterns. In the North American and European boreal regions, the SOS in most parts exhibited a delayed trend, while most of Siberia showed an advancing trend, and the significant changes were concentrated in specific areas. However, for the EVI_max, most of the boreal region showed a significant increasing trend. Similarly, the SOS trend varied in different vegetation types, but the EVI_max showed an increasing trend in every vegetation type. The SOS showed a nonsignificant advancing trend in deciduous needle forests and open shrublands and a nonsignificant delayed trend in evergreen needle forests and grasslands. However, the EVI_max showed a significant increasing trend in most vegetation types except deciduous needle forest (Figure 4a,b).

Unlike the different spatial patterns observed in SOS and EVI_max, the trends of SCD and SCE shared almost the same spatial pattern, but the extents of the trends varied. The SCD and SCE demonstrated clear spatial variability throughout the boreal region. A visible increasing SCD and advancing SCE was found in most of the Eurasia boreal region except for some regions near the Sea of Okhotsk and some regions in the Scandinavian and Russian taiga. However, the SCD and SCE trends were non-uniform in the North American boreal region. A visible increasing SCD and advancing SCE was found in Northern Cordillera forests and mid-continental Canadian forests, however, a visible decreasing SCD and delayed SCE was found in most of the northeast and northwest North American boreal region. Interestingly, for different vegetation types, the temporal trends in SCD and SCE were not as consistent as the spatial patterns. The SCD showed a nonsignificant decreasing trend in all vegetation types, however, a nonsignificant advancing SCE trend was found in evergreen needle forests and grasslands (Figure 4c,d).

3.3. Response of Vegetation dynamics to Snow Cover Phenology

Figure 5 shows the spatial distribution of the partial correlation coefficients between vegetation dynamics and snow cover phenology. Interestingly, the partial correlations between SOS and snow cover phenology showed two very different patterns, but the EVI_max and snow cover phenology exhibited almost the opposite spatial distributions (

Table A4. Percentages of the partial correlation coefficient between the start of growing season and snow cover duration (SCD_SOS), start of growing season and snow cover end date (SCE_SOS), annual maximum EVI and snow cover duration (SCD_EVImax), and annual maximum EVI and snow cover end date (SCE_EVImax) in the entire boreal region (All), North American boreal region (NA), and Eurasian boreal region (EA).

	Partial Correlation	All	NA	EA
SCD_EVImax	Positive	47.05%	42.11%	49.28%
	Negative	52.95%	57.89%	50.72%
SCD_SOS	Positive	56.32%	45.26%	61.32%
	Negative	43.68%	54.74%	38.68%
SCE_EVImax	Positive	57.71%	61.34%	56.06%
	Negative	42.29%	38.66%	43.94%
SCE_SOS	Positive	62.88%	68.41%	60.34%
	Negative	37.12%	31.59%	39.66%

). SOS and SCD were positively correlated in 56.32% of the boreal region, but more pronounced in Eurasia (61.32%) than in NA (45.26%). Additionally, the most significant positive correlation was found in the northeast and the northern region of East Siberia, while significant negative correlations were scattered throughout North America (Figure 5a). The positive correlation between SOS and SCE accounted for 62.88% of the boreal region. However, unlike the spatial distribution of the correlation between SOS and SCD, the correlation between SOS and SCE is higher in North America (68.41%) than in Eurasia (60.34%). Significant positive correlations were mainly concentrated in the northern part of east and west Siberia and the northwest part of North America and regions near James Bay. However, negative correlations were rare, mainly in northeastern Siberia (Figure 5b). For all vegetation types, the partial correlation between SOS and snow cover phenology showed no significant partial correlation (Figure 6a,b).

The EVI_max was negatively correlated with SCD throughout the boreal region (52.95%), Eurasian boreal region (50.72%), and North American boreal region (57.89%) (Figure 5c). However, the relationship between the EVI_max and SCE was almost the opposite, with a positive correlation demonstrated throughout the boreal region (57.71%), Eurasian boreal region (56.06%), and North American boreal region (61.34%) (Figure 5d). The spatial distribution of the significant partial correlation between the EVI_max and snow cover phenology was similar and was discretely distributed throughout the boreal region. Similarly, the partial correlation between the EVI_max and snow cover phenology of different vegetation types also showed the opposite pattern (Figure 6c,d). In addition, among the partial correlations between the EVI_max and SCD of different vegetation types, only mixed forests showed significant negative correlations, and among the partial correlations between the EVI_max and SCE, grasslands and mixed forests showed significant positive correlations.

3.4. Climate Controls on the Response of Vegetation Dynamics to Snow Cover Phenology

To study the control of the climate on the response of vegetation dynamics to snow cover phenology, we first used the annual temperature and precipitation data from 2000 to 2014 to obtain the annual average temperature and annual cumulative precipitation of these 14 years (Figure 7). As shown in Figure 7, the temperature gradient in the northern region clearly varied with latitude, and the coldest regions were located in eastern and northeastern Siberia. The trend in the precipitation gradient with latitude was not obvious. Eurasia showed high precipitation in the west and low precipitation in the east, while North America showed high precipitation in the east and low in the west.

Using the obtained temperature and precipitation gradients, we calculated the partial correlation coefficients between the vegetation dynamics and snow cover phenology under each hydrothermal gradient (Figure 8). As shown in Figure 8, under different hydrothermal gradients, the effects of snow cover phenology on vegetation dynamics showed clear aggregation and heterogeneity. SOS and SCD mainly showed positive partial correlations in different precipitation gradients below −10 °C. Above −10 °C, with the increase in the precipitation gradient, the partial correlation between SOS and SCD appeared to change from negative to positive and then negative. The significant partial correlations were mainly concentrated in areas with temperatures below −10 °C and where the temperature ranged from −5 to 0 °C and the precipitation was approximately 1000 mm (Figure 8a). However, SOS and SCE showed positive correlations in most areas; the two mainly showed significant positive partial correlations below 0 °C, negative partial correlations above 0 °C and below 1000 mm precipitation, but strong heterogeneity above 1000 mm precipitation (Figure 8b). Interestingly, the partial correlation between the EVI_max and snow cover phenology also showed opposite distribution patterns in different hydrothermal gradients. Below −10 °C, the EVI_max and SCD mainly showed positive partial correlations and significant positive correlations were mainly concentrated in the −10 °C region. The EVI_max and SCE mainly showed nonsignificant negative partial correlations below −10 °C. Above −10 °C, the EVI_max and SCD mainly exhibited negative partial correlations, while the EVI_max and SCE mainly exhibited positive partial correlations, and significant partial correlations were mainly concentrated in the region of −5 to 0 °C.

4. Discussion

4.1. Temporal Trends in Vegetation Dynamics and Snow Cover Phenology

For the temporal trend in the SOS, the major advancing SOS trends occurred in the Eurasian boreal region, while a reversed delayed SOS trend occurred in the Scandinavian and Russian taiga and some parts of North America. These results were consistent with the conclusions of Barichivich [54] and Wang [55]. Moreover, SOS also showed different trends in different vegetation types. The cause of the two different SOS trends may be the combined effects of a warmer November and December (late winter) and colder January and March (early spring), which delayed the spring phenology of earlier-unfolding plants and accelerated the spring phenology of later-unfolding plants, respectively [11,56]. However, there were significant increases in the EVI_max in different regions as well as in different vegetation types, which indicated that the increasing temperature enhanced the growth of vegetation in summer over the boreal region. For snow cover phenology, most advancing SCE trends were distributed in Eurasia, with a few distributed in the high latitudes of boreal NA (greater than 60° N); simultaneously, the major delayed SCE trends were distributed in the middle latitudes of boreal NA (less than 60° N). This result is consistent with the conclusion of Chen et al. [22,35]. The spatial distribution pattern of SCD was very similar to that of SCE, but the degrees of the changes in SCD and SCE in different vegetation types were very different. This difference was because the change in SCD was more complicated than the change in SCE, while SCD was affected by SCE, the SCO, and the snow depth [24]. Throughout the boreal region, SCD mainly showed a decreasing trend, which was consistent with the previous findings of Choi et al. [57], Chen et al. [22,35], and Peng et al. [24].

4.2. Influence of Snow Cover Phenology on Vegetation Dynamics

Previous studies suggested that spring vegetation phenology was principally affected by temperature [4,45,53,58] and many other factors, such as permafrost degradation [59], photoperiod [60], and winter chilling [53,61]. Our results also showed that the spring vegetation phenology was associated with snow cover phenology across the boreal region. This result was supported by Wang et al. [62] and Wang et al. [20], who showed that snow cover phenology had an important impact on SOS for the alpine plants over the Tibetan Plateau. In many parts of the Eurasian boreal region, SCD had a significant positive partial correlation with SOS. These results were consistent with previous conclusions that longer snow cover in winter delayed the arrival of spring phenology [53]. Moreover, in some parts of the North American boreal region, significant negative correlations between SCD and SOS could also be observed, and grasslands and mixed forests also exhibited nonsignificant negative partial correlations among different vegetation types. On the one hand, this result may be due to the strong correlation between SCD and winter chilling and the inconsistent demands of different vegetation types for winter chilling [54]. On the other hand, snow cover is also an effective insulator, which can protect the soil from exposure to harsh climatic conditions, such as low temperature and strong winds [43]. Therefore, the soil temperature in snow-covered areas was usually higher than that in snow-free areas [47], which was conducive to the SOS. Similar to SCD, SCE showed an inconsistent correlation with SOS over the boreal region, but the regions that showed significant positive correlations with SOS were larger than the areas that showed significant correlations between SCD and SOS. Compared to SCD, SCE had a more complex influence mechanism on the spring phenology of vegetation. In addition to affecting winter chilling, the change in SCE had a certain impact on the early spring temperature and soil moisture. On the one hand, the advancing trend of SCE increased the temperature in early spring and provided moisture for the growth of vegetation, which further promoted the arrival of vegetation in spring [21,63]. On the other hand, early snow cover melting may lead to more frost, and vegetation may choose to delay greening to protect against the damage caused by low temperatures [64].

Although snow cover and plant productivity were not directly related, snow cover changes can have a cascade effect on the drivers of known growth season productivity, such as soil moisture conditions [21]. Fortunately, snow manipulation experiments in the boreal region can help us understand how snow cover affects vegetation productivity. On the one hand, the melting time of the snow cover was closely related to the length of the vegetation growth period; on the other hand, the amount of snow was closely related to the soil moisture in the vegetation growing season [65]. Both processes had an important impact on vegetation productivity. Therefore, due to different vegetation types and climatic conditions, the relationship between the EVI_max and snow cover phenology showed great heterogeneity.

4.3. Impact of Climate in Vegetation Dynamics Response to Snow Cover Phenology

As indicated above, the response of vegetation dynamics to snow cover phenology showed strong heterogeneity in both spatial distribution and different vegetation types. By calculating the partial correlation between vegetation dynamics and snow cover phenology under different climatic gradients, we found that the partial correlation between them showed clear aggregation in different gradients. Interestingly, when the temperature was below −10 °C, there were positive partial correlations between SCD and SOS and SCD and EVI_max, which may be due to the shortening of SCD and the consequent increases in soil freezing and root mortality, resulting in delayed SOS and reduced vegetation growth [43]. When the temperature was greater than −10 °C, negative partial correlations began to appear under different precipitation gradients. Because SCD was closely related to winter chilling and less chilling accumulation delayed the SOS and further affected EVI_max [66], significant correlations between SCE and SOS as well as between SCE and EVI_max were mainly concentrated at −5 °C to 0 °C, with a significant positive partial correlation. This relationship may be due to the earlier SCE, which caused the spring phenology to advance and led to drier soil moisture, thereby reducing the EVI_max [47]. In other climate gradients, both positive and negative correlations were found, which were not significant. This result may be because the lag effects of snow cover phenology on the EVI_max in these regions were limited, and the EVI_max was more likely to be affected by the climate of the growing season.

5. Conclusions

Over the past few decades, the climate has undergone significant changes. Understanding the spatial-temporal changes in vegetation dynamics and snow cover phenology in the boreal region and their links is of great significance for global change research. Using remote sensing data from 2000 to 2014, we found that the annual maximum EVI (EVI_max) generally showed an increasing trend, but the start of growing season (SOS) varied in different regions and vegetation types, and the advancing trends were mainly concentrated in the Siberian boreal region. In the Eurasian boreal region, snow cover duration (SCD) mainly showed decreasing trends, snow cover end date (SCE) mainly showed advancing trends, while the spatial distribution of the two trends exhibited great heterogeneity in the North American boreal region. Our results also indicated that the snow cover phenology had a significant impact on the SOS and EVI_max, but the effect varied in different regions, vegetation types, and climate gradients. The uneven effect of snow cover phenology on vegetation dynamics emphasized the important role of snow cover in spring and summer vegetation physiology and the challenge of accurately characterizing the growth of vegetation with only temperature and precipitation in most current ecosystem models. Therefore, snow must be regarded as a non-negligible variable. The study of snow cover should not be limited to its feedback on air temperature, and more research should be conducted on other effects on soil temperature in winter and soil moisture in the growing season as well as further impacts on the structure and function of ecosystems. In the northern high latitude, a widespread decline in winter snow accumulation has been observed in recent decades and climate models unanimously project a decrease in the snow-to-rain ratio during the winter season in a warmer future [1,67]. For future research, it is very important that the dynamics of vegetation growth, changes in both snow cover phenology and snow accumulation, and their interactions under different vegetation types, regions, and climatic gradients are evaluated through a variety of models and different kinds of data sources, such as ground data, remote sensing data, and flux towers data.