Climate-mediated dynamics of the northern limit of paddy rice in China

Paddy rice agriculture plays an important role in food security and has a considerable influence on natural systems. In the context of climate change, understanding the nature and drivers of shifts in the northern limit of paddy rice (NLPR) is crucial for adaptation strategies and food security. However, quantitative studies on the effect of climate change on paddy rice distribution shifts have not been well performed. Here, we mapped the NLPR in China using Landsat imagery from 1984 to 2013, analyzed the latitudinal and elevational dynamics of the NLPR using Fishnet analysis, and explored the factors driving the changes in rice area across the NLPR regions using a linear regression model. Our results show that between 1984 and 2013, the NLPR shifted 24.93 km northward (the greatest movement was 88.01 km occurring at approximately 133° E) and elevational limits increased by 39.15 m (the greatest movement was 117.08 m occurring at approximately 129° E). While socioeconomic factors (e.g. benefits, policies, irrigation, and mulch) played significant roles in rice area changes, the changes in rice area across the NLPR regions had the strongest positive association with the increase in the previous temperature, indicating that rice cultivation in the NLPR regions has moved to higher latitudes over the 30 year study period to adapt to climate change. Our study highlighted that quantifying the interactions between climate change and crop production systems can facilitate a better understanding of the human responses to changes in the growing conditions in the face of climate change and ensuring regional and global food security.


Introduction
Paddy rice agriculture accounts for 11% of global arable land (Qian et al 2020) and is essential to global food security, contributing 19% of the world's daily energy supply and feeding 50% of its population (Elert 2014, Xia et al 2020. Rice accounts for over 40% of the total calorie intake in China and rice production is a source of employment and income for many farmers (You 2012). Paddy rice cultivation is also associated with a series of environmental impacts, including water resource depletion (Luo 2010) and greenhouse gas (GHG) emissions (Jiang et al 2017, Kritee et al 2018. For example, most paddy rice is grown in continuously or intermittently flooded fields and consumes a quarter to a third of the world's freshwater (Bouman et al 2007, Kuenzer and Knauer 2012. Furthermore, paddy rice planting accounts for approximately 11% of anthropogenic methane (CH 4 ) emissions under continuously flooded conditions (Carlson et al 2016, Jiang et al 2019 and enhances nitrous oxide (N 2 O) emissions under intermittently flooded conditions (Kritee et al 2018); both CH 4 and N 2 O are major GHGs contributing to climate change (Tan et al 2019). Therefore, understanding the spatio-temporal dynamics of rice cultivation is important for both global food security and environmental performance.
Paddy rice mainly grows in the tropical and subtropical regions of India and China, particularly in southern China (Zhang et al 2017). Over the past four decades, paddy rice production has dramatically expanded northeastward in China, increasing cultivation in mid to high latitude regions (Liu et al 2013, Fan et al 2018, Pan et al 2019. Climate change is expected to exacerbate the conversion of land to rice fields at higher latitudes. Paddy rice is a temperaturesensitive crop and climate change has a fundamental impact on rice planting in mid-high latitude regions by yielding an earlier transplanting date, as well as prolonging rice cultivation and decreasing cold damage (Ohta andKimura 2007, Seck et al 2012). For example, increasing temperature will advance the transplanting date of rice by 14-20 days and prolong safe rice cultivation by 25-30 days, further decreasing the risk of cold summer damage, which will provide more safe area for rice planting in northern Japan (Ohta and Kimura 2007). Research shows that the global average temperature has risen by approximately 0.8 • C since 1850 (Wheeler andBraun 2013, Sanchez et al 2014). However, temperature increases are unevenly distributed. In particular, mid to high latitude regions have experienced larger temperature increases (Lobell et al 2011), making them more suitable for rice cultivation . However, climate change impacts on locations of paddy rice and how, in turn, the rice planting system adapts to climate change remains highly uncertain.
The northern limit of vegetation, an important and widely used indicator by geographers and ecologists, is the northern frontier of the habitat at which vegetation grows, which indicates that beyond this northern limit, vegetation does not grow as the conditions (e.g. climate, topography, and other socioeconomic factors) are not suitable (Harsch et al 2009, Körner 2012). The indicator is valuable for characterizing environmental change (especially climate change) and for understanding the relationship between vegetation and the surrounding environment. In this way, the northern limit of paddy rice (NLPR) can be used to understand the interaction between climate and socioeconomic change and rice production and how rice production systems adapted to climate change. Therefore, an understanding of the NLPR is critical for managing both climate change and human adaptation.
Rice cultivation in Heilongjiang Province, Northeast China is the northernmost regions of paddy rice cultivation in China. Over the past five decades, Northeast China experienced a remarkable increase in surface annual mean temperature with a trend of about 0.38 • C/decade (Piao et al 2010, Liu et al 2012. Such climate change allows rice cultivation to expand to the north (Hu et al 2019). Furthermore, socioeconomic factors are thought to have impacts on the NLPR by altering rice planting behaviors. For example, increased agricultural inputs and higher profits (the net benefits to rice are one and a half times that of upland crops) (Pan et al 2019), improved technologies (e.g. greenhouse nurseries and irrigation) , and the implementation of relevant policies (e.g. the minimum purchase price for rice aiming to protect rice price) have resulted in expanded production in northeastern China. However, the quantitative relationship between paddy rice cultivation and climate and socioeconomic factors is not well understood. The distribution and shifts of the NLPR is also rarely quantified because of a lack of accurate time-series distribution data.
Remote sensing offers one means of mapping paddy rice distribution and its change over time and space (Weiss et al 2020). Many studies have sought to map rice distributions, including in southern China in 2002 (Xiao et al 2005), southeastern China between 2001 and 2013 (Qiu et al 2016), southeastern Asia in 2002 (Xiao et al 2006), and China and India between 2000 and 2015 (Clauss et al 2016, Zhang et al 2017. The spatiotemporal distribution of paddy rice has also been examined Xiao 2016, Boschetti et al 2017). However, previous work has tended to use low to medium resolution satellite data, such as Moderate Resolution Imaging Spectroradiometer (MODIS) products, which cover a large proportion of the global paddy rice planting area at a relatively low spatial resolution (250-500 m) and with limited temporal coverage (after 2000) (Yang et al 2007). These constraints limit the ability to assess dynamic change of the NLPR.
In comparison to MODIS-like data, Landsat imagery has finer resolution (30, 60, and 80 m) and longer duration (from 1976) (Yin et al 2018). Landsat-based maps are traditionally used for localto-regional-scale studies. For example, Landsat data were used to map the extent of paddy rice in the Mekong River Delta between 2000 and 2010 (Kontgis et al 2015), and in southern China between 1990 and 2015 (Jiang et al 2018). However, these studies were mainly focused on monitoring the overall distribution of paddy rice and less effort was made to quantitatively measure the NLPR and its shifts. The NLPR may experience substantial changes in response to climate change and socioeconomic development; however, knowledge of these potential impacts is limited.
To address these gaps, we aim to: (a) determine the distribution of paddy rice using Landsat imagery in northeastern China from 1984 to 2013; (b) detect the NLPR using a kernel density estimation algorithm (Peng et al 2016) based on the generated paddy rice maps; (c) investigate the spatiotemporal characteristics of the NLPR focusing on the latitude and elevation gradients; and (d) assess the contributions of climate and socioeconomic factors in driving changes in rice area across the NLPR regions. The study will help our understanding of climate and socioeconomic impacts on paddy rice planting and diagnose

Study area
The study area is located in the northern region of Heilongjiang  • N, figure 1), Northeast China. This is the northernmost area of paddy rice cultivation in China and also one of the most northern limits of rice cultivation in the world. The area has a humid monsoonal climate with mean annual air temperature from 5.8 • C to −3.9 • C (from south to north), accumulated temperature above 10 • C from 1600 • C·d to 3600 • C·d, precipitation from 500 to 800 mm (mainly in the summer), and approximately 140-170 frostfree days a year. Agriculture is predominantly single cropping and the major crops include rice, maize, and soybeans.

Data sources
Landsat imagery from 1984 to 2013 was used to map the distribution of paddy rice and extract the NLPR information. The period 1984-2013 was chosen because paddy rice expansion mainly occurred during this period according to China Statistical Yearbook, and then we selected 1984, 1990, 2000, and 2013 as analysis years based on the different stages of paddy rice changes. All cloud-free Landsat TM, Landsat ETM+, and Landsat OLI images (30 m resolution) during the growing season (May to September) in 1984September) in , 1990September) in , 2000September) in , and 2013 were obtained from the United States Geological Survey Center for Earth Resources Observation and Science (http:// earthexplorer.usgs.gov). Where cloud-free images in a target year were not available, images from the same month of the closest available year were used instead. A total of 64 standard level-1 terraincorrected images and other processing level images were used in the analysis. To determine the spatial dynamics of the NLPR based on altitude gradient, 30 m resolution elevation data were obtained from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (DEM) V2, which were downloaded for Northeast China from the Geospatial Data Cloud site of the Computer Network Information Center, Chinese Academy of Sciences (www.gscloud.cn).
To identify the factors affecting shifts in the NLPR, we tested the impact that socioeconomic and climate factors have on rice growing in the NLPR regions from 1984 to 2011 at the county level. Rice area data from 1984 to 2013 and socioeconomic variables covering the period from 1984 to 2011 in the NLPR regions (i.e. counties covered by the NLPR) were obtained from county level statistics books, including the China Statistical Yearbook  and the China Rural Statistical Yearbook (1984-2011), which are published by China's National Statistical Bureau (NBSC 2012a(NBSC , 2012b and the Agricultural Statistics of China (1984China ( -2008 published by the Ministry of Agriculture and Rural Affairs of China (MOA 2009). Historical weather observation data from 40 national meteorological stations covering the period from 1984 to 2013 (figure 1) were acquired from the China Meteorological Data Service Center (http:// data.cma.cn), including daily mean temperature and precipitation. As there are no meteorological stations in some counties in the NLPR regions, we recreated the meteorological observation data into annual 50 km resolution spatial datasets using the inverse distance weighted interpolation method (Shepard 1968).
Rice variety data from 1984 to 2013 were obtained from the Department of Crop Production, Ministry of Agriculture and Rural Affairs of China (www.zzys.moa.gov.cn/) for the ten main cultivated rice varieties according to planting area in each year . The corresponding accumulated temperature above 10 • C for each variety was determined. The accumulated temperature above 10 • C is an important and widely used indicator of thermal conditions in crop ecology that impacts crop variety and distribution (Dong et al 2009, Hou et al 2014.

Distribution mapping of paddy rice
We used the back-propagation neural network algorithm to extract information on paddy rice distribution from the Landsat images for the study period . This algorithm has been widely adopted as an effective means of error-minimization. It consists of a multi-layer feedforward network trained by error backpropagation, which can identify useful characteristics automatically through prior knowledge. The basic theory of the back-propagation algorithm is the minimization of output error through repeated iterations. When the output error of the forward-propagation algorithm does not meet the desired output error, the back-propagation algorithm reverts the error back to the original and adjusts the weighting of each neuron in the network, using repeated iterations until the output error no longer decreases with more iterations (Rumelhart et al 1986, Li et al 2017. Mapping of the distribution of paddy rice consisted of four steps. First, a field survey was carried out in 2013, the data from which were used as the input layer. For this, a set of ground-truth data for 196 quadrats (500 m × 500 m), including 71 rice quadrats and 125 non-rice quadrats, were used as the training sample. Second, optimized parameters for the activation function, contribution threshold, regulation rate, step size, root-mean-square error, number of hidden layers, and number of iterations with higher classification accuracy were logistic, being equal to 0.5, 0.5, 0.9, 0, 1, and 30, respectively. Third, filtering analysis was used to eliminate salt-and-pepper noise and small plots of the initial paddy rice maps from the back-propagation algorithm. Finally, a set of ground quadrats randomly located in fields and at points identified from high-resolution Google Earth images were used to assess the classification accuracy using a confusion matrix. We completed the distribution mapping using ENVI.

Extracting the NLPR using kernel density estimation
Kernel density estimation is a nonparametric statistical method for estimating the probability density of a random variable and has been widely used in ecology (Seaman and Powell 1996) and geography (Peng et al 2016) to map the density surface. The kernel density estimator for bivariate data is defined as follows: where f h (x) is the probability density estimator, n is the number of data points, h (h > 0) is a user-defined smoothing parameter or bandwidth, x is a coordinate vector of estimated points, and x i is a coordinate vector of sample points. K is a user-defined non-negative kernel function, which was taken as the quadratic Epanechnikov kernel described by Silverman (Silverman 1986): The bandwidth h is a free parameter which influences the estimation result that needs to be determined. Three different bandwidths (5, 10, and 15 km) were tested. The value of 10 km was finally selected with little fragmentation for estimating the density of paddy rice distributions. We converted final paddy rice maps from 1984 to 2013 obtained in section 2.3 from raster cells to point features. Then, according to the kernel density estimation algorithm, the bandwidth (radius parameter) and resolution of the output raster were set at 10 km and 30 m, respectively, consistent with the resolution of paddy rice maps based on Landsat imagery. The kernel density maps of paddy rice distribution from 1984 to 2013 were derived and contour maps were generated for each. These processes were finalized using ArcGIS 10.2. We used a percentage-based method to derive the NLPR by accounting for the paddy-rice-field area. This employed the selected contour (same kernel density estimator) of each kernel density map of paddy rice distribution, as follows: where P is the percentage of the paddy rice planting area accounting for the selected contour of kernel density maps among the total paddy rice planting area within the study area, A KDE is the area of paddy rice fields accounting for the selected contour of each kernel density map of paddy rice distribution, and A total is the total area of paddy rice. We adopted different kernel density map contours for the paddy rice distributions (KDE = 5, 10, 15, and 20) and then calculated P for each contour during different periods and the variance of P for each contour (figure S2 (available online at stacks.iop.org/ERL/16/ 064008/mmedia)). Finally, the NLPR was identified when the kernel density map contour of paddy rice distribution was 5 in equation (1) as the maximum percentage of paddy rice in different periods (more than 98%) and the minimum variance (0.15).

Spatio-temporal dynamics of the NLPR
To calculate the spatio-temporal dynamics of the NLPR, a 5 km × 5 km fishnet was generated over the study area with the X-directional lines removed (figure 2). A set of lines were constructed in the Y-direction at an interval of 5 km, which were numbered sequentially from west to east (MarkID) as i, i + 1…etc (figure 2). We adopted the NLPR in 1984 as the baseline. The shift distances were positive or negative where the detected movement direction was to the north or south, respectively, with shifts for the period 1984-1990, 1990-2000, and 2000-2013 being detected based on equation (4): where py is the preceding year, sy is the subsequent year, and Distance ipy-sy is the distance of the NLPR shift for each period at net line i. Directions were set to 1 (−1) if the movements of the NLPR from 1984 were to the north (south). The average rate of shift between 1984 and 2013 was also determined. The NLPR elevational shifts were calculated by overlaying the elevation data. The rates of change in elevation between 1984 and 2013 were calculated using the same principles as in equation (4).

Factors influencing the shifts of the NLPR
To assess the correlation factors for shifts of the NLPR in Northeast China, we developed a linear regression model using a county-level panel dataset from 1984 to 2011. However, the NLPR is hard to be included as the dependent variable directly. Thus, we chose rice area of 22 counties where the NLPR was located from 1984 to 2013 as an alternative to the dependent variable. Because we hypothesize that the northward expansion of rice is sensitive to changes in climate, climate change will be correlated with changes in rice area in regions that represent the northernmost limits of rice cultivation in China (Li et al 2015). These climate variables included the rice growing season accumulated temperature (GST), rice growing season precipitation (GSP). Socioeconomic variables were the per capita agricultural gross national product (AgGDP), rice yield, irrigation area (irrigation), plastic mulching consumption (mulch), and machinery consumption per arable land (energy) (table S2). To assess the extent to which each chosen variable affected rice relocation to north, we built a linear regression model with a log transform: where A i,t is the rice planting area in county i in year t; X i,j,t−1 represents the socioeconomic variables, which are the per capita AgGDP, rice yield, irrigation, mulch, energy, and per capita arable land for county i in year t − 1; Climate i,j,t−1 represents the climate factors, which include the GST and GSP for county i in year t − 1; D t is a set of time dummies, representing the effects of a few major policy reforms related to rice production during the study period; α j , β j , and δ t are the estimated parameters; and ε i,t is the error term. We adopted the variance inflation factor (VIF) to avoid collinearity (VIF < 5) ( We completed this process using the standard STATA package (www.stata.com). Figure 3 shows the distribution of paddy rice in 1984, 1990, 2000, and 2013. The accuracy assessment based on high-resolution Google Earth images indicates reasonable classifications, with an overall classification accuracy of 97% and a kappa coefficient of 0.92 (table 1). The producer and user accuracy for the paddy rice category were 96% and 93%, respectively (table 1). The extent of paddy rice cultivation increased dramatically between 1984 and 2013. The paddy rice planting area was 0.24 million ha in 1984, 0.45 million ha in 1990, and 2.65 million ha in 2013. By 2000, the paddy rice area had expanded northeast and was approximately 4.58 times the area of the baseline year. By 2013, the area of cultivation had increased to approximately 11 times that of the baseline. The majority of this expansion occurred in the northeast of the study region, around the Sanjiang Plain.

Landsat-derived paddy rice distributions and the NLPR for the period 1984-2013
The NLPR exacted from the paddy rice distribution maps are shown in figure 4. In 1984, the average latitude of the NLPR was approximately 47.34 • N, and increased to 47.36, 47.51, and 47.64 • N in 1990, 2000, and 2013, respectively. The Sanjiang Plain showed the most notable shifts in the NLPR.

Patterns of shifts in the NLPR
The distribution of the NLPR varied considerably over the past three decades. We quantified the variation based on latitude and elevation measures (figures 5 and S3  figure S3(b)) during 1984-1990.
of the socioeconomic factors are highly significant (P < 0.01 or P < 0.05). The improvements in the irrigation facilities have a strongly positive impact on changes in rice area, with the highest magnitude among all of the socioeconomic factors (each 1% increase in the irrigation was expected to increase the rice area by 0.6%). The magnitude of the previous yield impact is only second to those of the irrigation improvements. As expected, the introduction of greenhouse nursery technology can significantly motivate changes in rice area (each 1% increase in the mulch was expected to increase the rice area by 0.2%). Policy implementation also strongly impacted changes in rice area (P < 0.01); more intensive policies can have higher positive impacts (more intensive policies implemented after 2004 show a higher magnitude of effects). The per capita AgGDP with a one year lag is not significant but does have a positive impact. Nevertheless, the energy use per arable land has a significantly negative correlation with the rice cultivation (each 1% increase in the energy was expected to decrease the rice area by 0.3%).

Discussion
Paddy rice cultivation expanded dramatically in Northeast China over the past three decades in the context of climate change. Our study introduces the NLPR to assess the interaction between rice planting and climate in mid to high latitude regions. The NLPR provides information regarding the extent to which climate constrains the cultivation of paddy rice and how planting has responded to climate change. This information will help policy makers track changes in agricultural resources and optimize the ongoing management of agricultural systems to ensure that future expansion is socially, economically, and environmentally sustainable into climate change.
Our study indicates that increased temperature has the largest positive association with the changes in rice area across the NLPR regions. This is consistent with the findings of previous studies (Li et al 2015), which discussed the impact of rice relocation in China. Temperature is a determining factor for paddy rice planting in the mid to high latitude regions (i.e. our study area) because rice cannot be efficiently cultivated if the temperature is below a critical threshold (Bouman et al 2007, Zhang et al 2017. In previous decades, Northeast China has experienced a continuous warming trend at a rate of 0.4 • C/decade for the mean temperature and 76.2 • C· d/decade for the rice GST (figure S1). Such warming has altered regionalscale thermal conditions, increasing the more land area suitable for paddy rice cultivation in higher latitude regions (Yin et al 2016), thus allowing the rice expansion to north across the NLPR regions and the northward shifts of the NLPR.
While climate change has allowed the northward expansion of cultivation, this expansion has also relied on a certain set of socioeconomic constraints, including benefits, policies and physical inputs. Our results confirmed that the benefits of rice production (including the previous per capita AgGDP and yield) and related policies have positive impacts on rice expansion by affecting farmers' decision making. This suggests that increased benefits will motivate farmers to expand their rice growing area in the following year , Zhang et al 2017. These expansions will require essential adaptations, such as constructing more irrigation facilities and introducing greenhouse nursery technology (aiming to avoid chilling damage during the seeding phase and increase the yield) , as our results showed irrigation and mulch have strongly positive association with changes in rice area (table 2). However, we concluded that the energy input has an opposite impact, which is inconsistent with previous study (Hu et al 2019). The negative impact of energy is probably attributed to the constraints of topography on the distribution of paddy rice cultivation, resulting in significant spatial heterogeneities in the energy input and movement of the NLPR. Higher elevation positions, such as south of the Lesser Khingan Mountains, result in more difficult field management and water access for irrigation compared with the plain regions, which suggests the need for the more intensive use of energy inputs, but less rice expansion and shifts in the NLPR. For example, Tieli county, southwest of the Lesser Khingan Mountains characterized by smaller northward shifts of the NLPR (figure 5), increased its energy input per arable land four folds (from 0.88 in 1984 to 4.21 KW ha −1 in 2011) with a four-fold increase in the rice area (from 6000 to 31 759 ha). Nevertheless, Fujin county, at the center of Sanjiang Plain characterized by significant northward shifts of the NLPR (figure 5), only increased its energy input two folds (from 0.85 in 1984 to 2.59 KW ha −1 in 2011), but had 15-fold expansion in its rice area (from 5567 to 86 939 ha).
There are two major implications from these results. First, while changes in rice area across the NLPR regions are affected by the interaction between the climate and socioeconomic factors, temperature change has more important impacts on changes in rice area across the NLPR regions. Therefore, we assume that shifts of the NLPR same as rice area are associated with temperature change. Climate warming has allowed more land at higher latitude regions to be suitable for rice cultivation. Under the motivation of benefits, humans have adopted some field measures to expand the rice area to northern regions to adapt to climate warming. Changes in the global mean surface temperature will likely range from 0.3 • C to 0.7 • C from 2016 to 2035 relative to 1986(IPCC 2014, which indicates that there are probably more new areas that will be suitable for rice expansion in higher latitude regions, such that the NLPR will probably shift to higher latitude regions to adapt to climate warming. Such adaptation of rice production to continuing climate change in the future requires more timely policies support such as increasing the budget for irrigation construction, enhancing the rice price protective system, and improving yield by agronomic innovations and new rice variety in higher latitude regions (Zhang et al 2019b;table 2). Such shifts of the NLPR will have a positive contribution to increasing rice production and improving food security.
Second, although climate warming allows for more suitable areas for rice growing, we must also consider the critical environmental consequences during rice production adaptations to climate change, which yields more challenges for achieving regional Sustainable Development Goals set by the United Nations, such as water security. For example, irrigation is positively correlated with shifts of the NLPR as paddy fields have an extended flooding period (table 2). Rice relocation resulted in huge amount of surface and ground water consumption and increased water stress in Heilongjiang Province, which will even cause the water planetary boundary (safe operating space) of Heilongjiang Province to be overstepped in 2030 (Hu et al 2020). Strategies should therefore be taken to reduce water consumption. For example, the government should improve irrigation facilities and popularize water-saving agronomic management (e.g. alternating water-saving irrigation in post-transplanting stages) under the guidelines provided in the Comprehensive Agricultural Development Program (Bryan et al 2018, Zhang et al 2019b.
Some limitations of our study must be noted. Firstly, there may be some biases in the rice distribution maps that would impact the detection of the NLPR. However, we assume that these errors are randomly distributed and have little impact on the estimation of the density of paddy rice distribution. Second, the NLPR was defined where the kernel density estimator was 5, which covered more than 98% of the paddy rice area for each period. There are still very few sporadic areas of paddy rice that remain uncovered. We assume that mapping errors were marginal as we mainly focused on the northern boundary of paddy rice planting over a relatively large area. Additionally, we excluded paddy rice planting in other mid to high latitude regions. For example, we did not include paddy rice cultivation in Russia near the Chinese-Russian border because the paddy rice planting in this area distributes rather scattered near the north of the border with a small area (approximately 299 ha in 2013; www.gks.ru/). Third, the changes in rice area across the NLPR regions could not be explained completely by the factors we considered, indicating that other important factors were not included in this study as they are hard to quantify such as renewal of rice varieties. The introduced coldtolerant rice varieties decreased the required accumulated temperature above 10 • C by 93 • C·d from 1984 to 2013 (figure S4), which is also an important motivator for rice planting in mid to high latitude regions (Chen et al 2012, Song et al 2014. More work is needed to further explore the intertwined processes of climate and socioeconomic impacts on the cropping system.

Conclusions
We introduced the northern limit of cropping planting (e.g. the NLPR) as a critical indicator of agricultural shifts in China. Time-series distribution maps with a 30 m spatial resolution were used to evaluate the impact of climate change on paddy rice cultivation in Northeast China. Our results showed that the NLPR moved towards a higher latitude and altitude. Increased temperature was strongly associated with rice expansion across the NLPR regions, which provided more suitable growing conditions. However, farmers have adapted to rice cultivation in higher latitude and altitude regions by changing their production practices, including irrigation and agricultural intensification and by adopting new technologies under the policies support. Our study shows that climate change is fundamentally altering our relationship with agriculture. The northern limit of crop planting provides a new insight for crop systems worldwide to facilitate the understanding of the impact of climate change, and the adaptation and feedback within agricultural production systems, which will be crucial for ensuring the sustainable and equitable development of food systems in the future.

Data availability statement
The data generated and/or analysed during the current study are not publicly available for legal/ethical reasons but are available from the corresponding author on reasonable request.