Phenological Changes of Soybean in Response to Climate Conditions in Frigid Region in China over the Past Decades

Plant phenology becoming a focus of current research worldwide is a sensitive indicator of global climate change. To understand observed soybean phenology and explore its climatic determinants in frigid region (Northeast China and northeast in Inner Mongolia), we studied the phenological changes of soybean [Glycine max (L.) Merr.] for the frigid region during 1981–2017, then analyzed the contribution of major causal climate factors to phenology based on multiple stepwise regression. Altogether, the average temperature from sowing to maturity (WGP) was significant increasing, accumulated precipitation and sunshine hours were decreasing. More than 50% of observations showed delays in sowing, emergence and maturity stage and short durations of sowing to flowering (VGP), flowering to maturity (RGP) and sowing to maturity (WGP). The late sowing was getting the following phenological timing backward, but the flowering and maturity delaying trends were much less than that of sowing timing due to the warming accelerated growth of soybean. Detailed analysis indicated mean temperature and accumulated precipitation of the 1–3 months immediately preceding the mean emergence, flowering and maturity dates influenced the phenological timing in higher latitude areas (HLJ and FL), while in JL and LN, accumulated precipitation and sunshine hours(replacing mean temperature) were the climatic determinants. These results brought light the importance of research and policy to support strategies for adaptation to local condition under the climate change.


Introduction
With global warming in the last half century (Stocker et al., 2014), it is an indisputable fact in China, especially in the frigid region (Northeast provinces and northeast of Inner Mongolia in China) located in high latitude and cold areas (Du et al., 2018;Feng et al., 2018;Zhai et al., 2017). The increasing number of studies draw extensive attention of ecological consequence of global climate change. The climate warming accelerates the growth and development of crops, leading to earlier/later phenology or shorter/longer growing season (Bertin, 2008;Gordo & Sanz, 2010;Szabó et al., 2016;Vitasse et al., 2018). The shifts of plant phenology and growing season, being part of an ecosystem, reflect vegetation production and energy exchange in soil-vegetation-atmosphere system. That in turn influencing climate system is contributing to climate change (Richardson et al., 2013).
The phenology is one of the crucial indicator reflecting changes in climate and environment Qin et al., 2018;Wang et al., 2016;Xiao et al., 2015). Moreover, the climate has an immediate influence on agricultural production (Novikova et al., 2020;Richardson et al., 2013). Recently, interest in phenology shift of natural vegetation has increased. Several studies have reported phenological changes in response to climate change (Gordo & Sanz, 2010;Richardson et al., 2013;Szabó et al., 2016;Vitasse et al., 2018;Xu et al., 2020aXu et al., , 2020bYue et al., 2015). Studies have also used crop models driven by climate factors to simulate phenological events (Larijani, 2011;Shuai Zhang & Tao, 2013;Yue et al., 2015;Huang et al., 2017). The studies provide strong evidence that impact of climate change on phenological phases in vegetation, forest and crops. In recent 1 3 decades, spatiotemporal changes of crop phenology and its shifts in response to climate change have documented (Bai & Xiao, 2020;Fatima, 2021;Lizaso et al., 2018;Tao et al., 2012;Xiao et al., 2015Xiao et al., , 2016. Many studies have shown that an earlier or later onset of phenology due to climate change for rice (Bai & Xiao, 2020), maize (Wang et al., 2016;Xiao et al., 2016), soybean (Fatima et al., 2021), wheat (Xiao et al., 2015) and barley (Banihashemi et al., 2020). Thus, it is necessary to study the impact of climate change on soybean phenology to recognize the adaptation strategies in China.
The frigid region, the origin and top producers of soybean, is a base of cash crop in China (Xu et al., 2020a(Xu et al., , 2020b. In 2010, soybean-planting area in this region reached 4 million hectares about 48% of total areas in China. The yield of soybean achieved 7 million tons that accounted for nearly 47% total output (Hu et al., 2017). A certain degree of increasing temperature contributed larger yield in Heilongjiang Provinces (Jiang, 2011). The increasing temperature accelerated soybean growth progress, leading to much shorter whole growth period (Qiu et al., 2018). However, climate warming and drying can cause reduction of soybean yield (Hou et al., 2011;Qu et al., 2014). Bhatia and other researchers (Bhatia et al., 2008;Kessler et al., 2020;Mall et al., 2004) found that delaying planting could mitigate effect of high temperature stress on soybean production, increasing climatic potential productivity for soybean. Nevertheless, a great number of studies focus on effect of one climatic factor, the increasing temperature, on phenology or yield of soybean. In fact, the growth and development of crops are not only influenced by temperature, but also by precipitation, insolation and other climatic factors (Cai & Fu, 2018;Goldblum, 2009;Li et al., 2016). It could be increasing uncertainty if any factors ignored. Therefore, it is essential to understand the phenological variation in the frigid region in China and explore the driving climatic factors to help the government conduct agricultural activity and reduce the negative effects of climate change. Manual recording is the most intuitive and accurate method to obtain vegetation phenology, which provides important direct data for the study of phenological characteristics and its relationship with environmental factors. Since the 1980s, a lot of phenological data of soybean has been accumulated by agrometeorological observations, which lay a foundation for studying phenological changes and the effects of environmental changes on phenology of soybean. In this study, we investigated phenological characteristics and drivers of phenological shifts based on a set of long-term phenological records of soybean and climatic data of frigid region during 1981-2017. The study aims to (i) quantify characteristics of phenological changes of soybean in the frigid region, and (ii) to determine which climatic factors exert the greatest effect on the phenological shifts of soybean in recent decades.

Description of the Frigid Region in China
The studied frigid region located in northeast China, high latitude and cold weather (Fig. 1). The region included Heilongjiang Province (HLJ), Jilin Province (JL), Liaoning Province (LN) and 4 Leagues in northeast of Inner Mongolia (FL), including Hulun Buir, Hinggan League, Tongliao and Chifeng. The frigid region, located mainly North of 38.6°-53.7°N and East of 115.5°-135.1°E, had a vast geographic area covered with black and fertile soil (Kang et al., 2016;Zhuo et al., 2019;Shaoliang Zhang et al., 2020). Monsoon climate, temperate characteristics was long cold in winter and warm rainy in summer. Sowing time of soybean generally was in April to May and maturity time was in September to October.

Data Collection
To search for potential shifts in phenological events of soybean (sowing, emergence, flowering and maturity), we viewed that agricultural data from 33 agro-meteorological stations (Table 1) in the frigid region in China. We observed soybean phenology when 50% of plants reached the specific stage. For statistical analysis, the phenological event date was recorded as the days of year (DOY). The 1st January recorded DOY 1, 2nd January as DOY 2 and so forth, the 31st December as DOY 365. With the regard to climatic data, daily mean temperature, precipitation and sunshine hours were also gathered. All climatic records and phenological data of soybean based on the phenological observations were taken from the China Meteorological Administration (CMA). In this study, we defined vegetative growing period (VGP) of soybean to be the period from sowing stage to flowering stage, and reproductive growing period (RGP) to be the period from flowering stage to maturity stage. In addition, the period from sowing stage to maturity stage was explained as a whole growth period (WGP) of soybean.

Data Analysis
To identify the combination of the climatic factors and the effective periods, which actually affected phenology of the soybean, we used the method of Gordo and Szabó 1 3 (Gordo & Sanz, 2010;Szabó et al., 2016). Consequently, a series of monthly, bimonthly and trimonthly means (temperature) and sums (precipitation and sunshine hours) variables were calculated, which could influence the phenological events of soybean. 5 monthly mean temperature (T 1 ) variables relative to whole growth period of soybean from the May (T 15 ) to September (T 19 ) were created. The numbers 5 to 9 represented the months of May to September. Likewise, we defined the 4 bimonthly (T 2 ) variables, with T 26 meaning the mean temperature of the May to June, T 27 being the mean temperature of the June to July and the T 29 being the mean temperature of the August to September. 3 trimonthly temperature (T 3 ) variables were created by the similar way, with T 37 (mean temperature of current year's May to July), T 38 (June to August) and T 39 (July to September). In the same way as that of the T variables created, we defined precipitation and sunshine hours sum variables: 10 monthly (P 15 -P 19 , S 15 -S 19 ), 8 bimonthly (P 26 -P 29 , S 26 -S 29 ) and 6 trimonthly (P 37 -P 39 , S 37 -S 39 ) precipitation and sunshine hours sum variables. For example, the bimonthly variable (P 26 ) was the sum of precipitation from May to June, P 37 with the total of precipitation from May to July. A linear equation was defined to describe quantitatively variation tendency of climate and phenology of soybean for the time series.
where Y is the observed climatic factors or phenological date; b is the linear regression slope. It described the magnitude and significance of the trends in climatic variables and phenological date over the time and the sign of b reflects direction of change. The positive /negative sign indicates increase /decrease or delay/advance with time; x is the year. In this study, statistical significance is determined using the two-tailed t test analysis.
As a next step, based on the climatic variables computed before, we seek for statistical relationship between the phenological data and climatic variables. After this relationship step, we created separate linear egression models for each phenological stage (emergence, flowering and maturity) in the frigid region (HLJ, JL, LN and FL). Due to the intercorrelation of the predictors, we used stepwise selection to sequentially add variables from the set of the potential predictors until getting an optimum regression model. All data preparation was calculated in EXCEL and statistical analysis steps were performed in SPSS.

Climate Conditions
We found a significant increase of mean temperature (0.5-0.9 °C decade −1 ) during VGP for the frigid region in the study period. During the RGP, there were 2 parts (HLJ, FL) with the significant increase of mean temperature. During the RGP, all the trends of mean temperature were positive, even the nonsignificant ones. The strongest increase of mean temperature (0.8 °C decade −1 ) during the WGP was found in FL. Altogether, the mean temperature showed an increase for the frigid region during the WGP (Table 2).
Unlike the temperature, precipitation increased during the VGP in HLJ and LN (4.8 mm decade −1 and 5.5 mm decade −1 ), while decreased in JL and FL. For the frigid region, the precipitation reduced during RGP. In addition, the significant decrease of precipitation (− 18.6 mm decade −1 ) was during the RGP in HLJ. The precipitation during the WGP decreased in the study area with nonsignificant ones. The strongest decreased of precipitation was in FL during WPG.
For sunshine hours, we found significant decrease (− 21.0 to − 64.8 h decade −1 ) during VGP for the frigid region in the study period. Except for LN (− 8.9 h decade −1 ), the Table 2 Trends of climate factors in different growth period of soybean in the frigid region Trends are expressed as linear coefficients (positive values mean increase and negative values mean decrease), significance level symbols: ***0 < P < 0.001; **0.001 < P < 0.01; *0.01 < P < 0.05. Relationships found to be significant above α = 0.05 are highlighted in bold 1 3 sunshine hours increased during RGP in HLJ, JL and FL. The sunshine hours exhibited decrease during WGP in the frigid region. There were 2 parts (JL, LN) with a significant decrease of total sunshine hours. Given all that, we knew mean temperature increased significantly, precipitation and sunshine hours decreased during WGP in the frigid region in China. As climatic conditions during the WGP of soybean became warmer and drier, so did the phenology of soybean with the combination of climatic factors (temperature, precipitation and sunshine hours).

Spatiality of Phenological Phases
Phenological events of soybean showed an obvious spatial pattern for the study period ( Fig. 2A-D The periods of growth stages for the frigid region was plotted in Fig. 2E-H. The time duration in days from sowing to emergence stage (S-E) was 10-22 days. The VGP was 53-75 days. In addition, the RGP was 60-84 days. The WGP was 122-150 days. There was clear spatial patterns in VGP and WGP. They remarkably decreased with the increasing latitude. Whereas RGP was not relevant to latitude. However, mean length of RGP, RGP in HLJ and JL was longer 5 and 12 days than it was in LN and FL, respectively.

Temporality of Phenological Phases
The postponed tendency of soybean phenology predominated in different phases in the frigid region. 63.2%, 68.4% and 60.5% of the 33 sites had delayed trends respectively in sowing, emergence and maturity stage (Fig. 3). 50.0% of sites updated in flowering date. There was a significant (P < 0.05) delay up to 10.9 days decade −1 in sowing stage of soybean in 11 sites and a significant advance up to 8.2 days decade −1 in 5 sites. Emergence stage of soybean delayed significantly up to 6.1 days decade −1 in 10 sites and advanced 7.4 days decade −1 in 5 sites. Unlike the early two stages of soybean, the number of sites with early trend was same as the number of sites with postponed trend in flowering stage. Flowering stage of 4 sites advanced up to 5.0 days decade −1 and 6 sites delayed up to 6.4 days decade −1 significantly. Similar to sowing and emergence stage, delayed maturity was noted in the majority of sites. Moreover, there were 6 sites with delayed trend and 6 sites with early trend significantly.
Phenological phases of soybean corresponded with the duration of different soybean growth stages (Fig. 3E-H). 81.6%, 78.9%, 55.3% and 81.6% of sites had short a duration for time of S-E, VGP, RGP and WGP. From S-E, there was significant (P < 0.05) for 16 sites with shortened time and 1 sites with prolonged time. The time of VGP shortened significantly up to 9.9 days decade −1 for 14 stations. Unlike VGP, the time of RGP prolonged significantly by up to 4.9 days decade −1 for 5 sites and shortened notably on average 3.0 days decade −1 for 6 sites. The time of WGP shortened on average 3.8 days decade −1 in 26 sites, in addition, significant for 17 stations.
For average changes in parts of the region, we found that phenological stages (sowing, emergence, flowering and maturity) of soybean advanced in HLJ and delayed in JL and LN. There was a delayed trend of emergence and advanced trends of flowering and maturity in FL. Except the maturity trend in LN, the delayed trends gradually decreased with growth of soybean. E.g. in JL, the trend of sowing stage was 4.6 days decade −1 , emergence trend was 2.5 days decade −1 , then the trend reduced to 0.17 days decade −1 in the maturity. Although phenological variety in frigid region, the shortened duration of soybean growth stages showed up consistently except the RGP in HLJ (Fig. 4).

Phenology and Climate
11 significant relationships were discovered between the date of the emergence and the meteorological variables studied in the four different parts. 21 meteorological variables had a significant effect on flowering events and 33 meteorological variables had a bearing on the maturity events of soybean (Table 3). There were 36 relationships linked with temperature, whereas precipitation and sunshine hours were found 10 and 19 relationships, respectively. For the phenological events, the most predictors were meteorological variables of 1 to 3 months preceding the dates of the emergence, flowering or maturity events. We also found 1 significant temperature predictor of the 3 months (T 36 ) that delayed influence of the past emergence event in HLJ. Compared with the preceding meteorological variables, the delayed ones, which connected with phenological events, were much weaker and less. The coefficients for most temperature and sunshine hour variables were generally negative, showing that warm year or many sunshine hours generally induced earlier phenological events of soybean. On the other hand, coefficients for precipitation variables were usually positive, that meant more precipitation generally induced later flowering and maturity events. There were 33(4 for emergence, 13 1 3 for flowering, 16 for maturity) out of the 65 climatic predictors studied in HLJ that were included in any of significant relationships identified among regions and phenology. In addition, 20 significant relationships (3 for emergence, 3 for flowering and 14 for maturity) were found in FL. However, there were 5 and 7 significant relationships were found in JL and LN, respectively. Coefficients for meteorological variables varied among phenological events. The variables influencing emergence was primarily temperature variables studied. For flowering and maturity events, more precipitation and sunshine variables also were involving in growth of soybean. So according to the analysis of single predictor and phenological event, it was difficult to estimate the influence of climatic change on soybean growth. We took the climatic predictors, temperature variables, precipitation variables and sunshine hour variables that were studied, as the potential predictor. Then the potential predictors were analyzed with multivariate stepwise regression to explore influence of the climatic variables on soybean phenology. In the end, we created climate models predicting phenology (emergence, flowering and maturity) of soybean (Table 4).
Except mean temperature of April to May (T 25 ), sunshine hours (S 14 ) in April was also the significant climatic variable influencing the emergence of soybean in HLJ. The two climatic variables accounted for about 1/2 influence on emergence date, of which the leading role was T 25 with significance of 0.0001. The date of emergence advanced with rising temperature of April to May. In JL and LN, we also found that temperature and sunshine hours (T 14 and S 25 ) affected emergence of soybean, respectively. However, precipitation in April (P 14 ) played a key role in emergence in FL. The leading climatic variable, with significance of 0.06, accounted for 21percent influence of the event in this region (Table 4).
Flowering stages were affected by the temperature of May to June (T 26 ) and May to July (T 37 ) in HLJ and FL, separately. Two temperature variables took above 25 percent influence on flowering event. Unlike in HLJ and FL, precipitation and sunshine hours influenced the flowering event in JL and LN. Sunshine hours of August (S 18 ) and precipitation of April (P 14 ) were the dominant climatic factors, which clarified 27 percent of flowering event in JL. For LN, besides of temperature (T 17 ) and sunshine hours (S 36 ), precipitation (P 26 ) of May to June took part in the regression. Except JL, temperature factors dominated the influence of flowering event in HLJ, LN and FL, meaning higher temperature could advance the flowering event of soybean.
Mean temperature of June to July (T 27 ), precipitation of September (P 19 ) and April to May (P 25 ) were the climatic factors influencing the maturity event of soybean in HLJ. All factors accounted for 48 percent maturity of the changes. The T 27 was the principal factor, of which partial correlation coefficient passed the statistical test. In JL, we found that sunshine hours of September (S 19 ) and precipitation of August (P 18 ), collectively explaining 24 percent of maturity variation, dominated the maturity event. Sunshine hours of April (S 14 ) and mean temperature of July to August (T 28 ) affected maturity event in LN, the two climatic factors jointly accounted for 29 percent of the maturity period as well. Unlike the other three parts, there were three climatic factors (S 29 , T 27 and P 38 ) influencing the maturity event in FL. However, the leading roles were sunshine hours of August to September and mean temperature of June to July, with significance of 0.01. The three climatic factors altogether accounted for 47 percent of maturity change. Except in JL, warming within the 2 to 3 months before the maturity stage accelerated soybean ripening.

Discussion
Shifts in soybean phenological events in the frigid region were particularly poorly recorded, with a few exceptions coming from the works of Jiang et al., (2011), Liu et al., (2007 and Qiu et al., (2018). Based on datasets for the frigid region, the high latitude and cold area of China, so far this paper presents the most thorough of the phenology and their climatic drivers.
Although soil conditions in parts of the frigid region are suitable for soybean (Kang et al., 2016;Zhuo et al., 2019), phenology of soybean significantly vary with climatic conditions. Sowing and emergence events decreased with increasing latitude from south to north while flowering and maturity did not change significantly with the latitude. Moreover, for the duration of the growth stages, the time from sowing to flowering and sowing to maturity decreased with the increasing latitude, while the flowering to maturity did not show increase or decrease with increasing latitude.
Throughout the HLJ, the majority of phenology during 1980 to 2008 revealed advancing trends, especially abnormal advance in the west of HLJ (Jiang et al., 2011). Similar insights into the delayed phenology in JL by Qiu et al. (2018), we also discovered that delaying trends in 8 of 9 cases for emergence period, 7 of 9 cases for flowering and 8 of 9 cases in LN for maturity period. The delaying phenology by the number of more cases in LN was similar with the results of Liu et al., (2007). We also detected that shifts of phenology in FL were not similar to the tendencies in other parts. Moreover, shorter duration of WGP was consistent in most of studied cases. This is consistent with the previous reports about crops, for example, wheat (Tao et al., 2012), maize (Lizaso et al., 2018), rice (Bai & Xiao, 2020) and typical temperate grassland (Xu et al., 2020a(Xu et al., , 2020b. Nonetheless, differences in phenological changes exist among regions and studies due to differences in methodology, study periods and areas. There was no doubt that warming increased the likelihood of early sowing, while increased the probability of drought. Meanwhile, sowing depended not only on climatic conditions, but also on soil moisture, soil fertility, or sowing management practices. Therefore, sowing with water (Szabó et al., 2016)or delaying sowing date (Kessler et al., 2020;Mall et al., 2004), main cultivation modes, could partially counteract the negative effects of warming climate in the semiarid region (e.g. west of LN, west of JL and FL) Zhen et al., 2020). Few analyses about climate change that would result on sowing event were conducted in this paper. As noted earlier, delayed trends in flowering and maturity stages behind the sowing decreased in JL and LN.
Similarly to the results exploring the relationship between phenology and climatic variables, we detected variables in relation to temperature (P < 0.01) were the most influential factors of the phenological events. Moreover, temperature and sunshine hours were negatively correlated with soybean phenology. Warming climate with more insolation might reduce water availability (Dai et al., 2004). It is difficult for crops to effectively utilize the increased heat and carbon dioxide resources to improve the photosynthetic production, with negative impacts on plant growth (Du et al., 2018). In addition, this may accelerate chlorophyll degradation to increase plant maturity (Xiao et al., 2016;Zhao et al., 2015). The correlation between precipitation and soybean phenology was positive over the study region. This is in line with studies of soybean (Goldblum, 2009) and maize (Wang et al., 2016). In most growth stages, the effective climatic variables were 1-3 months preceding or containing the mean month of phenological events. This result was approvingly supported by other studies as well (Moradi et al., 2013;Szabó et al., 2016). Beyond question, temperature is the most important among these climatic factors influencing the phenological events. Even though spring phenology are known to be advanced to rise temperature in the Northeast of China (Huai et al., 2020;Qin et al., 2018), we did not detect phenology of soybean advanced consistently.
Climate warming in frigid region in China has become an indisputable fact in future (Chu et al., 2017;Hou et al., 2019). The rise is about 2-3 °C in mean annual temperature, the accumulated temperature ≥ 10 °C has increased, the start date ≥ 10 °C has advanced by the end of this century. It seems that the increase of heat resources prolongs the growing season of crops. However, due to the warming, the phenology shifts and the growing period gets shorter, which might accelerate the rate of grain filling and reduces the time of photosynthesis (Bai & Xiao, 2020;Lizaso et al., 2018;Tao et al., 2012). The yield might increase but the grain quality would decrease (Guo, 2015;Zhao et al., 2020). Precipitation in this region remains essentially unchanged with a slight increase. Nevertheless, the imbalance with water resources may have impacts on agricultural productivity. Furthermore, climate warming might lead to extreme weather and climate events, increasing the frequency of agricultural meteorological disasters such as drought, flood, high temperature and low temperature. The risk of meteorological disaster for agriculture production might increase (Bucchignani et al., 2017;Zhao & Dai, 2017). Therefore, when adjusting the crop distribution in the frigid region, we should not only fully consider the advantages of the increase in thermal resource, but also the precipitation and insolation resources. In order to ensure the safety of grain production, we should consider the allocation of insolation, heat and

Conclusions
The study investigated the observed climate-induced changes in soybean phenology and discussed the climatic factors influencing phenology for soybean in the frigid region of China. On the spatial distribution, phenology in sowing and emergence showed a gradual delay from south to north. The VGP and WGP were shorter from south to north as well.  Temperature  T 14  ----T 15 - Phenology in sowing, emergence and maturity was delayed at most representative sites in frigid region. The number of sites with delayed flowering event was equal to the number of sites with early flowering event. Due to climate change, S-E, VGP, RGP and WGP were shorter in most sites. Soybean phenology was negatively correlated with temperature and sunshine hours, while positively correlated with precipitation. Furthermore, phenology and causes of phenology of soybean in the frigid region were also studied with multiple factors combination statistical method. Climate variables of the 1-3 months preceding the emergence, flowering and maturity dates influenced the phenological events in the frigid region. The connection between soybean phenology and climate factors offered some insight into the future adaptation of climate change. Moreover, this research can help conduct agricultural production in the frigid region of China.

Conflicts of interest
The authors did not have any conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.