Modeling and Mapping Forest Floor Distributions of Common Bryophytes Using a LiDAR-Derived Depth-to-Water Index ()
1. Introduction
Data available for mapping natural vegetation distributions at high-resolution have become increasingly accessible and important for environmental research, monitoring, and impact assessments [1] . Predictive vegetation mapping of bryophytes, however, is still largely unexplored [2] . Where bryophyte-environment relationships have been studied, the focus has either been placed on the micro-scale (e.g., [3] [4] [5] ) or on large regional to global scales (e.g., [6] [7] [8] [9] [10] ). This research differs in perspective by looking at bryophyte distributions along the forest floor at the landscape scale, at 1 m resolution. Earlier landscape studies were typically done at coarser scales (e.g., [2] [3] [11] ). The objectives of this article are:
1. To determine how the occurrences of specific moss species vary by forest floor and location conditions across select New Brunswick locations, based on 1-m2 plot surveys along upland/wetland transects, and by variations in soil wetness, slope, canopy closure, and litter depth. Other variables refer to aspect, ecosite classes, vegetation type (poor to rich), upland versus lowland, microsite topography (mound, flat, pit), soil type (organic, mineral), ecoregion and sampling location.
2. To transform the resulting presence/absence patterns into species-specific occurrence probability models by way of logistic regression analysis.
3. To validate and to apply the models so obtained for landscape-wide moss- occurrence probability mapping.
The means to determine how changes in soil moisture regimes would affect upland-to-wetland distributions of individual moss species was enabled through the increasing availability of high-resolution bare-earth digital elevation models (DEMs). Through DEM raster derivations involving slope, aspect, flow direction, flow accumulation, and flow network, two important soil moisture indicators emerged: (i) the terrain wetness index (TWI), and (ii) the cartographic depth-to-water index (DTW). The reported utilities of these indices favor DTW in terms of achieving greater conformance with field-measured variables such as soil drainage regime, soil type, various soil physical and chemical properties, soil trafficability, and vegetation type by soil moisture regime preference [19] [20] [21] . The difference between DTW and TWI pertains to the greater dependence of TWI on DEM resolution and smoothing [20] . Therefore, this article focuses on DTW for the purpose of forest floor bryophyte mapping at 1-m resolution.
2. Methods
Plot-based surveys were conducted across New Brunswick by selecting 12 areas for which 1-m resolution LiDAR-derived digital elevation models have become available (Figure 1; LiDAR: Light Detection and Ranging). These areas from ridge tops through valleys and wetlands, and were reasonably representative of the physiographic and climatic conditions across New Brunswick (Table 1). Additional criteria for location selection involved (i) accessibility, (ii) ownership (New Brunswick crownlands), (iii) low to no levels of industrial/residential developments, (iv) minimal human influence (e.g., not affected by dissecting roads), and (v) overall representation of New Brunswick ecoregions.
Sample plots were centered on forest floor and walkable wetland surfaces locations, as opposed to potential plots on large boulders and logs. Species composition and relative abundance of all ground flora were recorded for 980 1-m2 quadrats. Mosses were identified using hand lenses, and samples were brought back to the lab for identification under a microscope as needed. Relative abundance was recorded for each species based on percent area coverage when viewed from above [22] . Plot observations were done twice through independent viewing. Each plot was photographed for further validation purposes. The sampling locations were assigned to two groups: eleven for modelling moss presence and absence of specific moss species across New Brunswick (950 quadrats), and one for model validation (301 quadrats; University of New Brunswick Forest, Fredericton).
Figure 1. Moss survey locations across New Brunswick.
Table 1. Surveyed areas, by ecoregion, location, number of plots (n), climate conditions (1981-2010 Climate Normals, Environment Canada), and elevation.
Forest floor sampling plots were spaced >10 m apart at each location from upland to wetlands [23] , and were representatively chosen across 8 DTW classes ranging from very wet (DTW < 0.1 m) to dry (DTW > 12 m). Each plot determination involved specifying:
・ Plot location by longitude and latitude.
・ Moss species, by presence and abundance (all species); species were identified on-site using hand lenses, followed by microscopic indoor confirmation using a microscope, and Ireland and Hanes (1982) as authoritative species identification guide.
・ Leaf litter layer (L) depths.
・ Micro-topography classes (mound versus pit) [24] [25] .
・ Canopy closure.
・ Tree species composition by vegetation type (VT [26] ), with VT ranging from 1 (predominantly ericaceous species associated with poor soil conditions) to 4 (tolerant hardwood tree species associated with rich site conditions).
・ Ecosite category (eight classes: bog, fen, freshwater marsh, shrub wetland, forested wetland, ecotone, or riparian zone) [26] .
・ Wetland versus upslope location, with 1 denoting wetland, 0 otherwise.
These plot-specific assignments were further checked by cross-referencing with (i) GIS data layers pertaining to forest inventory and wetland cover and type, and (ii) aerial photographs [27] . LiDAR-derived 1-m resolution digital elevation models (DEMs) were used to generate the cartographic depth-to-water index (DTW) for each site (Figure 2). This index determines the difference in elevation between the ground surface and the nearest open-water features, such as flow channels and water pools. As such, it also emulates the gradation in soil drainage from very poor, poor, imperfect, moderate, well and excessively well [20] [28] . Table 2 provides a summary of all transect and plot determined numerical and categorical variables, by map-versus plot-based categories. The map-based variables were used for province-wide moss distribution mapping, whereas the plot-determined variables were intended to reveal finer resolution and therefore not-yet projectable moss-distribution differences.
All the data so assembled were entered into a single spreadsheet, which each column identified by plot, location, species and variable, and each row referring to plot-specific observations. Species absence/presence was also noted through binary coding (1 for present, and 0 for absent). The data so compiled were then used for (i) generating the probability response curves (Figure 3) for each species, (ii) multivariate logistic regression analyses and classification, and (iii) model-based presence-absence probability projections.
The probability response curves were calculated with the statistical programming environment R 3.1.2 [29] , using the eHOF package v. 1.5.7 [30] , with all model parameters (a, b, c, d) obtained through non-linear maximum likelihood estimation procedures [31] ). Selection of the most adequate model type was done using the Akaike Information Criterion (AIC [32] [33] ). For this analysis,
Figure 2. Panels show side-by-side surface image and depth-to-water map for each of the 12 survey locations across New Brunswick, also showing wetland outlines and the individual upland-to-wetland plot locations.
Table 2. Variables associated with each moss sampling plot.
Figure 3. Theoretical species absence/presence probability response models, ranked by increasing complexity: (I) no response, (II) sigmoidal, (III) sigmoidal with plateau, (IV) unimodal symmetric, and (V) unimodal skewed; direction of x is reversible (from [37] ).
species with less than 50 observations (n < 50) needed to be omitted [34] [35] [36] . This narrowed the analysis to seven species: Bazzania trilobata (BT), Dicranum polysetum (DP), Hylocomium splendens (HS), Polytrichum commune (PC), Pleurozium schreberi (PS), Sphagnum fuscum (SF), and Sphagnum girgensohnii (SG).
The binary presence/absence were also evaluated for each moss species in relation to (i) all the variables in Table 2 and (ii) the best-fitted HOF response models. This was done using Equation (1) and the logistic multivariate regression process (SPSS Statistics for Windows, Version 2.2.0; IBM Corp., 2013) as follows:
(1)
with p as the species-specific occurrence probability, α and βi as parameters, and xi (i = 1, 2, … k) as explanatory regression variables [38] [39] .
The presence/absence data of species with n > 50 was analysed by way of logistic regression, using all the variables in Table 2 as potential moss presence/absence predictor variables. The least significant predictor variables with α < 0.05 (Wald test [40] ) were eliminated, one step at a time. Within a set of closely correlated predictor variables, only the more significant variables were retained. The model coefficients so obtained were checked for possible collinearity using the variance inflation factor measure (VIF [41] ). Akaike’s Information Criterion (AIC) scores were used to select among the best-fitted models [32] , with the goodness of fit evaluated using the likelihood ratio test and two pseudo R2 measures: Cox and Snell R2 and Nagelkerke R2 [32] [39] [42] .
2.1. Model-Based Classification Including Validation.
Model Performance Receiver Operating Characteristic (ROC, SPSS) curves served as a measure of the classification accuracy of logistic regression models for the most common moss species, i.e., DP, PS, and SG. The area under the resulting ROC curves (AUC) determines the discriminatory ability of the model, and varies from 0.5 (no better than chance) to 1.0 (perfect discrimination [43] ). Models with AUC values >0.7 are considered to have high discriminatory power [38] . The classification accuracy was further evaluated using the proportional-by-chance accuracy rate, which is calculated by summing the squared proportion that each group (present/absent as observed) of the total sample across all sampling locations [44] . For model development and validation purpose, sampling locations were assigned to two groups: eleven for model development (980 quadrats), and one for model validation (325 quadrats; University of New Brunswick Forest, Fredericton). The logistic regression models generated from the former served to predict the plot-based moss occurrences for the latter. The resulting classification accuracy was evaluated as described above.
2.2. Model Projections and Mapping
The best-fitted regression models were used to map the probability of DP, PS, and SG presence across the landscape at each of the 12 sampling locations using ArcGIS mapping procedures [45] .
3. Results
Occurrence frequencies varied widely by moss species, from common to uncommon (Table 3). Climacium dendroides, Ptilium crista-castrensis, Rhytidiadelphus triquetrus, Sphagnum squarrosum, Aulacomnium palustre and Sphagnum wulfianum were locally abundant in terms of percent cover at only a few locations. In contrast, PS and SG were most common, appearing respectively in 328 and 301 of the 980 plots. DP was next, with 233 plot observations. SF, HS, BT, and PC occurred within 63 to 137 plots. The plots in Figure 4 show how the more abundant moss species were related plot-specific variations in DTW, canopy closure, slope and forest litter depth, with Sphagnum spp. occurring more frequently than the other species on wet, flat and open locations with little to no litter fall. The extent of percent moss cover per plot (i) by upland versus wetland and (ii) by mound versus pit locations on is also shown in Figure 4.
The HOF modelled presence/absence response curves to the continuous variables log10DTW, slope, % canopy closure and forest litter depth followednon-linear trends (Figure 5), with DP, HS, PC and PS displaying unimodal bell-shaped responses (HOF models IV and V). In contrast, SG and SF decreased monotonically (HOF models II and III) from wet to dry, as to be expected [47] [48] . The resulting best-fitted model parameters, optima and niche values are listed in Table 4, with R2 values ranging from 0.16 to 0.99. Only 5 of the 28 resulting models were not statistically significant with P > 0.1. Among these, HS had the lowest R2 values in response to log10DTW and forest litter depth. After excluding the HS response model to log10DTW, the best-fitted R2 and P values followed this sequence: DTW > Litter depth > Slope > Canopy closure.
Sources: [13] [24] [46] .
Figure 4. Mean values and their corresponding 95% confidence intervals for depth-to-water (a), percent slope (b), percent canopy closure (c), and leaf litter depth (d) associated with the occurrences of 7 moss species across the 11 sampling locations used for model development. Also shown: bar diagrams for comparing mean percent moss cover per plot of the same 7 species by upland versus wetland (top right), and by mound versus pit (bottom right) locations.
Figure 5. Best-fitted HOF presence probability response curves for the more frequent moss species of this study (n > 50), by DTW, slope, canopy closure, and forest litter depth.
The best-fitted logistic regression results for the species with n > 50 are compiled and summarized in Table 5 and Table 6. Of note are: (i) the consistently negative occurrence responses to (i) increasing litter depth, (ii) the positive versus negative effects of mound versus pit, upland versus wetland, and increasing versus decreasing log10DTW, and (iii) the varying effect of canopy closure. In one case (SG), slope (20 m focal mean) produced better results than log10DTW. Dropping one or more of the variables in Table 4 in favour of other somewhat correlated variables (e.g., log10DTW versus slope or wetland/upland location; canopy closure versus litter depth, SW or MX) increased the significance of some of those variables, but weakened the overall goodness of fit. SW or MX) increased the significance of some of these variables, but weakened the overall goodness of fit.
The best-fitted logistic regression results (logistic regression coefficients, significance levels and odds ratios) for DP, PS and SG using only map instead of plot-based predictor variables (i.e., log10 DTW, mapped classes for wetland/upland and forest type; Table 2) are listed in Table 7. Among these predictor variables, PS and DP increased while SG decreased with increasing log log10DTW. The wetland specification also discriminated against DP and PS. Without the wetland variable, log log10 DTW would be the strongest predictor
Table 4. Moss presence probability models (HOF models II to V, Figure 5) for the more frequent species of this study (n > 50): best-fitted parameters, goodness of fitness indicators (R2 and P values), with Model III maximum response position for DTW, slope, canopy closure, and forest litter depth.
variable for each of these species. All three species were additionally favored to occur under SW coverage, with MX also contributing to the presence of PS, only weakly so for DP, but not at all for SG.
Details about the goodness-of-fit of the logistic regression models in Table 7 are presented in Table 8. This table shows that the best-fitted models (i) are highly significant by way the likelihood ratio tests (p < 0.001), (ii) confirm that the results are not due to chance alone (AUC > 0.78), and (iii) account for the overall variations in moss absence and presence although the pseudo R2 indices remain low.
Table 5. Best-fitted logistic regression results for seven moss species (CC: Canopy Closure; SW: Softwoods; MX: Mixedwoods).
Also shown in Table 8 are the goodness-of-fit results for model validation, to indicate that the models derived from the 11 sampling locations across New Brunswick predict the DP, PS and SG occurrences for the validation location
Table 6. Summary of the best-fitted positive (+) and negative (−) trends in Table 5.
Table 7. Best-fitted binary multivariate logistic regression parameters including their statistical significance and odds ratios to predict the moss presence/absence probability responses as affected by log10DTW, wetland versus upland location, and forest type (SW, MX).
Table 8. Goodness-of-fit comparisons (Likelyhood Ratio Test, Cox & Snell and Nagelkerke R2 analogues, and AUC) for the logistic regression models in Table 7 for DP, PS and SG.
Table 9. Observed and predicted frequencies for bryophyte presence by logistic regression analysis (LR) across the 11 New Brunswick locations that were used for model development, and the Fredericton location used for model validation (Figure 1).
more or less equally well. The corresponding classification results are listed in Table 9 by number of cases correctly and incorrectly predicted. The accuracy ranged from 67% to 83%, with false positives and negatives amounting to 18% to 42%. The overall classification accuracy exceeded the probability by chance alone at 76.8%, with the significance level reached at 75.3%. The average level of probability by chance alone was 60.3%.
4. Discussion
4.1. DTW
The species response frequencies to the DTW-modelled variations in soil moisture regime are generally in agreement with the following expectations (Figure 4): hydric mosses such as SF are commonly found on wetlands and on poorly drained soils where DTW is low (i.e., < 25 cm). In comparison, SG spreads widely on fens, bogs and forested wetlands, but also grows in wet microsites (pits) on moderately well-drained upland forests [49] . In contrast, BT, DP, HS, PS and PC prefer moist to dry sites on uplands and on mounds (Figure 5). In combination, bryophytes often grow intermixed with each other, because of similar survival requirements [12] [50] . Hence, many of their response curves in Figure 5 overlap to varying degrees. For example, DP and PS typically occur together, with optimal DTW responses at 4 to 5 m. This also reflects their common occurrences on dry woodland towards the margins of swamps, and also on stumps [13] . While all of the species gravitated towards optimal DTW locations, HS and BT did only weakly so. This could be due to two reasons: (i) these species are found on suitable microhabitats across a wide range of substrate conditions [17] ; (ii) their occurrence frequencies, however, are low and widely spread across the DTW from 0.6 to 25 m.
4.2. Slope & Aspect
Since slope is in part correlated with DTW (slope = 1.57 DTW; R2 = 0.51), there is a general similarity between moss species occurrence with respect to changing slope and DTW, with SG gravitating towards flat locations, as observed elsewhere [51] , with Sphagnum mosses generally present only on the lowest and wettest locations [52] . The other species prefer upland locations, and here they are more frequent on mounds than in pits. Although aspect has the potential to affect microclimates and plant species distribution [53] , no such dependence was detected in this study, possible due to small sample size, and the preponderance of diffuse light under forest canopies
4.3. Ecosite
While there are clear differences in moss occurrences by uplands versus wetlands (Figure 5, top), this was much less the case when the occurrence frequencies were analysed by ecosite type. In part, this may be due to the sampling unevenness by ecosite type. Typically, regression results across categorical variables need to be well and evenly represented by sample size [40] . Bryophytes species, however, are not evenly distributed across varying wetland and upland conditions [24] [54] [55] . For this study, upland and wetland locations were essentially equally represented (n = 497 versus 483, respectively), but the wetland locations were split into nine categories, with n per category ranging from 14 (meadows) to 129 (fens).
4.4. Vegetation Type
The softwood stands surveyed generally had abundant bryophyte mats, thereby rendering SW to be the dominant vegetation type variable to predict the occurrence of DP, PS and SG. In general, SG is associated with wet coniferous forests dominated by black spruce (Picea mariana), tamarack (Larix laricina), and balsam fir (Abies balsamea). In contrast, PS and DP are mostly associated with upland forests dominated by red spruce (Picea rubens), eastern hemlock (Tsuga canadensis), balsam fir, and white pine (Pinus strobus). PS and DP were also correlated with the MX variable, which represents early to late successional mixedwood forests dominated by red maple (Acer rubrum), white birch (Betula papyrifera), and balsam fir, and later successional stages composed of yellow birch (Betula alleghaniensis), red spruce or eastern hemlock. Within these forests, ground flora is generally herbaceous and accompanied with bryophyte diversity [25] . In pure hardwood stands, ground flora is dominated by shrubs and ferns. Bryophytes, when present, are generally constrained to grow on exposed coarse woody debris, likely due to moss-suppressing leaf litter cover elsewhere [25] .
4.5. Leaf Litter Depth, and Mound Versus Pit Effects
Bryophytes growing on the forest floor have been shown to be negatively affected by deciduous leaf cover [17] [56] . In the present study, all species occurred less frequently with increasing litter depth above 0.5 - 1 cm. The effect was almost as significant as the moss occurrence relationships with DTW.
Some bryophyte species such as DP and PC appear to be more adapted to growing through broadleaf litter [11] . Of the species included in this study, none were frequently found where the depth of the leaf litter was greater than 6 cm thick, although some species fared better than others in the > 6 cm zone. In part, this could be due to their growth form. Studies have found that mosses with upright (acrocarpous) growth forms such as DP and PC respond better to litter burial than creeping mosses with prostrate (pleurocarpous) growth forms [56] . In addition, bryophyte community composition is closely related to microhabitat [57] , thereby leading to the mound versus pit related leaf litter distributions shown in Figure 5. In contrast, macro-habitat features related to macro changes in topography, aspect, and canopy closure tend to be less influential.
4.6. Canopy Closure
Forest canopy closure affects levels of light and precipitation reaching the forest floor, creating a complex microclimatic gradient, with differing levels of light, temperature, and humidity [5] . Canopy closure helps explain some variation in bryophyte occurrence, particularly for Sphagnum mosses, which were typically found in wetlands with little to no tree canopy, but remain constantly moist in full sun. In dry and fully exposed environments, bryophytes become metabolically inactive or desiccate altogether [58] . Aside from Sphagnum species, most bryophytes require at least small amounts of shade [59] . Heavy shading, however, reduces moss occurrences [5] [56] [59] [60] [61] [62] , with bryophyte mortality increasing with increasing aspen leaf coverage, mainly due to shading and lack of phytochemical responses [59] , with light/shade effects thought to become growth-limiting once all other resource needs are met [63] .
4.7. Multiple Regression Analyses
The best-fitted multivariate logistic regression trends in Table 5 show that litter depth, mound location, DTW, wetland versus upland type and coniferous versus deciduous forest type are the dominant positive or negative predictor variables. Generally, wetland and coniferous forest types have already been shown to influence the occurrence of bryophyte species across landscapes [2] [3] [11] . For the non-coniferous forest types, only mixedwood (MX) remained as a predictor variable for PS. Among the somewhat correlated terrain variables, DTW and wetland versus upland remained as significant predictor variables, while slope also became significant predictor variable for SG. Only three species were either positively (BT, HS) or negatively (DP) affected by canopy closure, while increasing litter depth had a negative effect on all 7 species.
In total, only up to 30% of the plot-by-plot moss occurrences could be quantified by way of the above procedures and analyses. The unexplained sampling variations would, for the most part, be due to local variations in substrate and vegetation mix conditions as these vary locally on account of stochastically varying processes pertaining to litterfall, shading, canopy throughfall, and vegetation growth and competition [64] [65] . The mix in resolution and the approximate nature of some of the predictor variables would also contribute to the unexplained variations, as follows:
・ The values for canopy closure, wetland and forest type are based on assuming uniform conditions across each sampling location. Hence, these data do not reflect the meter-by-meter differences in vegetation type and canopy shading across the forest floor.
・ In contrast, the DEM-derived DTW variable and the field determined litter depth and mound versus pit location are matching the 1 m2 sampling procedure quite closely. However, using slope as a predictor variable produced best results for SG but only so after 20 m focal smoothing.
・ Sampling was done unevenly across vegetation type, i.e., mostly but not exclusively on softwood sites. Similarly, bogs were more frequently sampled than fens, while steep riparian zones were less frequently sampled. Riparian moss species such as Climacium dendroides and Aulacomnium palustre were therefore under-sampled. More equally distributed samples could have been obtained by stratifying by more than one environmental variable, e.g. forest type or wetland type in addition to stratifying by DTW zones alone.
・ There is a 3 to 5 meter inaccuracy of determining the exact sampling location using the hand-held GPS device.
4.8. Model Validation
Testing the best-fitted presence/absence models for DP, PS, and Sphagnum moss against the independent data for the University of New Brunswick Forest (Fredericton) produced an overall classification accuracy of 67% to 83%. For DP and PS, absence was better predicted than presence. Sphagnum had a particularly high percentage of false negatives (40%, Table 9).
The above model results produced with the mappable predictor variables Table 2 are similar to the results produced with the plot-based predictor variables in Table 2 for all of the sampling locations in Table 1. This similarity is illustrated
Figure 6. Presence/absence observations (dots, top portion of each panel) and the corresponding modelled presence probability for each dot at St. Stephen (top) and at Dorn Ridge (bottom), all overlaid on the mapped probability of presence for Dichranum polysetum (left), Pleurozium schreberi (middle) and Sphagnum girgensohnii (right).
Figure 7. Presence/absence observations (dots) overlaid on the predicted presence probability map for Dichranum polysetum (top), Pleurozium schreberi (middle) and Sphagnum girgensohnii (bottom) for the sampling plots within the University of New Brunswick Forest in Fredericton (validation location).
in Figure 6 and Figure 7 for three of these locations: St. Stephen, Dorn Ridge and Fredericton. Generally, data-derived models tend to perform less well on data that are not part of model training. For example, [66] found that only a half of data-trained models performed reasonably well when tested against independent data.
5. Conclusions
The question of whether changes in bryophyte distribution could be modelled and mapped was explored using univariate HOF models and a series of logistic regression models, with depth-to-water, slope, forest canopy closure, and depth of the leaf litter layer as species-specific presence/absence probability predictors. These probabilities varied by species, as expected. Along the DTW gradient, Bazzania trilobata, Dicranum polysetum, Polytrichum commune, Hylocomium splendens, and Pleurozium schreberi had a greater probability of occurrence towards the drier end of the water table gradient (well-drained forested land), whereas Sphagnum fuscum and Sphagnum girgensohnii had the reverse trend. The depth-to-water index and slope worked best as a predictor for the Sphagnum species, due to their affinity for poorly drained, flat areas. Canopy closure appears to affect the presence probability of Bazzania trilobata, Dicranum polysetum, Hylocomium splendens, and Sphagnum girgensohnii. All species had considerably lower occurrences with increased depth of deciduous leaf litter.
The above results suggest that the presence/absence probabilities of bryophyte species can be predicted using a combination of environmental variables, with the depth-to-water index as a prominent predictor variable for at least some of the species. This research therefore extends some of the existing knowledge regarding common bryophyte responses to environmental gradients pertaining to variations in shading, exposure, moisture and nutrient availabilities, and substrate quality and conditions. While this work is limited to New Brunswick, the sampling effort reached across several upland/wetland combinations under maritime climate. This suggests that the approach taken and the models so produced could be useful for estimating the presence/absence probability of common moss species across New Brunswick and elsewhere under similar climate conditions.
Acknowledgements
This research was in part supported by the Environmental Trust Fund of New Brunswick and by the Natural Sciences and Engineering Research Council of Canada (NSERC: Discovery Grant, Canada Graduate Scholarship). Special thanks go to the Forest Watershed Research Centre at UNB for support in processing the LiDAR-derived wet areas maps (Jae Ogilvie, Doug Hiltz) and field sampling (Matthew Bolton, Benno Künzig, Friedrich Johan Wüthrich).