Quantitative Analysis of Winter Wheat Growth and Yields Responding to Climate Change in Xinjiang, China

: The knowledge of climate change effects on variations of winter wheat yields are crucial for productions. Our objectives were to investigate the relationship between yield-related indices of winter wheat and the related climatic variables (selected using variance inﬂation factors) at the 20 sites of Xinjiang, China over 1981–2017. The background of climate and yield changes was analyzed from temporal and spatial respects. The number of independent climatic variables was selected with the variance inﬂation factor method to remove the multicollinear feature. The Pearson correlation was conducted between the ﬁrst difference values of climatic variables and yield-related indices of winter wheat (namely plant height, growth period duration, 1000-kernel weight, kernel number per ear, biomass and yield) to ﬁnd the key climatic variables that impacted winter wheat growth and yields. The multi-variate linear and nonlinear functions were established step by step using the selected key climatic variables. The best function was determined for each site (signiﬁcant for p < 0.05). From the results, there were general wetter and warmer trends of the climatic variables. Correspondingly, shortened winter wheat phenology and increased growth and yields were observed for most sites. Still, the climatic trends had mixed effects on winter wheat yields. The effects of precipitation, mean air temperature and relative humidity on plant height and growth period duration agreed well. Different sites had different major climatic drivers for winter wheat growth or yields, and the best functions of growth and yields could be linearly or nonlinearly, mostly described by multi-variate functions. The winter wheat growth or yield indices were also found to be closely connected with the soil water content status at the eight sites. The relationship between winter wheat growth or yield and climate provided useful references for forecasting crop production and for projecting the impact of future climate changes.


Introduction
Climate change has resulted in both positive and negative effects on global crop production, with the latter being more common than the former [1]. The theme for World Food Day is "Climate is changing. Food and agriculture must be too." Climate change threatens global crop production and security [2,3]. Sustainably improving crop production is urgently required to meet people's living demand. Nevertheless, the ongoing influences of climate change will increase the risk of meeting this demand for maintaining steady crop production [4]. Understanding the relationship between climatic variables and historical crop production is of importance for assessing the sustainability of the agricultural production.
Wheat is one of the most prevalent grain crops which supports 35-40% people in the world [5]. Wheat is the second food crops following rice in China. The wheat production accounts for approximately 29% of total food production in China [6]. In recent years, climate change has posed huge risks to wheat grain production. With growing population and enhanced living standards following economic development, it is much important to maintain and increase wheat production. Therefore, to quantity the influence of climate change on winter wheat yields is highly necessary for guiding regional farmland production practice and for taking measures to adapt to the global climate change.
Agricultural production is highly dependent on climate, which affects the quantity and quality of harvested crops [7]. Several researchers have reported the effects of climate change on yields of wheat or other crops. For example, Lobell and Field [8] found that growing season temperatures and precipitation explained 30% or more of annual variations in global average wheat yields. Lobell et al. [9] showed that global corn and wheat production decreased by 3.8% and 5.5%, respectively, due to climate change. Ray et al. [10] analyzed the relationship between climate and crop yields and highlighted the roles in temperature, precipitation and their interaction to explain yield variability. Nicholls [11] estimated the impact of climate trends on wheat yields in Australia since 1952-1992 and found that the minimum temperature was the main influencing factor in the related climate variables that contributed to increased wheat yield by 30-50%. Lobell et al. [12] found that the increased wheat yield by 25% in Mexico over 1980-2001, mostly due to climate trends in the northwest United States, especially the decline in nighttime temperature during the growing season. Tao et al. [13] found that China's rice, wheat, corn and soybean production changed by 3.2 × 10 5 , −1.2 × 10 5 , −21.2 × 10 5 and 0.7 × 10 5 t decade − 1 , respectively, due to the climate change. Tao et al. [14] concluded that the trends in mean temperature, precipitation and sunshine duration over 1980-2008 reduced wheat, maize and soybean yields by 1.27%, 1.73% and 0.41%, respectively, while they increased rice yields by 0.56%. Zhang and Huang [15] revealed overall negative effects of the average temperature and diurnal temperature range on wheat and maize yields over 1980-2008 in China. Xiong et al. [16] found that the increasing growing-season temperature had negative responses to wheat, maize and soybean yields in China. Wei et al. [3] estimated the effects of temperature and precipitation on the yields of main grain crops based on provincial data from 1980 to 2008 in China. They found that the contribution rates of temperature change to wheat and yield growth were 1.3% and 0.4%, respectively, and the effect of precipitation change was very small. Yang et al. [17] reported that rising temperatures had a positive impact on the productions of wheat, maize and rice in China. In addition, Bhatt et al. [18] reported that the warming trends of average temperature in the past few decades had a negative impact on yields of rice, corn and wheat in the Koshi River basin of Nepal. Wang et al. [19] found that recent climatic changes have increased wheat yields by 8.5% to 21.2% in four different wheat belts of New South Wales in Australia over 1922-2000. Although studies of climate change effects on wheat yields have been conducted at different scales with much useful and referable results, their applications are confined to varying regions, which require updated relationships between crop yields and climatic variables. This brings difficulties to applying the previous research results to the other places.
To reveal the effects of climate variability on winter wheat yields is a challenging topic since crop yields data are precious and lacking in most places. There are usually two methods to reveal the effects of climate yield relationships: statistical models and crop models [20]. Statistical method is a useful technique to reveal the effects [11,12,21]. Most research focused on the effects of temperature, precipitation and solar radiation rather than wind speed or relative humidity [19,[22][23][24][25][26]. Though temperature, precipitation and solar radiation has significant effects on winter wheat yields, whether the statistical method is a better strategy or not than crop model remains an unanswered question. It is still a matter of importance to study the effects of the statistical method [27][28][29][30]. Statistical models can be obtained quickly and simply. Since climate variabilities often have a non-monotonic effect on winter wheat yields (e.g., increased temperatures will have a positive effect on yields at cool temperatures; however, it will have a negative effect on yields at high temperatures). Therefore, there are linear and nonlinear relationships in statistical methods [7,8]. Nonlinear relationships have a polynomial, cubic spline and interactive. In recent decades, machine learning methods have attracted more and more attention. They have been applied in many fields, such as digital water metering [31], climate extremes and crop yields [29]. Machine learning methods can be used to study the nonlinear and hierarchical relationships between predictors and responses, including artificial neural network-based methods, regression tree-based methods, stochastic-based methods and hybrid methods [32]. These methods generally perform better than traditional linear regression models.
The best relationship between crop growth-yield and climatic variables are not so easy to determine, since the growth-yield data for different sites are not as easy to collect as the climatic variables. Our overall objective was to determine the best quantitative relationship between winter wheat growth-yields and climatic variables at the selected 20 sites of Xinjiang, China. Our specific objectives were (i) to investigate the spatiotemporal variation characteristics in the related climatic variables during the winter wheat growth periods of 1981-2017, (ii) to identify the major climatic drivers of winter wheat and (iii) to establish the best quantitative relationships between yield-related indices of winter wheat and the key climatic variables. This research would provide useful references for further prediction of regional wheat production.

Study Area and Data Sets
Xinjiang Uygur Autonomous Region is located in northwestern China (73 • 40 -96 • 23 E, 34 • 22 -49 • 10 N). Xinjiang is a typical arid and semi-arid region and belongs to a temperate continental dry climate. It is characterized by 'three mountains (i.e., the Kunlun Mountains in the south, the Tianshan Mountains in central Xinjiang, and the Altai Mountains in the northeast) and two basins (i.e., Zhungaer and Tulufan)' [24,33,34]. Winter wheat is one of the three major food crops in Xinjiang.
Daily weather, growth and yield data during the winter wheat growth period from 20 agricultural meteorological stations in Xinjiang during 1981-2017 are collected from the China Meteorological Data Sharing Network with strict quality control. Site distribution and the corresponding geographical information is mapped in Figure 1. The elevations of the selected sites range between 469 and 8464 m. The related climatic variables include precipitation (P re ), mean temperature (T mean ), maximum temperature (T max ), minimum temperature (T min ), wind speed at 2 m (U 2 ), sunshine hour (S un ) and relative humidity (RH). Detailed geographical, regional and long-term mean meteorological information during wheat growth period is presented in Table 1. The P re , T mean , U 2 , S un and RH during the winter wheat growth period range from 23-489 mm, 6.6-12.

The Multicollinearity Test among the Selected Climate Variables
Crop yield is influenced by both climatic and non-climatic factors, which include improved varieties, fertilization development, biocide application, cultivation, irrigation, etc. To remove the non-climatic effects and only evaluate the impacts of climate change on winter wheat growth and yields, the first-difference (∆) of the variable has been effectively used [11,23,35].
Before selecting the best models of yields, the multicollinearity should be analyzed for the selected climate variables to remove the variable that is highly dependent to another. The variance inflation factor (VIF) is useful in assessing the multicollinearity between the studied climatic variables [36]: where R 2 is the coefficient of determination between pairs of climatic variables. The variable which has R 2 > 0.9 or VIF > 10 is removed from the selected climatic variables because of high collinearity. The multicollinearity test is repeated several times until the collinearity relationship is not shown. Finally, the climatic variables which are independent from each other for establishing models of winter wheat growth and yields are kept. In this study, T max is removed from the first run and T min is removed from the second run. After the third run, the variables without any multicollinearity are determined, namely P, T mean, U 2 , S un and RH.

Trend Test
The modified Mann-Kendall test is used here to detect the temporal trends of the climatic and the yield-related variables [37][38][39]. This method accounts for the effects of autocorrelation in the studied time series. The modified Mann-Kendall statistic Z m considers the influences of auto-correlation on original Mann-Kendall statistic (Z) using a correction factor n s 1 [40], described as follows: where: where r j is the jth-order of serial auto-correlation coefficient.
If Z (or Z m ) is positive (negative), the time series have an upward (downward) trend. The null hypothesis is rejected if |Z m | ≥ Z 1−β/2 at a confidence level of β. For β = 0.05, when |Z m | ≥ 1.96, the series have a significant trend [40].

The Best Linear Equations
The Pearson correlation is conducted between pairs of the first-order difference of growth indices (including Y ie , B, G, H, K and W) and climate variables (P re , T mean , U 2 , S un and RH) during winter wheat growth seasons. From this procedure, the sensitive climate variables that influence winter wheat yields are primarily shown. The equation that had the largest Pearson correlation coefficient is selected as a best linear function for winter wheat growth-yield indicators.
The stepwise single and multi-variate linear regression is then conducted to find the best linear functions of winter wheat growth-yield indices, described by the following equation [11]: where ∆Y i is first-order difference of winter wheat growth indices (Y ie , B, G, H, L or W), ∆X i is first-order difference of climate variables (P re , T mean , U 2 , S un or RH) and a i is the regression coefficients (i = 1 to 37 years). The intercept is forced through the origin to avoid trend effects [11]. The stability of the regression functions is estimated by the adjusted coefficient of determination (R adj 2 ). The larger the R adj 2 values, the more stable the model at a significant level [7]. The t-test is used to examine the significance of the functions at the confidence levels of 0.01 or 0.05.

The Selection of the Best Nonlinear Equations
The climatic variables probably have nonlinear and non-monotonic effects on winter wheat growth and yields caused partially by climate variability [7,41]. Considering this nonlinear feature, the multivariate nonlinear regressions are conducted for the winter wheat growth-yield indices and the climatic variables. The nonlinear components include cubic, quadratic and interactive items. If the linear equation includes n climatic variables, 2n + C 2 n of nonlinear variables are added to the basis of the original n linear variables, then a total of 3n + C 2 n variables are used to improve the linear equation. The 1-, 2-, . . . , to 3n + C 2 n -variable regression procedure contains C 1 3n+C2n , C 2 3n+C2n , C 3 3n+C2n , . . . , to C 3n+C2n 3n+C2n (Equation (5)), total ∑C i 3n+C2n equations for each growth-yield indicator will be obtained. The multi-variate regression model [30] is described as: Among all the available nonlinear equations, the equation with the largest R adj 2 is chosen as the best multivariate nonlinear equation of the wheat growth indicators. Its performance is compared with the former obtained best linear equation. Finally, the best equation of the winter wheat growth indices is determined with the largest R adj 2 . The regression analysis and other statistical analysis are performed via R 3.5.3 language software. The spatial distributions of the studied variables are mapped in ArcGIS 10.3 software.

Variations of the Climatic Variables
Five climatic variables of P re , T mean , U 2 , S un and RH passed the multicollinearity test and were selected as the dependent variables for describing winter wheat production. The variations of the selected climatic variables during winter wheat-growing season of 1981-2017 at the studied 20 sites of Xinjiang are plotted in Figure 2. P re ranged from 28.    The trends of the growth-and yield-related parameters over 1981-2017 are mapped for the 20 sites in Figure 5. Different trends were observed for different sites. Few sites showed significant trends. Trends in kernel number per ear increased at all 20 sites. Plant height at 16 and growth period duration at 19 sites had decreasing trends, while 1000-kernel weight, biomass and yield at 18 sites showed increasing trends. General shortened winter wheat phenology and increased growth and yields were observed in Xinjiang.

Pearson Correlation Analysis
The matrix of Pearson correlation coefficients (r) between yield-related indices and climatic variables at each site were obtained (Tables 2 and 3). The results demonstrated that (1) ∆H was positively correlated to ∆P re and ∆RH but negatively correlated to ∆U 2 , ∆T mean and ∆S un at most of the sites. (2) ∆G had negative correlations with ∆T mean but positive correlations with ∆P re or ∆RH at most of the sites. However, the correlations between ∆G and ∆U 2 or ∆S un were not consistent for different sites. The effects of precipitation, mean temperature and relative humidity on plant height and growth period duration agreed well. (3) Correlations between the first differences of 1000-kernel weight (or kernel number per ear, yield and biomass) and climatic variables differed with sites. Values of r changed with sites and variables.

The Best Multivariate Linear Regression Equations
For selecting the best multivariate linear functions, the yield at the site Wusu is taken for example. The first difference of yield (∆Y ie ) was linear correlated with each of the selected five climate variables (i.e., P re , T mean , U 2 , S un and RH) ( Table 4). A linear function of ∆Y ie that performed the best and had the largest R adj 2 among the five was selected (∆Y ie = −8.4∆P re ). The multivariate linear functions of ∆Y ie correlated with 2, 3, 4 and 5 climate variables (22 functions in total) were obtained step by step, of which the equations are given in Table 5. The best 2-, 3-, 4-, 5-variate linear functions that had the largest R adj 2 values in each group were selected (bold in Table 5). Through further comparing four better equations, one 2-variate equation was selected (∆Y ie = −12.2∆P re − 623∆T mean ) as the best linear function of ∆Y ie . However, this best relationship performed not very good with R adj 2 value of 0.29. Nonlinear relationship should be tried to see if the performance could be improved.
Four variables 28∆S un 0.28 * * and ** mean significance at confidence levels of 0.05 and 0.01, respectively.

The Best Multivariate Linear Regression Equations
The performance of the single and multi-variate nonlinear regression is shown in Table 6, also taking the site Wusu as an example. It was shown that both 2-and 3-order polynomial and interactive regression functions had improved performance (R adj 2 > 0.34) compared to all of the linear functions. Finally, the equation ∆Y ie = −12.1∆P re − 604∆T mean − 4.04 ∆P re × ∆T mean was considered the best function for describing the relationship of ∆Y ie with climate variables, considering its largest and relatively simple format. A similar procedure was conducted for finding the best functions of the other growthand yield-related variables at all sites. The site-specific best functions at all sites are presented in Table 7. Unfortunately, not all of the yield-related variables could be described as functions of climatic variables, and the best functions were obtained for only less than half the studied sites. In addition, the largest R adj 2 values for plant height, growth period duration, 1000-kernel weight, kernel number per ear, biomass and yield were 0.41, 0.27, 0.27, 0.50, 0.41 and 0.40, respectively. The performance of these relationships showed that climate change only explained less than 50% of the variability in winter wheat growth. The obtained equations could be used in prediction of winter wheat yields but with not too high accuracy. Still, these equations are useful, since yield prediction is important.

Roles of Soil Water Storage (Content) on Winter Wheat Growth and Yields
The yield and growth of winter wheat should also be related to soil water content, since irrigation is also a key technique when soil water is deficient. However, since soil moisture data were not as complete as weather data, the best functions of the growth and yield indices did not include the effects of soil water content. Still, its effects on winter wheat growth and yields could be discussed using some limited data. The volumetric soil water content (θ v ) at a depth of 0-10, 10-20, 20-30, 30-40 and 40-50 cm of eight sites during 1981-2012 were collected for assessing its effects on winter wheat growth and yields. The mean θ v during winter wheat growth period at different depths were correlated to the winter wheat yield-related indices for the eight sites and the Pearson correlation coefficients are presented in Table 8. It was shown that θ v had complicated varying influences on winter wheat growth and yields. This result agreed well with Li et al. [30]; their results showed that cotton stalk weight was positively correlated with soil moisture at depths of 0-10 to 40-50 cm at 11, 9, 10, 9 and 10 sites, while negative correlations were also detected between soil moisture at different depths and cotton yield indices. Previous studies have emphasized the importance of soil moisture on crop yields [42,43]. The plant height, growth period duration, 1000-kernel weight and kernel number per ear could be positively or negatively related to θ v at certain depth. This occurred because (i) θ v values were averaged for the entire growth period, (ii) the depths played different roles and (iii) the other factors including wheat variety and soil properties also affected θ v . Overall, θ v at different depths explained less than 61% of variability in the winter wheat growth and yields.

Discussion
In this research, the winter wheat growth period duration reduced during the research period (Figure 5b). Previous studies agree with our results very well [44][45][46]. The trends of winter wheat's 1000-kernel weight, kernel number per ear, biomass and yield in Xinjiang increased at most sites. Adaptions such as delaying sowing date and changing cultivars would be considered. Zhou et al. [47] reported that the genetic improvement in grain yield was mainly contributed to increased grain weight per spike and reduced plant height, their results also supported our results in Figure 5a,c.
The impacts of climate change on crop yields can vary generally depending on the assessment technique (e.g., crop models or statistical models), crop varieties, climate data and crop model initial settings; therefore, the impacts of climate change on winter wheat growth and yields can be different [19,[48][49][50]. Research using crop models without changing crop variety and management measures have shown that climate change had a negative impact on crop yields over the past three decades because of reduced growth duration [21,[51][52][53]. In this research, the site-level data showed that first-difference climate variables had both negative and positive effects on the first-difference of yield-related indices (Tables 2 and 3). One possible explanation is that in different growth periods, winter wheat had changing temperature thresholds and optimum temperatures [27,54]. The second reason is that the studied sites had different topographies. The third is that wheat growth processes (e.g., photosynthesis) occurred exclusively during the day, and the temperature during the day may not always increase during night [55].
The largest R adj 2 values for the best functions of plant height, growth period duration, 1000-kernel weight, kernel number per ear, biomass and yield were 0.41, 0.27, 0.27, 0.50, 0.41 and 0.40, respectively, all <0.50. This was due to the complexity (non-linearity) of the climate system and the complexity of crop growth and yield processes. Our research showed that the significance of climate change only explained less than 50% of the variability in winter wheat growth. Previous studies obtained different ranges of R adj 2 . For example, Lobell et al. [12] obtained an R adj 2 value of 0.56 when using the multivariate linear regression functions between first differences of wheat yield and climatic settings in Mexico. The highest R 2 value of regression model between wheat yield and climate over 1988-2002 was 0.47 in the southwestern plains of Australia [19]. Lobell and Field [8] established a regression model between yield and climate first-differences over 1988-2002 with an R 2 value of 0.41 for global scale of wheat. Liu et al. [5] and Lobell and Burke [21] obtained a large R 2 value of 0.37 in the key wheat-growing region of China during 1960-2009. Since the advances in genotype and agricultural management measures such as fertilization and irrigation also affect and contribute to winter wheat yields, it was reasonable that climate change could not explain more than 56% of the winter wheat growth and yield.
Further research is, therefore, needed to better understand the roles of other factors, such as genotype and agronomic management, on winter wheat production.

Conclusions
Five climate variables (namely P re , T mean, U 2 , S un and RH) without any multicollinearity were selected for further finding the best functions of winter wheat yield-related indices at the 20 stations in Xinjiang, China. Most sites were wetter, warmer, received less sunshine and experienced lower wind speed and humidity. Correspondingly, winter wheat plant height and growth period duration decreased, and the kernel number per ear, biomass and yield increased at most of the sites.
The Pearson correlation coefficients between yield-related indices and climatic variables ranged between −0.66 and 0.57. In general, the yield-related indices had strong or weak positive correlations with P and RH, but had negative correlations with T mean , U 2 and S un . From the Pearson correlations, some key climatic variables that affected winter wheat crop yields were selected and used for further multi-variate linear and nonlinear regressions. The best functions of the wheat growth-and yield-related indices which could be linear or nonlinear were finally obtained step by step. The adjusted coefficient of deter-mination varied between 0.09 to 0.50, which implied that climate change only explained at largest 50% of the variability in winter wheat growth and yields. Some other factors, such as genotype and agronomic management should be considered in future study. Appendix A