Spatiotemporal patterns of remotely sensed phenology and their response to climate change and topography in subtropical bamboo forests during 2001-2017: a case study in Zhejiang Province, China

ABSTRACT

Vegetation phenology has long been adapted to environmental change and is highly sensitive to climate change. Shifts in phenology also affect feedbacks of vegetation to environmental factors such as topography and climate by influencing spatiotemporal fluctuations in productivity, carbon fixation, and the carbon water cycle. However, there are limited studies which explores the combined effects of the climate and terrain on phenology. Bamboo forests exhibit the outstanding phenological phenomena and play an important role in maintaining global carbon balance in climate change. Therefore, the interaction mechanisms of climate and topography on bamboo forest phenology were analyzed in Zhejiang Province, China during 2001–2017. The partial least squares path model was applied to clarify the interplay between the climate and terrain impacts on phenology under land cover/use change. The results revealed that the average start date of the growing season (SOS) significantly advanced by 0.81 days annually, the end date of the growing season (EOS) was delayed by 0.27 days annually, and the length of the growing season (LOS) increased by 1.08 days annually. There were obvious spatial differences in the partial correlation coefficients between the climate factors and phenological metrics. Although the SOS, EOS and LOS were affected by different climatic factors, precipitation was the dominant factor. Due to the sensitivity of the SOS and EOS to precipitation, a 100 mm increase in regional annual precipitation would cause the average SOS to advance by 0.18 days and the EOS to be delayed by 0.12 days. Regarding the terrain factors affecting climate conditions, there were clear differences in the influences of different altitudes, slopes and aspect gradients on bamboo forest phenology. This study further showed that topographic factors mainly affected the interannual variations in phenological metrics under land cover/use change by affecting precipitation. This study clarified the spatial pattern of bamboo forest phenology and the interactive mechanisms between vegetative phenology and environmental conditions, as this information is crucial in assessing the impact of phenological change on the carbon sequestration potential of bamboo forests.

KEYWORDS:

• Bamboo forest
• phenology
• climate change
• topographic effects

1. Introduction

Vegetation phenology refers to regular changes in vegetation (leaf spreading, flowering, deciduous leaves etc.) with cyclical changes in seasonal climate (Zhang et al. 2003; Broich et al. 2014; Wang, Luo, and Shafeeque 2019; Li et al. 2021c). Vegetation phenology is particularly susceptible to climate change and can serve as a clear signal of global warming over extended time scales (Sun et al. 2004; Liu et al. 2015; Zhao et al. 2022). Climate change affects phenology and controls seasonal changes in vegetation growth through phenology. Phenology regulates vegetation photosynthesis and other physiological ecological processes and has a strong impact on spatiotemporal changes in ecosystem productivity and the carbon cycle (Sacks, Schimel, and Monson 2007; Han et al. 2018; Richardson et al. 2010; Ahlström et al. 2015; Wu et al. 2021; Li et al. 2021a). Therefore, studying the response mechanisms of vegetation phenology and climate change is important for understanding material and energy exchanges between climate and vegetation and for accurately evaluating vegetation productivity and the global carbon budget.

Significant changes in plant phenology (the start and end dates and length of the growing season, SOS/EOS/LOS) have resulted from climatic warming, and the most notable changes are an earlier SOS in spring, a delayed EOS in the fall, and an extended LOS at the regional and global level (Yuan et al. 2020; Shen et al. 2019; Yang et al. 2022; Zu et al. 2018; Julien and Sobrino 2009; Zhao et al. 2015). Previous studies have shown that vegetation phenology is greatly influenced by climatic factors (temperature, precipitation, and radiation) (Xia et al. 2019; Suepa et al. 2016; Ge et al. 2021b; Gao and Zhao 2022; Li et al. 2021b). Phenology has a different sensitivity to climate in different regions. Temperature is a major climatic variable that controls the SOS and EOS in temperate forests (Xia et al. 2019; Gao and Zhao 2022; Jeong et al. 2011; Zhao et al. 2015; Ding et al. 2016). Gao and Zhao (2022) indicated that surface temperature has a substantial impact on advanced spring phenology in the north temperate zone. Additionally, precipitation has a significant role in influencing phenology in a complex manner in northern temperate forests (Piao et al. 2006; Fu et al. 2018).

Recent studies have also found that phenology is affected not only by climate change but also by topographic conditions (Yuan et al. 2020; Vitasse, Signarbieux, and Fu 2018; Xia et al. 2019). Topography influences the growing environment of vegetation and the climate change within forests. Vitasse, Signarbieux, and Fu (2018) indicated that high temperatures might encourage phenological progress at lower elevations. Yuan et al. (2020) found that the SOS had a nonlinear delay with increasing elevation. Therefore, the combination of climatic and geographic variables affects vegetation phenology, and the responses and feedbacks to climate change are inconsistent across different regions.

The most efficient technique used to assess how large-scale vegetation phenology responds to climate change and topography through the application of remote sensing data (Wang et al. 2019a; Xu et al. 2020; Silveira et al. 2019; Zhao et al. 2022; Piao et al. 2019). Remote sensing-based research on vegetation phenology has frequently ignored the effects of changing land cover/use, which weakens and obscures the understanding of how climate change and topography affect plant phenology (Gao et al. 2019; Wang et al. 2017; Wang, Luo, and Shafeeque 2019; Wang et al. 2022). Therefore, it is unclear whether the phenological tendency is consistent or varies with climate change and different topographic conditions under a background of land cover/use change.

Response of vegetation phenology to climate change and topographic conditions can exhibit complex (Wang et al. 2022; Zani et al. 2020; Zu et al. 2018). The partial correlation coefficient can better reflect the relationship between the two variables by eliminating the influence of other variables. The structural equation model can clearly analyze the role of individual variables on the whole and the relationship between individual variables, and can effectively solve the problem of multivariable complex collinearity (Jakobowicz and Derquenne 2007; Chin and Newsted 1999; Sarstedt, Ringle, and Hair 2022). Recently the partial correlation coefficient and structural equation model have been widely applied on phenology researches (Li et al. 2020; Mao et al. 2022; Jiang et al. 2020; Wang et al. 2022).

Although an increasing number of studies have reported how climate change and topographic conditions affect the phenology of plants in temperate forests (Zhao et al. 2022; Iler et al. 2017; Li et al. 2022; Xia et al. 2019; Wu et al. 2021), there have been few reports on tropical and subtropical forests. Tropical and subtropical forests have diverse species compositions and frequent human interference, resulting in the effects of their phenology on environmental conditions being more complex than those in temperate forests. Bamboo forests form a significant component of subtropical forest ecosystems in China. Bamboo forests are distinguished by on- and off-years in contrast to the phenological traits of other forest types (Mao et al. 2016; Mei et al. 2020). The “on-year” refers to bamboo shoots being exposed at the surface in mid-March and early April and then growing rapidly above the ground until reaching their full height in April and May. The “off-year” refers to the leaves of old bamboos turn yellow gradually in mid-March and early April, and then the leaves fall off, after which new leaves grow in April and May of the next year (Chen 2010; Mei et al. 2020; Li et al. 2021c). Therefore, it is crucial to examine the processes by which the phenological features of bamboo forests respond to topographic and climatic changes.

Bamboo forests, with an area of 6.41 million ha, account for 2.91% of China’s forest area, which is mostly dispersed among 15 provinces in Southeast China (Cui et al. 2019; Kang et al. 2022). Bamboo forests are highly effective at sequestering carbon, maintaining the global carbon balance, and adapting to climate change according to previous research (Li et al. 2019a; Mao et al. 2017; Kang et al. 2022; Mao et al. 2020, 2022). Therefore, to better understand how bamboo forest ecosystems respond and adapt to climate change, this study clarified the interactive response mechanism of bamboo forest phenology to climate change and geographical conditions under the condition of land cover/use change, which is of great scientific significance. The subject of the present study was the bamboo forest resource of Zhejiang Province. The phenological metrics (SOS, EOS and LOS) and distribution information for these bamboo forests were extracted using MODIS products from 2001 to 2017, and features of the geographical evolution of the bamboo forest phenology as well as the connection between the phenological metrics and environmental variables were examined. With reference to land cover/use change, the following were the primary goals of this study: (1) to analyze the spatiotemporal variation in the relationships among phenology and climatic and topographical factors and (2) to investigate the impacts of geographical and climatic conditions on the phenological metrics of bamboo forests under land cover/use changes.

2. Study area and data

2.1. Study area

Zhejiang Province, which is situated on China’s southeast coast, has a typical subtropical monsoon climate (Figure 1a). The climate is controlled by the subtropical anticyclone and has high temperature, rain and air humidity in summer. The area is controlled by the northern Mongolian high pressure, with sunny, cold, dry, and low-temperature weather and little rain in winter. The terrain inclines from the southwest to the northeast, and the average elevation is between 300 and 500 m. The central region has hills and basins, while the northeast region is mostly plains. The bamboo resources are very rich, with 0.83 million ha, and most are distributed in the mountainous and hilly areas of Zhejiang Province (Figures 1(b,c)).

Figure 1. (a) Meteorological stations, (b) bamboo abundance distribution and (c) elevation data of Zhejiang Province.

2.2. Data

2.2.1. Climate data

China’s Meteorological Data Center (http://data.cma.cn/), which has 824 observation stations (Figure 1), provided the daily climatic data, including daytime temperature (Tmax, unit: 0.1°C), nighttime temperature (Tmin, unit: 0.1°C), precipitation (Pre, unit: 0.1 mm), and solar radiation (Rad, unit: 0.01 W m⁻²). Daily climate data with a 500 m spatial resolution were used to match the MODIS LAI data with the inverse distance weighted method based on these 824 meteorological observation stations. The spatial distribution of climate data obtained by inverse distance interpolation method are basically consistent with the actual situation of China in 2001–2017. The climate data for Zhejiang Province were obtained from the national climate data. The annual mean climate data values were calculated from the averages or sums of the daily data and used to analyze the spatiotemporal impact mechanism of climate on phenology. The interannual change values of climate data were calculated from the averages of the regional annual mean climate data values and used to analyze the interaction mechanism of climate and topography on phenology under land cover/use change.

2.2.2. Topographic data

The DEM data were provided by the Geospatial Data Cloud site, Computer Network Information Center, Chinese Academy of Sciences (http://www.gscloud.cn). The local average approach described by Shang et al. (2013) was used to resample the DEM with the original spatial resolution of 30 m to 500 m. The altitude values were obtained from the DEM. The “slope” and “aspect” functions of the “Spatial Analyst” module in ArcGIS v. 10.6 were used to derive the slope and aspect from the DEM.

To explore the influence of topographical factors (e.g. altitude, slope and aspect) on bamboo forest phenological metrics, we calculated the mean value of the phenological metrics as the altitude changed from 50 to 1500 m at 50-m intervals and as the slope changed from 0° to 50° at 1°-intervals. The aspect was divided into 8 cardinal directions by rotating clockwise from the north at an intersection angle of 45° according to Bindajam et al. (2020), and we calculated the mean values of the bamboo forest phenological metrics in these eight directions.

3. Methods

Figure 2 shows the detailed research process and includes the following 5 steps: (1) environmental data collection and preprocessing, including collection of meteorological station data and interpolated data with 500 m spatial resolution, and collection of DEM data images and calculated altitude, slope and aspect dataset; (2) bamboo forest abundance information was extracted by coupling the maximum likelihood method and least square mixed pixel decomposition method using MODIS NDVI and reflectance data; (3) MODIS LAI time series data were assimilated using the particle filter algorithm, and then spatiotemporal distribution of bamboo forest phenology was extracted using a logistic regression equation based on assimilated MODIS LAI and bamboo forest abundance information; (4) The spatiotemporal patterns of bamboo forest phenology were analyzed by using Mann-Kendall test and Sen’s slope estimator; (5) the impact mechanism of climate and topographical factors on phenology was analysis based on linear trend analysis, partial correlation coefficient and structural equation model.

Figure 2. The flowchart of spatiotemporal patterns of bamboo forest phenology and its response to climate change and topography.

3.1. Extracting bamboo forest information from remote sensing data

In this study, land surface reflectance (MOD09A1) and NDVI time series data (MOD13Q1) were employed to extract the bamboo forest abundance distribution. Previous studies have reported the extraction of bamboo forest abundance information by coupling the maximum likelihood method and the fully constrained least square mixed pixel decomposition method using MODIS NDVI and reflectance data (Li et al. 2019a, 2018); the kappa coefficients were >80%, and the estimated bamboo forest area was highly consistent with the continuous forest resource inventory data of Zhejiang Province (R² >0.80). Because the annual change in the bamboo forest distribution was not too large and the workload was relatively low, we mapped the bamboo forest distribution every two years based on this method during 2001–2018.

3.2. Extraction of phenological metrics from assimilated MODIS LAI

The MODIS leaf area index (LAI) products (MOD15A2) were used to extract the bamboo forest phenological metrics from 2001 to 2017. However, the MODIS LAI time series data cannot continuously monitor vegetation growth processes due to atmospheric conditions, including clouds, aerosols, water vapor, and ozone, which affect the estimation accuracy of the phenological metrics (Vermote et al. 1997; Chen, Deng, and Chen 2006; Fang et al. 2008; Xiao, Wang, and Wang 2008). We used a particle filter algorithm to assimilate the MODIS LAI time series, and the R² and RMSE values between the observed and assimilated LAI were 4.49% higher and 58.54% lower, respectively, than those between the observed and smoothed MODIS LAI (Appendix S1). Therefore, the bamboo forest phenological metrics were extracted by combining the assimilated MODIS LAI time series data based on the particle filter algorithm and a logistic regression equation (Li et al. 2021c).

We used the maximum assimilated LAI of a year to divide the growth process into ascending and descending sections (Figure 3a). The rate of change (ROC) (see Equation (2))(2) $ROC=b^{3}cz\{\frac{3z\left(1-z\right){\left(1+z\right)}^3\left[2{\left(1+z\right)}^3+b^2{c}^{2}z\right]}{{\left[{\left(1+z\right)}^4+{\left(bcz\right)}^2\right]}^{2.5}}-\frac{{\left(1+z\right)}^2\left(1+2z-5z^2\right)}{{\left[{\left(1+z\right)}^4+{\left(bcz\right)}^2\right]}^{1.5}}\}$(2) of the curvature in the ascending and descending sections was calculated using the logistic regression equation (see Equation 1(1) $f\left(t\right)=\frac{c}{1+\mathrm{e}\mathrm{x}\mathrm{p}\left(a+bt\right)}+d$(1) ). The SOS and LOS were determined by the maximum ROC value in the two sections (Figure 3a). The LOS was determined by subtracting the SOS from the EOS. The R² and RMSE values between the ground-observed and predicted SOS, EOS and LOS using the assimilated LAI were, on average, 65.74% higher and 54.90% lower than those using the smoothed MODIS LAI, respectively (Figure 3b). Therefore, the assimilated LAI time series data could accurately extract the phenological metrics of the bamboo forests.

Figure 3. (a) Phenological metrics of the bamboo forest extracted based on the assimilated LAI time series, (b) comparison between the ground-observed and predicted phenological metrics using the assimilated LAI and smoothed MODIS LAI. ROC represents the curvature value of the LAI time series data fitted by the logistic regression equation. PF_Pre_SOS, PF_Pre_EOS and PF_Pre_LOS represent the predicted SOS-, EOS- and LOS-based assimilated LAI time series, respectively. MODIS_Pre_SOS, MODIS_Pre_EOS and MODIS_Pre_LOS represent the predicted SOS, EOS and LOS based on the smoothed MODIS LAI, respectively.

(1) $f\left(t\right)=\frac{c}{1+\mathrm{e}\mathrm{x}\mathrm{p}\left(a+bt\right)}+d$(1)

(2) $ROC=b^{3}cz\{\frac{3z\left(1-z\right){\left(1+z\right)}^3\left[2{\left(1+z\right)}^3+b^2{c}^{2}z\right]}{{\left[{\left(1+z\right)}^4+{\left(bcz\right)}^2\right]}^{2.5}}-\frac{{\left(1+z\right)}^2\left(1+2z-5z^2\right)}{{\left[{\left(1+z\right)}^4+{\left(bcz\right)}^2\right]}^{1.5}}\}$(2)

where t is the day of the year; f(t) is the LAI value at time t; a, b, c, and d are fitting parameters; $z=e^{a+bt}$; and ROC is the curvature value.

3.3. Statistical analysis

3.3.1. Nonparametric trend analysis

Time series data trend analysis methods can be divided into parametric and nonparametric methods (Gocic and Trajkovic 2013). Although parametric methods are usually used to evaluate the change trend of time series data, these methods require data to conform to a normal distribution. Nonparametric methods (e.g. the Mann-Kendall test and Sen’s slope estimator) are more efficient than parametric methods, do not require samples to follow a certain distribution, and are not affected by a few outliers (Zhu et al. 2022; Zu et al. 2018). Therefore, this study used the Mann-Kendall test and Sen’s slope estimator to analyze the trend change in phenological time series.

(1) Mann-Kendall test

The Mann-Kendall test (Kendall and Kendall 1949; Mann 1945) uses the standard normal test statistic z to detect the change trend of phenology. The standard normal test statistic z is calculated as follows:

(3) $z=\frac{s-1}{\sqrt{Var\left(s\right)}},ifs00,ifs=0\frac{s+1}{\sqrt{Var\left(s\right)}},ifs0$(3)

(4) $s=\sum _{i=1}^{n-1}\sum _{j=i+1}^{n}sgn\left(x_j-x_i\right)$(4)

where s represents the Mann-Kendall test statistic values of phenological time series; n is the number of years; x_i and x_j represent the phenological value of the ith year and jth year (i < j); and Var(s) represents the variance in the Mann-Kendall test statistic s.

Positive z values indicate increasing trends, while negative z values indicate decreasing trends. The change trend of the phenological time series was divided into five levels according to the z value and significance level α = 0.05. Specifically, z > 0 and p < 0.05 represent a significant increasing trend; z > 0 and p > 0.05 represent a nonsignificant increasing trend; z = 0 represents a stable trend; z < 0 and p < 0.05 represent a significant decreasing trend; z < 0 and p > 0.05 represent a nonsignificant decreasing trend.

(2) Sen’s slope estimator

The Sen’s slope estimator first calculates the slope of the trend in a sample of N pairs of data and then ranks them from smallest to largest; the median of the slope is then selected as Sen’s slope (Sen 1968; Gocic and Trajkovic 2013). Sen’s slope (Q_med) is calculated as follows:

(5) $Q_{med}=Q_{\frac{N+1}{2}},ifNisodd\frac{Q_{\frac{N}{2}}+Q_{\frac{N+2}{2}}}{2},ifNiseven$(5)

(6) $Q{\hspace{0.17em}}_k^{\hspace{0.17em}}=\frac{x_j-x_i}{j-i}(k=1,2,\dots ,N;1\le i<j\le n)$(6)

where Q_i represents the slope of ith pairs of phenological data; n is the number of years; and x_i and x_j represent the phenological value of the ith year and jth year.

3.3.2. Linear trend analysis

The slope of the linear regression model (y = ax+b) is typically employed to describe the sensitivity of climate change to phenology. The slope (a) and the uncertainty (σ) of the linear regression model are calculated as follows:

(7) $a=\frac{n\sum _{i=1}^n{x}_i{y}_i-\sum _{i=1}^n{x}_i\sum _{i=1}^n{y}_i}{n\sum _{i=1}^n{x}_i^2-{\left(\sum _{i=1}^n{x}_i\right)}^2}$(7)

(8) $\sigma =\sqrt{\frac{1}{n-2}(\sum _{i=1}^n{(y_i-\bar{y})}^2-a\sum _{i=1}^n(x_i-\bar{x})(y_i-\bar{y}))}$(8)

where a is the sensitivity value of climate change to phenology; σ represents the uncertainty of the linear regression model; n is the number of years; x_i represents the ith year of climate change, and $\bar{x}$ represents the mean of x_i; and y_i represents the phenological value in the ith year, and $\bar{y}$ represents the mean of y_i.

3.3.3. Driving force analysis

(1) Partial correlation coefficient

The partial correlation coefficient (PCC) was employed to assess the reaction of the phenological metrics (SOS, LOS and EOS) to Pre, Rad, Tmax, and Tmin to explore the response of vegetation phenology to a single climatic component. The PCC ($\rho$) is calculated as follows (Wu et al. 2015; Ge et al. 2021a):

(9) ${\rho }_{12\bullet 3\text{~}p}=\sqrt{\frac{R_{1(2,3,\dots ,p)}^2-R_{1(3,\dots ,p)}^2}{1-R_{1(3,\dots ,p)}^2}}$(9)

where ${\rho }_{12\bullet 3\sim p}$ represents the PCC of variables 1 and 2; $R_{1(2,3,\dots ,p)}^2$ is the coefficient of determination of the regression analysis for variable 1 and variables (2, 3, … , p); and $R_{1(3,\dots ,p)}^2$ represents the coefficient of determination for variable 1 and variables (3, … , p) in the regression analysis. For instance, the PCC of the SOS and the Pre can be calculated with the control variables Rad, Tmax, and Tmin. ${\rho }_{12\bullet 3\text{~}p}$ refers to the PCC of the SOS and Pre variables; $R_{1(2,3,\dots ,p)}^2$ is the coefficient of determination of the regression of the SOS variable and the Pre, Rad, Tmax and Tmin variables; and $R_{1(3,\dots ,p)}^2$ represents the coefficient of determination of the regression of the SOS variable and the Rad, Tmax and Tmin variables.

(2) Structural equation model

We constructed a path model to clarify the influence mechanism of climate change on phenology. The path model is a special structural equation model without latent variables (LVs) that can analyze the linear relationship between multiple independent and dependent variables and address more complex variable relationships than regression analysis. The absolute value of the path coefficient may indicate the significance of the independent variable in relation to the dependent variable.

The partial least squares path model (PLS-PM) was used to examine the direct and indirect impacts of the driving elements of phenology to fully understand how topographical and climatic change factors affected phenological metrics under land cover/use change. The PLS-PM model is generally called the PLS method of structural equation modeling, which is used to analyze the linear statistical relationship between multiple variables and can effectively solve the problem of multivariable complex collinearity. The PLS-PM model does not demand that the variables adhere to specific distributional assumptions, and it can still be used in cases of multiple correlations of variables or small sample sizes (Jakobowicz and Derquenne 2007; Chin and Newsted 1999; Sarstedt, Ringle, and Hair 2022). This method can overcome problems such as complex distributions of variable data and correlations between variables.

The two components of the PLS-PM model are measurements and structural models (Fernandes et al. 2018; Li et al. 2020). The link between the measured variable (MV) and the LV is described by the measurement model using the weight and correlation coefficient. The link between the LVs in the structural model is described by the path coefficient (PC). The LVs (including phenological factors, climate factors, topographic factors and abundance factors) cannot be directly observed but are indirectly described by a block of observable variables called MVs (including the SOS, EOS, LOS, Pre, Rad, Tmax, Tmin, altitude, slope, aspect and abundance). For example, the phenological factor of the LV can be calculated from the weighted SOS, EOS and SOS values; the topographic factor of the LV can be calculated from weighted altitude, slope and aspect values. The level of impact of each MV over its associated LV is indicated by its weight. The correlation coefficient measures how closely each MV and its accompanying LV are correlated. The amount of MV information contained in the LV increases with the correlation coefficient. The PLS-PM model was implemented in the “plspm” package in R (version 4.2.0).

4. Results

4.1. Spatiotemporal patterns of bamboo forest phenological metric

Figure 4 displays the geographical distribution of the average annual values and the Sen’s slope of the bamboo forest phenological metrics, and the significant spatial trends of the phenology based on the Mann-Kendall test are shown in Appendix S2. Most of the bamboo forest in Zhejiang Province is in the northwest and center of the province, with 49.84% of the total distribution represented by an average SOS between 90 and 100 days (Figure 4a). Area with an average SOS <70 days accounted for 4.82% (Figure 4a). A total of 20.51% and 14.60% of the pixels had significant (p < 0.05) increasing (delaying) and decreasing (advancing) trends, respectively (Appendix S2). As shown in Figure 4b, the temporal trends for the SOS indicated a delay from north to south. Areas with changes in the SOS between −2 and 0 days per year accounted for 64.38% and were mainly dispersed in the north and the center. Areas with changes in the SOS >0 days occupied 26.77% and were mainly distributed in southern Zhejiang Province.

Figure 4. The spatial patterns of the average annual (a) SOS, (c) EOS, and (e) LOS and the Sen’s slope of (b) SOS, (d) EOS, and (f) LOS in Zhejiang Province.

The average EOS showed a delay from north to south (Figure 4c). An EOS between 310 and 320 days accounted for 40.97% of the area, mainly in the west and south. A total of 26.68% and 9.78% of the pixels had significant (p < 0.05) delayed and advancing trends, respectively (Appendix S2). Sen’s slope of the EOS increased from north to south, as shown in Figure 4d. Areas with changes in the EOS between 0 and 2 days accounted for 74.17% and were mainly in the north and center of Zhejiang Province.

The average LOS increased from north to south (Figure 4e). Areas with an LOS greater than 220 days accounted for 38.09% and were mainly in the northeast and south. Areas with an LOS between 200 and 220 days accounted for 37.60%. A total of 31.45% and 8.62% of the pixels had significant (p < 0.05) increasing and decreasing trends, respectively (Appendix S2). The LOS increased from north to south (Figure 4f). An estimated 61.90% of the northern and central research locations had LOS changes between 0 and 2 days per year.

The average SOS increased significantly by 0.81 days per year (p < 0.01) in the 17 years of the study, indicating that the SOS advanced significantly (Figure 5), while the EOS and LOS significantly increased by 0.27 and 1.08 days per year, respectively (p < 0.01), from 2001–2017, indicating that the EOS and LOS were significantly delayed (Figure 5).

Figure 5. The change trends in the average SOS, EOS and LOS in Zhejiang Province from 2001–2017. The shadow area represents the 95% confidence intervals of the linear fitting curve of phenological annual change.

4.2. Influence of climate change on bamboo forest phenology

The PCCs of the phenological metrics and climate factors (Pre, Rad, Tmax and Tmin) had obvious spatial heterogeneities (Figure 6), and the significance values of the PCCs between the phenological and climatic variables are shown in Table 1. The PCC between the SOS and Pre changed from negative to positive from north to south (Figure 6a). The SOS displayed a negative PCC with Pre in 53.32% of the bamboo distribution area (3.16% of them were significant with p < 0.05), with partial correlations (< −0.4) in 25.43% of all pixels located in the northern areas of Zhejiang Province. Positive PCCs >0.4 were found in 21.66% of the pixels. The PCCs between the EOS and LOS and the Pre exhibited similar spatial patterns (Figures 6(b,c)). There were positive partial correlations between the EOS and LOS and the Pre in 58.10% (6.23% of which were significant at p < 0.05) and 56.71% (4.49% of which were significant at p < 0.05) of all pixels, respectively, which were mainly distributed in northern Zhejiang Province. Negative correlations (< −0.4) were found in 21.72% and 21.09% of the pixels, respectively.

Figure 6. Spatial variations in PCCs between (a) SOS and Pre, (b) EOS and Pre, (c) LOS and Pre, (d) SOS and Rad, (e) EOS and Rad, (f) LOS and Rad, (g) SOS and Tmax, (h) EOS and Tmax, (i) LOS and Tmax, (j) SOS and Tmin, (k) EOS and Tmin, and (l) LOS and Tmin in Zhejiang Province from 2001–2017.

Table 1. The percentages of significant PCCs between the phenological metrics and climate variables (+ represents positive significance; - represents negative significance).

Proportion of significance		Pre	Rad	Tmax	Tmin
SOS	+	2.17%	1.75%	2.50%	2.01%
	-	3.16%	2.97%	2.00%	2.71%
EOS	+	6.23%	4.24%	2.81%	2.78%
	-	2.29%	2.05%	2.95%	2.65%
LOS	+	4.49%	3.57%	2.11%	2.68%
	-	2.18%	1.65%	2.66%	2.14%

The partial correlations between the SOS and Rad exhibited obvious differences from north to south (Figure 6d). The SOS exhibited negative PCCs with Rad in 55.38% of the bamboo distribution area (2.97% of the pixels were statistically significant; p < 0.05), with PCCs (< −0.4) in 26.65% of all the pixels in northwestern Zhejiang Province. Positive correlations (>0.4) were found in 19.52% of the pixels. The partial correlations between the EOS and LOS and the Rad exhibited similar spatial patterns (Figures 6(e,f)). There were positive PCCs >0.4 between the EOS, LOS and Rad, accounting for ~50% of the areas. Negative PCCs between the EOS, LOS and Rad were found on the southwestern and southeastern coasts, and negative correlations (< −0.4) were found in 19.78% and 18.64% of the pixels, respectively.

The SOS had positive PCCs with Tmax in 51.38% of the pixels (2.50% of which were significant at p < 0.05), with PCCs >0.4 in 23.03% of the pixels located in northwestern Zhejiang Province (Figure 6g). Negative PCCs < −0.4 were found in 21.93% of the pixels. The partial correlations between the EOS, LOS and Tmax exhibited similar spatial patterns (Figures 6(h,i)). The positive PCCs between the EOS, LOS and Tmax in 48.74% (2.81% of which were significant at p < 0.05) and 48.56% (2.11% of which were significant at p < 0.05) of the pixels were found in northeastern and southeastern Zhejiang Province, with partial correlations >0.4 accounting for 22.44% and 21.90%, respectively. Negative PCCs with < −0.4 between the EOS and LOS and Tmax were found in 23.77% and 23.54% of the pixels, respectively.

The SOS displayed a negative PCC with Tmin in 53.70% of the pixels (2.71% of which showed significant correlation at p < 0.05), and PCCs >0.4 were found in 25.22% of the pixels located in northwestern and southwestern Zhejiang Province (Figure 6j). Positive partial correlations of >0.4 were found in 21.26% of the pixels. The positive PCCs between the EOS and Tmin in 50.16% of the pixels were relatively scattered (Figure 6k). Negative PCCs between the EOS and Tmin were discovered in the middle part, and negative correlations of < −0.4 accounted for 23.31% of the pixels (Figure 6k). The positive PCCs between the LOS and Tmin in 52.85% of the pixels were primarily found in the northeastern and southeastern regions, with the PCCs >0.4 accounting for 25.02% (Figure 6l). Negative PCCs between the LOS and Tmin were found mainly in the northeastern and southeastern regions, and negative correlations of < −0.4 were found in 21.76% of the pixels (Figure 6l).

4.3. The sensitivity of climate change to phenological variations

The spatial uncertainties in the sensitivities of the SOS and EOS to climate factors are shown in Appendix S3. The uncertainties in the sensitivities of SOS and EOS to climate factors (Pre, Rad, Tmax and Tmin) were basically within 10 days. However, the uncertainties in the sensitivity of the SOS to climate factors were slightly higher than those of the EOS to climate factors.

Figure 7 depicts the sensitivities of the SOS and EOS to climate factors (Pre, Rad, Tmax and Tmin) in Zhejiang Province. The sensitivity of the SOS to Pre showed that a 100 mm increase in regional annual precipitation would cause the SOS to advance by 0.18 days. The negative sensitivity values of the SOS to Pre accounted for 55.85% across the study area, and a sensitivity ranging from −0.04 to 0 days mm⁻¹ accounted for 49.41% of the pixels. Positive sensitivities ranging from 0 to 0.03 days mm⁻¹ occupied 37.60% of the pixels (Figure 7a). The sensitivity of the EOS to Pre indicated that a 100 mm increase in the regional annual precipitation would delay the EOS by 0.12 days. The positive sensitivity values of the EOS to Pre accounted for 57.32% of the pixels across the study area, and sensitivities ranging from 0 to 0.02 days mm⁻¹ accounted for 49.41% of the pixels. The sensitivity values of the EOS to Pre mainly ranged from −0.01–0.01 days mm⁻¹, accounting for 64.49% of the distribution area of the bamboo forest (Figure 7b).

Figure 7. Spatial sensitivity change of (a) SOS to Pre, (b) EOS to Pre, (c) SOS to Rad, (d) EOS to Rad, (e) SOS to Tmax, (f) EOS to Tmax, (g) SOS to Tmin, and (h) EOS to Tmin in Zhejiang Province.

The sensitivity of the SOS to Rad showed that an increase in the average solar radiation of 1 W m⁻² would cause the SOS to advance by 0.03 days (Figure 7c). The sensitivity values of the SOS to Rad in the northeast and northwest areas varied from 0 to 2 days/(W m⁻²) for 36.78% of the pixels (Figure 7c). The sensitivity of the EOS to Rad showed that an increase in the average solar radiation of 1 W m⁻² would cause an EOS delay of 0.08 days. The sensitivity values ranged from −1 to 0 days/(W m⁻²) in 31.01% of the pixels and from 0 to 1 days/(W m⁻²) in 38.57% of the pixels, with most sites found within 28°–30° N in Zhejiang Province (Figure 7d).

According to the sensitivity of the SOS to Tmax, a 1°C increase in the average Tmax would cause the SOS to advance by 0.34 days. The negative sensitivity values of the SOS to Tmax accounted for 50.23% of the pixels, and the sensitivity ranged from −10 to 0 days/℃, accounting for 34.24% of the pixels. The sensitivity values ranged from −50 to −10 days/°C, accounting for 15.99% of the pixels, and from 10 to 50 days/°C, accounting for 14.48% of the pixels, and these sites were in southern Zhejiang Province. The sensitivity values of the SOS to Tmax between −5 and 5 days/°C accounted for 43.92% of the pixels, and they were mostly located in northwestern, northeastern and western Zhejiang Province (Figure 7e). The sensitivity of the EOS to Tmax showed that an increase in the average Tmax of 1°C would advance the EOS by 0.30 days. The negative sensitivity values of the EOS to Tmax accounted for 50.24% of the pixels across the study area, and the sensitivity ranged from −5 to 0 days/℃, accounting for 31.32% of the pixels. The sensitivity values ranged from −5 to 5 days/°C in 64.43% of the pixels, which were mostly located within 28°–31.2° N in Zhejiang Province (Figure 7f).

The sensitivity of the SOS to Tmin revealed that a 1°C increase in the average Tmin would cause the SOS to advance by 1.56 days. The sensitivity values ranged from −50 to −15 days/°C and from 15 to 50 days/°C, accounting for 16.60% and 13.20% of the pixels, respectively, which were primarily located in southern Zhejiang Province. The sensitivity values of the SOS to Tmin between −7.5 and 7.5 days/°C accounted for 43.85% of the pixels across the study area, and they were mostly located in northwestern and northeastern Zhejiang Province (Figure 7g). The positive sensitivity values of the EOS to Tmin accounted for 50.17% of the pixels across the study area, and the sensitivity ranged from 0 to 7.5 days/℃, accounting for 34.26% of the pixels. The sensitivity values ranged from −7.5 to 7.5 days/°C in 64.49% of the pixels, and they were mostly located within 28°–31.2° N in Zhejiang Province (Figure 7h).

Figure 8 shows the path coefficients of the phenological metrics and climate factors. Pre and Tmax had significant effects on the SOS, with scores of −0.13 and 0.14, respectively. The path coefficient values (absolute values) of the SOS to climate factors ranked in the order of Tmax > Pre > Rad > Tmin, indicating that Tmax and Pre had the greatest impacts on the SOS and Tmin had the smallest impact. Similarly, the path coefficient values (absolute values) of the EOS to climate factors ranked in the order of Pre > Rad > Tmin > Tmax, indicating that Pre had the greatest impact on the EOS and Tmax had the smallest impact. The path coefficient values (absolute values) of the LOS to climate factors ranked in the order of Pre > Tmax > Rad > Tmin, indicating that Pre had the greatest impact on the LOS and Tmin had the smallest impact. These results showed that Pre was the dominant factor affecting the spatiotemporal variation in bamboo forest phenology and that the impacts of radiation and temperature on phenology cannot be ignored.

Figure 8. Path analysis of the influence of meteorological factors on phenological metrics based on the path model.

*** represents p<0.01.

4.4. Topographical factors impact bamboo forest phenology

Topographical factors (altitude, slope and aspect) affect the phenology of bamboo forests due to the varying terrain of Zhejiang Province. The average SOS, EOS and LOS were calculated at 100-m altitude intervals, and the temporal trends varied across the increasing altitude gradient. There were discrepancies in the phenological trends at different altitudes (Figures 9a–c). The SOS was delayed at <350 m, and this trend increased with altitude, while the LOS was shortened. The advance in the SOS occurred between 350 and 950 m and increased with altitude, while the LOS was lengthened. Increases in altitude caused the SOS to be delayed when the altitude was greater than 950 m, while the LOS was shortened. The EOS was postponed by 0.41 days with every 100-m increase in altitude. Moreover, the PCCs between the SOS, EOS and Pre along with altitude were highest, while the LOS had the highest partial correlation coefficient with Tmax (Figure 10a), which indicated that altitude mainly affected phenological changes by affecting Pre and Tmax.

Figure 9. Change trends in (a) SOS, (b) EOS, and (c) LOS with altitude, (d) SOS, (e) EOS, (f) LOS with slope, (g) SOS, (h) EOS, and (i) LOS with aspect in Zhejiang Province. The shadow area of (a)-(f) represents the 95% confidence interval (CI) of the fitting curve.

Figure 10. The PCCs between phenology and climatic factors along with (a) altitude, (b) slope and (c) aspect in Zhejiang Province.

The average SOS, EOS and LOS were calculated for every 1° of slope, and the temporal trends varied across the increasing slope gradient. There were discrepancies in the phenological trends at different slopes (Figures 9(d–f)). The SOS was delayed at 0-5°, and this delay increased with slope, while the trends of the EOS and LOS were opposite to that of SOS. The advancing SOS trend between 5° and 25° increased with slope, while the LOS was lengthened. The delay in the EOS at 5-15° increased with slope, reaching the highest value of 312 days at 15°, but the advancing trend weakened at 15-25°. There was no obvious trend between the SOS, EOS and LOS and slope at slopes >25°. In addition, the SOS, EOS and LOS had the highest partial correlation coefficients (absolute values) with Tmin, Rad and Pre, respectively (Figure 10b). Therefore, there were obvious differences in how the climatic factors affected the phenological changes with changes in slope.

The average SOS, EOS and LOS were calculated, and the trends varied with aspect (Figures 9(g–i)). The mean SOS tended to be advanced from the north, northeast, east, southeast to south, while the EOS and LOS tended to be delayed and lengthened, respectively. The SOS tended to be delayed in the south, southwest, and west to the northwest, while the EOS and LOS tended to be advanced and shortened, respectively. In addition, the SOS, EOS and LOS had the highest partial correlation coefficients (absolute values) with Pre along with aspect (Figure 10c). Therefore, Pre was the main climatic factor controlling phenological change with aspect.

4.5. Drivers of interannual phenological metric variation

Figure 11 shows the direct and indirect impacts of topographic factors, climatic factors, and bamboo abundance factors on the phenological changes in bamboo forests. The direct effects of climate, topography and abundance factors on the phenological factors were −0.21, 0.49 and 0.47, respectively. The topographic factors indirectly affected the phenological factor by influencing the climate factor; specifically, the indirect influence factor was −0.30, and the total influence factor was 0.19. The abundance factor indirectly affected the phenological factor by influencing the terrain and climate factors; specifically, the indirect influence factor was 0.33, and the total influence factor was 0.80. This result indicated that land cover/use change was the main factor affecting interannual phenological variations. In addition, the SOS, LOS and EOS had higher correlation coefficients (absolute values) with Pre (Figure 12). Thus, topography and abundance mainly affected interannual phenological variations by affecting precipitation and temperature. The SOS and LOS had significantly higher correlation coefficients (absolute values) with altitude, slope and abundance; in contrast, the EOS had a significantly higher correlation coefficient with altitude and slope (Figure 12), which indicated that altitude and slope were the main factors affecting interannual phenological metric variations under land cover/use changes.

Figure 11. Analysis of influences on phenological driving factors based on the PLS-PM model. Circles represent the LVs. Rectangles represent the MVs. Arrows represent the links between MVs and associated LVs, as well as among related LVs, while arrow labels are the correlation coefficients and PCs that quantify those links. Solid arrows and dashed arrows indicate direct and indirect impacts, respectively.

Figure 12. The correlation coefficient between phenological factors and driving factors based on the Pearson correlation coefficient method.

*, ** and *** represent p<0.1, p<0.05 and p<0.01, respectively.

5. Discussion

Bamboo forest is called “the second largest forest in the world” (FAO 2005, 2010) and is a well-known, abundant, environmentally friendly, and renewable resource that makes it a great alternative to wood (Li et al. 2019b). Therefore, it is crucial to comprehend the phenological variations in bamboo forests and their changing correlations with environmental factors to adapt to climate change. With global warming, the phenology of global vegetation exhibits an advancing SOS and a delayed EOS, but different types of plants have different phenologies depending on their location (Myneni et al. 1997; Nezlin, Kostianoy, and Li 2005; Julien and Sobrino 2009; Xia et al. 2019; Stuble, Bennion, and Kuebbing 2021). Myneni et al. (1997) used advanced very-high-resolution radiometer (AVHRR) data to study the vegetation phenology at latitudes of 45–47° N from 1981 to 1991 and found that the SOS advanced by 0.8 days annually and that the EOS was delayed by 0.4 days annually. Julien and Sobrino (2009) used the GIMMS NDVI to monitor evolution patterns in global vegetation phenology from 1981 to 2003 and discovered a yearly advance in the SOS of 0.38 days and a delay in the EOS of 0.45 days. We found a 0.27-day annual delay rate for the EOS that was less than that found by Myneni et al. (1997) and Julien and Sobrino (2009), and it was basically consistent with that of Xia et al. (2019), indicating a 0.22-day annual delay trend in the EOS from 2001 to 2017. These results indicated that the change rate of the EOS is gradually slowing, possibly due to a gradual or an intermittent slowdown of global warming.

The uncertainties in bamboo forest phenology extraction based on remote sensing may be caused by the LAI data assimilation algorithm. The special growth characteristics of bamboo are different from those of other forest types in spring. Bamboo grows swiftly from shoots to new bamboo stalks in approximately 50 days during on-years. Old bamboo leaves deteriorate and fall off, while new leaves emerge during off-years (Chen 2010). Thus, LAI data assimilation based on remote sensing data every 8 days may insufficiently depict the phenological characteristics of on- and off-years, so it is impossible to accurately depict the growth process of bamboo forests, resulting in spatiotemporal uncertainty in phenological metric extraction. In addition, the LAI data assimilation algorithm does not take into account the impacts of management activities, extreme climate change and other factors on the growth of bamboo forests, which could also result in uncertainty in phenology estimation, leading to uncertainty in the spatial trends of bamboo forest phenology.

The spatial trends of bamboo forest phenology have obvious heterogeneity (Figure 4) and are possibly affected by climate, terrain factors and bamboo forest abundance (Wu et al. 2021; Xia et al. 2019; Stuble, Bennion, and Kuebbing 2021; Silveira et al. 2019; Ren, Qin, and Ren 2019; Chen et al. 2018). Many studies have shown that climate factors are the dominant factors affecting phenological changes (Wang et al. 2019b; Xu, Lu, and Yu 2005; Wang, Luo, and Shafeeque 2019; Wang et al. 2019a; Ge et al. 2021b). The PCC was employed to examine the effects of climate conditions on the SOS, EOS, and LOS to illustrate how phenology responds to climate change. We discovered that there were glaring discrepancies in the geographical patterns of climatic variable impacts on the SOS, EOS, and LOS (Figure 6). The effects of climatic conditions on phenology varied within the same location; for instance, Pre had a negative influence on the SOS but a positive impact on the EOS and LOS in northwestern Zhejiang Province (Figure 6). We also discovered that the PCCs and path coefficients (absolute value) between Pre and phenology were higher, demonstrating that Pre was the dominant factor affecting the spatiotemporal variation in bamboo forest phenology (Figures 6 and 8). This result conflicted with those of Piao et al. (2006), Song et al. (2011), Wang et al. (2019b) and Xia et al. (2019), who found that precipitation was not the dominant factor affecting phenology; however, our result was in line with that of Verger et al. (2016), who found that vegetation phenology and rainfall were obviously correlated at low latitudes.

Additionally, we found that the uncertainties in the sensitivity of the SOS to climate factors were slightly higher than those of the sensitivity of the EOS to climate factors (Appendix S3), which may be because bamboo forest is sensitive to climate at the beginning of the growing season. This is because bamboo shoots emerge and grow at a very rapid rate in the spring, reaching the maximum water and temperature demand period of the bamboo growth process (Xu and Tan 1997; Zhang 1995; Li et al. 2019b). The climate demand of bamboo forests differs in different growth stages. Bamboo shoot tips need certain climatic conditions when exposed to the surface. Previous studies have shown that the minimum amount of precipitation and temperature suitable for bamboo growth are 337 mm and 8°C in spring (Li et al. 2019b; Chen 2009). The suitable temperature for the rapid growth of bamboo culm is 17–19°C, and the suitable temperature for the growth of branches and leaves is 18–25°C; meanwhile, a seasonal precipitation distribution is required (Chen 2009). Thus, the amount of climatic conditions, whether daily, weekly, or monthly, has a great influence on the development of bamboo forests. These may be the reasons for the relatively large uncertainty in the sensitivity of the SOS to climate factors.

Terrain conditions led to differences in the climate conditions of the study area, causing different altitudes, slopes and aspect gradients to show different phenological characteristics (Figure 9). Previous studies have shown that altitude has important impacts on the SOS (Colombo et al. 2009; Wang, Zhang, and Rodman 2021; Hua, Ohlemüller, and Sirguey 2022). Our results also indicated that the altitude gradient had a substantial impact on the bamboo forest phenology, and the shift trend for the phenology that considered topographic parameters was not constant, which concurred with the outcomes of Yuan et al. (2020) and Xia et al. (2019). Altitude and aspect are connected to precipitation and temperature, and slope is associated with erosion and land management (Li et al. 2017). The effects of altitude, slope and aspect on bamboo forest phenology were quite different (Figure 9). Phenology had a significant correlation with altitude and slope, and altitude and slope had the highest correlation coefficients with Pre and Tmin (Figure 12). Therefore, topographic factors affect and restrict the interannual change in bamboo forest phenology to a certain extent through precipitation and temperature.

Bamboo forest phenology could also be affected by abundance. First, a phenology pixel may contain an unknown composition of vegetation species or land cover/use types because of the low spatiotemporal resolution of MODIS products, resulting in a possibly inaccurate depiction of the phenological changes in bamboo forests in remote sensing-based phenology. Thus, a low abundance would cause the uncertainty of vegetation phenology estimates to a certain extent (Chen et al. 2018). Second, the influences of climate and topographical factors on phenology were obviously different at different class levels of bamboo forest abundance (Appendix S4). The bamboo forest abundance increased from 0–0.8 to 0.8–1, and the relatively important climate factors influencing the SOS changed from Tmax and Pre to Tmin and Pre. Pre and Rad had a great influence on the EOS at an abundance value of 0.2–1. In terms of topographical factors, the bamboo forest abundance increased from 0.2–0.6 to 0.6–1, and the relatively important topographic factor influencing the SOS changed from altitude to slope. These results indicate that bamboo forest abundance had a greater impact on phenological variation through an indirect influence on climate and topography. Additionally, bamboo forests incur frequent business activities (Mao et al. 2016), and human activities, extreme drought and rainfall, and other natural disasters will cause great fluctuations in the interannual changes in bamboo forest phenology. However, this study ignored how human management (e.g. fertilization, felling) and extreme climate affect phenology. Therefore, how to quantify the different impacts of extreme climate, terrain, human activities and land cover/use type changes on bamboo forest phenology should be addressed in future research.

6. Conclusions

Bamboo forest phenology was analyzed to determine the impacts of climatic and topographical variables on the SOS, EOS and LOS of bamboo forests under land cover/use change during 2001–2017. The results revealed that the spatial distribution and change trends of the bamboo forest phenology had obvious heterogeneity. The average annual advance in the SOS, delay in the EOS and lengthening of the LOS were 0.81 days, 0.27 days and 1.08 days, respectively. There were obvious spatial differences in the partial correlation coefficients between climate factors and phenological metrics. Precipitation was the main factor affecting the bamboo forest phenology across the distribution area of the bamboo forest, and the average advance in the SOS and delay in the EOS were 0.18 days and 0.12 days, respectively, for every 100-mm increase in regional annual precipitation. Terrain conditions caused differences in climate conditions in the study area, with different altitudes, slopes and aspect gradients showing different phenological characteristics. Topography and abundance mainly affected interannual phenological variations by affecting precipitation and temperature. Therefore, this study clarified the interactive response mechanism of bamboo forest phenology to climate change and geographical conditions under the condition of land cover/use change, thus providing a theoretical basis for further evaluating the carbon sink potential of bamboo forest. However, the study ignored influence mechanism of human management and extreme climate on phenology. We will accurately evaluate the interactive effects and contributions of extreme climate, terrain, human activities and land cover/use type changes on bamboo forest phenology in future research.

Acknowledgments

The research was supported by the National Natural Science Foundation of China (No. 32201553, 32171785, U1809208, 31901310), Talent launching project of scientific research and development fund of Zhejiang A & F University (No. 2021LFR029).

Disclosure statement

The authors declare no possible conflicts of interest.

Data availability statement

The MODIS, Climate and DEM datasets used during this study are openly available from NASA’s website (https://ladsweb.modaps.eosdis.nasa.gov/search/), the National Meteorological Science Data Center of China (http://data.cma.cn/), and the Geospatial Data Cloud site, the Computer Network Information Center, and the Chinese Academy of Sciences (http://www.gscloud.cn), respectively.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15481603.2022.2163575

Additional information

Funding

The work was supported by the National Natural Science Foundation of China [32201553]; National Natural Science Foundation of China [31901310]; National Natural Science Foundation of China [U1809208]; Talent launching project of scientific research and development fund of Zhejiang A & F University [2021LFR029]; National Natural Science Foundation of China [32171785].