Landscape Patterns of Rare Vascular Plants in the Lower Athabasca Region of Alberta, Canada

Predicting habitat for rare species at landscape scales is a common goal of environmental monitoring, management, and conservation; however, the ability to meet that objective is often limited by the paucity of location records and availability of spatial predictors that effectively describe their habitat. To address this challenge, we used an adaptive, model-based iterative sampling design to direct four years of rare plant surveys within 0.25 ha plots across 602 sites in northeast Alberta, Canada. We used these location records to model and map rare plant habitats for the region using a suite of geospatial predictors including airborne light detection and ranging (LiDAR) vegetation structure metrics, land cover types, soil pH, and a terrain wetness model. Our results indicated that LiDAR-derived vegetation structural metrics and land cover were the most important individual factors, but all variables contributed to predicting the occurrence of rare plants. For LiDAR variables, rarity was negatively related to maximum canopy height, but positively related to canopy relief ratio. Rarity was therefore more likely in places with shorter canopy heights and greater structural complexity. This included fens, which overall had the highest rates of rare plant occurrence. Model-based allocation of sampling led to detections of uncommon species at nearly all sites, while the rarest species in the region were detected at an average encounter rate of 8%. Landscape predictions of rare plant habitat can improve our understanding of environmental limits in rarity, guide local management decisions and monitoring plans, and provide regional tools for assessing impacts from resource development.


Introduction
A key tool in biodiversity conservation is mapping where species are most likely to occur. This allows for the understanding of patterns in species occurrence, identification of biodiversity hotspots, and provides needed visualizations for managers. However, creating these products and understanding broad patterns in the spatial distribution of taxa require field data that can be difficult and costly to acquire. Furthermore, spatial models of plant richness, occurrence, or abundance require environmental predictors that have been collected over broad spatial extents at an adequate resolution [1]. Sampling and modeling of rare species (i.e., those of high conservation value) is complicated by the fact that they are often cryptic in nature, associated with uncommon habitats, or are largely data-deficient [2]. These factors contribute to knowledge gaps on rare species including information on where they are most likely to occur, how they may respond to disturbances, and where best to allocate limited conservation resources. adaptive, model-based sampling design that defines rare plant habitat for a suite of rare species in the Lower Athabasca region of northeast Alberta. Specifically, our objectives were to: (1) identify the landscape factors that most affect presence of rare vascular plants regionally; and (2) predict (map) rare plant habitat across the Lower Athabasca region. Understanding habitat relationships and the development of spatial tools and maps for the area provides an opportunity to improve land use decision-making, regional monitoring, and stewardship.

Background to Initial Model Building, Field Methods, and Data Sources
We adopted a model-based approach to surveying and predicting rare plant habitat within the Lower Athabasca region (south of Lake Athabasca) in northeast Alberta, Canada ( Figure 1) using the following three steps: (1) creating an initial model of potential rare plant habitat using existing data on rare plant species occurrences and readily available landscape predictors in the study region, in particular Ducks Unlimited Enhanced Wetland Classification (DU-EWC) land cover types [24] and terrain-derived measures from a 50 m digital elevation model (see Section 2.4 for more details); (2) stratifying field sampling based on model predictions and more heavily weighting visits to sites that were predicted to have rare plants, while still including predicted non-habitats to ensure widely distributed sampling and testing of the model; and (3) re-developing the rare plant model using new field data to improve predictions and the locations of subsequent sampling sites (i.e., we iterated this process for each year of sampling). Further details of each of these three steps are given below. The definitions of rare plants followed sub-national status ranks (S-ranks) assigned by the Alberta Conservation Information Management System (ACIMS), following the globally applied NatureServe methodology [25]. Specifically, we used any species ranked S1 and S2 as 'rare' and those ranked S3 as 'vulnerable' or uncommon.
Field sampling took place across a four-year period from 2012 to 2015, and in each year, steps 2 and 3, as described above, were repeated to improve sampling efficiency. A detailed field sampling methodology is provided elsewhere [26,27]. Briefly, plot size was 0.25 ha (50 × 50 m), which was selected to align with the plot sized used by a provincial-scale biodiversity monitoring program [28]. Experienced observers (i.e., individuals with multiple years of survey experience, n = 18 across all years) who underwent region-specific training were asked to survey the plot by walking parallel lines without time constraints. Surveys were completed during the growing season (mid-June to mid-August). Observers worked in pairs, but surveyed plots independently. Although the emphasis was rare plants, observers recorded the presence of all vascular plants within plots to fully describe assemblages. Given the large plot size, no effort was made to estimate cover or abundance of common species. Sites were constrained to be within 2.5 km of roads to ensure accessibility, but in one year of sampling helicopter support was used to access remote sites. Although the geographic extent of sampling was therefore limited in some areas, this did not appear to affect the environmental space of what was sampled, which was of primary concern here. In total, 602 sites were sampled over the four-year period at approximately similar sample sizes per year.
Site selection for field sampling in 2012 was weighted toward places that were predicted to have rare species based on species distributions predicted using maximum entropy models [20] for 66 regionally rare plants that used regional sources of existing occurrence data (~14,000 ACIMS records;~1100 records from other sources). Once initial field data were collected in 2012, we adjusted our approach to use only our field observations to ensure consistent protocols, a consistent plot size (0.25 ha), and recent data. We also focused on S1 and S2 ranked species, since nearly all sites encountered at least one S3 species. Rather than focus on individual species with limited numbers of observations, we instead coded any plots with at least one S1 or S2 species as 'present' (1) and those sites without a S1 or S2 species as 'absent' (0). This emphasized our objective of identifying potential rare plant habitat and not necessarily numerous models of individual species. Logistic regression [29] Forests 2020, 11, 699 4 of 15 was then used to model the probability of a site having a rare plant (S1 and/or S2) present using STATA 14 software [30]. Again, the DU-EWC layer was used in all years, with additional terrain variables included as predictors. A map of potential rare plant habitat was then estimated annually from the best supported model, determined using Akaike's information criteria (AIC) comparison [31], and sampling sites for subsequent years were directed by map predictions (described further in Section 2.2).
Forests 2020, 11, x FOR PEER REVIEW 4 of 16 best supported model, determined using Akaike's information criteria (AIC) comparison [31], and sampling sites for subsequent years were directed by map predictions (described further in Section 2.2).

Figure 1.
Location of Alberta in western Canada and the study area within North America (a); and terrain, elevation, major water features, and place names (towns) for the Lower Athabasca region (south of Lake Athabasca) (b).

Final Landscape Map of Rare Plant Habitat
The final model predicting rare plant habitat compared all field observations from four years of sampling (n = 602 plots) to the DU-EWC land cover classification (Figure 2a), a terrain-derived moisture index (2b), soil characteristics (soil pH) (2c), and LiDAR-derived vegetation structure metrics (2d). LiDAR-derived metrics represented vertical vegetation structure because initial tests indicated that these measures were stronger predictors of plant rarity than horizontal measures such as canopy cover. Climate variables were not considered since the study region is small relative to variations in climate.
The DU-EWC land cover classification included a number of detailed wetland classes (Table 1), for example, graminoid rich fen. Separation of wetland types (e.g., treed bog, graminoid poor fen) was needed given the prevalence and importance of these land cover types in the Lower Athabasca, where some areas are 60-70% wetland (primarily peatlands and much of it treed peatlands) [32]. Anthropogenic habitats (clear cuts, agriculture, industry developments) as well as aquatic habitats dominated by open water were not considered and thus masked in the final model by assuming absence of terrestrial rare plant species.

Final Landscape Map of Rare Plant Habitat
The final model predicting rare plant habitat compared all field observations from four years of sampling (n = 602 plots) to the DU-EWC land cover classification (Figure 2a), a terrain-derived moisture index (2b), soil characteristics (soil pH) (2c), and LiDAR-derived vegetation structure metrics (2d). LiDAR-derived metrics represented vertical vegetation structure because initial tests indicated that these measures were stronger predictors of plant rarity than horizontal measures such as canopy cover. Climate variables were not considered since the study region is small relative to variations in climate.
The DU-EWC land cover classification included a number of detailed wetland classes (Table 1), for example, graminoid rich fen. Separation of wetland types (e.g., treed bog, graminoid poor fen) was needed given the prevalence and importance of these land cover types in the Lower Athabasca, where some areas are 60-70% wetland (primarily peatlands and much of it treed peatlands) [32]. Classification: Protected A Terrain-derived variables included a moisture index estimated from a 50-m digital elevation model (DEM) [33] using the Compound Topographic Index (CTI) method [34,35]. Although a one meter LiDAR-based depth-to-water (DTW) metric [36] was available for part of the study region, comparisons within that area indicated that the 50 m DEM-derived CTI model performed as well or better than the more detailed DTW predictions, and thus CTI was subsequently used in all models. Soil characteristics of the region were described by soil pH obtained from Soil Landscapes of Canada version 3.2 [37]. Although other soil variables were available, they were either highly correlated with soil pH or did not correlate well with rare plant locations. LiDAR-derived metrics was available for most, but not all, parts of the study area (see gray areas in Figure 1d) effectively representing crown lands outside of the Cold Lake Air Weapons Range. Models using LiDAR data therefore represented a subset of plots (469 sites). We first developed models that covered the area of LiDAR extent (n = 469) and then applied the same model structure, minus LiDAR variables, for the full dataset (n = 602), thus allowing us to fill in those areas without LiDAR data with similar models. LiDAR point cloud metrics were summarized for the region at the same scale of the plot (50 m raster resolution) using FUSION software [38]. Final LiDAR-derived variables used for rare plant models included canopy relief ratio (CRR) [39,40] and maximum canopy height (95 th centile, shown in Figure 1b) as the other variables were highly correlated (r > |0.7|) and/or had little predictive value based on univariate tests. CRR is simply the mean height minus the minimum height divided by the maximum height minus the minimum height. Final predictor variables had a 50-m raster cell size to ensure that they corresponded to the plot size and were checked to ensure low between-variable correlations (r < |0.7|) to minimize problems with multicollinearity. Non-linear responses for continuous variables such as soil pH were considered and included if initial univariate tests demonstrated an improved model fit. Finally, AIC was used to compare combinations of variables as single-factor models (i.e., Soils [S], Terrain wetness [T], Land cover [L], and Vegetation structure [V]), two-factor models that had combinations of these single factors (note that single factor models could include more than one variable such as non-linear variables or multiple measures of the same factor such as vegetation structure), three-factor models that combined three of the four factors, and finally a single global model that included all four factors.  Table 1. Ducks Unlimited Enhanced Wetland Classification (DU-EWC) land cover types considered in models of rare plant habitat (occurrence) in the Lower Athabasca region of northeast Alberta, Canada (source: Ducks Unlimited). Note that some classes were combined in the analysis given their similarity and lower sample sizes. Number of plots (n) and percent coverage within study area by land cover class is provided. Terrain-derived variables included a moisture index estimated from a 50-m digital elevation model (DEM) [33] using the Compound Topographic Index (CTI) method [34,35]. Although a one meter LiDAR-based depth-to-water (DTW) metric [36] was available for part of the study region, comparisons within that area indicated that the 50 m DEM-derived CTI model performed as well or better than the more detailed DTW predictions, and thus CTI was subsequently used in all models. Soil characteristics of the region were described by soil pH obtained from Soil Landscapes of Canada version 3.2 [37]. Although other soil variables were available, they were either highly correlated with soil pH or did not correlate well with rare plant locations.
LiDAR-derived metrics was available for most, but not all, parts of the study area (see gray areas in Figure 1d) effectively representing crown lands outside of the Cold Lake Air Weapons Range. Models using LiDAR data therefore represented a subset of plots (469 sites). We first developed models that covered the area of LiDAR extent (n = 469) and then applied the same model structure, minus LiDAR variables, for the full dataset (n = 602), thus allowing us to fill in those areas without LiDAR data with similar models. LiDAR point cloud metrics were summarized for the region at the same scale of the plot (50 m raster resolution) using FUSION software [38]. Final LiDAR-derived variables used for rare plant models included canopy relief ratio (CRR) [39,40] and maximum canopy height (95th centile, shown in Figure 1b) as the other variables were highly correlated (r > |0.7|) and/or had little predictive value based on univariate tests. CRR is simply the mean height minus the minimum height divided by the maximum height minus the minimum height. Final predictor variables had a 50-m raster cell size to ensure that they corresponded to the plot size and were checked to ensure low between-variable correlations (r < |0.7|) to minimize problems with multicollinearity. Non-linear responses for continuous variables such as soil pH were considered and included if initial univariate tests demonstrated an improved model fit. Finally, AIC was used to compare combinations of variables as single-factor models (i.e., Soils [S], Terrain wetness [T], Land cover [L], and Vegetation structure [V]), two-factor models that had combinations of these single factors (note that single factor models could include more than one variable such as non-linear variables or multiple measures of the same factor such as vegetation structure), three-factor models that combined three of the four factors, and finally a single global model that included all four factors.
Models were fit for both datasets (with and without LiDAR). The lowest AIC was considered the most supported model. We then used the reported coefficients to map predictions using the ArcGIS (version 10.5, ESRI, 2016) raster calculator. Model performance and the predictive accuracy of final models were reported as percent deviance explained (pseudo-R 2 ) and Area-Under-the-Curve Receiver Operating Characteristic (AUC-ROC). Although ecological models often have poor explanatory power (2-5% R 2 , [41]), we considered models with pseudo-R 2 > 0.2 as being reasonably good in explanatory power and predictive accuracy assessed by AUC-ROC values based on the criteria of ROC < 0.7 being poor model accuracy, 0.7-0.9 good model accuracy, and > 0.9 high model accuracy [42,43]. As airborne LiDAR data describing vegetation structure were not available across the entire study area (Figure 1d), but we were interested in using the best available information (predictions), the final map of potential rare plant habitat used the LiDAR-based model where available, with the remaining areas filled in with predictions from the non-LiDAR model using the Conditional tool in ArcGIS (see the grey hillshade area in Figure 2d).

General Patterns in Plant Rarity
A total of 59 observations of 17 rare vascular plant species (i.e., those with S1 or S2 status, Appendix A, Table A1) were detected within 47 of the 602 plots (~8% of sites with at least one S1 or S2 species; 39 plots when limited to the extent of available LiDAR data) (Figure 3a). Vulnerable species (i.e., those ranked S3) were located at nearly 100% of plots. Encounter (detection) rate of rare plants (S1 and S2) varied substantially among land cover classes with the highest encounter rates found in  (Figure 3b), although some of these classes had lower sample sizes. In particular, graminoid and treed fens had high encounter rates of rare plants, at >10% of sites, more than 3-fold higher than that of upland deciduous and upland conifer (non-pine dominated) stands which averaged <3%. Indeed, the highest encounter rate of rare species (e.g., Carex heleonastes, Malaxis paludosa, and Juncus stygius) was found within graminoid poor fens at 23% of sites. Interestingly, although upland pine forest is considered to be species-poor in vascular plant composition, it contains a number of rare species (e.g., Spiranthes lacera, Carex umbellata, and Carex adusta) that reach presence in nearly 10% of survey sites.
Forests 2020, 11, x FOR PEER REVIEW 7 of 16 was found within graminoid poor fens at 23% of sites. Interestingly, although upland pine forest is considered to be species-poor in vascular plant composition, it contains a number of rare species (e.g., Spiranthes lacera, Carex umbellata, and Carex adusta) that reach presence in nearly 10% of survey sites.

Landscape Predictors of Plant Rarity
The most supported model (lowest AIC) predicting rare plant presence for S1 and S2 species included all individual landscape factors related to soils, terrain wetness, land cover, and vegetation structure (Table 2). When considering individual (single) factors, LiDAR-derived vertical vegetation structure, as measured by maximum canopy height (95 th centile) and canopy relief ratio, was more supported (∆AIC > 4) than soils, terrain wetness, or land cover alone. Land cover alone was ranked second, terrain wetness third, and last, soil pH that was similar to the null model suggesting no support when considered individually. Although land cover was less supported than vegetation structure, this could be due to classification errors in the land cover product or simply from greater model complexity since model accuracy (ROC) was similar between the two and fit was marginally higher for the more complicated land cover model ( Table 2). The most supported two-factor model included vegetation structure and land cover and this model was more supported than vegetation structure alone (single factor model). Interestingly, when considering the three combined landscape factors (more supported in AIC than the 2-factor model), soil pH was added to vegetation structure and land cover before that of terrain wetness, despite being initially neutral when alone. The final adjusted global model contained all four landscape factors, suggesting that the additional model

Landscape Predictors of Plant Rarity
The most supported model (lowest AIC) predicting rare plant presence for S1 and S2 species included all individual landscape factors related to soils, terrain wetness, land cover, and vegetation structure (Table 2). When considering individual (single) factors, LiDAR-derived vertical vegetation structure, as measured by maximum canopy height (95th centile) and canopy relief ratio, was more supported (∆AIC > 4) than soils, terrain wetness, or land cover alone. Land cover alone was ranked second, terrain wetness third, and last, soil pH that was similar to the null model suggesting no support when considered individually. Although land cover was less supported than vegetation structure, this could be due to classification errors in the land cover product or simply from greater model complexity since model accuracy (ROC) was similar between the two and fit was marginally higher for the more complicated land cover model ( Table 2). The most supported two-factor model included vegetation structure and land cover and this model was more supported than vegetation structure alone (single factor model). Interestingly, when considering the three combined landscape factors (more supported in AIC than the 2-factor model), soil pH was added to vegetation structure and land cover before that of terrain wetness, despite being initially neutral when alone. The final adjusted global model contained all four landscape factors, suggesting that the additional model complexity of including all four factors was important in predicting rare plant presence. Overall, the final (global) model had good model fit (pseudo-R 2 of 0.228) and model accuracy (ROC = 0.841; Table 2). Table 2. Comparison of candidate models describing the presence of rare plants (S1 or S2 rank) within 0.25 ha plots (n = 469) in the Lower Athabasca region based on soils (S), terrain (T), land cover (L), and vertical vegetation structure (V). AIC values in bold font represent the most-supported model (lower AIC) within single, two, three, and four-factor sets. Model complexity is represented by number of parameters (K). ROC represents predictive accuracy based on Receiver Operating Characteristic, while model fit (R 2 ) was measured by percent deviance explained. When examining model coefficients (β), the ranking of land cover type coefficients was similar to those illustrated in Figure 3b as raw encounter rates. For example, graminoid poor fens had the highest adjusted predicted presence of rare plants, being 205-times (odds ratio, exp(β)) more likely to encounter rare plants in this land cover type than the reference category (most common land cover type) of deciduous forests (Table 3). Non-linear responses of plant rarity to soil pH and terrain wetness (CTI) were supported with peak occurrence of rare plants at moderate levels of both variables. Finally, for vegetation structure, canopy relief ratio (CRR) was positively related to rare plant occurrences, whereas vegetation height (95th centile) was negatively related to rare plant occurrences (Table 3). Parameters included in the model without LiDAR metrics were similar to those with LiDAR variables (Table 3), while still maintaining moderately good model fit (pseudo-R 2 = 0.182) and similar overall model accuracy (ROC = 0.812). Predictions of S1 and S2 ranked vascular plant habitat are illustrated as patchy patterns of habitat throughout the region at both local and meso-scales, reflecting the importance of different land cover types, vegetation structure, soils, and terrain wetness (Figure 4). Some of the larger blocks of unsuitable habitat in the center and north of the study area appear to be due to soil limitations (soil pH), whereas localized variation of low, moderate, and high rare plant habitat appears to be more related to stand type and history (vegetation structure), land cover, and terrain wetness. Variation in hotspots of rare plant habitat throughout the area illustrate differential conservation values and threats across the region. Some notable sites with higher probability rare plant habitat occurrence include the southern parts of the Birch Mountains, the area around Winfred Lake east of Conklin, and the Marguerite River Wildland along the Saskatchewan border east of Fort McKay (see Figure 1 for place names). Table 3. Logistic regression coefficients (β), standard errors (SE), and significance (P) of the most-supported (AIC) model (with and without LiDAR data) describing the probability of a S1 or S2 vascular plant being present within 0.25 ha plots in the Lower Athabasca region of Alberta, Canada. Land cover variables are in comparison to the reference category of deciduous forest. Classification: Protected A Figure 4. Predicted distribution of rare vascular plant (S1 or S2 conservation rank) habitat within the Lower Athabasca of northeast Alberta, Canada. Prediction classes are based on model sensitivity (true positive rate), specificity (true negative rate), and optimal threshold classification probability where sensitivity and specificity are equalized. Both 'unlikely' and 'low' classes had model probabilities lower than the optimal threshold; the 'unlikely' class had a model specificity ≥ 0.9, while the 'low' class had model specificity < 0.9. In contrast, the 'moderate' and 'high' classes both had model probabilities higher than the optimal threshold probability, with the 'moderate' class having sensitivity < 0.9 and the 'high' class having a model sensitivity > 0.9 indicating higher confidence in a rare plant being present. Oil and gas lease boundaries, which represent areas wherein rights have been granted to companies by the crown for exploration and resource development, current to 2015, are shown. . Predicted distribution of rare vascular plant (S1 or S2 conservation rank) habitat within the Lower Athabasca of northeast Alberta, Canada. Prediction classes are based on model sensitivity (true positive rate), specificity (true negative rate), and optimal threshold classification probability where sensitivity and specificity are equalized. Both 'unlikely' and 'low' classes had model probabilities lower than the optimal threshold; the 'unlikely' class had a model specificity ≥0.9, while the 'low' class had model specificity <0.9. In contrast, the 'moderate' and 'high' classes both had model probabilities higher than the optimal threshold probability, with the 'moderate' class having sensitivity <0.9 and the 'high' class having a model sensitivity >0.9 indicating higher confidence in a rare plant being present. Oil and gas lease boundaries, which represent areas wherein rights have been granted to companies by the crown for exploration and resource development, current to 2015, are shown.

Discussion
Rare plant habitat in the Lower Athabasca region is best related to land cover type, soil pH, terrain, and vegetation structure metrics. When considering landscape predictors of rarity, we found that not only was the DU-EWC effective in predicting the occurrence of rare species, but so were LiDAR-derived vegetation structure metrics [23], particularly vegetation height (95th percentile) and canopy relief ratio (CRR). In fact, when considered individually, the LiDAR-derived vegetation metrics were similar to or better than land cover in predicting rare plant habitat. This suggests that remote-sensing based proxies of rare plant habitat may be used to not only predict current habitat, but also potentially used for monitoring change. The value of model-directed sampling for rare species was also clearly demonstrated in our study; we encountered "vulnerable" (S3) species at nearly 100% of plots and encountered rare species, here including "critically imperiled" (S1) and "imperiled" (S2) species at 3.6% and 5.3% of sites, respectively (see Appendix A, Table A1).
For the LiDAR metric CRR, higher values (i.e., more vertically complex vegetation with more emphasis of structure above mean height) correlated with increasing probability of rare plant occurrence. This may relate to forest age, where older stands may have more variation in their canopy height, or canopy patchiness may be reflective of variable stand history due to processes such as fire. Vegetation height had an inverse relationship with rare plant occurrence, where shorter canopies were more related to the occurrence of rare plants. This is logical given that fens were found to support more rare species occurrences, and these land-cover types have shorter canopy heights, with either stunted trees due to saturated growing conditions or shrubs and graminoids as the dominant vegetation layer. Rare plant occurrence was highest at intermediate values of both terrain wetness and pH. This is interesting given that some notable rare species in the region occur in 'extreme' habitats such as fens (e.g., Malaxis paludosa, Carex heleonastes, and Juncus stygius), where pH and wetness are high, or in jack pine forests (e.g., Spiranthes lacera and the Athabasca sand-plain endemic Lechea intermedia var. depauperata) where pH is high and wetness is very low, although many species included in modeling occurred at more mesic sites. Undoubtedly, models developed for individual species would demonstrate such idiosyncrasies in niche space.
Previous work that used remotely sensed data to predict rare plant habitat had similar predictive success to what was observed here, but more often focused on single species or very small groups [6,12], and very rarely on suites of rare species in a specific landcover type [6]. For example, de Queiroz et al. (2012) found that LANDSAT data in combination with LiDAR had stronger predictive ability for locating gyposphytes given the presence of gypsum mounds. Working also in Alberta, Nijland et al. (2015) combined a LiDARr-derived wetness index (DTW) and a Landsat-derived vegetation index to model site productivity [44]. In Brazil, da Silveira et al. (2018) combined Sentinel-2 and LiDAR metrics to better understand relationships between soils and vegetation [45]. In Florida, a LiDAR-derived digital terrain model was used in combination with other remote sensing products to map plant communities, with the management goal of monitoring change in a community strongly associated with an endangered plant species [46]. Our results further demonstrate the value of LiDAR, in combination with other variables representing important aspects of species niches, in both detecting new rare plant occurrences through field surveys and identifying important habitat associations. In addition, we note that rare plant species may respond at even more local scales than what was considered here (0.25 ha), as shown in previous examples that have linked rare plant occurrence to microtopographic features at scales <10 m [6,12]. Data obtained at an even finer resolution therefore could improve model predictions. However, at our scale of analysis, clear patterns between land cover types (e.g., fens, pine forests) and vegetation structure were observed and these relationships provide a basis for understanding regional patterns in rarity and conserving rare species.
Further to this, recent work has demonstrated that some rare plant species in our study area occur in species-poor sites like pine forests and poor fens, and thus approaches to conserving the most diverse communities will not satisfy the conservation of rare species and principles of complementary need to be considered [47]. Moreover, these habitats are scattered across the region, and mapping areas that are likely to contain rare plants shows a highly patchy distribution, particularly over areas that are leased to oil and gas companies for oil and natural gas extraction. This overlap is of concern given the known effects of oil and gas development on habitats, and in particular, the challenges associated with reclaiming peatlands [16,20], which were highlighted as important habitats for rare species in this study. This demonstrates the real-world applicability of regional mapping tools in species conservation and management; models employing spatial data such as LiDAR can guide management at regional to local scales and potentially aid in minimizing impacts to important rare plant habitats. Here, we show an overlap of areas where the probability of rare plant occurrence is high with oil and gas lease boundaries (i.e., areas where resource extraction and resulting land alternation are actively occurring). At the lease level, our rare plant model could be used to guide pre-development surveys, avoid areas of high rare plant habitat value, and plan for mitigation strategies.
Finally, in late 2015, as part of an effort to update the conservation status (ranks) of species in Canada, ACIMS published updated conservation status ranks for all species known to occur in the province based on observations and data collected since the time of their last assessment including data from this project. This resulted in a change in species status for some rare species in our study region. Given that our goal was to demonstrate the applicability and efficacy of employing remotely sensed variables in adaptive, model-based sampling, and that all field sampling and models were based on conservation status ranks published prior to 2015, we do not discuss how this change in status may change model predictions. However, this raises an important point that rarity status is in large part a product of our present understanding of species occurrence, distribution, and population trends, and thus are not static for any species over long-time frames, particularly for regions that lack significant sampling and thus knowledge gaps about rarity. Future applications of this method should consider the best available data of conservation status at the time of model building and sampling, as we have done here.

Conclusions
Although the majority of applications mapping rare plant habitat create landscape predictions for single species [3,4,6], we modeled the occurrence of any rare species (S1 & S1) within our study region, and demonstrated that environmental factors and LiDAR metrics can successfully predict the occurrence of a suite of species. Specifically, we found rare plant habitat to be associated with short and varied canopies, certain land cover types, especially fens, and areas with intermediate pH and terrain wetness. Although single species applications are valuable for modeling the realized niche of a species across space [3,12], multi-species approaches such as this allow for regional conservation planning and identification of landscape patterns in rarity by highlighting habitats that are most likely to house rare species [3]. We consider the application of an adaptive, model-based sampling method to be an efficient and effective approach to identifying rare elements and the availability and use of LiDAR products provide a unique opportunity to monitor rare plant habitat. Our success in detecting rare and vulnerable species using this approach may create a feedback loop wherein as more detections of rare species are made, we better understand niche associations, allowing for more accurate predictions, or conversely better understand that some species are less rare than originally perceived, particularly in poorly sampled landscapes. Therefore, use of this approach for sampling rare species and the use of its final products (landscape predictions of rarity) should be considered in regional environmental assessments and lease-level management, particularly since spatially explicit information on conservation values for plants is often lacking, especially when compared to wildlife.  Acknowledgments: Ducks Unlimited provided in-kind support by making available the Enhanced Wetland Classification of the Lower Athabasca layer.

Conflicts of Interest:
The authors declare no conflicts of interest. Project funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.
Appendix A Table A1. Rare (S1 and S2-ranked) vascular plant species encountered across 602 field plots (0.25 ha) in northeastern Alberta, Canada. Rare species were encountered at 54 sites. The number of plots where each species was detected is reported along with the species conservation status rank as of 2014.