Populus euphratica Phenology and Its Response to Climate Change in the Upper Tarim River Basin, NW China

: Quantifying the phenological variations of Populus euphratica Olivier ( P. euphratica ) resulting from climate change is vital for desert ecosystems. There has previously been great progress in the inﬂuence of climate change on vegetation phenology, but knowledge of the variations in P. euphratica phenology is lacking in extremely arid areas. In this study, a modiﬁed method was proposed to explore P. euphratica phenology and its response to climate change using 18-year Global Land Surface Satellite (GLASS) leaf area index (LAI) time series data (2000–2017) in the upper Tarim River basin. The start of the growing season (SOS), length of the growing season (LOS), and end of the growing season (EOS) were obtained with the dynamic threshold method from the reconstructed growth time series curve by using the Savitzky–Golay ﬁltering method. The grey relational analysis (GRA) method was utilized to analyze the inﬂuence between the phenology and the key climatic periods and factors. Importantly, we also revealed the positive and negative effects between interannual climate factors and P. euphratica phenology using the canonical correlation analysis (CCA) method, and the interaction between the SOS in spring and EOS in autumn. The results revealed that trends of P. euphratica phenology (i.e., SOS, EOS, and LOS) were not signiﬁcant during the period from 2000–2017. The spring temperature and sunshine duration (SD) controlled the SOS, and the EOS was mainly affected by the temperature and SD from June–November, although the impacts of average relative humidity (RH) and precipitation (PR) on the SOS and EOS cannot be overlooked. Global warming may lead to SOS advance and EOS delay, and the increase in SD and PR may lead to earlier SOS and later EOS. Runoff was found to be a more key factor for controlling P. euphratica phenology than PR in this region. P. euphratica phenology; (2) identify the key meteorological factors and key periods that affect the phenology of P. euphratica , and (3) conﬁrm whether the beginning and the end of the P. euphratica growing season inﬂuence each other. The results provide insights into respectively. The typical variable V 1 had positive correlations with average RH with the typical load coefﬁcient of 0.24. V 1 had negative correlations with PR, average ST, average AT, EA, and SD, however, with typical load coefﬁcients of − 0.26, − 0.22, − 0.63, − 0.87, − 0.17, respectively. These results generally indicated that the SOS had negative correlations with PR, SD, and temperature (i.e., AT, ST, and EA), and a positive correlation with RH. The EOS had positive correlations with PR and temperature, and negative correlations with relative humidity. factors, factors on vegetation phenology change, which has not been adequately considered in current vegetation phenology research and ecosystem models. This study also revealed that global P. with in of the Tarim The results provide insights into quantitative assessment of spatio-temporal variations in P. euphratica phenology under climate change and provide a scientiﬁc reference for future management of the desert ecosystem.


Introduction
Vegetation phenology is the research of events in the plant life cycle, and how these events adapt to climate changes [1,2], which makes it the most intuitive and sensitive biological indicator of climate change. Climate change, such as modified precipitation, increasing temperature, and increasing frequency of extreme weather events, will change the vegetation phenology [3], which will have extensive effects on plant community structure, plant distribution, energy cycles, vegetation ecosystems, and primary production of vegetation [4][5][6][7][8]. Altered phenological timing, furthermore, may have negative effects on species interactions in the food chain, such as migration, survival, reproduction, and occupying feeding habits [9,10]. Therefore, a better understanding of how vegetation phe-nology responds to climate change can improve the accuracy of phenological simulation models and provide a scientific reference for future management of vegetation ecosystems.
Plant growth, including leaf unfolding, flowering, fruiting, and leaf senescence, has been recorded to reflect the phenological characteristics of vegetation based on ground monitoring sites of phenology. However, this method has some limitations in the spatiotemporal scale and biome scale [11]. Therefore, the remote sensing technique for phenological monitoring, which is used to simulate various logistic regression models through long-term series remote sensing data and to extract vegetation phenology (i.e., the start of the growing season (SOS), length of the growing season (LOS), and end of the growing season (EOS)) from suitable models, has been widely applied to explore the spatial-temporal variations in vegetation phenology at the regional and global scales [12][13][14][15][16][17]. However, there is a wide range of mixed pixels in satellite remote sensing data, and the spatial resolution of remote sensing data will affect the accuracy of vegetation phenological feature inversion, especially in the extreme arid areas where the vegetation distribution is more sparse and the coverage is lower. In addition, previous studies on vegetation phenology have focused on plant phenology change trends [18][19][20], phenology changes in different vegetation types [21], temperature [22][23][24], rainfall [25][26][27], aspect [28], altitude [23], solar radiation [29], and other single climatic factors affecting vegetation phenology. Therefore, there is still a need to demonstrate the relationship of the preseason, interannual, and multi-climatic factors with phenology, and what kind of interaction exists between the SOS and EOS.
Arid and semi-arid ecosystems are important parts of Earth's climate system, covering 41% of the Earth's land surface [30][31][32]. These regions are characterized by harsh environments and little rainfall, especially in extremely arid areas, where the average annual rainfall is often <60-100 mm [33]. The Tarim River basin is located in the hinterland of Eurasia, an extremely arid area. This region is highly vulnerable and sensitive to climate change due to the particular features comprising mountains and a mosaic of innermountain basins [34]. Populus euphratica (Salicaceae), an ancient, precious, endangered species, is a constructive tree in extremely arid areas [34]. P. euphratica is a vital part of the extreme drought ecosystem in this region due to its advantageous characteristics such as cold resistance, heat resistance, drought resistance, wind-sand resistance, and salt-alkali resistance [35]. In the upper Tarim River basin, since there it is less affected by the Ecological Water Conveyance Project than the lower Tarim River, the results on the response of P. euphratica phenology to climate change are more believable. However, research in this field mainly focuses on the relationship between vegetation and groundwater [34,36], changes in rainfall and temperature under climate change and global warming, its impact on water resources [37][38][39][40][41], and the impact of water on physical and chemical indicators of vegetation growth [30,32,42,43]. To date, few, if any, studies have revealed the dynamics of P. euphratica phenology and its response to climate change in arid areas.
The novelty of this study is that a modified method was proposed to explore P. euphratica phenology and its response to climate change using 18-year Global Land Surface Satellite (GLASS) leaf area index (LAI) time series data (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017) in the upper Tarim River basin. A sandwich spatial sampling strategy, which can reduce the accuracy error caused by many mixed pixels in satellite remote sensing data, was used to extract the LAI information of P. euphratica from the GLASS LAI time series data, and the phenological information (i.e., SOS, LOS, EOS) was obtained with the dynamic threshold method from the reconstructed growth time series curve by using the Savitzky-Golay filtering method. The grey relational analysis (GRA) and canonical correlation analysis (CCA) methods were utilized to explore the relationship between phenology and multi-climatic factors. Specifically, the main objectives of this study were to (1) quantitatively reveal the change trends of P. euphratica phenology; (2) identify the key meteorological factors and key periods that affect the phenology of P. euphratica, and (3) confirm whether the beginning and the end of the P. euphratica growing season influence each other. The results provide insights into quantitative assessment of interannual variations of P. euphratica phenology under climate change and provide a scientific reference for future management of the desert ecosystem.

Study Area
Through field investigation, P. euphratica was found to be mainly distributed along the upper reaches of the Tarim River and the Yerqiang River. Therefore, this investigation mainly focused on the upper Tarim River and the Yeerqiang River basin, which was divided into two regions ( Figure 1b): Bachu (BC) and Shaya (SY). The study area experiences a typical temperate arid continental climate, with long sunshine duration, strong annual evaporation, and low annual rainfall. The extreme minimum temperature is −25 • C, the extreme maximum temperature is 39.4 • C, and the annual average temperature is 10.4 • C [34]. The annual average precipitation is 30-50 mm, and the annual potential evaporation is 2000-2900 mm [43]. the end of the P. euphratica growing season influence each other. The results provide insights into quantitative assessment of interannual variations of P. euphratica phenology under climate change and provide a scientific reference for future management of the desert ecosystem.

Study Area
Through field investigation, P. euphratica was found to be mainly distributed along the upper reaches of the Tarim River and the Yerqiang River. Therefore, this investigation mainly focused on the upper Tarim River and the Yeerqiang River basin, which was divided into two regions ( Figure 1b): Bachu (BC) and Shaya (SY). The study area experiences a typical temperate arid continental climate, with long sunshine duration, strong annual evaporation, and low annual rainfall. The extreme minimum temperature is −25 °C, the extreme maximum temperature is 39.4 °C, and the annual average temperature is 10.4 °C [34]. The annual average precipitation is 30-50 mm, and the annual potential evaporation is 2000-2900 mm [43].

Data
The Global Land Surface Satellite (GLASS) leaf area index (LAI) has a temporal resolution of eight days and spans the period from 1981-2017. The LAI product from 1982-1999 was developed from AVHRR reflectance with a spatial resolution of 0.05°, and that from 2000-2017 with a spatial resolution of 1 km was derived from MODIS reflectance [44]. The GLASS LAI data product is a national key project, and its dataset is superior to other similar products in terms of spatial and temporal accuracy [45]. This product feature has significant advantages in the continuity of its time series and the integrity of its spatial scope [46]. Moreover, the GLASS LAI dataset has passed the verification analysis of the measurement data of the common international site, with the results demonstrating that the verification accuracy is higher than other similar products [47,48]. Since P. euphratica has a small distribution area and low vegetation coverage, this study used the GLASS 1 km resolution data, which were from the National Earth System Science Data Sharing Infrastructure, National Science & Technology Infrastructure of China (http://www.geodata.cn (14 October 2019)). A total of eighteen years of data from 2000-2017 were obtained, and all the LAI values were the original values of the GLASS LAI dataset, which were equal to 10 times the actual values.

Data
The Global Land Surface Satellite (GLASS) leaf area index (LAI) has a temporal resolution of eight days and spans the period from 1981-2017. The LAI product from 1982-1999 was developed from AVHRR reflectance with a spatial resolution of 0.05 • , and that from 2000-2017 with a spatial resolution of 1 km was derived from MODIS reflectance [44]. The GLASS LAI data product is a national key project, and its dataset is superior to other similar products in terms of spatial and temporal accuracy [45]. This product feature has significant advantages in the continuity of its time series and the integrity of its spatial scope [46]. Moreover, the GLASS LAI dataset has passed the verification analysis of the measurement data of the common international site, with the results demonstrating that the verification accuracy is higher than other similar products [47,48]. Since P. euphratica has a small distribution area and low vegetation coverage, this study used the GLASS 1 km resolution data, which were from the National Earth System Science Data Sharing Infrastructure, National Science & Technology Infrastructure of China (http://www.geodata.cn (14 October 2019)). A total of eighteen years of data from 2000-2017 were obtained, and all the LAI values were the original values of the GLASS LAI dataset, which were equal to 10 times the actual values.
This study used climatic data, including sunshine duration (SD), average air temperature (AT), average land surface temperature (ST), effective accumulated temperature (EA), average relative humidity (RH), and precipitation (PR). These data were derived from http://data.cma.cn/ (14 July 2021) from 1999-2018. Annual runoff data were derived . Additionally, in this study, based on Landsat TM/OLI long time series remote sensing data with 30 m resolution from 1990-2015, the distribution data of P. euphratica for 1990, 1995, 2000, 2005, 2010, and 2015 were extracted using the object-oriented classification method. These data were then used to select typical sample sites of P. euphratica.

Methods
The technique flowchart of this study is shown in Figure 2, the research steps were as follows. (1) P. euphratica samples were selected and extracted from LAI time series data using the sandwich spatial sampling method. The growth time series curve of P. euphratica was then reconstructed using the Savitzky-Golay filtering method. (2) The phenological features of P. euphratica, including the SOS, EOS, and LOS, were extracted using the dynamic threshold method. (3) The grey relational analysis (GRA) method was used to analyze the correlations between climatic factors and the SOS during three periods (spring, SOS preseason, and annual), and the EOS during three periods (autumn, EOS preseason, and annual). (4) The canonical correlation analysis (CCA) method was used to explore the positive and negative effects between phenology (SOS, EOS) and interannual climatic factors (SD, AT, ST, EA, RH, and PR). This study used climatic data, including sunshine duration (SD), average air temperature (AT), average land surface temperature (ST), effective accumulated temperature (EA), average relative humidity (RH), and precipitation (PR). These data were derived from http://data.cma.cn/ (14 July 2021) from 1999-2018. Annual runoff data were derived from the Statistical Bulletin of the Ministry of Water Resources of the People's Republic of China (http://www.mwr.gov.cn/sj/#tjgb (14 July 2021)). Additionally, in this study, based on Landsat TM/OLI long time series remote sensing data with 30 m resolution from 1990-2015, the distribution data of P. euphratica for 1990, 1995, 2000, 2005, 2010, and 2015 were extracted using the object-oriented classification method. These data were then used to select typical sample sites of P. euphratica.

Methods
The technique flowchart of this study is shown in Figure 2, the research steps were as follows. (1) P. euphratica samples were selected and extracted from LAI time series data using the sandwich spatial sampling method. The growth time series curve of P. euphratica was then reconstructed using the Savitzky-Golay filtering method. (2) The phenological features of P. euphratica, including the SOS, EOS, and LOS, were extracted using the dynamic threshold method. (3) The grey relational analysis (GRA) method was used to analyze the correlations between climatic factors and the SOS during three periods (spring, SOS preseason, and annual), and the EOS during three periods (autumn, EOS preseason, and annual). (4) The canonical correlation analysis (CCA) method was used to explore the positive and negative effects between phenology (SOS, EOS) and interannual climatic factors (SD, AT, ST, EA, RH, and PR).

Extract LAI Features
Based on the distribution data of P. euphratica in 1990, 1995, 2005, 2010, and 2015, the regional composition sample data of P. euphratica, which fully met the 1 km 2 area criterion, were automatically identified by a 1 km × 1 km square window. The typical sample points were selected from high-resolution imagery (e.g., Google Earth imagery, Version 2018) by

Extract LAI Features
Based on the distribution data of P. euphratica in 1990, 1995, 2005, 2010, and 2015, the regional composition sample data of P. euphratica, which fully met the 1 km 2 area criterion, were automatically identified by a 1 km × 1 km square window. The typical sample points were selected from high-resolution imagery (e.g., Google Earth imagery, Version 2018) by the Sandwich Spatial Sampling tool, which is developed based on spatial optimization sampling strategy of Trinity theory. The LAI information was extracted from each sample using the GLASS LAI products from 2000-2017.

Filtering and Reconstruction of LAI Curves
Although the GLASS LAI products have been processed to remove noise, the time series information of the extracted sample points still exhibits the "pseudo" peak phenomenon. Thus, the Savitzky-Golay (S-G) filtering method was utilized to remove the noise and reconstruct the LAI curves of P. euphratica. S-G filtering can be used with the time series index data of different time scales and is not limited by sensor differences, so the noise in the original data can be removed using the minimum root mean square error [49].
where Y j is the resultant LAI value, Y j is the original LAI value, C i is the coefficient for the LAI value of the filter, and N is the number of convoluting integers and is equal to the smoothing window size (2 m + 1).

Dynamic Threshold Method
The dynamic threshold method is widely utilized to extract phenological features of vegetation because it has better spatio-temporal applicability [50,51]. In this method, the SOS is defined as the time when the vegetation LAI increases by a certain percentage, while the EOS is defined as the time when the vegetation LAI decreases by a certain percentage. The threshold used in this method is not a specific value, but a dynamic ratio. When the dynamic threshold method was proposed to extract phenological features, Jönsson and Eklundh [52] suggested that the threshold at the SOS and EOS is 0.2, which has been subsequently used by many scholars in phenological research. However, the choices of threshold values depend on different plant species and different regions, so we had to proceed cautiously and rigorously when selecting thresholds. In this study, we combined remote sensing data with 2016 ground-based monitoring data and used the following procedure illustrated by the flow chart in Figure 3. The SOS and EOS values were obtained using various dynamic thresholds (i.e., A 1 , A 2 ) ranging from 0.1-0.4 at 0.02 intervals. Then, the results were verified by using the ground-based monitoring data (Table 1). In terms of these tests, we opted to use the thresholds of A = 0.2 and A 2 = 0.3 to extract the SOS and EOS of P. euphratica phenology in the study area from 2000 to 2017, respectively.

Grey Relational Analysis (GRA)
The relationship between many factors falls within a gray area, due to the limitations of people's understanding and the complexity of geographical phenomena, and it is difficult to objectively measure the degree of relevance by the relative correlation coefficient method [53]. The grey relational analysis method makes up for the above shortcomings. GRA is not limited by sample type and probability distribution, so it is a new method to investigate uncertain problems under limited information and data [54,55]. Details of the principles and computation steps of the PNPI were presented in Li et al. [53]. In this study, the SOS is used as the reference factor sequences from 2000 to 2017, and the comparative factor sequence consists of climatic factors (i.e., SD, AT, ST, EA, RH, and PR) during three periods (spring, SOS preseason, and annual) from 2000 to 2017. The EOS from 2000 to 2017 is used as the reference factor sequence, and the comparative factor sequence consists of climatic factors during three periods (autumn, EOS preseason, and annual). The correlations between climatic factors and the P. euphratica phenology during different periods were explored.

Grey Relational Analysis (GRA)
The relationship between many factors falls within a gray area, due to the limitations of people's understanding and the complexity of geographical phenomena, and it is difficult to objectively measure the degree of relevance by the relative correlation coefficient method [53]. The grey relational analysis method makes up for the above shortcomings. GRA is not limited by sample type and probability distribution, so it is a new method to investigate uncertain problems under limited information and data [54,55]. Details of the principles and computation steps of the PNPI were presented in Li et al. [53]. In this study, the SOS is used as the reference factor sequences from 2000 to 2017, and the comparative factor sequence consists of climatic factors (i.e., SD, AT, ST, EA, RH, and PR) during three periods (spring, SOS preseason, and annual) from 2000 to 2017. The EOS from 2000 to 2017 is used as the reference factor sequence, and the comparative factor sequence consists of climatic factors during three periods (autumn, EOS preseason, and annual). The correlations between climatic factors and the P. euphratica phenology during different periods were explored.

Canonical Correlation Analysis (CCA)
Canonical correlation is a statistical analysis method used to study the correlation between two groups of variables and was first proposed by Hotelling in 1935 [56]. It is also a dimensionality reduction technology, which is widely used in various fields [57]. Vegetation phenology is not only related to climatic factors, but also the physiological characteristics of vegetation. The change in climatic factors may not lead to one-way phenological changes due to the interaction effects between the SOS in spring and the EOS in autumn [8,58]. It is not rigorous to use simple climatic data to describe these changes. Therefore, the first set of variables (X1) consists of the EOS and SOS from 2000 to 2017 in this study, and the interannual climatic factors (SD, AT, ST, EA, RH, and PR) during the period of 2000-2017 as the second set of variables (Y1). The CCA method not only revealed the relationship between SOS and EOS, but also explored the correlations between the P. euphratica phenology and interannual climate factors. The principle and computation steps are detailed in Hotelling [56].

Canonical Correlation Analysis (CCA)
Canonical correlation is a statistical analysis method used to study the correlation between two groups of variables and was first proposed by Hotelling in 1935 [56]. It is also a dimensionality reduction technology, which is widely used in various fields [57]. Vegetation phenology is not only related to climatic factors, but also the physiological characteristics of vegetation. The change in climatic factors may not lead to one-way phenological changes due to the interaction effects between the SOS in spring and the EOS in autumn [8,58]. It is not rigorous to use simple climatic data to describe these changes. Therefore, the first set of variables (X 1 ) consists of the EOS and SOS from 2000 to 2017 in this study, and the interannual climatic factors (SD, AT, ST, EA, RH, and PR) during the period of 2000-2017 as the second set of variables (Y 1 ). The CCA method not only revealed the relationship between SOS and EOS, but also explored the correlations between the P. euphratica phenology and interannual climate factors. The principle and computation steps are detailed in Hotelling [56].

Interannual Variations of P. euphratica Phenology
Interannual changes and P. euphratica phenology (SOS, LOS, EOS) are shown in Figure 5 in the upper Tarim River basin from 2000-2017. Generally, the upper Tarim River basin showed advancing SOS from 2000-2017 and the EOS of P. euphratica was advanced, but a slight lengthening in LOS was observed in the study area. However, it was found that the interannual variation trends of the EOS, LOS, and SOS derived from GLASS LAI

Interannual Variations of P. euphratica Phenology
Interannual changes and P. euphratica phenology (SOS, LOS, EOS) are shown in Figure 5 in the upper Tarim River basin from 2000-2017. Generally, the upper Tarim River basin showed advancing SOS from 2000-2017 and the EOS of P. euphratica was advanced, but a slight lengthening in LOS was observed in the study area. However, it was found that the interannual variation trends of the EOS, LOS, and SOS derived from GLASS LAI data were not significant in the study time series. It is reasonable that a longer study time series for GLASS LAI data should be considered to research P. euphratica phenology under global climate change. Overall, the spatial heterogeneity of the P. euphratica phenology was higher, and the changing trend of P. euphratica phenology varied in different regions (BC and SY). BC showed delayed EOS, but in SY it was significantly advanced. The P. euphratica phenology Overall, the spatial heterogeneity of the P. euphratica phenology was higher, and the changing trend of P. euphratica phenology varied in different regions (BC and SY). BC showed delayed EOS, but in SY it was significantly advanced. The P. euphratica phenology varied among different years in the same region, and the spatial heterogeneity of P. euphratica phenology in the same year and region (BC or SY) was higher. This is reasonable since climate, which is variable and complex, fundamentally determines the phenology and interannual change trends of phenology.

Effects of Key Climatic Periods and Factors
The GRA between phenology (SOS and EOS) and climatic factors (SD, AT, ST, EA, RH, and PR) is presented in Figure 6. The SOS of P. euphratica was mainly affected by spring climatic factors, while the EOS of P. euphratica was mainly affected by the EOS preseason (June-November) meteorological factors. The SOS of P. euphratica was mainly affected by the climatic factors (0.81) in spring (March-May) (Figure 6a), and less affected by the annual (January-December) climatic factors (0.75) or the SOS preseason (December-May) climatic factors (0.67) (Figure 6c,e). Moreover, the influence of spring climatic factors in the different regions (BC and ALE) on the SOS of P. euphratica was greater than that of the annual climatic factors or the SOS preseason climatic factors.
varied among different years in the same region, and the spatial heterogeneity of P. euphratica phenology in the same year and region (BC or SY) was higher. This is reasonable since climate, which is variable and complex, fundamentally determines the phenology and interannual change trends of phenology.

Effects of Key Climatic Periods and Factors
The GRA between phenology (SOS and EOS) and climatic factors (SD, AT, ST, EA, RH, and PR) is presented in Figure 6. The SOS of P. euphratica was mainly affected by spring climatic factors, while the EOS of P. euphratica was mainly affected by the EOS preseason (June-November) meteorological factors. The SOS of P. euphratica was mainly affected by the climatic factors (0.81) in spring (March-May) (Figure 6a), and less affected by the annual (January-December) climatic factors (0.75) or the SOS preseason (December-May) climatic factors (0.67) (Figure 6c,e). Moreover, the influence of spring climatic factors in the different regions (BC and ALE) on the SOS of P. euphratica was greater than that of the annual climatic factors or the SOS preseason climatic factors.   (Figure 6b,d). Moreover, the influence of autumn climatic factors on the EOS of P. euphratica was greater than that of the annual climatic factors or the EOS preseason meteorological factors in the BC and SY regions. It indicated that spring climatic factors (March-May) influence the SOS of P. euphratica more than climatic factors in other periods (i.e., annual and preseason), and the EOS is more affected by the variations in climatic factors from June to November.
Different climatic factors had different effects on phenology ( Figure 6). The effects of average ST (0.91), average AT (0.90), average EA (0.89), and SD (0.86) on the SOS of P. euphratica were greater than those of average RH (0.82) and precipitation (0.47). The effects of average AT (0.86), SD (0.86), average EA (0.85), and average ST (0.84) on the EOS of P. euphratica were greater than those of average RH (0.79) and precipitation (0.39). Therefore, we could infer that temperature (i.e., AT, ST, and EA) and SD controlled P. euphratica phenology, although the impacts of RH and PR on P. euphratica phenology cannot be overlooked in the upper Tarim River basin. between the EOS and the annual (January-December) climatic factors, (e) GRA between the SOS and the SOS preseason (December-May) climatic factors, (f) GRA between the EOS and the EOS preseason (June-November) climatic factors.

Exploration of the Relationships between SOS, EOS, and Interannual Climate
The EOS of P. euphratica was mainly affected by the EOS preseason (June-November) meteorological factors (0.82) (Figure 6f), and less affected by the annual (January-December) climatic factors (0.78) or the meteorological factors (0.77) in autumn (September-November) (Figure 6b,d). Moreover, the influence of autumn climatic factors on the EOS of P. euphratica was greater than that of the annual climatic factors or the EOS preseason meteorological factors in the BC and SY regions. It indicated that spring climatic factors (March-May) influence the SOS of P. euphratica more than climatic factors in other periods (i.e., annual and preseason), and the EOS is more affected by the variations in climatic factors from June to November.
Different climatic factors had different effects on phenology ( Figure 6). The effects of average ST (0.91), average AT (0.90), average EA (0.89), and SD (0.86) on the SOS of P. euphratica were greater than those of average RH (0.82) and precipitation (0.47). The effects of average AT (0.86), SD (0.86), average EA (0.85), and average ST (0.84) on the EOS of P. euphratica were greater than those of average RH (0.79) and precipitation (0.39). Therefore, we could infer that temperature (i.e., AT, ST, and EA) and SD controlled P. euphratica phenology, although the impacts of RH and PR on P. euphratica phenology cannot be overlooked in the upper Tarim River basin.   The canonical redundancy analysis (Table 2) revealed that the first pair of canonical variables U 1 could explain 53.14% of the intragroup variation and 8.97% of the other group variation. The typical variable V 1 , however, could explain 22.51% of the intragroup variation and 21.20% of the other group variation. Therefore, we could infer that the SOS and EOS of P. euphratica affect each other. The phenology (SOS, EOS) of P. euphratica is not only related to meteorological factors but also the physiological characteristics of P. euphratica. Therefore, it was found that the SOS of P. euphratica was strongly related to EOS with a high canonical correlation coefficient of 0.63, implying that a later SOS may generally be accompanied by a later EOS. Note: X 1 , X 2 , Y 1 , Y 2 , canonical variable.

Validation of P. euphratica Phenology
Remote sensing monitoring of vegetation phenology is based on vegetation growth characteristics (LAI, EVI, NDVI, and other features), and the vegetation phenology is simulated using a mathematical model. This method is different from the traditional ground observation phenology method [11], but previous studies showed that the vegetation phenology phase determined by remote sensing data was closely related to those obtained from the ground phenology monitoring points [59]. This study was based on one year's worth of data from the ground phenology monitoring points of P. euphratica, and established a mathematical model of P. euphratica phenology ( Table 1). The phenological features of P. euphratica, including the SOS, EOS, and LOS, were extracted using the dynamic threshold method. To verify the accuracy of the research results and the applicability of the method, based on the GLASS LAI and higher resolution of MODIS NDVI (MOD09Q1 V6) data in 2006, 2010, and 2014, 25 samples were randomly selected to compare and analyze the results of key phenological phase extraction of P. euphratica through the process flow chart steps in Figure 2. As shown in Figure 8, The results revealed that using the GLASS LAI and MODIS NDVI to extract the phenology of P. euphratica yielded higher correlations (SOS: R 2 = 0.861, p < 0.001; EOS: R 2 = 0.824, p < 0.001). Overall, this method could quickly identify and extract the key phenological period of P. euphratica in a large area, and this study's findings are reliable and have scientific reference value for the study of vegetation phenology, especially in extremely arid areas. tion and 21.20% of the other group variation. Therefore, we could infer that the SOS and EOS of P. euphratica affect each other. The phenology (SOS, EOS) of P. euphratica is not only related to meteorological factors but also the physiological characteristics of P. euphratica. Therefore, it was found that the SOS of P. euphratica was strongly related to EOS with a high canonical correlation coefficient of 0.63, implying that a later SOS may generally be accompanied by a later EOS.

Validation of P. euphratica Phenology
Remote sensing monitoring of vegetation phenology is based on vegetation growth characteristics (LAI, EVI, NDVI, and other features), and the vegetation phenology is simulated using a mathematical model. This method is different from the traditional ground observation phenology method [11], but previous studies showed that the vegetation phenology phase determined by remote sensing data was closely related to those obtained from the ground phenology monitoring points [59]. This study was based on one year's worth of data from the ground phenology monitoring points of P. euphratica, and established a mathematical model of P. euphratica phenology ( Table 1). The phenological features of P. euphratica, including the SOS, EOS, and LOS, were extracted using the dynamic threshold method. To verify the accuracy of the research results and the applicability of the method, based on the GLASS LAI and higher resolution of MODIS NDVI (MOD09Q1 V6) data in 2006, 2010, and 2014, 25 samples were randomly selected to compare and analyze the results of key phenological phase extraction of P. euphratica through the process flow chart steps in Figure 2. As shown in Figure 8, The results revealed that using the GLASS LAI and MODIS NDVI to extract the phenology of P. euphratica yielded higher correlations (SOS: R 2 = 0.861, p < 0.001; EOS: R 2 = 0.824, p < 0.001). Overall, this method could quickly identify and extract the key phenological period of P. euphratica in a large area, and this study's findings are reliable and have scientific reference value for the study of vegetation phenology, especially in extremely arid areas. In addition, phenological research of P. euphratica in the Tarim River basin has been rare, and other researchers [60][61][62][63][64][65][66] believe that air temperature is the main factor affecting plant phenology in arid and semi-arid areas, and the accumulated temperature should therefore be the condition used to judge P. euphratica phenology (Table 3). For example, Zhao et al. [61,66], Liu et al. [64,65], Wu et al. [63], and Zhang et al. [62] monitored the phenology of P. euphratica by taking the date the daily average temperature reached ≥5 • C as the SOS of P. euphratica, and the date the daily average temperature reached ≤5 • C as the EOS in Hexi Corridor. These results showed that the SOS of P. euphratica was early March, which was earlier than the findings of this study. The EOS of P. euphratica was late October, which was consistent with the findings of this study. However, our study showed that the temperature (i.e., AT, ST, and EA) and SD were strongly related to phenology (SOS, EOS) in the upper Tarim River basin. Moreover, vegetation phenology is affected by many factors such as geographical factors and climatic factors. Therefore, using a single accumulated temperature as the criterion for judging the phenology of P. euphratica cannot directly reflect this cycle, and this method is not rigorous. Table 3. Canonical redundancy comparative study on the phenology of P. euphratica [60][61][62][63][64][65][66].

Study Area SOS (Mean) EOS (Mean) Method Data
Lower Tarim

Relationships between Phenology of P. euphratica and Climatic Factors
Previous researchers showed that the SOS of vegetation was advanced and EOS was delayed with climate warming [67][68][69][70]. In this study, however, the SOS and EOS of P. euphratica had an earlier trend, but the change was not significant in the upper Tarim River basin from 2000-2017. Wu et al. [71] reported that there was no significant advance trend in SOS from 1982-2016 in six temperate areas, and Zhang et al. [11] found that there was no advanced trend of spring phenology and the rate of earlier SOS was slow with variations in climate from 2006-2017. Zhang et al.'s [11] study results implied that the trend of EOS was not significant from 2006-2017, which was consistent with the findings of Li et al. [72] and Shen et al. [8]. It is reasonable that the SOS displays a later trend and EOS shows an earlier trend as the microclimate of this region is different. In addition, vegetation phenology is regional, which is closely related to species, geographical factors, and climatic factors. The response of vegetation phenology to climate change is not universally identical in different regions [73]. Therefore, it is necessary to explore the variations in P. euphratica phenology and its response to climate change in extremely arid areas.
In this study, we found that temperature (i.e., AT, ST, and EA) and SD controlled P. euphratica phenology, although the impacts of RH and PR on P. euphratica phenology cannot be overlooked in the upper Tarim River basin. Importantly, Figure 9 revealed that the SOS had negative correlations with PR, SD, and temperature (i.e., AT, ST, and EA), but a positive correlation with RH, while the relationships between the EOS and meteorological factors were the opposite. Air temperature is one of the main factors affecting the phenology of P. euphratica, but precipitation has the least effect on phenology compared to other climatic factors (i.e., AT, EA, ST, SD, RH) in the study area. The temperature increase may lead to earlier SOS, which is attributed to the fact that increasing temperature can promote the chemical reaction rate of vegetation to meet the thermal requirement for budburst and leaf unfolding, and it can accelerate the rate of earlier SOS [8,11,74,75]. positive correlation with RH, while the relationships between the EOS and meteorological factors were the opposite. Air temperature is one of the main factors affecting the phenology of P. euphratica, but precipitation has the least effect on phenology compared to other climatic factors (i.e., AT, EA, ST, SD, RH) in the study area. The temperature increase may lead to earlier SOS, which is attributed to the fact that increasing temperature can promote the chemical reaction rate of vegetation to meet the thermal requirement for budburst and leaf unfolding, and it can accelerate the rate of earlier SOS [8,11,74,75]. Meanwhile, as reported in previous studies (e.g., Che et al. [73]; Dreesen et al. [76]; Anderegg et al. [77]) higher temperatures in spring and autumn will increase surface evapotranspiration and decrease water availability, preventing plant photosynthesis and growth, increasing plant death and chlorophyll degradation, and resulting in earlier EOS. This study found that a temperature increase may lead to earlier SOS and later EOS, which was similar to the findings of Wang et al. [67] and Tao et al. [78], who reported that a temperature decrease may lead to later SOS and earlier EOS among different vegetation types. In addition, the results of this study revealed that the SOS of P. euphratica was strongly related to EOS with a high canonical correlation coefficient of 0.63, indicating that a later SOS may generally be accompanied by a later EOS. This result was similar to the findings of Wu et al. [58], who reported that SOS was positively correlated with EOS from 1997-2014.

Which has a Stronger Influence on Phenology, Runoff, or Precipitation?
Water is a key climatic factor restricting vegetation growth in arid areas [79,80]. The phenological characteristics of P. euphratica in the Tarim River basin are affected by rainfall and runoff. Rainfall directly increases soil moisture and provides direct absorption by the roots of P. euphratica [81,82]. On the other hand, runoff compensates the groundwater so the groundwater depth becomes shallow and is absorbed by the roots of P. euphratica Meanwhile, as reported in previous studies (e.g., Che et al. [73]; Dreesen et al. [76]; Anderegg et al. [77]) higher temperatures in spring and autumn will increase surface evapotranspiration and decrease water availability, preventing plant photosynthesis and growth, increasing plant death and chlorophyll degradation, and resulting in earlier EOS. This study found that a temperature increase may lead to earlier SOS and later EOS, which was similar to the findings of Wang et al. [67] and Tao et al. [78], who reported that a temperature decrease may lead to later SOS and earlier EOS among different vegetation types. In addition, the results of this study revealed that the SOS of P. euphratica was strongly related to EOS with a high canonical correlation coefficient of 0.63, indicating that a later SOS may generally be accompanied by a later EOS. This result was similar to the findings of Wu et al. [58], who reported that SOS was positively correlated with EOS from 1997-2014.

Which has a Stronger Influence on Phenology, Runoff, or Precipitation?
Water is a key climatic factor restricting vegetation growth in arid areas [79,80]. The phenological characteristics of P. euphratica in the Tarim River basin are affected by rainfall and runoff. Rainfall directly increases soil moisture and provides direct absorption by the roots of P. euphratica [81,82]. On the other hand, runoff compensates the groundwater so the groundwater depth becomes shallow and is absorbed by the roots of P. euphratica [34,43]. Based on the annual runoff and precipitation data, this study used the grey correlation method to explore the grey correlation grade between P. euphratica phenology and annual rainfall and runoff (Figure 9). The results demonstrated that the grey correlation grade between runoff and phenology is higher than that between rainfall and phenology. Precipitation is scarce in arid areas, and P. euphratica growth is not sensitive to precipitation events, which is consistent with previous research results [83,84]. However, runoff affects the growth of P. euphratica by replenishing the groundwater table, which is then absorbed by the roots of P. euphratica, hence there is a certain lag effect [85]. Therefore, the grey correlation grade of rainfall and runoff for the phenology of P. euphratica is low, but runoff is still the dominant factor, which is easy to regulate and control. The regulation of runoff and the stable control or reduction in shallow groundwater depth may weaken the impact of global climate change on the phenology of P. euphratica.

Limitations
Through field investigation, the concentrated distribution of P. euphratica has been found to be decreased [86,87]. P. euphratica has gradually declined due to climate changes and disturbance by human activities, and its vegetation coverage is low [34,85]. The retrieval accuracy of the spatio-temporal patterns of P. euphratica phenology using 1 km × 1 km GLASS LAI remote sensing data is not high. Moreover, considering that there is a wide range of mixed pixels in satellite remote sensing data, the spatial resolution of remote sensing data will affect the accuracy of P. euphratica phenological feature inversion. Therefore, to solve the above two problems, this study utilized a grid spatial sampling strategy to select typical sample areas in order to obtain P. euphratica phenological information, thereby reducing the impact of the above two problems to explore the changing trend of P. euphratica phenology in arid areas and its response to climate change. This paper has provided an approach in which the satellite product data are used to reveal the interannual variations of the phenology of the vegetation, plant community, and a specific tree species, of which the spatial aggregation and coverage are low.
Yang et al. [88] found that the preseasons of EOS were defined as 1-6 months before early November, and the preseason temperature is the major control factor of the dormancy onset date. Fu et al. [89] defined May to October as the EOS preseason, and then revealed the influences of climate change and spring phenology on autumn phenology from 1982-2015. Ren et al. [90] found that the phenology was most correlated with climatic factors in the period of about three months to phenology (SOS, EOS). However, the key periods that affect the phenology may vary among different areas [91][92][93]. In this study, the climatic factors could be divided into five temporal periods, which have a certain subjectivity. We should use a suitable method in the future, such as linear correlation analysis, partial correlation analysis, and Pearson correlation analysis, to determine the optimal periods, and make the climatic factors in the period have a strong correlation with phenology.
Since the study area mainly consists of modern alluvial plain and delta plain, the topography is characterized by monotonous flatness and a slight slope of the entire plain, and the range of elevation is small, the impact of topographical factors on P. euphratica phenology was ignored. Due to the limitation of vegetation coverage and aggregation degree of P. euphratica in arid areas, this study needs a longer time series and higher spatial resolution remote sensing data to reveal the phenological characteristics of P. euphratica. Meanwhile, this study established a model based on a single year of ground monitoring data, and the accuracy of the simulation was insufficient, indicating that further improvement is needed in simulation and validation. Thus, these questions await further investigation.

Conclusions
Vegetation phenology is very sensitive to climate change in extremely arid regions, and is thus the best indicator of climate change as well as an essential indicator of ecosystem function. Using the Global Land Surface Satellite leaf area index time series data from 2000 to 2017, the phenological information (i.e., SOS, LOS, and EOS) was obtained with the dynamic threshold method from the reconstructed growth time series curve by using the Savitzky-Golay filtering method. The interannual variation trends of P. euphratica phenology were not significant from 2000 to 2017 because there was no significant trend in annual climatic factors. P. euphratica phenology was remarkably changed with respect to its high spatial heterogeneity. The spring climatic factors had a great influence on the SOS of P. euphratica, while the EOS was mainly affected by the EOS preseason (June-November) meteorological factors. The temperature and SD controlled P. euphratica phenology, although the impacts of RH and PR on P. euphratica phenology cannot be overlooked in the upper Tarim River basin, implying that it is not easy to characterize phenological change with simple climate data. This difficulty is due to the interaction between the SOS in spring and the EOS in autumn, as well as the comprehensive influence of multiple climatic factors on vegetation phenology change, which has not been adequately considered in current vegetation phenology research and ecosystem models. This study also revealed that global warming may lead to SOS advancement and EOS delay, and the increase in SD may lead to earlier SOS and later EOS. Due to the scarcity of rainfall in extremely arid regions, vegetation growth primarily depends on groundwater recharge, so the phenology of P. euphratica has a higher correlation with runoff than rainfall in the upper reaches of the Tarim River. The results provide insights into quantitative assessment of spatio-temporal variations in P. euphratica phenology under climate change and provide a scientific reference for future management of the desert ecosystem.