Spatial Evaluation of Greenhouse Gas Fluxes in a Sasa (Dwarf Bamboo) Invaded Wetland Ecosystem in Central Hokkaido, Japan

: To evaluate the effect of vegetation change on greenhouse gas (GHG) budget from a wetland ecosystem, the CO 2 , CH 4 and N 2 O budgets from whole area (21.5 ha) of the Bibai Wetland, where dwarf bamboo ( Sasa ) or Ilex has invaded into original Sphagnum dominated vegetation, located in Hokkaido, Japan were estimated. The original Sphagnum -dominated vegetation was changed from a sink to a source of CO 2 by invasion of short-Sasa (50 cm > height), while the invasion of tall-Sasa (50 cm < height < 150 cm) or Ilex increased CO 2 uptake. Annual CH 4 emission was decreased by the invasion of Sasa or Ilex . The annual N 2 O emission was slightly increased by invasion of Ilex only. These GHG budgets were correlated with the environmental factors related to the water table depth. The distribution of vegetation and environmental factors was estimated from satellite image bands, and the GHG budget of the entire wetland was estimated. The whole wetland area was considered to be a sink for GHG ( − 113 Mg CO 2 -eq y − 1 ) and CO 2 uptake by tall-Sasa occupied 71% of the GHG budget. The vegetation change due to the lowering of the water table depth currently increases the rate of carbon accumulation in the ecosystem by about 5 times.


Introduction
The area of peat-forming wetlands is estimated to be 230-250 Mha globally [1,2], which is about 1.5% of the total terrestrial area.Wetland ecosystems in their natural state generally accumulate large amounts of carbon (C) in the form of peat over long periods of time, as plant production exceeds organic matter decomposition due to their wet environment [3].The amount of carbon stored in northern peatlands is estimated to be 450 Pg (petagram, 10 15 g), which is one third of the global soil C storage [1].While peatlands absorb carbon dioxide (CO 2 ) and play an important role as a sink for C, they also produce and emit methane (CH 4 ) through the anaerobic decomposition of organic matter.CH 4 emission from natural wetlands are considered to be an important source of CH 4 [4].The release of nitrous oxide (N 2 O) has been shown to be small in wetlands under natural conditions, but can be a large source when subjected to anthropogenic drainage and drying [5].
Greenhouse gas (GHG) fluxes from wetlands are highly spatio-temporally variable and are influenced by a variety of factors, including environmental factors such as air and soil temperature, moisture condition, and solar radiation, as well as soil chemical properties and vegetation.The net ecosystem exchange of CO 2 (NEE) in wetlands is determined by the balance between CO 2 uptake by photosynthesis and CO 2 emission from plant and soil respiration.Although wetland ecosystems exhibit CO 2 uptake in the long-term, fluctuations in photosynthesis and respiration can cause interannual variability in NEE over a considerable range from uptake to release [6].Instantaneous values of photosynthesis are generally related to photosynthetically available radiation (PAR), while respiration is related to temperature (soil and air temperature) [7,8].Spatial variations in annual NEE values have been explained by differences in vegetation, groundwater level, and plant biomass [7,9].CH 4 fluxes from wetlands have been measured mainly in northern peatlands and reported to have large spatio-temporal variability.Temporal variations of CH 4 fluxes are due to air temperature [10], soil temperature, groundwater level, and atmospheric pressure [11].Spatial variations in CH 4 fluxes have been reported to occur due to differences in groundwater level [12,13], plant biomass and vegetation characteristics [14], temperature and topography [15].Spatio-temporal variations in N 2 O have been reported to be due to lowering of the groundwater table [16][17][18][19], but there are few examples of measurements and the details are unknown.
Scaling up plot-level GHG flux measurements to regional and global scales is important for understanding the interaction between climate change and GHG emissions from wetlands.However, scaling up GHG fluxes from wetlands with large spatio-temporal variability is challenging.As a method to evaluate GHG emissions from wetlands over a wide area, classifications of wetland components (topography and vegetation) and estimation of area cover using satellite images have been used.A combination of land surface classification using Landsat satellite imagery and flux measurements has been used to estimate basin-level CO 2 and CH 4 budgets in the tundra of northeastern Europe in Russia [20] and landscape-level CH4 budgets in the wetland-upland complex ecosystem of Canada [21] estimates have been reported.
Wetland ecosystems are based on a delicate balance of environmental factors such as moisture conditions, and are considered to be sensitive to environmental changes.In addition to indirect effects such as climate change, wetland environments are also greatly affected by direct anthropogenic disturbances such as drainage for agricultural use.The anaerobic peat soils in wetlands rapidly become aerobic when drained.This results in increased CO 2 emission due to accelerated organic matter decomposition [22], decreased CH 4 emission [17,18], and depending on the type of wetland, increased N 2 O emission [17][18][19].
About 80% of Japan's peatlands are distributed in the Hokkaido region (about 200,000 ha), a northern part of Japan with a cool climate.It is thought that almost all this peatland was originally a wetland, but most of it has disappeared by now due to agricultural development since 1870.The area of Ishikari Mire Complex, the largest in Hokkaido, which used to be spread out in the Ishikari River basin, has decreased to 0.2% in 1983 (119 ha) compared to 1870 (50,690 ha) [23].Even in wetlands that have avoided disappearance due to agricultural development, changes in the wetland environment have occurred due to indirect effects of development.Typical examples are the influx of sediment into the wetland (floodplain mire) and following expansion of alder (Alnus japonica) forests due to human activities around the wetland [24,25], and the drying of the wetland and following invasion of Sasa (dwarf bamboo) [26][27][28].These changes in the wetland environment may alter its GHG budget.
In this study, to assess the impact of vegetation change triggered by drainage due to agricultural land development on the GHG budget of wetland ecosystems, we conducted a spatial estimation of GHG fluxes in a Sasa invaded wetland ecosystem.

Study Site
This study was conducted at the Bibai Wetland (43 • 19 N, 141 • 48 E), located in Ishikari peatland, central Hokkaido, northern part of Japan (Figure 1).In the surrounding area of Bibai Wetland, most of the wetlands have been disappeared due to agricultural reclamation since 1870.At present, about 22 ha of the Bibi wetland is preserved, but it is surrounded by agricultural land and related open ditches on all sides, and it has been pointed out that the wetland is drying up [29].In the central part of the wetland, where we conducted this survey, 3.5 to 4.5 m of peat accumulated on the mineral soil (Figure 2).In the peat profile, low moor peat, (minerotrophic fen type), transitional peat and high moor peat (ombrotrophic bog type) are deposited in order from the bottom [30].The high moor peat consisting of undecomposed residues of peat moss (Sphagnum spp.), sedges (Carex spp.), and vines (Vaccinium oxycoccus), is deposited at nearly 2 m.The peat at a depth of 10-30 cm has a high ash content and low T-C content due to the contamination of volcanic tephra that erupted nearly 300 years ago [30].The original vegetation of the wetland was mainly sphagnum moss (Sphagnum) and sedge plants, which covered the entire wetland in the early 1960s [29], but as the wetland became drier, dwarf bamboo (Sasa) invaded from the surrounding areas.The current vegetation consists of about 1 ha of Sphagnum papillosum, Moliniopsis japonica, Carex middendorffii, and Eriophorum vaginatum in the center of the wetland.Surrounding the Sphagnum-dominated vegetation, there is a vegetation of mainly Sasa palmata and Gale belgica var.tomentosa.Rhus trichocarpa is found in the more dried eastern part of the wetland, and Ilex crenata var.paludosa in the southwestern part of the wetland (Figure S1).The mean annual temperature at this region is 7.1 • C and the mean annual precipitation is 1156 mm.The snow season is usually from mid-November to mid-April.This study was conducted from January 2004 to November 2006.GHG flux measurement plots were set up from the area of Sphagnum-Moliniopsis vegetation in the central part of the wetland to the area dominated by Sasa and Ilex (Figure 1).Corresponding to the vegetation change, three to seven measurement plots were established in the Sphagnum-Moliniopsis site (hereafter "Sphagnum site"), the short-Sasa site (height of Sasa < 50 cm), the tall-Sasa site (50 cm < height of Sasa < 130 cm), and the Ilex site (from June 2005), respectively (Figure 1 and Figure S1).The measurement plots were changed in June 2005 within each vegetation site.GHG (CO 2 , CH 4 and N 2 O) fluxes from ecosystems were measured using a closed chamber method.In order to measure the emission and absorption of GHG by vegetation, a transparent acrylic chamber (30 × 30 cm (small chamber) or 60 × 60 cm (large chamber) at the bottom, 50-130 cm high depending on the vegetation) was used and left the vegetation inside the chamber during the measurement [31][32][33][34].For small chambers, chamber collars (30 × 30 cm at the bottom) equipped with a groove on the top were inserted into the peat surface to a depth of 5 cm before the beginning of this study.The groove was filled with water, and the chamber was inserted into the groove during each gas flux measurement.For large chambers, the transparent acryl chamber consisted of a rectangular main body (60 × 60 cm at the bottom, 50-130 cm in height) equipped with a groove on the top, and a detachable lid was used.The main bodies were inserted into the peat surface to a depth of 5 cm before the beginning of this study, and remained in place throughout the measurement period.As well as the small chamber, the groove of the large chamber was filled with water, and the lid was inserted into the groove during each gas flux measurement.Boardwalks were set up near the measurement plots to prevent disturbance during the measurement.
GHG fluxes were measured at 0 and 4 min after sealing the chamber for CO 2 (net ecosystem exchange of CO 2 : NEE), 0, 10 and 20 min for CH 4 and 0 and 20 min for N 2 O, respectively.Gas samples were collected in Tedlar bags (CO 2 ) and vacuum vials (CH 4 and N 2 O) and brought back to the laboratory for analysis.The frequency of surveys was 2-4 times per month during the growing season (non-snow season) and about once every month or two during the snow season.GHG fluxes from the snow surface were measured during the snow season.Simultaneously with the GHG flux measurements, soil temperature at a depth of 3 cm, water table depth, and photosynthetically active radiation (PAR) by a photon sensor (PAR-01, PREDE, Tokyo, Japan) were measured.
The concentration of CO 2 was analyzed using an infrared gas analyzer (ZFP5 or ZFP9, Fuji Electric, Tokyo, Japan), and CH 4 and N 2 O were analyzed using a gas chromatograph equipped with a flame ionization detector and an electron capture detector, respectively (GC-8A and GC-14B, respectively, Shimadzu, Kyoto, Japan).The GHG fluxes were calcu-lated by linear regression of the measurement time and the change in GHG concentration in the chamber.Positive and negative values of GHG fluxes indicate the net emission from the ecosystem to the atmosphere and the net uptake from the atmosphere to the ecosystem, respectively.

Total Ecosystem Respiration and Gross Photosynthesis
Net ecosystem exchange of CO 2 , i.e., CO 2 fluxes from the whole ecosystem (plant and soil), consisted of CO 2 emission by respiration of the whole ecosystem and CO 2 uptake by photosynthesis.Therefore, the relationship is described as Equation (1): where R TOT is the total ecosystem respiration, and P G is the gross photosynthetic CO 2 uptake.In terms of NEE, CO 2 uptake is a negative value and CO 2 emission to the atmosphere is a positive value.In the equation, both R TOT and P G are presented by absolute values.
After the measurement of NEE by the transparent chamber, the R TOT was immediately measured under light shielding condition (100% shading) with an opaque cover.P G was estimated by subtracting NEE by R TOT according to Equation (1).
To establish a PAR-P G curve, shading nets (60, 70, and 90% shading) were placed over the transparent chambers, and NEE measurements at low PAR level were also performed about once a month during the growing season.A photon sensor was installed at 120 cm above the ground in the center of the wetland and the PAR was continuously measured at 5-min intervals using a data logger (F-80, M.C.S., Sapporo, Japan).The soil temperature (3 cm) of each vegetation site was also measured continuously at 1-h intervals using a data logger.

Calculation of Cumulative GHG budgets
To calculate the cumulative NEE value during the growing season, estimation curve of T 3cm (soil temperature at a depth of 3 cm)-R TOT and PAR-P G were established at each measurement plots: where T 3 cm is the soil temperature at a depth of 3 cm, a and b are regression constants.
A PAR-P G curve is generally expressed by the following Equation (3) (e.g., [7]): where Q is the maximum P G under light-saturated conditions, K is the light intensity at which 50% of the maximum P G and a and b are regression constants.In order to account for seasonal changes in the curve, Q was assumed to vary linearly with the variation of the average temperature for two weeks before and after the measurement date (totally about one month), and was converted to "Q = c × (average temperature for one month), where c is a constant".This gives Equation (4) as follows: Using Equation ( 4), c and K were determined by non-linear regression (Kaleida graph, Synergy Softare).The c and K were calculated for each plot, and since the number of low PAR measurements was not sufficient for some of the short-Sasa and tall-Sasa plots, the average of constant K for other plots in the same vegetation site was used to create the equation.
Using the continuous values of soil temperature (3 cm) and PAR, continuous values of R TOT and P G at 1-h intervals were calculated from Equations ( 2) and (4), respectively.NEE was then calculated using Equation ( 1) and integrated over the growing season.NEE during the snow season and the cumulative values of CH 4 and N 2 O were calculated by a linear interpolation.The cumulative period of GHG budget was two years, from January to December in 2004 and from June 2005 to June 2006, and the averaged value was indicated as the annual GHG budget of each plot.The net GHG budget was calculated by multiplying the annual budgets of CO 2 , CH 4 , and N 2 O by the global warming potential (GWP, CO 2 :CH 4 :N 2 O = 1:34:298) [4], which takes radiative forcing into account.

Physico-Chemical Properties of Peat and Dissolved Organic Carbon (DOC) Concentration of Groundwater at the Flux Measurement Plots
In August 2006, bulk peat samples were collected from the surface layer of 0-10 cm at GHG flux measurement plots (20 plots in total).The samples were dried at 105 • C for 48 h to determine the moisture content.A part of the samples was air-dried and then ground for analysis.The degree of humification of the peat was determined according to Kaila's method [35,36], and the total carbon and nitrogen contents (T-C and T-N, respectively) were measured with an N/C analyzer (SUMIGRAPH NC-1000, Sumika Chemical Analysis Service, Osaka, Japan).
As with the peat samples, a stainless-steel tube was inserted into the peat at the flux measurement plots and groundwater was collected at a depth of 50 cm.The dissolved organic carbon (DOC) concentration of the groundwater was measured using a TOC Analyzer (SHIMADZU TOC-5000A, Shimadzu, Kyoto, Japan).

Observing the Spatial Distributions of Vegetation and Environmental Factors
Grid points were set at 50 × 50 m intervals on the southern half of the wetland (Figure 1), and latitude and longitude were measured using GPS.In July 2006, the vegetation coverage within a radius of about 15 m from the grid points was recorded and classified into five vegetation communities: Sphagnum-Moliniopsis, Moliniopsis-Gale, Sasa, Sasa-Gale, Ilex.Sasa height (8 replicates at each point) and water table depth were measured 5 times between September 2005 and September 2006 and averaged.Peat and groundwater were also collected and analyzed in August 2006 from the grid points as well as the gas flux measurement plots indicated in Section 2.5., except for peat C and N contents.

Satellite Imagery
The IKONOS image taken in September 2002 was used in this study.The ground resolution of IKONOS is 4 m.The four bands obtained by the satellite are Band 1, blue (445-516 nm); Band 2, green (506-595 nm); Band 3, red (632-698 nm); Band 4, near infrared (757-853 nm).The pixel values of the bands at each studied plot (grid points and GHG flux measurement plots) were obtained.Normalized vegetation index (NDVI) and ratio vegetation index (RVI) as vegetation indices indicating the amount and vigor of plants were obtained from the following Equations ( 5) and ( 6)

Estimating the Distribution of Vegetation
Using the vegetation classification surveyed at all plots as the objective variable and the satellite pixel values (Band 1-4, NDVI and RVI) as explanatory variables, discriminant analysis was conducted on all combinations of explanatory factors.Vegetation discrimination using Mahalanobis' distance was conducted for the combinations of factors with the highest discriminant validity (Band 1, 2, and 4).The accuracy rate (the ratio of the number of correctly discriminated samples to the number of samples used in the discriminant analysis) was calculated for each vegetation community and for the entire vegetation.Using the Sasa and Sasa-Gale communities as the Sasa site, the Sphagnum-Moliniopsis and Moliniopsis-Gale communities as the Sphagnum site, and the Ilex community as the Ilex site, the area of each vegetation site was calculated.

Estimating the Distributions of Environmental Factors
Multiple regression was performed with environmental factors as the objective variable and satellite pixel values (Band 1-4, NDVI and RVI) as the explanatory variables.The obtained regression equation was used to estimate the distribution of Sasa height and environmental factors in the entire wetland.The estimation of Sasa height was conducted in the Sasa site (Sasa and Sasa-Gale communities).The water table depth, moisture content of surface peat and DOC concentration of groundwater were estimated by multiple regression for the Sphagnum, Sasa and Ilex sites, respectively.Multiple regression was performed using 80% of the total data.Since a multiple regression was not able to be conducted for the Ilex site due to the small number of plots, a single regression was conducted and the factor with the highest regression coefficient was selected.

Verification of Estimation Accuracy
In order to verify the accuracy of the vegetation classification and Sasa height, a field survey was conducted in September 2006 at different plots from the grid points.The vegetation classification was verified at 26 sites, and the accuracy rate of the classification was calculated.Sasa height was verified at 12 sites, and the measured values were compared with the estimated values.For the water table depth, moisture content of the surface peat, and groundwater DOC, the remaining 20% of the data that were not used for multiple regression were used to compare the measured and estimated values.

Statistical Analysis
Statistical analyses were performed using Excel statistics 2012 for Windows (SSRI, Tokyo, Japan).Differences in the environmental factors, GHG budgets and bands among the vegetation sites were compared using a Tukey's HSD test (p < 0.05).

Spatio-Temporal Variation of GHG Fluxes
Seasonal variations of environmental factors and GHG fluxes during the measurement period are shown in Figure 3 and Figure S2.
The soil temperature (3 cm) in the tall-Sasa and Ilex sites were relatively lower compared to that in the Sphagnum site.During periods of low precipitation, the water table depth in the Ilex site was lower than in the other sites, with a maximum drop of nearly 50 cm.
The daytime NEE in the Sphagnum and Ilex sites showed an uptake of CO 2 throughout the growing season, and the uptake increased during the summer.Daytime NEE in the short-and tall-Sasa sites showed emission of CO 2 just after the snow melting, thereafter, tended to uptake CO 2 from the beginning of July.The NEE in the tall-Sasa site tended to show greater absorption of CO 2 than that in the short-Sasa site.The NEE in the Ilex site tended to be absorbed CO 2 throughout the growing season and remained higher absorption than that of other vegetation site.Large CO 2 uptake was found just before the snow season in the tall-Sasa and Ilex sites.The CH 4 flux in the Sphagnum site tended to show emission throughout the growing season, and tended to increase in the summer.Though the CH 4 fluxes in the short-and tall-Sasa sites remained at a lower level compared to that in the Sphagnum site, sometimes episodic high CH 4 emissions were found.The CH 4 fluxes in the Ilex site varied between slight emission and uptake.The N 2 O fluxes in the studied wetland varied between slight emission and uptake, except for slightly higher emission from the Ilex site.During the snow season, there was little emission or uptake of any GHG in all vegetation sites.R TOT were significantly positively correlated with soil temperature (3 cm) and water table depth in all vegetation sites (Table S1).The CH 4 fluxes were significantly positively correlated with soil temperature (3 cm) in the Sphagnum and short-Sasa sites.There was also a significant positive correlation between CH 4 flux and water table depth in the Sphagnum site.The temporal variation in CH 4 fluxes in the tall-Sasa and Ilex sites could not be explained by soil temperature (3 cm) and water table depth.N 2 O fluxes were not significantly correlated with soil temperature of (3 cm) and water table depth in all vegetation sites.

Estimating Cumulative NEE
There were significant exponential correlations between the soil temperature (3 cm) and R TOT , and significant regression expressions with R 2 = 0.673 to 0.941 (p < 0.05) were obtained in all plots (Table S2 and Figure S3).Based on the regression of the P G with PAR and the average monthly temperature, significant estimation curves for 11 of the 23 plots were established (Table S2).For the plots where estimation curves could not be established, the constant K was set to the average value of the same vegetation site, and the constant c was determined to establish the estimation curve.The correlation coefficients (R) of the estimation curves for P G in the Sphagnum site ranged from 0.373 to 0.877 (Table S2), and the curves corresponded relatively well to the measured values compared to other vegetation sites (Figure S4).High P G were sometimes measured at low PAR in the Ilex site, and some plots deviated from the estimation curves.

Annual GHG Budgets
The annual GHG budgets in each site are shown in Table 1.The R TOT in the Ilex site was significantly higher than those in the other three sites.The P G in the Sphagnum and short-Sasa sites were similar, and the P G in the tall-Sasa and Ilex sites were significantly higher than those in the other two sites.The annual NEE of the Sphagnum and short-Sasa sites were negative and positive (−49.7 ± 40.2 and 124 ± 42.9 g C m −2 y −1 , respectively), indicating net CO 2 uptake and emission, respectively.The annual NEE of the tall-Sasa and Ilex sites showed higher CO 2 uptake (−240 ± 79.5 and −990 ± 90.9 g C m −2 y −1 , respectively) than that of the Sphagnum site, and the annual NEE of the Ilex site was significantly different from those of the other three sites.The annual CH 4 emission was lower in the order of the Sphagnum > short-Sasa > tall-Sasa > Ilex sites, and was significantly lower in the tall-Sasa and Ilex sites than in the Sphagnum site.The annual N 2 O emission from the Ilex site was significantly higher than that in the other three sites.
The net GHG budget in each vegetation site was positive (emission) in the Sphagnum and short-Sasa sites, but negative (uptake) in the tall-Sasa and Ilex sites.The N 2 O emissions accounted for only a small percentage of net GHG budget in all sites (<1%).In absolute values, CO 2 and CH 4 accounted for 34.6 and 65.0% of the net GHG budget in the Sphagnum site, 68.1 and 31.6% in the short-Sasa site, 89.2 and 10.6% in the tall-Sasa site, and 99.3 and 0.142% in the Ilex site, respectively.The contribution of CO 2 to the net GHG budget was larger at the sites invaded by Sasa and Ilex.

Relationships between Annual GHG Budgets and Environmental Factors
The physico-chemical properties of surface peat and groundwater DOC concentrations at the flux measurement plots are shown in Table 2.The moisture content in the Ilex site was significantly lower than that in the other three sites.Kaila's degree of humification, which is used as an indicator of peat decomposition, was significantly lower in the Sphagnum site than in the other three sites.The T-C and T-N contents and C/N ratio of peat in the Ilex site tended to be lower than those in the other three sites.Although there was no significant difference, the groundwater DOC concentration in the Sphagnum site tended to be lower than those in the other three sites.The averaged water table depth was positively and significantly correlated with P G , R TOT and N 2 O, and negatively correlated with NEE and CH 4 (Table S3).Annual NEE in the Sasa sites (short-and tall-Sasa) showed a significant negative correlation (R = −0.743,p < 0.01) with Sasa height (Figure 4), indicating that net CO 2 uptake increased as Sasa grew.The annual CH 4 emission including all vegetation was significantly negatively correlated with groundwater DOC (R = −0.618,p < 0.01).The annual N 2 O emission was significantly negatively correlated with the moisture content of the surface peat (R = −0.608,p < 0.01).

Estimated Distributions of Vegetation and Environmental Factors Using Satellite Imagery
The water table depth of the Sasa and Ilex communities tended to be lower than that of the Sphagnum-Moliniopsis community in the measurements at the grid points on the southern half of the wetland (Table S4).The Sasa height and moisture content of surface peat were significantly correlated with water table depth (R = 0.552, p < 0.01 and R = −0.754,p < 0.01, respectively) (Figure 5 and Table S5).There was a positive, but not significant, correlation between groundwater DOC and water table depth (R = 0.231, p = 0.08).The satellite band values at the grid points and flux measurement plots for each vegetation community are shown in Table S6.The band values (1-4) of the Sasa community tended to be higher than those of the other communities.Discriminant analysis using the vegetation communities at the study site (Sphagnum-Moliniopsis, Moliniopsis-Gale, Sasa, Sasa-Gale, Ilex) as the objective variable and IKONOS bands 1, 2, 3, 4, NDVI and RVI as explanatory variables showed that the discriminant contributors were bands 1, 2, and 4 (Table S7).Accuracy rate of vegetation classification by discriminant analysis was lowest in the Sasa-Gale community (66.7%), but was 100% in the Sphagnum-Moliniopsis and Moliniopsis-Gale communities.Overall, a 93.6% of accuracy rate was obtained.Figure 6 shows the distribution of vegetation in the entire wetland estimated by the discriminant analysis.Multiple regression analysis was conducted by vegetation site using Sasa height, water table depth, moisture content of surface peat and groundwater DOC as objective variables and IKONOS band 1, 2, 3, 4, NDVI and RVI as explanatory variables and significant regression equations were established for each (Table S8).For Sasa height, all data were used for the analysis.For other environmental factors, 80% of all data were used to create regression equations, and the remaining 20% of data were used to verify the accuracy of the estimates (see below).Since multiple regressions could not be conducted for the Ilex site due to the small number of data, a single regression was conducted and those with high correlation coefficients were selected.Based on the above regression equation, Figure 7 shows the distributions of Sasa height in the Sasa sites and environmental factors in the entire wetland estimated from the satellite image bands.In order to verify the estimation accuracy of the satellite band, the vegetation at 26 sites was surveyed separately from the grid points.Although the accuracy rate for the Sphagnum-Moliniopsis community was low, the accuracy rate for the other communities was high, and the overall accuracy rate was 76.9% (Table S9).Bubier et al. (2005) [21] and Heikkinen et al. (2004) [20] reported a vegetation classification accuracy of 80% and 84%, respectively, which are comparable to the values obtained in this study.Environmental factors were relatively well estimated by satellite bands, but the estimated water table depth and groundwater DOC in the Sasa site had a narrower range than the measured values (Figure S5).

Estimation of GHG Budget from the Entire Wetland
The area of each vegetation site estimated using satellite bands was 3.86 ha for the Sphagnum, 0.122 ha for the short-Sasa (height < 50 cm), 17.3 ha for the tall-Sasa (50 cm < height < 150 cm), and 0.210 ha for the Ilex.The net GHG budget for the entire wetland (21.5 ha) was calculated using the GHG budgets for each vegetation site (Table 1), the regression equation between GHG budgets and environmental factors (Figure 5), and the distribution of vegetation and environmental factors estimated using the satellite bands (Figures 6 and 7).The regression equations of Sasa height-NEE (for Sasa sites), groundwater DOC-CH 4 , and the moisture content of surface peat -N 2 O were used to estimate the respective GHG budget.The NEE of the Sphagnum and Ilex sites were calculated by multiplying the annual NEE of each site by the area.
The NEE of CO 2 for the entire wetland was estimated to be −42.6Mg C y −1 , CH 4 to be 0.929 Mg C y −1 , and N 2 O to be 1.72 kg N y −1 (Table 3).The net GHG budget of the entire wetland was −113 Mg CO 2 -eq y −1 , and it served as a net GHG sink.The contribution of CO 2 uptake by the tall-Sasa site to the total GHG uptake of the wetland was 71% in absolute value, accounting for the majority of the total.The net ecosystem exchange of CO 2 in soil-plant ecosystems under natural vegetation has been reported to be −154 to −18.9 g C m −2 y −1 [37] and −105 to −53.0 g C m −2 y −1 [38] (i.e., uptake) in wetlands in western Finland.Heikkinen et al. (2002) [7] reported that NEE varied from −43 (uptake, intermediate fen) to 3 g C m −2 y −1 (emission, dry hummock) in the Russian tundra depending on the vegetation and water regime at the site.In the present study, the NEE in the Sphagnum site (−49.7 ± 40.2 g C m −2 y −1 ) was similar to those reported values, and the NEE in the tall-Sasa and Ilex sites (−240 ± 80 and −990 ± 91 g C m −2 y −1 , respectively) were lower (i.e., higher CO 2 uptake) than those reported values.
The NEE of CO 2 over the growing season in forest ecosystems has been reported to be about −350 g C m −2 in a Canadian poplar forest [39] and −122 g C m −2 in a Finnish birch forest (Leaf area index: LAI = 2.5) [40].The annual value of NEE was reported to be −525 g C m −2 y −1 for a deciduous mixed forest (LAI = 4.9) in Tennessee, USA [41].The NEE value of tall-Sasa site (LAI = 2.6) in this study was comparable to these reported values.CH 4 emissions have been reported ranging from 0.1 to 24.7 g C m −2 y −1 [37] and 13.1 to 36.0 g C m −2 y −1 [38] in wetlands of western Finland, and 2.0 to 16.0 g C m −2 y −1 in the Russian tundra [7].The CH 4 emission of 11.2 ± 1.9 g C m −2 y −1 from the Sphagnum site in this study was comparable to these reported values.Martikainen et al. (1993) [17] reported that N 2 O emissions were <4 mg N m −2 y −1 from natural wetlands and <4-143 mg N m −2 y −1 from drained wetlands in Finland.The N 2 O emission of 5.05 ± 4.02 mg N m −2 y −1 for the Sphagnum site in this study was comparable to that reported for the natural wetlands.Even the highest N 2 O emission in this study, in the Ilex site (47.2 ± 9.5 mg N m −2 y −1 ), was lower than the emission from the drained wetland reported previously.

Spatial Variation of GHG Budgets
Differences in cumulative NEE between sites have been reported to be due to groundwater level, differences in vegetation, and plant biomass [7,9].The cumulative NEE in this study also showed differences due to vegetation (Table 1).The invasion of shrubby Ilex into the original Sphagnum community in our study area increased CO 2 uptake.The NEE in the Sasa site was significantly negatively correlated with Sasa height (Figure 4), suggesting that the differences in NEE among the plots in the Sasa site were due to differences in plant biomass.Increasing Sasa biomass increases both the respiratory emission of CO 2 and photosynthetic uptake of CO 2 , but the effect of increasing photosynthesis was greater than that of ecosystem respiration, which is thought to have increased negative NEE (i.e., CO 2 uptake).
The emission pathway of CH 4 from wetland to the atmosphere consists of diffusion of dissolved CH 4 in groundwater along the concentration gradient, transport via aeration tissue (aerenchyma) of aquatic vascular plants, and ascent of CH 4 -containing air bubbles accumulated in peat (i.e., ebullition) [42].Diffusion may be the dominant pathway for CH 4 emission from peatlands dominated by Sphagnum [43,44].While it has been reported that most of the CH 4 transported by diffusion is oxidized by CH 4 oxidizing bacteria in the unflooded layer of the peat profile [44,45].Vascular plants, such as sedges, transport CH 4 in peat to the atmosphere through their aerenchyma without undergoing CH 4 oxidation in the aerobic layer, resulting in increased CH 4 fluxes [10,[46][47][48][49][50].Factors that have been reported to control the spatial variation of CH 4 emission from wetlands include groundwater level [12,13], plant biomass [14], and differences in temperature and topography [15].
In this study, the annual CH 4 emission showed a significant negative correlation with water table depth and a significant positive correlation with soil temperature (Table S3).The lowering of the water table depth, which was prominent in the tall-Sasa and Ilex sites, could cause the aerobic peat layer to thicken.As a result, CH 4 emission was considered to be low due to the lower CH 4 production and enhanced CH 4 oxidation.The high CH 4 emission in the Sphagnum site (Table 2), where the water table depth is slightly higher than in the short-Sasa site, suggests that the contribution of sedge-mediated CH 4 emission was significant.Possible reasons for decrease in CH 4 emission associated with Sasa invasion include reduced gas exchange at the peat surface [32] and the decrease in surface temperature due to ground cover by Sasa (Figure 3 and Figure S2) and following decrease in CH 4 production activity.
Although peatlands are considered to emit less N 2 O in their natural state (i.e., wetland), increased N 2 O emission due to drainage or lowering of the water table caused by climate change has been pointed out [5,17,18].In this study, the annual N 2 O emission was significantly positively correlated with the water table depth (Table S3), and N 2 O emission was slightly higher in the Ilex site where the water table depth was lower than in the Sphagnum site (Table 1).Maljanen et al. (2012) [51] reported that the N 2 O emission from boreal drained peatlands was greater when the C/N ratio of the surface peat was less than 20.They also reported that low soil pH, high nitrate availability and water table depth and addition of mineral soil were also associated with high N 2 O emissions.In this study site, the C/N ratio of the peat at depths of 40 to 150 cm was high (>50).In addition, the peat deeper than 40 cm was almost saturated with groundwater throughout the year (Figure 3), suggesting that N 2 O production in this zone was very limited.In contrast, the C/N ratio of the peat near the surface was lower than 20 (13.5-19.3 for 0-10 cm, Table 2), and the oxygen supply could be sufficient due to long duration without being saturated with groundwater.Therefore, the N 2 O production potential in the surface layer is considered to be high.The C/N ratio of the surface peat in the Ilex site with high N 2 O emission was the lowest among the sites (Table 2).The low T-C content, i.e., low organic matter content in the Ilex site indicates that the peat decomposition is more advanced than in the other sites.In the Ilex site, the surface peat dried out, and the supply of oxygen could accelerate the mineralization of nitrogen and following nitrification, resulting in the emission of N 2 O.
The annual N 2 O emission was also related to R TOT , P G and NEE (Table S3).This means that there was a positive correlation between plant biomass (mainly Sasa) and N 2 O emission.It has been reported that the rate of decomposition of organic matter and GHG release in peat soils is limited not only by the availability of oxygen, but also by the availability of mineral nutrients (especially N and P) for microbial growth [52,53].In this study, the lowering of the water table depth could accelerate the decomposition of peat and the mineralization and release of nutrients that had been stored in the organic matter pool [54,55], which could increase plant biomass and increase substrate and microbial activity followed by increasing N 2 O emission.
The study site is characterized by the presence of volcanic tephra mixed with surface peat and high mineral content (Figure 2 and Table 2).It has been suggested that the addition of mineral soil to drained peatland enhances N 2 O emission [51].In this study, no significant N 2 O emission was observed except for the Ilex site.However, in the future, when the peat is further dried at other sites and the environment becomes more favorable for N 2 O production and emission, the mineral mixed with the peat may promote N 2 O emission.

Sasa Invasion into the Sphagnum-Dominated Wetland
The vegetation in the studied wetland was determined by the water table depth, which was maintained at a high level in the Sphagnum community and lowered in the Sasa and Ilex communities (Figure 3 and Table S4).Sasa height was also correlated with water table depth (Table S4).In the studied field, these environmental changes are thought to have resulted in increased CO 2 uptake due to vegetation invasion and increased plant biomass, as well as changes in surface gas exchange such as decreased CH 4 emission and increased N 2 O emission.Fujimoto et al. (2006) [56] reported that the evapotranspiration rate in the Sphagnum-dominated community was higher than that in the Sasa-dominated community, and that 80-100% of rainfall was consumed by evapotranspiration in the Sphagnum-dominated community and 70-80% of rainfall in the Sasa-dominated community.In the Sarobetsu Mire in northern Hokkaido, where the invasion of Sasa into the Sphagnum vegetation is expanding similarly to the present study site, it has been reported that hydrological conditions such as distance from natural ditches and groundwater flow influenced by drainage and inflow conditions have a significant effect on the expansion of Sasa distribution [26,28].In addition to the above, considering the fact that the averaged water table depth is lower near the open ditches (Figure 7), the difference in water table depth for each vegetation in the study area is considered to be strongly influenced by the open ditches.
Several studies have shown that climate, hydrology, and nutrient status collectively control the broad distribution of vegetation in northern peatlands (e.g., [57][58][59]).Nordbakken (1996a, b) [60,61] investigated fine scale variation in vegetation species cover in relation to groundwater level in peatlands in Norway.While many of these studies have estimated species cover, recent studies on carbon and CH 4 exchange and nutrient cycling in peatlands require quantification of plant biomass, such as photosynthetic organs (e.g., [62]).Bubier et al. (2006) [63] investigated the relationship between species and environmental factors in order to improve carbon modeling studies, and reported that groundwater level was significantly positively correlated with leaf biomass and total biomass of vascular plants.If the relationship between groundwater level and vegetation biomass, soil environment, and GHG fluxes is common, as found in this study, then estimating broad-scale fluxes from satellite-based estimates of vegetation distribution and biomass may be an effective method in peatland ecosystems that are sensitive to environmental changes.

The Past, Present, and Future of the Bibai Wetland
As a result of the invasion of Sasa and Ilex into the original Sphagnum-diminated vegetation in the Bibai Wetland, the uptake of CO 2 increased, the emission of CH 4 decreased, and the emission of N 2 O slightly increased in the Ilex site (Table 1).Net GHG budget changed from emission in the Sphagnum site to uptake in the tall-Sasa and Ilex sites.Minkkinen et al. (2002) [64] predicted that draining and reforesting peatlands in Finland would increase CO 2 emission from the soil but greatly increase carbon stocks in tree biomass, reduce CH 4 emission and increase N 2 O emission to a lesser extent, resulting in reduced net GHG emission.The changes in the GHG budget with vegetation change in the wetland shown in this study are consistent with this report.
However, it has also been reported that net ecosystem production (NEP) decreases in older forests in the boreal and temperate zones [65].Therefore, there is a possibility that the current high CO 2 uptake in the tall-Sasa and Ilex sites will decrease in the future.In the present Bibai Wetland, CO 2 uptake has increased by 5 and 20 times, respectively, due to the invasion of Sasa and Ilex into the original Sphagnum-dominated community (Table 1), but at the same time progressing peat decomposition due to the lowering of the water table depth is also presented (Table 2 and Table S4).Furthermore, the NEE in the site was 124 ± 43 g C m −2 y −1 , suggesting net CO 2 emission.Therefore, it is expected that carbon emission from the peat was dominant in the process of invasion and establishment of short-Sasa into the original Sphagnum-diminated vegetation, leading to the present tall-Sasa.
Before 1960s, most of the wetland had been occupied by Sphagnum-dominated vegetation [29].The carbon accumulation rate of the entire wetland (21.5 ha) at the time estimated by current carbon exchange rate (CO 2 uptake by NEE-CH 4 emission) in the Sphagnum site was 8.26 Mg C y −1 .In contrast, the current carbon accumulation rate of the entire wetland was 41.7 Mg C y −1 (Table 3), which was about five times higher than the estimated rate for Shsagnum-dominated wetland due to the invasion of Sasa and Ilex, which increased carbon fixation.However, in this study, we did not separate peat decomposition (microbial respiration) from plant respiration, nor did we measure the supply of plant matter to the peat as litter.Therefore, the amount of peat decomposition due to the invasion of Sasa, and the contribution of carbon fixed by plants to new peat formation is not clear.A further understanding of the carbon cycle between the atmosphere, plants, and peat is needed.

Conclusions
The net GHG budget of the entire wetland where dwarf bamboo (Sasa) has invaded the original Sphagnum-dominated vegetation was estimated by combining the measurement of GHG fluxes by vegetation and the estimation of the distribution of environmental factors using satellite imagery.The net GHG budget for the entire wetland was estimated to be −113 Mg CO 2 -eq y −1 , acting as a net source of GHG.The CO 2 uptake by the tall-Sasa accounted for 71% of the total, indicating the importance of the accuracy enhancement on estimation for biomass and CO 2 balance of Sasa.Consequently, drawdown of water table followed by invasion of Sasa and Ilex into Sphagnum-dominated wetland increased carbon accumulation to the wetland ecosystem and changed it from a source to a sink for GHG at present.
The site studied here is an example that lowering of water table of soils with high organic matter content will not always lead to carbon loss from the ecosystem because the CO 2 uptake by developing vegetation can more than compensate the possible C loss from the soil.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10.3390/atmos12040448/s1, Figure S1: Vegetation of the study field (Bibai Wetland), Table S1: Correlation coefficients between GHG fluxes and environmental factors; Figure S2: Seasonal changes in air temperature, precipitation and snow depth (a: line, black and grey bars, respectively), soil temperature at 3cm depth (b), water table depth in a Sphagnum plot (c) and PAR (photosynthetically active radiation) (d), Table S2: Parameters for estimating total ecosystem respiration (R TOT ) and gross photosynthesis (P G ); Figure S3: Relationships between soil temperature at 3cm depth and total ecosystem respiration at each measurement plot, Table S3: Correlation coefficients between annual GHG budgets and environmental factors; Figure S4: Relationships between photosynthetically active radiation (PAR) and gross photosynthesis / monthly air temperature at each measurement plot, Table S4: Water table depth, physico-chemical properties of surface peat layer (0-10 cm) and ground water DOC (dissolved organic carbon) at grid points; Figure S5: Relationships between measured values of environmental factors and their estimated values using satellite bands, Table S5: Correlation coefficients among environmental factors at grid points; Table S6: Satellite image band values for each vegetation community; Table S7: Accuracy rate of vegetation classification by discriminant analysis; Table S8: Coefficients of multiple regression analysis of environmental factors using satellite image bands; Table S9: Verification of vegetation classification and Sasa height estimation using satellite bands.
Author Contributions: A.K., F.T., O.N. and R.H. conceived and designed the experiments; A.K., F.T. and O.N. performed the field measurements; A.K. analyzed the data; M.T. contributed analysis of satellite imagery; A.K. and F.T. wrote the paper; R.H. gave many constructive comments on this manuscript.All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Figure 1 .
Figure 1.Location of the study field (Bibai Wetland) and survey points.

Figure 2 .
Figure 2. Peat C and N content (a,b), C/N ratio (c), and ash content (d) by depth in the Bibai wetland.The high ash content and low T-C content of the peat at a depth of 10-30 cm is due to the contamination of volcanic tephra.

Figure 3 .
Figure 3. Seasonal changes in daily air temperature, precipitation and snow depth ((a) line, black and grey bars, respectively), soil temperature at 3 cm depth (b), water table depth (c), Net Ecosystem Exchange (NEE) (d) and CH 4 (e) and N 2 O (f) fluxes.Error bars indicate minimum and maximum value.Positive fluxes indicate net emissions to the atmosphere and negative fluxes indicate net uptake from the atmosphere.

Figure 4 .
Figure 4. Relationships between environmental factors and annual GHG budgets.** indicates significant at p < 0.01.

Figure 5 .
Figure 5. Relationships between averaged water table depth and Sasa height, ground water DOC or moisture content of surface peat at grid points.** indicates significant at p < 0.01.

Figure 6 .
Figure 6.Spatial distribution of vegetation communities in the Bibai Wetland estimated using a satellite image.

Figure 7 .
Figure 7. Spatial distributions of environmental factors in the Bibai Wetland estimated using a satellite image.

Table 1 .
Greenhouse gas (GHG) budgets at each vegetation site., net ecosystem exchange of CO 2 ; P G : gross photosynthesis; R TOT : total ecosystem respiration.Both R TOT and P G are presented by absolute values.For NEE, CH 4 and N 2 O, positive values indicate net emissions to the atmosphere and negative values indicate net uptake from the atmosphere.Numbers represent the average ± standard error.Numbers within a line followed by different letters differ significantly among the sites.(p < 0.05, Tukey-Kramer test). NEE

Table 2 .
Physico-chemical properties of surface peat layer (0-10 cm) and ground water dissolved organic carbon (DOC) at flux measurement plots.
Numbers represent the average ± standard error.Numbers within a line followed by different letters differ significantly among the sites.(p < 0.05, Tukey-Kramer test).

Table 3 .
GHG budget from the Bibai Wetland estimated by spatial distributions of environmental factors.