Characterizing Boreal Peatland Plant Composition and Species Diversity with Hyperspectral Remote Sensing

: Peatlands, which account for approximately 15% of land surface across the arctic and boreal regions of the globe, are experiencing a range of ecological impacts as a result of climate change. Factors that include altered hydrology resulting from drought and permafrost thaw, rising temperatures, and elevated levels of atmospheric carbon dioxide have been shown to cause plant community compositional changes. Shifts in plant composition a ﬀ ect the productivity, species diversity, and carbon cycling of peatlands. We used hyperspectral remote sensing to characterize the response of boreal peatland plant composition and species diversity to warming, hydrologic change, and elevated CO 2 . Hyperspectral remote sensing techniques o ﬀ er the ability to complete landscape-scale analyses of ecological responses to climate disturbance when paired with plot-level measurements that link ecosystem biophysical properties with spectral reﬂectance signatures. Working within two large ecosystem manipulation experiments, we examined climate controls on composition and diversity in two types of common boreal peatlands: a nutrient rich fen located at the Alaska Peatland Experiment (APEX) in central Alaska, and an ombrotrophic bog located in northern Minnesota at the Spruce and Peatland Responses Under Changing Environments (SPRUCE) experiment. We found a strong e ﬀ ect of plant functional cover on spectral reﬂectance characteristics. We also found a positive relationship between species diversity and spectral variation at the APEX ﬁeld site, which is consistent with other recently published ﬁndings. Based on the results of our ﬁeld study, we performed a supervised land cover classiﬁcation analysis on an aerial hyperspectral dataset to map peatland plant functional types (PFTs) across an area encompassing a range of di ﬀ erent plant communities. Our results underscore recent advances in the application of remote sensing measurements to ecological research, particularly in far northern ecosystems. creating the ﬁnal classiﬁcation analysis. These results support prior research in peatlands E.S.K. and Writing – original Writing – review J.R.R., E.S.K.,


Introduction
The far northern regions of the globe are warming at twice the global average, resulting in significant effects on the vegetation composition of arctic tundra and boreal ecosystems. Prominent among the changes to plant communities include the pan-arctic northward expansion of shrubs into tundra, and the spread of wetland vegetation as a result of permafrost thaw and ground subsidence [1][2][3][4]. Boreal forest ecosystems are also shifting from dominance by conifers to deciduous species as a result of increased fire prevalence [5]. The relationship between climate and plant community composition has important implications for biogeochemical cycles. For example, research indicates that elevated temperatures in arctic environments increase the rate of nitrogen mineralization, promoting the encroachment of woody plants into the tundra [6]. Increased shrub cover leads to an increase in the rate of carbon uptake, but with uncertain consequences for the overall arctic carbon cycle [4,7,8]. Shifts in arctic and boreal plant community composition and biogeochemical function can have significant impacts on wildlife [9,10] and the livelihoods of local communities [11,12]. Despite these wide-ranging impacts, resource and logistical challenges present an obstacle to the study of remote ecosystems. These challenges make the examination of factors such as plant community compositional change difficult. Indeed, given the extent and remoteness of northern ecosystems, remote sensing approaches to assessing vegetation change are necessary for understanding the ecological effects of climate change over large spatial scales.
Peatlands cover approximately 15% of land cover in the northern hemisphere and factor significantly in the global carbon cycle [13][14][15]. Peatlands are highly heterogeneous ecosystems characterized by thick organic soils and dominated by a variety of plant functional types (PFTs) depending on landscape position and hydrology [16,17]. PFTs refer to the physiological, morphological, structural, and chemical attributes that determine how plants allocate resources, take up carbon, and reproduce themselves [18,19]. In this study, we use the PFT framework as a trait-based method of vascular plant-species classification based on shared structure and function [18]. Changes in mean annual temperature and hydrologic conditions in peatlands can lead to shifts in the composition of PFTs, ultimately affecting ecosystem productivity [20][21][22][23]. Prior research in peatlands has indicated that warming and drying affect PFT composition by promoting dominance by woody shrubs [23][24][25]. Conversely, raised water tables from increased precipitation or permafrost thaw lead to an increase in graminoid species such as sedges [20][21][22]. Lowered water tables and increased shrub dominance have been associated with decreased annual gross primary productivity (GPP) relative to elevated water tables [22,26]. Studies have also demonstrated a link between the graminoid cover and methane emissions rates from wetland ecosystems [27,28]. Given their importance in regulating ecosystem function, mapping the distribution of PFTs could be an important element of predicting ecosystem responses to environmental change [29,30].
High spatial resolution hyperspectral imagery is emerging as a promising method for conducting remote assessments of PFT distribution [19,31]. Hyperspectral data leverages a far greater number of spectral bands as compared to multi-spectral data, ultimately allowing for the differentiation of vegetation characteristics with a high degree of fidelity [32]. For example, prior research has demonstrated that variation in the structural, physiological, and chemical characteristics of particular species or PFTs interact uniquely with different regions of the electromagnetic spectrum, allowing them to be distinguished from other species [33,34]. Species-specific properties that affect reflectance include unique leaf pigments such as chlorophyll and anthocyanin [35][36][37][38], leaf nitrogen content [39,40], and leaf water content [41]. Recent research has shown that phylogenetic distance both within and among species is correlated with spectral dissimilarity, a relationship that scales from leaf-level to plot-level canopy observations [42]  . Leveraging this type of information could allow for species to be mapped over large areas through classification analysis of hyperspectral data, improving our ability to track changes in species' distributions and diversity across broad spatial extents. Although mapping PFTs rather than individual species may not allow for accurate estimates of species diversity, it still affords valuable information regarding the structure and function of ecosystems Remote Sens. 2019, 11, 1685 3 of 22 as a whole [29]. The use of hyperspectral remote sensing as a tool for mapping PFTs could provide ecologists with the ability to track ecosystem responses to climate change [43]. Furthermore, these types of analyses could ultimately help inform management strategies and environmental policy approaches focused on ecosystem conservation and management in the face of climate change.
In addition to tracking compositional changes using remote sensing, a number of approaches to characterizing species and functional diversity using hyperspectral methods have also been investigated. In low diversity environments such as arctic and boreal ecosystems, the reduction of species richness may have amplified effects on ecosystem productivity due to low functional redundancy in plant types [44][45][46][47]. For example, peatland studies have found a loss of species in arctic systems in response to the warming effects of carbon cycling rates [45,48]. Thus, tracking changes in the vascular functional diversity of boreal peatlands is important for understanding the potential impacts of climate change on the function of these ecosystems [19]. Recent research using hyperspectral methods has linked spectral variation among pixels with species diversity by treating a stack of pixels as a reflectance signature that captures the optical characteristics of vegetation across hundreds of spectral bands [49,50]. Remote sensing estimation of diversity has previously been accomplished by linking spectral heterogeneity among pixels with species richness, a method based on the spectral variation hypothesis that unique species will present with strongly differentiated spectral properties [42,[51][52][53][54]. Several different statistical approaches have been used in the calculation of spectral heterogeneity. A simple method of estimating spectral heterogeneity from a multi-or hyperspectral dataset is comparing the coefficient of variation of multiple measurements [55][56][57][58]. Another approach is calculating Euclidian distances among different scans or pixels derived from multivariate eigenvector methods such as principal component analysis [52]. The further a pixel lands from the centroid of the data, the greater the spectral heterogeneity of that pixel [51,52]. This approach has been correlated with alpha diversity metrics that account for both species richness and evenness, such as the Shannon index, as well as with metrics such as functional diversity or the number of different plant types present in a pixel [42,59]. Together, these studies indicate the ability of multiple different remote sensing approaches to characterize the ecological diversity of a landscape. Ultimately, remotely sensed estimates of diversity are integral to mapping essential biodiversity variables (EBVs) -structural and functional characteristics of ecosystems that predict biodiversity levels at a global scale [60][61][62].
The objective of this research was to assess the ability of remote sensing approaches for characterizing and tracking changes in the plant communities of peatland ecosystems in response to global change drivers. We collected the vegetation and hyperspectral data at two different peatland ecosystems characteristic of boreal and arctic regions across the northern hemisphere. Both sites are located within large-scale ecosystem manipulation experiments examining the effects of a suite of climate change drivers on peatland ecosystems. We analyzed our data to determine which elements of vascular plant functional composition and structure influence the reflectance spectra. We also examined relationships between plot-level species diversity and spectral variation using multiple methodological approaches. Predicated on the relationships that we demonstrated between near-earth spectral reflectance and plant functional cover, we leveraged aerial hyperspectral imagery to map the distribution of PFTs across a boreal peatland in interior Alaska, USA.

Study Sites
We collected data at two sites representing typical northern boreal and arctic peatland ecosystems. The first site was at the Alaska Peatland Experiment (APEX) at the Bonanza Creek Long Term Ecological Research Station, located 30 miles west of Fairbanks, Alaska (64.82 • N, 147.87 • W) ( Figure 1). The site is located in a sedge-and shrub-dominated rich fen on the floodplain of the Tanana River. APEX was initiated in 2005 as a long-term manipulation to study the impact of altered hydrology on peatland ecology and biogeochemistry. APEX consists of three water table manipulation treatments: a raised, a lowered, and a control treatment, each 120 m 2 in size ( Figure 2, left panel). The lowered treatment plot is drained by a trench that borders the plot. Water from the trench is then pumped into the raised water table treatment. The water table in the control plot is not manipulated. We performed field data  collection within the APEX water table manipulation treatment

Vegetation Cover Sampling
Species and PFT cover were measured differently at the two study sites owing to differences in their respective data collection protocols (Table 1). In particular, non-vascular plants at APEX (which accounted for less than 10% of total cover) were recorded where present, whereas, at SPRUCE, mosses carpeted the soil surface and were therefore considered to have 100% cover across the site

Vegetation Cover Sampling
Species and PFT cover were measured differently at the two study sites owing to differences in their respective data collection protocols (Table 1). In particular, non-vascular plants at APEX (which accounted for less than 10% of total cover) were recorded where present, whereas, at SPRUCE, mosses carpeted the soil surface and were therefore considered to have 100% cover across the site The second study site was at the Spruce and Peatland Response Under Changing Environments (SPRUCE) project located at the Marcell Experimental Forest in northern Minnesota (47.30 • N, 93.29 • W) ( Figure 1). SPRUCE was initiated in 2014 as a climate change experiment to explore the effects of increased temperatures and elevated CO 2 levels on ecological and biogeochemical processes in peatlands at the southern edge of the boreal region ( Figure 2, right panel). The experiment is located within an ombrotophic bog dominated by Sphagnum mosses and black spruce (Picea mariana) with an understory of ericaceous shrubs, graminoids, and forbs. The SPRUCE experiment is a regression-based model in which ecological changes are measured in response to a broad range of temperatures [63]. Temperatures and CO 2 levels are manipulated in 10 chambered treatment plots ( Figure 2, right panel). An additional three unenclosed control plots were also tracked in this experiment. The chambers consist of octagonal enclosures 12.8 m in diameter and 7 m tall. The enclosures penetrate completely through the peat soils down to mineral material below. Peat soils average 2.5 m in depth across the site. Ambient temperatures and CO 2 levels are maintained in one control chamber while another control chamber has ambient temperatures and elevated CO 2 . The remaining eight treatment chambers are maintained at +2.25, +4.5 +6.75, and +9 • C above ambient temperatures to simulate a range of possible warmer climate scenarios. Carbon dioxide levels are maintained at ambient levels in half of the temperature treatment chambers and elevated to between 800 and 900 parts per million in the other half of the chambers. Each chamber, as well as the unchambered control enclosures, contain three 2-m 2 vegetation plots, where sampling was conducted for this experiment.

Vegetation Cover Sampling
Species and PFT cover were measured differently at the two study sites owing to differences in their respective data collection protocols (Table 1). In particular, non-vascular plants at APEX (which accounted for less than 10% of total cover) were recorded where present, whereas, at SPRUCE, mosses carpeted the soil surface and were therefore considered to have 100% cover across the site and were not included in the sampling protocol. At APEX, species cover was recorded using a point bar vegetation survey method [64]. Species were recorded along two perpendicular transects positioned across a series of replicate plots in each treatment. A metal frame was positioned above the canopy 1 m above the surface of the peat. A laser pointer was inserted into 10 evenly spaced holes in the metal frame along the transect. Species and height above the ground were recorded for all interceptions of the laser beam. This was repeated for each "hit" from the canopy to the surface of the peat. If the target was not photosynthetic it was recorded as either "standing dead", indicating that it was non-photosynthetic biomass but still standing, or as "litter", indicating it was part of the thatch forming at the surface of the peat. For the lowered and raised treatments, four replicate plots were sampled, and in the control treatment five replicates were sampled. Each species was categorized by PFT, and percent cover for each type was calculated based on the percentage of hits of each PFT relative to the entire sample (see Table 2 for species list). In addition, we calculated the leaf area index (LAI) for each point along each transect as the number of living layers of vegetation present between the canopy and the soil surface [25].  Vegetation cover at SPRUCE was sampled in mid-July across the range of temperature and carbon dioxide treatments. Three 2-m 2 plots were sampled in each of the 13 chambers for a total of 33 plots in the final analysis. Each plot was sampled using a 1-x 2-m frame that was set onto PVC pipes installed in the experiment to facilitate repeated data collection in the exact same area. The frame was divided into 50 cells 20 x 20 cm in size. The presence of all species was recorded within each cell of the grid, and species were later grouped into PFTs (Table 2). If a species or PFT was present in half of the cells it was considered to occupy 50% of the plot. Our sampling method at SPRUCE did not allow for the calculation of the LAI.

Spectral Sampling
To link spectral reflectance with functional composition and species diversity, sampling of spectral reflectance was performed at APEX and SPRUCE in the same plots where vegetation cover data were collected. Spectral reflectance measurements were taken at APEX using an Analytical Spectral Devices Fieldspec Pro that measured reflectance in 1-nm bandwidths between 300 and 2500 nm ( Figure 3). Data were collected during peak growing season on 29 June 2016. Scans were performed during a 1-h window on each side of solar noon, which occurred at approximately 12:30 p.m. local time. The sun at the time of data collection was at an azimuth angle of approximately 151 • east of north. Three scans were performed at each plot at the same location to acquire a local average. The instrument was positioned 1 m above the peat surface, and captured a field of view of 44 cm in diameter (Equation (1)) [65].
Data were processed such that water and atmospheric absorption bands that had a relative reflectance of greater than one or less than zero were automatically excluded, as well as areas that exhibited noise from the instrument. This resulted in three regions of the electromagnetic spectrum being included in the analysis ( Figure 3). Thirteen plots were sampled for vegetation and spectral reflectance among the three APEX water table manipulation treatments: five plots in the control treatment and four plots each of the lowered and raised plots. Remote Sens. 2018, 10, x FOR PEER REVIEW 8 of 25 Spectral data were collected at SPRUCE on 22 September 2016 under clear sky conditions beginning at 1:00 p.m. local time. The sun at the time of data collection was at an azimuth angle of 178° east. Although it was past peak growing season, senescence had only just begun at the site, and experimentally warmed plots were still near peak growing conditions. Data collection was performed using a PP Systems UniSpec-DC spectroradiometer that detected incoming and reflected solar radiation in 3-4 nm bandwidths between 310 and 1100 nm ( Figure 3). Three scans were performed above each vegetation plot in the same sampling location to acquire a local average. Details on data collection and processing are provided in McPartland et al. 2018 [25]. Reflectance was calculated for the average of the three scans following the methods of Harris et al. (2014) [16] and Wang et al. (2016Wang et al. ( & 2018 [49,50]. Bands at the beginning and end of the spectrum were eliminated due to noise, resulting in a data range of 400-1000 nm used in the analysis ( Figure 3). Data for SPRUCE chamber number 4 (+4.5 °C, elevated CO2) were also excluded because reflectance data for that plot presented as a significant outlier, with values that were far outside of the typical reflectance distribution collected at the other plots. This was likely due to changing sky conditions between the time that the white reference scan was taken, and the time of measurement.

Aerial Hyperspectral Data Collection
We used a hyperspectral dataset collected in July of 2014 at the Bonanza Creek Long Term Ecological Research Station, where the APEX field site is located. (Figure 2, right panel). The aerial dataset included the APEX study site and surrounding peatland ecosystems (Figure 4). The aerial image encompassed a variety of peatland cover types including black spruce and tamarack bogs, sedge meadows, tussock grasslands, and birch and willow shrub communities. Airborne hyperspectral data were collected by the remote sensing company SpecTir. The instrument used in the data collection was a ProSpecTir VNIR SWIR dual sensor push broom imaging spectrometer that collected data at 360 channels from approximately 400 to 2450 nm at a spectral resolution of approximately 2.9 nm in the visible and near-infrared, and 8.5 nm in the short-wave infrared portions of the electromagnetic spectrum. The imagery was collected approximately 200 m above ground surface, resulting in a spatial resolution of 1 m 2 [66]. The study area is approximately 2 km 2 in size. Data were processed using the Spectral Angle Mapper algorithm in the Python-based Spectral Classification Plug-in [67]. Ground control was performed by SpecTir and participants in this study. More detailed information on the geolocation and post-processing of the aerial hyperspectral data is provided in Anderson et al. (in review) [68]. Spectral data were collected at SPRUCE on 22 September 2016 under clear sky conditions beginning at 1:00 p.m. local time. The sun at the time of data collection was at an azimuth angle of 178 • east. Although it was past peak growing season, senescence had only just begun at the site, and experimentally warmed plots were still near peak growing conditions. Data collection was performed using a PP Systems UniSpec-DC spectroradiometer that detected incoming and reflected solar radiation in 3-4 nm bandwidths between 310 and 1100 nm ( Figure 3). Three scans were performed above each vegetation plot in the same sampling location to acquire a local average. Details on data collection and processing are provided in McPartland et al. 2018 [25]. Reflectance was calculated for the average of the three scans following the methods of Harris et al. (2014) [16] and Wang et al. (2016Wang et al. ( & 2018 [49,50]. Bands at the beginning and end of the spectrum were eliminated due to noise, resulting in a data range of 400-1000 nm used in the analysis (Figure 3). Data for SPRUCE chamber number 4 (+4.5 • C, elevated CO 2 ) were also excluded because reflectance data for that plot presented as a significant outlier, with values that were far outside of the typical reflectance distribution collected at the other plots. This was likely due to changing sky conditions between the time that the white reference scan was taken, and the time of measurement.

Aerial Hyperspectral Data Collection
We used a hyperspectral dataset collected in July of 2014 at the Bonanza Creek Long Term Ecological Research Station, where the APEX field site is located. (Figure 2, right panel). The aerial dataset included the APEX study site and surrounding peatland ecosystems (Figure 4). The aerial image encompassed a variety of peatland cover types including black spruce and tamarack bogs, sedge meadows, tussock grasslands, and birch and willow shrub communities. Airborne hyperspectral data were collected by the remote sensing company SpecTir. The instrument used in the data collection was a ProSpecTir VNIR SWIR dual sensor push broom imaging spectrometer that collected data at 360 channels from approximately 400 to 2450 nm at a spectral resolution of approximately 2.9 nm in the visible and near-infrared, and 8.5 nm in the short-wave infrared portions of the electromagnetic spectrum. The imagery was collected approximately 200 m above ground surface, resulting in a spatial resolution of 1 m 2 [66]. The study area is approximately 2 km 2 in size. Data were processed using the Spectral Angle Mapper algorithm in the Python-based Spectral Classification Plug-in [67]. Ground control was performed by SpecTir and participants in this study. More detailed information on the geolocation and post-processing of the aerial hyperspectral data is provided in Anderson et al. (in review) [68].

Analysis of Vegetation Composition and Species Diversity
For both sites, we calculated the percent cover of each PFT as a fraction of the total observations made at the plot level. The total number of plots included in our analysis was 13 across the three treatments at APEX, and 33 plots across 13 chambers at SPRUCE. At APEX, the effect of water table treatment on PFT cover was analyzed using a one-way ANOVA with a post-hoc pairwise means comparison to determine whether significant differences existed among treatments. At SPRUCE, we used a linear mixed effects modeling framework to evaluate the response of PFT cover to warming and elevated CO 2 . Temperature was added to the model as a continuous variable, and CO 2 level was included as a factor. We included treatment chamber in the model as a random effect and evaluated model fit using the Akaike information criterion, as well as both marginal and conditional R 2 , which describe the variance explained by the fixed effects and full model, respectively [69]. Linear mixed effects modeling was done using the lme4 package in R.
We calculated species diversity using the Shannon diversity index which incorporates both species richness and evenness into a single measure of alpha diversity: In this equation, Shannon diversity (H') is calculated at the plot-level in which p i represents the proportion of the population represented by species i [70]. Shannon diversity was calculated for each plot using the R Vegan package for community ecology [71]. We examined whether significant differences in species diversity among plots existed using a one-way ANOVA with a post-hoc pairwise means comparison at APEX, and linear mixed effects modeling at SPRUCE. We assessed significance based on an alpha level of 0.05 in the ANOVA, and evaluated model coefficients in the linear mixed effects model.

Spectral Data Analysis
For both sites, we performed principal component analyses (PCA) on the full reflectance spectra to determine which PFTs contributed most strongly to variation among spectral scans (Figures S1 and S2) [72,73]. This analysis reduced the dimensionality of the full dataset while preserving variation among plots within each experiment [74,75]. We were then able to determine the strength of the relationship between spectral variation and PFT cover. We performed the analyses using the R vegan package [71], and used permutation tests to determine the significance of different cover types on the spectral PCA results. This modeled the relationship between the PCA axes of the spectral data and PFT cover as a linear relationship. By this method, significance values are assigned by automatically comparing the sample data with a randomly permuted dataset. We considered different PFTs to be significant drivers of spectral variation at a p-value of less than or equal to 0.05.
We also examined the relationship between Shannon diversity and spectral variation through several measures of community and spectral heterogeneity following similar methods to those used by Wang et al. (2016) [49], who found a positive relationship between plot-level plant species diversity and the coefficient of variation of the associated spectral scans. Based on their findings, we predicted that greater coefficients of variation calculated for our average per-plot spectra would correlate with higher Shannon diversity index values [49]. We also performed a similar analysis based on research done by Rocchini et al. (2010) [53] that theorized that greater community heterogeneity would be correlated with greater spectral heterogeneity. In this analysis, we used measured heterogeneity as centroid distances derived from a principal component analysis where heterogeneity was measured as the distance from the center point of the PCA [53,76]. Basic linear regression was used to assess the relationship between diversity and spectral coefficients of variation, as well as between spectral heterogeneity and community heterogeneity measured using Euclidian distances from the global mean. Histograms were used to determine whether the data were normally distributed, and we used Anderson-Darling tests to reject the null hypothesis of uniformity in the distribution of our data ( Figure S3).

Hyperspectral Image Analysis
For the APEX field site, we used the aerial imagery to perform a supervised classification analysis using a random forests (RF) model to classify PFTs across the landscape [77]. We designated four different classes representative of dominant site PFTs. The classes were coniferous forests, low shrub, tussock meadow, and graminoid fen ( Figure S4). Before performing the analysis, we manually identified approximately one hundred training points within each class using a high-resolution aerial photograph collected in July of 2016 using an eBee senseFly unmanned aerial system (UAS) mounted with a Canon S110 camera (Figure 4, left panel). Although identifying training points in the field would have been ideal, we are confident that we were able to select training points that accurately reflected the ascribed vegetation cover class based on the image and our knowledge of the site. The UAS was flown at a height of approximately 92 m above the ground. The resulting imagery had a 3-cm ground resolution and was mosaicked together and georeferenced using a series of high-resolution GPS points that were collected via total station at the time of data collection. The UAS image was mosaicked using the Agisoft Photscan software platform. RF analyses were completed using the randomForest R package [78] and data were processed using the rgdal R package [79]. Standard model parameters were used in the RF model. Five hundred trees were created, and the number of variables considered at each node was the square root of the number of parameters. We used a model selection approach, as implemented in the rfUtilities R package [80], to find the most parsimonious and explanatory set of predictors. We used the overall out-of-bag accuracy and Cohen's kappa coefficient to evaluate model performance [81]. Variable importance was assessed using the mean decrease in accuracy and mean decrease in impurity upon iterated variable permutation, as implemented in the randomForest R package; variables that are more important will tend to cause higher mean decreases in both accuracy and node impurity when they permuted [77]. All associated mapmaking was performed using Quantum GIS [82].
Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 25 distributed, and we used Anderson-Darling tests to reject the null hypothesis of uniformity in the distribution of our data ( Figure S3).

Hyperspectral Image Analysis
For the APEX field site, we used the aerial imagery to perform a supervised classification analysis using a random forests (RF) model to classify PFTs across the landscape [77]. We designated four different classes representative of dominant site PFTs. The classes were coniferous forests, low shrub, tussock meadow, and graminoid fen ( Figure S4). Before performing the analysis, we manually identified approximately one hundred training points within each class using a high-resolution aerial photograph collected in July of 2016 using an eBee senseFly unmanned aerial system (UAS) mounted with a Canon S110 camera (Figure 4, left panel). Although identifying training points in the field would have been ideal, we are confident that we were able to select training points that accurately reflected the ascribed vegetation cover class based on the image and our knowledge of the site. The UAS was flown at a height of approximately 92 m above the ground. The resulting imagery had a 3cm ground resolution and was mosaicked together and georeferenced using a series of highresolution GPS points that were collected via total station at the time of data collection. The UAS image was mosaicked using the Agisoft Photscan software platform. RF analyses were completed using the randomForest R package [78] and data were processed using the rgdal R package [79]. Standard model parameters were used in the RF model. Five hundred trees were created, and the number of variables considered at each node was the square root of the number of parameters. We used a model selection approach, as implemented in the rfUtilities R package [80], to find the most parsimonious and explanatory set of predictors. We used the overall out-of-bag accuracy and Cohen's kappa coefficient to evaluate model performance [81]. Variable importance was assessed using the mean decrease in accuracy and mean decrease in impurity upon iterated variable permutation, as implemented in the randomForest R package; variables that are more important will tend to cause higher mean decreases in both accuracy and node impurity when they permuted [77]. All associated mapmaking was performed using Quantum GIS [82].

Response of Plant Functional Composition and Species Diversity to Experimental Maniuplation
Analysis of the vegetation cover data indicated that plant communities at both sites had diverged as a result of the experimental manipulation of water tables, temperatures, and CO 2 levels. At APEX, the most striking trend was the overall higher percent cover of litter and reduced cover of nearly all other PFTs in the lowered water table treatment as compared with the control and raised water table treatments (Table 3, Figure 5). Sedges, equisetum, forbs, and shrubs all showed trends towards higher abundance in the control treatment, and mosses were significantly more abundant in the control than other treatment (Table 3). At SPRUCE, we found a trend towards greater shrub cover with warming, coinciding with a decrease in forb cover ( Figure 6, Table 4). Graminoid species remained fairly consistent across the treatments ( Figure 6).
Our analysis of diversity also showed a trend towards greater alpha diversity in the control treatment as compared to the raised or lowered treatments at APEX (ANOVA Diversity F 2,10 = 3.52, p = 0.07) (Figure 7). We did not find any trends in plot-level diversity with warming or CO 2 level at the SPRUCE site (Table 3).

Response of Plant Functional Composition and Species Diversity to Experimental Maniuplation
Analysis of the vegetation cover data indicated that plant communities at both sites had diverged as a result of the experimental manipulation of water tables, temperatures, and CO2 levels. At APEX, the most striking trend was the overall higher percent cover of litter and reduced cover of nearly all other PFTs in the lowered water table treatment as compared with the control and raised water table treatments (Table 3, Figure 5). Sedges, equisetum, forbs, and shrubs all showed trends towards higher abundance in the control treatment, and mosses were significantly more abundant in the control than other treatment (Table 3). At SPRUCE, we found a trend towards greater shrub cover with warming, coinciding with a decrease in forb cover ( Figure 6, Table 4). Graminoid species remained fairly consistent across the treatments ( Figure 6).
Our analysis of diversity also showed a trend towards greater alpha diversity in the control treatment as compared to the raised or lowered treatments at APEX (ANOVADiversity F2,10 = 3.52, p = 0.07) (Figure 7). We did not find any trends in plot-level diversity with warming or CO2 level at the SPRUCE site (Table 3).    Figure 6. Percent PFT cover in response to warming at the SPRUCE project. Linear mixed effects models were fit to the PFT data with temperature as a continuous variable, CO2 level as a factor, and chamber included as random effect.  Figure 6. Percent PFT cover in response to warming at the SPRUCE project. Linear mixed effects models were fit to the PFT data with temperature as a continuous variable, CO 2 level as a factor, and chamber included as random effect.

Relationship between Community Composition and Spectral Response
The results of the principal component analyses (PCAs) of our plot-level spectral datasets indicated that variation in plot-level vegetation characteristics was related to differentiation in

Relationship between Community Composition and Spectral Response
The results of the principal component analyses (PCAs) of our plot-level spectral datasets indicated that variation in plot-level vegetation characteristics was related to differentiation in spectral reflectance at both sites (Figure 8). At APEX, results from our permutation tests indicated that spectral reflectance varied by treatment (r 2 = 0.45, p = 0.005), with the most variation in reflectance explained by plot-level species diversity (Table 5). Variation in the spectra at APEX was closely associated with the percent cover of mosses and forbs. Mosses and forbs together accounted for less than 10% on average of the cover at APEX, but had a strong effect on spectral reflectance, particularly in the control and raised water table treatment plots. The lowered water table treatment was associated with a greater litter cover, coinciding with lower plot-level species diversity (Figure 8, left panel). Overall, the first principal component explained 79% of the variation in the dataset, driven primarily by greater diversity in the control and raised water table treatments, and higher litter cover in the lowered water table treatment ( Figure 8, left panel, Table 5).
PCA results from the reflectance data at SPRUCE showed significant spectral differences across temperature (r 2 = 0.28, p = 0.01) and CO 2 treatments (r 2 = 0.19, p = 0.01) (Figure 8, right panel). Like APEX, reflectance spectra responded to differences in forb cover (Table 6). Unlike APEX, shrub cover also emerged as a significant vector in the SPRUCE data. Overall, the first principal component explained 96% of total variance among spectral scans. Since we did not find any strong differences in diversity among treatments, diversity did not emerge as a significant vector in the analysis. These results suggest that the relative cover of shrubs versus forbs associated with temperature and CO 2 level is the main driver of spectral differences among treatments at SPRUCE.

Relationship between Species Diversity and Spectral Variation
We found that plot-level species diversity was a significant predictor of spectral variation among plots at the APEX site, but not at the SPRUCE site (Tables 5 and 6). Based on this result we examined the relationship between species diversity and spectral diversity at the APEX site. The relationship between Shannon diversity and spectral CV showed a positive relationship (R 2 adj = 0.36, F 2,11 = 4.23, p = 0.06, RMSE = 0.30) (Figure 9, left panel). The relationship between community and spectral heterogeneity showed a positive relationship (R 2 adj = 0.21, F 2.11 = 4.11, p = 0.07, RMSE = 0.78) (Figure 9, right panel). It should be noted that differences in heterogeneity are relative to a global mean value. Points closest to zero along both axes are less heterogeneous than those further from the global mean.

Relationship between Species Diversity and Spectral Variation
We found that plot-level species diversity was a significant predictor of spectral variation among plots at the APEX site, but not at the SPRUCE site (Tables 5 and 6). Based on this result we examined the relationship between species diversity and spectral diversity at the APEX site. The relationship between Shannon diversity and spectral CV showed a positive relationship (R 2 adj = 0.36, F2,11 = 4.23, p = 0.06, RMSE = 0.30) (Figure 9, left panel). The relationship between community and spectral heterogeneity showed a positive relationship (R 2 adj = 0.21, F2.11 = 4.11, p = 0.07, RMSE = 0.78) (Figure 9, right panel). It should be noted that differences in heterogeneity are relative to a global mean value. Points closest to zero along both axes are less heterogeneous than those further from the global mean.

Hyperspectral Image Analysis-Mapping of PFTs
The results of the mapping project indicate that the random forests supervised classification analysis was highly successful at distinguishing among cover types ( Figure 10). The random forests model selection algorithm was set to separate cover into four classes based on the entire hyperspectral

Hyperspectral Image Analysis-Mapping of PFTs
The results of the mapping project indicate that the random forests supervised classification analysis was highly successful at distinguishing among cover types ( Figure 10). The random forests model selection algorithm was set to separate cover into four classes based on the entire hyperspectral dataset of 360 spectral bands. The overall model had an out-of-bag error rate of 9.85% and a Cohen's kappa coefficient of 0.869. Forest and shrub classes were classified with the greatest accuracy, and graminoid fen was the least accurate (Table 7). We found several spectral bands to be consistently important in making these classifications, as estimated by both the mean decrease in accuracy as well as the mean decrease in node impurity due to variable permutation (Table 8). Nearly all these bands fell within the 1800-2500 nm range representing the short-wave infrared region of the electromagnetic spectrum.

Hyperspectral Remote Sensing of Peatland Response to Climate Change
By leveraging two large ecosystem manipulation experiments we have demonstrated that hyperspectral remote sensing can capture the response of peatland functional composition and species diversity to hydrologic change and increasing temperatures. Prior research has already demonstrated that hydrology, temperature, and CO 2 affect peatland ecosystem structure and function [22,23,26,[83][84][85]. Remote sensing studies in peatlands have also demonstrated the efficacy of hyperspectral methods in mapping peatland flora [17,86,87]. However, we have not encountered any study that has used hyperspectral remote sensing to track peatland ecological response to the diverse suite of global change drivers explored here. Our field results from the APEX site indicate that long-term hydrologic change in rich fen peatlands leads to shifts in species diversity and ecosystem productivity. The dominant source of variation within the principal component analysis of the APEX spectral data appeared largely driven by the greater diversity and overall cover associated with the control and raised plots, and a higher percent cover of litter associated with the lowered plots, which describes a gradient of ecosystem GPP [47]. At the SPRUCE experiment, we found that species composition varied among treatments. The results of our vegetation cover sampling indicated an increasing trend in shrub cover, and an associated decreasing trend in forb cover with warming. These results align with previous research linking warmer air temperatures and increases in soil nitrogen availability with an increased shrub growth in far northern ecosystems [4,88,89]. In our analysis of the SPRUCE spectral data, shrub and forb cover both emerged as statistically significant vectors within the multi-variate analysis of reflectance. Our results support the feasibility of tracking changes in PFT cover, in particular of shrubs, using hyperspectral methods.
We did not track non-vascular moss cover at SPRUCE because they were ubiquitous across the site. However, prior research indicates that the moisture status of peat-forming mosses can strongly influence reflectance spectra [85,90]. Further research at SPRUCE would characterize the role of moss moisture status on hyperspectral reflectance properties. However, focusing solely on vascular plant cover, our results highlight the potential to use hyperspectral data to identify and track a range of different climate change-induced impacts to northern boreal peatlands.

Remote Sensing of Boreal Peatland Species Diversity
In low-diversity ecosystems such as peatlands, the loss of species diversity may have cascading effects on ecosystem function and productivity [44]. The low functional redundancy of peatlands is such that the loss of a species may also represent the loss of an entire plant type, thereby negatively effecting ecological productivity [47]. We demonstrate that drought in peatlands negatively affects diversity, particularly of mosses and herbaceous species. Specifically, variation and heterogeneity in spectral reflectance increased with increasing Shannon diversity among plots. Increasing spectral heterogeneity has been previously used as a proxy for diversity in studies that have mapped peatland ecological diversity using remote sensing [91][92][93]. This type of study is based on the species-area relationship in which species richness increases with scale across a heterogeneous landscape [94]. These studies have extrapolated species richness by calculating indices such as the normalized difference vegetation index (NDVI), then calculating heterogeneity for pixel aggregates of varying spatial scales [53,55]. Where these approaches have modeled diversity as a function of area, our approach captures diversity directly by comparing the spectral heterogeneity among plots with varying levels of diversity [53]. Alternative approaches have attempted to map diversity through the identification of unique species [54,[95][96][97]. However, this approach is challenging in non-forested systems in which individual organisms do not occupy an entire pixel. Our results indicate remote sensing may be used to estimate relative levels of diversity directly by comparing the spectral heterogeneity of measurements performed over a variety of different canopy types. We suggest that these results offer the ability to track diversity in systems with mixed or heterogeneous community assemblages that include trees, shrubs, graminiods, or other plant types [98].

Hyperspectral Characterization and Mapping of Plant Functional Types
Our results are in conversation with a growing body of research leveraging remotely sensed data to track vegetation change in peatland ecosystems [17,86,87]. The results of our analyses of both our plot-level and aerial datasets demonstrate the efficacy of using hyperspectral measurements to characterize the distribution of peatland PFTs. Through our close-range hyperspectral data collection we saw that the cover of PFTs was a driver of spectral variation. Despite relatively weak trends in our vegetation cover data at both sites, the reflectance spectra were highly responsive to subtle changes in the vegetation cover resulting from treatment. In particular, at SPRUCE, differences in the relative cover of shrubs versus forbs had a strong effect on spectral reflectance. This is likely due to the differences in stature, leaf structure, and foliar chemistry between high-growing woody shrubs compared with ground-layer herbaceous PFTs such as forbs [36,99]. We also predict that collecting data in the early fall improved our ability to detect variation among treatments because carotenoids and anthocyanins present in leaf tissues around senescence may have increased optical diversity, allowing for greater distinction among PFTs. Future research could examine the role of seasonality in improving the efficacy of hyperspectral remote sensing in making species-level distinctions.
Our plot-level findings prompted us to apply hyperspectral remote sensing to map peatland plant functional types over a larger spatial scale. The cover types represented within the aerial hyperspectral dataset included the PFTs that were present in our field measurements. We found the land-cover classification analysis was successful at distinguishing among different PFTs. This was particularly true for forest and shrub cover, which had the lowest error rate in the classification analysis. The percent cover of shrubs was also one of the strongest sources of variation in the field-based measurements at SPRUCE. We predict that the forested and shrub-covered regions of our study area were the most successfully mapped because their stature and leaf structure are clearly distinguishable from other PFTs, whereas graminiod fen and tussock grasses were less optically distinct from each other.
Our results indicate that reflectance in the short-wave infrared regions of the electromagnetic spectrum are particularly important in mapping peatland PFTs. Variable importance indicators such as model accuracy indicated that several SWIR bands between 1800 and 2500 nm were instrumental in creating the final classification analysis. These results support prior research in peatlands that has shown SWIR reflectance to be particularly sensitive to changes in peatland plant cover and moisture status [85]. Given the importance of changing shrub cover within the greater arctic and boreal landscape, these results support the application of hyperspectral measurements to track changes in the distribution of shrubs [4].
In arctic and boreal ecosystems where the environment is undergoing rapid change, the methods described here could prove invaluable to developing assessments of shifts in species composition over time. Due to the large spatial extent and low population density of the far north, examining changes to the landscape across broad spatial scales and over time must largely be achieved using remote sensing techniques. The use of hyperspectral imagery is integral to this approach because of the ability to select spectral bands that best capture variation across community types [90]. The high spatial resolution of aerial over satellite data is also key to this approach because it allows for fine-scale distinctions among cover classes. There are multiple research programs that collect hyperspectral data and make it publicly available to scientists. For example, NASA's Airborne Visible/Infrared Imaging Spectrometer collected aerial hyperspectral imagery over the APEX field site in 2017 and 2018 as part of the Arctic and Boreal Vulnerability Experiment (ABoVE) [100]. With the launch of NASA's Hyperspectral InfraRed Imager (HySPIRI) planned for 2021, global satellite hyperspectral imagery will become available to the remote sensing research community within the next several years [101]. The National Ecological Observatory Network also recently launched its hyperspectral Aerial Observation Platform with the mission of bridging spatial scales from organisms to landscapes [102]. A 2017 survey published by the National Academies of Science, Engineering and Medicine regarding the future of earth observation placed a high priority on hyperspectral data acquisition for the purposes of tracking diversity and changes in PFT cover (National Academies Press 2017). Our results establish statistical relationships between remotely sensed hyperspectral measurements, species diversity, and plant functional cover that could be used in future studies to identify and track community-and ecosystem-level changes in arctic and boreal regions.

Conclusions
The results from our study demonstrate that warming and hydrologic change lead to detectable changes in species composition and ecosystem productivity in boreal peatlands. We have further demonstrated the utility and efficacy of applying remotely sensed data toward characterizing and tracking changes in plant functional cover over time. Remote sensing is emerging as an effective tool to support ecological research, given its growing ability to characterize key elements of ecosystem structure and function such as leaf area, foliar chemistry, species diversity, and, more recently, the cover of different species and PFTs. We have shown that remote sensing can be used to detect the effects of a range of global change drivers on the productivity and functional composition of two types of peatland ecosystems representative of extensive areas of the northern boreal environment. We have also shown the feasibility of using hyperspectral remote sensing to map the distribution of PFTs over large spatial extents. In light of rapid global change, this approach could provide valuable insight into the changing ecology of the far north.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/11/14/1685/s1, Figure S1: Mean treatment reflectance spectra for the APEX site. Five plots are averaged in the control, and four are averaged in the lowered and raised plots. Each plot-level average is comprised of three scans taken in succession, Figure S2: Mean treatment reflectance spectra for the SPRUCE site. Each spectra represents six plots across two chambers. Each plot-level average is comprised of three scans taken in succession, Figure S3: Distribution of data included in the linear modeling of the relationship between hyperspectral and plot-level characteristics. We also performed Anderson-Darling tests to support our use of these data in our models. We were able to reject the null hypothesis that the data were uniformly distributed in all cases, Figure S4: Spectral response curves for the cover types represented in the aerial hyperspectral dataset to demonstrate differentiation among cover types.