Climate Change and Anthropogenic Impacts on Wetland and Agriculture in the Songnen and Sanjiang Plain, Northeast China

Influences of the increasing pressure of climate change and anthropogenic activities on wetlands ecosystems and agriculture are significant around the world. This paper assessed the spatiotemporal land use and land cover changes (LULCC), especially for conversion from marshland to other LULC types (e.g., croplands) over the Songnen and Sanjiang Plain (SNP and SJP), northeast China, during the past 35 years (1980–2015). The relative role of human activities and climatic changes in terms of their impacts on wetlands and agriculture dynamics were quantitatively distinguished and evaluated in different periods based on a seven-stage LULC dataset. Our results indicated that human activities, such as population expansion and socioeconomic development, and institutional policies related to wetlands and agriculture were the main driving forces for LULCC of the SJP and SNP during the past decades, while increasing contributions of climatic changes were also found. Furthermore, as few studies have identified which geographic regions are most at risk, how the future climate changes will spatially and temporally impact wetlands and agriculture, i.e., the suitability of wetlands and agriculture distributions under different future climate change scenarios, were predicted and analyzed using a habitat distribution model (Maxent) at the pixel-scale. The present findings can provide valuable references for policy makers on regional sustainability for food security, water resource rational management, agricultural planning and wetland protection as well as restoration of the region.


Introduction
As some of the world's most productive ecosystems with their fragility and sensitivity to global changes, wetlands provide a diversity services vital for human well-being, such as health, safety and welfare, and perform a non-replaceable role in the processes of carbon and nitrogen storage, climate regulation, flood abatement, water supply and bio diversity conservation (e.g., [1][2][3]), which are now recognized as "the kidneys of the landscape", "ecological supermarkets" [4], and "act like a sponge" [5].Unfortunately, with the impact of the increasing anthropogenic and climatic changes, the degradation and loss of wetlands takes place worldwide, and is more rapid than that of other ecosystems [6].A recent review of 189 wetland assessments reported that long-term loss of global natural wetlands area averaged between 54% and 57%, but loss may have been as high as 87% since 1700 [7].Wetlands are continuing to decline globally, both in area and quality [8], particularly in Asia [7] and many midand high-latitude areas [9].Additionally, it has been frequently reported that these changes being Given the specific geography of the SNP and SJP, a case study was carried out on climate change and anthropogenic impacts on wetlands and agriculture of the region, with the definition of "wetland(s)" utilized hereafter comprised of natural or man-made ones.Natural wetlands are further divided into two categories, marshland and waterbody, while the man-made ones mainly refer to the paddy field areas, which will be described in detail in the Materials Section.The objectives of this study are mainly focused on: (1) monitoring and assessing the spatiotemporal LULCC, especially for conversion from marshland to other LULC types (e.g., croplands); (2) quantitatively distinguishing the relative role of human activities from climatic changes in terms of their impact on marshland conversions; and (3) modeling how predicted climate change scenarios will spatially and temporally impact wetlands and croplands in the study region.The reminder of this study is organized as follows: Section 2 presents the study region and the methodologies and materials utilized.The results and discussions are given in Sections 3 and 4, respectively, and the concluding remarks of the study follow in Section 5.

Study Region
The SNP, located in the central part of northeast China, lies within 42°30′ N-49°30′ N and 121°30′ E-129° E (Figure 1).The region has a total area of approximately 223,492 km 2 and belongs to the drainage of the Songhua River, as well as one of three major regions of soda saline-alkali soil in the world [31].Elevation ranges from 101 m to 1632 m above sea level, with most below 200 m.It is dominated by a semi-humid continental monsoon type climate.Annual rainfall varies geographically from about 350 to 800 mm with 70-80% of precipitation occurring in June to September.The mean monthly temperature varies from −26 to −16 °C in January to 21-23 °C in July [31].Geographically, the SJP features much in common with the SNP.It lies within 43 • 30 N-48 • 30 N and 129 • E-135 • 30 E (Figure 1), spanning a total area of about 108,829 km 2 and is an alluvial plain where three major regional rivers, namely the Heilong River, Songhua River and Ussuri River, flow through the region.Three Ramsar wetlands, namely Honghe, Sanjiang and Xingkai Lake National Nature Reserve, are distributed within the SJP [32].Before the 1950s, the region was once the largest marshland complex in China [33], but is presently a very important national grain production base after 60 years of agricultural development, primarily through transforming marshlands to arable lands.
In short, due to the large-scale anthropogenic activities (e.g., agricultural development), approximately 70-80% of the freshwater wetlands has been reclaimed during the past several decades in these two typical plains [26,34,35].Intensive anthropogenic activities also resulted enormous changes in ecological environment as increasingly frequent drought, decreased river runoff, decline of groundwater, and so on.For example, continuous overexploitation of groundwater has caused serious environmental concerns, as significant decrease of groundwater level has been found extensively in groundwater irrigated region [31,36].Meanwhile, a significant increase of air temperature with a rate of 0.3 • C per decade was detected [37], and the precipitation presented a decreasing trend in the rainy season but inverse in the dry season [38].Furthermore, significant warming has been predicted to occur under different scenarios by the end of the 21st century in the region [29,30].

LULCC Detection-Especially for the Conversion from Marshland to Other LULC Types
Understanding LULCC dynamics at regional scale based on remote sensing and GIS is not only essential to determine the coupled driving mechanisms of climatic changes and anthropogenic activities on the LULCC, but also helpful for adjusting anthropogenic activities, optimizing resources allocation for realizing regional sustainable land uses.In this study, the following approaches were adopted to explore the patterns and processes of LULCC in the SJP and SNP. (

1) Statistics on time series of LULC area and its ratio
The ratio of one LULC type in a study area is calculated as: The change of specific LULC type in the study region is calculated as: where A i and A T are the area of the i LULC type and the study region, respectively; A ic is the change areas of the i type in a certain period; and A it1 and A it2 represent the total area of the LULC type i at times t1 and t2 [39].
(2) Dynamic degree To identify systematic transitions within a matrix of LULCC, the transitions relative to the sizes of the categories must be interpreted [40].Dynamic degree (DD) [41] of a single land use type is thus introduced to analyze timely area difference of the LULC type: where A a and A b refer to area of a certain LULC type at the beginning and the ending of a period, respectively.T refers to time interval (e.g., years).

Quantitatively Distinguishing the Impacts of Anthropogenic Activities and Climatic Changes on Conversion and Loss of Marshland
Trajectory analysis can recover the history of LULCC and relate the spatiotemporal pattern of such changes to natural and anthropogenic factors [23,26,42], which can also be utilized to effectively analyze the trends of LULCC over time [43,44].This study referenced this approach and the trajectory codes were adopted to express the change trajectory of time series of LULCC as follows: where n indicates the number of time nodes and (G1) i , (G2) i , and (Gn) i are the codes of the LULC types in parcel i, at a particular time.In this study, codes 1-8 were used to represent the paddy field, dry farmland, marshland, waterbody, forest land, grassland, settlement and other unused land, respectively.Additionally, following [26], we aggregated LULCC types into three categories according to Y i : human induced changes, naturally induced changes, and unchanged land.The LULCC types were defined as human-induced as long as at least one of the observed LULCC types is recognized as human-induced, as listed in Table 1, no matter when it occurs, can distinguish human induced changes from naturally induced changes.Five periods of LULCC, i.e., 1980-2015, 1990-2015, 1995-2015, 2000-2015 and 2005-2015, in this study, were selected to study the variations of human and naturally induced contributions.

Maxent Modeling
Maxent (v.3.4.1; http://biodiversityinformatics.amnh.org/open_source/maxent/) is a machine learning software package attempting to simulate and predict how the distribution of a selected species will be modified in response to given environment changes (mostly climate and land use changes) [45] with a maximum entropy approach for species distribution simulations.It has proven useful and is widely used in habitats (e.g., wetlands and agriculture) distribution modeling [46][47][48][49][50].In this study, the model was developed by running 10 replicates with randomly splitting the distribution data into two subsets: 75% for calibrating and training the models and the reminder for testing and evaluating the model performance [51].To evaluate model projection accuracy, statistics of the area under curve (AUC) of the receiver operating characteristic (ROC) plot was used [52].AUC is a threshold independent statistic measure of model performance, varies from 0 to 1 with values >0.9 and <0.5 referring to model performance excellent and worse than a random, while values in the ranges 0.5-0.6,0.6-0.7,0.7-0.8 and 0.8-0.9indicate fail, poor, fair, and good, respectively [53].A jackknife approach [54] was used to examine the sensitivity of each individual environmental variable to the model.The immediate outputs of Maxent model include the predicted habitat map ranked with successive scale from 0 to 1.To get an effective approach to assess the suitability of the predicted target distributions and their changes, following the treatment of uncertainty for the IPCC AR5 (2013) by [55], we grouped the suitability of the predicted target distributions into three levels, including low, medium and high suitability (see Table 2).Additionally, to examine the predicted suitability changes, we overlaid the future predictions onto the current baseline prediction of suitability of target distribution and then grouped new layer of suitability values into five categories, including low suitable/no change (code 11), medium suitable/no change (code 22), high suitable/no change (code 33), increasing suitability (codes 12, 13 and 23) and decreasing suitability (codes 31, 32 and 21).

LULC Dataset
LULCC matrixes for both the SJP and SNP were generated from the seven sets of LULC maps (1980, 1990, 1995, 2000, 2005, 2010 and 2015) covering a period of 35 years derived from National Land Cover Dataset (NLCD) [56].These seven sets of LULC maps were generated through visual interpretation and digitalization on the fine resolution remotely sensed imaginaries (here the 1 km resolution maps were used) by the Chinese Academy of Sciences.For convenience of analysis, the original twenty-five LULC types were grouped into eight major LULC categories including: paddy field, dry farmland, marshland, waterbody, forest land, grassland, settlement and other unused land, as described by [57], among which, geographically, marshland is a type of natural wetland that are frequently or continually inundated with water, characterized by emergent herbaceous vegetation that is adapted to saturated soil conditions, and often found at the edges of lakes and rivers [58], while waterbody refers to the cover with obvious fresh water surface, such as rivers and lakes.Here, the cultivated lands were further sub-divided into paddy field and dry farmland.For the Maxent modeling, this paper focused on the marshland, waterbody, paddy field and dry farmland, which involves both the wetlands and agricultural land.The presence of the four LULC types were randomly created within the un-changed areas since 2000 for the individual cover.As a result, a total of 547, 159, 451 and 3002 points for four LULC types were generated within the geographic locations of the two plains.Additionally, extensive field surveys by using Unmanned Aerial Vehicle (UAV) were carried out to collect the ground truth of different wetland types (see Figure 2) in the SJP in June 2017 for gaining an intuitive knowledge of the wetlands in the specific region.

Environmental Data for the Maxent Modeling
To run the Maxent model, a set of thirty-three environmental variables, which are grouped into three categories and considered as important factors for the distribution of marshland, waterbody, paddy field and dry farmland, including nineteen bioclimatic variables, eight climatic parameters and seven topographic indexes (Table 3), was adopted for process.The nineteen bioclimatic variables, often used in habitat distribution modeling and generated from the monthly average temperature and rainfall to represent annual trends, seasonality and extremes or limiting environmental factors, were downloaded and selected from WorldClim database from the website (http://www.worldclim.org/).Please refer to [59] and [60] for the detailed explanation of WorldClim database.According to some previous studies [47,49,50], for prediction of wetlands distributions, eight climatic variables, including growing season temperature (GST), growing season precipitation (GSP), annual biological temperature (ABT) [61], coldness index (CI) [62,63], warmness index (WI) [62,63], dryness index (DI) [37], humidity index (HI) [64] and potential evapotranspiration ratio (PER) [61], were calculated based on the corresponding meteorological dataset.It should be noted that the growing season is defined as the period from May to September [50].Additionally, wetlands vary in their source and degree of wetness, and as a first-order control on spatial variation of hydrological conditions, topography affects the spatial distribution of soil moisture and groundwater flow and thus influences distribution of wetlands [65].Hence, seven terrain parameters including elevation, slope, topographic wetness index (TWI) [66], transform of aspect (Trasp), relief amplitude (RA) and elevation stand deviation (ESD), generated from 30 arc-second DEM data distributed with HydroSHED [67] database by using ArcGIS Spatial Analyst Tools, were selected in addition to the climatic data based on a priori expectations of their roles in determining locations for developing wetlands and croplands.

Environmental Data for the Maxent Modeling
To run the Maxent model, a set of thirty-three environmental variables, which are grouped into three categories and considered as important factors for the distribution of marshland, waterbody, paddy field and dry farmland, including nineteen bioclimatic variables, eight climatic parameters and seven topographic indexes (Table 3), was adopted for process.The nineteen bioclimatic variables, often used in habitat distribution modeling and generated from the monthly average temperature and rainfall to represent annual trends, seasonality and extremes or limiting environmental factors, were downloaded and selected from WorldClim database from the website (http://www.worldclim.org/).Please refer to [59,60] for the detailed explanation of WorldClim database.According to some previous studies [47,49,50], for prediction of wetlands distributions, eight climatic variables, including growing season temperature (GST), growing season precipitation (GSP), annual biological temperature (ABT) [61], coldness index (CI) [62,63], warmness index (WI) [62,63], dryness index (DI) [37], humidity index (HI) [64] and potential evapotranspiration ratio (PER) [61], were calculated based on the corresponding meteorological dataset.It should be noted that the growing season is defined as the period from May to September [50].Additionally, wetlands vary in their source and degree of wetness, and as a first-order control on spatial variation of hydrological conditions, topography affects the spatial distribution of soil moisture and groundwater flow and thus influences distribution of wetlands [65].Hence, seven terrain parameters including elevation, slope, topographic wetness index (TWI) [66], transform of aspect (Trasp), relief amplitude (RA) and elevation stand deviation (ESD), generated from 30 arc-second DEM data distributed with HydroSHED [67] database by using ArcGIS Spatial Analyst Tools, were selected in addition to the climatic data based on a priori expectations of their roles in determining locations for developing wetlands and croplands.2070)) (hereafter RCP4550, RCP4570, RCP8550 and RCP8570) was preferred over a single GCM due to large uncertainties in climate conditions associated with climate model results.The five GCMs were selected based on their historical performance at the global and regional simulations [69,70].

Other Supplementary Data
It is widely recognized that wetlands play an important role in the hydrological cycle, and in return, the hydrological conditions pose a vital impact on the developing and degradation of wetlands.Some specific hydro-meteorological data, therefore, were used to provide auxiliary information for exploring wetland dynamics in the study region, i.e., a long-term monthly precipitation (and temperature) dataset gridded in 0.5 • × 0.5 • for 1961-2015 generated and provided by the China Meteorological Administration (CMA) [71,72] based on 2472 meteorological stations across China, and an evapotranspiration (ET) simulation dataset, spanning 1982-2013, inversed from remotely sensed products based on the National Oceanic and Atmospheric Administration-Advanced Very High Resolution Radiometer satellite (NOAA-AVHRR) [73].In addition, in situ groundwater level record from 2005 to 2014 for two stations Hailun and Sanjiang in the SNP and SJP run by the Chinese Ecosystem Research Network (CERN) were collected.Crop yield data at province-level including wheat, maize, soybean and rice for the period of 1949-2015 in the Heilongjiang Province were downloaded from a set of nationwide crop statistical data from the website (http://zzys.agri.gov.cn/).

Spatiotemporal Variations of LULCC in the SJP and SNP
LULC in the SJP and SNP changed dramatically during the study period of 1980-2015.The temporal changes of main LULC types, such as marshlands and cultivated lands (including paddy field and dry farmland) over the SJP (and SNP) are illustrated in Table 4.The most notable changes of LULC in SJP were a decline in area of marshland and an increase in cultivated land.In 1980, marshland accounted for approximately 16.84% of the estimated total area of about 18,247 km 2 for the study area, while, by 2015, the total area of wetland was estimated to decrease substantially by about 47% to 9720 km 2 .Meanwhile, associated with the decrease of marshland in area during different periods, the cultivated land increased by about two-fifth from 39,049 km 2 in 1980 to 54,983 km 2 in 2015 with an annual increasing rate of 0.57-0.75% in the SJP.With respect to the SNP, similar changes occurred with the marshland area ratio decreasing from 9.48% to 8.01% and the area of the cultivated land increased from 54.08% to 59.69%.However, although the area of marshland and cultivated land changed in a smaller proportion in the SNP, they still accounted for a large area of the SNP relative to that of the SJP.Furthermore, regarding the changes of cultivated land, the paddy field had the largest increase in area with 13,870 km 2 (about 283.64% increase relative to the year 1980), followed by dry farmland (2064 km 2 , 6.04%) in the SJP.Similarly, in the SNP, the paddy field had the largest increase (6744 km 2 , 70.63%), followed by the dry farmland (5796 km 2 , 5.21%).Finally, the cultivated land took more than 50% of the whole area of both the SJP and SNP by 2015, which became the substrate and the dominant landscape of the two plains.As for the spatial distributions of marshland and other LULC types, as shown in Figures S1 and S2, a variety of the marshland constant landscape fragmentation now is mainly distributed within several large natural wetland reserves in the two plains, especially in the SJP.Furthermore, analysis of the seven periods of LULCC indicated that, during the past 35 years, different LULC types had their specific variation characteristics.For example, although the area of paddy field in the SJP was continuously increasing, except 1990-1995, a fast increase period was found during 1995-2000 (1.2% increase per year).Meanwhile, the area of dry farmland changed slightly, with characterization of increase from 1980 to 1995, but decrease gradually during the period of 1995-2015, and generally the area fluctuated from 31.52% in 1980 to 33.42% in 2015.
Statistics on LULCC (Table 4) and dynamic degree of certain LULC types (Table 5) indicated that LULC changed with time obviously in the SJP and SNP during the period 1980-2015.For the main LULC type within the SJP, paddy field area decreased during 1990-1995 and increased continuously with an accelerated rate after 2000, but the highest increase rate occurred from 1995 to 2000.The dry farmland area increased greatly at first before 1995 and then decreased at a relatively stable rate in the last two decades.Correspondingly, the marshland area reduced year by year with the decreasing rate declining gradually.As for the SNP, a similar change pattern of marshland conversion from 1980 to 2015 was observed, except the increasing period of 1995-2000, while the areal variations of paddy field and dry farmland presented large differences between the SNP and SJP.Firstly, the paddy field area increased during 1980-1995, then decreased at a declining rate from 1995 to 2005, and right after followed another increasing period in the recent decade, which revealed an irregular change with time over the SNP.In terms of the dry farmland area, it featured an increasing trend at first , but then decreasing gradually since 2005.However, as this study mainly focused on the changes of the wetlands and agricultural land, the analysis of the LULCC for the other types was beyond our scope for discussion, but several general points should be noted: (1) forest and grassland decreased significantly in both plains, with the largest decreasing rates found during the period of 1990-2000, which was considerably higher than that of the other three periods; and (2) the waterbody area in the SJP generally displayed a slightly changes from 1980 to 2015, while a large continuously decreasing trend in the SNP during 1990-2005 was observed.1 The ratio (%) is the area ratio of one LULC type in a study area calculated by Equation (1). 2 Cultivation index is the proportion of cultivation land area in the area of whole region.

The Relative Role of Anthropogenic Activities and Climatic Changes in Terms of Their Impacts on Marshland
Quantifying spatiotemporal impacts of anthropogenic activities and climatic changes on marshland changes was conducted, as shown in Table 6, and Figures 3 and 4. In the past 35 years, marshland changes attributed to the human induced changes took a predominate proportion, accounting for about 69.12% and 56.23% in the SJP and SNP, respectively, while the naturally induced changes accounted for only 8.44% and 23.70%, correspondingly.Theoretically, about 1,273,609 (8 7 −7 7 ) kinds of possible trajectories could occur, however, some trajectories never show up and some others took too small parts in all trajectories to be accounted.In fact, altogether, about 1047 kinds of trajectories existed in the SJP and the top 200 took an area about 89.60% of the region; correspondingly, about 1745 kinds of trajectories in the SNP and the top 200 took an area about 84.83% of the whole region.Among the transformations related to marshland conversion, the transformation of "marshland→dry farmland" accounted for the largest proportion about 54.40% of marshland changes in the SJP, followed by the conversion of "marshland→paddy field" with a proportion of 14.68%.However, the largest transformation (code "3222222") took about 7.97% of marshland area changes, which indicated that once the marshland was transformed to other LULC types, it would be difficult to recover.Similar to the SJP, "marshland→dry farmland" took the predominate transformation, accounted for about 35.00% of marshland area changes in the SNP, followed by the conversion of "marshland→paddy field" with a proportion of 14.30%.Therefore, the main driving force of LULCC related to marshland conversion in both plains can be attributed to anthropogenic activities, particularly the agricultural activities.Moreover, characteristic analysis on the trajectories suggested that mutual transformations between marshland and grassland as well as between marshland and waterbody occurred frequently in both plains, and these mutual transformations were mainly governed by the local precipitations.
force of LULCC related to marshland conversion in both plains can be attributed to anthropogenic activities, particularly the agricultural activities.Moreover, characteristic analysis on the trajectories suggested that mutual transformations between marshland and grassland as well as between marshland and waterbody occurred frequently in both plains, and these mutual transformations were mainly governed by the local precipitations.Table 6 lists the areas simulated in the five different periods for three types of LULCC; the corresponding area change ratio is presented in Figure 3. From this study, the following understandings concerning with marshland conversion and degradation are worth paying attention: (1) the proportion of the human-induced LULCC in the SJP was higher than that in the SNP, whereas the contributions of anthropogenic impacts on marshland conversion in the both plains were gradually diminished through five different periods, in which the most notable period occurred during 1980-2000; (2) the contribution of marshland changes related to climatic changes in LULCC in the both plains had increased gradually, particularly during the period of 1980-1995; and (3) although the naturally induced changes took a minor proportion for both plains in recent years (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015), the ratio of the unchanged marshland increased considerably, which implied that on the one hand, with the promotion of protection and restoration policy for marshland over the region, the degradation and loss of marshland has been slowed down gradually in the both plains in the recent decade, on the other hand, these results gave us confidence by using the Maxent model to explore the impacts of climate changes.
For three types of LULCC, in terms of spatial distribution, as simulation results illustrate in Figure 4, during five different periods, notable differences can be recognized between the two plains.The human induced changes were mainly concentrated in the northeast and southeast part of the SJP but scattered widely along the rivers on the right bank of the SNP, and the naturally induced changes usually appeared on the left bank of the SNP, especially in some wetland natural reserves.Generally, the decreasing rate for the human induced changes was higher in the SJP than that in the SNP.Table 6 lists the areas simulated in the five different periods for three types of LULCC; the corresponding area change ratio is presented in Figure 3. From this study, the following understandings concerning with marshland conversion and degradation are worth paying attention: (1) the proportion of the human-induced LULCC in the SJP was higher than that in the SNP, whereas the contributions of anthropogenic impacts on marshland conversion in the both plains were gradually diminished through five different periods, in which the most notable period occurred during 1980-2000; (2) the contribution of marshland changes related to climatic changes in LULCC in the both plains had increased gradually, particularly during the period of 1980-1995; and (3) although the naturally induced changes took a minor proportion for both plains in recent years (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015), the ratio of the unchanged marshland increased considerably, which implied that on the one hand, with the promotion of protection and restoration policy for marshland over the region, the degradation and loss of marshland has been slowed down gradually in the both plains in the recent decade, on the other hand, these results gave us confidence by using the Maxent model to explore the impacts of climate changes.
For three types of LULCC, in terms of spatial distribution, as simulation results illustrate in Figure 4, during five different periods, notable differences can be recognized between the two plains.The human induced changes were mainly concentrated in the northeast and southeast part of the SJP but scattered widely along the rivers on the right bank of the SNP, and the naturally induced changes usually appeared on the left bank of the SNP, especially in some wetland natural reserves.Generally, the decreasing rate for the human induced changes was higher in the SJP than that in the SNP.

Maxent Modeling Performance and Sensitivity of Environment Variables
In this study, the AUC statistic was used to assess the predicting accuracy of the Maxent model with presence-only data.The closer the AUC is to 1, the better the model performance is, and a high AUC value indicates that sites with high suitability values tend to be areas of known presence or potential distribution.Overall, the Maxent model performed well with mean AUC values ranging from 0.711 to 0.962, as shown in Table 7, at predicting the potential distributions of the four LULC types under the five climate scenarios.In which, the model performance evaluated as "excellent" (AUC ≥ 0.9) for the waterbody and paddy field, "good" (AUC = 0.8-0.9)for the marshland and "fair" (AUC = 0.7-0.8)for the dry farmland with relatively poor but acceptable performance.In brief, the different performances between dry farmland and the other LULC types might be attributed to that, unlike the paddy field (rice predominates according to [74]), marshland and waterbody mostly consist of unique LULC type with a high coincidence of the sample points in the study region, however, the majority of the dry farmland was planted with different crops, for instance, wheat, maize and soybean [74], but were not further subdivided in present study.
Permutation importance and percent contribution, measured as 0-1 in the Maxent model by using the Jackknife test, were selected to evaluate the importance of variables and help us to understand sensitivity of variables on changes of the different LULC types, as shown in Table 8.The permutation importance is up to the final performance of the model rather than the path used in an individual run and therefore is better for evaluating the importance of a particular variable [75].

Maxent Modeling Performance and Sensitivity of Environment Variables
In this study, the AUC statistic was used to assess the predicting accuracy of the Maxent model with presence-only data.The closer the AUC is to 1, the better the model performance is, and a high AUC value indicates that sites with high suitability values tend to be areas of known presence or potential distribution.Overall, the Maxent model performed well with mean AUC values ranging from 0.711 to 0.962, as shown in Table 7, at predicting the potential distributions of the four LULC types under the five climate scenarios.In which, the model performance evaluated as "excellent" (AUC ≥ 0.9) for the waterbody and paddy field, "good" (AUC = 0.8-0.9)for the marshland and "fair" (AUC = 0.7-0.8)for the dry farmland with relatively poor but acceptable performance.In brief, the different performances between dry farmland and the other LULC types might be attributed to that, unlike the paddy field (rice predominates according to [74]), marshland and waterbody mostly consist of unique LULC type with a high coincidence of the sample points in the study region, however, the majority of the dry farmland was planted with different crops, for instance, wheat, maize and soybean [74], but were not further subdivided in present study.
Permutation importance and percent contribution, measured as 0-1 in the Maxent model by using the Jackknife test, were selected to evaluate the importance of variables and help us to understand sensitivity of variables on changes of the different LULC types, as shown in Table 8.The permutation importance is up to the final performance of the model rather than the path used in an individual run and therefore is better for evaluating the importance of a particular variable [75].3.18 3.34 1 Obtained by changing the variable values between presence and background points to observe the changes of AUC. 2 represents how much the variables contributed to the model judged by the path selected for a particular simulation.
For different LULC types, environmental variables own varying effects on their changes.Among them, just as PER for dry farmland, elevation is the most powerful determining variable in the potential distribution of marshland, waterbody and paddy field and its permutation importance was much higher than others, while climatic variables were of moderate importance.However, although the terrain variables might play an important role in determining the distribution of these three LULC types predicted, it is notable that the model only based on these terrain variables did not perform as well as that by integration of climatic variables.Bio13 (precipitation of the wettest month) and Bio4 (temperature seasonality) owned the highest permutation importance for the marshland, while Bio4 and Bio3 (isothermally) played a vital role for waterbody, which indicated that changes of both temperature and humidity conditions will result in marshland changes, and waterbody changes were mostly controlled by the temperature conditions besides terrain effects.For cultivated lands, particularly the dry farmland, climatic variables (e.g., PER, HI and DI) have vital impacts, which performed a different way relative to others; however, it is not surprising given that crops planted in the dry farmland and paddy field, in other words, agricultural activities, might mostly be affected by the local weather (particularly the extreme weather) rather than terrain conditions [76].It should be noted that the variations of these variables are period specific upon the data availability.Annual mean and seasonally averaged temperature increased significantly over the past five decades with an increasing rate about 0.35 • C and 0.33 • C per decade for the SJP and SNP, respectively.The increment ratio was the highest in winter, followed by autumn, and temperature tended to increase more quickly in the SJP than that in the SNP, except spring.As for annual precipitation, both plains varied differently with no significant trend.However, it can be divided into three variation periods, i.e., the decreasing period of 1961-1980, the increasing period of 1981-1998, and the recent increasing period since 2000.Both periods with increasing trends might be mostly caused by the increase of precipitation in spring and winter.For the estimated actual evapotranspiration based on the AVHRR satellite observations, significant increasing trends in annual, summer and autumn were found in the SNP, whereas similar trends were detected in annual, winter and autumn in the SJP.It should be noted that, although the temperature and precipitation increased more significantly in the SJP compared to that of the SNP, the higher increasing trend of evapotranspiration for the SNP might be attributed to the compensation of groundwater and runoff for lack of precipitation of the region to a large extent [77].Additionally, as two major grain crops in the Heilongjiang Province where the SJP is located (about 90% grain field of the whole province distributed within the SJP), the grain yields of corn and rice varied in similar phases, especially for the period of 1995-2005, which exhibited good consistency with the transpiration from dry farmland to paddy field of the region.The groundwater levels varied in inverse ways as observed in the Sanjiang and Hailun stations, which probably attributed to the different agricultural development levels in both plains, especially the larger area of groundwater-fed paddy field in the SJP [36], where the groundwater storage has been continuously decreased in recent decade relative to that of the SNP.
Remote Sens. 2018, 10, 356 14 of 25 particularly the dry farmland, climatic variables (e.g., PER, HI and DI) have vital impacts, which performed a different way relative to others; however, it is not surprising given that crops planted in the dry farmland and paddy field, in other words, agricultural activities, might mostly be affected by the local weather (particularly the extreme weather) rather than terrain conditions [76].It should be noted that the variations of these variables are period specific upon the data availability.Annual mean and seasonally averaged temperature increased significantly over the past five decades with an increasing rate about 0.35 °C and 0.33 °C per decade for the SJP and SNP, respectively.The increment ratio was the highest in winter, followed by autumn, and temperature tended to increase more quickly in the SJP than that in the SNP, except spring.As for annual precipitation, both plains varied differently with no significant trend.However, it can be divided into three variation periods, i.e., the decreasing period of 1961-1980, the increasing period of 1981-1998, and the recent increasing period since 2000.Both periods with increasing trends might be mostly caused by the increase of precipitation in spring and winter.For the estimated actual evapotranspiration based on the AVHRR satellite observations, significant increasing trends in annual, summer and autumn were found in the SNP, whereas similar trends were detected in annual, winter and autumn in the SJP.It should be noted that, although the temperature and precipitation increased more significantly in the SJP compared to that of the SNP, the higher increasing trend of evapotranspiration for the SNP might be attributed to the compensation of groundwater and runoff for lack of precipitation of the region to a large extent [77].Additionally, as two major grain crops in the Heilongjiang Province where the SJP is located (about 90% grain field of the whole province distributed within the SJP), the grain yields of corn and rice varied in similar phases, especially for the period of 1995-2005, which exhibited good consistency with the transpiration from dry farmland to paddy field of the region.The groundwater levels varied in inverse ways as observed in the Sanjiang and Hailun stations, which probably attributed to the different agricultural development levels in both plains, especially the larger area of groundwater-fed paddy field in the SJP [36], where the groundwater storage has been continuously decreased in recent decade relative to that of the SNP.Concerning climate changes, the normalized regionally averaged TOP 5 climatic variables derived by dividing each climatic variable with the corresponding maximum of the variable value under current and four future climatic scenarios for each of the four LULCs are illustrated in Figure 6.Overall, although the TOP 5 climatic variables for predicting the distributions of the four LULCs are important, the role each variable played in the LULCC is different.For instance, although the Bio13 and CI for marshland; GST, HI and Bio5 for paddy field; PER, HI and Bio13 for dry farmland; and Bio13 for waterbody in both plains varied considerably under different climatic change scenarios, the TOP 5 climatic variables remained almost unchanged, which reflected the complex dynamics governing the LULCCs of the region.Moreover, the averaged values of climatic variables at plain scale may just give us a qualitative understanding of climate changes, the studies at much finer scales (e.g., pixel scale) are needed to assess the impacts of climate changes on LULCCs, and this will be presented in the following section.Concerning climate changes, the normalized regionally averaged TOP 5 climatic variables derived by dividing each climatic variable with the corresponding maximum of the variable value under current and four future climatic scenarios for each of the four LULCs are illustrated in Figure 6.Overall, although the TOP 5 climatic variables for predicting the distributions of the four LULCs are important, the role each variable played in the LULCC is different.For instance, although the Bio13 and CI for marshland; GST, HI and Bio5 for paddy field; PER, HI and Bio13 for dry farmland; and Bio13 for waterbody in both plains varied considerably under different climatic change scenarios, the TOP 5 climatic variables remained almost unchanged, which reflected the complex dynamics governing the LULCCs of the region.Moreover, the averaged values of climatic variables at plain scale may just give us a qualitative understanding of climate changes, the studies at much finer scales (e.g., pixel scale) are needed to assess the impacts of climate changes on LULCCs, and this will be presented in the following section.

Predicted Wetlands and Agricultural Lands Suitability under Different Climatic Change Scenarios
The area of predicted suitability of wetlands and agricultural lands at different levels under the current and four future climatic change scenarios over the SJP and SNP are presented in Figures 7  and 8, respectively.The detailed statistics and spatial distributions of these predictions are shown in Tables S2 and S3 and Figure S3, respectively.In general, the relative changes between low suitability and the sum of medium-and high-suitability represent the basic suitability variation characteristics of a specific LULC type.Results showed that predicted suitable distribution area of marshland will increase under the four climate change scenarios over these two plains, especially the medium suitability, but the trend and extent of projected changes differ between plains and scenarios.Under RCP4550, RCP8550, RCP4570 and RCP8570 scenarios, the total area of suitable distributions of marshland in the SJP (SNP) is predicted to increase by 4.98% (20.62%), 12.08% (30.86%), 9.11% (24.59%) and 13.04% (25.11%), respectively.Similar to marshland, the predicted suitability changes of waterbody area over these two plains under the four future climatic scenarios almost tripled relative to that of marshland (see Tables S2 and S3), except RCP4570 in the SJP, which will decline slightly (−1.73%).In terms of the agricultural land, under RCP4.5 scenario, paddy field in the SJP will lose by 3.81% in the 2050s, but increase by 3.60% in the 2070s, while it will decrease by 3.82% and 2.95% under RCP8550 and RCP8570 scenarios, respectively.However, statistical results reflect significant differences in the SNP, as the predicted paddy field area will increase (at least 22.16%) under these future climatic scenarios.For the dry farmland, area of its potential suitable distribution will remain unchanged or increase (mainly in the SJP, at most 3.17%) or decrease (mainly in the SNP, at most 1.38%) slightly.

Predicted Wetlands and Agricultural Lands Suitability under Different Climatic Change Scenarios
The area of predicted suitability of wetlands and agricultural lands at different levels under the current and four future climatic change scenarios over the SJP and SNP are presented in Figures 7 and 8, respectively.The detailed statistics and spatial distributions of these predictions are shown in Tables S2  and S3 and Figure S3, respectively.In general, the relative changes between low suitability and the sum of medium-and high-suitability represent the basic suitability variation characteristics of a specific LULC type.Results showed that predicted suitable distribution area of marshland will increase under the four climate change scenarios over these two plains, especially the medium suitability, but the trend and extent of projected changes differ between plains and scenarios.Under RCP4550, RCP8550, RCP4570 and RCP8570 scenarios, the total area of suitable distributions of marshland in the SJP (SNP) is predicted to increase by 4.98% (20.62%), 12.08% (30.86%), 9.11% (24.59%) and 13.04% (25.11%), respectively.Similar to marshland, the predicted suitability changes of waterbody area over these two plains under the four future climatic scenarios almost tripled relative to that of marshland (see Tables S2 and S3), except RCP4570 in the SJP, which will decline slightly (−1.73%).In terms of the agricultural land, under RCP4.5 scenario, paddy field in the SJP will lose by 3.81% in the 2050s, but increase by 3.60% in the 2070s, while it will decrease by 3.82% and 2.95% under RCP8550 and RCP8570 scenarios, respectively.However, statistical results reflect significant differences in the SNP, as the predicted paddy field area will increase (at least 22.16%) under these future climatic scenarios.For the dry farmland, area of its potential suitable distribution will remain unchanged or increase (mainly in the SJP, at most 3.17%) or decrease (mainly in the SNP, at most 1.38%) slightly.

Processes of LULCC over the Two Plains and Their Driving Forces
In general, the conversion of marshland to cultivated land predominated the long-term LULCC over the SJP and SNP.However, characteristic analysis on the long-term LULCC over these two plains suggested the complexity and discrepancies of spatiotemporal variations between these two plains.Based on previous studies in the region [22,23,26,27,78] together with main findings of this study, LULCC in different periods were further analyzed with attempt to explore the primary driving forces for the LULCC over the region.It should be noted firstly that, given the availability of continuous long-term LULC observation with a unified classification system in the study region, we mainly focused on the period after the reform and opening-up policy in 1978.
Firstly, with respect to the quantitation of the relative role of human activities from climatic changes in terms of their impacts on marshland conversions derived from this study, further in consideration of the results from [26], indicated that, on the one hand, comparing with natural factors, human activities played a predominate role (89.67%) in the conversion of marshland to the other LULC types during 1954-2013 over the SJP, especially the dramatic marshland degradation taken place in the earlier period from 1954 to 1980.On the other hand, the role human activities played tended to decreasing from 89.67% to 10.30%, while the role of climatic changes predominated the conversion from marshland to the other LULC types, which demonstrated the reliability of the fundamental assumption in predicting the effects of future climate change on LULCC.However,

Processes of LULCC over the Two Plains and Their Driving Forces
In general, the conversion of marshland to cultivated land predominated the long-term LULCC over the SJP and SNP.However, characteristic analysis on the long-term LULCC over these two plains suggested the complexity and discrepancies of spatiotemporal variations between these two plains.Based on previous studies in the region [22,23,26,27,78] together with main findings of this study, LULCC in different periods were further analyzed with attempt to explore the primary driving forces for the LULCC over the region.It should be noted firstly that, given the availability of continuous long-term LULC observation with a unified classification system in the study region, we mainly focused on the period after the reform and opening-up policy in 1978.
Firstly, with respect to the quantitation of the relative role of human activities from climatic changes in terms of their impacts on marshland conversions derived from this study, further in consideration of the results from [26], indicated that, on the one hand, comparing with natural factors, human activities played a predominate role (89.67%) in the conversion of marshland to the other LULC types during 1954-2013 over the SJP, especially the dramatic marshland degradation taken place in the earlier period from 1954 to 1980.On the other hand, the role human activities played tended to decreasing from 89.67% to 10.30%, while the role of climatic changes predominated the conversion from marshland to the other LULC types, which demonstrated the reliability of the fundamental assumption in predicting the effects of future climate change on LULCC.However,

Processes of LULCC over the Two Plains and Their Driving Forces
In general, the conversion of marshland to cultivated land predominated the long-term LULCC over the SJP and SNP.However, characteristic analysis on the long-term LULCC over these two plains suggested the complexity and discrepancies of spatiotemporal variations between these two plains.Based on previous studies in the region [22,23,26,27,78] together with main findings of this study, LULCC in different periods were further analyzed with attempt to explore the primary driving forces for the LULCC over the region.It should be noted firstly that, given the availability of continuous long-term LULC observation with a unified classification system in the study region, we mainly focused on the period after the reform and opening-up policy in 1978.
Firstly, with respect to the quantitation of the relative role of human activities from climatic changes in terms of their impacts on marshland conversions derived from this study, further in consideration of the results from [26], indicated that, on the one hand, comparing with natural factors, human activities played a predominate role (89.67%) in the conversion of marshland to the other LULC types during 1954-2013 over the SJP, especially the dramatic marshland degradation taken place in the earlier period from 1954 to 1980.On the other hand, the role human activities played tended to decreasing from 89.67% to 10.30%, while the role of climatic changes predominated the conversion from marshland to the other LULC types, which demonstrated the reliability of the fundamental assumption in predicting the effects of future climate change on LULCC.However, marshland loss and fragmentation driven by the combined effects of anthropogenic activities and climatic changes resulted in the area of marshland of the two plains tended to be reducing in general over the whole study period.
Thus far, in most studies, the LULCC, especially the dynamics of wetlands and agriculture were closely related to shifts in industry policies [23,26,34,79].Hence, we further explored the LULCC characteristics and their driving forces in terms of these aspects.Since the founding of the People's Republic of China, reclamation has been an important reason for marshland loss and degradation [26], mainly because of the rapid increase of population and the concomitant requirement of grain [34].Promoted by the policy of the agricultural modernization in the early 1980s, advanced agricultural machinery was introduced, and large modernized farms were built to produce more food [22].Since 1992, the market-directed economic system has replaced the former planned economic [80] when land reclamation was further enhanced.Since more profit could be yield from rice growing compared with dry farming during 1992-1995, the agricultural activities were much focused on the conversion from dry farmland to paddy field.With the importance of wetlands in ecological functions were widely recognized in the late 1990s, national ecological projects such as "farmland back to wetland", "construction of ecological province for Heilongjiang Province", etc. were adopted [81], during which, the wetlands got certain level of protection.However, profit oriented economic activities along with the increased stress from food security during the period of 2000-2015, resulted in marshland reclamation moved into a new rapid growth period.The conflict among food security, regional economic development and wetland protection till now remained as an unsettled issue in this region.
Meanwhile, climate change also affected the wetlands and agricultural changes in both plains.The noticeable issues have motivated many climate change studies in the region [22,23,26,34,38,78,81].In this study, the close relationships among historical climate changes, wetlands and agricultural activities were recognized, too.Examples include: (1) through affecting some eco-hydrological processes of, e.g., vegetation growing, permafrost degenerated and evapotranspiration, the continue increasing of temperature (as illustrated in Figure 5a) has stimulated the transition from marshland to croplands, because temperature is a vital climatic variable that controls the distributions of crops in the study region; and (2) except the influences of planting structural adjustment (e.g., conversion between paddy field and dry farmland), fundamentally, hydrothermal conditions over both plains took the predominate role in the conversions from wetland to agricultural land.However, increasing temperature favored the crop growth on one hand, but also increased the risk of drought, which resulted in water resource redistributions by changing irrigation from surface water to groundwater on another, especially for the SJP [36].Gradual decreases of groundwater level and storage, as shown in Figure 5d, are only the reflection of this issue, which indirectly caused the loss and degradation of marshland.
In short, because of profit oriented economic activities along with the increased stress from food security, population, economic and social development as well as climate changes, several issues related to ecological unbalance were observed in the study region, such as marshland degradation, groundwater depletion, and increase of extreme climatic events, that should be paid more attention.Consequently, this also reflected the necessity and values of further prediction of suitability of LULC under climate changes.

Suitability Changes of Predicted Wetlands and Agricultural Lands Distributions under Future Climate Changes and Its Implications
Pioneer studies exploring spatiotemporal changes of marshland as well as its related issues, such as reclamation, wetlands hydrology and reallocation of water resources over the SJP and SNP, can be found (e.g., [36,78]).However, the focus of these studies mainly concentrated on historical large-scale changes, e.g., basin-, plain-, and country-scale.Little attention has been paid to future risk of both wetlands and agriculture at the pixel-scale.Figures 9 and 10 display the temporal and spatial variations of the predicted wetlands and agricultural land suitability changes under the four future climate change scenarios relative to the "current" baseline, respectively.Tables S4 and S5 present the detailed statistics on the predicted wetlands and agriculture suitability changes for the SJP and SNP, respectively.In this study, the suitability was predicted for a specific LULC based on the topography and other geographical factors as the backgrounds, and mainly concerned with the impacts of climate change and its variability that is described in Section 2.3.2.For instance, in terms of the suitability of marshland, the change characteristics predicted under the future climate scenarios of RCP4550 and RCP4570, suggested that, on the one hand, distribution areas of all of the marshland fragments of the low-, medium-and high-suitable/no change exhibited decreasing trends, while the area of the fragments of increasing suitability was larger than that of the decreasing suitability, while, on the other hand, due to the varied changing levels of some climatic variables, i.e., Bio13, Bio4 and CI, spatiotemporal discrepancies for these variables existed over the region.Similar phenomenon was also found for the GST and HI in the prediction the paddy field distribution under RCP4550 and RCP8570.This result directly bridged the connections between changes of climatic factors and suitability distribution of LULC at the pixel-scale, and helps better understand the relationships between climatic changes and LULCC, which favor the policy making on sustainable management of ecosystems in the region.detailed statistics on the predicted wetlands and agriculture suitability changes for the SJP and SNP, respectively.In this study, the suitability was predicted for a specific LULC based on the topography and other geographical factors as the backgrounds, and mainly concerned with the impacts of climate change and its variability that is described in Section 2.3.2.For instance, in terms of the suitability of marshland, the change characteristics predicted under the future climate scenarios of RCP4550 and RCP4570, suggested that, on the one hand, distribution areas of all of the marshland fragments of the low-, medium-and high-suitable/no change exhibited decreasing trends, while the area of the fragments of increasing suitability was larger than that of the decreasing suitability, while, on the other hand, due to the varied changing levels of some climatic variables, i.e., Bio13, Bio4 and CI, spatiotemporal discrepancies for these variables existed over the region.Similar phenomenon was also found for the GST and HI in the prediction the paddy field distribution under RCP4550 and RCP8570.This result directly bridged the connections between changes of climatic factors and suitability distribution of LULC at the pixel-scale, and helps better understand the relationships between climatic changes and LULCC, which favor the policy making on sustainable management of ecosystems in the region.

Uncertainty Analysis
This study qualitatively and quantitatively assessed the effects of human activities and natural climatic changes on LULCC over the SJP and SNP, northeast China, during the past 35 years by using several LULCC detection methods and a trajectory analysis approach.The Maxent model was utilized to predict and evaluate the future suitability of main LULC types under different climate change scenarios.The utilization of a long-term seven-stage LULC dataset under the consistent classification system over the past 35 years guaranteed the quality of the important factors in the predictions.For example, for the LULCC analysis, always remember that there is not a single universally definition of wetland, and, unfortunately, precisely complete LULC data cannot be acquired because of the different definitions and techniques employed by the various assessments.For the predictions of Maxent model, on the one hand, although the LULC classification used here demonstrated the accuracy of the modeling approach, there were still some insufficiencies.For instance, the further division the dry farmland type, which may need more complicated progress.On the other hand, the selection of environment factors might also become a source of uncertainties, because there may exist overfitting results, causing by strong relationship between factors sometimes, which should be checked at first [82].Additionally, as for future projections, IPCC global climate

Uncertainty Analysis
This study qualitatively and quantitatively assessed the effects of human activities and natural climatic changes on LULCC over the SJP and SNP, northeast China, during the past 35 years by using several LULCC detection methods and a trajectory analysis approach.The Maxent model was utilized to predict and evaluate the future suitability of main LULC types under different climate change scenarios.The utilization of a long-term seven-stage LULC dataset under the consistent classification system over the past 35 years guaranteed the quality of the important factors in the predictions.For example, for the LULCC analysis, always remember that there is not a single universally definition of wetland, and, unfortunately, precisely complete LULC data cannot be acquired because of the different definitions and techniques employed by the various assessments.For the predictions of Maxent model, on the one hand, although the LULC classification used here demonstrated the accuracy of the modeling approach, there were still some insufficiencies.For instance, the further division the dry farmland type, which may need more complicated progress.On the other hand, the selection of environment factors might also become a source of uncertainties, because there may exist overfitting results, causing by strong relationship between factors sometimes, which should be checked at first [82].Additionally, as for future projections, IPCC global climate models tell us unambiguously that warming trend will continue [29], but uncertainties about its extent and pace are large [69,70].To reach a more definitive conclusion, future work must improve (e.g., regional climate simulations) [83,84].

Conclusions
In the past 35 years, anthropogenic activities, such as population expansion, socioeconomic development, and institutional policies related to wetlands and agriculture, were the main driving forces for LULCC over the SJP and SNP, while the increasing contributions of climatic change dominate the future suitable distributions of main LULC types in the study region in general.From this perspective, analyses in this paper emphasized the significance of rational utilization of wetlands and agricultural lands, which improved our understanding of the implications of wetland and agricultural ecosystem variations in undeveloped regions of northeast China.Additionally, the main findings of the predicted wetlands and agriculture suitability changes in this study contributed to policy-making by allowing for the reasonable use of land and environmental protection associated with global change and sustainability in the typical region.Therefore, under the backgrounds of ensuring food security, revitalization of old industrial base in northeast China, and the building of the Belt and Road policies in China, some actions must be taken to promote the rational use of wetland resources to balance agricultural development with wetland ecosystem sustainability for the region, which can refer to this study as a scientific basis.

Supplementary Materials:
The following are available online at www.mdpi.com/2072-4292/10/3/356/s1, Figure S1: Spatial changes of the four LULC types over the SJP, Figure S2: Spatial changes of the four LULC types over the SNP, Figure S3: Suitability distributions predicted for the wetlands and agriculture over the SJP and SNP under different climate scenarios with the Maxent model, Table S1: Bioclimatic variables, Table S2: Predicted suitable area (in square kilometers) for marshland, paddy field, dry farmland and waterbody under the current and four future climatic scenarios over the SJP, Table S3: Predicted suitable area (in square kilometers) for marshland, paddy field, dry farmland and waterbody under the current and four future climatic scenarios over the SNP, Table S4: Predicted suitable area changes (in square kilometers) for marshland, paddy field, dry farmland and waterbody under the current and four future climatic scenarios over the SJP, Table S5: Predicted suitable area changes (in square kilometers) for marshland, paddy field, dry farmland and waterbody under the current and four future climatic scenarios over the SNP.

Figure 1 .
Figure 1.Location map of the SNP and SJP showing digital elevation model (DEM), major rivers and country boundaries.

Figure 1 .
Figure 1.Location map of the SNP and SJP showing digital elevation model (DEM), major rivers and country boundaries.

Figure 3 .
Figure 3. Areal ratio variations of LULCC types over the SJP and SNP for the five different periods.

Figure 3 .
Figure 3. Areal ratio variations of LULCC types over the SJP and SNP for the five different periods.

Figure 4 .
Figure 4. Spatial variations of LULCC types for five different periods over the SJP and SNP.

Figure 4 .
Figure 4. Spatial variations of LULCC types for five different periods over the SJP and SNP.

3. 4 .
Figure 5 illustrates the variation trend lines of hydro-meteorological variables, including: (a) annual mean air temperature; (b) precipitation; (c) evapotranspiration; and (d) the observed groundwater depth variations of the Sanjiang and Hailun station, and statistics in grain yields in the Heilongjiang Province.It should be noted that the variations of these variables are period specific upon the data availability.Annual mean and seasonally averaged temperature increased significantly over the past five decades with an increasing rate about 0.35 • C and 0.33 • C per decade for the SJP and SNP, respectively.The increment ratio was the highest in winter, followed by autumn, and temperature tended to increase more quickly in the SJP than that in the SNP, except spring.As for annual precipitation, both plains varied differently with no significant trend.However, it can be divided into three variation periods, i.e., the decreasing period of 1961-1980, the increasing period of 1981-1998, and the recent increasing period since 2000.Both periods with increasing trends might be mostly caused by the increase of precipitation in spring and winter.For the estimated actual evapotranspiration based on the AVHRR satellite observations, significant increasing trends in annual, summer and autumn were found in the SNP, whereas similar trends were detected in annual, winter and autumn in the SJP.It should be noted that, although the temperature and precipitation increased more significantly in the SJP compared to that of the SNP, the higher increasing trend of evapotranspiration for the SNP might be attributed to the compensation of groundwater and runoff for lack of precipitation of the region to a large extent[77].Additionally, as two major grain crops in the Heilongjiang Province where the SJP is located (about 90% grain field of the whole province distributed within the SJP), the grain yields of corn and rice varied in similar phases, especially for the period of 1995-2005, which exhibited good consistency with the transpiration from dry farmland to paddy field of the region.The groundwater levels varied in inverse ways as observed in the Sanjiang and Hailun stations, which probably attributed to the different agricultural development levels in both plains, especially the larger area of groundwater-fed paddy field in the SJP[36], where the groundwater storage has been continuously decreased in recent decade relative to that of the SNP.

Wetlands and Agriculture 3 . 4 . 1 .
Figure 5 illustrates the variation trend lines of hydro-meteorological variables, including: (a) annual mean air temperature; (b) precipitation; (c) evapotranspiration; and (d) the observed groundwater depth variations of the Sanjiang and Hailun station, and statistics in grain yields in the Heilongjiang Province.It should be noted that the variations of these variables are period specific upon the data availability.Annual mean and seasonally averaged temperature increased significantly over the past five decades with an increasing rate about 0.35 °C and 0.33 °C per decade for the SJP and SNP, respectively.The increment ratio was the highest in winter, followed by autumn, and temperature tended to increase more quickly in the SJP than that in the SNP, except spring.As for annual precipitation, both plains varied differently with no significant trend.However, it can be divided into three variation periods, i.e., the decreasing period of 1961-1980, the increasing period of 1981-1998, and the recent increasing period since 2000.Both periods with increasing trends might be mostly caused by the increase of precipitation in spring and winter.For the estimated actual evapotranspiration based on the AVHRR satellite observations, significant increasing trends in annual, summer and autumn were found in the SNP, whereas similar trends were detected in annual, winter and autumn in the SJP.It should be noted that, although the temperature and precipitation increased more significantly in the SJP compared to that of the SNP, the higher increasing trend of evapotranspiration for the SNP might be attributed to the compensation of groundwater and runoff for lack of precipitation of the region to a large extent[77].Additionally, as two major grain crops in the Heilongjiang Province where the SJP is located (about 90% grain field of the whole province distributed within the SJP), the grain yields of corn and rice varied in similar phases, especially for the period of 1995-2005, which exhibited good consistency with the transpiration from dry farmland to paddy field of the region.The groundwater levels varied in inverse ways as observed in the Sanjiang and Hailun stations, which probably attributed to the different agricultural development levels in both plains, especially the larger area of groundwater-fed paddy field in the SJP[36], where the groundwater storage has been continuously decreased in recent decade relative to that of the SNP.

Figure 5 .
Figure 5.The observed trends of hydro-meteorological variables, variations of groundwater depth and changes in grain yields in the Heilongjiang Province.(a) the variations of annual mean air temperature observed between 1961 and 2013 over the SJP (in orange yellow color) and SNP (in light blue color), expressed as deviation from the mean during that period (in line).The dotted line is a fit to the data.The inset displayed trends in seasonal temperature (°C per year) during the period 1961-2013.A p value of <0.05 was taken as statistically significant.Levels of significance were denoted as: * 0.01 < p < 0.05, ** p < 0.01; (b,c) Same as (a) but for precipitation (1961-2013) and evapotranspiration (1982-2013) variations; (d) The variations of grain (wheat, corn, soybean and rice) fields of the Heilongjiang Province during the period 1949-2015, and the groundwater depth observations of the Sanjiang and Hailun station distributed in the SJP and SNP, respectively.

Figure 5 . 25 Figure 5 .
Figure 5.The observed trends of hydro-meteorological variables, variations of groundwater depth and changes in grain yields in the Heilongjiang Province.(a) the variations of annual mean air temperature observed between 1961 and 2013 over the SJP (in orange yellow color) and SNP (in light blue color), expressed as deviation from the mean during that period (in line).The dotted line is a fit to the data.The inset displayed trends in seasonal temperature ( • C per year) during the period 1961-2013.A p value of <0.05 was taken as statistically significant.Levels of significance were denoted as: * 0.01 < p < 0.05, ** p < 0.01; (b,c) Same as (a) but for precipitation (1961-2013) and evapotranspiration (1982-2013) variations; (d) The variations of grain (wheat, corn, soybean and rice) fields of the Heilongjiang Province during the period 1949-2015, and the groundwater depth observations of the Sanjiang and Hailun station distributed in the SJP and SNP, respectively.

Figure 6 .
Figure 6.Variations of the TOP 5 normalized climatic environment variables under different climatic scenarios for each of the four LULCs of the: (a) SJP; and (b) SNP.The normalized values were calculated by dividing each variable with the maximum of the variable under different climatic change scenarios.

Figure 6 .
Figure 6.Variations of the TOP 5 normalized climatic environment variables under different climatic scenarios for each of the four LULCs of the: (a) SJP; and (b) SNP.The normalized values were calculated by dividing each variable with the maximum of the variable under different climatic change scenarios.

Figure 7 .
Figure 7. Area changes of suitability predicted for the wetlands and agriculture under climate change scenarios over the SJP.

Figure 8 .
Figure 8. Area changes of suitability predicted for the wetlands and agriculture under climate change scenarios over the SNP.

Figure 7 . 25 Figure 7 .
Figure 7. Area changes of suitability predicted for the wetlands and agriculture under climate change scenarios over the SJP.

Figure 8 .
Figure 8. Area changes of suitability predicted for the wetlands and agriculture under climate change scenarios over the SNP.

Figure 8 .
Figure 8. Area changes of suitability predicted for the wetlands and agriculture under climate change scenarios over the SNP.

Figure 9 .
Figure 9. Predicted suitable area changes for marshland, paddy field, dry farmland and waterbody under the four future climatic scenarios over the (a) SJP and (b) SNP relative to the "current" baseline.Note: LN, low suitable/no change; MN, medium suitable/no change; HM, high suitable/no Change; IS, increasing suitability; DS, decreasing suitability.

Figure 9 .
Figure 9. Predicted suitable area changes for marshland, paddy field, dry farmland and waterbody under the four future climatic scenarios over the (a) SJP and (b) SNP relative to the "current" baseline.Note: LN, low suitable/no change; MN, medium suitable/no change; HM, high suitable/no Change; IS, increasing suitability; DS, decreasing suitability.

Figure 10 .
Figure 10.Suitability variations predicted for the wetlands and agriculture over the SJP and SNP under four future climate scenarios relative to the baseline (1970-2000) with the Maxent model.

Figure 10 .
Figure 10.Suitability variations predicted for the wetlands and agriculture over the SJP and SNP under four future climate scenarios relative to the baseline (1970-2000) with the Maxent model.

Table 2 .
Classification of the predicted changes of environmental suitability.

Table 3 .
Environment variables used in Maxent model.
Abbreviation of all the equations for calculating environmental variables, where T G is the mean temperature of the growing season; P G is the accumulated precipitation of the growing season; T i is the monthly mean temperature from January to December; P s and T s are the mean temperature and precipitation of the summer months (June to August), respectively; P is the annual precipitation; α is the local upslope area draining through a certain point per unit contour length and tan β is the local slope; Aspect represents the slope aspect ( • ).3Please refer to the website (https://hydrosheds.cr.usgs.gov/index.php).It should be noted that the bioclimatic and climatic variables were downloaded or generated based on the grid data of monthly precipitation and monthly mean temperature from the WorldClim database and all environmental variables were processed with GIS software and resampled to a grid resolution of 1 km × 1 km.

Table 4 .
Major LULC statistics in the SJP and SNP during the period 1980-2015.

Table 6 .
Areas simulated in the five different periods for three types of LULCC (km 2 ).

Table 7 .
The AUC for the four LULC types predicted with different climate scenarios by the Maxent model.

Table 7 .
The AUC for the four LULC types predicted with different climate scenarios by the Maxent model.

Table 8 .
Environmental variables (i.e., TOP 3 topographic variables (normal) and TOP 5 climatic variables (bold)) and their averaged percent contribution and percent permutation importance for each of the LULC types in the Maxent model.Variables were arranged in order as permutation importance from the highest to the lowest.