Impacts of Land Use Changes on Wetland Ecosystem Services in the Tumen River Basin

: Climate change and global rapid agricultural expansion have drastically reduced the area of wetlands globally recently, so that the ecosystem functions of wetlands have been impacted severely. Therefore, this study integrated the land use data and the integrated valuation of ecosystem services and tradeo ﬀ s (InVEST) model to evaluate the impacts of the land-use change (LUC) on wetland ecosystem services (ES) from 1976 to 2016 in the Tumen River Basin (TRB). Results reveal that the area of wetlands in TRB had decreased by 22.39% since 1976, mainly due to the rapid conversion of wetlands to dry ﬁelds and construction lands, and the LUC had induced notable geospatial changes in wetland ES consequently. A marked decrease in carbon storage and water yield was observed, while the habitat quality was enhanced slightly. Speciﬁcally, the conversion of rivers and paddy ﬁelds to ponds and reservoirs were the main reasons for the increase in habitat quality and caused the habitat quality to increase by 0.09. The conversion of marshes to lakes, paddy ﬁelds, grasslands, dry ﬁelds, and artiﬁcial surfaces were the key points for the decline in carbon storage; the conversion of marshes to lakes (5.38 km 2 ) and reservoir ponds (1.69 km 2 ) were the dominant factors driving the losses of water yield. According to our results, we should center on the conservation of wetlands and rethink the construction of the land use. The ﬁndings are expected to provide a theoretical reference and basis for promoting environmental protection in TRB and the construction of ecological civilization in border areas.


Introduction
Ecosystem services (ES) refer to the benefits that humans receive from ecosystems [1]. Wetlands are the areas of marsh, peatland or water, whether natural or artificial, permanent or temporary, with water that is static or flowing, fresh, brackish or salt, including areas of marine water whose depth at low tide does not exceed six meters [2]. Wetlands are the land use type that provides the most ES per unit area; varieties of ES provided by wetlands are of great significance for improving human welfare and promoting regional sustainable development.
Under the background of the industrial revolution and global warming, wetland degradation has been a global and tricky problem [3][4][5][6][7]. For instance, approximately 50% of wetlands have disappeared due to human interference since 1970 [8][9][10][11]. Moreover, nearly 40% of the intertidal wetlands in China have disappeared since 1990, and the losses of waterbird habitats in eastern China were as high as 19.4%

Data
Remote sensing data, meteorological data, soil data and geographic information system (GIS) auxiliary data were employed in this study. The soil data were obtained from the Harmonized world soil database (HWSD, http://westdc.westgis.ac.cn/data/tag/key/HWSD). The GIS auxiliary data included the Space shuttle radar topographic mission (SRTM) dataset from the joint survey of the National Aeronautics and Space Administration (NASA) and the National Geospatial-Intelligence Agency (NIMA) of the United States Department of Defense. Subwatershed boundaries were extracted from digital elevation model (DEM) data. The meteorological data were obtained from the China Meteorological Data Sharing Service Center (https://data.cma.cn/). The land use data were obtained from Landsat 1975/1976 Landsat MSS, and 2016/2017/2018 Landsat ETM+/OLI were downloaded from the United States Geological Survey (USGS). Moreover, the images with less than 10% cloud cover were given priority, although some images had cloud cover of more than 10%. We also employed multi-year remote sensing images to cut down the cloud effect (Table S1 in Supplementary Materials). Before wetland extraction, filling in the missing values caused by unscanned gaps on Landsat ETM+ data, radiometric, atmospheric, topographic and geometric corrections were both done through the ENVI 5.1 software package [26]. The Landsat ETM+ and OLI have a PAN band with a spatial resolution of 15 m and eight multispectral bands with 30 m resolution, so Ehlers fusion algorithm (pixel-level fusion) was applied in sharpening multispectral bands (bands 1 to 7) with a PAN band to increase the spatial resolution and reduce mixed pixels [27]. We then made image composites with bands blue, green, red and near-infrared as the optimal setting for wetland mapping.

Wetland Extraction
An object-based classification approach, known for its advantages in delineating object boundaries and reducing "salt and pepper" effects, was used for land use interpretation [28]. Firstly, we extracted wetland information with the multi-scale segmentation algorithm, which is the most common method of the image segmentation module. Referencing the polygon objects hand-digitized based on empirical values and images, the image was segmented at three scale thresholds of 100

Data
Remote sensing data, meteorological data, soil data and geographic information system (GIS) auxiliary data were employed in this study. The soil data were obtained from the Harmonized world soil database (HWSD, http://westdc.westgis.ac.cn/data/tag/key/HWSD). The GIS auxiliary data included the Space shuttle radar topographic mission (SRTM) dataset from the joint survey of the National Aeronautics and Space Administration (NASA) and the National Geospatial-Intelligence Agency (NIMA) of the United States Department of Defense. Subwatershed boundaries were extracted from digital elevation model (DEM) data. The meteorological data were obtained from the China Meteorological Data Sharing Service Center (https://data.cma.cn/). The land use data were obtained from Landsat 1975/1976 Landsat MSS, and 2016/2017/2018 Landsat ETM+/OLI were downloaded from the United States Geological Survey (USGS). Moreover, the images with less than 10% cloud cover were given priority, although some images had cloud cover of more than 10%. We also employed multi-year remote sensing images to cut down the cloud effect (Table S1 in Supplementary Materials). Before wetland extraction, filling in the missing values caused by unscanned gaps on Landsat ETM+ data, radiometric, atmospheric, topographic and geometric corrections were both done through the ENVI 5.1 software package [26]. The Landsat ETM+ and OLI have a PAN band with a spatial resolution of 15 m and eight multispectral bands with 30 m resolution, so Ehlers fusion algorithm (pixel-level fusion) was applied in sharpening multispectral bands (bands 1 to 7) with a PAN band to increase the spatial resolution and reduce mixed pixels [27]. We then made image composites with bands blue, green, red and near-infrared as the optimal setting for wetland mapping.

Wetland Extraction
An object-based classification approach, known for its advantages in delineating object boundaries and reducing "salt and pepper" effects, was used for land use interpretation [28]. Firstly, we extracted wetland information with the multi-scale segmentation algorithm, which is the most common method of the image segmentation module. Referencing the polygon objects hand-digitized based on empirical values and images, the image was segmented at three scale thresholds of 100 (Level 1), 80 (Level 2) and 50 (Level 3), respectively. Specifically, objects created at Level 1 were used to detect changes for large areas of vegetation and non-vegetation. Objects at Level 2 were created to identify changes for wetland, farmland and grassland. Finally, those at Level 3 for construction land and bare land. Meanwhile, the shape factor was set to 0.2, the compactness parameter was set at 0.4, and the smoothness was 0.6. Then, the land use classification decision tree was established according to the reference [25]. It is worth noting that the wetland patch maps from 2016 were provided by various counties and municipalities in the TRB, and we also combined 480 field sampling points that come from Google Earth's field measurements and high-resolution remote sensing data with the establishment of stratified random sampling method [29] to verify the classification accuracy. Additionally, the verification points from 1976 were selected from the wetland patches that were virtually unchanged from 1976 until 2016, according to the wetland data in 2016 and fieldwork (e.g., forest swamp in the mountains, natural river). This process led to the final land use data required for this study ( Figure S1). Furthermore, all land use maps were uniformly resampled to 250 m after the classification and improved with field wetland patch data provided by various forestry bureaus in TRB to keep the same resolution with other data employed in InVEST. ArcGIS software was used to perform conversion analysis between various classes, providing data support for the subsequent analysis of the impacts of LUC on ES changes.

Assessment of Wetland ES
With reference to the research results of millennium ecosystem assessment(MEA) [3] and in combination with the availability of data, this article selected three ES types as the key wetland ES in this area: habitat quality, carbon storage and water production.
Habitat Quality. The habitat quality module generates habitat quality maps that combine information on land use and biodiversity threat factors. The principle is to obtain the distribution of habitat quality by combining the sensitivity of the landscape types and the intensity of external threats and use the habitat index to reflect the quality of habitats [30]. The specific calculation formula is as follows: where Q xj represents the quality of habitat in grid cell x in LULC type j; D xj represents the total threat level in grid cell x with land use or habitat type j; H j represents the habitat suitability of land use type j; k is a half-saturation constant and usually taken as half of the maximum value, and z is a normalized constant and usually set to 2.5 [31]. Referring to the studies of Sharp et al. [30], Sallustio et al. [32] and Bao et al. [33], parameters for the various land types were obtained (Tables S2 and S3). Carbon Storage. The carbon storage and sequestration model, commonly used to estimate the amount of carbon currently stored in a landscape or the amount of carbon stored over time, is mainly composed of four carbon pools: aboveground biomass, belowground biomass, soil and dead organic matter. The calculation formula is as follows: where C tol is the total carbon storage in the basin (t), C above is the aboveground carbon storage (t), C below is the belowground carbon storage (t), C soil is the soil carbon storage (t) and C dead is the dead organic matter carbon storage (t).
where C j represents four different carbon pools, C ij represents the carbon density of the corresponding land use type i of carbon pool j, S ij represents the area of corresponding land use type i of carbon pool j, and n represents the number of land-use types. Carbon density data were obtained from the studies and field surveys of Xiang et al. [25] and Yan et al. [34] (Table S4). Water Yield. The water yield model is based on the water cycle and the Budyko curve. The calculation formula is as follows: Y xj is the annual water yield for each grid cell x on landscape j; P x is the annual precipitation on grid cell x and AET xj is the actual annual evapotranspiration for grid cell x on landscape j. The calculation of AET xj P x was performed using the Budyko water-energy coupled balance equation proposed by Fu [35] and Zhang et al. [36]. The calculation formula is as follows: R xj is the drying index (ratio of actual evapotranspiration to precipitation) of grid cell x when the land use type is j and x is the ratio of the corrected annual available vegetation to the expected precipitation. The calculation formulas are as follows: k is the vegetation evapotranspiration coefficient, ET 0 is potential evapotranspiration (mm/a), Z is the Zhang coefficient [31] and AWC x is the effective moisture content of the soil; the calculation formula is as follows: RA is solar radiation at the top of the atmosphere (astronomical radiation) (MJ · m −2 · d −1 ), T av is the average daily maximum temperature and daily minimum temperature ( • C), TD is the average difference between daily maximum temperature and daily minimum temperature ( • C), Max Soil Depth x is the maximum thickness of the soil layer, Root Depth x is the root depth (Table S5) and PAWC x is the available moisture content for vegetation. The datasets used for the water yield model input included climate data, digital elevation models, soil data, land use data and watershed boundary data [37,38].

Evaluation of LUC Data Accuracy
Through the combination of field data and high-resolution remote sensing data downloaded from satellites, woodland, grassland, farmland, wetland, construction land and bare land were finally extracted. The research results showed that the overall accuracy of the result in 2016 was 92.40; the Kappa coefficient was 0.91. The overall accuracy in 1976 was 88.55; the Kappa coefficient was 0.87 (Table 1). The lower accuracy of the result in 1976 could be ascribed to (1) less spectral bands for classification and (2) low quality of remote sensing images from 1976, while the 2016 Landsat OLI/ETM+ data contained a panchromatic band with a resolution of 15 m, which improved the resolution of the remote sensing image data and made the wetland classification accuracy in 2016 the highest. Moreover, it is difficult to classify the wetland accurately in land use data. Field wetland patch data provided by various forestry bureaus in TRB were collected to improve the wetland classification accuracy, which was 83.72-86.14% in the work of Zheng et al. [18]. Hence, the overall accuracy of the classification results was higher than 88%, and the Kappa coefficient was greater than 0.85, which met the classification requirements. Therefore, the results extracted in this study could objectively and accurately reflect the LUC process in TRB.

Changes in Wetland Area from 1976 to 2016
Remarkable and continuous wetland losses were examined in the TRB from 1976 to 2016. The wetland area reduced from 1161.80 km 2 in 1976 to 901.62 km 2 in 2016, with a decrease of 22.39% and an average annual reduction of 6.50 km 2 ( Figure 2). The most extensive loss was observed in BEH, whose wetland area decreased from 245.26 km 2 in 1976 to 158.57 km 2 in 2016, i.e., a reduction of 86.69 km 2 , accounting for 33.32% of the total wetland area loss in the area. At the same time, the wetland reduction in LDH was the least. From 1976 to 2016, LDH wetland area decreased by 1.22 km 2 , accounting for only 0.47% of the total wetland reduction in the area. In addition, the wetland area in GPG increased from 25.10 km 2 in 1976 to 34.18 km 2 in 2016, a 9.078 km 2 increase.
Sustainability 2020, 11, x FOR PEER REVIEW 6 of 15 0.85, which met the classification requirements. Therefore, the results extracted in this study could objectively and accurately reflect the LUC process in TRB.

Changes in Wetland area from 1976 to 2016
Remarkable and continuous wetland losses were examined in the TRB from 1976 to 2016. The wetland area reduced from 1161.80 km 2 in 1976 to 901.62 km 2 in 2016, with a decrease of 22.39% and an average annual reduction of 6.50 km 2 ( Figure 2). The most extensive loss was observed in BEH, whose wetland area decreased from 245.26 km 2 in 1976 to 158.57 km 2 in 2016, i.e., a reduction of 86.69 km 2 , accounting for 33.32% of the total wetland area loss in the area. At the same time, the wetland reduction in LDH was the least. From 1976 to 2016, LDH wetland area decreased by 1.22 km 2 , accounting for only 0.47% of the total wetland reduction in the area. In addition, the wetland area in GPG increased from 25.10 km 2 in 1976 to 34.18 km 2 in 2016, a 9.078 km 2 increase. In detail, the river area in TRB reduced the most in the past forty years. The loss of the river has reached 157.62 km 2 from 1976 to 2016, accounting for 57.93% of the total wetland loss in the area ( Table 2). Lake area, by contrast, exhibited the least reduction. The reduction was only 17.47% from 1976 to 2016, accounting for 0.59% of the total wetland loss in the area. Additionally, the increase of the pond and reservoir zones was observed significantly, from 8.60 km 2 in 1976 to 76.69 km 2 in 2016, an increase of 68.09 km 2 . In detail, the river area in TRB reduced the most in the past forty years. The loss of the river has reached 157.62 km 2 from 1976 to 2016, accounting for 57.93% of the total wetland loss in the area ( Table 2). Lake area, by contrast, exhibited the least reduction. The reduction was only 17.47% from 1976 to 2016, accounting for 0.59% of the total wetland loss in the area. Additionally, the increase of  For a specific assessment of wetland area change among all subwatersheds in the TRB, the changes for the past forty years over all the subwatersheds were compared in Table 3. The results of temporal changes in wetland area from 1976 to 2016 showed a decreasing trend, except in the reservoir/pond and lake landscapes. Particularly, the reservoir/pond in HCH showed an increase of 37.25%, while the lake among all the subwatersheds remained virtually unchanged. Furthermore, the area of swamp wetland declined seriously with a change percentage of more than 40% among GPG, GYH and JXQ, and the greatest reduction was shrub swamp (45.13%), forest swamp (49.96%) and marsh (63.22%), respectively. River degradation was obviously observed with the amount of 34.46 km 2 , 41.38 km 2 in BEH and GYH, respectively. Paddy exhibited the largest decrease in BEH, with a decrease of 47.08%.

Changes in Wetland ES
In 2016, the average habitat quality in the TRB was 0.81. For a specific assessment of subwatersheds, the habitat quality in GPG was the highest, 0.91, while it was the lowest in MJH, 0.31 ( Figure 3A). Moreover, the carbon storage, on average, was 56.06 Tg (1 Tg = 10 6 t), and the subwatershed GYH showed the highest carbon storage, 18.05 Tg, accounting for 32.19% of the total carbon storage in the area, while YQG had the lowest carbon storage, 0.66 Tg, accounting for 1.17% of the total carbon storage in the TRB ( Figure 3B). The total amount of water yield was 217.33 × 10 5 m 3 in 2016. HCH had the highest water yield, 44.36 × 10 5 m 3 , accounting for 20.41% of the total water yield in the TRB, and MJH had the lowest water yield, 4.41 × 10 5 m 3 , accounting for 2.03% of the total water yield in the area ( Figure 3C). The pattern of these ES was largely related to the distribution of wetland types. of the total carbon storage in the TRB ( Figure 3B). The total amount of water yield was 217.33 × 10 5 m 3 in 2016. HCH had the highest water yield, 44.36 × 10 5 m 3 , accounting for 20.41% of the total water yield in the TRB, and MJH had the lowest water yield, 4.41 × 10 5 m 3 , accounting for 2.03% of the total water yield in the area ( Figure 3C). The pattern of these ES was largely related to the distribution of wetland types. The results of temporal changes in wetland ES from 1976 to 2016 showed a decreasing trend, except habitat quality, which was improved slightly. Particularly, the habitat quality in the TRB increased from 0.81 to 0.90, an increase of 11.11%. The habitat quality in LDH was enhanced by 0.17, while that in HQH and MJH both increased only 0.01 ( Figure 4A  A striking reduction was identified in GYH for carbon storage with an amount of 4.82 Tg, accounting for 25.09% of the total carbon storage loss in the area, and the carbon storage in YQG declined slightly with the change of about 0.11 Tg, accounting for 0.57% of the total carbon storage loss in the TRB ( Figure 4B). A significant decrease of 190.28 × 10 5 m 3 from 1976 to 2016 was observed. Particularly, the amount of water yield in BEH showed the largest reduction, 70.58 × 10 5 m 3 , The results of temporal changes in wetland ES from 1976 to 2016 showed a decreasing trend, except habitat quality, which was improved slightly. Particularly, the habitat quality in the TRB increased from 0.81 to 0.90, an increase of 11.11%. The habitat quality in LDH was enhanced by 0.17, while that in HQH and MJH both increased only 0.01 ( Figure 4A). The amount of carbon storage declined from 75.29 Tg in 1976 to 56.06 Tg in 2016, a decrease of 25.54% and an average annual reduction of 1.37 Tg. of the total carbon storage in the TRB ( Figure 3B). The total amount of water yield was 217.33 × 10 5 m 3 in 2016. HCH had the highest water yield, 44.36 × 10 5 m 3 , accounting for 20.41% of the total water yield in the TRB, and MJH had the lowest water yield, 4.41 × 10 5 m 3 , accounting for 2.03% of the total water yield in the area ( Figure 3C). The pattern of these ES was largely related to the distribution of wetland types. The results of temporal changes in wetland ES from 1976 to 2016 showed a decreasing trend, except habitat quality, which was improved slightly. Particularly, the habitat quality in the TRB increased from 0.81 to 0.90, an increase of 11.11%. The habitat quality in LDH was enhanced by 0.17, while that in HQH and MJH both increased only 0.01 ( Figure 4A  A striking reduction was identified in GYH for carbon storage with an amount of 4.82 Tg, accounting for 25.09% of the total carbon storage loss in the area, and the carbon storage in YQG declined slightly with the change of about 0.11 Tg, accounting for 0.57% of the total carbon storage loss in the TRB ( Figure 4B). A significant decrease of 190.28 × 10 5 m 3 from 1976 to 2016 was observed. Particularly, the amount of water yield in BEH showed the largest reduction, 70.58 × 10 5 m 3 , A striking reduction was identified in GYH for carbon storage with an amount of 4.82 Tg, accounting for 25.09% of the total carbon storage loss in the area, and the carbon storage in YQG declined slightly with the change of about 0.11 Tg, accounting for 0.57% of the total carbon storage loss in the TRB ( Figure 4B). A significant decrease of 190.28 × 10 5 m 3 from 1976 to 2016 was observed. Particularly, the amount of water yield in BEH showed the largest reduction, 70.58 × 10 5 m 3 , accounting for 37.14% of the total reduction in water yield in the area, and STH had the smallest reduction, 1.20 × 10 5 m 3 , accounting for 0.63% of the total reduction in water yield in the area. In addition, HQH showed an increasing trend in water yield, increasing by 2.30 × 10 5 m 3 ( Figure 4C).
For a specific assessment of ES in different types of wetlands in the TRB, the changes for the past forty years in the three ES over all the subwatersheds were compared in Figure 5. In terms of forty years in the three ES over all the subwatersheds were compared in Figure 5. In terms of habitat quality, ponds and reservoirs were obviously enhanced with the increasing index of about 0.48 ( Figure 5A), and forest swamps exhibited the largest decrease with 0.06.
In terms of carbon storage, the ES in rivers and marshes dramatically decreased by 9.06 Tg and 6.84 Tg, respectively, accounting for 47.13% and 35.56% of the total carbon storage loss in the area ( Figure 5B). The carbon storage in ponds and reservoirs showed an increasing trend, with a total increase of 4.01 Tg. In terms of water yield, paddy fields exhibited the largest reduction, 126.60 × 10 5 m 3 , accounting for 66.53% of the total loss of water yield in the area. Ponds and reservoirs showed an increasing trend in water yield, with a total increase of 9.75 × 10 5 m 3 ( Figure 5C).

Discussion
Remarkable LUC-driven sharp conversions of ES were observed in this region during the period of 1976 to 2016. Although the conversion of wetlands, including shrub swamps and rivers, to dry fields and artificial surfaces were the main reasons for the degradation of habitat quality in the TRB, making the habitat quality reduced by approximately 0.08-0.17, the conversion of rivers and paddy fields to ponds and reservoirs enhanced the habitat quality by approximately 0.15. In terms of carbon storage, the ES in rivers and marshes dramatically decreased by 9.06 Tg and 6.84 Tg, respectively, accounting for 47.13% and 35.56% of the total carbon storage loss in the area ( Figure 5B). The carbon storage in ponds and reservoirs showed an increasing trend, with a total increase of 4.01 Tg. In terms of water yield, paddy fields exhibited the largest reduction, 126.60 × 10 5 m 3 , accounting for 66.53% of the total loss of water yield in the area. Ponds and reservoirs showed an increasing trend in water yield, with a total increase of 9.75 × 10 5 m 3 ( Figure 5C).

Discussion
Remarkable LUC-driven sharp conversions of ES were observed in this region during the period of 1976 to 2016. Although the conversion of wetlands, including shrub swamps and rivers, to dry fields and artificial surfaces were the main reasons for the degradation of habitat quality in the TRB, making the habitat quality reduced by approximately 0.08-0.17, the conversion of rivers and paddy fields to ponds and reservoirs enhanced the habitat quality by approximately 0.15.
The conversion of marshes to lakes (0.03 Tg), paddy fields (0.04 Tg), grasslands (0.01 Tg), dry fields (0.24 Tg) and artificial surfaces (0.03 Tg) induced a significant loss in carbon storage, respectively, accounting for 4.55%, 6.06%, 1.52%, 36.36% and 4.55% of the total loss of carbon storage. In addition, the conversion of rivers to dry fields and artificial surfaces should also be noted, which induced a remarkable loss of carbon storage by 0.09 Tg and 0.03 Tg, respectively.
The conversion of farmland and construction land to reservoirs/ponds, rivers and paddy fields and the conversion of forests to reservoirs/ponds and rivers also cannot be ignored as they were the dominant forces driving the losses of water yield in the TRB.
Although the conversion of rivers (93.75 km 2 , 19.13 km 2 ) and paddy fields (162.56 km 2 , 46.75 km 2 ) to farmlands and construction lands increased the water yield by 199.69 × 10 5 m 3 , the reduction in water yield due to the conversion of non-wetland to wetland was greater than the former. Meantime the reduction in water yield caused by the conversion of marshes to lakes (5.38 km 2 ) and reservoirs and ponds (1.69 km 2 ) were also significant, and the water yield reductions were 2.81 × 10 5 m 3 and 0.76 × 10 5 m 3 , respectively.
For a specific assessment of the policy-driven conversions of land use and the LUC-driven wetland ES change, we know that the change in habitat quality was mainly related to development trends in the agricultural economy in the TRB. Typically, China, Russia, North Korea, South Korea and Mongolia jointly launched the Tumen River Regional Cooperative Development project in 1990, which gradually increased the pressure on local agricultural and socio-economic development. Meanwhile, 11.31 km 2 , 294.13 km 2 and 70.19 km 2 of wetlands were converted into paddy fields, farmland and artificial surfaces, respectively, for meeting the demand of the growing population and the increasing grain supply-demand. Moreover, climate change has also been recognized as one of the key drivers of wetland degradation globally [39,40].
The precipitation in the TRB showed a decreasing trend [41], which caused many rivers and lakes to dry up during the periods 1976-1988, 1992-1994 and 2005-2015. They were thereby converted to agricultural land, which reduced carbon storage and the water yield. In addition, 17.06 km 2 of other types of wetlands were converted into reservoirs or ponds, which increased habitat quality.
The main reason for the decrease in carbon storage of the wetland ecosystem in the TRB was that 22.51 km 2 of swamp were drained and cultivated as agricultural land in recent years. The low soil and water retention capacity of crops caused serious losses in soil organic carbon ( Figure S2).
The construction of large hydrological infrastructures, irrigation facilities, roads, ditches, bridges and culverts changed the soil organic carbon content of the underlying surface and reduced the diversity of wetland vegetation in the TRB [42]. It even changed the hydrological conditions of the wetlands, which exacerbated the degradation of marshes and caused a significant decline in carbon storage.
In order to actively maintain and improve wetland ES, the State Forestry Administration of the People's Republic issued some documents, for example, 'wetland ecological benefit compensation pilot project' in 2009, so clearing the policy-driven LUC and the responses of wetland to these policies, which has enhanced diverse wetland ES, are especially necessary. Consequently, 13 local and national forest parks, nature reserves and wetland parks (Table 4), i.e., Gayahe National Wetland Park, Lanjia Grand Canyon National Forest Park, Jilin Hunchun National Nature Reserve, were established for protecting the wetlands in the TRB. We found that these conservation policies in the TRB contributed significantly to the increase in wetland ES. Consequently, habitat quality was slightly enhanced in the natural reserve in MJH, which was 0.51 higher than that across the whole MJH. In addition, the carbon storage of wetlands in the natural reserve in GYH in 2016 was 0.75Tg, accounting for 4.17% of the amount across the entire GYH. The water yield in the natural reserves in GPG, JXQ and HCH accounted for 84.96%, 55.34% and 20.53% of the total amount in each subwatershed because of the establishment of Tumen River of the Legal Protected National Forest Park (established in 2002), Jilin Hunchun Amur Tiger National Nature Reserve (established in 2005) and Tumen River National Forest Park (established in 1997), accounting for 66.73%, 26.89% and 40.86% of the total wetland area in each subwatershed, respectively. Meanwhile, we also found that some wetlands outside the natural reserve have not been effectively protected. Therefore, we suggest wetlands that have high ecosystem services values, especially outside the protected areas, should be taken as the key conservation areas in the next step so that we can slow down the distinct degradation that negative policy-driven LUC had induced on wetlands.
LUC has been identified as one of the most influential factors for the change of ES [4,43]. Therefore, it is necessary to take the protection and restoration of wetlands into consideration from the perspective of LUC [44]. Moreover, maintaining and improving the wetland ecosystem functions and services is also a key step for wetland management [45][46][47]. We propose that wetland nature reserves should be established in BEH and JXQ, making wetlands in these subbasins can become well developed. As the salmon spawning backflow route in MJH, a national natural reserve area for conserving fish and other aquatic animals should be established at the same time. Finally, we suggest that all wetland conservation areas should be extracted based on wetland ES in the future.
There also are some limitations to this study. Although the land use data used in this paper was combined with the wetland patch field survey data from the counties and cities from 2016, the wetland area obtained by remote sensing still included some uncertainties in terms of accuracy. Moreover, the impacts of LUC on ES were not evaluated at a mechanism level in-depth. However, these restrictions did not have had a fundamental impact on our results, which combined wetland patch survey data and multi-source remote sensing data to evaluate the impacts of LUC on wetland ES in the TRB at multiple scales. Therefore, we will integrate the optical images with different spectral characteristics and the microwave remote sensing data with field survey methods to characterize wetland spatial information accurately in the future. Second, make efforts to integrate ecological process-based models and evaluate the impacts of LUC on wetland ES from the mechanism perspective. Although we quantitatively evaluated wetland ES of different patch areas on a large scale and calculated the spatial land use differences in ES correctly with InVEST, and while the explicit model results demonstrated the response of ES to LUC will still find some issues. For example, the model application, based on a simplified statistical analysis of ecosystem processes, still presents some limitations that cannot be ignored. In the "habitat quality" model, all threats are additive, which makes the quality value consequently reduce. However, the real impact of multiple threats is much greater than the sum of individual threat levels linearly in most ecosystems [48]. As for the modeling approaches in "carbon storage", it does not include greenhouse gas emissions related to LUC, which will affect the net carbon sequestration capacity [49]. Moreover, the "water yield" model, based on a simplification of the hydrological cycle, does not consider the base flow. The model results are reliable only at subwatershed to watershed scale, where the model is developed [50]. While considering the aim is to compare LUC effects on wetland ES at basin level, rather than to provide a fairly accurate estimate, so these modeling limitations can be acceptable when the ES is modeled over wide temporal ranges. Generally, the findings provide a theoretical reference and basis for promoting environmental protection in the TRB and the construction of ecological civilization in border areas. We will continue to improve data accuracy in the future, making models that integrate more ecological mechanisms, which will be a valuable tool for the description of precise, punctual ES.

Conclusions
We analyzed the impacts of LUC on wetland ES in the TRB from 1976 to 2016. We found that the area of wetlands in the TRB showed a rapid diminishing trend. The marsh, in particular, had been cut into half. Moreover, among the eleven subwatersheds, the wetland area in BEH decreased the most, while LDH decreased the least.
As a result of LUC, wetland ES showed a decreasing trend in the past forty years, except for habitat quality, which showed an increase of 11.11% due to the conversion of rivers and paddy fields to ponds and reservoirs. Moreover, the habitat quality in LDH increased the most. The conversion of marshes to lakes, paddy fields, grasslands, dry fields and artificial surfaces were the main reasons for the decline in carbon storage, causing a decrease in carbon storage of approximately 19.23 Tg. Carbon storage in GYH showed the largest reduction, while YQG had the lowest carbon storage reduction. The conversion of non-wetland to wetland was the main reason for the decline in water yield in the TBR basin, a decrease of 190.28 × 10 5 m 3 from 1976 to 2015, in which the BEH had the largest reduction.
Supplementary Materials: The following are available online at http://www.mdpi.com/2071-1050/12/23/9821/s1, Figure S1: Wetland in the TRB. (A) is the wetland use data in 1976, and (B) is the wetland use data in 2016., Figure S2: Status of drainage and reclamation of marsh wetlands in the TRB, Table S1: The information list  of landsat image data: title, Table S2: Threators, Table S3: Sensitivity Table, Table S4: Land use and carbon density(t/ha), Table S5: Land use and carbon density(t/ha).
Author Contributions: Conceptualization, W.Z. and D.Z.; methodology, Y.Z. and R.J.; software, Y.Z. and X.Z.; writing-original draft preparation, Y.Z.; writing-review and editing, D.Z. and R.J.; supervision, project administration, W.Z. All authors have read and agreed to the published version of the manuscript.