Modelling Lichen Abundance for Woodland Caribou in a Fire-Driven Boreal Landscape

: Woodland caribou ( Rangifer tarandus caribou ) are reliant on Cladonia spp. ground lichens as a major component of their diet and lichen abundance could be an important indicator of habitat quality, particularly in winter. The boreal forest is typiﬁed by large, stand-replacing forest ﬁres that consume ground lichens, which take decades to recover. The large spatial extent of caribou ranges and the mosaic of lichen availability created by ﬁres make it challenging to track the abundance of ground lichens. Researchers have developed various techniques to map lichens across northern boreal and tundra landscapes, but it remains unclear which techniques are best suited for use in the continuous boreal forest, where many of the conﬂicts amongst caribou and human activities are most acute. In this study, we propose a two-stage regression modelling approach to map the abundance (biomass, kg / ha) of Cladonia spp. ground lichens in the boreal forest. Our study was conducted in Woodland Caribou Provincial Park, a wilderness-class protected area in northwestern Ontario, Canada. We used ﬁeld sampling to characterize lichen abundance in 109 upland forest stands across the local time-since-ﬁre continuum (2–119 years-since-ﬁre). We then used generalized linear models to relate lichen presence and lichen abundance to forest structure, topographic and remote sensing attributes. Model selection indicated ground lichens were best predicted by ecosite, time-since-ﬁre, and canopy closure. Lichen abundance was very low ( < 1000 kg / ha) across the time-since-ﬁre continuum in upland forest stands with dense tree cover. Conversely, lichen abundance increased steadily across the time-since-ﬁre continuum in upland forest stands with sparse tree cover, exceeding 3000 kg / ha in mature stands. We interpolated the best lichen presence and lichen abundance models to create spatial layers and combined them to generate a map that provides a reasonable estimation of lichen biomass ( R 2 = 0.39) for our study area. We encourage researchers and managers to use our method as a basic framework to map the abundance of ground lichens across ﬁre-prone, boreal caribou ranges. Mapping lichens will aid in the identiﬁcation of suitable habitat and can be used in planning to ensure habitat is maintained in adequate supply in areas with multiple land-use objectives. We also encourage the use of lichen abundance maps to investigate questions that improve our understanding of caribou ecology.


Introduction
The boreal ecotype of woodland caribou (Rangifer tarandus caribou) have evolved to occupy a niche unexploited by other northern ungulates [1]. Caribou tend to select low-productivity forests where potential high-quality habitats to target for protection or areas of overlap between humans and wildlife that present a high risk of conflict [23]. In conjunction with spatial predictions of predation risk, forage layers can be used to study the tradeoffs between nutrition and predator avoidance experienced by prey species [26,27].
In this study, we create a predictive model of lichen abundance in the boreal forest of Woodland Caribou Provincial Park, in Ontario, Canada and interpolate this model to create a spatial prediction (map) of lichen biomass. We use a regression modelling approach, first conducting field sampling within the study area to parameterize relationships between lichen abundance and environmental conditions (forest type, time-since-fire, canopy closure). We then relate lichen presence and lichen abundance to remote sensing and GIS data and use an a priori model selection procedure to identify the best explanatory variables. We interpolate the top lichen presence and abundance models across the study area and combine them to generate a map predicting lichen biomass (kg/ha). We show that our approach is straightforward and could be applied in other boreal caribou ranges with site-specific field data. Lichen maps could help managers develop more effective conservation strategies for woodland caribou. Managers could use lichen maps to track the availability of this important food resource over time and ensure a constant supply of lichen-rich habitat through resource or fire management planning. Paired with GPS collar locations, lichen maps could be used to identify the quantity of lichen in stands selected by caribou, aiding in the delineation of suitable habitat patches based on available forage resources.

Study Area
Our study area encompasses Woodland Caribou Provincial Park, a 5000 km 2 wilderness-class protected area in northwestern Ontario, Canada ( Figure 1) [28]. The park is a part of Pimachiowin Aki, a World Heritage Site that has received international recognition for its intact boreal forest and Indigenous cultural heritage [29]. The region is characteristic of the continuous boreal forest and is characterized by rolling terrain of bedrock outcroppings and numerous small lakes. Elevation varies from 309 m to 430 m above sea level and the park is situated on a plateau slightly elevated above the surrounding area, causing sparse conifer and dense conifer ecosites to compose a large proportion of the study area [30]. Sparse conifer (ecosite B012) occurs primarily on bedrock outcroppings where soils are very shallow (<15 cm) and moisture, nutrient availability, and plant diversity are low [12]. The overstory is dominated by jack pine (Pinus banksiana Lamb.) and the understory plant community consists primarily of Cladonia spp. ground lichens and velvet-leaf blueberry (Vaccinium myrtilloides Michx.). Dense conifer (ecosite B049) dominates upland sites with deeper, rocky soils (>15 cm) and nutrient and moisture conditions are more favourable for plant growth compared to sparse conifer [12]. A mixed overstory of black spruce (Picea mariana (Mill.) BSP) and jack pine characterizes dense conifer ecosites and the understory plant community consists primarily of feathermosses (e.g., Pleurozium schreberi (Brid.) Mitt.) and herbaceous plants (e.g., bunchberry, Cornus canadensis L.)). Small peatlands supporting black spruce and tamarack (Larix laricina (Du Roi) K. Koch) form in bedrock depressions and support an understory plant community dominated by Sphagnum spp. mosses and ericaceous shrubs (e.g., Labrador tea (Ledum groenlandicum Oeder)) [28]. There are no roads or resource development activities, historic or current, within Woodland Caribou Provincial Park. Development is limited to portage trails, campsites and several fly-in fishing camps. Large, frequent forest fires persist as an integral component of the local ecosystem due to a dry, continental climate [28]. The average annual area burned in the park over the last 30 years (1985-2015) is 0.6%, above the average for northern protected areas in Canada [31]. The study area is home to woodland caribou belonging to the Owl-Flinstone and Atikaki-Berens ranges in Manitoba and the Sydney and Berens ranges in Ontario [32].

Methods Overview
We combined field sampling with spatial environmental covariates to generate a map of lichen biomass for our study area ( Figure 2). First, we conducted vegetation surveys to quantify lichen cover and canopy closure in sparse conifer and dense conifer ecosites. We used conversion factors to estimate the stand-level lichen biomass (kg/ha) of each sampling location. Second, we derived nine environmental covariates from remote sensing and GIS data. Third, we assigned our field observations and environmental covariates to the GPS waypoint of each sampling location. We then used generalized linear models to predict lichen presence and lichen biomass as a function of a priori hypotheses built from our environmental covariates. We used model selection to identify the best candidate model and interpolated each top model to generate lichen presence and lichen biomass maps, which we combined to generate a final lichen biomass map for the study area. There are no roads or resource development activities, historic or current, within Woodland Caribou Provincial Park. Development is limited to portage trails, campsites and several fly-in fishing camps. Large, frequent forest fires persist as an integral component of the local ecosystem due to a dry, continental climate [28]. The average annual area burned in the park over the last 30 years (1985-2015) is 0.6%, above the average for northern protected areas in Canada [31]. The study area is home to woodland caribou belonging to the Owl-Flinstone and Atikaki-Berens ranges in Manitoba and the Sydney and Berens ranges in Ontario [32].

Methods Overview
We combined field sampling with spatial environmental covariates to generate a map of lichen biomass for our study area ( Figure 2). First, we conducted vegetation surveys to quantify lichen cover and canopy closure in sparse conifer and dense conifer ecosites. We used conversion factors to estimate the stand-level lichen biomass (kg/ha) of each sampling location. Second, we derived nine environmental covariates from remote sensing and GIS data. Third, we assigned our field observations and environmental covariates to the GPS waypoint of each sampling location. We then used generalized linear models to predict lichen presence and lichen biomass as a function of a priori hypotheses built from our environmental covariates. We used model selection to identify the best candidate model and interpolated each top model to generate lichen presence and lichen biomass maps, which we combined to generate a final lichen biomass map for the study area.

Field Sampling
To quantify lichen abundance, we conducted vegetation surveys at 109 sampling locations within and adjacent to Woodland Caribou Provincial Park from June-August 2018. We selected sampling locations based on time-since fire, stratified into decadal classes ( Figure 3; range = 2-119 years post-fire). We confirmed time-since-fire at sampling locations using an increment bore. Due to access constraints and the dominance of upland conifer in the study area, we constrained sampling to dense conifer and sparse conifer ecosites. Within each time-since-fire class we selected an equal number of sampling locations in each ecosite using a forest inventory map.

Field Sampling
To quantify lichen abundance, we conducted vegetation surveys at 109 sampling locations within and adjacent to Woodland Caribou Provincial Park from June-August 2018. We selected sampling locations based on time-since fire, stratified into decadal classes ( Figure 3; range = 2-119 years post-fire). We confirmed time-since-fire at sampling locations using an increment bore. Due to access constraints and the dominance of upland conifer in the study area, we constrained sampling to dense conifer and sparse conifer ecosites. Within each time-since-fire class we selected an equal number of sampling locations in each ecosite using a forest inventory map.

Field Sampling
To quantify lichen abundance, we conducted vegetation surveys at 109 sampling locations within and adjacent to Woodland Caribou Provincial Park from June-August 2018. We selected sampling locations based on time-since fire, stratified into decadal classes ( Figure 3; range = 2-119 years post-fire). We confirmed time-since-fire at sampling locations using an increment bore. Due to access constraints and the dominance of upland conifer in the study area, we constrained sampling to dense conifer and sparse conifer ecosites. Within each time-since-fire class we selected an equal number of sampling locations in each ecosite using a forest inventory map.   We accessed sampling locations by canoe and portage within the park and by truck in adjacent areas. At each sampling location, we established a start point 25 m from the edge of the mapped ecosite boundary and used a fiberglass tape to establish a 50 m transect oriented in a primary or secondary compass direction ( Figure 4). We placed a 1 m 2 quadrat at the 5 m, 15 m, 25 m, 35 m and 45 m marks of the transect to conduct five vegetation surveys per sampling location. We recorded the xy coordinates of each sampling location at the 25 m mark of the transect with a handheld GPS unit (accuracy ± 5 m). We spaced sampling locations a minimum of 100 m apart to reduce spatial autocorrelation. We accessed sampling locations by canoe and portage within the park and by truck in adjacent areas. At each sampling location, we established a start point 25 m from the edge of the mapped ecosite boundary and used a fiberglass tape to establish a 50 m transect oriented in a primary or secondary compass direction ( Figure 4). We placed a 1 m 2 quadrat at the 5 m, 15 m, 25 m, 35 m and 45 m marks of the transect to conduct five vegetation surveys per sampling location. We recorded the xy coordinates of each sampling location at the 25 m mark of the transect with a handheld GPS unit (accuracy ± 5 m). We spaced sampling locations a minimum of 100 m apart to reduce spatial autocorrelation. For each 1 m 2 quadrat, a single observer visually estimated the percent cover of each of the six most common Cladonia spp. ground lichens in the region (Table 1) and used a concave spherical densiometer to estimate the canopy closure above the quadrat. We recorded lichen cover for the sampling location by taking the average of the total lichen cover of each quadrat. Similarly, we recorded a single canopy closure value for each sampling location as the average value from the five quadrats. To derive estimates of lichen biomass, we multiplied the cm 2 area of the quadrat covered by each lichen species by its corresponding cover-to-biomass conversion factor (developed by  Table 1) [33]. We validated the conversion factors for use in our study area using destructive sampling (Appendix A). We estimated stand-level lichen biomass (kg/ha) for each For each 1 m 2 quadrat, a single observer visually estimated the percent cover of each of the six most common Cladonia spp. ground lichens in the region (Table 1) and used a concave spherical densiometer to estimate the canopy closure above the quadrat. We recorded lichen cover for the sampling location by taking the average of the total lichen cover of each quadrat. Similarly, we recorded a single canopy closure value for each sampling location as the average value from the five quadrats. To derive estimates of lichen biomass, we multiplied the cm 2 area of the quadrat covered by each lichen species by its corresponding cover-to-biomass conversion factor (developed by McMullin et al. (2011); Table 1) [33]. We validated the conversion factors for use in our study area using destructive sampling (Appendix A). We estimated stand-level lichen biomass (kg/ha) for each sampling location by adding the biomass estimates for each quadrat, converting from g to kg and multiplying by 2000 (see Appendix B for example calculation). We assigned the stand-level estimates of lichen cover, canopy closure and lichen biomass to the GPS waypoint of each sampling location for use in spatial modelling. Table 1. Cover-to-biomass (g/cm 2 ) conversion factors for the six most common Cladonia spp. ground lichens found in northwestern Ontario, Canada

Environmental Covariates
We selected nine environmental covariates (Table 2) supported by the literature to generate spatial models of lichen presence and lichen biomass [9,15,18,21,34,35]. The details of how the covariate layers were created are found in Appendix C. Note that the canopy closure layer was generated using a generalized linear model with forest inventory data, field measurements, time-since-fire and normalized difference vegetation index (NDVI; Appendix C). We converted all polygon datasets to rasters with 30 m pixels in ArcGIS 10.5 [36]. We resampled all covariate layers to have matching 30 m pixels and subsequently assigned values of each of covariate to the GPS waypoint of each sampling location using the mask() and extract() functions in the raster package in R version 3.6.0 [37,38]. This enabled us to subsequently relate lichen presence and biomass to forest structure, topographic and remote sensing attributes.

Spatial Modelling
We used our environmental covariates to generate a set of seven candidate models (Table 3) based on a priori hypotheses. Our base model includes ecosite, canopy closure and time-since-fire, which we anticipated would be the strongest predictors of lichen abundance. Each additional candidate model built on the base model by adding a topographic or remote sensing covariate. Covariates in the same candidate model have a Pearson's correlation coefficient < |0.6| to reduce collinearity within candidate models. We included a statistical null model (intercept) to assess the robustness of our candidate models. Table 3. Name and structure of candidate models used to predict lichen presence (0,1) and lichen abundance (biomass, kg/ha) as a function of forest structure, topographic and remote sensing attributes. Covariates within the same model have a Pearson's correlation coefficient < |0.6|. TSF = time-since-fire, Canopy = canopy closure, NDVI = normalized difference vegetation index, NDMI = normalized difference moisture index, SWIR2 = short-wave infrared reflectance (details in Appendix C).

Model Name Model Structure
Base Lichen~TSF + Ecosite + Canopy Elevation Lichen~TSF + Ecosite + Canopy + Elevation All Topography Lichen~TSF To generate a raster with cell values representing lichen biomass (kg/ha), we used our candidate models to conduct a two-stage modelling approach [24]: (1) lichen presence, (2) lichen abundance. We first used generalized linear models (family = binomial, link = logit) to identify the candidate model best explaining lichen presence (0 = absent, 1 = present). Lichen was considered present at sampling locations with >1% lichen cover (n = 87). We ranked competing models using Akaike's Information Criterion corrected for a small sample size (AIC c [45]) and considered the model with the lowest AIC c score as the top model. We interpolated this top model across the study area to create a raster with cell values representing probability of occurrence (0-1) for ground lichens. We used model-based interpolation as defined by Elith and Leathwick (2009) [46], implemented using the predict() function in the raster package of R version 3.6.0 [37,38]. We then created a binary layer where lichen is predicted to be absent (0) or present (1) in each pixel. We used the point on the receiver operator criterion (ROC) curve closest to the top left corner of the graph (0.71) as our presence threshold [47]. Lichen was classified as present (1) in cells with a probability of occurrence > 0.71 and absent (0) in cells with a probability of occurrence ≤0.71. We conducted k-fold cross-validation (k = 100; 60% training, 40% testing) to assess the accuracy of the lichen presence raster based on the mean area under the curve (AUC) statistic [48].
Once we generated the lichen presence raster, we used generalized linear models (family = Gamma, link = log) to identify the candidate model best explaining lichen biomass (kg/ha). We identified the top model as the candidate model with the lowest AIC c score and interpolated it across the study area to create a raster with pixel values representing lichen biomass (kg/ha). We multiplied this new layer by the lichen presence raster to create a layer that only predicts biomass in pixels where lichen is predicted to be present. We assessed the accuracy of this final lichen abundance raster by running a simple linear regression (R 2 ) between observed and predicted lichen biomass at each sampling location.

Lichen Biomass
Preliminary analysis of the field data revealed that post-fire lichen recovery differed markedly between sparse conifer and dense conifer ecosites ( Figure 5). Ground lichens were essentially absent from burns 0-19 years old in both ecosites and dense conifer supported low lichen abundance across the time sequence. Twenty years after fire, lichen biomass began to increase quickly in sparse conifer, reaching a median of 2648 kg/ha 40-49 years post-fire and leveling off thereafter. Mature sparse conifer ecosites supported approximately 2000-3700 kg/ha of ground lichens.
Preliminary analysis of the field data revealed that post-fire lichen recovery differed markedly between sparse conifer and dense conifer ecosites ( Figure 5). Ground lichens were essentially absent from burns 0-19 years old in both ecosites and dense conifer supported low lichen abundance across the time sequence. Twenty years after fire, lichen biomass began to increase quickly in sparse conifer, reaching a median of 2,648 kg/ha 40-49 years post-fire and leveling off thereafter. Mature sparse conifer ecosites supported approximately 2,000-3,700 kg/ha of ground lichens.

Spatial Modelling
The top candidate model for predicting lichen presence included ecosite, time-since-fire and canopy closure ( Table 4). The average AUC score from the k-fold cross-validation (k = 100) for the lichen presence model was 0.80, indicating 'good' model fit [48]. Table 4. Ranking of candidate models used to predict lichen presence as a function of forest structure, topographic and remote sensing attributes (Table 3). Models with a lower Akaike Information Criterion score (AICc) better describe the data. k = number of fixed effects (+1 for intercept) and wi = Akaike weight. SWIR2 = short-wave infrared reflectance, NDVI = normalized difference vegetation index, NDMI = normalized difference moisture index (Appendix C).

Spatial Modelling
The top candidate model for predicting lichen presence included ecosite, time-since-fire and canopy closure ( Table 4). The average AUC score from the k-fold cross-validation (k = 100) for the lichen presence model was 0.80, indicating 'good' model fit [48]. Table 4. Ranking of candidate models used to predict lichen presence as a function of forest structure, topographic and remote sensing attributes (Table 3). Models with a lower Akaike Information Criterion score (AIC c ) better describe the data. k = number of fixed effects (+1 for intercept) and w i = Akaike weight. SWIR2 = short-wave infrared reflectance, NDVI = normalized difference vegetation index, NDMI = normalized difference moisture index (Appendix C). Beta coefficients from the model describe the direction and magnitude of the effect of a covariate on the response variable. For example, in the top lichen presence model, probability of occurrence is positively associated with time-since-fire, increasing~1.6% per year since fire (Table 5). In the top model, lichen presence is negatively associated with dense conifer ecosites and there is a weak positive association between lichen presence and canopy closure (Table 5). The top candidate model for predicting lichen abundance was the same as lichen presence, including ecosite, time-since-fire and canopy closure ( Table 6). Table 6. Ranking of candidate models used to predict lichen abundance (biomass; kg/ha) as a function of forest structure, topographic and remote sensing attributes (Table 3). Models with a lower Akaike Information Criterion score (AICc) better describe the data. k = number of fixed effects (+ 1 for intercept) and w i = Akaike weight. SWIR2 = short-wave infrared reflectance, NDVI = normalized difference vegetation index, NDMI = normalized difference moisture index (Appendix C). The top lichen abundance model indicates lichen biomass is positively associated with time-since-fire, increasing 1.3% per year since fire (Table 7). Lichen biomass is negatively associated with dense conifer ecosites. There is a weak negative association between lichen biomass and canopy closure, with biomass decreasing by 0.4% per unit increase in canopy closure (Table 7).  Figure 6 displays the post-fire recovery of lichen biomass in sparse conifer and dense conifer ecosites as predicted by the top lichen abundance model. Note the shallow slope of the curve for dense conifer-lichen biomass is never predicted to exceed~1000 kg/ha. By comparison, lichen biomass increases quite steadily in sparse conifer ecosites, reaching~2000 kg/ha 50 years post-fire and exceeding 3000 kg/ha in stands 80-100 years post-fire ( Figure 6). The final lichen biomass map is displayed in Figure 7. Simple linear regression between observed and predicted lichen biomass at sampling locations (R 2 = 0.39) indicates our model performs to a similar standard as previous studies that created forage abundance layers for ungulates [18,25,26].  The final lichen biomass map is displayed in Figure 7. Simple linear regression between observed and predicted lichen biomass at sampling locations (R 2 = 0.39) indicates our model performs to a similar standard as previous studies that created forage abundance layers for ungulates [18,25,26]. The final lichen biomass map is displayed in Figure 7. Simple linear regression between observed and predicted lichen biomass at sampling locations (R 2 = 0.39) indicates our model performs to a similar standard as previous studies that created forage abundance layers for ungulates [18,25,26].

Discussion
We mapped lichen biomass across Woodland Caribou Provincial Park using a spatial modelling approach that can provide a framework to generate lichen biomass maps for resource management and ecological research in Canada's boreal forest. By relating our field observations of lichen abundance to forest structure, topographic and remote sensing attributes, we were able to identify environmental features useful in predicting ground lichens. We found that time-since-fire and ecosite were important predictors of ground lichens. Probability of occurrence and biomass of ground lichens was negatively associated with dense conifer ecosites ( Table 5, Table 7) and such stands demonstrated low lichen abundance (<1000 kg/ha) across the local time-since-fire continuum (Figure 6). Conversely, sparse conifer ecosites supported very low lichen abundance in the first 20 years after fire, but lichen biomass increased steadily from 20-50 years post-fire ( Figure 5). Mature sparse conifer (≥ 70 years old) supported approximately three times more lichen biomass than dense conifer of the same age ( Figure 6). The lichen abundance model appears to overestimate lichen biomass in young sparse conifer stands (0−19 years post-fire; Figure 6) relative to what was observed in the field ( Figure 5). Similarly, the model appears to exaggerate the accumulation of lichen biomass in older stands (≥ 50 years old). These discrepancies could be due to the unbalanced sampling design we employed, as we focused most of our sampling effort on middle-aged stands due to a secondary objective to test post-fire lichen recovery. This resulted in few observations at the young (0−19 years post-fire, n = 12) and older (50−119 years post-fire, n = 20) portions of the local post-fire continuum (Figure 3). In addition, our field observations suggest lichen biomass may follow a non-linear pattern with time-since-fire in sparse conifer ecosites ( Figure 5). We were unable to fully capture this trend in our analysis because generalized linear models assume a linear relationship between the response variable and the predictor variables. Other model types such as generalized additive models can improve predictions of non-linear trends [46]. Species distribution models [46] such as those developed through MaxEnt [49], provide a highly flexible workflow for mapping the distribution of plants, and have been used to map the presence of lichens [50]. Future research could incorporate these approaches to generate lichen maps for caribou conservation.
In our study, lichen presence was positively associated with canopy closure (Table 5). Conversely, lichen biomass was negatively associated with canopy closure (Table 7). Lichen growth is typically maximized at intermediate levels of canopy closure (~40%) [51], beyond which the growth of mosses is promoted at the expense of lichens [6]. Thus, lichens may require a minimum amount canopy closure to be present at a site but experience reduced growth at high levels of canopy closure, perhaps explaining the opposing responses of lichen presence and biomass observed here. In the oldest stands we sampled (70−119 years old), high mortality of mature trees created large gaps in the canopy and increased sun exposure at ground level. This promoted the growth of juniper shrubs (Juniperus communis L.), which often covered the ground lichens, possibly reducing access to foraging caribou. We had limited observations in overmature conifer stands (n = 14) and suggest future work should measure lichen biomass and caribou habitat selection in mature (50−70 years old) and overmature stands (≥70 years old) to estimate the optimal renewal period for caribou habitat. This information is essential to develop effective fire response and resource management plans that consider caribou conservation.
Most previous studies quantifying lichen over large areas only used remote sensing [9,18,21,22] or environmental [13,14] data. We anticipated that combining forest structure and topographic attributes with remote sensing attributes would provide the best results. Contrary to our expectations, models with only forest structure and/or topographic attributes were just as predictive as models including remote sensing attributes. The candidate models for both lichen presence and lichen abundance had small differences between AICc scores (∆AIC C < 2) ( Table 4; Table 6), which would typically indicate support for multiple candidate models [52]. However, the best candidate model for lichen presence and lichen abundance was the base model, the most parsimonious of the candidate set, only containing ecosite (i.e., dense conifer vs. sparse conifer), time-since-fire and canopy closure as covariates. The penalty weight assigned by AIC c to more complex models [52] indicates that the additional topographic and remote sensing covariates did not improve explanatory power over the base model. The lack of support for candidate models with remote sensing covariates could arise from multiple sources. First, environmental and remote sensing covariates are often correlated. We controlled for collinearity within models but, because we were interested in predicting lichen abundance rather than inferring ecological relationships, we did not account for correlation amongst models. Second, the coarse spatial resolution of Landsat imagery can cause trees to mask the spectral signature of ground lichens [14]. Keim et al. (2017) [9] report an R 2 = 0.74 for a lichen map generated using QuickBird satellite imagery (2.5 m pixels) and LiDAR data (1 m pixels). They found that QuickBird imagery predicted lichens better than SPOT (6 m pixels) and Landsat (30 m pixels) imagery in their study area in the continuous boreal forest of northeastern Alberta [9]. We suggest that the continued incorporation of finer resolution satellite [9] or UAV imagery [53] may help improve the accuracy of lichen mapping in years to come.
The lichen abundance map we generated in this study (Figure 7) highlights the patchy distribution of lichens on the landscape, which is driven primarily by the prevalence of fire in the study area. Lichen-rich forest is relatively restricted on the landscape, only occurring in mature, sparse conifer ecosites (≥50 years post-fire), where biomass often exceeds 3000 kg/ha ( Figure 5). Historically, many researchers used habitat type as a proxy for lichen abundance [54][55][56]. However, some studies have explicitly measured lichen availability and suggest caribou preferentially select stands with ≥3000 kg/ha of ground lichens as winter habitat [8,57,58]. Given that animal nutrition is necessarily related to the amount of available food, ecologists and land managers should strive to understand caribou's use of lichen biomass across time and space. In particular, identifying the use of lichen biomass during different seasons could be used to delineate nutritionally important habitat patches.
The relatively low accuracy of our map (R 2 = 0.39) is unsurprising given we used relatively coarse spatial covariates (30 m pixels) to model the presence and abundance of lichens, which are responding to environmental conditions at a very small scale (i.e., microsite). However, we feel that our lichen map provides a reasonable estimation of lichen biomass across our study area and suggest our modelling approach provides a useful framework for researchers to apply and improve in future lichen mapping projects. Most previous research mapped lichen cover [9,13,18], but we suggest lichen biomass is more biologically relevant than lichen cover as biomass is more closely related to animal energetics and fitness [26]. We stress the importance of validating biomass conversion factors and landscape covariates for new study areas, as growing conditions for lichen may vary. For example, in northern Alberta, peatlands are a dominant landscape feature. Previous studies indicate peatlands support much lower lichen abundance than upland sites [35], however raised 'islands' of drier peat within bogs can provide better conditions for lichen growth and support locally abundant ground lichens [9,11]. In other parts of the boreal forest, sandy areas dominated by jack pine support thick mats of ground lichens [33]. Integrating abiotic information, such as substrate type, groundwater depth and terrain ruggedness into spatial models may improve lichen predictive mapping, especially when the study area spans multiple biophysical regions. We incorporated data from numerous sources, data types and spatial resolutions to map the abundance of lichens in our study (Table 2). Researchers must be cognizant of the vintage of the source data in each layer they incorporate in their modelling framework to ensure temporal consistency. We suggest future research should focus on incorporating multiple sources of information, including time-since-fire and attributes derived from high resolution satellite imagery (e.g., spectral values, landcover type, forest structure [59]). This will improve spatio-temporal consistency and repeatability. We also encourage researchers to ensure they are selecting the most appropriate model for predicting the distribution of lichens and suggest utilizing generalized additive models [48] may be of particular utility to address some of the deficiencies of this study. Once a lichen abundance map has been generated, we encourage researchers to conduct independent validation using additional field sampling to improve certainty in their spatial predictions.

Conclusions
In this study, we propose a modelling framework for predicting the abundance of ground lichens in the boreal forest. We show that ecosite, time-since-fire and canopy closure are important drivers of lichen presence and abundance. We encourage researchers to use and improve our modelling framework to generate spatial predictions of lichen across caribou ranges. There is an increasing emphasis in wildlife ecology on including more biologically relevant variables in habitat selection analyses [60]. Quantifying nutritional landscapes can help researchers and managers measure how food availability changes with succession and varies by habitat type. Explicitly measuring selection for forage abundance can aid in the identification of high-quality habitat and to ensure continuous availability through resource planning and fire response. Mapping forage resources can also be used to test hypotheses, such as the effect of forage abundance on fitness or the tradeoff between nutrition and predator avoidance. Lichen abundance maps should be applied by researchers to help improve our understanding of caribou foraging ecology and support better conservation and resource management decisions.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A
To validate the cover-to-biomass conversion factors [33] for our study area, we collected a subsample of the lichen material in one 1 m 2 quadrat at 34 randomly selected sampling locations. Within each selected quadrat, we placed a 25 cm × 25 cm subplot and recorded the percent cover of all six lichen species. We then collected all thallus material of each Cladonia spp. ground lichen in the subplot, placing each species in a separate, labelled paper bag.
We air-dried our lichen samples after returning from the field to prevent mold and decomposition. We later cleaned the lichens of debris (moss, needles, etc.) and dried each sample in a biomass oven at 60 • C for 24 h. We weighed the dried samples using a digital scientific balance (measured in grams to two decimal places) and recorded a g/cm 2 value for each sample by dividing the weight of the dried sample (g) by the area it covered in the subplot (cm 2 ). We derived a cover-to-biomass conversion factor (g/cm 2 ) for each lichen species by taking the average g/cm 2 for all samples of each species.
We compared our conversion factors to McMullin's using two-tailed T-tests. We considered conversion factors not statistically different at an α-level = 0.05. The validation procedures could not be performed for C. stellaris or C. stygia because their rarity precluded them from being present in the destructive samples. There was considerable overlap in the conversion factors for each lichen species ( Figure A1). In the two-tailed T-test for each lichen species, the p-values (all ≥ 0.43) exceeded the α-level = 0.05, indicating the conversion factors developed by McMullin do not differ significantly from those recorded in this study. We therefore concluded that McMullin's conversion factors were appropriate for our study area and applied them to our subsequent analyses.  Figure A1). In the two-tailed T-test for each lichen species, the p-values (all ≥ 0.43) exceeded the αlevel = 0.05, indicating the conversion factors developed by McMullin do not differ significantly from those recorded in this study. We therefore concluded that McMullin's conversion factors were appropriate for our study area and applied them to our subsequent analyses. Figure A1. Mean biomass (g/cm²) recorded for each Cladonia spp. ground lichen destructively sampled by McMullin [33] and Silva (this study). Error bars represent ±1 SE.

Appendix B
To estimate lichen biomass in each 1 m 2 quadrat, we visually estimated the percent cover of each lichen species and converted to proportions (Table A1). We multiplied each proportion by 10,000 (10,000 cm 2 = 1 m 2 ) to determine the square centimetre area covered by each lichen species in the quadrat. We multiplied the square centimetre area covered by each lichen species by its conversion factor to derive a biomass estimate (g/m 2 ). We added the biomass of all species present in the quadrat to determine a biomass estimate for the quadrat (Table A1).

Appendix B
To estimate lichen biomass in each 1 m 2 quadrat, we visually estimated the percent cover of each lichen species and converted to proportions (Table A1). We multiplied each proportion by 10,000 (10,000 cm 2 = 1 m 2 ) to determine the square centimetre area covered by each lichen species in the quadrat. We multiplied the square centimetre area covered by each lichen species by its conversion factor to derive a biomass estimate (g/m 2 ). We added the biomass of all species present in the quadrat to determine a biomass estimate for the quadrat (Table A1). We repeat this procedure (Table A2) for all five 1 m 2 quadrats along the transect, resulting in five estimates of lichen biomass per sampling location (Table A2). To determine the stand-level lichen biomass for the sampling location (kg/ha), we add the biomass estimates for the five quadrats (g/5 m 2 ). We then divide by 1000 (1000 g = 1 kg) to convert to kilograms and multiply the result by 2000 (5 m 2 × 2000 = 10,000 m 2 = 1 ha), resulting in a stand-level biomass estimate (Table A2; kg/ha).

Appendix C
We used nine environmental covariates to construct a set of candidate models for predicting lichen presence and lichen abundance: ecosite, canopy closure, time-since-fire, elevation, slope, blue reflectance, short-wave infrared (SWIR2) reflectance, normalized difference vegetation index (NDVI) and normalized difference moisture index (NDMI). The pre-processing details for these datasets are described in the following sections.

C.1.1. Ecosite
We created an ecosite layer from the primary ecosite attribute (PRI_ECO) [61] for each polygon in the Forest Resource Inventory datasets for Woodland Caribou Provincial Park (2009) and the surrounding Forest Management Units: Kenora (2015), Red Lake (2013) and Whiskey Jack (2015) [39]. We grouped the 68 ecosites present in our study area into eleven broad categories: sparse conifer, dense conifer, anthropogenic, bog, fen, hardwood swamp, mixedwood, rock, shrubland and upland mixed conifer (Table A3). We assigned the simplified forest classification to each inventory dataset, merged them together, clipped them to the study area and created an ecosite raster in ArcGIS 10.5 [36]. We only sampled sparse conifer (ecosite B012) and dense conifer (ecosite B049) in our study, so the ecosite variable used for modelling was a factor with two levels: 1 = sparse conifer, 2 = dense conifer. The beta coefficient for ecosite in the lichen presence and lichen abundance models indicates the effect of dense conifer relative to sparse conifer. Unsampled ecosites were assigned 'NoData' in the lichen presence and abundance rasters.  We created our own time-since-fire layer using the fire perimeters captured into two provincial GIS polygon datasets: FiresByDecade   [40] and Fire Disturbance Area (1960-2013) [41]. We clipped the two datasets to the extent of the study area, merged them and converted the new layer to a raster format in ArcGIS 10.5 [36]. We calculated time-since-fire by subtracting the fire year from 2014 (the study year for producing the lichen map). Areas unaffected by fire since 1929 were assigned a uniform value of 100.

C.1.3. Canopy Closure
We created a canopy closure layer from the overstorey crown closure attribute (OCCLO) [61] for each polygon in the Forest Resource Inventory datasets for Woodland Caribou Provincial Park (2009) and the surrounding Forest Management Units: Kenora (2015), Red Lake (2013) and Whiskey Jack (2015) [39]. We merged the individual inventory datasets together, clipped them to the study area and created a single canopy closure raster in ArcGIS 10.5 [36].
Simple linear regression indicated poor agreement (R 2 adj = 0.17) between our canopy closure layer and our field observations. To improve the accuracy of our canopy closure layer, we used generalized linear models (family = Gamma, link = logit) to predict our field observations as a function four environmental covariates: canopy closure derived from the inventory datasets (OCCLO), ecosite (1 = sparse conifer, 2 = dense conifer), time-since-fire (TSF) and normalized difference vegetation index (NDVI) ( Table A4). Table A4. Name and structure of candidate models used to create the canopy closure layer for modelling lichen presence and abundance. Obs = canopy closure recorded in the field; OCCLO = canopy closure derived from the inventory datasets; TSF = time since fire (years), NDVI (normalized difference vegetation index). We ranked the candidate models by AIC c score [36] and considered the model with the lowest AIC c score as the best of the candidate set. The model with the lowest AIC c score included canopy closure derived from the inventory datasets, ecosite, time-since-fire and NDVI (Table A5). Table A5. Ranking of candidate models used to create the canopy closure layer for modelling lichen presence and abundance. k = number of fixed effects (+ 1 for intercept), w i = Akaike weight. TSF = time-since-fire, NDVI = normalized difference vegetation index. The model summary for the top model is presented in Table A6. We interpolated this top model across the study area using the raster package in R version 3.6.0 [37,38] to create the canopy closure layer we used to model lichen presence and abundance. Simple linear regression indicated the new canopy closure layer showed greater agreement with our field observations (R 2 adj = 0.40).

. Elevation & Slope
Elevation (metres above sea level) was obtained directly from a provincial digital elevation model [44]. Slope values were derived from the digital elevation model using ArcMap 10.5 [36].

C.1.5. Remote Sensing Covariates
We derived our remote sensing covariates from the spectral bands of two Landsat 8 Surface Reflectance datasets (captured July 31, 2014) [42,43]. The individual spectral bands used in this study are: (normalized difference moisture index) [20].